頭部姿態(tài)估計(jì)

mediatv 2020-01-14

展開全文

基本思想（update)

通過(guò)Dlib獲得當(dāng)前人臉的特征點(diǎn)，然后通過(guò)(1)修改模型的幾何形狀參數(shù)和(2)旋轉(zhuǎn)平移模型，進(jìn)行擬合，計(jì)算標(biāo)準(zhǔn)模型求得的特征點(diǎn)與Dlib獲得的特征點(diǎn)之間的差，使用Ceres不斷迭代優(yōu)化，最終得到最佳的(1)模型幾何形狀參數(shù)和(2)旋轉(zhuǎn)和平移參數(shù)。

使用環(huán)境

系統(tǒng)環(huán)境：Ubuntu 18.04
使用語(yǔ)言：C++
編譯工具：CMake + VSCode

第三方工具

Dlib：用于獲得人臉特征點(diǎn)
Ceres：用于進(jìn)行非線性優(yōu)化
CMinpack：用于進(jìn)行非線性優(yōu)化 （OPTIONAL）
HDF5：用于讀取.h5格式數(shù)據(jù)（new）

源代碼

https://github.com/Great-Keith/head-pose-estimation/tree/master/cpp

基礎(chǔ)概念

旋轉(zhuǎn)矩陣(update)

頭部的任意姿態(tài)可以轉(zhuǎn)化為6個(gè)參數(shù)（yaw, roll, pitch, tx, ty, tz），前三個(gè)為旋轉(zhuǎn)參數(shù)，后三個(gè)為平移參數(shù)。
平移參數(shù)好理解，原坐標(biāo)加上對(duì)應(yīng)的變化值即可；旋轉(zhuǎn)參數(shù)需要構(gòu)成旋轉(zhuǎn)矩陣，三個(gè)參數(shù)分別對(duì)應(yīng)了繞y軸旋轉(zhuǎn)的角度、繞z軸旋轉(zhuǎn)的角度和繞x軸旋轉(zhuǎn)的角度。

具體代碼實(shí)現(xiàn)我們可以通過(guò)Dlib已經(jīng)封裝好的API，rotate_around_x/y/z(angle)。該函數(shù)返回的類型是dlib::point_transform_affine3d，可以通過(guò)括號(hào)進(jìn)行三維的變形，我們將其封裝成一個(gè)rotate函數(shù)使用如下：

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是我自己定義的一個(gè)三維點(diǎn)坐標(biāo)類型，因?yàn)镈lib中并沒有提供，而使用OpenCV中的cv::Point3f會(huì)與dlib::point定義起沖突。定義如下：

typedef dlib::vector<double, 3> point3f;
typedef dlib::point point2d;

[NOTE] Dlib中的dlib::vector不是std::vector，注意二者區(qū)分。

為了嘗試使用分析式的優(yōu)化求解，我們可以選擇不采用dlib的API進(jìn)行旋轉(zhuǎn)，而是通過(guò)直接將其數(shù)學(xué)式子列出，如下：

/* We use Z1Y2X3 format of Tait–Bryan angles */
double c1 = cos(roll  * pi / 180.0), s1 = sin(roll  * pi / 180.0);
double c2 = cos(yaw   * pi / 180.0), s2 = sin(yaw   * pi / 180.0);
double c3 = cos(pitch * pi / 180.0), s3 = sin(pitch * pi / 180.0);

for(std::vector<point3f>::iterator iter=landmarks3d.begin(); iter!=landmarks3d.end(); ++iter) {
    double X = iter->x(), Y = iter->y(), Z = iter->z(); 
    iter->x() = ( c1 * c2) * X + (c1 * s2 * s3 - c3 * s1) * Y + ( s1 * s3 + c1 * c3 * s2) * Z + tx;
    iter->y() = ( c2 * s1) * X + (c1 * c3 + s1 * s2 * s3) * Y + (-c3 * s1 * s2 - c1 * s3) * Z + ty;
    iter->z() = (-s2     ) * X + (c2 * s3               ) * Y + ( c2 * c3               ) * Z + tz; 
}

這樣同樣也完成了旋轉(zhuǎn)過(guò)程，并且方便進(jìn)行求導(dǎo)。
[注] 但目前還是選擇使用數(shù)學(xué)求解的方法，因?yàn)榉治鍪角蠼饨?jīng)常會(huì)出現(xiàn)不收斂的情況，原因未知。

LM算法

這邊不進(jìn)行贅訴，建議跟著推導(dǎo)一遍高斯牛頓法，LM算法類似于高斯牛頓法的進(jìn)階，用于迭代優(yōu)化求解非線性最小二乘問(wèn)題。在該程序中使用Ceres/CMinpack封裝好的API（具體使用見后文）。

三維空間到二維平面的映射

針孔模型

根據(jù)針孔相機(jī)模型我們可以輕松的得到三維坐標(biāo)到二維坐標(biāo)的映射：
$x = f_{x} \frac{X}{Z} + c_{x}$
$y = f_{y} \frac{Y}{Z} + c_{y}$
[NOTE] 使用上角標(biāo)來(lái)表示3維坐標(biāo)還是2維坐標(biāo)，下同。
其中 $f_{x}, f_{y}, c_{x}, c_{y}$ 為相機(jī)的內(nèi)參，我們通過(guò)OpenCV官方提供的Calibration樣例進(jìn)行獲?。?br> 例如我的電腦所獲得的結(jié)果如下：

從圖中矩陣對(duì)應(yīng)關(guān)系可以獲得對(duì)應(yīng)的參數(shù)值。

#define FX 1744.327628674942
#define FY 1747.838275588676
#define CX 800
#define CY 600

去除z軸的直接映射（new）

在程序中使用直接去除z軸的直接映射先進(jìn)行一次迭代求解，用該解來(lái)作為真正求解過(guò)程的初始值，詳細(xì)內(nèi)容見下方初始值的選定。
$x = X$
$y = Y$

BFM模型的表示(new)

BFM模型是通過(guò)PCA的方法得到的，可參見之前博文：BFM模型介紹及可視化實(shí)現(xiàn)（C++）
需要了解的是，雖然人臉幾何有上萬(wàn)個(gè)頂點(diǎn)，但我們可以通過(guò)199個(gè)PC系數(shù)來(lái)進(jìn)行表示，這也就是我們要求得的最終參數(shù)之一。

具體步驟

獲得模型的特征點(diǎn)(new)

根據(jù)BFM模型介紹及可視化實(shí)現(xiàn)（C++）中的介紹，我將其中對(duì)bfm模型的操作的代碼整合到了該項(xiàng)目當(dāng)中，并且將其中的數(shù)據(jù)類型根據(jù)實(shí)際使用dlib進(jìn)行的改進(jìn)：
*　使用point3f（dlib::vector<double, 3>）/ point2d（dlib::vector<int, 2>）替換自己寫的vec類：使得整體操作更為流暢和統(tǒng)一；

使用dlib::matrix替換std::vector<...>/ std::vector<std::vector<...>>：統(tǒng)一向量和矩陣，方便運(yùn)算；讀入數(shù)據(jù)也更加簡(jiǎn)單；在矩陣乘法等計(jì)算量比較大的運(yùn)算當(dāng)中，dlib/opencv有對(duì)此進(jìn)行優(yōu)化，速度會(huì)大大提升。
[注] 盡管如此，在模型幾何形狀的迭代上速度依舊很慢，在我的電腦上完成一次收斂的交替迭代需要40min左右。也是因?yàn)檫@個(gè)時(shí)間關(guān)系，如果要進(jìn)行real-time的頭部姿態(tài)捕捉的話，可以考慮先拍攝一張用戶照片作為輸入，后續(xù)的迭代中不對(duì)幾何形狀進(jìn)行優(yōu)化，在DVP其視頻上的處理就是這樣的。

模型獲得之后，我們根據(jù)GitHub上的這個(gè)朋友github.com/anilbas/BFMLandmarks給出的68個(gè)特征點(diǎn)下標(biāo)獲得對(duì)應(yīng)的特征點(diǎn)。

獲得標(biāo)準(zhǔn)模型的特征點(diǎn)(deprecated)

該部分可見前一篇文章：BFM使用 - 獲取平均臉模型的68個(gè)特征點(diǎn)坐標(biāo)
我們將獲得的特征點(diǎn)保存在文件 landmarks.txt 當(dāng)中。

使用Dlib獲得人臉特征點(diǎn)

該部分不進(jìn)行贅訴，官方有給出了詳細(xì)的樣例。
具體可以參考如下樣例：

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（通過(guò)這個(gè)樣例可以學(xué)習(xí)OpenCV如何調(diào)用攝像頭）

其中使用官方提供的預(yù)先訓(xùn)練好的模型，下載地址：http:///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就存放著我們通過(guò)Dlib獲得的當(dāng)前人臉的特征點(diǎn)二維點(diǎn)序列。

[NOTE] 在最后CMake配置的時(shí)候，需要使用Release版本（最重要），以及增加選項(xiàng)USE_AVX_INSTRUCTIONS和USE_SSE2_INSTRUCTIONS/USE_SSE4_INSTRUCTIONS，否則因?yàn)镈lib的檢測(cè)耗時(shí)較長(zhǎng)，使用攝像頭即時(shí)擬合會(huì)有嚴(yán)重的卡頓。

使用Ceres進(jìn)行非線性優(yōu)化(update)

Ceres的使用官方也提供了詳細(xì)的樣例，在此我們使用的是數(shù)值差分的方法，可參考：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;

[注] 具體的方法使用了Ridders（ceres::RIDDERS），而不是向前差分（ceres::FORWARD）或者中分（ceres::CENTRAL）,因?yàn)橛煤髢烧哌M(jìn)行處理的時(shí)候，LM算法 $β_{k + 1} = β_{k} - (J^{T} J + λ I)^{- 1} J^{T} r)$ 的更新項(xiàng)為0，無(wú)法進(jìn)行迭代，暫時(shí)沒有想到原因，之前這里也被卡了很久。
[注] 源代碼中還有使用了CMinpack的版本，該版本不可用的原因也是使用了封裝最淺的lmdif1_調(diào)用（返回結(jié)果INFO=4），該版本下使用的向前差分，如果改為使用lmdif_對(duì)其中的一些參數(shù)進(jìn)行調(diào)整應(yīng)該是可以實(shí)現(xiàn)的。

HeadPoseCostFunctor的構(gòu)建 - 頭部姿態(tài)的6個(gè)參數(shù)(update)

CostFunctor的構(gòu)建是Ceres，也是這個(gè)程序，最重要的部分。首先我們需要先把想要計(jì)算的式子寫出來(lái)：
$Q = \sum_{i}^{L A N D M A R K_N U M} ‖ q_{i} - p_{i} ‖^{2}$
$Q = \sum_{i}^{L A N D M A R K_N U M} ‖ q_{i} - \prod (R (y a w, r o l l, p i t c h) P_{i} + T (t_{x}, t_{y}, t_{z})) ‖^{2}$ 。
其中：

LANDMARK_NUM：該程序中為68，因?yàn)镈lib算法獲得的特征點(diǎn)數(shù)為68；；
$q_{i}$ ：通過(guò)Dlib獲得的2維特征點(diǎn)坐標(biāo)，大小為68的vector<dlib::point>
$p_{i}$ ：經(jīng)過(guò)一系列變換得到的標(biāo)準(zhǔn)模型的2維特征點(diǎn)坐標(biāo)，大小為68的vector<dlib::point>，具體計(jì)算方法是通過(guò) $p_{i}^{2 d} = M a p (R (y a w, r o l l, p i t c h) (p_{i}^{3 d}) + T (t_{x}, t_{y}, t_{z}))$ ；
$p_{i}$ ：標(biāo)準(zhǔn)模型的三維3維特征點(diǎn)坐標(biāo)，大小為68的vector<point3f>；
$R (y a w, r o l l, p i t c h)$ ：旋轉(zhuǎn)矩陣；
$T (t_{x}, t_{y}, t_{z})$ ：平移矩陣；
$\prod ()$ ：3維點(diǎn)轉(zhuǎn)2維點(diǎn)的映射，如上所描述通過(guò)相機(jī)內(nèi)參獲得。
$‖ \cdot ‖$ ：因?yàn)槭莾蓚€(gè)2維點(diǎn)的差，我們使用歐幾里得距離(的平方)來(lái)作為2點(diǎn)的差。
Ceres當(dāng)中的CostFunctor只需要寫入平方以內(nèi)的內(nèi)容，因此我們?nèi)缦聵?gòu)建：

struct HeadPoseCostFunctor {
public:
    HeadPoseCostFunctor(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 */
};

其中的參數(shù)x是一個(gè)長(zhǎng)度為6的數(shù)組，對(duì)應(yīng)了我們要獲得的6個(gè)參數(shù)。

ShapeCostFunctor的構(gòu)建 - 模型幾何的199個(gè)參數(shù)(new)

同HeadPoseCostFunctor的構(gòu)建，不同的是求解的參數(shù)改為幾何形狀的PC個(gè)數(shù)．

double head_pose[6] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
double shape_pc_basis[N_PC] = {0.f};

Problem head_pose_problem, shape_pc_problem;
CostFunction* head_pose_cost_function = new NumericDiffCostFunction<HeadPoseNumericCostFunctor, ceres::RIDDERS, LANDMARK_NUM, 6>(new HeadPoseNumericCostFunctor(obj_detection));
head_pose_problem.AddResidualBlock(head_pose_cost_function, NULL, head_pose);
CostFunction* shape_pc_cost_function = new NumericDiffCostFunction<ShapeNumericCostFunctor, ceres::RIDDERS, LANDMARK_NUM, N_PC>(new ShapeNumericCostFunctor(obj_detection));
shape_pc_problem.AddResidualBlock(shape_pc_cost_function, NULL, shape_pc_basis);

通過(guò)交替求解，即先求得頭部姿態(tài)參數(shù)的最優(yōu)解，以此為基礎(chǔ)去求幾何形狀參數(shù)的最優(yōu)解，以此往復(fù)，直到二者都達(dá)到收斂。偽代碼如下：

while(true) {
    Solve(head_pose);
    Solve(shape_pc_basis);
    if(two problems are both convergence)
        break;
}

通過(guò)該方法對(duì)同一張圖片進(jìn)行測(cè)試，記錄最終cost，如下：

	初始值	僅頭部姿態(tài)參數(shù)	交替迭代
cost(數(shù)量級(jí))	1e9	1e5	4e3

初始值的選定(update)

之前并沒有多考慮這個(gè)因素，在程序中除了第一幀的初始值是提前設(shè)置好的以外，后續(xù)的初始值都是前一幀的最優(yōu)值。
后面的表現(xiàn)都很好，但這第一幀確實(shí)會(huì)存在紊亂的情況（后來(lái)發(fā)現(xiàn)是超出最大迭代次數(shù)50后依舊NO_CONVERGENCE）。
因?yàn)閷?duì)于這些迭代優(yōu)化方法，初始值的選擇決定了會(huì)不會(huì)陷入局部最優(yōu)的情況，因此考慮使用一個(gè)粗估計(jì)的初始值。

粗暴的估計(jì)

通過(guò)例如鼻子的傾斜角度、嘴巴的傾斜角度、頭部的寬度等等直觀上的數(shù)據(jù)，進(jìn)行一個(gè)粗步估計(jì)。但影響效果一般。

使用直接映射進(jìn)行估計(jì)

不使用針孔相機(jī)模型，而是使用直接去掉z軸來(lái)進(jìn)行三維到二維的映射，在此基礎(chǔ)上進(jìn)行迭代，將迭代結(jié)果作為真正迭代的初始值。
影響效果一般，往往在這個(gè)迭代模型上已經(jīng)下降了3個(gè)數(shù)量級(jí)，但換作真正迭代的時(shí)候計(jì)算初始值的cost又會(huì)恢復(fù)到1e9級(jí)別。

后續(xù)考慮

了解到了DVP中初值使用POSIT算法（改進(jìn)版SoftPOSIT算法），后續(xù)考慮嘗試這種方法進(jìn)行估計(jì)。

測(cè)試結(jié)果(deprecated)

臉部效果：
only_face
輸出工作環(huán)境：
work_place

本站是提供個(gè)人知識(shí)管理的網(wǎng)絡(luò)存儲(chǔ)空間，所有內(nèi)容均由用戶發(fā)布，不代表本站觀點(diǎn)。請(qǐng)注意甄別內(nèi)容中的聯(lián)系方式、誘導(dǎo)購(gòu)買等信息，謹(jǐn)防詐騙。如發(fā)現(xiàn)有害或侵權(quán)內(nèi)容，請(qǐng)點(diǎn)擊一鍵舉報(bào)。

小男孩‘自慰网亚洲一区二区,亚洲一级在线播放毛片,亚洲中文字幕av每天更新,黄aⅴ永久免费无码,91成人午夜在线精品,色网站免费在线观看,亚洲欧洲wwwww在线观看