CV学习-Morphing

Morphing - 人脸渐变

概述

影像变形（Morphing），是由一张图像流畅地变成另一张图像的视觉效果，最常见的应用是由一张人脸影像变化到另一张人脸影像。
影像变形最基本的做法是借由交叉融合（cross-dissolving）的方法来做，也就是把两张影像在变化时间上作
线性的内插。这样的方法固然简单，但由于两张影像形状不同（例如两个人脸型不同)，会产生鬼影的现象。为了解决鬼影的现象，必须在交叉融合之前进行改变影像形状的弯曲转换（warping）这个动作。因此，整个影像变形的过程，就是翘曲加上交叉融合。

具体效果如下：

过程

Warping是由一张影像渐进到另一张影像的效果，作法是计算在影像变化时间内所有影格的中间影像，最后再把这一连串的影像接续拨放就达到变形的效果，以下步骤为对于某一个影格的中间影像的产生方法：

• 根据使用的特征指定方法（网格、点、线段），将两张图片的对应特征作线性内插，得到中间影像的特征资讯
• 把两张影像都弯曲转换到中间影像
• 把这两张弯曲转换过后的影像作交叉融合，也就是把这两张影像对应的颜色取平均，如此便可以得到在这个影格的中间影像

具体实现（by C++）

测试数据

第三方库准备

1. CImg库用来对图片进行读写以及像素级别的操作。Cimg下载
   在使用的地方添加头文件#include “CImg.h”，引用命名空间using namespace cimg_library
2. Dlib库用来获取标准68个人脸特征点。Dlib下载

实现步骤

基本的数据结构

点

#ifndef POINT_H#define POINT_H// 3D vectorclass Point {public:    Point() {}    Point(double fx, double fy, double fz)        :x(fx), y(fy), z(fz) {    }    Point operator - (const Point& v) const {        return Point(x - v.x, y - v.y, z - v.z);    }    // Dot product    double Dot(const Point& v) const {        return x * v.x + y * v.y + z * v.z;    }    // Cross product    Point Cross(const Point& v) const {        return Point(            y * v.z - z * v.y,            z * v.x - x * v.z,            x * v.y - y * v.x);    }

边、三角形以及他们的容器

typedef std::vector<Point> PointArray;//定义点类的vector容器  typedef struct  {      int left;      int right;      int count;//边的计数，如果计数为0，则删除此边  }Edge;//定义边类  typedef std::vector<Edge> EdgeArray;//定义边类的vector容器  typedef struct  {      int v[3];//三角形的三个顶点      Edge s[3];//三角形的三条边      double xc;//三角形外接圆圆心的x坐标      double yc;//三角形外接圆圆心的y坐标      double r;//三角形外接圆的半径    double** m1; //到原图的映射矩阵    double** m2; //到目标图的映射矩阵}Triangle;//定义三角形类  typedef std::vector<int> intArray;//定义int类的vector容器

获取标准68个人脸特征点

如果在VS中使用Dlib以及它的GUI可以参考这篇Blog
但是我们实验中不需要使用它的GUI，简单的调用可以参考如下过程：

1. 将下载好的Dlib库解压，复制用C++写的src文件夹到项目工程目录下
   

2. 将头文件引入项目中
   

3. 将all_console.cpp和all_gui.cpp的内容注释掉，因为我们不需要使用它的调试和gui功能，否则报错。

4. 在源文件中添加source.cpp,否则不会编译头文件。

   

   

5. 检测代码(模型文件下载地址shape_predictor_68_face_landmarks.dat.bz2)

Delaunay三角剖分

获取img1和img2的特征点后，用插值法获得中间图片的特征点位置，加上图片四个角顶点，然后对三张图片的特征点进行Delaunay三角剖分（一共68+4=72个顶点）。结果会将中间图片分成分成很多三角形块。由于Delaunay良好的性质，中间图片某个三角形的顶点一定与生成这个顶点的img1和img2中的特征点一一对应，所以为了计算方便我们可以考虑建立索引：

    用vector保存每张图片的特征点，并不再修改。把每个点在vector中的下标index当作索引，任何操作都是对下标进行操作，然后再从vector里返回点的坐标值。

中间图片特征点生成

    std::vector<Point> points;    CImg<unsigned char> img(img1.width(), img1.height(), 1, 3);    img.fill(0);    for (int i = 0; i < points1.size(); i++) {        float x, y;        x = (1 - alpha) * points1[i].x + alpha * points2[i].x;        y = (1 - alpha) * points1[i].y + alpha * points2[i].y;        Point temp = { x, y, 0 };        points.push_back(temp);    }

三角化算法:

Bowyer-Watson算法的基本步骤是：
1. 构造一个超级三角形，包含所有散点，放入三角形链表。
2. 将点集中的散点依次插入，在三角形链表中找出外接圆包含插入点的三角形（称为该点的影响三角形），删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来，完成一个点在Delaunay三角形链表中的插入。
3. 根据优化准则对局部新形成的三角形优化。将形成的三角形放入Delaunay三角形链表。
4. 循环执行上述第2步，直到所有散点插入完毕.

具体代码见工程代码“Delaunay.h”

确定像素值

因为img1、img2、中间图的三角块是一一对应的，中间图的每一个三角块中的像素值由img1和img2中的对应的三角块确定，所以每个三角块之间存在着Afine Projection的变换。我们需要求出每个三角块的坐标变换矩阵。求解的方法可以采用矩阵求逆的方式：

 A,B,C三个点的坐标构成3*3的方阵A=[ax,bc,cx; ay,by,cy; 1,1,1]，由Afine Projection的性质可以确定M=[a,b,c; d,e,f; 0,0,1],  投影矩阵A'= [a'x,b'x,c'x; a'y,b'y,c'y; 1,1,1]。 根据A'= MA, 可以推出M = A'*inv(A).

至于整个确定像素值的过程，我采用的是思路比较简单的做法，预先处理计算完中间图到img1和img2的所有变换矩阵得到m1, m2，存储在中间图片的三角序列里。然后遍历中间图片所有点，先判断它在哪个三角形内部，然后求出三角形的投影变换矩阵,当然更好的方法是直接遍历三角形中的点。

那么基本的问题如下：
1. 如何确定一个点在三角形内部？可以参考这个链接判断点在三角形内部
2. 矩阵求逆方法？完整代码参考工程代码”Mymatrix.h”

#define MYDIM 3// 计算每一行double getA(double** arcs, int n) {    if (n == 1) {        return arcs[0][0];    }    double ans = 0;    double** temp = getAMatrix(MYDIM, MYDIM);    int i, j, k;    for (i = 0; i<n; i++) {        for (j = 0; j<n - 1; j++) {            for (k = 0; k<n - 1; k++) {                temp[j][k] = arcs[j + 1][(k >= i) ? k + 1 : k];            }        }        double t = getA(temp, n - 1);        if (i % 2 == 0) {            ans += arcs[0][i] * t;        }        else {            ans -= arcs[0][i] * t;        }    }    deleteAMatrix(temp);    return ans;}//计算每一行每一列的每个元素所对应的余子式，组成A*void getAStart(double** arcs, int n, double** ans) {    if (n == 1) {        ans[0][0] = 1;        return;    }    int i, j, k, t;    double** temp = getAMatrix(MYDIM, MYDIM);    for (i = 0; i<n; i++) {        for (j = 0; j<n; j++) {            for (k = 0; k<n - 1; k++) {                for (t = 0; t<n - 1; t++) {                    temp[k][t] = arcs[k >= i ? k + 1 : k][t >= j ? t + 1 : t];                }            }            ans[j][i] = getA(temp, n - 1);  //此处顺便进行了转置            if ((i + j) % 2 == 1) {                ans[j][i] = -ans[j][i];            }        }    }    deleteAMatrix(temp);}//得到给定矩阵src的逆矩阵保存到des中。bool GetMatrixInverse(double** src, int n, double** des) {    double flag = getA(src, n);    double** t = getAMatrix(MYDIM, MYDIM);    if (0 == flag) {        std::cout << "原矩阵行列式为0，无法求逆。请重新运行" << std::endl;        return false;//如果算出矩阵的行列式为0，则不往下进行    }    else {        getAStart(src, n, t);        for (int i = 0; i<n; i++) {            for (int j = 0; j<n; j++) {                des[i][j] = t[i][j] / flag;            }        }    }    deleteAMatrix(t);    return true;}

双线性插值法确定像素值

unsigned char* Morphing::BilinearFilter(CImg<unsigned char> &src, Point * p) {    int xs = static_cast<int> (p->x);    int xl = xs + 1 > src.width() ? xs + 1 : xs;    int ys = static_cast<int> (p->y);    int yl = ys + 1 > src.height() ? ys + 1 : ys;    float ax = p->x - xs;    float by = p->y - ys;    float Si, Sj;    unsigned char *rgb = new unsigned char[3];    for (int i = 0; i < 3; i++) {        Si = src(xs, ys, 0, i) + ax* (src(xs, yl, 0, i) - src(xs, ys, 0, i));        Sj = src(xl, ys, 0, i) + ax* (src(xl, yl, 0, i) - src(xl, ys, 0, i));        rgb[i] = static_cast<unsigned char>(Si + by*(Sj - Si));    }    //cout << (int)rgb[0] << (int)rgb[1] << (int)rgb[2] << endl;    return rgb;}

执行过程

// alpha:变换率 outfile:保存图片路径void Morphing::doMorphing(float alpha, string outfile) {    // 异常排查    if (path1.empty() || path2.empty() ||         points1.empty() || points2.empty() ||        points1.size() != points2.size()) {        cout << "error!" << endl;        return;    }    // 生成中间图特征点    std::vector<Point> points;    CImg<unsigned char> img(img1.width(), img1.height(), 1, 3);    img.fill(0);    for (int i = 0; i < points1.size(); i++) {        float x, y;        x = (1 - alpha) * points1[i].x + alpha * points2[i].x;        y = (1 - alpha) * points1[i].y + alpha * points2[i].y;        Point temp = { x, y, 0 };        points.push_back(temp);    }    // 三角剖分    Delaunay delaunay = doTriangle(img, points);    // 预先算出所有三角形的映射矩阵    for (int i = 0; i < delaunay.m_Tris.size(); i++) {        int* v = delaunay.m_Tris[i].v;        double **A = makeAMatrix(points[v[0]], points[v[1]], points[v[2]]);        double **invA = getInverse(A);        double **U1 = makeAMatrix(points1[v[0]], points1[v[1]], points1[v[2]]);        double **U2 = makeAMatrix(points2[v[0]], points2[v[1]], points2[v[2]]);        delaunay.m_Tris[i].m1 = getMatrixProduct(U1, invA);        delaunay.m_Tris[i].m2 = getMatrixProduct(U2, invA);        deleteAMatrix(A);        deleteAMatrix(invA);        deleteAMatrix(U1);        deleteAMatrix(U2);    }    // 遍历中间图所有像素点    cimg_forXY(img, x, y) {        Point p(x, y, 0);        //cout << x << " " << y << endl;        for (int i = 0; i < delaunay.m_Tris.size(); i++) {            int* v = delaunay.m_Tris[i].v;            // 确定该点所在的三角形            bool in = PointinTriangle1(points[v[0]], points[v[1]], points[v[2]], p);            if (!in) continue;            double **M1 = delaunay.m_Tris[i].m1;            double **M2 = delaunay.m_Tris[i].m2;            // 矩阵乘法            double u1 = M1[0][0] * x + M1[0][1] * y + M1[0][2];            double v1 = M1[1][0] * x + M1[1][1] * y + M1[1][2];            Point *p1 = new Point(u1, v1, 0);            unsigned char* rgb1 = BilinearFilter(img1, p1);            double u2 = M2[0][0] * x + M2[0][1] * y + M2[0][2];            double v2 = M2[1][0] * x + M2[1][1] * y + M2[1][2];            Point *p2 = new Point(u2, v2, 0);            unsigned char* rgb2 = Biline            for (int k = 0; k < 3; k++) {                img(x, y, 0, k) = (1 - alpha) * rgb1[k] + alpha * rgb2[k];            }            delete p1, p2;            delete[] rgb1, rgb2;            break;        }    }    img.display();    img.save(outfile.c_str());}

结果演示

中间帧
alpha = 0.2

alpha = 0.5

alpha = 0.7

最后结果11帧

---

整个project代码参见github