根据变换矩阵的不同，可以分为三种变形方法，分别是：仿射变换、相似变换、刚性变换。其中刚性变换的效果是最好的，这边主要讲最简单的仿射变换，其余的只要替换一下变换矩阵即可实现。
一般的坐标变换
坐标变换一般分为：绕中心的旋转变换 + 平移变换。
l v ( v ) = ( v − p ∗ ) M + q ∗ l_v(v) = (v - p_*)M + q_* lv(v)=(v−p∗)M+q∗
一般的坐标变换只要一个方程，进行线性变换就可以了。
在多个控制点下的仿射变换
很显然，之前的公式并不适用于多个控制点，所以我们要利用最小二乘法了。其中 p p p为变换前的控制点， q q q为变换后的控制点。
每个控制点对输入的影响权重如下，其中 α \alpha α一般取1：
w i = 1 ∣ p i − v ∣ 2 α w_i = \frac{1}{|p_i -v|^{2\alpha}} wi=∣pi−v∣2α1
变换中心 p ∗ ， q ∗ p_*，q_* p∗，q∗计算如下：
p ∗ = ∑ ( w i p i ) ∑ w i p_* = \frac{\sum (w_ip_i)}{\sum w_i} p∗=∑wi∑(wipi) ， q ∗ = ∑ ( w i q i ) ∑ w i q_* = \frac{\sum (w_iq_i)}{\sum w_i} q∗=∑wi∑(wiqi)
我们要求的是，在该权重分配下，变换矩阵M是什么。
由于根据影响权重的不同，我们要求距离近的控制点影响要大一些，故最终的加权最小二乘表达式如下：
∑ w i ∣ p ˉ i M − q ˉ i ∣ 2 \sum w_i|\bar p_i M - \bar q_i|^2 ∑wi∣pˉiM−qˉi∣2
在线性代数中，我们求解形如 A x = b Ax=b Ax=b的正规方程最小二乘解，可得：
x = ( A T A ) − 1 A T b x = (A^TA)^{-1}A^Tb x=(ATA)−1ATb
故有 M = ( ∑ p i T w i p i ) − 1 ∑ w i p i T q i M = (\sum p_i^T w_i p_i)^{-1} \sum w_i p_i^T q_i M=(∑piTwipi)−1∑wipiTqi
有了旋转矩阵M的表达式后，我们得到变形的表达式：
l v ( v ) = ( v − p ∗ ) M + q ∗ = ( v − p ∗ ) ( ∑ p i T w i p i ) − 1 ∑ w i p i T q i + q ∗ l_v(v) = (v - p_*)M + q_* = (v - p_*)(\sum p_i^T w_i p_i)^{-1} \sum w_i p_i^T q_i + q_* lv(v)=(v−p∗)M+q∗=(v−p∗)(∑piTwipi)−1∑wipiTqi+q∗
CVfpoint MLS(CVfpoint t , CVVECTORS_F32 pList, CVVECTORS_F32 qList)
{
CVfpoint fv;
double A[2][2], B[2][2], M[2][2];
int num = pList->dataNum;
float *w = CV_malloc0(num * sizeof(float));
CVfpoint* p = pList->datap;
CVfpoint* q = qList->datap;
//计算各个控制顶点的权重，也就是计算点t到各个顶点的距离1/sqr(d)
for (int i = 0; i < num; i++)
{
float temp;
if (p[i].x != t.x || p[i].y != t.y)
temp = 1 / ((p[i].x - t.x) * (p[i].x - t.x) + (p[i].y - t.y) * (p[i].y - t.y));
else//如果t为控制顶点，那么需要把该控制顶点的权重设置为无穷大
temp = 100000;
w[i] = temp;
}
CVfpoint pc, qc;
//q为目标图像的控制点的位置，我们的目标是找到t在q中的对应位置
double px = 0, py = 0, qx = 0, qy = 0, tw = 0;
for (int i = 0; i < num; i++)
{
px += w[i] * p[i].x;//所有控制顶点p的加权位置
py += w[i] * p[i].y;
qx += w[i] * q[i].x;//所有控制顶点q的加权位置
qy += w[i] * q[i].y;
tw += w[i];//总权重
}
pc.x = px / tw;
pc.y = py / tw;
qc.x = qx / tw;
qc.y = qy / tw;
//初始化
for (int i = 0; i < 2; i++)
for (int j = 0; j < 2; j++)
{
A[i][j] = 0;
B[i][j] = 0;
M[i][j] = 0;
}
for (int i = 0; i < num; i++)
{
double P[2] = { p[i].x - pc.x,p[i].y - pc.y };//P
double PT[2][1];//P^T
PT[0][0] = p[i].x - pc.x;
PT[1][0] = p[i].y - pc.y;
double Q[2] = { q[i].x - qc.x,q[i].y - qc.y };//Q
double T[2][2];
//T = P^T * P : 列乘行=矩阵
T[0][0] = PT[0][0] * P[0];
T[0][1] = PT[0][0] * P[1];
T[1][0] = PT[1][0] * P[0];
T[1][1] = PT[1][0] * P[1];
//A = wi * (P^T * P)
for (int k = 0; k < 2; k++)
for (int j = 0; j < 2; j++)
{
A[k][j] += w[i] * T[k][j];
}
// T = P^T * Q
T[0][0] = PT[0][0] * Q[0];
T[0][1] = PT[0][0] * Q[1];
T[1][0] = PT[1][0] * Q[0];
T[1][1] = PT[1][0] * Q[1];
//B = wi * (P^T * Q)
for (int k = 0; k < 2; k++)
for (int j = 0; j < 2; j++)
{
B[k][j] += w[i] * T[k][j];
}
}
CV_free0(w);
// 求 A^-1
double det = A[0][0] * A[1][1] - A[0][1] * A[1][0];
if (det < 0.0000001)
{
fv.x = t.x + qc.x - pc.x;
fv.y = t.y + qc.y - pc.y;
return fv;
}
double temp1, temp2, temp3, temp4;
temp1 = A[1][1] / det;
temp2 = -A[0][1] / det;
temp3 = -A[1][0] / det;
temp4 = A[0][0] / det;
A[0][0] = temp1;
A[0][1] = temp2;
A[1][0] = temp3;
A[1][1] = temp4;
// M = A^-1 * B
M[0][0] = A[0][0] * B[0][0] + A[0][1] * B[1][0];
M[0][1] = A[0][0] * B[0][1] + A[0][1] * B[1][1];
M[1][0] = A[1][0] * B[0][0] + A[1][1] * B[1][0];
M[1][1] = A[1][0] * B[0][1] + A[1][1] * B[1][1];
// V = t - p*
double V[2] = { t.x - pc.x,t.y - pc.y };
double R[2][1];
// R = M * V
R[0][0] = V[0] * M[0][0] + V[1] * M[1][0];//lv（x）总计算公式
R[1][0] = V[0] * M[0][1] + V[1] * M[1][1];
// R + q*
fv.x = R[0][0] + qc.x;
fv.y = R[1][0] + qc.y;
// Ax = b ==> x = (A^T*A)^-1 * A^T b
return fv;
}
相似变形
注意：是相似变形，和矩阵相似变换 P − 1 A P = B P^{-1}AP=B P−1AP=B不是一个概念。
加权最小二乘形式如下：
∑ w i ∣ ( p ˉ i − p ˉ i ⊥ ) M − ( q ˉ i − q ˉ i ⊥ ) ∣ 2 \sum w_i| \left(\begin{matrix} \bar p_i \\ -\bar p_i^⊥ \end{matrix} \right) M - \left(\begin{matrix} \bar q_i \\ -\bar q_i^⊥ \end{matrix} \right) |^2 ∑wi∣(pˉi−pˉi⊥)M−(qˉi−qˉi⊥)∣2
即为了求解方程组 P M = Q PM=Q PM=Q，类似的有 M = ( P T W P ) − 1 P T W Q M=(P^TWP)^{-1}P^TWQ M=(PTWP)−1PTWQ
故 ( P T W P ) − 1 = 1 ∑ w i p ˉ i p ˉ i T ( 1 0 0 1 ) (P^TWP)^{-1}=\frac{1}{\sum w_i \bar p_i \bar p_i^T}\left(\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix} \right) (PTWP)−1=∑wipˉipˉiT1(1001)
根据该条件得到相似变形的转换矩阵 M M M如下：
M = 1 μ s ∑ w i ( p ˉ i − p ˉ i ⊥ ) ( q ˉ i − q ˉ i ⊥ ) T M = \frac{1}{\mu_s}\sum w_i \left(\begin{matrix} \bar p_i \\ -\bar p_i^⊥ \end{matrix} \right) \left(\begin{matrix} \bar q_i & -\bar q_i^⊥ \end{matrix} \right)^T M=μs1∑wi(pˉi−pˉi⊥)(qˉi−qˉi⊥)T
刚性变形
刚性变形，也就是说变形不含任何缩放效果，我们进一步限制转换矩阵 M M M使其满足 M T M = I M^TM = I MTM=I，这样可以得到刚性变形。
我们知道，矩阵元素每一列的模长为拉伸系数 μ s λ \mu_sλ μsλ，取第一列可得：
μ s λ = ( p ˉ i q ˉ i T ) 2 + ( p ˉ i ( q ˉ i ⊥ ) T ) 2 \mu_sλ =\sqrt{(\bar p_i\bar q_i^T)^2+(\bar p_i(\bar q_i^⊥)^T)^2} μsλ=(pˉiqˉiT)2+(pˉi(qˉi⊥)T)2
故可得，刚性变形的 M M M矩阵解如下
M = 1 λ M 相 似 变 形 = 1 ( ∑ w i p ˉ i q ˉ i T ) 2 + ( ∑ w i p ˉ i ( q ˉ i ⊥ ) T ) 2 ∑ w i ( p ˉ i − p ˉ i ⊥ ) ( q ˉ i − q ˉ i ⊥ ) T M=\frac{1}{λ } M_{相似变形} = \frac{1}{\sqrt{(\sum w_i\bar p_i\bar q_i^T)^2+(\sum w_i\bar p_i(\bar q_i^⊥)^T)^2}}\sum w_i \left(\begin{matrix} \bar p_i \\ -\bar p_i^⊥ \end{matrix} \right) \left(\begin{matrix} \bar q_i & -\bar q_i^⊥ \end{matrix} \right)^T M=λ1M相似变形=(∑wipˉiqˉiT)2+(∑wipˉi(qˉi⊥)T)2 1∑wi(pˉi−pˉi⊥)(qˉi−qˉi⊥)T
参考资料
注：所有参考资料的公式都有误，我重新推导的才是对的。
最终项目：基于移动最小二乘的图像变形
卡通图像变形算法
图像处理（十九）基于移动最小二乘的图像变形