头部姿态估计 - OpenCV/Dlib/Ceres

基本思想(update)

通过Dlib获得当前人脸的特征点,然后通过(1)修改模型的几何形状参数和(2)旋转平移模型,进行拟合,计算标准模型求得的特征点与Dlib获得的特征点之间的差,使用Ceres不断迭代优化,最终得到最佳的(1)模型几何形状参数和(2)旋转和平移参数。

使用环境

系统环境:Ubuntu 18.04
使用语言:C++
编译工具:CMake + VSCode

第三方工具

Dlib:用于获得人脸特征点
Ceres:用于进行非线性优化
CMinpack:用于进行非线性优化 (OPTIONAL)
HDF5:用于读取.h5格式数据(new)

源代码

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

基础概念

旋转矩阵(update)

头部的任意姿态可以转化为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;
typedef dlib::point point2d;

[NOTE] Dlib中的dlib::vector不是std::vector,注意二者区分。

为了尝试使用分析式的优化求解,我们可以选择不采用dlib的API进行旋转,而是通过直接将其数学式子列出,如下:

/* 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; 
}

这样同样也完成了旋转过程,并且方便进行求导。
[注] 但目前还是选择使用数学求解的方法,因为分析式求解经常会出现不收敛的情况,原因未知。

LM算法

这边不进行赘诉,建议跟着推导一遍高斯牛顿法,LM算法类似于高斯牛顿法的进阶,用于迭代优化求解非线性最小二乘问题。在该程序中使用Ceres/CMinpack封装好的API(具体使用见后文)。

三维空间到二维平面的映射

针孔模型

根据针孔相机模型我们可以轻松的得到三维坐标到二维坐标的映射:
\(x=f_x\frac{X}{Z}+c_x\)
\(y=f_y\frac{Y}{Z}+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

去除z轴的直接映射(new)

在程序中使用直接去除z轴的直接映射先进行一次迭代求解,用该解来作为真正求解过程的初始值,详细内容见下方初始值的选定
\(x=X\)
\(y=Y\)

BFM模型的表示(new)

BFM模型是通过PCA的方法得到的,可参见之前博文:BFM模型介绍及可视化实现(C++)
需要了解的是,虽然人脸几何有上万个顶点,但我们可以通过199个PC系数来进行表示,这也就是我们要求得的最终参数之一。

具体步骤

获得模型的特征点(new)

根据BFM模型介绍及可视化实现(C++)中的介绍,我将其中对bfm模型的操作的代码整合到了该项目当中,并且将其中的数据类型根据实际使用dlib进行的改进:
* 使用point3fdlib::vector<double, 3>)/ point2ddlib::vector<int, 2>)替换自己写的vec类:使得整体操作更为流畅和统一;

  • 使用dlib::matrix替换std::vector<...>/ std::vector<std::vector<...>>:统一向量和矩阵,方便运算;读入数据也更加简单;在矩阵乘法等计算量比较大的运算当中,dlib/opencv有对此进行优化,速度会大大提升。
    [注] 尽管如此,在模型几何形状的迭代上速度依旧很慢,在我的电脑上完成一次收敛的交替迭代需要40min左右。也是因为这个时间关系,如果要进行real-time的头部姿态捕捉的话,可以考虑先拍摄一张用户照片作为输入,后续的迭代中不对几何形状进行优化,在DVP其视频上的处理就是这样的。

模型获得之后,我们根据GitHub上的这个朋友github.com/anilbas/BFMLandmarks给出的68个特征点下标获得对应的特征点。

获得标准模型的特征点(deprecated)

该部分可见前一篇文章:BFM使用 - 获取平均脸模型的68个特征点坐标
我们将获得的特征点保存在文件 landmarks.txt 当中。

使用Dlib获得人脸特征点

该部分不进行赘诉,官方有给出了详细的样例。
具体可以参考如下样例:

其中使用官方提供的预先训练好的模型,下载地址: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_INSTRUCTIONSUSE_SSE2_INSTRUCTIONS/USE_SSE4_INSTRUCTIONS,否则因为Dlib的检测耗时较长,使用摄像头即时拟合会有严重的卡顿。

使用Ceres进行非线性优化(update)

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;

[注] 具体的方法使用了Ridders(ceres::RIDDERS),而不是向前差分(ceres::FORWARD)或者中分(ceres::CENTRAL),因为用后两者进行处理的时候,LM算法\(\beta_{k+1}=\beta_k-(J^TJ+\lambda I)^{-1}J^Tr)\)的更新项为0,无法进行迭代,暂时没有想到原因,之前这里也被卡了很久。
[注] 源代码中还有使用了CMinpack的版本,该版本不可用的原因也是使用了封装最浅的lmdif1_调用(返回结果INFO=4),该版本下使用的向前差分,如果改为使用lmdif_对其中的一些参数进行调整应该是可以实现的。

HeadPoseCostFunctor的构建 - 头部姿态的6个参数(update)

CostFunctor的构建是Ceres,也是这个程序,最重要的部分。首先我们需要先把想要计算的式子写出来:
\(Q=\sum_i^{LANDMARK\_NUM} \|q_i-p_i\|^2\)
\(Q=\sum_i^{LANDMARK\_NUM} \|q_i-\prod(R(yaw, roll, pitch)P_i+T(t_x, t_y, t_z))\|^2\)
其中:

  • LANDMARK_NUM:该程序中为68,因为Dlib算法获得的特征点数为68;;
  • \(q_i\):通过Dlib获得的2维特征点坐标,大小为68的vector<dlib::point>
  • \(p_i\):经过一系列变换得到的标准模型的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\):标准模型的三维3维特征点坐标,大小为68的vector<point3f>;
  • \(R(yaw, roll, pitch)\):旋转矩阵;
  • \(T(t_x, t_y, t_z)\):平移矩阵;
  • \(\prod()\):3维点转2维点的映射,如上所描述通过相机内参获得。
  • \(\|·\|\):因为是两个2维点的差,我们使用欧几里得距离(的平方)来作为2点的差。
    Ceres当中的CostFunctor只需要写入平方以内的内容,因此我们如下构建:
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 */
};

其中的参数x是一个长度为6的数组,对应了我们要获得的6个参数。

ShapeCostFunctor的构建 - 模型几何的199个参数(new)

同HeadPoseCostFunctor的构建,不同的是求解的参数改为几何形状的PC个数.

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);

通过交替求解,即先求得头部姿态参数的最优解,以此为基础去求几何形状参数的最优解,以此往复,直到二者都达到收敛。伪代码如下:

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

通过该方法对同一张图片进行测试,记录最终cost,如下:

初始值 仅头部姿态参数 交替迭代
cost(数量级) 1e9 1e5 4e3

初始值的选定(update)

之前并没有多考虑这个因素,在程序中除了第一帧的初始值是提前设置好的以外,后续的初始值都是前一帧的最优值。
后面的表现都很好,但这第一帧确实会存在紊乱的情况(后来发现是超出最大迭代次数50后依旧NO_CONVERGENCE)。
因为对于这些迭代优化方法,初始值的选择决定了会不会陷入局部最优的情况,因此考虑使用一个粗估计的初始值。

粗暴的估计

通过例如鼻子的倾斜角度、嘴巴的倾斜角度、头部的宽度等等直观上的数据,进行一个粗步估计。但影响效果一般。

使用直接映射进行估计

不使用针孔相机模型,而是使用直接去掉z轴来进行三维到二维的映射,在此基础上进行迭代,将迭代结果作为真正迭代的初始值。
影响效果一般,往往在这个迭代模型上已经下降了3个数量级,但换作真正迭代的时候计算初始值的cost又会恢复到1e9级别。

后续考虑

了解到了DVP中初值使用POSIT算法(改进版SoftPOSIT算法),后续考虑尝试这种方法进行估计。

测试结果(deprecated)

脸部效果:
only_face
输出工作环境:
work_place

posted @ 2019-07-26 23:28  Bemfoo  阅读(5278)  评论(0编辑  收藏  举报