三維網格補洞算法（Radial Basis Function）

　　在逆向工程中，由於設備或模型的原因，我們獲取得到的三維模型數據往往並不完整，從而使得生成的網格模型存在孔洞，這對后續的模型分析會造成影響。下面介紹一種基於徑向基函數（RBF：Radial Basis Function）的三角網格補洞方法。

Step 1：檢測孔洞邊界

　　三角網格是由一系列頂點（V）以及由這些頂點所構成的三角面片（F）所組成，由三角面片可以得到網格的邊（E）。通常一條邊連接兩個三角面片，這種邊稱為網格內部邊，而如果某條邊僅連接一個三角面片，那么稱這條邊為網格邊界邊，所有的邊界邊按順序連接之后就形成了網格的孔洞。

Step 2：初始化網格

　　為了使孔洞填充簡單、健壯，可以采用最小角度法進行網格修補，具體步驟如下：

　　（1）得到孔洞邊界點信息，計算邊界邊長度的平均值l。

　　（2）計算每個邊界點的兩條相鄰邊的夾角大小。

　　（3）找出具有最小夾角的邊界點，計算它的兩個相鄰邊界點的距離s，判斷s < 2×l是否成立：若成立，則按下圖所示增加一個三角形，若不成立則增加兩個三角形。

　　（4）更新邊界點信息。

　　（5）判斷孔洞是否修補完整，若未完整，轉（2），否則結束。

Step 3：最小二乘網格

　　初始化補洞得到的網格質量不是很好，我們需要優化網格頂點的位置，優化的條件是：

其中d_i為頂點v_i的1環鄰域頂點數。

　　上式可以用一個線性方程組來描述：LV = 0，其中L是Laplace矩陣，具體形式為，矩陣L的秩等於n – k，n為網格頂點的數目，k為網格連通區域的個數，如果網格是全連通的，那么矩陣L的秩是n – 1。因此如果我們要使方程有解，需要加入m（m ≥ k）個控制頂點坐標v_s作為方程的邊界條件，在實際中我們將初始化網格的邊界頂點作為控制頂點。

　　其實上述線性方程組的求解等價於如下能量函數最小化的求解：

　　能量函數的第一部分是使得網格頂點盡量光滑，即每個頂點位於其1環鄰域頂點的中心，第二部分是為了控制頂點的位置滿足要求。

　　最小二乘網格的優勢是能夠生成高質量的光滑網格，生成過程僅需要網格的拓撲連接關系和少數控制點的坐標信息。

Step 4：徑向基函數（RBF）隱式曲面

　　徑向基函數是一個僅依賴於離控制點c距離的函數，即h(x,c) = h(||x - c||)，距離通常是歐式距離，也可以是其他形式距離。徑向基函數網絡是一個三層BP網絡，其可以表示為多個基函數的線性組合：

　　下面列出2個最常用的基函數形式：

　　（1）Gaussian：h(r) = e^-(εr)2

　　（2）Polyharmonic spline：h(r) = r^k，k = 1、3、5、… 或h(r) = r^kln(r)，k = 2、4、6、…

　　利用徑向基函數網絡並通過給定控制點x_i和對應的值f_i之后，可以求解得到網絡中的系數λ_i，因此徑向基函數網絡能夠解決空間散亂數據點的平滑插值問題，函數的零等值面就是我們要求的曲面。

　　在實際求解時函數f(x)表達式中通常會增加一個一次多項式P(x)，如下：

，其中：

　　將控制點位置x_i和值f_i代入函數f(x)可以得到如下方程，求解之后即可確定隱式曲面f(x)。

　　控制點x_i可分為三類：

　　（1）邊界控制點：曲面通過的點，f(x_i) = 0

　　（2）外部控制點：將邊界控制點沿法向正方向移動一小段距離而得到的控制點，取f(x_i) = -1

　　（3）內部控制點：將邊界控制點沿法向負方向移動一小段距離而得到的控制點，取f(x_i) = 1

　　計算時我們選擇的基函數為h(r) = r³，其實此時隱式曲面函數滿足Δ³f = 0，將基函數代入后得到隱式曲面的表達式如下：

function options = RBF_Create(x, y, type)
    % x: input data; y: input data value
    options = [];
    options.nodes = x;
    switch type
        case 'linear'
            options.phi = @rbfphi_linear;
        case 'cubic'
            options.phi = @rbfphi_cubic;
    end    
    
    [n, dim] = size(x);
    A = zeros(n,n);
    for i = 1:n
        r = normrow(bsxfun(@minus, x(i,:), x));
        temp = feval(options.phi, r);
        A(i,:) = temp';
        A(:,i) = temp;
    end
    % Polynomial part
    P = [ones(n,1) x];
    A = [ A      P
          P' zeros(dim+1, dim+1)];
    B = [y; zeros(dim+1, 1)];
    
    coeff = A\B;
    options.coeff = coeff;
end

% Radial Base Functions
function u = rbfphi_linear(r)
    u = r;
end
function u = rbfphi_cubic(r)
    u = r.*r.*r;
end

Step 5：牛頓插值

　　為了得到插值隱式曲面的網格，我們需要把最小二乘網格的頂點投射到隱式曲面上，這里我們采用牛頓迭代法：

　　最小二乘網格的頂點位置作為迭代初始條件，當||x_k+1 – x_k|| < ε時，迭代停止，ε為設定的誤差。上式中梯度表達式如下：

效果：

 本文為原創，轉載請注明出處：http://www.cnblogs.com/shushen

 

 

參考文獻：

[1] Olga Sorkine and Daniel Cohen-Or. 2004. Least-Squares Meshes. In Proceedings of the Shape Modeling International 2004 (SMI '04). IEEE Computer Society, Washington, DC, USA, 191-199.

[2] Greg Turk and James F. O'Brien. 1999. Shape transformation using variational implicit functions. In Proceedings of the 26th annual conference on Computer graphics and interactive techniques (SIGGRAPH '99). ACM Press/Addison-Wesley Publishing Co., New York, NY, USA, 335-342.

[3] 袁天然. 三角網格模型光順簡化和縫補技術的研究及應用[D]. 南京航空航天大學, 2007.