CN106157321B - 基于平面表面高动态范围图像的真实点光源位置测算方法 - Google Patents

Abstract

本发明涉及一种基于平面表面高动态范围图像的真实点光源位置测算方法，其特征在于：采用彩色数码相机与三角架相连，彩色数码相机与计算机相连，计算机中存有若干即62个圆，并将其编号dt记为1,2,3,...,62的半径为10mm大正圆随机分布，做为随机点云标定模板，将随机点云标定模板印刷于210mm*297mm纯白色表面光滑不透光pvc板上，记作随机点云标记板；其具有高效、简便等优点，更适合在实际用使用；其高效、快捷地估计出实际空间中点光源的位置并且为后期在三维建模软件中将其虚拟化做准备。

Description

基于平面表面高动态范围图像的真实点光源位置测算方法

技术领域

本发明涉及一种基于平面表面高动态范围图像的真实点光源位置测算方法，属于增强现实技术领域。

背景技术

目前，为了加强增强现实中虚拟对象的沉浸感，大多数的光照环境重构技术基本都是对测算设备或者场景有较高要求，并且为了减少对其他应用的影响，光照环境重构操作过程相对独立，使得整个增强现实系统的完成过程显得十分繁琐。

基于增强现实的立体电影制作过程中，决定真实场景、角色同虚拟场景、角色合成效果真实感的一个关键技术是光照一致性技术，包括明暗、反射、阴影等方面。早期基于增强现实的立体电影制作过程中，光照效果的解决方案多是虚拟布光和后期调整，这种方法只能提供一种近似的逼真感，并且这种方法需要花费大量的后期制作时间和精力，如果更换了环境，则需要对场景进行重新布光和调整。

基于图像的光照，其核心技术是记录真实场景光照信息载体的获取。经典的球镜模型是在2002年Kanbara等人将一个二维方形板和一个球镜组合在一起，形成一个三维标志物，能够同时完成几何注册和光照注册，以及更适用于复杂场景的改进方法，将随机点云标定技术应用其中，从而增强了该模型的鲁棒性和估测范围；美国特拉华大学的Wei Zhou等人，针对于以往光照估计中要求单一特定类型光源的局限，提出了一种场景光照模型的统一化框架，对不同类型的光源建立了更普遍适用的光源模型；Willianms提出的阴影映射算法，首先提取从光源到场景物体所有轮廓边，然后得到由光源点及每个轮廓边生成的半无限四边形，这些四边形构成物体的阴影体。

目前真实光照恢复方法主要有三种：第一种是在已知物体几何形状以及物体表面为Lambert漫反射的前提下，从单幅图像中根据物体表面明暗信息恢复出光源的数目和方向；第二种是通过已知几何形状的物体在场景中投射出的阴影区域，计算区域内的光能辐射分布来恢复出场景中的光照分布；第三种是为了模拟真实世界中的复杂光照，利用纯镜面反光球计算场景中某一点处的环境映照图，然后通过高动态范围的恢复算法计算场景的辐射度图。当获得了光照信息之后，其中虚拟光生成方法分为：一种是先用经验光照明模型和场景几何来恢复物体材料的光照属性，再根据这些属性生成新光照条件或新视点下的同一场景的图像。该方法可以实现漫游和交互，但比较复杂；另一种是基于不同视点、光照方向的采样图像，采用合适的基函数对采样图像进行插值、拟合等处理，从而获得以视点和光照方向为变量的重光照函数，包含了间接光照和环境光照，具有全局性，但全局光照计算复杂度较高。

发明内容

本发明的目的在于提供一种基于平面表面高动态范围图像的真实点光源位置测算方法，测算出真实点光源位置的测算模型，具有高效、简便等优点，更适合在实际用使用；其高效、快捷地估计出实际空间中点光源的位置并且为后期在三维建模软件中将其虚拟化做准备。

本发明的技术方案是这样实现的：一种基于平面表面高动态范围图像的真实点光源位置测算方法，其特征在于：采用彩色数码相机与三角架相连，彩色数码相机与计算机相连，计算机中存有62个圆，并将其编号dt记为1,2,3,...62,这62的半径为10mm大正圆随机分布，做为随机点云标定模板，将随机点云标定模板印刷于210mm*297mm纯白色表面光滑不透光pvc板上，记作随机点云标记板；

具体的测算步骤如下：

步骤1、将三角架调至适当高度，并调整三脚架的云台部分，使得通过彩色数码相机具体步骤如下的取景器目镜能够看到随机点云标记板的全部，将彩色数码相机调至M档，拍摄五张仅曝光时间不同的标记板照片，分别记为Pic0、Pic1、Pic2、Pic3、Pic4，对应的曝光时间依次记为Ev0、Ev1、Ev2、Ev3、Ev4，五张标记板照片的图像分辨率均为pix1*pix2；

步骤2、使用MALAB R2012b中的makehdr函数和Ev0、Ev1、Ev2、Ev3、Ev4，将Pic0、Pic1、Pic2、Pic3、Pic4合成为一张随机点云标记板高动态范围图像，记为HDRI，存储于计算机中；

步骤3、在随机点云标记板上建立以随机点云标记板的几何中心为原点O、沿纸面向上为Z轴正方向的右手坐标系，记为世界坐标系；

步骤4、使用MALAB R2012b将HDRI上所有像素点的亮度值存储到pix1*pix2的数组中，记为lum，其中lum(i,j)表示HDRI上坐标为(i,j)像素点的亮度值；在HDRI中随机点云标记板区域内的白色部分任意选取五个像素点，依次分别记为point0、point1、point2、point3、point4，并记它们在HDRI中的像素坐标依次分别为

(x₀,y₀)、(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)，同时，易获取point0、point1、point2、point3、point4所对应的亮度值，依次分别记为

lum(x₀,y₀)、lum(x₁,y₁)、lum(x₂,y₂)、lum(x₃,y₃)、lum(x₄,y₄)；

步骤5、使用OpenCV1.0中的cvThreshold函数，对随机点云标定模板进行二值化；使用OpenCV1.0中的cvFindContours函数，在二值化后的HDRI中找出各圆形的轮廓的点集；使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个圆形的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFindContours函数，对每个一维数组进行圆拟合，得到每个圆的圆心坐标，记为

(O_{dt_x},O_{dt_y},1)，

并从rotatedrect数据结构中提取出每个圆的二次曲线方程记为：

Q_dt(O_{dt_x},O_{dt_y})＝A_q_dt*O_{dt_x} ²+B_q_dt*O_{dt_x}*O_{dt_y}+C_q_dt*O_{dt_y} ²+D_q_dt*O_{dt_x}+E_q_dt*O_{dt_y}+F_q_dt＝0

，其中A_q_dt、B_q_dt、C_q_dt、D_q_dt、E_q_dt、F_q_dt分别为编号为dt的圆二次曲线方程的多项式系数，令

并将P_q_dt归一化记为P_qe_dt。

步骤6、将随机点云标定模板中的每一个圆的P_qe_dt后文称做“中心圆”与其最近的4个圆的P_qe_dt后文称做“邻近圆”构成一个按照中心圆编号为次序的PN*5矩阵，记为YU，其中YU[dt][0]记录中心圆的P_qe_dt，YU[dt][ii]记录其最近4个圆的P_qe_dt，ii＝1,2,3,4，如图4所示为一组，计算出每个中心圆曲线与其邻近圆曲线间的PP_{dt_ii}做为边长，其中

PP_{dt_ii}＝Tr(YU[dt][ii]/YU[dt][0])；

将4个邻近圆P_qe_dt编号为邻近圆曲线1、邻近圆曲线2、邻近圆曲线3、邻近圆曲线4。定义Tr(YU[dt][2]/YU[dt][1])为邻近圆曲线1与邻近圆曲线2的边长，定义Tr(YU[dt][3]/YU[dt][2])为邻近圆曲线2与邻近圆曲线3的边长，定义Tr(YU[dt][4]/YU[dt][3])为邻近圆曲线3与邻近圆曲线4的边长，定义Tr(YU[dt][1]/YU[dt][4])为邻近圆曲线4与邻近圆曲线1的边长。

根据余弦定理计算方法如下：

将每个中心圆的四个余弦值存储为62*4的数组，记为YXJ1，其中YXJ1(dt,0)、YXJ1(dt,1)、YXJ1(dt,2)、YXJ1(dt,3)表示随机点云标定模板中编号为dt的圆的4个相关余弦值，记为

[cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]；

步骤7、使用OpenCV1.0中的cvThreshold函数，对HDRI进行二值化；使用OpenCV1.0中的cvFindContours函数，找出每一个可见椭圆的轮廓点集；随后，需使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个可见椭圆的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFitEllipse函数，对每个一维数组进行椭圆拟合，得到每个椭圆的圆心在世界坐标系下的坐标，记为

(T_{dt′_x},T_{dt′_y},1)，

并从rotatedrect数据结构中提取出每个可见椭圆的二次曲线方程记为：

Q′_dt′(T_{dt′_x},T_{dt′_y})＝A′_q_dt′*T_{dt′_x} ²+B′_q_dt′*T_{dt′_x}*T_{dt′_y}+C′_q_dt′*T_{dt′_y} ²+D′_q_dt′*T_{dt′_x}+E′_q_dt′*T_{dt′_y}+F′_q_dt′＝0

其中A'_q_dt′、B'_q_dt′、C'_q_dt′、D'_q_dt′、E'_q_dt′、F'_q_dt′分别为编号为dt′的椭圆二次曲线方程的多项式系数，令

并将P′_q_dt′归一化记为P′_qe_dt′，其中dt′表示HDRI中可见椭圆的编号；

步骤8、将HDRI中的每一个可见的椭圆的P′_qe_dt′(后文称做“中心椭圆”)与其最近的4个椭圆的P′_qe_dt′(后文称做“邻近椭圆”)构成一个按照中心椭圆编号为次序的TN*5矩阵，记为TU，其中TU[dt’][0]记录中心椭圆的P′_qe_dt′，TU[dt’][ii]记录其最近4椭圆的P′_qe_dt′，ii＝1,2,3,4，如图5所示为一组，计算出每个中心椭圆与其邻近椭圆间的PP′_{dt′_ii}做为边长，其中

PP′_{dt′_ii}＝Tr(TU[dt′][ii]/TU[dt′][0])；

将4个邻近椭圆P′_qe_dt′编号为邻近椭圆曲线1、邻近椭圆曲线2、邻近椭圆曲线3、邻近椭圆曲线4。定义Tr(TU[dt'][2]/TU[dt'][1])为邻近椭圆曲线1与邻近椭圆曲线2的边长，定义Tr(TU[dt'][3]/TU[dt'][2])为邻近椭圆曲线2与邻近椭圆曲线3的边长，定义Tr(TU[dt'][4]/TU[dt'][3])为邻近椭圆曲线3与邻近椭圆曲线4的边长，定义Tr(TU[dt'][1]/TU[dt'][4])为邻近椭圆曲线4与邻近椭圆曲线1的边长；

根据余弦定理计算方法如下：

将每个中心椭圆处的四个余弦值存储为TN*4的数组，记为YXJ2[TN][4]，其中YXJ2[dt’]表示摄影球图片中可见的编号为dt’的椭圆的4个相关余弦值可记为[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}]

，TN表示在HDRI中观察到的椭圆圆点的总数；

步骤9、将数组YXJ2与YXJ1进行比较，找出其中相匹配的点，其满足：|[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}][cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]^-1|≤0.01

按照YXJ1中各点顺序将YXJ2中与YXJ1中相匹配的点存储于TN*4的数组中，记为YXJ3。

步骤10、将YXJ1中每个圆形的圆心的坐标(O_{dt_x},O_{dt_y},1)和YXJ3中对应的椭圆圆心的坐标(T_{dt′_x},T_{dt′_y},T_{dt′_z})传递给OpenCV1.0中的摄像机标定函数CalibrateCamera2，便可得彩色数码相机2的3*3旋转矩阵R和3*1平移矩阵T，从而可以计算得到point0、point1、point2、point3、point4在世界坐标系中的坐标(W_{N_1},W_{N_2},0)，N＝0,1,2,3,4：

令

由(x_N,y_N,0,1)^T＝A*(W_{N_1},W_{N_2},0,1)^T

解出采样点point0、point1、point2、point3、point4在世界坐标系中的坐标。

步骤11、将point0、point1、point2、point3、point4中的任意四个点为一组，分为不重复的五组，并将这五组采样点所对应的在世界坐标系中的坐标和亮度值存储到一个5*4*4数组中，记为g_pi，其中g_pi[n][p][q]q＝0,1,2表示第n组采样点中第p个采样点在世界坐标系中的坐标，g_pi[n][p][3]表示第n组采样点中第p个采样点的亮度值；

步骤12、分别假设n＝0、1、2、3、4，对g_pi[n][p][3]执行步骤13、14、15操作；

步骤13、根据对实际点光源的位置，估计出大致测算范围M，使得欲测算的点光源能够位于以随机点云标记板所在平面为底面的M*M*M的立方体内，该立方体记为C。设该场景下环境光的亮度值为lum_amb，则可针对每组采样点利用牛顿迭代法解如下方程组：

设最终待求点光源在图2所示坐标系中的坐标为(L_x,L_y,L_z),

，该方程组记为“三球方程组”；

步骤14、每组采样点在解步骤13中涉及的方程组时所选取的初始值可采用如下方法确定：

，该式记为“最小距离式”。其中L_{x_f},L_{y_f},L_{z_f}为步骤13中牛顿迭代法解方程组时以x_f,y_f,z_f为初始值的估算结果。

步骤15、选取使得最小距离式成立的体中心点所在较小立方体作为新的待分割立方体，重复操作步骤14，直到最小距离式的结果值小于0.001或者待分割立方体的边长不大于0.01，重复步骤14终止，至此，求得每组采样点所估算出的关于lum_amb的点光源坐标，记为

(L′_{n_x},L′_{n_y},L′_{n_z})；

步骤16、以步骤13中涉及的lum_amb为自变量，求解出满足下式的lum_amb值：

步骤17、根据步骤16中所得的lum_amb值，依照步骤14、15的操作求解三球方程组，并将五个三球方程组的结果坐标取平均值。

通过以上步骤便可快速、高精度地估计实际空间中点光源的位置并且为后期在三维建模软件中将其虚拟化做准备。

本发明的积极效果是提出了一种基于平面表面高动态范围图像中的照度信息，测算出真实点光源位置的测算模型。该方法操作简单，对输入参数要求较低，仅需一副平面表面的高动态范围图像，获得点光源的空间位置；同时，为了增强鲁棒性，使用了点云标定，克服了传统的方形标定板无法标定大角度位置的缺陷；可以与各种增强现实的应用相结合，对沉浸感的增强起到了很大的帮助，对光照的重构提供了参考。

附图说明

图1是一种基于单幅图像中的点光源信息，生成虚拟点光源姿态的方法所需设备构成图，其中：1为A4随机点云标记板，2为彩色数码相机，3为三角架，4为计算机，此图是说明书摘要附图。

图2是随机点云标记板1坐标系图，其中：1为A4随机点云标记板，此图是说明书摘要附图。

图3为本方法中使用的A4随机点云标记图，此图是说明书摘要附图。

图4是随机点云标记板1中五个圆点组示意图，此图是说明书摘要附图。

图5是随机点云标记板1中五个圆点组透视示意图，此图是说明书摘要附图。

图6是估计立方体空间的分割示意图，此图是说明书摘要附图。

具体实施方式

下面结合附图对本发明做进一步的描述：如图1所示，一种基于平面表面高动态范围图像的真实点光源位置测算方法，其特征在于：采用彩色数码相机2与三角架3相连，彩色数码相机2与计算机4相连，计算机4中存有62个圆，并将其编号dt记为1,2,3,...62,这62的半径为10mm大正圆随机分布，做为随机点云标定模板，将随机点云标定模板印刷于210mm*297mm纯白色表面光滑不透光pvc板上，记作随机点云标记板1；

具体的测算步骤如下：

步骤1、将三角架3调至适当高度，并调整三脚架3的云台部分，使得通过彩色数码相机2具体步骤如下的取景器目镜能够看到随机点云标记板1的全部，将彩色数码相机2调至M档，拍摄五张仅曝光时间不同的标记板照片，分别记为Pic0、Pic1、Pic2、Pic3、Pic4，对应的曝光时间依次记为Ev0、Ev1、Ev2、Ev3、Ev4，五张标记板照片的图像分辨率均为pix1*pix2。

步骤2、使用MALAB R2012b中的makehdr函数和Ev0、Ev1、Ev2、Ev3、Ev4，将Pic0、Pic1、Pic2、Pic3、Pic4合成为一张随机点云标记板1高动态范围图像，记为HDRI，存储于计算机4中。

步骤3、在随机点云标记板1上建立以随机点云标记板1的几何中心为原点O、沿纸面向上为Z轴正方向的右手坐标系，记为世界坐标系。

步骤4、使用MALAB R2012b将HDRI上所有像素点的亮度值存储到pix1*pix2的数组中，记为lum，其中lum(i,j)表示HDRI上坐标为(i,j)像素点的亮度值。在HDRI中随机点云标记板1区域内的白色部分任意选取五个像素点，依次分别记为point0、point1、point2、point3、point4，并记它们在HDRI中的像素坐标依次分别为

(x₀,y₀)、(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)，

同时，易获取point0、point1、point2、point3、point4所对应的亮度值，依次分别记为

lum(x₀,y₀)、lum(x₁,y₁)、lum(x₂,y₂)、lum(x₃,y₃)、lum(x₄,y₄)。

步骤5、使用OpenCV1.0中的cvThreshold函数，对随机点云标定模板进行二值化；使用OpenCV1.0中的cvFindContours函数，在二值化后的HDRI中找出各圆形的轮廓的点集；使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个圆形的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFindContours函数，对每个一维数组进行圆拟合，得到每个圆的圆心坐标，记为

(O_{dt_x},O_{dt_y},1)，

并从rotatedrect数据结构中提取出每个圆的二次曲线方程记为：

Q_dt(O_{dt_x},O_{dt_y})＝A_q_dt*O_{dt_x} ²+B_q_dt*O_{dt_x}*O_{dt_y}+C_q_dt*O_{dt_y} ²+D_q_dt*O_{dt_x}+E_q_dt*O_{dt_y}+F_q_dt＝0

其中A_q_dt、B_q_dt、C_q_dt、D_q_dt、E_q_dt、F_q_dt分别为编号为dt的圆二次曲线方程的多项式系数，令

并将P_q_dt归一化记为P_qe_dt。

步骤6、将随机点云标定模板中的每一个圆的P_qe_dt后文称做“中心圆”与其最近的4个圆的P_qe_dt后文称做“邻近圆”构成一个按照中心圆编号为次序的PN*5矩阵，记为YU，其中YU[dt][0]记录中心圆的P_qe_dt，YU[dt][ii]记录其最近4个圆的P_qe_dt，ii＝1,2,3,4，如图4所示为一组，计算出每个中心圆曲线与其邻近圆曲线间的PP_{dt_ii}做为边长，其中

PP_{dt_ii}＝Tr(YU[dt][ii]/YU[dt][0])。

将4个邻近圆P_qe_dt编号为邻近圆曲线1、邻近圆曲线2、邻近圆曲线3、邻近圆曲线4。定义Tr(YU[dt][2]/YU[dt][1])为邻近圆曲线1与邻近圆曲线2的边长，定义Tr(YU[dt][3]/YU[dt][2])为邻近圆曲线2与邻近圆曲线3的边长，定义Tr(YU[dt][4]/YU[dt][3])为邻近圆曲线3与邻近圆曲线4的边长，定义Tr(YU[dt][1]/YU[dt][4])为邻近圆曲线4与邻近圆曲线1的边长。

根据余弦定理计算方法如下：

将每个中心圆的四个余弦值存储为62*4的数组，记为YXJ1，其中YXJ1(dt,0)、YXJ1(dt,1)、YXJ1(dt,2)、YXJ1(dt,3)表示随机点云标定模板中编号为dt的圆的4个相关余弦值，记为

[cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]。

步骤7、使用OpenCV1.0中的cvThreshold函数，对HDRI进行二值化；使用OpenCV1.0中的cvFindContours函数，找出每一个可见椭圆的轮廓点集；随后，需使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个可见椭圆的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFitEllipse函数，对每个一维数组进行椭圆拟合，得到每个椭圆的圆心在世界坐标系下的坐标，记为

(T_{dt′_x},T_{dt′_y},1)，

并从rotatedrect数据结构中提取出每个可见椭圆的二次曲线方程记为：

Q′_dt′(T_{dt′_x},T_{dt′_y})＝A′_q_dt′*T_{dt′_x} ²+B′_q_dt′*T_{dt′_x}*T_{dt′_y}+C′_q_dt′*T_{dt′_y} ²+D′_q_dt′*T_{dt′_x}+E′_q_dt′*T_{dt′_y}+F′_q_dt′＝0

其中A'_q_dt′、B'_q_dt′、C'_q_dt′、D'_q_dt′、E'_q_dt′、F'_q_dt′分别为编号为dt′的椭圆二次曲线方程的多项式系数，令

并将P′_q_dt′归一化记为P′_qe_dt′，其中dt′表示HDRI中可见椭圆的编号。

步骤8、将HDRI中的每一个可见的椭圆的P′_qe_dt′(后文称做“中心椭圆”)与其最近的4个椭圆的P′_qe_dt′(后文称做“邻近椭圆”)构成一个按照中心椭圆编号为次序的TN*5矩阵，记为TU，其中TU[dt’][0]记录中心椭圆的P′_qe_dt′，TU[dt’][ii]记录其最近4椭圆的P′_qe_dt′，ii＝1,2,3,4，如图5所示为一组，计算出每个中心椭圆与其邻近椭圆间的PP′_{dt′_ii}做为边长，其中

PP′_{dt′_ii}＝Tr(TU[dt′][ii]/TU[dt′][0])。

将4个邻近椭圆P′_qe_dt′编号为邻近椭圆曲线1、邻近椭圆曲线2、邻近椭圆曲线3、邻近椭圆曲线4。定义Tr(TU[dt'][2]/TU[dt'][1])为邻近椭圆曲线1与邻近椭圆曲线2的边长，定义Tr(TU[dt'][3]/TU[dt'][2])为邻近椭圆曲线2与邻近椭圆曲线3的边长，定义Tr(TU[dt'][4]/TU[dt'][3])为邻近椭圆曲线3与邻近椭圆曲线4的边长，定义Tr(TU[dt'][1]/TU[dt'][4])为邻近椭圆曲线4与邻近椭圆曲线1的边长。

根据余弦定理计算方法如下：

将每个中心椭圆处的四个余弦值存储为TN*4的数组，记为YXJ2[TN][4]，其中YXJ2[dt’]表示摄影球图片中可见的编号为dt’的椭圆的4个相关余弦值可记为[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}]

，TN表示在HDRI中观察到的椭圆圆点的总数。

步骤9、将数组YXJ2与YXJ1进行比较，找出其中相匹配的点，其满足：|[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}][cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]^-1|≤0.01

按照YXJ1中各点顺序将YXJ2中与YXJ1中相匹配的点存储于TN*4的数组中，记为YXJ3。

步骤10、将YXJ1中每个圆形的圆心的坐标(O_{dt_x},O_{dt_y},1)和YXJ3中对应的椭圆圆心的坐标(T_{dt′_x},T_{dt′_y},T_{dt′_z})传递给OpenCV1.0中的摄像机标定函数CalibrateCamera2，便可得彩色数码相机2的3*3旋转矩阵R和3*1平移矩阵T，从而可以计算得到point0、point1、point2、point3、point4在世界坐标系中的坐标(W_{N_1},W_{N_2},0)，N＝0,1,2,3,4：

令

由(x_N,y_N,0,1)^T＝A*(W_{N_1},W_{N_2},0,1)^T

解出采样点point0、point1、point2、point3、point4在世界坐标系中的坐标。

步骤11、将point0、point1、point2、point3、point4中的任意四个点为一组，分为不重复的五组，并将这五组采样点所对应的在世界坐标系中的坐标和亮度值存储到一个5*4*4数组中，记为g_pi，其中g_pi[n][p][q]q＝0,1,2表示第n组采样点中第p个采样点在世界坐标系中的坐标，g_pi[n][p][3]表示第n组采样点中第p个采样点的亮度值。

步骤12、分别假设n＝0、1、2、3、4，对g_pi[n][p][3]执行步骤13、14、15操作。

步骤13、根据对实际点光源的位置，估计出大致测算范围M，使得欲测算的点光源能够位于以随机点云标记板1所在平面为底面的M*M*M的立方体内，该立方体记为C。设该场景下环境光的亮度值为lum_amb，则可针对每组采样点利用牛顿迭代法解如下方程组：

设最终待求点光源在图2所示坐标系中的坐标为(L_x,L_y,L_z),

，该方程组记为“三球方程组”。

步骤14、每组采样点在解步骤13中涉及的方程组时所选取的初始值可采用如下方法确定：

由M*M*M立方体的各边中点，将M*M*M的待分割立方体分割为等大的八个较小的立方体，按照方位顺序将这八个较小立方体编号，并记为C_f，其中f＝0,1,2,…,7，此八个较小立方体的体中心点在世界坐标系中的空间坐标记为(x_f,y_f,z_f)，并计算如下式子：

，该式记为“最小距离式”。其中L_{x_f},L_{y_f},L_{z_f}为步骤13中牛顿迭代法解方程组时以x_f,y_f,z_f为初始值的估算结果。

步骤15、选取使得最小距离式成立的体中心点所在较小立方体作为新的待分割立方体，重复操作步骤14，直到最小距离式的结果值小于0.001或者待分割立方体的边长不大于0.01，重复步骤14终止。至此，求得每组采样点所估算出的关于lum_amb的点光源坐标，记为

(L′_{n_x},L′_{n_y},L′_{n_z})。

步骤16、以步骤13中涉及的lum_amb为自变量，求解出满足下式的lum_amb值：

步骤17、根据步骤16中所得的lum_amb值，依照步骤14、15的操作求解三球方程组，并将五个三球方程组的结果坐标取平均值。

通过以上步骤便可快速、高精度地估计实际空间中点光源的位置并且为后期在三维建模软件中将其虚拟化做准备。

Claims (1)

1.一种基于平面表面高动态范围图像的真实点光源位置测算方法，其特征在于：采用彩色数码相机与三角架相连，彩色数码相机与计算机相连，计算机中存有62个圆，并将其编号dt记为1,2,3,...62,这62个半径为10mm大正圆随机分布在一个平面上，做为随机点云标定模板，将随机点云标定模板印刷于210mm*297mm纯白色表面光滑不透光pvc板上，记作随机点云标记板；

具体的测算步骤如下：

步骤1、将三角架调至适当高度，并调整三脚架的云台部分，使得通过彩色数码相机的取景器目镜能够看到随机点云标记板的全部，将彩色数码相机调至M档，拍摄五张仅曝光时间不同的标记板照片，分别记为Pic0、Pic1、Pic2、Pic3、Pic4，对应的曝光时间依次记为Ev0、Ev1、Ev2、Ev3、Ev4，五张标记板照片的图像分辨率均为pix1*pix2；

步骤2、使用MATLAB R2012b中的makehdr函数和Ev0、Ev1、Ev2、Ev3、Ev4，将Pic0、Pic1、Pic2、Pic3、Pic4合成为一张随机点云标记板高动态范围图像，记为HDRI，存储于计算机中；

步骤3、在随机点云标记板上建立以随机点云标记板的几何中心为原点O、沿纸面向上为Z轴正方向的右手坐标系，记为世界坐标系；

步骤4、使用MATLAB R2012b将HDRI上所有像素点的亮度值存储到pix1*pix2的数组中，记为lum，其中lum(i,j)表示HDRI上坐标为(i,j)像素点的亮度值，在HDRI中随机点云标记板区域内的白色部分任意选取五个像素点，依次分别记为point0、point1、point2、point3、point4，并记它们在HDRI中的像素坐标依次分别为(x₀,y₀)、(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)，

同时，易获取point0、point1、point2、point3、point4所对应的亮度值，依次分别记为lum(x₀,y₀)、lum(x₁,y₁)、lum(x₂,y₂)、lum(x₃,y₃)、lum(x₄,y₄)；

步骤5、使用OpenCV1.0中的cvThreshold函数，对随机点云标定模板图像进行二值化；使用OpenCV1.0中的cvFindContours函数，在二值化后的模板图像中找出各圆形的轮廓的点集；使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个圆形的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFitEllipse函数，对每个一维数组进行圆拟合，得到每个圆拟合的圆心坐标，记为(O_{dt_x},O_{dt_y},1)，

并从rotatedrect数据结构中提取出每个圆拟合的二次曲线方程记为：

Q_dt(O_{dt_x},O_{dt_y})＝A_q_dt*O_{dt_x} ²+B_q_dt*O_{dt_x}*O_{dt_y}+C_q_dt*O_{dt_y} ²+D_q_dt*O_{dt_x}+E_q_dt*O_{dt_y}+F_q_dt＝0，

其中A_q_dt、B_q_dt、C_q_dt、D_q_dt、E_q_dt、F_q_dt分别为编号为dt的圆二次曲线方程的多项式系数，令

并将P_q_dt归一化记为P_qe_dt；

步骤6、将随机点云标定模板中的每一个圆的P_qe_dt,后文称做“中心圆”与其最近的4个圆的P_qe_dt后文称做“邻近圆”构成一个按照中心圆编号为次序的PN*5矩阵，记为YU，其中YU[dt][0]记录中心圆的P_qe_dt，YU[dt][ii]记录中心圆的最近4个邻近圆的P_qe_dt，ii＝1,2,3,4,计算出每个中心圆曲线与其邻近圆曲线间的PP_{dt_ii}做为边长，其中

PP_{dt_ii}＝Tr(YU[dt][ii]/YU[dt][0])；

将4个邻近圆P_qe_dt编号为邻近圆曲线1、邻近圆曲线2、邻近圆曲线3、邻近圆曲线4；定义Tr(YU[dt][2]/YU[dt][1])为邻近圆曲线1与邻近圆曲线2的边长，定义Tr(YU[dt][3]/YU[dt][2])为邻近圆曲线2与邻近圆曲线3的边长，定义Tr(YU[dt][4]/YU[dt][3])为邻近圆曲线3与邻近圆曲线4的边长，定义Tr(YU[dt][1]/YU[dt][4])为邻近圆曲线4与邻近圆曲线1的边长；

根据余弦定理计算方法如下：

将每个中心圆的四个余弦值存储为62*4的数组，记为YXJ1，其中YXJ1(dt,0)、YXJ1(dt,1)、YXJ1(dt,2)、YXJ1(dt,3)表示随机点云标定模板中编号为dt的圆的4个相关余弦值，记为[cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]；

步骤7、使用OpenCV1.0中的cvThreshold函数，对HDRI进行二值化；使用OpenCV1.0中的cvFindContours函数，找出每一个可见椭圆的轮廓点集；随后，需使用OpenCV1.0中的cvCvtSeqToArray函数，依次将每一个可以看到的椭圆的轮廓的点集转换为一维数组，一维数组中元素的类型都为CvPoint；最后，使用OpenCV1.0中的cvFitEllipse函数，对每个一维数组进行椭圆拟合，得到每个椭圆的圆心在世界坐标系下的坐标，记为(T_{dt′_x},T_{dt′_y},1)，

并从rotatedrect数据结构中提取出每个可见椭圆的二次曲线方程记为：

Q′_dt′(T_{dt′_x},T_{dt′_y})＝A′_q_dt′*T_{dt′_x} ²+B′_q_dt′*T_{dt′_x}*T_{dt′_y}+C′_q_dt′*T_{dt′_y} ²+D′_q_dt′*T_{dt′_x}+E′_q_dt′*T_{dt′_y}+F′_q_dt′＝0

其中A'_q_dt′、B'_q_dt′、C'_q_dt′、D'_q_dt′、E'_q_dt′、F'_q_dt′分别为编号为dt′的椭圆二次曲线方程的多项式系数，令

并将P′_q_dt′归一化记为P′_qe_dt′，其中dt′表示HDRI中可见椭圆的编号；

步骤8、将HDRI中的每一个可见的椭圆的P′_qe_dt′后文称做“中心椭圆”与其最近的4个椭圆的P′_qe_dt′后文称做“邻近椭圆”构成一个按照中心椭圆编号为次序的TN*5矩阵，记为TU，其中TU[dt’][0]记录中心椭圆的P′_qe_dt′，TU[dt’][ii]记录中心椭圆的最近4个邻近椭圆的P′_qe_dt′，ii＝1,2,3,4,计算出每个中心椭圆与其邻近椭圆间的PP_d′_t′__ii做为边长，其中

PP′_dt′__ii＝Tr(TU[dt′][ii]/TU[dt′][0])；

将4个邻近椭圆P′_qe_dt′编号为邻近椭圆曲线1、邻近椭圆曲线2、邻近椭圆曲线3、邻近椭圆曲线4；定义Tr(TU[dt'][2]/TU[dt'][1])为邻近椭圆曲线1与邻近椭圆曲线2的边长，定义Tr(TU[dt'][3]/TU[dt'][2])为邻近椭圆曲线2与邻近椭圆曲线3的边长，定义Tr(TU[dt'][4]/TU[dt'][3])为邻近椭圆曲线3与邻近椭圆曲线4的边长，定义Tr(TU[dt'][1]/TU[dt'][4])为邻近椭圆曲线4与邻近椭圆曲线1的边长；

根据余弦定理计算方法如下：

将每个中心椭圆处的四个余弦值存储为TN*4的数组，记为YXJ2[TN][4]，TN表示在HDRI中观察到的椭圆圆点的总数；其中YXJ2[dt’]表示摄影球图片中可见的编号为dt’的椭圆的4个相关余弦值记为[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}]

步骤9、将数组YXJ2与YXJ1进行比较，找出其中相匹配的点，其满足：

|[cosα′_{dt′_0},cosα′_{dt′_1},cosα′_{dt′_2},cosα′_{dt′_3}][cosα_{dt_0},cosα_{dt_1},cosα_{dt_2},cosα_{dt_3}]^-1|≤0.01；

按照YXJ1中各点顺序将YXJ2中与YXJ1中相匹配的点存储于TN*4的数组中，记为YXJ3；

步骤10、将YXJ1中每个圆形的圆心的坐标(O_{dt_x},O_{dt_y},1)和YXJ3中对应的椭圆圆心的坐标(T_{dt′_x},T_{dt′_y},T_{dt′_z})传递给OpenCV1.0中的摄像机标定函数CalibrateCamera2，便可得彩色数码相机的3*3旋转矩阵R和3*1平移矩阵T，从而可以计算得到point0、point1、point2、point3、point4在世界坐标系中的坐标(W_{N_1},W_{N_2},0)，N＝0,1,2,3,4：

令

由(x_N,y_N,0,1)^T＝A*(W_{N_1},W_{N_2},0,1)^T解出采样点point0、point1、point2、point3、point4在世界坐标系中的坐标；

步骤11、将point0、point1、point2、point3、point4中的任意四个点为一组，分为不重复的五组，并将这五组采样点所对应的在世界坐标系中的坐标和亮度值存储到一个5*4*4数组中，记为g_pi，其中g_pi[n][p][q]q＝0,1,2表示第n组采样点中第p个采样点在世界坐标系中的坐标，g_pi[n][p][3]表示第n组采样点中第p个采样点的亮度值；

步骤12、分别假设n＝0、1、2、3、4，对g_pi[n][p][3]执行步骤13、14、15操作；

步骤13、根据对实际点光源的位置，估计出大致测算范围M，使得欲测算的点光源能够位于以随机点云标记板所在平面为底面的M*M*M的立方体内，该立方体记为C，设该场景下环境光的亮度值为lum_amb，则针对每组采样点利用牛顿迭代法解如下方程组：

设最终待求点光源在坐标系中的坐标为(L_x,L_y,L_z),

，

该方程组记为“三球方程组”；

步骤14、每组采样点在解步骤13中涉及的方程组时所选取的初始值可采用如下方法确定：

由M*M*M立方体的各边中点，将M*M*M的待分割立方体分割为等大的八个较小的立方体，按照方位顺序将这八个较小立方体编号，并记为C_f，其中f＝0,1,2,…,7，此八个较小立方体的体中心点在世界坐标系中的空间坐标记为(x_f,y_f,z_f)，并计算如下式子：

，

该式记为“最小距离式”，其中L_{x_f},L_{y_f},L_{z_f}为步骤13中牛顿迭代法解方程组时以x_f,y_f,z_f为初始值的估算结果；

步骤15、选取使得最小距离式成立的体中心点所在较小立方体作为新的待分割立方体，重复操作步骤14，直到最小距离式的结果值小于0.001或者待分割立方体的边长不大于0.01，重复步骤14终止，至此，求得每组采样点所估算出的关于lum_amb的点光源坐标，记为(L′_{n_x},L′_{n_y},L′_{n_z})；

步骤16、以步骤13中涉及的lum_amb为自变量，求解出满足下式的lum_amb值：

步骤17、根据步骤16中所得的lum_amb值，依照步骤14、15的操作求解三球方程组，并将五个三球方程组的结果坐标取平均值；

通过以上步骤1到步骤17便可快速、高精度地估计实际空间中点光源的位置并且为后期在三维建模软件中将其虚拟化做准备。

Images