前言

之前寫(xiě)過(guò)徑向基函數(shù)(RBF)神經(jīng)網(wǎng)絡(luò)做分類(lèi)或者擬合。然后挖了個(gè)坑說(shuō)在《Phase-Functioned Neural Networks for Character Control》里面提到了用于做地形編輯，所以這篇博客就是解析一下如何用RBF做網(wǎng)格編輯系統(tǒng)。

參考博客：

Noe’s tutorial on deforming 3D geometry using RBFs

基于參考博客的人臉網(wǎng)格編輯code

有很多網(wǎng)格變形算法的python包PyGem

《Real-Time Shape Editing using Radial Basis Functions》

簡(jiǎn)介

為了更好理解什么是網(wǎng)格變形，直接看第二個(gè)參考博客的效果：

圖中實(shí)現(xiàn)的效果就是藍(lán)色的點(diǎn)是黃色人臉模型的幾個(gè)關(guān)鍵點(diǎn)，綠色的點(diǎn)是另一個(gè)沒(méi)被顯示的人臉模型的對(duì)應(yīng)人臉3D關(guān)鍵點(diǎn)，使用RBF變形算法，基于藍(lán)色到綠色的變換關(guān)系，將人臉變形到綠色關(guān)鍵點(diǎn)上。至于這個(gè)功能的用途請(qǐng)自行探索，后續(xù)有時(shí)間的話(huà)，搞不好我也會(huì)基于這個(gè)技術(shù)做個(gè)好玩的東東出來(lái)。

理論

其實(shí)理論基本就是參考博客的內(nèi)容，不過(guò)那個(gè)博客里面對(duì)徑向基函數(shù)的描述不是特別清晰。特別注意的是從PyGem工具包提供的RBF形變函數(shù)來(lái)看，RBF網(wǎng)格編輯會(huì)有個(gè)半徑參數(shù)，可以控制形變影響區(qū)域，其實(shí)就是在計(jì)算徑向基的時(shí)候，加了個(gè)額外的系數(shù)。

RBF插值

徑向基函數(shù)插值的作用是能夠在一些2D/3D離散點(diǎn)之間進(jìn)行平滑插值。假設(shè)在3維空間中有M個(gè)離散點(diǎn)$x_i$，RBF就能提供整個(gè)離散空間中的平滑插值函數(shù)。函數(shù)就是M個(gè)徑向基函數(shù)$g(r_i)$的結(jié)果之和，其中$r_i$是估算點(diǎn)和原始點(diǎn)的距離：
$$F(\mathbf{x})=\sum _{i=1}^M{a}_{i}g(||\mathbf{x}–{\mathbf{x}}_i||)+c_0+c_{1}x+c_{2}y+c_{3}z,\mathbf{x}=(x,y,z)?(1)$$
其中$a_i$是常量系數(shù)，后面四項(xiàng)$c_0$到$c_3$是一次多項(xiàng)式系數(shù)，這些項(xiàng)的就是無(wú)法單獨(dú)使用徑向基函數(shù)完成的一個(gè)仿射變換。

根據(jù)已有的M個(gè)離散值，可以構(gòu)建M個(gè)函數(shù)$F(x_i,y_i,z_i)=F_i$，然后基于此，構(gòu)建M+4個(gè)線(xiàn)性方程組$GA=F$，其中$F=(F_1,F_2,\dots ,F_M,0,0,0,0)$,$A=(a_1,a_2,\dots ,a_M,c0,c1,c2,c3)$，以及G是一個(gè)$(M+4)\times (M+4)$的矩陣
$$\mathbf{G}=?\begin{array}{cccccccccc}{g}_{11}& g_{12}& ?& ?& ?& g_{1M}& 1& x_1& y_1& z_1\\ g_{21}& g_{22}& ?& ?& ?& g_{2M}& 1& x_2& y_2& z_2\\ ?& ?& ?& ?& ?& ?& ?& ?& ?& ?\\ ?& ?& ?& ?& ?& ?& ?& ?& ?& ?\\ ?& ?& ?& ?& ?& ?& ?& ?& ?& ?\\ g_{M1}& g_{M2}& ?& ?& ?& g_{MM}& 1& x_M& y_M& z_M\\ 1& 1& ?& ?& ?& 1& 0& 0& 0& 0\\ x_1& x_2& ?& ?& ?& x_M& 0& 0& 0& 0\\ y_1& y_2& ?& ?& ?& y_M& 0& 0& 0& 0\\ z_1& z_2& ?& ?& ?& z_M& 0& 0& 0& 0\end{array}?$$
其中$g_{ij}=g(||x_i?x_j||)$，g有很多不同的選擇，同時(shí)也會(huì)導(dǎo)致不同的解。在參考博客中，使用了shift log function：
$$g(t)=\sqrt{\log (t^2+k^2)},k^2\ge 1$$
可以直接設(shè)置$k=1$，然后根據(jù)式(1)方程組，得到系數(shù)向量A。

其實(shí)上面的$g$就是傳說(shuō)中的徑向基函數(shù)了，根據(jù)scipy的文檔發(fā)現(xiàn)常用的徑向基函數(shù)有：

'multiquadric': sqrt((r/self.epsilon)**2 + 1) 'inverse': 1.0/sqrt((r/self.epsilon)**2 + 1) 'gaussian': exp(-(r/self.epsilon)**2) 'linear': r 'cubic': r**3 'quintic': r**5 'thin_plate': r**2 * log(r)

這些函數(shù)在PyGem里面有獨(dú)立的實(shí)現(xiàn)，只不過(guò)PyGem里面加入了半徑系數(shù)r控制插值影響區(qū)域，比如對(duì)gaussian的實(shí)現(xiàn)：

def gaussian_spline(X, r=1):result = np.exp(-(X * X) / (r * r))return result

詳細(xì)可自行查看PyGem源碼

偏移插值(Interpolating displacements)

上一節(jié)的插值是給了一系列的離散點(diǎn)，求解這些位于這些離散點(diǎn)間的值。而偏移插值的意思又是什么捏？假設(shè)M個(gè)3D坐標(biāo)點(diǎn)對(duì)應(yīng)的形變信息都知道了，也就是有另一個(gè)3D偏移向量$u_i$，將對(duì)應(yīng)索引的$x_i$進(jìn)行了偏移，此時(shí)就可以將$x_i$視為控制點(diǎn)，這些點(diǎn)被移動(dòng)到了$x_i+u_i$，RBF插值就是能夠計(jì)算其余剩下的3D點(diǎn)對(duì)應(yīng)的偏移量。

設(shè)${\mathbf{x}}_i=(x_i,y_i,z_i)$和${\mathbf{u}}_i=(u_i^x,u_i^y,u_i^z)$，那么同樣可以列出線(xiàn)性方程組：

$$\mathbf{G}{\mathbf{A}}_x=(u_1^x,u_2^x,\dots ,u_M^x,0,0,0,0{)}^T$$

$$\mathbf{G}{\mathbf{A}}_y=(u_1^y,u_2^y,\dots ,u_M^y,0,0,0,0{)}^T$$

$$\mathbf{G}{\mathbf{A}}_z=(u_1^z,u_2^z,\dots ,u_M^z,0,0,0,0{)}^T$$

其中$\mathbf{G}$和$\mathbf{A}$與上一節(jié)描述的差不多，只不過(guò)求解目標(biāo)從$F$變成了偏移量$u$。

求解完畢以后，需要對(duì)任意點(diǎn)x進(jìn)行偏移量的計(jì)算，這一點(diǎn)文章沒(méi)講，但是從源碼中分析發(fā)現(xiàn)就是$F$的計(jì)算公式：
$$F(\mathbf{x})=\sum _{i=1}^M{a}_{i}g(||\mathbf{x}–{\mathbf{x}}_i||)+c_0+c_{1}x+c_{2}y+c_{3}z,\mathbf{x}=(x,y,z)?(1)$$
這也難怪，畢竟偏移插值也是RBF理論，所以在即使求解目標(biāo)不同，但是形式相同的情況下，插值函數(shù)還是一樣的。

核心代碼解析

C++

求解階段

第二個(gè)參考博客的代碼就是第一個(gè)參考博客的c++實(shí)現(xiàn)，其中有兩個(gè)核心與上述理論對(duì)應(yīng)，代碼基于RBF插值理論實(shí)現(xiàn)了開(kāi)頭的人臉變形；

所以細(xì)節(jié)一：計(jì)算控制點(diǎn)的偏移量

// move control points for (unsigned int i = 0; i<controlPointPosX.size(); i++ ) {controlPointDisplacementX[i] = (destinationPointPosX[i]-controlPointPosX[i]);controlPointDisplacementY[i] = (destinationPointPosY[i]-controlPointPosY[i]);controlPointDisplacementZ[i] = (destinationPointPosZ[i]-controlPointPosZ[i]); }

源碼里面還有個(gè)系數(shù)(-cosf(morphtime)/2+0.5)不用管它，單純?yōu)榱孙@示變換過(guò)程設(shè)置的時(shí)間。

細(xì)節(jié)二計(jì)算形變系數(shù)，也就是公式1的代碼如下：

RBFInterpolator(vector<real> x, vector<real> y, vector<real> z, vector<real> f) {successfullyInitialized = false; // default value for if we end init prematurely.M = f.size();// all four input vectors must have the same length.if ( x.size() != M || y.size() != M || z.size() != M )return; ColumnVector F = ColumnVector(M + 4);P = Matrix(M, 3);Matrix G(M + 4,M + 4);// copy function valuesfor (unsigned int i = 1; i <= M; i++)F(i) = f[i-1];F(M+1) = 0; F(M+2) = 0; F(M+3) = 0; F(M+4) = 0;// fill xyz coordinates into P for (unsigned int i = 1; i <= M; i++){P(i,1) = x[i-1];P(i,2) = y[i-1];P(i,3) = z[i-1];}// the matrix below is symmetric, so I could save some calculations Hmmm. must be a todofor (unsigned int i = 1; i <= M; i++)for (unsigned int j = 1; j <= M; j++){real dx = x[i-1] - x[j-1];real dy = y[i-1] - y[j-1];real dz = z[i-1] - z[j-1];real distance_squared = dx*dx + dy*dy + dz*dz;G(i,j) = g(distance_squared);}//Set last 4 columns of Gfor (unsigned int i = 1; i <= M; i++){G( i, M+1 ) = 1;G( i, M+2 ) = x[i-1];G( i, M+3 ) = y[i-1];G( i, M+4 ) = z[i-1];}for (unsigned int i = M+1; i <= M+4; i++)for (unsigned int j = M+1; j <= M+4; j++)G( i, j ) = 0;//Set last 4 rows of Gfor (unsigned int j = 1; j <= M; j++){G( M+1, j ) = 1;G( M+2, j ) = x[j-1];G( M+3, j ) = y[j-1];G( M+4, j ) = z[j-1];}Try { Ginv = G.i(); A = Ginv*F;successfullyInitialized = true;}CatchAll { cout << BaseException::what() << endl; } }

這里調(diào)用的時(shí)候，$x,y,z$就是控制點(diǎn)的原始坐標(biāo)，$f$就是偏移量，比如對(duì)三個(gè)坐標(biāo)分別建立插值函數(shù)：

RBFX = RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementX ); RBFY = RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementY ); RBFZ = RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementZ );

整個(gè)過(guò)程就是先利用距離和徑向基函數(shù)構(gòu)建$(M+4)\times (M+4)$維度的$G$

for (unsigned int i = 1; i <= M; i++)for (unsigned int j = 1; j <= M; j++){real dx = x[i-1] - x[j-1];real dy = y[i-1] - y[j-1];real dz = z[i-1] - z[j-1];real distance_squared = dx*dx + dy*dy + dz*dz;G(i,j) = g(distance_squared);}

然后構(gòu)建$(M+4)$維度的$F$，再去求解$GA=F$。

推斷階段

當(dāng)任意的頂點(diǎn)過(guò)來(lái)以后，需要調(diào)用如下代碼：

interpolate(real x, real y, real z) {if (!successfullyInitialized)return 0.0f;real sum = 0.0f;// RBF partfor (unsigned int i = 1; i <= M; i++){real dx = x - P(i,1);real dy = y - P(i,2);real dz = z - P(i,3);real distance_squared = dx*dx + dy*dy + dz*dz;sum += A(i) * g(distance_squared);} //affine partsum += A(M+1) + A(M+2)*x + A(M+3)*y + A(M+4)*z;return sum; }

大致意思就是，需要將被估算頂點(diǎn)與M個(gè)已知形變頂點(diǎn)求基于距離的徑向基函數(shù)的加和，然后使用系數(shù)乘一下，就得到了新的坐標(biāo)，調(diào)用方法如下

void deformObject(TriangleMesh* res, TriangleMesh* initialObject, RBFInterpolator rbfX, RBFInterpolator rbfY, RBFInterpolator rbfZ) {for (unsigned int i = 0; i < res->getParticles().size(); i++){Vector3 oldpos = initialObject->getParticles()[i].getPos();Vector3 newpos;newpos[0] = oldpos[0] + rbfX.interpolate(oldpos[0], oldpos[1], oldpos[2]);newpos[1] = oldpos[1] + rbfY.interpolate(oldpos[0], oldpos[1], oldpos[2]);newpos[2] = oldpos[2] + rbfZ.interpolate(oldpos[0], oldpos[1], oldpos[2]);res->getParticles()[i].setPos(newpos);} }

python

求解階段

在PyGem代碼中有實(shí)現(xiàn)，由于python對(duì)矩陣的操作更加簡(jiǎn)潔，所以看起來(lái)很清晰

def _get_weights(self, X, Y):npts, dim = X.shapeH = np.zeros((npts + 3 + 1, npts + 3 + 1))H[:npts, :npts] = self.basis(cdist(X, X), self.radius)#, **self.extra)H[npts, :npts] = 1.0H[:npts, npts] = 1.0H[:npts, -3:] = XH[-3:, :npts] = X.Trhs = np.zeros((npts + 3 + 1, dim))rhs[:npts, :] = Yweights = np.linalg.solve(H, rhs)return weights

其中X，Y就是對(duì)應(yīng)控制點(diǎn)形變前后的坐標(biāo)位置，可以發(fā)現(xiàn)這里并沒(méi)有按照常理出牌，正常情況應(yīng)該就是對(duì)形變前后的差值做求解，不過(guò)事實(shí)證明，并不影響最終結(jié)果。

其中self.basis對(duì)應(yīng)的是幾種徑向基函數(shù)的一個(gè):

__bases = {'gaussian_spline': RBFFactory.gaussian_spline,'multi_quadratic_biharmonic_spline': RBFFactory.multi_quadratic_biharmonic_spline,'inv_multi_quadratic_biharmonic_spline': RBFFactory.inv_multi_quadratic_biharmonic_spline,'thin_plate_spline': RBFFactory.thin_plate_spline,'beckert_wendland_c2_basis': RBFFactory.beckert_wendland_c2_basis,'polyharmonic_spline': RBFFactory.polyharmonic_spline }

求解方式也是沒(méi)有用$G^{?1}F$，而是直接用np.linalg.solve求解了。

推斷階段

def __call__(self, src_pts):self.compute_weights()H = np.zeros((src_pts.shape[0], self.n_control_points + 3 + 1))H[:, :self.n_control_points] = self.basis(cdist(src_pts, self.original_control_points), self.radius)#**self.extra)H[:, self.n_control_points] = 1.0H[:, -3:] = src_ptsreturn np.asarray(np.dot(H, self.weights))

與求解階段一樣，先計(jì)算待預(yù)測(cè)坐標(biāo)點(diǎn)與控制點(diǎn)的距離，然后輸入到徑向基函數(shù)中，最后利用系數(shù)乘一下，就得到了待預(yù)測(cè)坐標(biāo)點(diǎn)的形變后位置。

實(shí)驗(yàn)

簡(jiǎn)單網(wǎng)格形變

藍(lán)色的標(biāo)記為規(guī)規(guī)整整的原始網(wǎng)格，黑色圓點(diǎn)為選取的形變點(diǎn)，紅色三角為目標(biāo)形變點(diǎn)，最終整個(gè)網(wǎng)格的形變結(jié)果就是紅色三角形了。

人頭形變

將c++兩個(gè)人頭拿過(guò)來(lái)，將控制點(diǎn)進(jìn)行著色后，使用meshlab可視化看看

在python中將PyGem的代碼摳出來(lái)，測(cè)試一下半徑影響：

當(dāng)形變半徑影響為0.5時(shí)候

當(dāng)形變半徑影響為1.0時(shí)候，就是標(biāo)準(zhǔn)的RBF插值了

將0.5和1.0的放一起對(duì)比一下：

可以發(fā)現(xiàn)半徑的設(shè)置的確對(duì)變形結(jié)果有一定影響。半徑較小的時(shí)候，鼻子都被捏塌了，半徑大點(diǎn)就好點(diǎn)，實(shí)際人臉變形的過(guò)程中，我們是不會(huì)使用這么簡(jiǎn)單的幾個(gè)點(diǎn)的，肯定會(huì)加入更多的關(guān)鍵點(diǎn)。

后記

好像理論非常的簡(jiǎn)單，就是一個(gè)方程組的求解，只不過(guò)方程組的建立與徑向基函數(shù)有關(guān)系。這玩意好像可以用來(lái)做捏臉什么的，以后有機(jī)會(huì)試試。

python版本的代碼可以去PyGem上找到，或者本博客的代碼在我的github上能找到，關(guān)注微信公眾號(hào)查看公眾號(hào)簡(jiǎn)介有，或者CSDN左側(cè)欄有g(shù)ithub網(wǎng)址。

總結(jié)

以上是生活随笔為你收集整理的径向基函数RBF三维网格变形的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。