高階Laplace曲面形變算法（Polyharmonic Deformation）

　　數學上曲面的連續光滑形變可以通過最小化能量函數來建模得到，其中能量函數用來調節曲面的拉伸或彎曲程度，那么能量函數最小化同時滿足所有邊界條件的最優解就是待求曲面。

　　能量函數通常是二次函數形式：

其中S_*代表關於曲面參數u和v的k階偏導。

　　對於上述優化問題的求解方法，通常利用變分法得到對應的Euler-Lagrange方程，然后求解該方程得到最優解。對於二次能量函數形式，其對應的Euler-Lagrange方程為如下多階調和方程：

　　例如對於，那么對應的Euler-Lagrange方程就是2階Laplace方程Δ²x = 0。上式中Δ代表Laplace算子，邊界區域δΩ的限制條件b_j保證了系統存在非平凡解。

　　對於k = 1，上述方程描述的是滿足曲面面積最小的表面；對於k = 2，方程描述的是滿足曲面彎曲最小的表面，對於k = 3，方程描述的是滿足曲面曲率變化最小的表面。

　　方程的邊界條件可以選擇Dirichlet邊界條件（第一類邊界條件），即指定曲面邊界點的位置，或者Neumann 邊界條件（第二類邊界條件），即指定曲面邊界點的法向方向。

　　如果曲面用三角網格形式表示，那么我們在方程中需要使用離散Laplace算子，這是一個線性算子，具體形式在“網格形變算法”中有所解釋，其表達式如下：

其中α_ij和β_ij為邊e_ij對應的2個對角，N₁(p_i)表示頂點p_i的1環鄰域點，A(p_i)表示頂點p_i的Voronoi面積。

function L = Laplace_Matrix(V, F)
    nV = size(V,1);
    L = sparse(nV,nV);
    for i = 1:3
       i1 = mod(i-1,3)+1;
       i2 = mod(i  ,3)+1;
       i3 = mod(i+1,3)+1;
       Vec1 = V(F(:,i2),:) - V(F(:,i1),:);
       Vec2 = V(F(:,i3),:) - V(F(:,i1),:);

       cotangle = 0.5 .* dot(Vec1, Vec2, 2) ./ norm(cross(Vec1, Vec2, 2));

       L = L + sparse(F(:,i2), F(:,i3), cotangle, nV, nV);
    end
    L = L + L';
    L = -spdiags(full( sum(L,2) ), 0, nV, nV) + L;
end

　　利用離散線性Laplace算子，方程可以變化為稀疏線性方程組：

其中p = (p₁, … , p_p)是計算域內的待求點，對應下圖中黃色點集；h = (h₁, … , h_p)是計算域內的控制點，對應下圖中紅色和藍色點集，它們代表了方程中的邊界條件。需要注意的是為滿足求解得到的曲面C^k連續，方程中需要待求點外k+1環的控制點集作為邊界條件。

 

 

 

 

圖：k階能量方程使得曲面邊界滿足C^k-1連續

左：k = 1，中：k = 2，右：k = 3

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

 

 

參考文獻：

[1] Mario Botsch and Leif Kobbelt. 2004. An intuitive framework for real-time freeform modeling. ACM Trans. Graph. 23, 3 (August 2004), 630-634.