基本思想
通過Dlib獲得當前人臉的特徵點,然後通過旋轉平移標準模型的特徵點進行擬合,計算標準模型求得的特徵點與Dlib獲得的特徵點之間的差,使用Ceres不斷迭代優化,最終得到最佳的旋轉和平移引數。
使用環境
系統環境:Ubuntu 18.04
使用語言:C++
編譯工具:CMake
第三方工具
Dlib:用於獲得人臉特徵點
Ceres:用於進行非線性優化
CMinpack:用於進行非線性優化 (OPTIONAL)
原始碼
https://github.com/Great-Keith/head-pose-estimation
基礎概念
旋轉矩陣
頭部的任意姿態可以轉化為6個引數(yaw, roll, pitch, tx, ty, tz),前三個為旋轉引數,後三個為平移引數。
平移引數好理解,原座標加上對應的變化值即可;旋轉引數需要構成旋轉矩陣,三個引數分別對應了繞y軸旋轉的角度、繞z軸旋轉的角度和繞x軸旋轉的角度。
具體程式碼實現我們可以通過Dlib已經封裝好的API,rotate_around_x/y/z(angle)
。該函式返回的型別是dlib::point_transform_affine3d
,可以通過括號進行三維的變形,我們將其封裝成一個rotate函式使用如下:
void rotate(std::vector<point3f>& points, const double yaw, const double pitch, const double roll)
{
dlib::point_transform_affine3d around_z = rotate_around_z(roll * pi / 180);
dlib::point_transform_affine3d around_y = rotate_around_y(yaw * pi / 180);
dlib::point_transform_affine3d around_x = rotate_around_x(pitch * pi / 180);
for(std::vector<point3f>::iterator iter=points.begin(); iter!=points.end(); ++iter)
*iter = around_z(around_y(around_x(*iter)));
}
[NOTE] 其中point3f是我自己定義的一個三維點座標型別,因為Dlib中並沒有提供,而使用OpenCV中的cv::Point3f會與dlib::point定義起衝突。定義如下:
typedef dlib::vector<double, 3> point3f;
[NOTE] Dlib中的dlib::vector不是std::vector,注意二者區分。
LM演算法
這邊不進行贅訴,建議跟著推導一遍高斯牛頓法,LM演算法類似於高斯牛頓法的進階,用於迭代優化求解非線性最小二乘問題。在該程式中使用Ceres/CMinpack封裝好的API(具體使用見後文)。
三維空間到二維平面的對映
根據針孔相機模型我們可以輕鬆的得到三維座標到二維座標的對映:
\(X^{2d}=f_x(\frac{X^{3d}}{Z^{3d}})+c_x\)
\(Y^{2d}=f_y(\frac{Y^{3d}}{Z^{3d}})+c_y\)
[NOTE] 使用上角標來表示3維座標還是2維座標,下同。
其中\(f_x, f_y, c_x, c_y\)為相機的內參,我們通過OpenCV官方提供的Calibration樣例進行獲取:
例如我的電腦所獲得的結果如下:
從圖中矩陣對應關係可以獲得對應的引數值。
#define FX 1744.327628674942
#define FY 1747.838275588676
#define CX 800
#define CY 600
[NOTE] 本程式不考慮外參。
具體步驟
獲得標準模型的特徵點
該部分可見前一篇文章:BFM使用 - 獲取平均臉模型的68個特徵點座標
我們將獲得的特徵點儲存在檔案 landmarks.txt
當中。
使用Dlib獲得人臉特徵點
該部分不進行贅訴,官方有給出了詳細的樣例。
具體可以參考如下樣例:
- https://github.com/davisking/dlib/blob/master/examples/face_landmark_detection_ex.cpp
- https://github.com/davisking/dlib/blob/master/examples/webcam_face_pose_ex.cpp(通過這個樣例可以學習OpenCV如何呼叫攝像頭)
其中使用官方提供的預先訓練好的模型,下載地址:http://dlib.net/files/shape_predictor_68_face_landmarks.dat.bz2
具體在程式碼中使用如下:
cv::Mat temp;
if(!cap.read(temp))
break;
dlib::cv_image<bgr_pixel> img(temp);
std::vector<rectangle> dets = detector(img);
cout << "Number of faces detected: " << dets.size() << endl;
std::vector<full_object_detection> shapes;
for (unsigned long j = 0; j < dets.size(); ++j) {
/* Use dlib to get landmarks */
full_object_detection shape = sp(img, dets[j]);
/* ... */
}
其中shape.part
就存放著我們通過Dlib獲得的當前人臉的特徵點二維點序列。
[NOTE] 在最後CMake配置的時候,需要使用Release
版本(最重要),以及增加選項USE_AVX_INSTRUCTIONS
和USE_SSE2_INSTRUCTIONS
/USE_SSE4_INSTRUCTIONS
,否則因為Dlib的檢測耗時較長,使用攝像頭即時擬合會有嚴重的卡頓。
使用Ceres進行非線性優化
Ceres的使用官方也提供了詳細的樣例,在此我們使用的是數值差分的方法,可參考:https://github.com/ceres-solver/ceres-solver/blob/master/examples/helloworld_numeric_diff.cc
Problem problem;
CostFunction* cost_function = new NumericDiffCostFunction<CostFunctor, ceres::RIDDERS, LANDMARK_NUM, 6>(new CostFunctor(shape));
problem.AddResidualBlock(cost_function, NULL, x);
Solver::Options options;
options.minimizer_progress_to_stdout = true;
Solver::Summary summary;
Solve(options, &problem, &summary);
std::cout << summary.BriefReport() << endl;
這裡我直接使用了數值差分的方法(NumericDiffCostFunction
),而不是使用自動差分(AutoDiffCostFunction
),是因為自動差分的CostFunctor是通過Template實現的,利用Template來實現Jacobian矩陣的計算使用的同一個結構,這樣的話下方旋轉矩陣就不能直接通過呼叫Dlib提供的三維座標旋轉介面,而是要將整個矩陣拆解開來實現(這邊暫時沒有細想到底能不能實現),因此出於簡便,使用數值差分,在準確性上是會受到影響的。
並且注意到,具體的方法使用了Ridders(ceres::RIDDERS
),而不是向前差分(ceres::FORWARD
)或者中分(ceres::CENTRAL
),因為用後兩者進行處理的時候,LM演算法\(\beta_{k+1}=\beta_k-(J^TJ+\lambda I)^{-1}J^Tr)\)的更新項為0,無法進行迭代,暫時沒有想到原因,之前這裡也被卡了很久。
[NOTE] 原始碼中還有使用了CMinpack的版本,該版本不可用的原因也是使用了封裝最淺的lmdif1_
呼叫(返回結果INFO=4),該版本下使用的向前差分,如果改為使用lmdif_
對其中的一些引數進行調整應該是可以實現的。
CostFunctor的構建
CostFunctor的構建是Ceres,也是這個程式,最重要的部分。首先我們需要先把想要計算的式子寫出來:
\(Q=\sum_i^{LANDMARK\_NUM} \|q_i^{2d}-p_i^{2d}\|^2\)
\(Q=\sum_i^{LANDMARK\_NUM} \|q_i^{2d}-Map(R(yaw, roll, pitch)p_i^{3d}+T(t_x, t_y, t_z))\|^2\)。
其中:
- LANDMARK_NUM:該程式中為68,因為Dlib演算法獲得的特徵點數為68;;
- \(q_i^{2d}\):通過Dlib獲得的2維特徵點座標,大小為68的vector<dlib::point>
- \(p_i^{2d}\):經過一系列變換得到的標準模型的2維特徵點座標,大小為68的vector<dlib::point>,具體計算方法是通過\(p_i^{2d}=Map(R(yaw, roll, pitch)(p_i^{3d})+T(t_x, t_y, t_z))\);
- \(p_i^{3d}\):標準模型的三維3維特徵點座標,大小為68的vector<point3f>;
- \(R(yaw, roll, pitch)\):旋轉矩陣;
- \(T(t_x, t_y, t_z)\):平移矩陣;
- \(Map()\):3維點轉2維點的對映,如上所描述通過相機內參獲得。
- \(\|·\|\):因為是兩個2維點的差,我們使用歐幾里得距離來作為2點的差。
Ceres當中的CostFunctor只需要寫入平方以內的內容,因此我們如下構建:
struct CostFunctor {
public:
CostFunctor(full_object_detection _shape){ shape = _shape; }
bool operator()(const double* const x, double* residual) const {
/* Init landmarks to be transformed */
fitting_landmarks.clear();
for(std::vector<point3f>::iterator iter=model_landmarks.begin(); iter!=model_landmarks.end(); ++iter)
fitting_landmarks.push_back(*iter);
transform(fitting_landmarks, x);
std::vector<point> model_landmarks_2d;
landmarks_3d_to_2d(fitting_landmarks, model_landmarks_2d);
/* Calculate the energe (Euclid distance from two points) */
for(int i=0; i<LANDMARK_NUM; i++) {
long tmp1 = shape.part(i).x() - model_landmarks_2d.at(i).x();
long tmp2 = shape.part(i).y() - model_landmarks_2d.at(i).y();
residual[i] = sqrt(tmp1 * tmp1 + tmp2 * tmp2);
}
return true;
}
private:
full_object_detection shape; /* 3d landmarks coordinates got from dlib */
};
其中的引數x是一個長度為6的陣列,對應了我們要獲得的6個引數。
初始值的選定
當前並沒有多考慮這個因素,在landmark-fitting-cam程式中除了第一幀的初始值是提前設定好的以外,後續的初始值都是前一幀的最優值。
後面的表現都很好,但這第一幀確實會存在紊亂的情況。
因此後續優化可以考慮使用一個粗估計的初始值,因為對於這些迭代優化方法,初始值的選擇決定了會不會陷入區域性最優的情況。
測試結果
臉部效果:
輸出工作環境: