《Image Deformation Using Moving Least Squares》-- 阅读笔记

参考资料：
https://zhuanlan.zhihu.com/p/103212151
https://blog.csdn.net/hjimce/article/details/46550001

一.背景

《Surface Defect Saliency of Magnetic Tile》一文中提到过一种基于移动最小二乘图像变形的数据增强方法

该方法可以通过操纵少数点完成网格的变形，实现图像变形，可以生成更多的形状不同的缺陷图像，增大数据量。

二. 详细步骤

三.算法实现

具体讲仿射变换的变形算法实现。

**问题：**原图像的各个控制顶点坐标p，原图像上的像素点v的坐标。变形后图像的控制顶点位置q，求v在变形后图像中对应位置f(v)。

总计算公式：

上面中lv(x)和f(v)是同一个函数。因为x就是我们输入的原图像的像素点坐标v。

因此我们的目标就是要知道p*，q*，变换矩阵M。这样输入一个参数v，我们就可以计算出它在变换后图像中的位置了。

下面的工作就是求出3式中的参数

1. 权重w的计算方法

也就是计算v到控制顶点pi的距离倒数作为权重，参数a一般取值为1。

//计算各个控制顶点的权重，也就是计算点t到各个顶点的距离1/sqr(d)
while(iter!=p.end())
{
	double temp;
	if(iter->x!=t.x || iter->y!=t.y)
	    temp=1/((iter->x-t.x)*(iter->x-t.x)+(iter->y-t.y)*(iter->y-t.y));
	else//如果t为控制顶点，那么需要把该控制顶点的权重设置为无穷大
		temp=MAXNUM;
	w.push_back(temp);
	iter++;
	}

2. q*,p*的计算公式

也就是计算控制顶点pi和qi的加权求和重心位置。

double px=0,py=0,qx=0,qy=0,tw=0;
while(iterw!=w.end())
{
px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
py+=(*iterw)*(iter->y);
qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
qy+=(*iterw)*(iterq->y);
tw+=*iterw;//总权重
iter++;
iterw++;
iterq++;
}
pc.x=px/tw;
pc.y=py/tw;
qc.x=qx/tw;
qc.y=qy/tw;

3. 仿射变换1矩阵M的计算公式

四. MLS 源代码

//输入原图像的t点，输出变形后图像的映射点f(v)
MyPoint CMLSDlg::MLS(const MyPoint& t)
{
	if(p.empty())//原图像的控制顶点p，与输入点t为同一副图像坐标系下
		return t;
	MyPoint fv;
	double A[2][2],B[2][2],M[2][2];
	iter=p.begin();
	w.erase(w.begin(),w.end());//初始化权重数据
	//计算各个控制顶点的权重，也就是计算点t到各个顶点的距离1/sqr(d)
	while(iter!=p.end())
	{
		double temp;
		if(iter->x!=t.x || iter->y!=t.y)
		    temp=1/((iter->x-t.x)*(iter->x-t.x)+(iter->y-t.y)*(iter->y-t.y));
		else//如果t为控制顶点，那么需要把该控制顶点的权重设置为无穷大
			temp=MAXNUM;
		w.push_back(temp);
		iter++;
	}
	vector<double>::iterator iterw=w.begin();
	vector<MyPoint>::iterator iterq=q.begin();//q为目标图像的控制点的位置，我们的目标是找到t在q中的对应位置
	iter=p.begin();
	//求p*，q*
	MyPoint pc,qc;
	double px=0,py=0,qx=0,qy=0,tw=0;
	while(iterw!=w.end())
	{
	px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
	py+=(*iterw)*(iter->y);
	qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
	qy+=(*iterw)*(iterq->y);
	tw+=*iterw;//总权重
	iter++;
    iterw++;
	iterq++;
	}
	pc.x=px/tw;
	pc.y=py/tw;
	qc.x=qx/tw;
	qc.y=qy/tw;
	iter=p.begin();
	iterw=w.begin();
	iterq=q.begin();
	//初试化
	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;
	}
	//计算pi-p*,qi-q*
	while(iter!=p.end())
	{
 
	double P[2]={iter->x-pc.x,iter->y-pc.y};
	double PT[2][1];
	PT[0][0]=iter->x-pc.x;
	PT[1][0]=iter->y-pc.y;
	double Q[2]={iterq->x-qc.x,iterq->y-qc.y};
	double T[2][2];
   
   //定义T储存piTpi相乘结果
    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];
 	//计算pjqj的类加
	for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
	{
	//在T乘上权重，赋给A数组
	A[i][j]+=(*iterw)*T[i][j];
	}
	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];
    
    
	for(int i=0;i<2;i++)
    for(int j=0;j<2;j++)
	{
	B[i][j]+=(*iterw)*T[i][j];
	}
 
	iter++;
    iterw++;
	iterq++;
	}
	//cvInvert(A,M);
    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
	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];
 	//算出x-p*
	double V[2]={t.x-pc.x,t.y-pc.y};
	double R[2][1];
  	//算出(x-p*)M
	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];
	//最后结果fv=(x-p*)M+q*
	fv.x=R[0][0]+qc.x;
	fv.y=R[1][0]+qc.y;
 
	return fv;
}

总结

图像变形算法，有正向映射和逆向映射，如果按照每个像素点，都通过上面的计算方法求取其对应变换后的像素点位置，那么其实计算量是非常大的，因为一幅图像的像素点，实在是太多了，如果每个像素点，都用上面的函数遍历过一遍，那计算量可想而知。
因此一般的变形算法是对待图像进行三角剖分：

然后只根据只对三角网格模型的顶点，根据变形算法，计算出三角网格模型每个顶点的新位置，最后再用三角形仿射变换的方法，计算三角形内每个像素点的值，得到变形后的图像，这样不仅速度快，同事解决了正向映射与逆向映射变形算法存在的不足之处，具体图像变形的正向和逆向映射存在的缺陷，可以自己查看相关的文献。