電腦視覺與模式識別（3）—— FaceMorphing

標籤：

        這是一個非常有意思的應用，可以將一個人的臉逐漸過渡為另一個人的臉。花了大概1天完成了最基本的功能，大概3天去完善它，可能還有不少bug等著我去修改，不過先把目前的進展記錄下來吧。

        圖形庫：CImg

        環境：Win10、C++11、VS2013

        一：

        二：

思路：

1)      分別對映像A和映像B採集控制點（座標）

2)      利用Delaunay演算法劃分三角形，保證映像A和映像B的三角形一一對應且不重合。

3)      計算每對三角形的過渡三角形，過渡量為0~1。

4)      計算三角形到三角形的仿射變換矩陣。

5)      利用變換矩陣分別將映像A和映像B變換到過渡三角形網格。

6)      利用過渡量計算合成圖的像素色彩。

        我覺得痛點主要在於 Delaunay演算法 和 Affine Transformation矩陣計算，其次是應用的魯棒性增強。由於目前還沒有學習Face Service的相關技術，所以還需要人為標定控制點，做一些“Dirty Work”。因此我也把工程分為了兩個：一是主工程，專門負責已經三角分割後的影像處理。二是輔助工程，專門負責控制點的採集以及三角分割。可以這麼說，主工程的核心是Morphing，輔助工程的核心是Delaunay。

        先看看輔助工程吧，事先定義好了常用的結構體後，就可以寫Delaunay演算法了。Delaunay分割的準則是：任何三角形ABC，不存在一點D，在其外接圓內。考慮到人臉的控制點數量一般在低兩位元，電腦平均每秒可計算百萬層級的資料，也就是支援每秒O(n^4)的計算量，那麼可以採用暴力演算法：

void Delaunay::buildDelaunayEx(const std::vector<Dot*>& vecDot, std::vector<Triangle*>& vecTriangle){int size = vecDot.size();if (size < 3)return;for (int i = 0; i < size - 2; ++i){for (int j = i + 1; j < size - 1; ++j){for (int k = j + 1; k < size; ++k){Dot* A = vecDot[i];Dot* B = vecDot[j];Dot* C = vecDot[k];Triangle* tri = new Triangle(*A, *B, *C);getTriangleCircle(*tri);bool find = true;for (int m = 0; m < size; ++m){Dot* P = vecDot[m];if (*P == *A || *P == *B || *P == *C) continue;if (pointInCircle(*P)){find = false;break;}}if (find){vecTriangle.push_back(tri);}}}}}

        我試了一下在控制點不超過40個時，是可以1秒跑完的。當然啦，這是一種最笨最簡單的方法。
        為了讓映像A和映像B的三角分割是一一對應的，我只對映像A進行了Delaunay分割，而映像B則根據映像A的三角構成來分割。因此，我會記錄映像A的每個點的序號，以及每個三角形三個點的序號組成。最終將這些序號以文本的方式輸出。

#include "CImg.h"#include "Delaunay.h"using namespace cimg_library;#include <iostream>#include <fstream>#include <string>using namespace std;int main(){// 輸入映像路徑string pathA, pathB;cout << "Image A path: "; cin >> pathA;cout << "Image B path: "; cin >> pathB;CImg<double> sourceA(pathA.c_str()), sourceB(pathB.c_str());CImgDisplay Adisp(sourceA, "Image A"), Bdisp(sourceB, "Image B");// 控制點數組vector<Dot*> vecDotA, vecDotB;// 映像角點預先置入數組vecDotA.push_back(new Dot(0, 0, 0));vecDotA.push_back(new Dot(0, sourceA.height() - 1, 1));vecDotA.push_back(new Dot(sourceA.width() - 1, sourceA.height() - 1, 2));vecDotA.push_back(new Dot(sourceA.width() - 1, 0, 3));vecDotB.push_back(new Dot(0, 0, 0));vecDotB.push_back(new Dot(0, sourceB.height() - 1, 1));vecDotB.push_back(new Dot(sourceB.width() - 1, sourceB.height() - 1, 2));vecDotB.push_back(new Dot(sourceB.width() - 1, 0, 3));// 點線顏色int color[3] = { 0, 255, 0 };// 點擊滑鼠擷取座標while (!Adisp.is_closed()){Adisp.wait();if (Adisp.button() & 1 && Adisp.mouse_y() >= 0){Dot* click = new Dot(Adisp.mouse_x(), Adisp.mouse_y(), vecDotA.size());sourceA.draw_circle(click->x, click->y, sourceA.width() / 40, color);sourceA.display(Adisp);vecDotA.push_back(click);}}while (!Bdisp.is_closed()){Bdisp.wait();if (Bdisp.button() & 1 && Bdisp.mouse_y() >= 0){Dot* click = new Dot(Bdisp.mouse_x(), Bdisp.mouse_y(), vecDotB.size());sourceB.draw_circle(click->x, click->y, sourceB.width() / 40, color);sourceB.display(Bdisp);vecDotB.push_back(click);}}// 三角形數組vector<Triangle*> vecTriA, vecTriB;// 對映像A進行Delaunay三角形分割 Delaunay().buildDelaunayEx(vecDotA, vecTriA);// 同步映像B的三角形分割 for (int i = 0; i < vecTriA.size(); ++i) { Dot* A = vecDotB[vecTriA[i]->a.value]; Dot* B = vecDotB[vecTriA[i]->b.value]; Dot* C = vecDotB[vecTriA[i]->c.value]; vecTriB.push_back(new Triangle(*A, *B, *C)); } CImg<double> targetA(pathA.c_str()), targetB(pathB.c_str()); CImgDisplay TAdisp(targetA), TBdisp(targetB);// 點擊滑鼠逐步顯示三角形int i = 0;while (!TAdisp.is_closed()){TAdisp.wait();if (TAdisp.button() && i < vecTriA.size()){targetA.draw_line(vecTriA[i]->a.x, vecTriA[i]->a.y, vecTriA[i]->b.x, vecTriA[i]->b.y, color);targetA.draw_line(vecTriA[i]->a.x, vecTriA[i]->a.y, vecTriA[i]->c.x, vecTriA[i]->c.y, color);targetA.draw_line(vecTriA[i]->b.x, vecTriA[i]->b.y, vecTriA[i]->c.x, vecTriA[i]->c.y, color);targetA.display(TAdisp);targetB.draw_line(vecTriB[i]->a.x, vecTriB[i]->a.y, vecTriB[i]->b.x, vecTriB[i]->b.y, color);targetB.draw_line(vecTriB[i]->a.x, vecTriB[i]->a.y, vecTriB[i]->c.x, vecTriB[i]->c.y, color);targetB.draw_line(vecTriB[i]->b.x, vecTriB[i]->b.y, vecTriB[i]->c.x, vecTriB[i]->c.y, color);targetB.display(TBdisp);++i;}}// 擷取映像命名, 改變命名為txt格式string filenameA = pathA.substr(0, pathA.find_last_of(".")) + ".txt";string filenameB = pathB.substr(0, pathB.find_last_of(".")) + ".txt";// 控制點和三角形寫入文本ofstream outputA(filenameA, ios::out), outputB(filenameB, ios::out);outputA << "[Points]" << endl;outputB << "[Points]" << endl;for (int i = 0; i < vecDotA.size(); ++i){outputA << vecDotA[i]->x << "," << vecDotA[i]->y << endl;outputB << vecDotB[i]->x << "," << vecDotB[i]->y << endl;}outputA << "[Triangles]" << endl;outputB << "[Triangles]" << endl;for (int i = 0; i < vecTriA.size(); ++i){outputA << vecTriA[i]->a.value << "," << vecTriA[i]->b.value << "," << vecTriA[i]->c.value << endl;outputB << vecTriB[i]->a.value << "," << vecTriB[i]->b.value << "," << vecTriB[i]->c.value << endl;}}

        接下來是主工程，有了輔助工程提供的文本資料以後，主工程需要計算每對三角的過渡三角形：

// 計算過渡三角形Triangle* middleTriangle(Triangle* A, Triangle* B, float rate){float ax = rate*(A->a.x) + (1 - rate)*(B->a.x);float ay = rate*(A->a.y) + (1 - rate)*(B->a.y);float bx = rate*(A->b.x) + (1 - rate)*(B->b.x);float by = rate*(A->b.y) + (1 - rate)*(B->b.y);float cx = rate*(A->c.x) + (1 - rate)*(B->c.x);float cy = rate*(A->c.y) + (1 - rate)*(B->c.y);return new Triangle(Dot(ax, ay), Dot(bx, by), Dot(cx, cy));}

        然後計算變換矩陣係數：

Matrix3x3 Warp::TriangleToTriangle(float u0, float v0, float u1, float v1, float u2, float v2,float x0, float y0, float x1, float y1, float x2, float y2){// |A|int detA;detA = u0*v1 + u1*v2 + u2*v0 - u2*v1 - u0*v2 - u1*v0;// A*int A11, A12, A13, A21, A22, A23, A31, A32, A33;A11 = v1 - v2;A21 = -(v0 - v2);A31 = v0 - v1;A12 = -(u1 - u2);A22 = u0 - u2;A32 = -(u0 - u1);A13 = u1*v2 - u2*v1;A23 = -(u0*v2 - u2*v0);A33 = u0*v1 - u1*v0;Matrix3x3 result;result.a11 = (float)(x0*A11 + x1*A21 + x2*A31) / detA;result.a21 = (float)(y0*A11 + y1*A21 + y2*A31) / detA;result.a31 = (float)(A11 + A21 + A31) / detA;result.a12 = (float)(x0*A12 + x1*A22 + x2*A32) / detA;result.a22 = (float)(y0*A12 + y1*A22 + y2*A32) / detA;result.a32 = (float)(A12 + A22 + A32) / detA;result.a13 = (float)(x0*A13 + x1*A23 + x2*A33) / detA;result.a23 = (float)(y0*A13 + y1*A23 + y2*A33) / detA;result.a33 = (float)(A13 + A23 + A33) / detA;return result;}

        計算過渡三角形變換矩陣的目的是，把原映像轉換為過渡映像，使得兩幅原映像在一致的三角網格下進行色彩疊加(Morph)。為了達到這個目的，需要把原映像的三角網格利用變換矩陣映射到過渡三角形網格上。每一對三角形有一個變換矩陣，存入動態數組中。

// 計算三角變換矩陣vector<Matrix3x3*> vecAB, vecBA;for (int i = 0; i < vecTriA.size(); ++i){Matrix3x3 HAM = Warp::TriangleToTriangle(*vecTriA[i], *vecTriM[i]);vecAB.push_back(new Matrix3x3(HAM));Matrix3x3 HBM = Warp::TriangleToTriangle(*vecTriB[i], *vecTriM[i]);vecBA.push_back(new Matrix3x3(HBM));}

         三角映射比較簡單了，判斷每個像素點是否在原三角形內，然後計算映射三角形內的座標，最後把該座標的像素在原圖上做一個雙線性插值計算，就可以作為該目標像素點的色彩值了。以映像A為例：

// 對原圖A進行三角變換cimg_forXY(targetA, x, y){bool isFind = false;for (int i = 0; i < vecTriB.size(); ++i){if (vecTriB[i]->pointInTriangle(Dot(x, y))){float tx = x*vecBA[i]->a11 + y*vecBA[i]->a12 + vecBA[i]->a13;float ty = x*vecBA[i]->a21 + y*vecBA[i]->a22 + vecBA[i]->a23;if (tx >= 0 && tx < width && ty >= 0 && ty < height){cimg_forC(sourceA, c)targetA(x, y, 0, c) = sourceA.linear_atXY(tx, ty, 0, c);}isFind = true;break;}}if (!isFind){cimg_forC(sourceA, c)targetA(x, y, 0, c) = sourceA.linear_atXY(x, y, 0, c);}}

        最後合成靶心圖表的像素色彩：

// 計算色彩插值cimg_forXYZC(target, x, y, z, c)target(x, y, z, c) = rate*targetB(x, y, z, c) + (1 - rate)*targetA(x, y, z, c);

       

         到此，應用已經基本成型。更多的測試還在進行中，我希望最終能通過Canvas製作成網頁版的變臉應用。

電腦視覺與模式識別（3）—— FaceMorphing