南宁本地有几家网站开发,discuz论坛 整合到网站,广州地址设计网站,网络浏览器前言
之前写过径向基函数(RBF)神经网络做分类或者拟合。然后挖了个坑说在《Phase-Functioned Neural Networks for Character Control》里面提到了用于做地形编辑#xff0c;所以这篇博客就是解析一下如何用RBF做网格编辑系统。
参考博客#xff1a;
Noe’s tutorial on d…前言
之前写过径向基函数(RBF)神经网络做分类或者拟合。然后挖了个坑说在《Phase-Functioned Neural Networks for Character Control》里面提到了用于做地形编辑所以这篇博客就是解析一下如何用RBF做网格编辑系统。
参考博客
Noe’s tutorial on deforming 3D geometry using RBFs
基于参考博客的人脸网格编辑code
有很多网格变形算法的python包PyGem
《Real-Time Shape Editing using Radial Basis Functions》
简介
为了更好理解什么是网格变形直接看第二个参考博客的效果 图中实现的效果就是蓝色的点是黄色人脸模型的几个关键点绿色的点是另一个没被显示的人脸模型的对应人脸3D关键点使用RBF变形算法基于蓝色到绿色的变换关系将人脸变形到绿色关键点上。至于这个功能的用途请自行探索后续有时间的话搞不好我也会基于这个技术做个好玩的东东出来。
理论
其实理论基本就是参考博客的内容不过那个博客里面对径向基函数的描述不是特别清晰。特别注意的是从PyGem工具包提供的RBF形变函数来看RBF网格编辑会有个半径参数可以控制形变影响区域其实就是在计算径向基的时候加了个额外的系数。
RBF插值
径向基函数插值的作用是能够在一些2D/3D离散点之间进行平滑插值。假设在3维空间中有M个离散点xix_ixiRBF就能提供整个离散空间中的平滑插值函数。函数就是M个径向基函数g(ri)g(r_i)g(ri)的结果之和其中rir_iri是估算点和原始点的距离 F(x)∑i1Maig(∣∣x–xi∣∣)c0c1xc2yc3z,x(x,y,z)⋯(1)F(\mathbf{x}) \sum_{i1}^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) \cdots \mathbf{(1)} F(x)i1∑Maig(∣∣x–xi∣∣)c0c1xc2yc3z,x(x,y,z)⋯(1) 其中aia_iai是常量系数后面四项c0c_0c0到c3c_3c3是一次多项式系数这些项的就是无法单独使用径向基函数完成的一个仿射变换。
根据已有的M个离散值可以构建M个函数F(xi,yi,zi)FiF(x_i,y_i,z_i)F_iF(xi,yi,zi)Fi然后基于此构建M4个线性方程组GAFGAFGAF其中F(F1,F2,…,FM,0,0,0,0)F(F_1,F_2,\ldots,F_M,0,0,0,0)F(F1,F2,…,FM,0,0,0,0),A(a1,a2,…,aM,c0,c1,c2,c3)A(a_1,a_2,\ldots,a_M,c0,c1,c2,c3)A(a1,a2,…,aM,c0,c1,c2,c3)以及G是一个(M4)×(M4)(M4)\times(M4)(M4)×(M4)的矩阵 G[g11g12∙∙∙g1M1x1y1z1g21g22∙∙∙g2M1x2y2z2∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙gM1gM2∙∙∙gMM1xMyMzM11∙∙∙10000x1x2∙∙∙xM0000y1y2∙∙∙yM0000z1z2∙∙∙zM0000]\mathbf{G} \left[\begin{array}{cccccccccc}g_{11} g_{12} \bullet \bullet \bullet g_{1M} 1 x_1 y_1 z_1 \\ g_{21} g_{22} \bullet \bullet \bullet g_{2M} 1 x_2 y_2 z_2 \\ \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \\ \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \\ \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \bullet \\ g_{M1} g_{M2} \bullet \bullet \bullet g_{MM} 1 x_M y_M z_M \\ 1 1 \bullet \bullet \bullet 1 0 0 0 0 \\ x_1 x_2 \bullet \bullet \bullet x_M 0 0 0 0 \\ y_1 y_2 \bullet \bullet \bullet y_M 0 0 0 0 \\ z_1 z_2 \bullet \bullet \bullet z_M 0 0 0 0\end{array}\right] G⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡g11g21∙∙∙gM11x1y1z1g12g22∙∙∙gM21x2y2z2∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙∙g1Mg2M∙∙∙gMM1xMyMzM11∙∙∙10000x1x2∙∙∙xM0000y1y2∙∙∙yM0000z1z2∙∙∙zM0000⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤ 其中gijg(∣∣xi−xj∣∣)g_{ij}g(||x_i-x_j||)gijg(∣∣xi−xj∣∣)g有很多不同的选择同时也会导致不同的解。在参考博客中使用了shift log function g(t)log(t2k2),k2≥1g(t)\sqrt{\log{(t^2k^2)}},k^2\geq1 g(t)log(t2k2),k2≥1 可以直接设置k1k1k1然后根据式(1)方程组得到系数向量A。
其实上面的ggg就是传说中的径向基函数了根据scipy的文档发现常用的径向基函数有
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)这些函数在PyGem里面有独立的实现只不过PyGem里面加入了半径系数r控制插值影响区域比如对gaussian的实现
def gaussian_spline(X, r1):result np.exp(-(X * X) / (r * r))return result详细可自行查看PyGem源码
偏移插值(Interpolating displacements)
上一节的插值是给了一系列的离散点求解这些位于这些离散点间的值。而偏移插值的意思又是什么捏假设M个3D坐标点对应的形变信息都知道了也就是有另一个3D偏移向量uiu_iui将对应索引的xix_ixi进行了偏移此时就可以将xix_ixi视为控制点这些点被移动到了xiuix_iu_ixiuiRBF插值就是能够计算其余剩下的3D点对应的偏移量。
设xi(xi,yi,zi)\mathbf{x}_i (x_i, y_i, z_i)xi(xi,yi,zi)和ui(uix,uiy,uiz)\mathbf{u}_i (u^x_i, u^y_i, u^z_i)ui(uix,uiy,uiz)那么同样可以列出线性方程组
GAx(u1x,u2x,…,uMx,0,0,0,0)T\mathbf{G} \mathbf{A}_x (u^x_1, u^x_2, \ldots, u^x_M, 0, 0, 0, 0)^TGAx(u1x,u2x,…,uMx,0,0,0,0)T
GAy(u1y,u2y,…,uMy,0,0,0,0)T\mathbf{G} \mathbf{A}_y (u^y_1, u^y_2, \ldots, u^y_M, 0, 0, 0, 0)^TGAy(u1y,u2y,…,uMy,0,0,0,0)T
GAz(u1z,u2z,…,uMz,0,0,0,0)T\mathbf{G} \mathbf{A}_z (u^z_1, u^z_2, \ldots, u^z_M, 0, 0, 0, 0)^TGAz(u1z,u2z,…,uMz,0,0,0,0)T
其中G\mathbf{G}G和A\mathbf{A}A与上一节描述的差不多只不过求解目标从FFF变成了偏移量uuu。
求解完毕以后需要对任意点x进行偏移量的计算这一点文章没讲但是从源码中分析发现就是FFF的计算公式 F(x)∑i1Maig(∣∣x–xi∣∣)c0c1xc2yc3z,x(x,y,z)⋯(1)F(\mathbf{x}) \sum_{i1}^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) \cdots \mathbf{(1)} F(x)i1∑Maig(∣∣x–xi∣∣)c0c1xc2yc3z,x(x,y,z)⋯(1) 这也难怪毕竟偏移插值也是RBF理论所以在即使求解目标不同但是形式相同的情况下插值函数还是一样的。
核心代码解析
C
求解阶段
第二个参考博客的代码就是第一个参考博客的c实现其中有两个核心与上述理论对应代码基于RBF插值理论实现了开头的人脸变形
所以细节一计算控制点的偏移量
// move control points
for (unsigned int i 0; icontrolPointPosX.size(); i )
{controlPointDisplacementX[i] (destinationPointPosX[i]-controlPointPosX[i]);controlPointDisplacementY[i] (destinationPointPosY[i]-controlPointPosY[i]);controlPointDisplacementZ[i] (destinationPointPosZ[i]-controlPointPosZ[i]);
}源码里面还有个系数(-cosf(morphtime)/20.5)不用管它单纯为了显示变换过程设置的时间。
细节二计算形变系数也就是公式1的代码如下
RBFInterpolator(vectorreal x, vectorreal y, vectorreal z, vectorreal 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(M1) 0; F(M2) 0; F(M3) 0; F(M4) 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, M1 ) 1;G( i, M2 ) x[i-1];G( i, M3 ) y[i-1];G( i, M4 ) z[i-1];}for (unsigned int i M1; i M4; i)for (unsigned int j M1; j M4; j)G( i, j ) 0;//Set last 4 rows of Gfor (unsigned int j 1; j M; j){G( M1, j ) 1;G( M2, j ) x[j-1];G( M3, j ) y[j-1];G( M4, j ) z[j-1];}Try { Ginv G.i(); A Ginv*F;successfullyInitialized true;}CatchAll { cout BaseException::what() endl; }
}这里调用的时候x,y,zx,y,zx,y,z就是控制点的原始坐标fff就是偏移量比如对三个坐标分别建立插值函数
RBFX RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementX );
RBFY RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementY );
RBFZ RBFInterpolator(controlPointPosX, controlPointPosY, controlPointPosZ, controlPointDisplacementZ );整个过程就是先利用距离和径向基函数构建(M4)×(M4)(M4)\times(M4)(M4)×(M4)维度的GGG
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);}然后构建(M4)(M4)(M4)维度的FFF再去求解GAFGAFGAF。
推断阶段
当任意的顶点过来以后需要调用如下代码
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(M1) A(M2)*x A(M3)*y A(M4)*z;return sum;
}大致意思就是需要将被估算顶点与M个已知形变顶点求基于距离的径向基函数的加和然后使用系数乘一下就得到了新的坐标调用方法如下
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代码中有实现由于python对矩阵的操作更加简洁所以看起来很清晰
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其中XY就是对应控制点形变前后的坐标位置可以发现这里并没有按照常理出牌正常情况应该就是对形变前后的差值做求解不过事实证明并不影响最终结果。
其中self.basis对应的是几种径向基函数的一个:
__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
}求解方式也是没有用G−1FG^{-1}FG−1F而是直接用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))与求解阶段一样先计算待预测坐标点与控制点的距离然后输入到径向基函数中最后利用系数乘一下就得到了待预测坐标点的形变后位置。
实验
简单网格形变
蓝色的标记为规规整整的原始网格黑色圆点为选取的形变点红色三角为目标形变点最终整个网格的形变结果就是红色三角形了。 人头形变
将c两个人头拿过来将控制点进行着色后使用meshlab可视化看看 在python中将PyGem的代码抠出来测试一下半径影响
当形变半径影响为0.5时候 当形变半径影响为1.0时候就是标准的RBF插值了 将0.5和1.0的放一起对比一下 可以发现半径的设置的确对变形结果有一定影响。半径较小的时候鼻子都被捏塌了半径大点就好点实际人脸变形的过程中我们是不会使用这么简单的几个点的肯定会加入更多的关键点。
后记
好像理论非常的简单就是一个方程组的求解只不过方程组的建立与径向基函数有关系。这玩意好像可以用来做捏脸什么的以后有机会试试。
python版本的代码可以去PyGem上找到或者本博客的代码在我的github上能找到关注微信公众号查看公众号简介有或者CSDN左侧栏有github网址。
