CN107608077A

CN107608077A - 一种多光源位置估计方法

Info

Publication number: CN107608077A
Application number: CN201710805388.6A
Authority: CN
Inventors: 李华; 韩成; 张超; 权巍; 蒋振刚; 薛耀红; 胡汉平; 徐超
Original assignee: Changchun University of Science and Technology
Current assignee: Changchun University of Science and Technology
Priority date: 2017-09-08
Filing date: 2017-09-08
Publication date: 2018-01-19
Anticipated expiration: 2037-09-08
Also published as: CN107608077B

Abstract

本发明涉及一种多光源位置估计方法，其特征在于具体步骤如下：记世界坐标系与摄像头坐标系重合；计算机利用Matlab2015中的rgb2gray函数将srcimage转化成灰度图像，利用公式，构造方程组，根据实对称矩阵的正交对角化分解原理，设置一个阈值，通过手动方式，构造方程组，构造矩阵P，表示第i个光源在世界坐标系中的坐标，通过以上步骤实现了多光源位置的估计。其能解决增强现实中基于多光源的虚实融合场景的无缝融合等问题，从而，使得虚拟对象叠加到真实场景中更具沉浸感，更加逼真。

Description

一种多光源位置估计方法

技术领域

本发明涉及一种多光源位置估计方法，属于增强现实技术领域。

背景技术

增强现实技术是将虚拟对象无缝渲染在真实环境中，它已成功应用在教育、医学、娱乐等众多领域。光照一致性是增强现实中的关键技术之一，是指真实和虚拟对象具有一致性的照明。一旦获取光源参数，就能在虚实融合场景中操纵虚拟对象的明暗及产生阴影。而多光源位置估计是复杂环境光照一致性的重点和难点。

增强现实中的光照检测算法可以分为基于图像光照(image-based lighting,IBL)和逆向渲染(inverse rendering)两大类光源恢复算法。基于图像光照技术是HDR(high dynamic range，高动态范围)技术与EM(environment map，环境映照)技术的结合，它使用HDREM代替虚拟光源来“照明”虚拟物体。HDR图像通常是拍摄若干张不同曝光度的周围环境的全景图像，然后根据相机响应曲线合成HDR图像。HDR图像能够较精确的展示周围环境的真实光照信息，最后利用环境映照技术将光照信息叠加到虚拟对象上。但它对虚拟对象位置具有严格限制，且无法生成虚拟对象的阴影，有时还需使用鱼眼镜头或者全景相机来获取图像。逆向渲染是根据光度模型(photometric model)从图像中追踪物体反射属性和光源信息。一种典型的实现方法是通过镜面反射校准物的明暗，使用Phone光照模型检测真实环境中点光源位置或方向信息。基于逆向渲染算法恢复光源信息无需配备昂贵的设配，只需一个普通的CMOS摄像头，而且操作简单。因此，逆向渲染算法在光照估计中具有广泛的应用价值。

Mark W.Powell等人通过三角化两个镜面反射球上的高光点来获取多光源位置。他们假设镜面反射球在图像上的投影是一个正圆，基于这个假设对两个镜面反射球进行定位；根据Phone光照模型和逆向光线追踪算法追踪球面高光处的反射光线方向；然后匹配两个镜面反射球上来自同一光源照明的高光点，高光点对与对应的光源构成一个三角形，最后通过正弦定理计算每个光源位置。Yanhong Liu等人提出一种实时追踪多光源位置的增强现实系统。他们使用SURF(Speeded Up Robust Features)算法进行相机位姿追踪；在自然特征平面标记中心上放置一个纯镜面反射球，基于逆向光线追踪原理，单幅图像恢复多光源方向，两幅图像定位多光源位置。

发明内容

为了准确地对复杂环境中多光源位置进行定位，本发明提出一种多光源位置估计方法，该方法可以有效地测算和恢复复杂场景中多个光源的位置，解决增强现实中基于多光源的虚实融合场景的无缝融合等问题，从而，使得虚拟对象叠加到真实场景中更具沉浸感，更加逼真。

本发明的技术方案是这样实现的：一种多光源位置估计方法，包括计算机1，摄像头2，第一个镜面反射球3，第二个镜面反射球4，实验台5，光源组6；其中计算机1和摄像头2通过电缆连接，第一个镜面反射球3和第二个镜面反射球4放置在实验台5上，第一个镜面反射球3的半径用Sr1表示，第一个镜面反射球4的半径用Sr2表示，光源组6由G个光源组成，分别称为第i个光源，i＝1,2,...,G；其特征在于具体步骤如下：

步骤1、记世界坐标系与摄像头2坐标系重合；利用计算机1控制摄像头2对准实验台5拍摄一幅包含第一个镜面反射球3和第二个镜面反射球4及第一个镜面反射球3和第二个镜面反射球4上分别有G个高亮点区域的彩色图像,该彩色图像用srcimage表示。

步骤2、计算机1利用Matlab2015中的rgb2gray函数将srcimage转化成灰度图像，该灰度图像用grayimage表示；设置一个阈值，用threshold1表示，利用Matlab2015中的im2bw函数和threshold1，将grayimage转化成二值图像，该二值图像用binaryimage1表示；利用Matlab2015中的bwperim函数提取binaryimage1中第一个镜面反射球3和第二个镜面反射球4的轮廓，并分别用contours1和contours2表示；

设contours1中有M的像素点，由这些像素点在图像坐标系中的坐标构成矩阵Q1，Q1(i,j)是矩阵Q1的第i行，第j列元素，i＝1,2,…,M,j＝1,2,3；Q1(i,1)为第i个像素点在图像坐标系中的X轴坐标，Q1(i,2)为第i个像素点在图像坐标系中的Y轴坐标，Q1(i,3)的值为1，i＝1,2,…,M。

设contours2中有N的像素点，由这些像素点在图像坐标系中的坐标构成矩阵Q2，Q2(i,j)是矩阵Q2的第i行，第j列元素，i＝1,2,…,N,j＝1,2,3；Q2(i,1)为第i个像素点在图像坐标系中的X轴坐标，Q2(i,2)为第i个像素点在图像坐标系中的Y轴坐标，Q2(i,3)的值为1，i＝1,2,…,N。

步骤3、利用公式

ZQ1＝Q1*K^-1，ZQ2＝Q2*K^-1

构成矩阵ZQ1、ZQ2，ZQ1(i,j)是矩阵ZQ1的第i行，第j列元素，i＝1,2,…,M,j＝1,2,3；ZQ2(i,j)是矩阵ZQ2的第i行，第j列元素，i＝1,2,…,N,j＝1,2,3；K为摄像头2的内参矩阵。

步骤4、构造方程组

A1*ZQ1(p,1)²+B1*ZQ1(p,1)*ZQ1(p,2)+C1*ZQ1(p,2)²+D1*ZQ1(p,1)+E1*ZQ1(p,2)+F1＝0，其中，p＝1,2,...,M，求解A1,B1,C1,D1,E1,F1，利用下式

构成矩阵

构造方程组

A2*ZQ2(q,1)²+B2*ZQ2(q,1)*ZQ2(q,2)+C2*ZQ2(q,2)²+D2*ZQ2(q,1)+E2*ZQ1(q,2)+F2＝0，其中，q＝1,2,...,N，求解A2,B2,C2,D2,E2,F2，利用下式

构成矩阵

步骤5、根据实对称矩阵的正交对角化分解原理，分别使用公式

对和进行分解，其中，U1和U2都是正交矩阵，C1和C2都是对角矩阵，C1(i,j)是矩阵C1的第i行，第j列元素，i＝1,2,3,j＝1,2,3，C2(i,j)是矩阵C2的第i行，第j列元素，i＝1,2,3,j＝1,2,3；设两个正常数，分别用r1和r2表示，r1和r2可根据下式计算得到。

步骤6、设第一个镜面反射球3的球心在世界坐标系下的坐标用向量SC11表示，根据公式

计算可得SC11。其中，u13是矩阵U1的第三列向量。

步骤7、设第二个镜面反射球4的球心在世界坐标系下的坐标用向量SC22表示，根据公式

计算可得SC22。其中，u23是矩阵U2的第三列向量。

步骤8、设置一个阈值，用threshold2表示，利用Matlab2015中的im2bw函数和threshold2，将grayimage转化成二值图像，该二值图像用binaryimage2表示。

步骤9、通过手动方式，根据属于同一光源产生的高亮点区域的原则，对binaryimage2中的高亮点区域进行匹配；在binaryimage2中第一个镜面反射球3上由第i个光源产生得到的高亮点区域，为第一个镜面反射球3第i高亮点区域，在binaryimage2中第二个镜面反射球4上由第i个光源产生得到的高亮点区域，为第二个镜面反射球4第i高亮点区域，i＝1,2,...G；使用质心算法，计算第一个镜面反射球3第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL1(i,1)和GL1(i,2),GL1(i,3)＝1,i＝1,2,...G，计算第二个镜面反射球4第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL2(i,1)和GL2(i,2),GL2(i,3)＝1,i＝1,2,...G。

步骤10、利用公式

H1＝GL1*K^-1，H2＝GL2*K^-1

构成矩阵H1、H2，其中，K为摄像头2的内参矩阵。

步骤11、利用公式

构造矩阵V1，V1(i,j)是矩阵V1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；利用公式

构造矩阵V2，V2(i,j)是矩阵V2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3。

步骤12、构造矩阵SS1和SS2；SS1(i,j)是矩阵SS1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3,4,5,6；SS2(i,j)是矩阵SS2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3,4,5,6；构造方程组

求解矩阵SS1；构造矩阵S1，S1(i,j)是矩阵S1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；如果SS1(i,3)<SS1(i,6)，那么S1(i,j)＝SS1(i,j)，i＝1,2,…,G,j＝1,2,3；否则，S1(i,1)＝SS1(i,4)，S1(i,2)＝SS1(i,5)，S1(i,3)＝SS1(i,6),i＝1,2,…,G。

构造方程组

求解矩阵SS2。构造矩阵S2，S2(i,j)是矩阵S2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；如果SS2(i,3)<SS2(i,6)，那么S2(i,j)＝SS2(i,j)，i＝1,2,…,G,j＝1,2,3；否则，S2(i,1)＝SS2(i,4)，S2(i,2)＝SS2(i,5)，S2(i,3)＝SS2(i,6),i＝1,2,…,G。

步骤13、利用公式

构造矩阵N1，N1(i,j)是矩阵N1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；利用公式

构造矩阵N2，N2(i,j)是矩阵N2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3。

步骤14、利用公式

L1(i,j)＝-2*(V1(i,1)*N1(i,1)+V1(i,2)*N1(i,2)+V1(i,3)*N1(i,3))*N1(i,j)+V1(i,j)

构造矩阵L1，L1(i,j)是矩阵L1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；利用公式

L2(i,j)＝-2*(V2(i,1)*N2(i,1)+V2(i,2)*N2(i,2)+V2(i,3)*N2(i,3))*N2(i,j)+V2(i,j)

构造矩阵L2，L2(i,j)是矩阵L2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3。

步骤15、构造矩阵a1，a2，a1(i,j)是矩阵a1的第i行，第j列元素，i＝1,2,…,G,j＝1；a2(i,j)是矩阵a2的第i行，第j列元素，i＝1,2,…,G,j＝1；构造矩阵P1，P2，P1(i,j)是矩阵P1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；P2(i,j)是矩阵P2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；构造方程组

求解矩阵a1，a2和矩阵P1，P2。

步骤16、构造矩阵P，P(i,j)是矩阵P的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；通过下式

求解矩阵P。(P(i,1),P(i,2),P(i,3))i＝1,2,…,G,表示第i个光源在世界坐标系中的坐标。

通过以上步骤实现了多光源位置的估计。

本发明的积极效果是适用于不同距离的多个点光源，能够有效地估计出光源的位置；能够精确定位镜面反射球在世界坐标系的球心坐标，在计算过程中，通过设定世界坐标系与相机坐标系重合消除了坐标系之间转换带来的误差；对于实验中来自同一照明的两条反射光线不能相交于一个3D空间点问题，采用与两条反射光线都垂直的最短线段并取中点作为估计的光源位置；最终获得较为精确的多光源位置估计值，，有效地解决多光源下虚实物体无缝融合等问题，从而，增强了融合场景的沉浸感和真实感。

附图说明

图1是一种多光源位置估计方法所需设备构成图。其中：1为计算机，2为摄像头，3为第一个镜面反射球，4为第二个镜面反射球，5为实验台，6为光源组。

图2实验彩色图像。

具体实施方式

下面结合附图对本发明做进一步的描述：如图1所示，一种多光源位置方法，包括计算机1，摄像头2，第一个镜面反射球3，第二个镜面反射球4，实验台5，光源组6；其中计算机1和摄像头2通过电缆连接，第一个镜面反射球3和第二个镜面反射球4放置在实验台5上，第一个镜面反射球3的半径用Sr1表示，第一个镜面反射球4的半径用Sr2表示，光源组6由G个光源组成，分别称为第i个光源，i＝1,2,...,G。其特征在于具体步骤如下：

步骤1、记世界坐标系与摄像头2坐标系重合。利用计算机1控制摄像头2对准实验台5拍摄一幅包含第一个镜面反射球3和第二个镜面反射球4及第一个镜面反射球3和第二个镜面反射球4上分别有G＝2个高亮点区域的彩色图像,该彩色图像用srcimage表示，如图2所示，彩色图像srcimage的尺寸为640*480。

步骤2、计算机1利用Matlab2015中的rgb2gray函数将srcimage转化成灰度图像，该灰度图像用grayimage表示；设置一个阈值threshold1＝100，利用Matlab2015中的im2bw函数和threshold1，将grayimage转化成二值图像，该二值图像用binaryimage1表示；利用Matlab2015中的bwperim函数提取binaryimage1中第一个镜面反射球3和第二个镜面反射球4的轮廓，并分别用contours1和contours2表示；

步骤3、利用公式

ZQ1＝Q1*K^-1，ZQ2＝Q2*K^-1

构成矩阵ZQ1、ZQ2，ZQ1(i,j)是矩阵ZQ1的第i行，第j列元素，i＝1,2,…,M,j＝1,2,3；ZQ2(i,j)是矩阵ZQ2的第i行，第j列元素，i＝1,2,…,N,j＝1,2,3；K为摄像头2的内参矩阵，即

步骤4、构造方程组

A1*ZQ1(p,1)²+B1*ZQ1(p,1)*ZQ1(p,2)+C1*ZQ1(p,2)²+D1*ZQ1(p,1)+E1*ZQ1(p,2)+F1＝0，其中，p＝1,2,...,M；求解A1,B1,C1,D1,E1,F1，利用下式

构成矩阵即

构造方程组

A2*ZQ2(q,1)²+B2*ZQ2(q,1)*ZQ2(q,2)+C2*ZQ2(q,2)²+D2*ZQ2(q,1)+E2*ZQ1(q,2)+F2＝0，其中，q＝1,2,...,N；求解A2,B2,C2,D2,E2,F2，利用下式

构成矩阵即

r1＝0.058053432442408，r2＝0.071989874769235

计算可得SC11＝(20.428455453126958,69.064549202249740,410.6574005211783)，其中，u13是矩阵U1的第三列向量。

计算可得SC22＝(111.5112625894254,-25.958894125159865,331.4576199603303)，其中，u23是矩阵U2的第三列向量。

步骤8、设置一个阈值threshold2＝250，利用Matlab2015中的im2bw函数和threshold2，将grayimage转化成二值图像，该二值图像用binaryimage2表示。

步骤9、通过手动方式，根据属于同一光源产生的高亮点区域的原则，对binaryimage2中的高亮点区域进行匹配；在binaryimage2中第一个镜面反射球3上由第i个光源产生得到的高亮点区域，为第一个镜面反射球3第i高亮点区域，在binaryimage2中第二个镜面反射球4上由第i个光源产生得到的高亮点区域，为第二个镜面反射球4第i高亮点区域，i＝1,2；使用质心算法，计算第一个镜面反射球3第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL1(i,1)和GL1(i,2),GL1(i,3)＝1,i＝1,2，计算第二个镜面反射球4第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL2(i,1)和GL2(i,2),GL2(i,3)＝1,i＝1,2。

步骤10、利用公式

H1＝GL1*K^-1，H2＝GL2*K^-1

构成矩阵H1、H2，其中，K为摄像头2的内参矩阵。经计算得

步骤11、利用公式

构造矩阵V1，V1(i,j)是矩阵V1的第i行，第j列元素，i＝1,2,j＝1,2,3；利用公式

构造矩阵V2，V2(i,j)是矩阵V2的第i行，第j列元素，i＝1,2,j＝1,2,3；经计算得

步骤12、构造矩阵SS1和SS2；SS1(i,j)是矩阵SS1的第i行，第j列元素，i＝1,2,j＝1,2,3,4,5,6；SS2(i,j)是矩阵SS2的第i行，第j列元素，i＝1,2,j＝1,2,3,4,5,6；构造方程组

求解矩阵SS1。构造矩阵S1，S1(i,j)是矩阵S1的第i行，第j列元素，i＝1,2,j＝1,2,3；如果SS1(i,3)<SS1(i,6),那么S1(i,j)＝SS1(i,j),i＝1,2,j＝1,2,3；否则，S1(i,1)＝SS1(i,4),S1(i,2)＝SS1(i,5),S1(i,3)＝SS1(i,6),i＝1,2；经计算得

构造方程组

求解矩阵SS2；构造矩阵S2，S2(i,j)是矩阵S2的第i行，第j列元素，i＝1,2,j＝1,2,3；如果SS2(i,3)<SS2(i,6),那么S2(i,j)＝SS2(i,j),i＝1,2,j＝1,2,3；否则，S2(i,1)＝SS2(i,4),S2(i,2)＝SS2(i,5),S2(i,3)＝SS2(i,6),i＝1,2；经计算得

步骤13、利用公式

构造矩阵N2，N2(i,j)是矩阵N2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；经计算得

步骤14、利用公式

L1(i,j)＝-2*(V1(i,1)*N1(i,1)+V1(i,2)*N1(i,2)+V1(i,3)*N1(i,3))*N1(i,j)+V1(i,j)

构造矩阵L1，L1(i,j)是矩阵L1的第i行，第j列元素，i＝1,2,j＝1,2,3；利用公式

L2(i,j)＝-2*(V2(i,1)*N2(i,1)+V2(i,2)*N2(i,2)+V2(i,3)*N2(i,3))*N2(i,j)+V2(i,j)

构造矩阵L2，L2(i,j)是矩阵L2的第i行，第j列元素，i＝1,2,j＝1,2,3；经计算得

步骤15、构造矩阵a1，a2，a1(i,j)是矩阵a1的第i行，第j列元素，i＝1,2,j＝1；a2(i,j)是矩阵a2的第i行，第j列元素，i＝1,2,j＝1；构造矩阵P1，P2；P1(i,j)是矩阵P1的第i行，第j列元素，i＝1,2,j＝1,2,3；P2(i,j)是矩阵P2的第i行，第j列元素，i＝1,2,j＝1,2,3；构造方程组

求解矩阵a1，a2和矩阵P1，P2；经计算得

求解矩阵P。(P(i,1),P(i,2),P(i,3))i＝1,2,表示第i个光源在世界坐标系中的坐标；经计算得

通过以上步骤实现了多光源位置的估计。

Claims

1.一种多光源位置估计方法，包括计算机1，摄像头2，第一个镜面反射球3，第二个镜面反射球4，实验台5，光源组6；其中计算机1和摄像头2通过电缆连接，第一个镜面反射球3和第二个镜面反射球4放置在实验台5上，第一个镜面反射球3的半径用Sr1表示，第一个镜面反射球4的半径用Sr2表示，光源组6由G个光源组成，分别称为第i个光源，i＝1,2,...,G；其特征在于具体步骤如下：

步骤1、记世界坐标系与摄像头2坐标系重合；利用计算机1控制摄像头2对准实验台5拍摄一幅包含第一个镜面反射球3和第二个镜面反射球4及第一个镜面反射球3和第二个镜面反射球4上分别有G个高亮点区域的彩色图像,该彩色图像用srcimage表示；

设contours1中有M的像素点，由这些像素点在图像坐标系中的坐标构成矩阵Q1，Q1(i,j)是矩阵Q1的第i行，第j列元素，i＝1,2,…,M,j＝1,2,3；Q1(i,1)为第i个像素点在图像坐标系中的X轴坐标，Q1(i,2)为第i个像素点在图像坐标系中的Y轴坐标，Q1(i,3)的值为1，i＝1,2,…,M；

设contours2中有N的像素点，由这些像素点在图像坐标系中的坐标构成矩阵Q2，Q2(i,j)是矩阵Q2的第i行，第j列元素，i＝1,2,…,N,j＝1,2,3；Q2(i,1)为第i个像素点在图像坐标系中的X轴坐标，Q2(i,2)为第i个像素点在图像坐标系中的Y轴坐标，Q2(i,3)的值为1，i＝1,2,…,N；

步骤3、利用公式

ZQ1＝Q1*K^-1，ZQ2＝Q2*K^-1

构成矩阵ZQ1、ZQ2，ZQ1(i,j)是矩阵ZQ1的第i行，第j列元素，i＝1,2,…,M,j＝1,2,3；ZQ2(i,j)是矩阵ZQ2的第i行，第j列元素，i＝1,2,…,N,j＝1,2,3；K为摄像头2的内参矩阵；

步骤4、构造方程组

<mrow> <mover> <mi>Q</mi> <mo>&OverBar;</mo> </mover> <mn>1</mn> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>A</mi> <mn>1</mn> </mrow> </mtd> <mtd> <mrow> <mi>B</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>D</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>B</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>C</mi> <mn>1</mn> </mrow> </mtd> <mtd> <mrow> <mi>E</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>D</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>E</mi> <mn>1</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>F</mi> <mn>1</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

构成矩阵

构造方程组

<mrow> <mover> <mi>Q</mi> <mo>&OverBar;</mo> </mover> <mn>2</mn> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>A</mi> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>B</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>D</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>B</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>C</mi> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>E</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>D</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>E</mi> <mn>2</mn> <mo>/</mo> <mn>2</mn> </mrow> </mtd> <mtd> <mrow> <mi>F</mi> <mn>2</mn> </mrow> </mtd> </mtr> </mtable> </mfenced> </mrow>

构成矩阵

<mrow> <mover> <mi>Q</mi> <mo>&OverBar;</mo> </mover> <mn>1</mn> <mo>=</mo> <mi>U</mi> <mn>1</mn> <mi>C</mi> <mn>1</mn> <mi>U</mi> <msup> <mn>1</mn> <mi>T</mi> </msup> <mo>,</mo> <mover> <mi>Q</mi> <mo>&OverBar;</mo> </mover> <mn>2</mn> <mo>=</mo> <mi>U</mi> <mn>2</mn> <mi>C</mi> <mn>1</mn> <mi>U</mi> <msup> <mn>2</mn> <mi>T</mi> </msup> </mrow>

计算可得SC11。其中，u13是矩阵U1的第三列向量；

计算可得SC22；其中，u23是矩阵U2的第三列向量；

步骤8、设置一个阈值，用threshold2表示，利用Mat lab2015中的im2bw函数和threshold2，将grayimage转化成二值图像，该二值图像用binaryimage2表示；

步骤9、通过手动方式，根据属于同一光源产生的高亮点区域的原则，对binaryimage2中的高亮点区域进行匹配；在binaryimage2中第一个镜面反射球3上由第i个光源产生得到的高亮点区域，为第一个镜面反射球3第i高亮点区域，在binaryimage2中第二个镜面反射球4上由第i个光源产生得到的高亮点区域，为第二个镜面反射球4第i高亮点区域，i＝1,2,...G；使用质心算法，计算第一个镜面反射球3第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL1(i,1)和GL1(i,2),GL1(i,3)＝1,i＝1,2,...G，计算第二个镜面反射球4第i高亮点区域的质心在图像坐标系中的X轴和Y轴的坐标分别为GL2(i,1)和GL2(i,2),GL2(i,3)＝1,i＝1,2,...G；

步骤10、利用公式

H1＝GL1*K^-1，H2＝GL2*K^-1

构成矩阵H1、H2，其中，K为摄像头2的内参矩阵；

步骤11、利用公式

构造矩阵V2，V2(i,j)是矩阵V2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；

求解矩阵SS1；构造矩阵S1，S1(i,j)是矩阵S1的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；如果SS1(i,3)<SS1(i,6)，那么S1(i,j)＝SS1(i,j)，i＝1,2,…,G,j＝1,2,3；否则，S1(i,1)＝SS1(i,4)，S1(i,2)＝SS1(i,5)，S1(i,3)＝SS1(i,6),i＝1,2,…,G；

构造方程组

求解矩阵SS2。构造矩阵S2，S2(i,j)是矩阵S2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；如果SS2(i,3)<SS2(i,6)，那么S2(i,j)＝SS2(i,j)，i＝1,2,…,G,j＝1,2,3；否则，S2(i,1)＝SS2(i,4)，S2(i,2)＝SS2(i,5)，S2(i,3)＝SS2(i,6),i＝1,2,…,G；

步骤13、利用公式

构造矩阵N2，N2(i,j)是矩阵N2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；

步骤14、利用公式

L1(i,j)＝-2*(V1(i,1)*N1(i,1)+V1(i,2)*N1(i,2)+V1(i,3)*N1(i,3))*N1(i,j)+V1(i,j)

L2(i,j)＝-2*(V2(i,1)*N2(i,1)+V2(i,2)*N2(i,2)+V2(i,3)*N2(i,3))*N2(i,j)+V2(i,j)

构造矩阵L2，L2(i,j)是矩阵L2的第i行，第j列元素，i＝1,2,…,G,j＝1,2,3；

求解矩阵a1，a2和矩阵P1，P2；

求解矩阵P；(P(i,1),P(i,2),P(i,3))i＝1,2,…,G,表示第i个光源在世界坐标系中的坐标；

通过以上步骤实现了多光源位置的估计。