CN106997605A - 一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法 - Google Patents

Abstract

本发明公开了一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法。用带IMU传感器的智能手机绕着裸脚拍摄360度获取脚型视频和对应的IMU数据，通过建立因子图，把摄像机姿态、三维点云坐标、IMU数据、视频轮廓等已有的各种数据做为因子，通过定义各个因子之间的约束误差方程并进行求解而得到摄像机参数，引入参考模型进行脚型重建。充分利用参考模型的点云分布、法向量、表面曲率等信息引导重建脚型的形变，使之在不偏离重建点云的情况下，尽量接近真实脚型，减小噪声点的干扰，从而实现一个方便地得到对应的脚型三维网格模型的过程。既降低了脚型数据的获取成本，又不失参数计算的精确性和鲁棒性。

Description

一种通过智能手机采集脚型视频和传感器数据获取三维脚型 的方法

技术领域

本发明涉及一种三维重建方法，尤其涉及一种基于带摄像机和传感器的智能手机通过采集脚型RGB视频和传感器数据并结合数据库中的参考脚型模型进行三维脚型数据获取的方法。

背景技术

近两年对结合传感器和三维重建算法的研究越来越多，产生了大量高质量的三维重建算法，迅速提升了基于移动终端的三维重建技术。对三维重建算法进行总结，可以将基于移动终端的三维重建算法大致分为如下几类：基于深度图融合的三维重建、基于IMU和视觉融合的三维重建、基于光照条件和激光扫描等其它传感器的三维重建。

基于深度图融合的三维重建，由深度传感器记录拍摄场景的深度信息，然后根据SFM算法或其它算法计算摄像机参数，并用该摄像机参数恢复三维空间坐标点并反投影到视频上，最小化反投影误差，如主要研究深度图的融合过程，即利用高斯模型或其它数据模型对同时摄像机参数和三维点云进行约束，过滤掉偏差较大的值。另一种算法是由建立体素，根据各个角度的深度图及摄像机参数计算当前位置的深度并进行融合，最后形成目标物体的体积网格。但这种方法运行时需要占用大量的内存，计算速度慢，不利用在线实时重建。

基于IMU和视觉融合的三维重建，比较成熟的算法中，主要考虑到视觉三维重建算法对纹理要求高的缺点，利用IMU单元的陀螺仪和加速计记录移动设备的朝向和加速度，并假设噪声服从高斯分布，借助Kalman滤波对IMU数据进行预测，得到视频对应的摄像机参数，并融入到视觉算法中，由计算的三维坐标点和二维特征点之间、两幅视频的二维特征点之间的约束关***一优化摄像机参数。在此基础上，Adrian等人加入了闭环检测优化姿态跟踪，即假设移动设备在拍摄过程中的移动成环形封闭，利用最后一帧视频和初始帧视频部分重叠的性质进一步约束摄像机参数以提高其精确度。

基于光照条件和激光扫描等其它传感器的三维重建，借助高级设备如激光传感器对目标重建物体进行扫描融合得到对应的三维模型。除此之外还有利用光源信息进行重建的算法，即利用光源传感器对目标重建物体从不同角度拍摄，或从同一角度用不同的光照拍摄，在光源信息已知的情况下通过拍摄视频RGB值及光源方向、物体表面法向量等的代数约束重建目标物体。此类方法尽管能够达到比较高的重建精度，但由于数据采集费用大，操作麻烦，不适用于一般性的重建中。

针对于脚型重建的科学研究则比较少。更多的是在医学研究上，用CT(ComputedTomography)电子计算机断层扫描视频对脚型进行三维重建以观察脚内部的病变情况，如用Laplacian表面形变框架对由CT视频重建出来的模型进行形变，但需要用户手工输入各项参数，而且仅限于医学研究领域，昂贵的CT拍摄设备及对用户的相关专业知识的高要求都使得基于CT视频的三维重建技术无法应用到普通用户的生活中。相对于室外场景，脚型的大小要求重建的误差足够小，因此很多室外场景的三维重建技术甚至室内的三维重建算法都因误差较大不适用于脚型重建。小误差的重建技术中激光扫描的重建效果最好，如结合轮廓、激光和视觉技术对脚型进行测量，得到的脚型参数误差可以达到5mm以内，或是对一系列利用结构光进行激光扫描重建三维脚型的技术进行对比，给出了利用激光扫描所能达到的精度范围。但激光扫描设备的昂贵也使得基于激光的重建技术无法普及。考虑到脚型模型在行走时与静止时的差异，给出了在行走状态下对脚型各项参数的测试方法，同样需要固定的大型设备。

在不需要大型昂贵设备的条件下，视觉重建算法占主导地位。如采用深度摄像机捆绑在移动小车上绕脚型拍摄大半圈进行脚型重建，拍摄过程中脚边放有标定板。深度相机虽然比激光设备便宜，却也并不是家庭常用设备。随着机器学习的发展与广泛应用，有越来越多的学者把三维重建研究技术和机器学习结合越来，最近的研究在通过脚型轮廓剔除外点后，利用GNG神经网络进行训练和学习，检测脚型模型中的标记得到网格图的参数，从而进行脚型重建。但网络学习需要大量的训练数据，数据收集和整理工作量巨大，重建结果依赖于数据和训练，相比于视觉和IMU融合的重建在误差精度上并没有太大的优势。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于带摄像机和传感器的智能手机通过采集脚型RGB视频和传感器数据并结合数据库中的参考脚型模型进行三维脚型数据获取的方法。

本发明的目的是通过以下技术方案来实现的：一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，该方法包括以下步骤：

(1)利用智能手机围绕裸脚旋转一圈采集脚型视频和传感器数据，将视频分解成视频帧，并取部分帧作脚型区域分割等预处理，同时对脚型数据库中的参考脚型模型求法向量等参数；

(2)根据多视角立体几何的约束关系，检测视频帧中的特征点并进行匹配，计算摄像机内外参数，恢复特征点对应的三维点云坐标；

(3)加入IMU数据，对步骤2得到的摄像机参数和三维点云坐标进行修正，包括修复多视角算法失败无法求得摄像机参数的情况和移动三维点云坐标以减小三维坐标计算误差；

(4)从已有的脚型数据库中选择相似的参考脚型模型，对齐重建脚型模型和参考脚型模型，并对重建脚型模型进行去噪、光滑、表面重建处理；

(5)根据重建脚型模型和参考脚型模型之间的局部和整体对齐关系，对参考脚型模型进行形变使其拟合重建脚型模型，输出最终的重建脚型模型结果。

进一步地，所述步骤1具体包括以下子步骤：

(1.1)用带有陀螺仪和加速计的智能手机绕脚360度拍摄视频一圈回到原点，拍摄时在脚底放置一张A4纸，拍摄过程中记录陀螺仪和加速计的数值。拍摄场景中的物体尺寸与真实尺寸之间的比例以A4纸作为参考。

(1.2)将采集到的视频分解成视频帧，从视频帧中等间隔地选择部分帧进行分割前背景区域，得到视频帧对应的关于脚型区域的mask视频。通过检测视频帧中A4纸的角点和边获取A4纸区域，根据RGB像素值分割A4纸及在A4纸区域内的部分脚型区域。同时，对整幅视频帧进行边缘检测，得到包括脚型边缘在内的视频帧边缘视频。然后由A4纸区域内的脚型区域向外搜索填充，最终得到初步估计的脚型区域及非脚型区域。建立GMM模型，通过k-mean算法分别对脚型区域像素和非脚型区域像素进行聚类，得到GMM模型中每个高斯模型的像素样本集，然后根据RGB值计算每个高斯模型的参数均值和协方差，每个高斯模型的权值定义为属于该高斯模型的像素个数与总像素个数的比值。对每个像素分配GMM模型中的高斯模型，以给定的脚型区域和非脚型区域为参考学习优化GMM模型，迭代分割，最终得到对应视频帧的脚型轮廓视频。

(1.3)用Kalibr工具箱对IMU进行标定，求出IMU到摄像机的旋转矩阵，利用该旋转矩阵变换摄像机坐标系和IMU坐标系。对于参考脚型模型，顶点v_i的法向量由该顶点和周围k个邻近点共同确定。记顶点v_i的第j个邻近点为(j∈1,…,k)。模型中心点为c，则法向量同时满足：

其中Face{v_i}为包含点v_i的面片，几何意义为顶点v_i的法向量由以该点为顶点的面片的法向量的权重和，方向为由中心点指向顶点外侧。最终得到每个顶点和每个面片的法向量。

进一步地，所述步骤2具体为：对每帧视频帧进行Brief和Freak特征点检测和提取，对视频帧两两进行特征匹配，计算单应矩阵，并按纯视觉方法计算基本矩阵、本质矩阵，对本质矩阵进行SVD分解求得初始的摄像机姿态参数，同时恢复特征点对应的三维点云坐标。

进一步地，所述步骤3具体包括以下子步骤：

(3.1)对IMU数据做如下预处理。由于传感器频率和摄像机频率不一致，假设IMU数据和摄像机在从时间t_i到t_j的间隔时间Δt里同步，状态参数x_i＝[R_i,p_i,v_i,b_i]记录了IMU数据的旋转R_i、位置p_i、速度v_i和传感器漂移偏差b_i，而传感器漂移偏差b_i又包括速度加速计和重力加速计的偏差。则视频帧对应的状态参数为IMU状态参数的子集。设白噪声θ＝[θ^gd,θ^ad]，θ^gd和θ^ad分别为重力加速计和速度加速计的噪声，角速度ω_i，加速度a_i，则有：

状态参数求导后仍然为偏差b的函数，则可以采用分两步进行优化的方法，第一步假设b已知，对状态参数求梯度最优化，第二步中考虑到传感器移动时b有变化，为减少计算，当b变成b+Δb时，梯度值只需要在原来的基础上增加由b变化带来的差值Δb即可。最后计算残差，得到IMU数据下的摄像机姿态。

(3.2)结合IMU数据和视觉数据建立因子图。其中因子节点包括由视觉方法计算得到的摄像机姿态x_i′，恢复的脚型三维点云坐标l_i、IMU状态参数x_i，偏差c_i等。记第k时刻IMU测量值为则因子图中各个约束可以表示为：IMU因子：

偏差因子：

f^bias(c_k+1,c_k)＝d(c_k+1-g(c_k))

先概率因子：

视觉投影因子

f^Proj(x,l)＝d(z-π(x,l))

其中z表示真实的视频投影坐标，V为所有顶点的集合，π为投影函数。

立体视觉因子：

f^Stereo(x,l)＝d(z^R-π^R(x,l))d(z^L-π^L(x,l))

其中L和R分别对应左、右摄像头。

闭环检测因子：

f^Close＝sig(s_ij)*f(x_i,u_ij)-x_j

其中sig()为sigmoid函数，s_ij记录传感器运动方向(正向、逆向)，u_ij记录摄像机运动方向(正向、逆向)，而f(x_i,u_ij)为运动约束模型，可以看成是f^IMU。

(3.3)把因子图中因子函数分解成一般形式并进行优化，得到优化后的各个因子的值。

进一步地，所述步骤4具体包括以下子步骤：

(4.1)根据轮廓约束从重建出的三维脚型点云中提取出脚型的三维点云数据，并补脚底点集，得到完整的脚型点云模型。

(4.2)从已有的脚型数据库中选择与重建脚型相似的脚型参考模型，先以坐标轴和包围盒为参考对重建模型和参考模型进行初步对齐，然后去掉增加的脚底点云及脚踝以上的点避免干扰，只保留重建点云模型中脚面上的点集，利用ICP算法进一步与参考模型对齐，最后重新加上脚底点云，并增加手工调整，确保重建脚型模型与参考脚型模型对齐。

(4.3)对重建脚型模型进行去噪和光滑处理，包括剔除重建脚型模型和参考脚型模型之间对应点对距离大于阈值的点，以及等密度去掉部分重建脚型点云。

(4.4)重建脚型模型的法向量计算和表面重建。先用邻近点计算法向量的初值，然后对点云中的点进行可见性判断，并根据输入点的八叉树选取空间中包含点云的单位立方体，对每个立方体中的点云，选取以目标点为圆心，半径h_i范围内的点进行平面拟合：

其中为单位法向量，x为范围内的邻近点，D为给定的范围阈值。最后以立方体中的点组成圈，判断每个面片的法向量夹角是否小于π/2，剔除大于该角度的面片，进行表面重建。

进一步地，所述步骤5具体包括以下子步骤：

(5.1)在参考脚型模型选择的区域中任意选定曲率线L，在L上取样本点：

L_s＝{(l_i,t_i),i＝1,…,numOfSample}

其中l_i,t_i分别表示该像素点的二维视频坐标和切线向量，numOfSample为样本点数。然后把这些样本点映射到参考脚型的曲面上得到：

其中l_i′,t_i′分别表示该投影点的三维空间坐标及沿L方向上的切线平面，而表示该样本点的法向量。把这些信息向曲面周围延伸，记x为该样本点的邻近点，即在与法向量和切平面都垂直的方向上移动一小步Δ＝{Δx,Δy}，定义一个新网格点表示原来的点移动Δ后所在的点。在曲率线的两边，即沿b_i和-b_i方向都进行这样的操作，直到生成一系列网格点覆盖选择区域，对这些网格点进行插值直到选择的视频区域中每个像素点h＝(i,j)都对应一个三维网格节点。然后建立目标能量函数，计算并输出深度映射值集合Z＝{z_h}，对目标能量函数运用梯度下降法输出区域的深度映射视频，在重建模型上选择对应的区域根据曲率进行修复。目标能量函数共包括四项,为：

第一项表示曲面光滑程度,其中z_xx(i,j)表示在像素点h＝(i,j)沿x方向的深度，同理，z_xy(i,j)和z_yy(i,j)分别为沿xy斜线方向和y方向的深度。

第二项F(z_h)满足x|f(x)＝0，表示模型曲面是函数F的零集。

第三项R(z_h)为同一像素在两幅不同的视频帧中的深度差的平方值，保证在重叠区域两个深度映射值的一致性。

第四项C(z_h)＝∑_i,j(U^THU)²为曲率线的变化，其中H为z_h的Hessian矩阵，U为曲率线上的点投影到两帧视频帧上对应的像素。

(5.2)对重建脚型模型进行整体变换操作。最终要求的目标模型S‘可以表示为由参考脚型模型S^*和重建脚型模型S中的点在变换θ和移动δ＝{d_i}的作用下得到的。即S‘＝T(S^*,S)+d。考虑点集的GMM概率密度函数：

gmm(S)＝∑p(x)＝1/(2πσ²)exp(-||x-μ||/(2σ²))

p(x)为邻近点x的概率，μ为均值，σ为方差。

代价函数设为两个模型之间的L2距离

d_L2＝(gmm(S^*)-gmm(T(S,θ)))²

其中T为对点云进行变形操作，利用梯度下降法最小化目标函数得到模型间的形变系数θ。

(5.3)对重建脚型模型进行局部调整。定义形变目标函数和几何体目标函数Ω(S^*,S)分别为：

其中t_j为点p_j的移动向量，N(j)为点p_j的邻域，e_jk为连接和的边，A_jk为包含j和k两点的面片的面积，A_j为只包含点j的面片的面积。整体目标函数定义为通过对模型中点的拉伸移动改变模型。为保证整个模型的光滑度，考虑光滑目标模型S‘＝T(S^*,S,θ)+d，两点集间每对对应点之间的移动d_i可以通过最小化光滑目标函数求得。光滑目标函数E_smooth如下：

其中权重系数ε由模型数据集确定，N为对应模型上点的邻域。最后，计算点云的PCA模型矩阵M，m是参数向量，定义形变代价函数：

其中θ和t分别为旋转和平移参数,s为缩放因子，R为旋转矩阵。对用户选择的每部分区域进行形变，最终得到调整后的重建脚型模型。

本发明的有益效果是：本发明用带IMU传感器的智能手机绕着裸脚拍摄360度获取脚型视频和对应的IMU数据，通过建立因子图，把摄像机姿态、三维点云坐标、IMU数据、视频轮廓等已有的各种数据做为因子，通过定义各个因子之间的约束误差方程并进行求解而得到摄像机参数，引入参考模型进行脚型重建。充分利用参考模型的点云分布、法向量、表面曲率等信息引导重建脚型的形变，使之在不偏离重建点云的情况下，尽量接近真实脚型，减小噪声点的干扰，从而实现一个方便地得到对应的脚型三维网格模型的过程。本发明方法使得在三维脚型重建过程中，当视觉重建算法因为特征点过少或其它原因失败时，能够用传感器预测的摄像机姿态进行修正，增强了重建的鲁棒性，同时结合参考脚型模型和传感器、脚型三维点、视频特征点等各项数据相互约束，能够有效地减小重建误差。实现了一个自动脚型模型生成***，能够低成本地重建出脚型三维模型并计算脚型的各项测量参数，解决了远程定制鞋子脚型参数测量的不便性。

附图说明

图1是本发明方法的总体流程图；

图2(a)是输入的原始视频数据；

图2(b)是输入的参考模型样例；

图3是步骤3)中所构建的因子图；

图4是所求得的最终脚型模型。

具体实施方式

下面结合附图和具体实施例对本发明作进一步相信说明。

本发明提供的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，该方法包括以下步骤：

(1)利用智能手机围绕裸脚旋转一圈采集脚型视频和传感器数据，将视频分解成视频帧，并取部分帧作脚型区域分割等预处理，同时对脚型数据库中的参考脚型模型求法向量等参数；

(2)根据多视角立体几何的约束关系，检测视频帧中的特征点并进行匹配，计算摄像机内外参数，恢复特征点对应的三维点云坐标；

(3)加入IMU数据，对步骤2得到的摄像机参数和三维点云坐标进行修正，包括修复多视角算法失败无法求得摄像机参数的情况和移动三维点云坐标以减小三维坐标计算误差；

(4)从已有的脚型数据库中选择相似的参考脚型模型，对齐重建脚型模型和参考脚型模型，并对重建脚型模型进行去噪、光滑、表面重建处理；

(5)根据重建脚型模型和参考脚型模型之间的局部和整体对齐关系，对参考脚型模型进行形变使其拟合重建脚型模型，输出最终的重建脚型模型结果。

进一步地，所述步骤1具体包括以下子步骤：

(1.1)用带有陀螺仪和加速计的智能手机绕脚360度拍摄视频一圈回到原点，拍摄时在脚底放置一张A4纸，拍摄过程中记录陀螺仪和加速计的数值。拍摄场景中的物体尺寸与真实尺寸之间的比例以A4纸作为参考。

(1.2)将采集到的视频分解成视频帧，从视频帧中等间隔地选择部分帧进行分割前背景区域，得到视频帧对应的关于脚型区域的mask视频。通过检测视频帧中A4纸的角点和边获取A4纸区域，根据RGB像素值分割A4纸及在A4纸区域内的部分脚型区域。同时，对整幅视频帧进行边缘检测，得到包括脚型边缘在内的视频帧边缘视频。然后由A4纸区域内的脚型区域向外搜索填充，最终得到初步估计的脚型区域及非脚型区域。建立GMM模型，通过k-mean算法分别对脚型区域像素和非脚型区域像素进行聚类，得到GMM模型中每个高斯模型的像素样本集，然后根据RGB值计算每个高斯模型的参数均值和协方差，每个高斯模型的权值定义为属于该高斯模型的像素个数与总像素个数的比值。对每个像素分配GMM模型中的高斯模型，以给定的脚型区域和非脚型区域为参考学习优化GMM模型，迭代分割，最终得到对应视频帧的脚型轮廓视频。

(1.3)用Kalibr工具箱对IMU进行标定，求出IMU到摄像机的旋转矩阵，利用该旋转矩阵变换摄像机坐标系和IMU坐标系。对于参考脚型模型，顶点v_i的法向量由该顶点和周围k个邻近点共同确定。记顶点v_i的第j个邻近点为(j∈1,…,k)。模型中心点为c，则法向量同时满足：

其中Face{v_i}为包含点v_i的面片，几何意义为顶点v_i的法向量由以该点为顶点的面片的法向量的权重和，方向为由中心点指向顶点外侧。最终得到每个顶点和每个面片的法向量。

进一步地，所述步骤2具体为：对每帧视频帧进行Brief和Freak特征点检测和提取，对视频帧两两进行特征匹配，计算单应矩阵，并按纯视觉方法计算基本矩阵、本质矩阵，对本质矩阵进行SVD分解求得初始的摄像机姿态参数，同时恢复特征点对应的三维点云坐标。

进一步地，所述步骤3具体包括以下子步骤：

(3.1)对IMU数据做如下预处理。由于传感器频率和摄像机频率不一致，假设IMU数据和摄像机在从时间t_i到t_j的间隔时间Δt里同步，状态参数x_i＝[R_i,p_i,v_i,b_i]记录了IMU数据的旋转R_i、位置p_i、速度v_i和传感器漂移偏差b_i，而传感器漂移偏差b_i又包括速度加速计和重力加速计的偏差。则视频帧对应的状态参数为IMU状态参数的子集。设白噪声θ＝[θ^gd,θ^ad]，θ^gd和θ^ad分别为重力加速计和速度加速计的噪声，角速度ω_i，加速度a_i，则有：

状态参数求导后仍然为偏差b的函数，则可以采用分两步进行优化的方法，第一步假设b已知，对状态参数求梯度最优化，第二步中考虑到传感器移动时b有变化，为减少计算，当b变成b+Δb时，梯度值只需要在原来的基础上增加由b变化带来的差值Δb即可。最后计算残差，得到IMU数据下的摄像机姿态。

(3.2)结合IMU数据和视觉数据建立因子图。其中因子节点包括由视觉方法计算得到的摄像机姿态x_i′，恢复的脚型三维点云坐标l_i、IMU状态参数x_i，偏差c_i等。记第k时刻IMU测量值为则因子图中各个约束可以表示为：

IMU因子：

偏差因子：

f^bias(c_k+1,c_k)＝d(c_k+1-g(c_k))

先概率因子：

视觉投影因子：

f^Proj(x,l)＝d(z-π(x,l))

其中z表示真实的视频投影坐标，V为所有顶点的集合，π为投影函数。

立体视觉因子：

f^Stereo(x,l)＝d(z^R-π^R(x,l))d(z^L-π^L(x,l))

其中L和R分别对应左、右摄像头。

闭环检测因子：

f^Close＝sig(s_ij)*f(x_i,u_ij)-x_j

其中sig()为sigmoid函数，s_ij记录传感器运动方向(正向、逆向)，u_ij记录摄像机运动方向(正向、逆向)，而f(x_i,u_ij)为运动约束模型，可以看成是f^IMU。

(3.3)把因子图中因子函数分解成一般形式并进行优化，得到优化后的各个因子的值。

进一步地，所述步骤4具体包括以下子步骤：

(4.1)根据轮廓约束从重建出的三维脚型点云中提取出脚型的三维点云数据，并补脚底点集，得到完整的脚型点云模型。

(4.2)从已有的脚型数据库中选择与重建脚型相似的脚型参考模型，先以坐标轴和包围盒为参考对重建模型和参考模型进行初步对齐，然后去掉增加的脚底点云及脚踝以上的点避免干扰，只保留重建点云模型中脚面上的点集，利用ICP算法进一步与参考模型对齐，最后重新加上脚底点云，并增加手工调整，确保重建脚型模型与参考脚型模型对齐。

(4.3)对重建脚型模型进行去噪和光滑处理，包括剔除重建脚型模型和参考脚型模型之间对应点对距离大于阈值的点，以及等密度去掉部分重建脚型点云。

(4.4)重建脚型模型的法向量计算和表面重建。先用邻近点计算法向量的初值，然后对点云中的点进行可见性判断，并根据输入点的八叉树选取空间中包含点云的单位立方体，对每个立方体中的点云，选取以目标点为圆心，半径h_i范围内的点进行平面拟合：

其中为单位法向量，x为范围内的邻近点，D为给定的范围阈值。最后以立方体中的点组成圈，判断每个面片的法向量夹角是否小于π/2，剔除大于该角度的面片，进行表面重建。

进一步地，所述步骤5具体包括以下子步骤：

(5.1)在参考脚型模型选择的区域中任意选定曲率线L，在L上取样本点：

L_s＝{(l_i,t_i),i＝1,…,numOfSample}

其中l_i,t_i分别表示该像素点的二维视频坐标和切线向量，numOfSample为样本点数。然后把这些样本点映射到参考脚型的曲面上得到：

其中l_i′,t_i′分别表示该投影点的三维空间坐标及沿L方向上的切线平面，而表示该样本点的法向量。把这些信息向曲面周围延伸，记x为该样本点的邻近点，即在与法向量和切平面都垂直的方向上移动一小步Δ＝{Δx,Δy}，定义一个新网格点表示原来的点移动Δ后所在的点。在曲率线的两边，即沿b_i和-b_i方向都进行这样的操作，直到生成一系列网格点覆盖选择区域，对这些网格点进行插值直到选择的视频区域中每个像素点h＝(i,j)都对应一个三维网格节点。然后建立目标能量函数，计算并输出深度映射值集合Z＝{z_h}，对目标能量函数运用梯度下降法输出区域的深度映射视频，在重建模型上选择对应的区域根据曲率进行修复。目标能量函数共包括四项,为：

第一项表示曲面光滑程度,其中z_xx(i,j)表示在像素点h＝(i,j)沿x方向的深度，同理，z_xy(i,j)和z_yy(i,j)分别为沿xy斜线方向和y方向的深度。

第二项F(z_h)满足x|F(x)＝0，表示模型曲面是函数F的零集。

第三项R(z_h)为同一像素在两幅不同的视频帧中的深度差的平方值，保证在重叠区域两个深度映射值的一致性。

第四项C(z_h)＝∑_i,j(U^THU)²为曲率线的变化，其中H为z_h的Hessian矩阵，U为曲率线上的点投影到两帧视频帧上对应的像素。

(5.2)对重建脚型模型进行整体变换操作。最终要求的目标模型S‘可以表示为由参考脚型模型S^*和重建脚型模型S中的点在变换θ和移动δ＝{d_i}的作用下得到的。即S‘＝T(S^*,S)+d。考虑点集的GMM概率密度函数：

gmm(S)＝∑p(x)＝1/(2πσ²)exp(-||x-μ||/(2σ²))

p(x)为邻近点x的概率，μ为均值，σ为方差。

代价函数设为两个模型之间的L2距离

d_L2＝(gmm(S^*)-gmm(T(S,θ)))²

其中T为对点云进行变形操作，利用梯度下降法最小化目标函数得到模型间的形变系数θ。

(5.3)对重建脚型模型进行局部调整。定义形变目标函数和几何体目标函数Ω(S^*,S)分别为：

其中t_j为点p_j的移动向量，N(j)为点p_j的邻域，e_jk为连接和的边，A_jk为包含j和k两点的面片的面积，A_j为只包含点j的面片的面积。整体目标函数定义为通过对模型中点的拉伸移动改变模型。为保证整个模型的光滑度，考虑光滑目标模型S‘＝T(S^*,S,θ)+d，两点集间每对对应点之间的移动d_i可以通过最小化光滑目标函数求得。光滑目标函数E_smooth如下：

其中权重系数ε由模型数据集确定，N为对应模型上点的邻域。最后，计算点云的PCA模型矩阵M，m是参数向量，定义形变代价函数：

其中θ和t分别为旋转和平移参数,s为缩放因子，R为旋转矩阵。对用户选择的每部分区域进行形变，最终得到调整后的重建脚型模型。

本发明方法的核心是根据输入的裸脚RGB视频和IMU数据求得视频帧对应的摄像机参数和三维脚型点云模型，并结合数据库中的参考模型与重建模型建立对应关系，定义能量函数并最小化之，从而优化重建模型。

以下利用一个我们实现的例子来描述具体流程的实施方式，步骤如下(见图1)：

1)用带摄像机和传感器的智能手机绕脚360度拍摄视频，同时打开IMU数据采集app记录IMU数据(见图2(a))；

2)对数据进行各项预处理。

3)加入IMU数据，构建因子图(见图3)并优化；

4)把重建得到的模型与数据库中的参考模型(见图2(b))进行对齐、形变等一系列操作，得到最终的脚型网格模型(见图4)。

Claims (6)

1.一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，该方法包括以下步骤：

(1)利用智能手机围绕裸脚旋转一圈采集脚型视频和传感器数据，将视频分解成视频帧，并取部分帧作脚型区域分割等预处理，同时对脚型数据库中的参考脚型模型求法向量等参数；

(2)根据多视角立体几何的约束关系，检测视频帧中的特征点并进行匹配，计算摄像机内外参数，恢复特征点对应的三维点云坐标；

(3)加入IMU数据，对步骤2得到的摄像机参数和三维点云坐标进行修正，包括修复多视角算法失败无法求得摄像机参数的情况和移动三维点云坐标以减小三维坐标计算误差；

(4)从已有的脚型数据库中选择相似的参考脚型模型，对齐重建脚型模型和参考脚型模型，并对重建脚型模型进行去噪、光滑、表面重建处理；

(5)根据重建脚型模型和参考脚型模型之间的局部和整体对齐关系，对参考脚型模型进行形变使其拟合重建脚型模型，输出最终的重建脚型模型结果。

2.根据权利要求1所述的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，所述步骤1具体包括以下子步骤：

(1.1)用带有陀螺仪和加速计的智能手机绕脚360度拍摄视频一圈回到原点，拍摄时在脚底放置一张A4纸，拍摄过程中记录陀螺仪和加速计的数值。拍摄场景中的物体尺寸与真实尺寸之间的比例以A4纸作为参考。

(1.2)将采集到的视频分解成视频帧，从视频帧中等间隔地选择部分帧进行分割前背景区域，得到视频帧对应的关于脚型区域的mask视频。通过检测视频帧中A4纸的角点和边获取A4纸区域，根据RGB像素值分割A4纸及在A4纸区域内的部分脚型区域。同时，对整幅视频帧进行边缘检测，得到包括脚型边缘在内的视频帧边缘视频。然后由A4纸区域内的脚型区域向外搜索填充，最终得到初步估计的脚型区域及非脚型区域。建立GMM模型，通过k-mean算法分别对脚型区域像素和非脚型区域像素进行聚类，得到GMM模型中每个高斯模型的像素样本集，然后根据RGB值计算每个高斯模型的参数均值和协方差，每个高斯模型的权值定义为属于该高斯模型的像素个数与总像素个数的比值。对每个像素分配GMM模型中的高斯模型，以给定的脚型区域和非脚型区域为参考学习优化GMM模型，迭代分割，最终得到对应视频帧的脚型轮廓视频。

(1.3)用Kalibr工具箱对IMU进行标定，求出IMU到摄像机的旋转矩阵，利用该旋转矩阵变换摄像机坐标系和IMU坐标系。对于参考脚型模型，顶点v_i的法向量由该顶点和周围k个邻近点共同确定。记顶点v_i的第j个邻近点为(j∈1,…,k)。模型中心点为c，则法向量同时满足：

其中Face{v_i}为包含点v_i的面片，几何意义为顶点v_i的法向量由以该点为顶点的面片的法向量的权重和，方向为由中心点指向顶点外侧。最终得到每个顶点和每个面片的法向量。

3.根据权利要求1所述的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，所述步骤2具体为：对每帧视频帧进行Brief和Freak特征点检测和提取，对视频帧两两进行特征匹配，计算单应矩阵，并按纯视觉方法计算基本矩阵、本质矩阵，对本质矩阵进行SVD分解求得初始的摄像机姿态参数，同时恢复特征点对应的三维点云坐标。

4.根据权利要求1所述的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，所述步骤3具体包括以下子步骤：

(3.1)对IMU数据做如下预处理。由于传感器频率和摄像机频率不一致，假设IMU数据和摄像机在从时间t_i到t_j的间隔时间Δt里同步，状态参数x_i＝[R_i,p_i,v_i,b_i]记录了IMU数据的旋转R_i、位置p_i、速度v_i和传感器漂移偏差b_i，而传感器漂移偏差b_i又包括速度加速计和重力加速计的偏差。则视频帧对应的状态参数为IMU状态参数的子集。设白噪声 和分别为重力加速计和速度加速计的噪声，角速度ω_i，加速度a_i，则有：

状态参数求导后仍然为偏差b的函数，则可以采用分两步进行优化的方法，第一步假设b已知，对状态参数求梯度最优化，第二步中考虑到传感器移动时b有变化，为减少计算，当b变成bΔb时，梯度值只需要在原来的基础上增加由b变化带来的差值Δb即可。最后计算残差，得到IMU数据下的摄像机姿态。

(3.2)结合IMU数据和视觉数据建立因子图。其中因子节点包括由视觉方法计算得到的摄像机姿态x′_i，恢复的脚型三维点云坐标l_i、IMU状态参数x_i，偏差c_i等。记第k时刻IMU测量值为则因子图中各个约束可以表示为：

IMU因子：

$f^{IMU}(x_{k+1},x_k,c_k)=d(x_{k+1}-h(x_k,c_k,z_k^{IMU}\left)\right)$

偏差因子：

f^bias(C_k+1,c_k)＝d(c_k+1-g(c_k))

先概率因子：

$f^{\mathrm{Pr}ior}\left(v\right)=d\left(v\right),\∀v\∈V$

视觉投影因子：

f^Proj(x,l)＝d(z-π(x,l))

其中z表示真实的视频投影坐标，V为所有顶点的集合，π为投影函数。

立体视觉因子：

f^Stereo(x,l)＝d(z^R-π^R(x,l))d(z^L-π^L(x,l))

其中L和R分别对应左、右摄像头。

闭环检测因子：

f^Close＝sig(s_ij)*f(x_i,u_ij)-x_j

其中sig()为sigmoid函数，s_ij记录传感器运动方向(正向、逆向)，u_ij记录摄像机运动方向(正向、逆向)，而f(x_i,u_ij)为运动约束模型，可以看成是f^IMU。

(3.3)把因子图中因子函数分解成一般形式并进行优化，得到优化后的各个因子的值。

5.根据权利要求1所述的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，所述步骤4具体包括以下子步骤：

(4.1)根据轮廓约束从重建出的三维脚型点云中提取出脚型的三维点云数据，并补脚底点集，得到完整的脚型点云模型。

(4.2)从已有的脚型数据库中选择与重建脚型相似的脚型参考模型，先以坐标轴和包围盒为参考对重建模型和参考模型进行初步对齐，然后去掉增加的脚底点云及脚踝以上的点避免干扰，只保留重建点云模型中脚面上的点集，利用ICP算法进一步与参考模型对齐，最后重新加上脚底点云，并增加手工调整，确保重建脚型模型与参考脚型模型对齐。

(4.3)对重建脚型模型进行去噪和光滑处理，包括剔除重建脚型模型和参考脚型模型之间对应点对距离大于阈值的点，以及等密度去掉部分重建脚型点云。

(4.4)重建脚型模型的法向量计算和表面重建。先用邻近点计算法向量的初值，然后对点云中的点进行可见性判断，并根据输入点的八叉树选取空间中包含点云的单位立方体，对每个立方体中的点云，选取以目标点为圆心，半径h_i范围内的点进行平面拟合：

$H_i=\{<\stackrel{\&RightArrow;}{n_i},x>-D=0,x\&Element;R^3,|\stackrel{\&RightArrow;}{n_i}|=1\}$

其中为单位法向量，x为范围内的邻近点，D为给定的范围阈值。最后以立方体中的点组成圈，判断每个面片的法向量夹角是否小于π/2，剔除大于该角度的面片，进行表面重建。

6.根据权利要求1所述的一种通过智能手机采集脚型视频和传感器数据获取三维脚型的方法，其特征在于，所述步骤5具体包括以下子步骤：

(5.1)在参考脚型模型选择的区域中任意选定曲率线L，在L上取样本点：

L_s＝{(l_i,t_i),i＝1,…,numOfSample}

其中l_i,t_i分别表示该像素点的二维视频坐标和切线向量，numOfSample为样本点数。然后把这些样本点映射到参考脚型的曲面上得到：

$L_s^{\&prime;}=\{(l_i^{\&prime;},\stackrel{\&RightArrow;}{n_i},t_i^{\&prime;}),i=1,...,numOfSample\}$

其中l′_i,t′_i分别表示该投影点的三维空间坐标及沿L方向上的切线平面，而表示该样本点的法向量。把这些信息向曲面周围延伸，记x为该样本点的邻近点，即在与法向量和切平面都垂直的方向上移动一小步Δ＝{Δx,Δy}，定义一个新网格点表示原来的点移动Δ后所在的点。在曲率线的两边，即沿b_i和-b_i方向都进行这样的操作，直到生成一系列网格点覆盖选择区域，对这些网格点进行插值直到选择的视频区域中每个像素点h＝(i,j)都对应一个三维网格节点。然后建立目标能量函数，计算并输出深度映射值集合Z＝{z_h}，对目标能量函数运用梯度下降法输出区域的深度映射视频，在重建模型上选择对应的区域根据曲率进行修复。目标能量函数共包括四项,为：

$E\left(Z\right)=\underset{h}{\mathrm{\&Sigma;}}S\left(z_h\right)+F\left(z_h\right)+R\left(z_h\right)+C\left(z_h\right)$

第一项表示曲面光滑程度,其中z_xx(i,j)表示在像素点h＝(i,j)沿x方向的深度，同理，z_xy(i,j)和z_yy(i,j)分别为沿xy斜线方向和y方向的深度。

第二项F(z_h)满足x|F(x)＝0，表示模型曲面是函数F的零集。

第三项R(z_h)为同一像素在两幅不同的视频帧中的深度差的平方值，保证在重叠区域两个深度映射值的一致性。

第四项C(z_h)＝∑_i,j(U^THU)²为曲率线的变化，其中H为z_h的Hessian矩阵，U为曲率线上的点投影到两帧视频帧上对应的像素。

(5.2)对重建脚型模型进行整体变换操作。最终要求的目标模型S‘可以表示为由参考脚型模型S^*和重建脚型模型S中的点在变换θ和移动δ＝{d_i}的作用下得到的。即S‘＝T(S^*,S)+d。考虑点集的GMM概率密度函数：

gmm(S)＝Σp(x)＝1/(2πσ²)exp(-||x-μ||/2σ²))

p(x)为邻近点x的概率，μ为均值，σ为方差。

代价函数设为两个模型之间的L2距离

d_L2＝(gmm(S^*)-gmm(T(S,θ)))²

其中T为对点云进行变形操作，利用梯度下降法最小化目标函数得到模型间的形变系数θ。

(5.3)对重建脚型模型进行局部调整。定义形变目标函数和几何体目标函数Ω(S^*,S)分别为：

$\mathrm{\&Omega;}(S^{*},S)=\underset{j\&Element;S}{\mathrm{\&Sigma;}}{\mathrm{\&omega;}}_j{A}_j\left|\right|p_j+t_j-p_j^{*}\left|\right|$

其中t_j为点p_j的移动向量，N(j)为点p_j的邻域，e_jk为连接和的边，A_jk为包含j和k两点的面片的面积，A_j为只包含点j的面片的面积。整体目标函数定义为通过对模型中点的拉伸移动改变模型。为保证整个模型的光滑度，考虑光滑目标模型S‘＝T(S^*,S,θ)+d，两点集间每对对应点之间的移动d_i可以通过最小化光滑目标函数求得。光滑目标函数E_smooth如下：

$\begin{array}{c}\min E_{smooth}=\underset{p_i\&Element;N\left(p_i\right)}{\mathrm{\&Sigma;}}\mathrm{\&epsiv;}\left(d_i-{\left(p_i^{*}-p_i\right)}^2\right)+\\ \underset{p_i\&Element;N\left(p_i\right)}{\mathrm{\&Sigma;}}\underset{p_j\&Element;N\left(p_j\right)}{\mathrm{\&Sigma;}}{\left(d_i-d_j\right)}^2+\underset{p\&Element;N\left(p\right)}{\mathrm{\&Sigma;}}\left|\right|n_p-n_p^{*}\left|\right|\end{array}$

其中权重系数ε由模型数据集确定，N为对应模型上点的邻域。最后，计算点云的PCA模型矩阵M，m是参数向量，定义形变代价函数：

$e(S,\mathrm{\&theta;},m,t)=\mathrm{\&Sigma;}\left|\right|s*R(Mm+p_i^{*}+t-p_i)\left|\right|$

其中θ和t分别为旋转和平移参数,为缩放因子，R为旋转矩阵。对用户选择的每部分区域进行形变，最终得到调整后的重建脚型模型。