MATLAB实现Catmull-Clark细分（CC细分）

终于调试好了Catmull-Clark细分（CC细分）的全部程序，将之前只适用于封闭四边形网格的程序进行了完善

主要一段代码来自于三维网格细分算法（Catmull-Clark subdivision & Loop subdivision）附源码，这个博主的很多篇博文都写的非常好，但是经常丢三落四的，像在这篇博文中他就用到了函数outline.m用来计算网格的边界，但是博主却没有给出outline函数，我自己重新编写了这个函数，并且能够成功执行^^，现在我贴出完整代码

 1 function [VV, FF, S] = CCSubdivision(V, F, iter)  
 2     % Catmull_Clark subdivision  
 3     if ~exist('iter','var')  
 4         iter = 1;  
 5     end  
 6     VV = V;  
 7     FF = F;  
 8       
 9     for i = 1:iter   
10        
11         nv = size(VV,1);  
12         nf = size(FF,1);   
13         
14         O = outline(FF);  
15           
16         original = 1:nv;  
17         boundary = O(:,1)';  
18         interior = original(~ismember(original, boundary));  
19           
20         no = length(original);  
21         nb = length(boundary);  
22         ni = length(interior);  
23   
24         %% Sv  
25         Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);  
26         [E, ~, idx] = unique(Etmp, 'rows');  
27           
28         Aeven = sparse([E(:,1) E(:,2)], [E(:,2) E(:,1)], 1, no, no);  
29         Aodd = sparse([FF(:,1) FF(:,2)], [FF(:,3) FF(:,4)], 1, no, no);  
30         Aodd = Aodd + Aodd';  
31           
32         val_even = sum(Aeven,2);  
33         beta = 3./(2*val_even);  
34           
35         val_odd = sum(Aodd,2);  
36         gamma = 1./(4*val_odd);  
37           
38         alpha = 1 - beta - gamma;  
39           
40         Sv = sparse(no,no);  
41         Sv(interior,:) = ...  
42             sparse(1:ni, interior, alpha(interior), ni, no) + ...  
43             bsxfun(@times, Aeven(interior,:), beta(interior)./val_even(interior)) + ...  
44             bsxfun(@times, Aodd(interior,:), gamma(interior)./val_odd(interior));  
45         Sboundary = ...  
46             sparse([O(:,1);O(:,2)],[O(:,2);O(:,1)],1/8,no,no) + ...  
47             sparse([O(:,1);O(:,2)],[O(:,1);O(:,2)],3/8,no,no);  
48         Sv(boundary,:) = Sboundary(boundary,:);  
49           
50         %% Sf  
51         Sf = 1/4 .* sparse(repmat((1:nf)',1 ,4), FF, 1);  
52         i0 = no + (1:nf)';  
53           
54         %% Se  
55         flaps = sparse([idx;idx], ...  
56                        [FF(:,3) FF(:,4);FF(:,4) FF(:,1);FF(:,1) FF(:,2);FF(:,2) FF(:,3)], ...  
57                        1);  
58         onboundary = (sum(flaps,2) == 2);  
59         flaps(onboundary,:) = 0;  
60           
61         ne = size(E,1);  
62         Se = sparse( ...  
63                 [1:ne 1:ne]', ...  
64                 [E(:,1); E(:,2)], ...  
65                 [onboundary;onboundary].*1/2 + ~[onboundary;onboundary].*3/8, ...  
66                 ne, ...  
67                 no) + ...  
68                 flaps*1/16;  
69           
70         %% new faces & new vertices  
71         i1 = no +   nf + (1:nf)';  
72         i2 = no + 2*nf + (1:nf)';  
73         i3 = no + 3*nf + (1:nf)';  
74         i4 = no + 4*nf + (1:nf)';  
75           
76         FFtmp = [i0 i4 FF(:,1) i1; ...  
77                  i0 i1 FF(:,2) i2; ...  
78                  i0 i2 FF(:,3) i3; ...  
79                  i0 i3 FF(:,4) i4];  
80   
81         reidx = [(1:no)'; no+(1:nf)'; no+nf+idx];  
82         FF = reidx(FFtmp);  
83           
84         S = [Sv; Sf; Se];  
85         VV = S*VV;  
86     end  
87  end  

其中outline函数如下

 

 1 function out = outline( FF )
 2 %OUTLINE Summary of this function goes here
 3 %   Detailed explanation goes here
 4 Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);
 5 [~, ~, idx] = unique(Etmp, 'rows');
 6 
 7 oriEtmp = [FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)];
 8 hh=sortrows([oriEtmp,idx],3);
 9 
10 x2=diff(sortrows(idx));
11 vector = all(x2==0, 2);
12 
13 index1=find(vector);
14 index2=index1+1;
15 index=[index1;index2];
16 
17 hh(index,:)=[];
18 out=hh(:, 1:2);
19 
20 end

对于任意四边形网格是适用的，下面是我们的测试代码

 1 [V,F]=obj__read('six.obj');
 2 V=V';F=F';
 3 iter=4;
 4 [VV, FF] = CCSubdivision(V, F, iter);
 5 %[VV, FF] = CCsub(V, F, iter);
 6 obj_write('six1.obj',VV',FF');
 7 
 8 [V,F]=obj__read('torus.obj');
 9 V=V';F=F';
10 iter=4;
11 [VV, FF] = CCSubdivision(V, F, iter);
12 %[VV, FF] = CCsub(V, F, iter);
13 obj_write('torus1.obj',VV',FF');

最后贴上细分效果

 

开放四边形网格
	六面开口盒子	四次细分
封闭四边形网格
	四边形框	四次细分

搞定哈哈^^

转载于:https://www.cnblogs.com/lafengxiaoyu/p/7575021.html