源地址：http://www.cnblogs.com/xiaotie/archive/2009/12/08/1619046.html
有用的地址：http://www.docin.com/p-406100176.html
很多时候，我们需要对一个图像的局部进行调整，这个调整必须是平滑的和可交互式的。Photoshop液化滤镜中向前变形工具就是这样一个工具，很好用。类似工具有美图秀秀(http://www.ddxia.com/view/123008460425144.html)的瘦脸功能。本文描述这类工具背后的原理与算法。
先以美图秀秀为例子，简单描述下向前变形功能。
首先，用鼠标控制一个圆形的选区。
然后，点击鼠标左键，向某个方向拖动，就可以产生光滑的向前变形图片：
通过这个工具，可对图片的局部进行调整，*度比较大，因此比较实用。
下面讲讲这类算法的原理。
上图中，阴影圆环代表一个半径为 rmax 的圆形选区。其中，C点是鼠标点下时的点，也就是圆形选区的圆心。鼠标从C拖到M，致使图像中的点U变换到点X。所以，关键问题是找到上面这个变换的逆变换——给出点X时，可以求出它变换前的坐标U（精确的浮点坐标），然后用变化前图像在U点附近的像素进行插值，求出U的像素值。如此对圆形选区内的每一个像素进行求值，便可得出变换后的图像。
Andreas Gustafsson 的 Interactive Image Warping 一文给出了这一逆变换公式:作者的官网为：http://www.gson.org/
论文下载地址为：http://www.gson.org/thesis/warping-thesis.pdf：
这个变形算法的特点是：
1 只有圆形选区内的图像才进行变形
2 越靠近圆心，变形越大，越靠近边缘的变形越小，边界处无变形
3 变形是平滑的
具体实现步骤如下：
1 对于圆形选区里的每一像素，取出其R,G,B各分量，存入3个Buff（rBuff, gBuff, bBuff）中（也即，三个Buff分别存储选区内的原图像的R,G,B三个通道的数值）
2 对于圆形选区里的每一个像素X，
2.1 根据上面的公式，算出它变形前的位置坐标精确值U
2.2 用插值方法，根据U的位置，和rBuff, gBuff, bBuff中的数值，计算U所在位置处的R,G,B等分量
2.3 将R，G，B等分量合成新的像素，作为X处的像素值
代码我就不贴了，真正对这功能有需求的，根据上面的文字可以很容易写出来——解决这类问题，重要的不是代码，而是思路和算法。
下面是我的实现演示：
上图中，左上角是原图，右下角是变形后的图。红色圆圈圈起来的是变形区域。可以看见，变形很光滑。我在上面的算法中引入了变形强度s(strength)，上图中strength=20。
引入strength，公式就得修改下，下面是我的修改版公式：
看看结果——
原图：
变形，strength=20：
变形，strength=120：
photoshop 与美图秀秀里这个功能可以连续的进行变形。我猜测，这个连续的变形是由一系列基础变形串联起来的，也就是，鼠标从M0拖到Mn位置，并不是只计算 M0->Mn这个变换，而是在鼠标轨迹上引入一系列中间点，M1,M2…Mn-1，然后，对图像进行 M0->M1,M1->M2,…,Mn-1->Mn等一系列变换。
下面jia_zhengshen原创。。。。并非原作者（小铁）所写。
这仅是我的猜想，我努力的实现了一下串联起来然后实现Photoshop的效果，但是我没有实现。并且这个算法有个天然的弊病只能在圆内进行变形，圆外由于上面提到的公式的特性相当陡，所以我认为并不是很实用。
下面贴一下我自己实现的代码，是使用opencv实现的。如果谁能实现Photoshop的变形，请告诉我哦。代码很乱，还望见谅。
#include<opencv\cv.h>
#include<opencv\highgui.h>
typedef struct warp
{
double origCenterX;
double origCenterY;
double nowCenterX;
double nowCenterY;
double radius;
bool isFirst;
IplImage *img;
IplImage *dst;
IplImage *temp;
IplImage *mask;
IplImage *temp2;
} _warp;
///warp *w ;
template<class T>
T biLine(double f00,double f01,double f10,double f11,double u,double v)//双线性插值算法。
{
return (T)(1-u)*(1-v)*f00+(1-u)*v*f01+u*(1-v)*f10 + u*v*f11;
}
int mapping(warp *w ,int x,int y, double *u,double *v)
{
double rmax2 = w->radius * w->radius;
double x_c2 = (x-w->origCenterX)*(x-w->origCenterX) + (y-w->origCenterY) *(y-w->origCenterY);
double m_c2 = (w->nowCenterX-w->origCenterX)*(w->nowCenterX-w->origCenterX) \
+(w->nowCenterY-w->origCenterY)*(w->nowCenterY-w->origCenterY);
double t =rmax2-x_c2+m_c2;
double f2;
if(t==0)
f2 = 0;
else
f2= ((rmax2-x_c2)/(rmax2-x_c2+m_c2))*((rmax2-x_c2)/(rmax2-x_c2+m_c2));
double m_cX = w->nowCenterX - w->origCenterX;
double m_cY = w->nowCenterY - w->origCenterY;
*u = x-f2*m_cX;
*v = y-f2*m_cY;
if(*u<0)
*u = 0;
if(*u>w->img->width)
*u = w->img->width-1;
if(*v<0)
*v = 0;
if(*v>w->img->height)
*v = w->img->height-1;
return 0;
}
int mapping2(warp *w ,int x,int y, double *u,double *v)
{
double fu,fv;
fu = x;
fv = y;
double dx = fu - w->origCenterX;
double dy = fv - w->origCenterY;
double rsq = dx*dx+dy*dy;
if(dx>0.0-w->radius
&& dx<w->radius
&& dy>0.0-w->radius
&& dy<w->radius
&& rsq< w->radius*w->radius
)
/*if(1)*/
{
double cmx = w->nowCenterX - w->origCenterX;
double cmy = w->nowCenterY - w->origCenterY;
//double msq = (dx-cmx)*(dx-cmx) +(dy-cmy)*(dy-cmy);
double msq = cmx*cmx + cmy*cmy;
double edge_dist = w->radius*w->radius -rsq;
double a = edge_dist/(edge_dist+msq);
a *= a;
fu -= a*cmx;
fv -= a*cmy;
*u = fu;
*v = fv;
}else
{
*u = x;
*v = y;
}
return 0;
}
int imProc(warp * w)
{
IplImage *img = w->img;
IplImage *dst = w->dst;
IplImage *temp = w->temp;
IplImage *mask = w->mask;
for(int i=0;i<img->height;i++)
{
uchar *ptr =(uchar*)(img->imageData+i*img->widthStep);
uchar *ptrWarp = (uchar *) (dst->imageData+i*dst->widthStep);
uchar *ptrTemp = (uchar *) (temp->imageData+i*temp->widthStep);
//uchar *ptrMask = (uchar *) (mask->imageData+i*mask->widthStep);
for(int j=0;j<img->width;j++)
{
double dx = j-w->nowCenterX;
double dy = i-w->nowCenterY;
if(dx> (0.0-w->radius)
&& dx<(w->radius)
&& dy>(0.0-w->radius)
&& dy<w->radius
&& (dx*dx+dy*dy) <= w->radius *w->radius
)
{
double u,v;
mapping2(w,j,i,&u,&v);
if(u>img->width)
u = img->width;
if(v>img->height)
v = img->height ;
uchar *ptr4 = (uchar *)(img->imageData+(int)v*img->widthStep);
//uchar *ptrMask = (uchar *)(mask->imageData+(int)v*img->widthStep);
//if(*(ptrMask+j) !
int iu= (int)u;
int iv = (int)v;
ptrWarp[3*j+0] = ptr4[3*iu+0];
ptrWarp[3*j+1] = ptr4[3*iu+1];
ptrWarp[3*j+2] = ptr4[3*iu+2];
ptrTemp[3*j+0] = ptr4[3*iu+0];
ptrTemp[3*j+1] = ptr4[3*iu+1];
ptrTemp[3*j+2] = ptr4[3*iu+2];
//double f00,f01,f10,f11;
double f00b,f00g,f00r;
double f01b,f01g,f01r;
double f10b,f10g,f10r;
double f11b,f11g,f11r;
int x00 = (int)u;
int y00 = (int)v;
uchar *ptr1 = (uchar *)(img->imageData+y00*img->widthStep);
f00b = ptr1[3*x00+0];
f00g = ptr1[3*x00+1];
f00r = ptr1[3*x00+2];
int x01 = (int)(u+0.5);
int y01 = (int )v;
uchar *ptr2 = (uchar *)(img->imageData+y01*img->widthStep);
f01b = ptr2[3*x01+0];
f01g = ptr2[3*x01+1];
f01r = ptr2[3*x01+2];
int x10 = (int)u;
int y10 = (int)(v+0.5);
uchar *ptr3 = (uchar *)(img->imageData+y10*img->widthStep);
f10b = ptr3[3*x10+0];
f10g = ptr3[3*x10+1];
f10r = ptr3[3*x10+2];
int x11 = (int)(u+0.5);
int y11 = (int)(v+0.5);
ptr4 = (uchar *)(img->imageData+y11*img->widthStep);
f11b = ptr4[3*x11+0];
f11g = ptr4[3*x11+1];
f11r = ptr4[3*x11+2];
//ptr4 = (uchar *)(img->imageData+(int)v*img->widthStep);
//int iu= (int)u;
//int iv = (int)v;
//ptrWarp[3*j+0] = ptr4[3*iu+0];
//ptrWarp[3*j+1] = ptr4[3*iu+1];
//ptrWarp[3*j+2] = ptr4[3*iu+2];
int a = biLine<int>(f00b,f01b,f10b,f11b,u-(int)u,v-(int)v);
//ptrWarp[3*j+0] = biLine<int>(f00b,f01b,f10b,f11b,u-(int)u,v-(int)v);
//ptrWarp[3*j+1] = biLine<int>(f00g,f01g,f10g,f11g,u-(int)u,v-(int)v);
//ptrWarp[3*j+2] = biLine<int>(f00r,f01r,f10r,f11r,u-(int)u,v-(int)v);
}//if
else
{
ptrWarp[3*j+0] = ptrTemp[3*j+0];
ptrWarp[3*j+1] = ptrTemp[3*j+1];
ptrWarp[3*j+2] = ptrTemp[3*j+2];
}
}//width
}//height
cvShowImage("ttt",dst);
//vWaitKey();
return 0;
}
int imProc2(warp * w)
{
IplImage *img = w->img;
IplImage *dst = w->dst;
IplImage *temp = w->temp;
IplImage *mask = w->mask;
for(int i=0;i<img->height;i++)
{
uchar *ptr =(uchar*)(img->imageData+i*img->widthStep);
uchar *ptrWarp = (uchar *) (dst->imageData+i*dst->widthStep);
uchar *ptrTemp = (uchar *) (temp->imageData+i*temp->widthStep);
uchar *ptrMask = (uchar *) (mask->imageData+i*mask->widthStep);
for(int j=0;j<img->width;j++)
{
double dx = j-w->nowCenterX;
double dy = i-w->nowCenterY;
if(dx> (0.0-w->radius)
&& dx<(w->radius)
&& dy>(0.0-w->radius)
&& dy<w->radius
&& (dx*dx+dy*dy) < w->radius *w->radius
)
/*if(1)*/
{
double u,v;
mapping2(w,j,i,&u,&v);
if(u>img->width)
u = img->width;
if(v>img->height)
v = img->height ;
if(u<0)
u =0;
if(v<0)
v = 0;
uchar *ptrMask = (uchar *)(mask->imageData+(int)v*mask->widthStep);
//if(*(ptrMask+j) ==0)
if(1)
{
//uchar *ptr4 = (uchar *)(img->imageData+(int)v*img->widthStep);
//uchar *ptr4 = (uchar *)(dst->imageData+(int)v*dst->widthStep);
uchar *ptr4 = (uchar *) (w->img->imageData+(int)v*img->widthStep);
int iu= (int)u;
int iv = (int)v;
ptrWarp[3*j+0] = ptr4[3*iu+0];
ptrWarp[3*j+1] = ptr4[3*iu+1];
ptrWarp[3*j+2] = ptr4[3*iu+2];
ptrTemp[3*j+0] = ptr4[3*iu+0];
ptrTemp[3*j+1] = ptr4[3*iu+1];
ptrTemp[3*j+2] = ptr4[3*iu+2];
//*(ptrMask+j) =1;
}else
{
uchar *ptr4 = (uchar *)(img->imageData+(int)v*img->widthStep);
int iu= (int)u;
int iv = (int)v;
ptrWarp[3*j+0] = ptr4[3*iu+0];
ptrWarp[3*j+1] = ptr4[3*iu+1];
ptrWarp[3*j+2] = ptr4[3*iu+2];
ptrTemp[3*j+0] = ptr4[3*iu+0];
ptrTemp[3*j+1] = ptr4[3*iu+1];
ptrTemp[3*j+2] = ptr4[3*iu+2];
//*(ptrMask+j) =1;
}
}//if
else
{
/*if( *(ptrMask+j)==0)
{
ptrWarp[3*j+0] = ptr[3*j+0];
ptrWarp[3*j+1] = ptr[3*j+1];
ptrWarp[3*j+2] = ptr[3*j+2];
}else*/
{
ptrWarp[3*j+0] = ptr[3*j+0];
ptrWarp[3*j+1] = ptr[3*j+1];
ptrWarp[3*j+2] = ptr[3*j+2];
}
}
}//width
}//height
//cvShowImage("ttt",dst);
//vWaitKey();
/*IplImage *ttt ;
ttt = w->temp;
w->temp = w->temp2;
w->temp2 = ttt;*/
cvCopy(dst,w->temp);
return 0;
}
void myMouseCallback(int event ,int x,int y,int flags,void *param)
{
warp *w = (warp*)param;
if(x<0 ||x>=w->img->width)
return ;
if(y<0 ||y>=w->img->height)
return ;
switch(event)
{
case CV_EVENT_MOUSEMOVE:
if(w->isFirst ==false)//确定了原来的圆心。
{
//w->nowCenterX = x;
//w->nowCenterY = y;
////std::cout<<w->nowCenterX<<std::endl;
//imProc2(w);
////;
////;//适当的处理
////if(fabs((double)w->origCenterX-x)> w->radius/5 )
//{
// w->origCenterX = x;
// w->origCenterY = y;
//}
}else//没有确定圆心。
{
;
}
break;
case CV_EVENT_LBUTTONDOWN:
w->origCenterX = x;
w->origCenterY = y;
w->isFirst= false;
break;
case CV_EVENT_LBUTTONUP:
w->isFirst = true;
w->nowCenterX = x;
w->nowCenterY = y;
//std::cout<<w->nowCenterX<<std::endl;
imProc2(w);
break;
}
//cvShowImage("dst",w->img);
}
int main()
{
IplImage *img = cvLoadImage("hello2.jpg");
IplImage *dst = cvCreateImage(cvGetSize(img),8,3);
IplImage *temp = cvCreateImage(cvGetSize(img),8,3);
IplImage *temp2 = cvCreateImage(cvGetSize(img),8,3);
IplImage *mask = cvCreateImage(cvGetSize(img),8,1);
cvZero(mask);
cvCopy(img,dst,NULL);
cvCopy(img,temp,NULL);
cvCopy(img,temp2,NULL);
warp *w = new warp;
w->img = img;
w->mask = mask;
w->dst = dst;
w->temp = temp;
w->temp2 = temp2;
w->isFirst = true;
w->radius = 100;
cvNamedWindow("hello");
cvSetMouseCallback("hello",myMouseCallback,w);
while(1)
{
int key = cvWaitKey(2);
switch(key)
{
case 'q':
case 27:
goto hear;
break;
case '1':
case '2':
case '3':
case '4':
case '5':
case '6':
case '7':
case '8':
case '9':
w->radius = (key-0x30)*10;
std::cout<<w->radius<<std::endl;
break;
case 's':
cvSaveImage("蛋蛋.jpg",w->dst);
}
cvShowImage("hello",w->dst);
//cvShowImage("dst",dst);
}
hear:
cvReleaseImage(&img);
return 0;
}