三維網格細分算法（Catmull-Clark subdivision & Loop subdivision）附源碼(轉載)

轉載：  http://www.javashuo.com/article/p-eecwaznt-me.html

下圖描述了細分的基本思想，每次細分都是在每條邊上插入一個新的頂點，能夠看到隨着細分次數的增長，折線逐漸變成一條光滑的曲線。曲面細分須要有幾何規則和拓撲規則，幾何規則用於計算新頂點的位置，拓撲規則用於肯定新頂點的鏈接關係。下面介紹兩種網格細分方法：Catmull-Clark細分和Loop細分。

Catmull-Clark subdivision：

　　Catmull-Clark細分是一種四邊形網格的細分法則，每一個面計算生成一個新的頂點，每條邊計算生成一個新的頂點，同時每一個原始頂點更新位置。下圖爲Catmull-Clark細分格式的細分掩膜，對於新增長的頂點位置以及原始頂點位置更新規則以下：

1.網格內部F-頂點位置：

　　設四邊形的四個頂點爲v₀、v₁、v₂、v₃，則新增長的頂點位置爲v = 1/4*(v₀ + v₁ + v₂ + v₃)。

2.網格內部V-頂點位置：

　　設內部頂點v₀的相鄰點爲v₁、v₂，…，v_2n，則該頂點更新後位置爲，其中α、β、γ分別爲α = 1 - β - γ。

3.網格邊界V-頂點位置：

　　設邊界頂點v₀的兩個相鄰點爲v₁、v₂，則該頂點更新後位置爲v = 3/4*v₀ + 1/8*(v₁ + v₂)。

4.網格內部E-頂點位置：

　　設內部邊的兩個端點爲v₀、v₁，與該邊相鄰的兩個四邊形頂點分別爲v₀、v₁、v₂、v₃和v₀、v₁、v₄、v₅，則新增長的頂點位置爲v = 1/4*(v₀ + v₁ + v_f¹ + v_f²) = 3/8*(v₀ + v₁) + 1/16*(v₂ + v₃ + v₄ + v₅)。

5.網格邊界E-頂點位置：

　　設邊界邊的兩個端點爲v₀、v₁，則新增長的頂點位置爲v = 1/2*(v₀ + v₁)。

效果：

function [VV, FF, S] = CC_subdivision(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

Loop subdivision：

　　Loop細分是一種三角形網格的細分法則，它按照1-4三角形分裂，每條邊計算生成一個新的頂點，同時每一個原始頂點更新位置。下圖爲Loop細分格式的細分掩膜，對於新增長的頂點位置以及原始頂點位置更新規則以下：

1.網格內部V-頂點位置：

　　設內部頂點v₀的相鄰點爲v₁、v₂，…，v_n，則該頂點更新後位置爲，其中。

2.網格邊界V-頂點位置：

　　設邊界頂點v₀的兩個相鄰點爲v₁、v₂，則該頂點更新後位置爲v = 3/4*v₀ + 1/8*(v₁ + v₂)。

3.網格內部E-頂點位置：

　　設內部邊的兩個端點爲v₀、v₁，相對的兩個頂點爲v₂、v₃，則新增長的頂點位置爲v = 3/8*(v₀ + v₁) + 1/8*(v₂ + v₃)。

4.網格邊界E-頂點位置：

　　設邊界邊的兩個端點爲v₀、v₁，則新增長的頂點位置爲v = 1/2*(v₀ + v₁)。

效果：