计算机视觉
图像处理

图像处理(十三)基于移动最小二乘的图像变形-Siggraph 2006

一、背景意义

写这篇博文是应为目前为止我看到了好多领域里的经典paper算法都有涉及到移动最小二乘(MLS)。可见这个算法非常重要,先来看一下它的相关经典应用:

1、图像变形。在图像处理领域paper:《Image Deformation Using Moving Least Squares》利用移动最小二乘的原理实现了图像的相关变形,而且这篇paper的引用率非常高,可以说是图像变形算法的经典算法,Siggraph上 面的paper。

利用移动最小二乘实现图像变形

2、点云滤波。利用MLS实现点云滤波,是三维图像学点云处理领域的一大应用,我所知道点云滤波经典算法包括:双边滤波、MLS、WLOP。

3、Mesh Deformation。用这个算法实现三角网格模型的变形应用也是非常不错的,相关的paper《3D Deformation Using Moving Least Squares

OK,接着我就以《Image Deformation Using Moving Least Squares》算法为例,进行讲解基于移动最小二乘的图像变形算法实现。

二、算法实现

在这里我没有打算将算法原理的推导过程,直接讲算法的实现步骤公式。

这篇paper根据变换矩阵的不同,可以分为三种变形方法,分别是仿射变换、相似变换、刚性变换。其中刚性变换的效果是最好的,我这边从简单的讲,只讲仿射变换的变形算法实现:

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

总计算公式为:

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

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

OK,只要我们知道上面公式中,各个参数的计算方法,我们就可以计算出变形后图像对应的坐标点f(v)了。

1、权重w的计算方法为:

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

这一步实现代码如下:

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

2、q*,p*的计算公式如下:

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

  1. double px=0,py=0,qx=0,qy=0,tw=0;
  2. while(iterw!=w.end())
  3. {
  4. px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
  5. py+=(*iterw)*(iter->y);
  6. qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
  7. qy+=(*iterw)*(iterq->y);
  8. tw+=*iterw;//总权重
  9. iter++;
  10.    iterw++;
  11. iterq++;
  12. }
  13. pc.x=px/tw;
  14. pc.y=py/tw;
  15. qc.x=qx/tw;
  16. qc.y=qy/tw;

3、仿射变换矩阵M的计算公式如下:

只要把相关的参数都带进去就可以计算了。

最后贴一些完整的MLS源代码:

  1. //输入原图像的t点,输出变形后图像的映射点f(v)
  2. MyPoint CMLSDlg::MLS(const MyPoint& t)
  3. {
  4.     if(p.empty())//原图像的控制顶点p,与输入点t为同一副图像坐标系下
  5.         return t;
  6.     MyPoint fv;
  7.     double A[2][2],B[2][2],M[2][2];
  8.     iter=p.begin();
  9.     w.erase(w.begin(),w.end());
  10.     //计算各个控制顶点的权重,也就是计算点t到各个顶点的距离1/sqr(d)
  11.     while(iter!=p.end())
  12.     {
  13.         double temp;
  14.         if(iter->x!=t.x || iter->y!=t.y)
  15.             temp=1/((iter->x-t.x)*(iter->x-t.x)+(iter->y-t.y)*(iter->y-t.y));
  16.         else//如果t为控制顶点,那么需要把该控制顶点的权重设置为无穷大
  17.             temp=MAXNUM;
  18.         w.push_back(temp);
  19.         iter++;
  20.     }
  21.     vector<double>::iterator iterw=w.begin();
  22.     vector<MyPoint>::iterator iterq=q.begin();//q为目标图像的控制点的位置,我们的目标是找到t在q中的对应位置
  23.     iter=p.begin();
  24.     MyPoint pc,qc;
  25.     double px=0,py=0,qx=0,qy=0,tw=0;
  26.     while(iterw!=w.end())
  27.     {
  28.     px+=(*iterw)*(iter->x);//所有控制顶点p的加权位置
  29.     py+=(*iterw)*(iter->y);
  30.     qx+=(*iterw)*(iterq->x);//所有控制顶点q的加权位置
  31.     qy+=(*iterw)*(iterq->y);
  32.     tw+=*iterw;//总权重
  33.     iter++;
  34.     iterw++;
  35.     iterq++;
  36.     }
  37.     pc.x=px/tw;
  38.     pc.y=py/tw;
  39.     qc.x=qx/tw;
  40.     qc.y=qy/tw;
  41.     iter=p.begin();
  42.     iterw=w.begin();
  43.     iterq=q.begin();
  44.     for(int i=0;i<2;i++)
  45.     for(int j=0;j<2;j++)
  46.     {
  47.     A[i][j]=0;
  48.     B[i][j]=0;
  49.     M[i][j]=0;
  50.     }
  51.     while(iter!=p.end())
  52.     {
  53.     double P[2]={iter->x-pc.x,iter->y-pc.y};
  54.     double PT[2][1];
  55.     PT[0][0]=iter->x-pc.x;
  56.     PT[1][0]=iter->y-pc.y;
  57.     double Q[2]={iterq->x-qc.x,iterq->y-qc.y};
  58.     double T[2][2];
  59.     T[0][0]=PT[0][0]*P[0];
  60.     T[0][1]=PT[0][0]*P[1];
  61.     T[1][0]=PT[1][0]*P[0];
  62.     T[1][1]=PT[1][0]*P[1];
  63.     for(int i=0;i<2;i++)
  64.     for(int j=0;j<2;j++)
  65.     {
  66.     A[i][j]+=(*iterw)*T[i][j];
  67.     }
  68.     T[0][0]=PT[0][0]*Q[0];
  69.     T[0][1]=PT[0][0]*Q[1];
  70.     T[1][0]=PT[1][0]*Q[0];
  71.     T[1][1]=PT[1][0]*Q[1];
  72.     for(int i=0;i<2;i++)
  73.     for(int j=0;j<2;j++)
  74.     {
  75.     B[i][j]+=(*iterw)*T[i][j];
  76.     }
  77.     iter++;
  78.     iterw++;
  79.     iterq++;
  80.     }
  81.     //cvInvert(A,M);
  82.     double det=A[0][0]*A[1][1]-A[0][1]*A[1][0];
  83.     if(det<0.0000001)
  84.     {
  85.         fv.x=t.x+qc.x-pc.x;
  86.         fv.y=t.y+qc.y-pc.y;
  87.         return fv;
  88.     }
  89.     double temp1,temp2,temp3,temp4;
  90.     temp1=A[1][1]/det;
  91.     temp2=-A[0][1]/det;
  92.     temp3=-A[1][0]/det;
  93.     temp4=A[0][0]/det;
  94.     A[0][0]=temp1;
  95.     A[0][1]=temp2;
  96.     A[1][0]=temp3;
  97.     A[1][1]=temp4;
  98.     M[0][0]=A[0][0]*B[0][0]+A[0][1]*B[1][0];
  99.     M[0][1]=A[0][0]*B[0][1]+A[0][1]*B[1][1];
  100.     M[1][0]=A[1][0]*B[0][0]+A[1][1]*B[1][0];
  101.     M[1][1]=A[1][0]*B[0][1]+A[1][1]*B[1][1];
  102.     double V[2]={t.x-pc.x,t.y-pc.y};
  103.     double R[2][1];
  104.     R[0][0]=V[0]*M[0][0]+V[1]*M[1][0];//lv(x)总计算公式
  105.     R[1][0]=V[0]*M[0][1]+V[1]*M[1][1];
  106.     fv.x=R[0][0]+qc.x;
  107.     fv.y=R[1][0]+qc.y;
  108.     return fv;
  109. }

调用方法示例:

  1.    int i=0,j=0;
  2. dImage=cvCreateImage(cvSize(2*pImage->width,2*pImage->height),pImage->depth,pImage->nChannels);//创建新的变形图像
  3. cvSet(dImage,cvScalar(0));
  4. MyPoint Orig=MLS(MyPoint(IR_X,IR_Y));
  5. int Orig_x=(int)(Orig.x)-(int)(pImage->width/2);
  6. int Orig_y=(int)(Orig.y)-(int)(pImage->height/2);
  7. for(i=0;i<pImage->height;i++)//遍历原图像的每个像素
  8. {
  9.     for(j=0;j<pImage->width;j++)
  10.     {
  11.         CvScalar color;
  12.         double x=j+IR_X;
  13.         double y=i+IR_Y;
  14.         MyPoint t=MLS(MyPoint(x,y));//MLS计算原图像(x,y)在目标图像的映射位置f(v)
  15.         int m=(int)(t.x);
  16.         int n=(int)(t.y);
  17.         m-=Orig_x;
  18.            n-=Orig_y;
  19.         color=cvGet2D(pImage,i,j);//像素获取
  20.         if(0<=m && dImage->width>m && 0<=n && dImage->height>n)
  21.         {
  22.         cvSet2D(dImage,n,m,color);
  23.         }
  24.     }
  25. }

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

因此一般的变形算法是对待图像进行三角剖分:

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

另外两种相似变换和刚性变换,可以自己查看M矩阵的计算公式,编写实现相关代码。

参考文献:

1、《Image Deformation Using Moving Least Squares》

2、《3D Deformation Using Moving Least Squares》

转载注明来源:CV视觉网 » 图像处理(十三)基于移动最小二乘的图像变形-Siggraph 2006

分享到:更多 ()
扫描二维码,给作者 打赏
pay_weixinpay_weixin

请选择你看完该文章的感受:

0不错 0超赞 0无聊 0扯淡 0不解 0路过

评论 4

评论前必须登录!