重读经典论文: Mean Value Coordinates for Closed Triangular Meshes
Mean Value Coordinates for Closed Triangular Meshes 发表于ACM SIGGRAPH 2003, 是图形学领域的经典论文，作者是之前组里毕业的居涛，论文只有6页，这在SIGGRAPH会议中并不多见，目前引用超过700次，在插值，参数化，变形等领域被广泛应用。
先看一下论文给的Teaser，这里主要是变形的一个应用。输入包括：待变形网格模型 + 控制网格（跟模型形状相似的稀疏网格），通过移动控制网格的顶点驱动模型形状发生改变，其实主要就是插值问题，将模型顶点表示为控制网格点的线性组合。
一、插值问题描述
本文所涉及的插值要做的事就是将2D多边形(3D多面体) PPP 内部的任意点 vvv 表示为PPP 的线性组合：
v=∑i=1nwipi(1)v = \sum^n_{i=1}w_ip_i \tag1v=i=1∑nwipi(1)
这里的 wiw_iwi 称为插值权重，满足 ∑i=1nwi=1\sum^n_{i=1}w_i = 1∑i=1nwi=1，插值问题的重点就是求解该权重。对于上图所示的变形问题，其实就是求出该权重并保持不变，然后修改控制网格: P→P′P \to P'P→P′, 最后，变形后的模型表示为：v′=∑i=1nwipi′v' = \sum^n_{i=1}w_ip'_iv′=∑i=1nwipi′ 。
二、2D多边形插值
这一部分主要探讨2D多边形内部插值，输入多边形 PPP, 内部一点 vvv, 输出对应的插值权重W={wi}W=\{w_i\}W={wi}。
如上图左所示，我们首先画一个以 vvv为圆心的单位圆，将 vvv 与各个多边形顶点连接的话，多边形每条边会对应到一段圆弧。
不失一般性，以 第一条边 p1p2p_1p_2p1p2 对应的圆弧为例，我们将vvv 指向这段圆弧上所有点的单位向量求和(积分)得到的向量 m1m_1m1 称为均值向量（mean vector）。 m1m_1m1 可以通过一组基 vp1=p1−v,vp2=p2−vvp_1=p_1-v,\ vp_2=p_2-vvp1=p1−v, vp2=p2−v 进一步表示为:
λ11(p1−v)+λ12(p2−v)=m1(2)\lambda^1_1(p_1-v) + \lambda^2_1(p_2-v) =m_1\tag2λ11(p1−v)+λ12(p2−v)=m1(2)
其中，λ11\lambda^1_1λ11 和 λ12\lambda^2_1λ12 分别表示 m1m_1m1 在这组基上的坐标值。另外，Floater曾在他的论文中给出了均值向量的表示形式：
Floater M S. Mean value coordinates[J]. Computer aided geometric design, 2003, 20(1): 19-27.
m1=tan[α/2]((p1−v)∣p1−v∣+(p2−v)∣p2−v∣)(3)m_1=\tan [\alpha / 2]\left(\frac{\left(p_{1}-v\right)}{\left|p_{1}-v\right|}+\frac{\left(p_{2}-v\right)}{\left|p_{2}-v\right|}\right)\tag3m1=tan[α/2](∣p1−v∣(p1−v)+∣p2−v∣(p2−v))(3)
其中，α\alphaα 是向量 vp1,vp2vp_1,vp_2vp1,vp2 之间的夹角。根据(2)(3)，我们就可以直接得出：
λ11=tan[α/2]∣p1−v∣,λ12=tan[α/2]∣p2−v∣(4)\lambda^1_1 = \frac{tan[\alpha/2]}{|p_1-v|}, \ \ \ \ \lambda^2_1= \frac{tan[\alpha/2]}{|p_2-v|}\tag4λ11=∣p1−v∣tan[α/2], λ12=∣p2−v∣tan[α/2](4)
那么，我们怎么根据上式导出插值权重呢？继续往下看。
注意（2）式只是针对第一条线段 p1p2p_1p_2p1p2 计算了均值向量 m1m_1m1，如果计算所有线段对应的均值向量，并对其求和，那么这些均值像之和为0向量。 相当于，圆心指向圆弧所有点构成的单位向量之和，必然为0，因为任意一个单位向量都存在一个方向相反的向量。也就是说：
{λ11(p1−v)+λ12(p2−v)=m1λ21(p2−v)+λ22(p3−v)=m2...λn1(pn−v)+λn2(p1−v)=mn\left\{ \begin{aligned} \lambda^1_1(p_1-v) + \lambda^2_1(p_2-v) =m_1\\ \lambda^1_2(p_2-v) + \lambda^2_2(p_3-v) =m_2\\ ...\\ \lambda^1_n(p_n-v) + \lambda^2_n(p_1-v) =m_n\\ \end{aligned} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧λ11(p1−v)+λ12(p2−v)=m1λ21(p2−v)+λ22(p3−v)=m2...λn1(pn−v)+λn2(p1−v)=mn
m1+m2+...+mn=0m_1 + m_2 + ... + m_n = 0m1+m2+...+mn=0
上述任意一组λi1\lambda^1_iλi1,λi2\lambda^2_iλi2均可通过（4）求解得到。因为任意 pi−vp_i - vpi−v 均出现两次（一个端点被两条边共享），因此重新整理上式为：
(λ11+λn2)(p1−v)+(λ12+λ21)(p2−v)+...+(λn−12+λn1)(pn−v)=0(5)( \lambda^1_1 + \lambda^2_n)(p_1-v) + (\lambda^2_1+\lambda^1_2)(p_2-v)+...+(\lambda^2_{n-1}+\lambda^1_n)(p_n-v) = 0\tag5(λ11+λn2)(p1−v)+(λ12+λ21)(p2−v)+...+(λn−12+λn1)(pn−v)=0(5)
我们令 k1=(λ11+λn2)k_1 =(\lambda^1_1 + \lambda^2_n)k1=(λ11+λn2), k2=(λ12+λ21)k_2=(\lambda^2_1+\lambda^1_2)k2=(λ12+λ21),…,kn=(λn−12+λn1)k_n = (\lambda^2_{n-1}+\lambda^1_n)kn=(λn−12+λn1), 上次可以进一步写为：
k1(p1−v)+k2(p2−v)+...+kn(pn−v)=0(6)k_1(p_1-v) +k_2(p_2-v)+...+k_n(p_n-v) = 0\tag6k1(p1−v)+k2(p2−v)+...+kn(pn−v)=0(6)
接下来，令 K=k1+k2+...+knK = k_1+k_2+...+k_nK=k1+k2+...+kn, w1=k1/K,w2=k2/K,...wn=kn/Kw_1=k_1/K, w_2=k_2/K, ... w_n=k_n/Kw1=k1/K,w2=k2/K,...wn=kn/K， 上式可以进一步改写为：
w1(p1−v)+w2(p2−v)+...+wn(pn−v)=∑i=1nwi(pi−v)=0(7)w_1(p_1-v) +w_2(p_2-v)+...+w_n(p_n-v) = \sum^n_{i=1}w_i(p_i-v)=0\tag7w1(p1−v)+w2(p2−v)+...+wn(pn−v)=i=1∑nwi(pi−v)=0(7)
我们发现 （7）跟（1）完全等价，因为 ∑i=1nwi=1\sum^n_{i=1}w_i = 1∑i=1nwi=1, v=∑i=1nwivv = \sum^n_{i=1}w_i vv=∑i=1nwiv。到此为止，我们得到了满足（1）式的插值权重也就是均值坐标。
三、3D多面体插值
在2D多边形中，我们通过将多边形投影到单位圆上求解均值坐标。容易想到，在3D多面体内部，我们可以将多面体投影到单位球来计算给定点的均值坐标。
如上图所示，我们将多面体投影到单位球面。如上图(3)所示，多面体上的任意三角面片对应单位球上的一个球面三角形，那么从vvv 指向该球面三角形上所有点的单位向量的积分称为该球面三角形的均值向量 mmm。文中给出了这一向量的计算方式（证明明详见论文）：
m=∑i=1312θini(8)m = \sum^3_{i=1}\frac{1}{2}\theta_in_i\tag8m=i=1∑321θini(8)
接下来，通过上图(5)所示的示例进行讲解，其均值向量可以表示为：
m=12(θ1n1+θ2n2+θ3n3)(9)m=\frac{1}{2}(\theta_1n_1 + \theta_2n_2 + \theta_3n_3)\tag9m=21(θ1n1+θ2n2+θ3n3)(9)
其中，θ1,θ2,θ3\theta_1,\theta_2,\theta_3θ1,θ2,θ3 分别是 (vp2,vp3),(vp1,vp3),(vp1,vp2)(vp_2,vp_3), (vp_1,vp_3),(vp_1,vp_2)(vp2,vp3),(vp1,vp3),(vp1,vp2) 的夹角，n1,n2,n3n_1,n_2,n_3n1,n2,n3 分别是 (vp2p3),(vp1p3),(vp1p2)(vp_2p_3), (vp_1p_3),(vp_1p_2)(vp2p3),(vp1p3),(vp1p2) 平面的法向。另一方面，球面三角形上的均值向量 mmm, 仍然可以通过一组基 (vp1,vp2,vp3)(vp_1,vp_2,vp_3)(vp1,vp2,vp3) 表示出来：
m=λ1(p1−v)+λ2(p2−v)+λ3(p3−v)(10)m =\lambda_1(p_1-v) + \lambda_2(p_2-v) + \lambda_3(p_3-v)\tag{10}m=λ1(p1−v)+λ2(p2−v)+λ3(p3−v)(10)
如何求解上述 λ\lambdaλ呢？因为n1n_1n1是平面 (vp2p3)(vp_2p_3)(vp2p3) 的法向量，所以 n1⋅(p2−v)=0n_1 \cdot (p_2-v) =0n1⋅(p2−v)=0, n1⋅(p3−v)=0n_1 \cdot (p_3-v) =0n1⋅(p3−v)=0，其它类似，那么有：
{n1m=n1λ1(p1−v)+0+0n2m=0+n2λ2(p2−v)+0n3m=0+0+n3λ3(p3−v)\left\{ \begin{aligned} n_1m =n_1\lambda_1(p_1-v) + 0+0\\ n_2m =0+n_2\lambda_2(p_2-v)+0 \\ n_3m =0+0+n_3\lambda_3(p_3-v)\\ \end{aligned} \right. ⎩⎪⎨⎪⎧n1m=n1λ1(p1−v)+0+0n2m=0+n2λ2(p2−v)+0n3m=0+0+n3λ3(p3−v)
所以：
{λ1=n1mn1(p1−v)λ2=n2mn2(p2−v)λ3=n3mn3(p3−v)(11)\left\{ \begin{aligned} \lambda_1 =\frac{n_1m}{n_1(p_1-v)} \\ \\ \lambda_2 =\frac{n_2m}{n_2(p_2-v)} \\ \\ \lambda_3 =\frac{n_3m}{n_3(p_3-v)} \\ \end{aligned} \right.\tag{11} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧λ1=n1(p1−v)n1mλ2=n2(p2−v)n2mλ3=n3(p3−v)n3m(11)
类似的，上式只是计算了多面体上一个三角面片对应的球面三角形的均值向量，以及对应的三个点的贡献。如果将所有三角面片投影到球面上，那么所有球面三角形的均值向量之和为0. 也就是：
∑T=1nλT1(pT1−v)+λT2(pT2−v)+λT3(pT3−v)=0(12)\sum^n_{T=1}\lambda^1_T(p^1_T-v) + \lambda^2_T(p^2_T-v) + \lambda^3_T(p^3_T-v) = 0\tag{12}T=1∑nλT1(pT1−v)+λT2(pT2−v)+λT3(pT3−v)=0(12)
其中，TTT表示遍历多面体 PPP 的所有三角形，pT1,pT2,pT3p^1_T,p^2_T,p^3_TpT1,pT2,pT3 表示某个三角形 TTT 的3个顶点。上式 λT1,λT2,λT3\lambda^1_T,\lambda^2_T,\lambda^3_TλT1,λT2,λT3 通过（11）式子计算。最后的最后，上式等价于：
∑i=1nwi(pi−v)=0\sum^n_{i=1}w_i(p_i-v)=0i=1∑nwi(pi−v)=0
这个过程跟2D的推导类似，最终归一化的权重 wiw_iwi 可以表示为：
wi=∑T∈N(i)∑j=13λTj⋅(pTi==pi)∑T=1n(λT1+λT2+λT3)(13)w_i =\frac{\sum_{T \in N(i)}\sum^3_{j=1}\lambda^j_T \cdot(p^i_T==p_i)}{\sum^n_{T=1}(\lambda^1_T+\lambda^2_T+\lambda^3_T)}\tag{13}wi=∑T=1n(λT1+λT2+λT3)∑T∈N(i)∑j=13λTj⋅(pTi==pi)(13)
也就是说，每个点的权重，跟所关联的所有三角形相关。 文中还有一些数值稳定性的优化方案，就不再描述。可以参看原文伪代码。
参考资料
[1] Floater M S. Mean value coordinates[J]. Computer aided geometric design, 2003, 20(1): 19-27.
[2] Tao Ju et al. Mean Value Coordinates for Closed Triangular Meshes
[3] 以及对应的slides： https://people.engr.tamu.edu/schaefer/teaching/689_Fall2006/talks/meanvalue.pdf