上期我們與大家探討了倍受業界歡迎的mesh deformation with ARAP energy [Sorkine and Alexa 2007]。它真正做到了我們一直所講的「在保持局部細節的前提下去編輯模型的大體形態」這一要求，而且易於實現，計算量較小，並且算法保證會收斂。但是，ARAP也存在諸如無法保持模型的體積，以及當模型很精細時解算會不太穩定等缺陷，畢竟這世上並不存在完美的工程技術或方法，只有特定情形下最合適的某種選擇。
本期我們給大家講一講浙江大學的黃勁教授在2006年CG頂會SIGGRAPH上發表的Subspace gradient domain mesh deformation [Huang et al. 2006]，它首先為待編輯的模型創建一個包裹在其外部但頂點數很少的control mesh，然後將定義在原mesh上的nonlinear least squares energy（其中包括保持體積不變的限制條件）投影到這個control mesh所在的subspace中，並用高斯 - 牛頓方法（Gauss-Newton method）去求解。由於是在變量數很少並且被平滑過的subspace中求解，該方法具有較好的穩定性，且算法收斂更快。
Control Mesh如上圖的馬和牛頭模型，由灰色三角形拼接而成的就是原始的待編輯模型，其外部包裹的透明網格就是新建的control mesh。control mesh在mesh deformation中就像木偶戲中的吊繩，通過對關鍵部位的控制提供了一種更加簡易的動作編輯方法。在工程問題的語境中，就是減少了自由度（degree of freedom），或降低了問題的維度（dimension reduction，因為變量數更少了）。
製作control mesh的方法有很多種。比如像[Sander et al. 2000]中描述的不斷刪減原始模型上的三角形和邊，並保證每次刪減後得到的新mesh在三維空間中都包裹著原始模型，而且還要和原始模型外形相近，直到模型足夠簡易（如下圖第二行）。或者也可以先將原始模型上的頂點都向模型外（朝該處的法向）移動一段距離，得到「膨脹」了的模型後再來刪減三角形和邊（這時可以更加容易地保證通過刪減得到的新模型在三維空間中包裹著原始模型，但「膨脹」了的模型可能會產生一些自交）。
*注1：上文提到的刪減模型的三角形和邊並且保持刪減前後模型外形相近這一問題叫做mesh simplification（如下圖，此時不涉及前文提到的包裹體積的限定條件），也是數字幾何處理（digital geometry processing）方向曾經的一個研究熱點，它能將三維掃描得到的過於精細的模型適當簡化以降低後續存儲與處理的功耗，不同精細程度的模型也可用於表示大型虛擬場景中不同遠近的物體以節約處理時間。個人覺得其中最重要的兩項工作是Progressive Meshes [Hoppe 1996]和QSlim [Garland and Heckbert 1997]。前者的主要貢獻是提供了一種支持隨意增減三角形和邊的數據結構，使得用戶可以根據視覺效果不斷調整所需的模型精細程度；後者的主要貢獻是定義了一種非常簡潔有效的度量，來近似描述每次刪減操作前後模型之間的差異。該問題用到的數學方法也是optimization，有興趣可以去瞧一瞧。
Subspace Projection有了control mesh，接下來我們需要定義原始模型和control mesh兩者形變間的關係，也就是原始模型將怎樣被control mesh的形變驅動。假設control mesh有個頂點，原始模型有個頂點，且。那麼將它們所有頂點的坐標分別壘在一起可用高維向量表示成和。考慮最簡單的線性映射，我們實際上就是要找到一個矩陣，使得
這樣每當control mesh的頂點坐標在維的subspace中求得後，原始模型的頂點坐標即可由上式求出。
其實上式就是在把原始模型上每個頂點的坐標都表示為control mesh上頂點坐標的線性組合。我們知道，在三維空間中，最少只需用四個不共面的點的線性組合就可表示任意一點的坐標，而這裡雖說，但肯定至少也是在左右這個量級上的，所以並不唯一，故我們尋求一個帶有一定的幾何特徵的儘可能平滑的矩陣：其平滑性表現在，對每一行，在空間中相近的control mesh上的頂點具有相近的係數，並且對每一列，在空間上相近的原始模型上的頂點在control mesh上同一頂點處的係數也要儘可能相近；其幾何特徵表現在，對每一行，距離待表示的原始模型上的頂點越近的control mesh的頂點，其係數越大。拓展到triangle mesh上的mean value coordinates [Ju et al. 2005]正好滿足這些特徵。之所以要使得該矩陣平滑，是因為它可以使得接下來會講到的優化求解過程更加穩定。
*注2：值得注意的一點是，本文所講的方法並不同於free-form mesh deformation [Sederberg and Parry 1986]那樣直接根據control mesh定義待優化的能量函數，而是仍然像前幾期所講的方法一樣根據原始模型定義能量函數（詳見下一章節），只不過在求解時將能量函數映射到了control mesh所在的subspace中對進行求解，從而縮小了解的搜索範圍。而像[Sederberg and Parry 1986]那樣在control mesh上定義的能量函數並不一定能很好地刻畫原始模型的形狀特徵。舉個最簡單的例子：如果能量函數中含有保持體積不變的項，那麼由[Sederberg and Parry 1986]方法得到的control mesh的體積確實能夠保持不變，但以它映射得到的原始模型體積很可能會變，因為映射並不具有保持映射前後模型體積變化量的特徵，這也是為什麼本文所講的方法得到的control mesh體積可能有所變化，但以它映射得到的原始模型的體積能夠保持原樣，如下圖對比：
帶有嚴格限定條件的非線性Least Squares能量函數注意，數學公式要來啦～！
[Huang et al. 2006]原文中的能量函數總共有4項，分別是局部細節保持限定、體積保持限定、基於視角投影的位置限定、和模型骨架限定。所謂限定，其實是指等式，[Huang et al. 2006]的思路是通過最小化Constrained Least Squares讓這些等式儘量或嚴格被滿足。方便起見，我們只講前兩項，並直接採用前兩期中普通的位置限定條件來處理用戶滑鼠輸入。那麼整個問題就是求解
其中包含局部細節保持限定和普通的位置限定，它們將儘量被滿足，因為這兩者都嚴格被滿足的情況基本不存在；為體積保持限定，它將嚴格被滿足。
由本系列的第一篇我們知道，描述局部細節的能量函數一定要具有對局部旋轉的容忍性（rotation invariance），這裡我們重新推導一種非線性的rotation invariant局部細節描述算子:：對原始模型上的一個頂點，假設和它直接相連的鄰點為（），那麼該點處的Laplacian coordinates（詳見本系列第一篇中的介紹）可以由該處的各鄰接三角形的法向的線性組合來表示：
其中係數有多種滿足條件的值。其實細心的讀者可能會發現，若把當成未知數，上式就是一個有無窮多個解的線性方程組，故我們可通過QR分解求出滿足上式的L2 norm最小的，並作為新的局部細節描述子的參數。由於局部旋轉不僅旋轉了，各三角形的法向量也經歷了相同的旋轉，所以一旦確定了，無論局部經歷何種旋轉，將上式中的變量換成形變後模型的坐標後依然成立。但是同理，局部的各種伸縮變換也會使得上式成立，故我們將頂點處的局部細節保持限定項定義為
其中為編輯後原始模型上頂點處的Laplacian coordinates，為編輯後模型的頂點坐標，是編輯後模型上點處鄰接三角形的法線的線性組合。這裡引入了normalization，可以使得局部的伸縮變換被避免掉。而開頭列出的能量函數中的就是計算Laplacian coordinates的Laplacian matrix（詳見本系列第一篇中的介紹），對應上式的左邊，對應上式的右邊。
至於保體積的限定，對原始模型這樣一個封閉的triangle mesh，其體積可表示為
即求出三維坐標系原點與每個三角形所組成的四面體的signed volume並求和。能量函數中的為形變前原始模型的體積。
我們可以看到，在這個問題中，除了用戶選取頂點的位置限定外，其它兩個限定都是非線性的，所以我們要求解的是一個Constrained Nonlinear Least Squares問題。
結合Line Search的Gauss-Newton方法這裡我們採用Gauss-Newton方法來求解nonlinear least squares問題，其思路和我們在上一期中提到的類似，就是通過不斷在局部進行近似求解來搜索最優值，即每一步都將高度非線性的能量函數用二次型能量來逼近，從而將問題轉化為一系列線性方程組求解問題。
對每一步的求解，我們可將在當前的X處進行泰勒展開並丟掉二階以上的項：，從而在局部把整個nonlinear least squares問題轉換成一個linear least squares問題：
注意到這裡還是有一個不太好對付的非線性的限定條件，故我們如法炮製，將其在局部近似為，並應用拉格朗日乘子法將這個帶有線性限定條件的linear least squares問題轉化為一個不帶限定條件的linear least squares問題，從而通過求解線性方程組得到每一步的搜索向量：
其中，分別為和的一階微分矩陣（相對於標量能量函數而言是二階微分信息），為拉格朗日乘子。於是我們就可以由編輯前的原始模型頂點坐標開始，每次求解一個linear least squares問題得到搜索向量，並以更新模型頂點坐標從而開始下一步，直至最終收斂（求得的趨近於時）。至於投影到subspace中求解，在數學上其實就是change of variables，即將本章節所有表達式中的替換成，對進行求解，然後再由重建出。故優化問題就變成了
每一步中搜索向量的求解式也就變成了：
由於在每一步中（比如第步），是依據被近似後的能量函數求得的，故只是該近似能量的最優解，而並不一定會保證原本的非線性能量函數值真的有所減小。為了確保新得到的確實更接近非線性能量的最優解，我們引入這個參數，即step size，它可以通過line search來確定。[Huang et al. 2006]中採用了最基本的基於Armijo Condition的backtracking line search，就是在每一步中將從開始不斷折半，直到檢測出用當前的計算得到的能夠使原本的非線性能量函數值有足夠的減小為止。這裡我就不詳細列舉公式了，有需要的話可以參考我當時實現這篇論文時寫的report。
*注3：求解nonlinear least squares問題常用的方法還有Levenberg–Marquardt method，它相較於Gauss-Newton方法而言在每一步的線性系統的係數矩陣中多加了一個damping項（均勻的擴大稀疏矩陣對角線上的係數），這樣可以使得求解更穩定，但收斂會慢一些。如果damping加得很大，掩蓋了原本的二階微分係數，Levenberg–Marquardt方法就會退化成gradient descent方法（只採用一階微分信息，在machine learning這種數據量極大的問題中應用很廣）。這些方法都屬於hill-climbing method，它們常用於求解連續能量函數的最優化問題，詳見convex optimization [Boyd and Vandenberghe 2004]。
Space Reduction在show結果之前我還想再多說一點，就是為什麼這樣的subspace mesh deformation會work，其實是有深層次的解釋的。雖然通過subspace projection我們搜索的範圍被限定小了，但它只是省去了對應於改變局部細節的那些自由度，而這些自由度我們正好不需要，因為我們就是要保持局部細節不變，只編輯大體形態。這在數學上，可以用原始模型的Laplacian matrix的奇異值分解來解釋：
假定編輯前後模型的頂點坐標分別為和，且其變化量為。如果用來表示Laplacian matrix，那麼編輯前後模型的Laplacian coordinate的變化可表示為。我們對進行奇異值分解：，並將用中的列向量表示為，其中為線性組合的係數，那麼編輯前後模型的Laplacian coordinates的變化量可由其L2 norm計算並表示為，其中為矩陣中的奇異值。在本系列的第一篇中我們說過，一個triangle mesh的Laplacian coordinates表徵了它的局部細節特徵，也就是我們要保持的信息，換句話說，我們希望它在經歷了編輯後的變化要儘量小，也就是要儘量小。那麼根據之前的奇異值分解，只要僅存在於對應奇異值很小的那些奇異向量所組成的空間中，我們的目的就達到了。而通過本文所講的subspace projection，正是把我們要搜索的限定在了這樣一個子空間中。
這種思想其實來源於信號處理，應用到這裡就是把一個三維物體的形體信息當成了一種信號來看待，算是數字幾何處理領域中比較獨特的一種觀念，詳見Spectral Geometry Processing with Manifold Harmonics [Vallet and Lévy 2008]：
效果比較與分析除了求解線性方程組，該方法在實現上較前兩期所講的方法而言主要難在非線性能量函數的求導上。像這種高度非線性的項，求導的最好方法就是應用鏈式法則了，即將它拆解成多個簡單函數的複合或四則運算，然後逐步求導。為了確保求導以及其實現的正確性，我在實現這個方法時還採用了finite difference的方法來進行驗證，詳見我的report和code（其中還包含了UBC數字幾何處理課程的編程作業）。
那麼現在先來看看原文[Huang et al. 2006]秀出的結果吧：
我們可以發現這裡能支持的模型姿態變化已經非常廣泛了，基本達到了skinning [Jacobson et al. 2014]的程度，而且相較之下還能將表面細節和體積保持的更好，只是運算時間會很長。事實上skinning也算是一種subspace mesh deformation，它的subspace就是由骨架上各頂點的坐標組成的。由於骨架中的骨骼是二維的線段，故這個subspace完全無法處理對體積的保持，所以會出現裹糖紙和關節坍縮等問題，但正因為它的subspace極其小，所以能夠做到實時。
原文還將算法與不進行subspace projection直接在完整的高維空間中求解的情況進行了比較：
我們可以由上圖的上半部分的模型編輯中間結果和能量函數值折線圖看到，在完整的高維空間中求解是非常不穩定而且會收斂到很不自然的結果的。
我在實現該方法時，還將它與我們上期講的ARAP [Sorkine and Alexa 2007]進行了比較，發現了本方法在保持體積和引入全局旋轉上明顯的優勢：
上圖說明本方法在保持編輯前後模型的體積上較ARAP有明顯優勢，但它的這種體積保持並不在任何時候都會產生很好的效果，因為被保持的體積是整個模型的體積，所以下圖的結果就這楊在我測試的過程中不經意地發生了：
馬的頭就這樣被拉長了...所以個人認為更好的方式是去引入模型的各組成部分（例如上圖馬模型可分為頭、軀幹、四肢）的語義信息，並分別保持各部位的體積。
我們上期說到過，如果沒有很好的局部旋轉初始化，ARAP會在有大改動的編輯中給出不太自然的結果，如下圖中的對比：
我們看到，ARAP收斂到了一個不太自然的局部最優，而本方法收斂到了最直接的解，對全局引入了合適的旋轉。（注意這裡兩者的求解都是一次操作完成的，並沒有用上期所講的那種滑鼠拖拽不斷求解的方法。）
結語再往後筆者就要進行deadline前最後的衝刺了，故將由origami dance在本系列的第四篇給大家介紹volumetric mesh deformation method，並把我們在這個系列中所介紹的這類geometric mesh deformation method做個總結。本系列的前三篇中所處理的模型都是surface mesh，它們只記錄了物體表面的形態，其體內是空的。而volumetric mesh，顧名思義，就是模型內部也是有離散採樣點的，這可以使形變過程中的體積保持更加可行。之所以稱這些方法為geometric method，是說這裡我們對於模型的編輯是瞬間完成的，從原始模型直接跳到編輯完的模型，沒有計算任何動態信息（如各離散採樣點的運動速度等）。這點有區別於基於物理規律來計算三維模型的彈性形變[Terzopoulos et al. 1987]，它是一個對加速度以時間來不斷積分的過程，和我們在電影工業中的流體模擬系列中所講的方法類似。另外值得注意的一點是，如果在物理模擬的框架中應用生物學物理模型（biomechanics）來模擬肌肉和肌腱的形變，得到的特效角色動畫將會非常真實，比如Musculotendon Simulation for Hand Animation [Shinjiro et al. 2008]。除了基於幾何和物理的方法之外，還可以用數據驅動（data-driven、example-based，比如[Jones et al. 2016]）以及計算共形幾何（Computational Conformal Geometry）[Gu and Yau 2008]的方法來進行mesh deformation，以後有機會我們可以再開新的系列一起探討。
參考文獻[Sorkine and Alexa 2007] Sorkine, O., & Alexa, M. (2007, July). As-rigid-as-possible surface modeling. In Symposium on Geometry processing (Vol. 4).
[Huang et al. 2006] Huang, J., Shi, X., Liu, X., Zhou, K., Wei, L. Y., Teng, S. H., ... & Shum, H. Y. (2006, July). Subspace gradient domain mesh deformation. In ACM Transactions on Graphics (TOG) (Vol. 25, No. 3, pp. 1126-1134). ACM.
[Sander et al. 2000] Sander, P. V., Gu, X., Gortler, S. J., Hoppe, H., & Snyder, J. (2000, July). Silhouette clipping. In Proceedings of the 27th annual conference on Computer graphics and interactive techniques (pp. 327-334). ACM Press/Addison-Wesley Publishing Co..
[Hoppe 1996] Hoppe, H. (1996, August). Progressive meshes. In Proceedings of the 23rd annual conference on Computer graphics and interactive techniques (pp. 99-108). ACM.
[Garland and Heckbert 1997] Garland, M., & Heckbert, P. S. (1997, August). Surface simplification using quadric error metrics. In Proceedings of the 24th annual conference on Computer graphics and interactive techniques (pp. 209-216). ACM Press/Addison-Wesley Publishing Co..
[Ju et al. 2005] Ju, T., Schaefer, S., & Warren, J. (2005, July). Mean value coordinates for closed triangular meshes. In ACM Transactions on Graphics (TOG) (Vol. 24, No. 3, pp. 561-566). ACM.
[Sederberg and Parry 1986] Sederberg, T. W., & Parry, S. R. (1986). Free-form deformation of solid geometric models. ACM SIGGRAPH computer graphics, 20(4), 151-160.
[Boyd and Vandenberghe 2004] Boyd, S., & Vandenberghe, L. (2004). Convex optimization. Cambridge university press.
[Vallet and Lévy 2008] Vallet, B., & Lévy, B. (2008, April). Spectral geometry processing with manifold harmonics. In Computer Graphics Forum (Vol. 27, No. 2, pp. 251-260). Blackwell Publishing Ltd.
[Jacobson et al. 2014] Jacobson, A., Deng, Z., Kavan, L., & Lewis, J. (2014). Skinning: Real-time shape deformation. In ACM SIGGRAPH.
[Terzopoulos et al. 1987] Terzopoulos, D., Platt, J., Barr, A., & Fleischer, K. (1987, August). Elastically deformable models. In ACM Siggraph Computer Graphics (Vol. 21, No. 4, pp. 205-214). ACM.
[Shinjiro et al. 2008] Sueda, S., Kaufman, A., & Pai, D. K. (2008, August). Musculotendon simulation for hand animation. In ACM Transactions on Graphics (TOG) (Vol. 27, No. 3, p. 83). ACM.
[Jones et al. 2016] Jones, B., Thuerey, N., Shinar, T., & Bargteil, A. W. (2016). Example-based plastic deformation of rigid bodies. ACM Transactions on Graphics (TOG), 35(4), 34.
[Gu and Yau 2008] Gu, X. D. (2008). Computational conformal geometry. S. T. Yau (Ed.). Somerville, Mass, USA: International Press.
_(:3」∠)_ _(・ω・」∠)_ _(:з)∠)_ ∠( ᐛ 」∠)＿ _(:зゝ∠)_
請毫不猶豫地關注我們：
我們的網站：https://graphicon.io
知乎專欄：GraphiCon圖形控
公眾號：GraphiCon
如果你有什麼想法，建議，或者想加入我們，你可以：
給我們發郵件：hi@graphicon.io
加入我們的QQ群：SIQGRAPH（342086343）
加入我們的slack群：GraphiCon
GraphiCon長期接受投稿，如果你想投稿給我們可以通過上面的方式聯繫我們！
本作品採用知識共享署名-非商業性使用-禁止演繹 4.0 國際許可協議進行許可。