影像處理（十九）基於移動最小二乘的映像變形-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。

這一步實現代碼如下：

[cpp]  view plain  copy //計算各個控制頂點的權重，也就是計算點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的加權求和重心位置。

[cpp]  view plain  copy 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、仿射變換矩陣M的計算公式如下:

只要把相關的參數都帶進去就可以計算了。

最後貼一些完整的MLS原始碼：

[cpp]  view plain  copy //輸入原映像的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();       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;       }       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[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];          for(int i=0;i<2;i++)       for(int j=0;j<2;j++)       {       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[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];          double V[2]={t.x-pc.x,t.y-pc.y};       double R[2][1];            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=R[0][0]+qc.x;       fv.y=R[1][0]+qc.y;          return fv;   }  

調用方法樣本：

[cpp]  view plain  copy    int i=0,j=0;   dImage=cvCreateImage(cvSize(2*pImage->width,2*pImage->height),pImage->depth,pImage->nChannels);//建立新的變形映像   cvSet(dImage,cvScalar(0));   MyPoint Orig=MLS(MyPoint(IR_X,IR_Y));   int Orig_x=(int)(Orig.x)-(int)(pImage->width/2);   int Orig_y=(int)(Orig.y)-(int)(pImage->height/2);      for(i=0;i<pImage->height;i++)//遍曆原映像的每個像素   {       for(j=0;j<pImage->width;j++)       {           CvScalar color;           double x=j+IR_X;           double y=i+IR_Y;           MyPoint t=MLS(MyPoint(x,y));//MLS計算原映像(x,y)在靶心圖表像的映射位置f(v)           int m=(int)(t.x);           int n=(int)(t.y);           m-=Orig_x;              n-=Orig_y;           color=cvGet2D(pImage,i,j);//像素擷取           if(0<=m && dImage->width>m && 0<=n && dImage->height>n)           {           cvSet2D(dImage,n,m,color);           }       }   }   映像變形演算法，有正向映射和逆向映射，如果按照每個像素點，都通過上面的計算方法求取其對應變換後的像素點位置，那麼其實計算量是非常大的，因為一幅映像的像素點，實在是太多了，如果每個像素點，都用上面的函數遍曆過一遍，那計算量可想而知。

因此一般的變形演算法是對待映像進行三角剖分：

然後只根據只對三角網格模型的頂點，根據變形演算法，計算出三角網格模型每個頂點的新位置，最後再用三角形仿射變換的方法，計算三角形內每個像素點的值，得到變形後的映像，這樣不僅速度快，同事解決了正向映射與逆向映射變形演算法存在的不足之處，具體映像變形的正向和逆向映射存在的缺陷，可以自己查看相關的文獻。

另外兩種相似變換和剛性變換，可以自己查看M矩陣的計算公式，編寫實現相關代碼。

本文地址：http://blog.csdn.net/hjimce/article/details/46550001     作者：hjimce     聯絡qq：1393852684   更多資源請關注我的部落格：http://blog.csdn.net/hjimce                原創文章，著作權，轉載請保留本行資訊。

參考文獻：

1、《Image Deformation Using Moving Least Squares》

2、《3D Deformation Using Moving Least Squares》

from: http://blog.csdn.net/hjimce/article/details/46550001