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

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

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

function [VV, FF, S] = CCSubdivision(V, F, iter)  
    % Catmull_Clark subdivision  
    if ~exist('iter','var')  
        iter = 1;  
    end  
    VV = V;  
    FF = F;  
      
    for i = 1:iter   
       
        nv = size(VV,1);  
        nf = size(FF,1);   
        
        O = outline(FF);  
          
        original = 1:nv;  
        boundary = O(:,1)';  
        interior = original(~ismember(original, boundary));  
          
        no = length(original);  
        nb = length(boundary);  
        ni = length(interior);  
  
        %% Sv  
        Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);  
        [E, ~, idx] = unique(Etmp, 'rows');  
          
        Aeven = sparse([E(:,1) E(:,2)], [E(:,2) E(:,1)], 1, no, no);  
        Aodd = sparse([FF(:,1) FF(:,2)], [FF(:,3) FF(:,4)], 1, no, no);  
        Aodd = Aodd + Aodd';  
          
        val_even = sum(Aeven,2);  
        beta = 3./(2*val_even);  
          
        val_odd = sum(Aodd,2);  
        gamma = 1./(4*val_odd);  
          
        alpha = 1 - beta - gamma;  
          
        Sv = sparse(no,no);  
        Sv(interior,:) = ...  
            sparse(1:ni, interior, alpha(interior), ni, no) + ...  
            bsxfun(@times, Aeven(interior,:), beta(interior)./val_even(interior)) + ...  
            bsxfun(@times, Aodd(interior,:), gamma(interior)./val_odd(interior));  
        Sboundary = ...  
            sparse([O(:,1);O(:,2)],[O(:,2);O(:,1)],1/8,no,no) + ...  
            sparse([O(:,1);O(:,2)],[O(:,1);O(:,2)],3/8,no,no);  
        Sv(boundary,:) = Sboundary(boundary,:);  
          
        %% Sf  
        Sf = 1/4 .* sparse(repmat((1:nf)',1 ,4), FF, 1);  
        i0 = no + (1:nf)';  
          
        %% Se  
        flaps = sparse([idx;idx], ...  
                       [FF(:,3) FF(:,4);FF(:,4) FF(:,1);FF(:,1) FF(:,2);FF(:,2) FF(:,3)], ...  
                       1);  
        onboundary = (sum(flaps,2) == 2);  
        flaps(onboundary,:) = 0;  
          
        ne = size(E,1);  
        Se = sparse( ...  
                [1:ne 1:ne]', ...  
                [E(:,1); E(:,2)], ...  
                [onboundary;onboundary].*1/2 + ~[onboundary;onboundary].*3/8, ...  
                ne, ...  
                no) + ...  
                flaps*1/16;  
          
        %% new faces & new vertices  
        i1 = no +   nf + (1:nf)';  
        i2 = no + 2*nf + (1:nf)';  
        i3 = no + 3*nf + (1:nf)';  
        i4 = no + 4*nf + (1:nf)';  
          
        FFtmp = [i0 i4 FF(:,1) i1; ...  
                 i0 i1 FF(:,2) i2; ...  
                 i0 i2 FF(:,3) i3; ...  
                 i0 i3 FF(:,4) i4];  
  
        reidx = [(1:no)'; no+(1:nf)'; no+nf+idx];  
        FF = reidx(FFtmp);  
          
        S = [Sv; Sf; Se];  
        VV = S*VV;  
    end  
 end  

其中outline函数如下

function out = outline( FF )
%OUTLINE Summary of this function goes here
%   Detailed explanation goes here
Etmp = sort([FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)],2);
[~, ~, idx] = unique(Etmp, 'rows');

oriEtmp = [FF(:,1) FF(:,2);FF(:,2) FF(:,3);FF(:,3) FF(:,4);FF(:,4) FF(:,1)];
hh=sortrows([oriEtmp,idx],3);

x2=diff(sortrows(idx));
vector = all(x2==0, 2);

index1=find(vector);
index2=index1+1;
index=[index1;index2];

hh(index,:)=[];
out=hh(:, 1:2);

end

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

[V,F]=obj__read('six.obj');
V=V';F=F';
iter=4;
[VV, FF] = CCSubdivision(V, F, iter);
%[VV, FF] = CCsub(V, F, iter);
obj_write('six1.obj',VV',FF');

[V,F]=obj__read('torus.obj');
V=V';F=F';
iter=4;
[VV, FF] = CCSubdivision(V, F, iter);
%[VV, FF] = CCsub(V, F, iter);
obj_write('torus1.obj',VV',FF');

最后贴上细分效果

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

搞定哈哈^^