图像变形算法之meitu瘦脸及放大眼睛（文末含代码）
导语：很多时候，我们需要对一个图像的局部进行调整，这个调整必须是平滑的和可交互式的。Photoshop液化滤镜中向前变形工具就是这样一个工具，很好用。类似工具有美图秀秀(http://xiuxiu.meitu.com/)的瘦脸功能。本文描述这类工具背后的原理与算法。
先以美图秀秀为例子，简单描述下向前变形功能。
首先，用鼠标控制一个圆形的选区。
然后，点击鼠标左键，向某个方向拖动，就可以产生光滑的向前变形图片：
通过这个工具，可对图片的局部进行调整，*度比较大，因此比较实用。
主要参考：交互式图像变形算法一文 。看过该文章后，自然了解瘦脸其实用的是图像局部平移变形，放大眼睛，用的是图像局部缩放变形。同时，文章还阐述了图像局部旋转变形的算法原理及思路。
瘦脸及放大眼睛的前提是需要检测到人脸，并提取特征点。谈到图像变形，最基础的思路是：由变形前坐标，根据变形映射关系，得到变形后坐标。这其中变形映射关系是最关键的，不同的映射关系，将得到不同的变形效果。平移、缩放、旋转，对应的是不同的映射关系，即不同的变换公式。当然实际在计算过程中，用的是逆变换，即由变形后坐标，根据逆变换公式反算变形前坐标，然后插值得到该坐标rgb像素值，将该rgb值作为变形后坐标对应的像素值。这样才能保证变形后的图像是连续、完整的。
下面讲讲这类算法的原理。
上图中，阴影圆环代表一个半径为 rmax 的圆形选区。其中，C点是鼠标点下时的点，也就是圆形选区的圆心。鼠标从C拖到M，致使图像中的点U变换到点X。所以，关键问题是找到上面这个变换的逆变换——给出点X时，可以求出它变换前的坐标U（精确的浮点坐标），然后用变化前图像在U点附近的像素进行插值，求出U的像素值。如此对圆形选区内的每一个像素进行求值，便可得出变换后的图像。
Andreas Gustafsson 的 Interactive Image Warping 一文给出了这一逆变换公式：
公式中，由于主要是像素点位置计算，因此涉及一些矢量运算，不过比较简单。其实上面公式就是逆变换公式了，x是变换后的位置，u是原坐标位置。整个计算在以c为圆心，r为半径的圆内进行。因为是交互式图像局部变形，所以c也可以看做鼠标点下时的坐标，而m为鼠标移动一段距离后抬起时的坐标，这样c和m就决定了变形方向。下面是示例代码，公式结合代码一起看，应该很快能弄明白。
这个变形算法的特点是：
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等一系列变换。
函数功能实现
void LocalTranslationWarp(Mat &img, int warpX, int warpY, int warpW, int warpH, int directionX, int directionY, double warpCoef)
{
RestrictBounds(warpX, warpY, warpW, warpH);
Mat imgCopy;
copyMakeBorder(img, imgCopy, 0, 1, 0, 1, BORDER_REPLICATE);
Point center(warpX + (warpW>>1), warpY + (warpH>>1));
double radius = (warpW < warpH) ? (warpW >> 1) : (warpH >> 1);
radius = radius * radius;
// 平移方向矢量/模
double transVecX = directionX - center.x;
double transVecY = directionY - center.y;
double transVecModel = transVecX*transVecX + transVecY*transVecY;
// 水平/垂直增量//映射后位置与原位置
double dx = 0, dy = 0, posX = 0.0, posY = 0.0, posU = 0.0, posV = 0.0;
// 点到圆心距离/平移比例
double distance = 0.0, ratio = 0.0;
// 插值位置
int startU = 0, startV = 0;
double alpha = 0.0, beta = 0.0;
int maxRow = warpY + warpH;
int maxCol = warpX + warpW;
uchar* pImg = NULL;
for (int i = warpY; i < maxRow; i++)
{
pImg = img.data + img.step * i;
for (int j = warpX; j < maxCol; j++)
{
posX = j;
posY = i;
dx = posX - center.x;
dy = posY - center.y;
distance = dx*dx + dy*dy;
if (distance < radius)
{
ratio = (radius - distance) / (radius - distance + transVecModel * warpCoef);
posU = posX - ratio * ratio * transVecX;
posV = posY - ratio * ratio * transVecY;
startU = (int)posU;
startV = (int)posV;
alpha = posU - startU;
beta = posV - startV;
BilinearInter(imgCopy, startU, startV, alpha, beta, pImg[3*j], pImg[3*j + 1], pImg[3*j + 2]);
}
}
}
}