发明名称 基于双目显微立体视觉的零件精确定位方法
摘要 本发明基于双目显微立体视觉的零件精确定位方法属于计算机视觉测量技术领域,涉及一种基于双目显微立体视觉的精密零件的精确定位方法。采用双目显微立体视觉系统,利用两个CCD摄像机采集被测零件的图像,立体显微镜将被测零件上的待测区域的图像信息放大;采用棋盘格标定板对两个CCD摄像机进行标定;采用Harris角点检测算法及亚像素提取算法进行特征点的提取。将提取后的特征点进行初匹配和匹配点对的校正,将特征点图像坐标输入到标定好的系统中得到特征点的空间实际坐标。本发明解决了由于待测目标区域小、定位精度要求高、非接触等问题所产生的测量难题。采用基于双目显微立体视觉的非接触式测量方法,很好的完成了精密零件的精确定位。
申请公布号 CN103247053B 申请公布日期 2015.10.14
申请号 CN201310182221.0 申请日期 2013.05.16
申请人 大连理工大学 发明人 刘巍;贾振元;屠先明;王福吉;王文强;赵凯
分类号 G06T7/00(2006.01)I 主分类号 G06T7/00(2006.01)I
代理机构 大连理工大学专利中心 21200 代理人 关慧贞
主权项 一种基于双目显微立体视觉的零件精确定位方法,其特征是,采用双目显微立体视觉系统,利用左右两个CCD摄像机(2、2’)采集被测零件(4)的图像位置信息,立体显微镜(3)将被测零件(4)上的待测区域(6)的图像位置信息进行放大;放大后的图像经过图像采集卡传到工业控制计算机(9)内,采用精度较高的棋盘格标定板对左右两个CCD摄像机进行标定;采用Harris角点检测算法及亚像素提取算法进行待测特征点的提取;将提取后的待测特征点进行初匹配和匹配点对的校正;将已经匹配好的特征点图像坐标输入到标定好的系统中得到特征点的空间实际坐标;测量方法的具体步骤如下:(1)左右两个CCD摄像机的标定左右两个CCD摄像机的标定包括摄像机内参数和摄像机外参数;摄像机内参数包括尺度因子α、β,主点坐标u<sub>0</sub>,v<sub>0</sub>,以及垂直因子γ;在标定过程中,需要先得到摄像机的五个内参数,在求得内参数矩阵的基础上,求解摄像机外参数矩阵;尺度因子是空间点在经过平移旋转变换之后,它在摄像机坐标系中的坐标与它在图像坐标系中的坐标的尺度关系,本发明对尺度因子的标定采用张氏摄像机标定法;根据摄像机的针孔模型可建立图像坐标系与世界坐标系的关系公式为:<maths num="0001" id="cmaths0001"><math><![CDATA[<mrow><msub><mi>Z</mi><mi>c</mi></msub><mfenced open='[' close=']'><mtable><mtr><mtd><mi>u</mi></mtd></mtr><mtr><mtd><mi>v</mi></mtd></mtr><mtr><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><mi>&alpha;</mi></mtd><mtd><mi>&gamma;</mi></mtd><mtd><msub><mi>u</mi><mn>0</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mi>&beta;</mi></mtd><mtd><msub><mi>v</mi><mn>0</mn></msub></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mn>0</mn></mtd><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mfenced open='[' close=']'><mtable><mtr><mtd><mi>R</mi></mtd><mtd><mi>t</mi></mtd></mtr></mtable></mfenced><mfenced open='[' close=']'><mtable><mtr><mtd><msub><mi>X</mi><mi>w</mi></msub></mtd></mtr><mtr><mtd><msub><mi>Y</mi><mi>w</mi></msub></mtd></mtr><mtr><mtd><msub><mi>Z</mi><mi>w</mi></msub></mtd></mtr><mtr><mtd><mn>1</mn></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000011.GIF" wi="1383" he="334" /></maths>其中:α和β就是需要标定的尺度因子,X<sub>w</sub>,Y<sub>w</sub>,Z<sub>w</sub>是空间中一点P的三维空间坐标,u,v是P点在图像上的图像坐标,R,t代表了摄像机坐标系相对于世界坐标系的旋转和平移矩阵;若令X′<sub>w</sub>=(X<sub>w</sub>,Y<sub>w</sub>,Z<sub>w</sub>)<sup>T</sup>,<img file="FDA0000741774530000012.GIF" wi="305" he="105" />Z<sub>c</sub>=s,则<maths num="0002" id="cmaths0002"><math><![CDATA[<mrow><mi>s</mi><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mi>w</mi></msub><mo>=</mo><mi>H</mi><msubsup><mi>X</mi><mi>w</mi><mo>&prime;</mo></msubsup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000013.GIF" wi="1384" he="82" /></maths>其中,H也称为单应性矩阵,且有H=K[r<sub>1</sub>,r<sub>2</sub>,t]                  (3)每一个位置的标定板图像都独立对应一个单应性矩阵,令H=[h<sub>1</sub>,h<sub>2</sub>,h<sub>3</sub>]则可由上式推得:[h<sub>1</sub> h<sub>2</sub> h<sub>3</sub>]=λK[r<sub>1</sub> r<sub>2</sub> t]        (4)其中,λ为任意的比例因子,K摄像机内参数矩阵,因r<sub>1</sub>和r<sub>2</sub>是单位正交向量,即<maths num="0003" id="cmaths0003"><math><![CDATA[<mrow><msup><msub><mi>r</mi><mn>1</mn></msub><mi>T</mi></msup><msub><mi>r</mi><mn>1</mn></msub><mo>=</mo><msup><msub><mi>r</mi><mn>2</mn></msub><mi>T</mi></msup><msub><mi>r</mi><mn>2</mn></msub><mo>=</mo><mn>1</mn></mrow>]]></math><img file="FDA0000741774530000021.GIF" wi="442" he="121" /></maths>且<maths num="0004" id="cmaths0004"><math><![CDATA[<mrow><msup><msub><mi>r</mi><mn>1</mn></msub><mi>T</mi></msup><msub><mi>r</mi><mn>2</mn></msub><mo>=</mo><mn>0</mn></mrow>]]></math><img file="FDA0000741774530000022.GIF" wi="256" he="121" /></maths>则:<maths num="0005" id="cmaths0005"><math><![CDATA[<mrow><msubsup><mi>h</mi><mn>1</mn><mi>T</mi></msubsup><msup><mi>K</mi><mrow><mo>-</mo><mi>T</mi></mrow></msup><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msub><mi>h</mi><mn>1</mn></msub><mo>=</mo><msubsup><mi>h</mi><mn>2</mn><mi>T</mi></msubsup><msup><mi>K</mi><mrow><mo>-</mo><mi>T</mi></mrow></msup><msup><mi>K</mi><mrow><mo>-</mo><mn>1</mn></mrow></msup><msub><mi>h</mi><mn>2</mn></msub><mo>=</mo><mn>1</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>5</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000023.GIF" wi="1384" he="103" /></maths>h<sub>1</sub>K<sup>‑T</sup>K<sup>‑1</sup>h<sub>2</sub>=0               (6)根据公式(5)公式(6),计算出立体显微镜左右光路的摄像机尺度因子;针对主点坐标,采用变倍率法对摄像机的主点进行标定;假设空间中任意一点P在摄像机坐标系中的坐标为x<sub>1</sub>,y<sub>1</sub>,z<sub>1</sub>,在不考虑非线性畸变以及图像坐标垂直度的情况下,该点的图像平面投影的坐标为:u<sub>1</sub>=rx<sub>1</sub>+u<sub>0</sub>                       (7)v<sub>1</sub>=ry<sub>1</sub>+v<sub>0</sub>其中u<sub>1</sub>,v<sub>1</sub>是该点的图像坐标,r是任一放大倍率,u<sub>0</sub>,v<sub>0</sub>是主点坐标;将式中的r消去,可得直线方程:<maths num="0006" id="cmaths0006"><math><![CDATA[<mrow><mfrac><mrow><msub><mi>u</mi><mn>1</mn></msub><mo>-</mo><msub><mi>u</mi><mn>0</mn></msub></mrow><msub><mi>x</mi><mn>1</mn></msub></mfrac><mo>=</mo><mfrac><mrow><msub><mi>v</mi><mn>1</mn></msub><mo>-</mo><msub><mi>v</mi><mn>0</mn></msub></mrow><msub><mi>y</mi><mn>1</mn></msub></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000024.GIF" wi="1384" he="157" /></maths>即在任何放大倍率下,点x<sub>1</sub>,y<sub>1</sub>,z<sub>1</sub>的图像坐标都在同一直线上,且该直线一定经过摄像机主点u<sub>0</sub>,v<sub>0</sub>;为了求取摄像机主点u<sub>0</sub>,v<sub>0</sub>,选定12个在图像投影平面上不重合的点,令这12个点分别在不同放大倍率下0.7、1、2、2.4、3、5、5.8、8进行投影,利用最小二乘法对获得的12条直线进行拟合,并用最小二乘法求取12条直线的交点坐标即摄像机主点u<sub>0</sub>,v<sub>0</sub>;由于垂直因子γ对成像精度影响不大,因此在初始参数标定时,将γ设为0,通过优化手段得到γ具体值;针对摄像机外参数旋转矩阵R和平移矩阵t;由公式(2)知,在标定完内参数矩阵K和单应性矩阵H后可以求得摄像机旋转矩阵R=[r<sub>1</sub>,r<sub>2</sub>,r<sub>3</sub>]和平移向量t如下:r<sub>1</sub>=λK<sup>‑1</sup>h<sub>1</sub>r<sub>2</sub>=λK<sup>‑1</sup>h<sub>2</sub>                           (9)r<sub>3</sub>=r<sub>1</sub>×r<sub>2</sub>t=λK<sup>‑1</sup>h<sub>3</sub>在初步标定完左右两个CCD摄像机的主参数后,下一步是对主参数和畸变参数的优化;为了简化优化参数,本发明采用四元数法将旋转矩阵<maths num="0007" id="cmaths0007"><math><![CDATA[<mrow><mi>R</mi><mo>=</mo><mo>[</mo><msub><mi>r</mi><mn>1</mn></msub><mo>,</mo><msub><mi>r</mi><mn>2</mn></msub><mo>,</mo><msub><mi>r</mi><mn>3</mn></msub><mo>]</mo><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><msub><mi>r</mi><mn>11</mn></msub></mtd><mtd><msub><mi>r</mi><mn>21</mn></msub></mtd><mtd><msub><mi>r</mi><mn>31</mn></msub></mtd></mtr><mtr><mtd><msub><mi>r</mi><mn>12</mn></msub></mtd><mtd><msub><mi>r</mi><mn>22</mn></msub></mtd><mtd><msub><mi>r</mi><mn>32</mn></msub></mtd></mtr><mtr><mtd><msub><mi>r</mi><mn>13</mn></msub></mtd><mtd><msub><mi>r</mi><mn>23</mn></msub></mtd><mtd><msub><mi>r</mi><mn>33</mn></msub></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA0000741774530000031.GIF" wi="710" he="255" /></maths>中的九个未知数缩减为四个,公式如下:<maths num="0008" id="cmaths0008"><math><![CDATA[<mrow><msub><mi>q</mi><mn>1</mn></msub><mo>=</mo><mfrac><msqrt><mn>1</mn><mo>+</mo><msub><mi>r</mi><mn>11</mn></msub><mo>+</mo><msub><mi>r</mi><mn>22</mn></msub><mo>+</mo><msub><mi>r</mi><mn>33</mn></msub></msqrt><mn>2</mn></mfrac><mo>,</mo><msub><mi>q</mi><mn>2</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>23</mn></msub><mo>-</mo><msub><mi>r</mi><mn>32</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>1</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>3</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>31</mn></msub><mo>-</mo><msub><mi>r</mi><mn>13</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>1</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>4</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>12</mn></msub><mo>-</mo><msub><mi>r</mi><mn>21</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>1</mn></msub></mrow></mfrac></mrow>]]></math><img file="FDA0000741774530000032.GIF" wi="1235" he="160" /></maths><maths num="0009" id="cmaths0009"><math><![CDATA[<mrow><msub><mi>q</mi><mn>2</mn></msub><mo>=</mo><mfrac><msqrt><mn>1</mn><mo>+</mo><msub><mi>r</mi><mn>11</mn></msub><mo>-</mo><msub><mi>r</mi><mn>22</mn></msub><mo>-</mo><msub><mi>r</mi><mn>33</mn></msub></msqrt><mn>2</mn></mfrac><mo>,</mo><msub><mi>q</mi><mn>1</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>23</mn></msub><mo>-</mo><msub><mi>r</mi><mn>32</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>2</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>3</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>12</mn></msub><mo>+</mo><msub><mi>r</mi><mn>21</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>2</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>4</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>31</mn></msub><mo>-</mo><msub><mi>r</mi><mn>13</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>2</mn></msub></mrow></mfrac></mrow>]]></math><img file="FDA0000741774530000033.GIF" wi="1232" he="161" /></maths>                         (10)<maths num="0010" id="cmaths0010"><math><![CDATA[<mrow><msub><mi>q</mi><mn>3</mn></msub><mo>=</mo><mfrac><msqrt><mn>1</mn><mo>-</mo><msub><mi>r</mi><mn>11</mn></msub><mo>+</mo><msub><mi>r</mi><mn>22</mn></msub><mo>-</mo><msub><mi>r</mi><mn>33</mn></msub></msqrt><mn>2</mn></mfrac><mo>,</mo><msub><mi>q</mi><mn>1</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>31</mn></msub><mo>-</mo><msub><mi>r</mi><mn>13</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>3</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>2</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>12</mn></msub><mo>-</mo><msub><mi>r</mi><mn>21</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>3</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>4</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>23</mn></msub><mo>-</mo><msub><mi>r</mi><mn>32</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>3</mn></msub></mrow></mfrac></mrow>]]></math><img file="FDA0000741774530000034.GIF" wi="1233" he="163" /></maths><maths num="0011" id="cmaths0011"><math><![CDATA[<mrow><msub><mi>q</mi><mn>4</mn></msub><mo>=</mo><mfrac><msqrt><mn>1</mn><mo>-</mo><msub><mi>r</mi><mn>11</mn></msub><mo>-</mo><msub><mi>r</mi><mn>22</mn></msub><mo>+</mo><msub><mi>r</mi><mn>33</mn></msub></msqrt><mn>2</mn></mfrac><mo>,</mo><msub><mi>q</mi><mn>1</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>12</mn></msub><mo>-</mo><msub><mi>r</mi><mn>21</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>4</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>2</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>31</mn></msub><mo>+</mo><msub><mi>r</mi><mn>13</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>4</mn></msub></mrow></mfrac><mo>,</mo><msub><mi>q</mi><mn>3</mn></msub><mo>=</mo><mfrac><mrow><msub><mi>r</mi><mn>23</mn></msub><mo>-</mo><msub><mi>r</mi><mn>32</mn></msub></mrow><mrow><mn>4</mn><msub><mi>q</mi><mn>4</mn></msub></mrow></mfrac></mrow>]]></math><img file="FDA0000741774530000035.GIF" wi="1237" he="163" /></maths>从公式(10)计算得到新的旋转矩阵:<maths num="0012" id="cmaths0012"><math><![CDATA[<mrow><mi>R</mi><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><mo>-</mo><mn>1</mn><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>2</mn><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>1</mn><mn>2</mn></msubsup></mtd><mtd><mn>2</mn><msub><mi>q</mi><mn>2</mn></msub><msub><mi>q</mi><mn>3</mn></msub><mo>-</mo><mn>2</mn><msub><mi>q</mi><mn>1</mn></msub><msub><mi>q</mi><mn>4</mn></msub></mtd><mtd><mn>2</mn><msub><mi>q</mi><mn>1</mn></msub><msub><mi>q</mi><mn>3</mn></msub><mo>+</mo><mn>2</mn><msub><mi>q</mi><mn>2</mn></msub><msub><mi>q</mi><mn>4</mn></msub></mtd></mtr><mtr><mtd><mn>2</mn><msub><mi>q</mi><mn>2</mn></msub><msub><mi>q</mi><mn>3</mn></msub><mo>+</mo><mn>2</mn><msub><mi>q</mi><mn>4</mn></msub><msub><mi>q</mi><mn>1</mn></msub></mtd><mtd><mo>-</mo><mn>1</mn><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>1</mn><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>3</mn><mn>2</mn></msubsup></mtd><mtd><mo>-</mo><mn>2</mn><msub><mi>q</mi><mn>1</mn></msub><msub><mi>q</mi><mn>2</mn></msub><mo>+</mo><mn>2</mn><msub><mi>q</mi><mn>3</mn></msub><msub><mi>q</mi><mn>4</mn></msub></mtd></mtr><mtr><mtd><mo>-</mo><mn>2</mn><msub><mi>q</mi><mn>1</mn></msub><msub><mi>q</mi><mn>3</mn></msub><mo>+</mo><mn>2</mn><msub><mi>q</mi><mn>4</mn></msub><msub><mi>q</mi><mn>2</mn></msub></mtd><mtd><mn>2</mn><msub><mi>q</mi><mn>1</mn></msub><msub><mi>q</mi><mn>2</mn></msub><mo>+</mo><mn>2</mn><msub><mi>q</mi><mn>3</mn></msub><msub><mi>q</mi><mn>4</mn></msub></mtd><mtd><mo>-</mo><mn>1</mn><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>1</mn><mn>2</mn></msubsup><mo>+</mo><mn>2</mn><msubsup><mi>q</mi><mn>4</mn><mn>2</mn></msubsup></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000036.GIF" wi="1615" he="246" /></maths>由于四元数q<sub>1</sub>,q<sub>2</sub>,q<sub>3</sub>,q<sub>4</sub>的平方和等于一,因此旋转矩阵的待优化参数简化为三个;将已求得的左右两个CCD摄像机标定参数值作为优化初始值,利用优化算法对摄像机主参数和畸变参数进行整体优化;对左右两个CCD摄像机主参数和畸变参数的优化是基于最大似然估计准则,对于给定的n个标定板图片提供的m×n个标定点坐标,左右两个CCD摄像机参数的最优化问题通过下式的最小化问题来表达:<maths num="0013" id="cmaths0013"><math><![CDATA[<mrow><mi>S</mi><mrow><mo>(</mo><mi>&theta;</mi><mo>)</mo></mrow><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><msubsup><mi>&epsiv;</mi><mi>i</mi><mn>2</mn></msubsup><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>m</mi></munderover><msup><mrow><mo>[</mo><msub><mi>y</mi><mi>i</mi></msub><mo>-</mo><msub><mi>m</mi><mi>i</mi></msub><mrow><mo>(</mo><msub><mi>C</mi><mi>j</mi></msub><mo>,</mo><msub><mi>p</mi><mi>j</mi></msub><mo>,</mo><msub><mi>X</mi><mi>i</mi></msub><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000041.GIF" wi="1526" he="135" /></maths>其中j是指参与计算的第j个摄像机,i是指第j个摄像机获取的第i个点,X<sub>i</sub>是输入的空间点坐标,y<sub>i</sub>是空间第i点的图像坐标,C<sub>j</sub>是固定不变的摄像机参数向量,其长度为n<sub>0</sub>,p<sub>j</sub>是需要调整的摄像机参数向量,其长度为n<sub>1</sub>,n<sub>0</sub>+n<sub>1</sub>即为摄像机所有参数向量长度,m<sub>i</sub>(C<sub>j</sub>,p<sub>j</sub>,X<sub>i</sub>)是摄像机的成像方程;本发明采用基于LM算法的光束平差法来解决上式的最小化优化问题;利用光束平差法的编写优化程序时除了要提供左右两个CCD摄像机的标定参数初值外还要提供一定数量的和形式的三维空间坐标点以及这些三维坐标点的相应图像坐标,因此为了精确得到整个视场空间内高精密标定板的角点坐标,采用数控机床z轴进给带动显微镜纵向移动的方法;并采用激光干涉仪(测量精度可达到0.1um)对z轴进行测量得到z轴移动的实际距离,保证所采集角点坐标点的z向实际精度;在光束平差法中,优化算法的迭代因子是由迭代参数的一阶导数组成的雅克比矩阵构成的,因此需要在计算过程中引入待成像方程中调整参数的雅克比矩阵;<maths num="0014" id="cmaths0014"><math><![CDATA[<mrow><mi>J</mi><mo>=</mo><mo>[</mo><mfrac><mi>df</mi><mrow><mi>dp</mi><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></mfrac><mo>,</mo><mfrac><mi>df</mi><mrow><mi>dp</mi><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow></mfrac><mo>,</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>&CenterDot;</mo><mo>,</mo><mfrac><mi>df</mi><mrow><mi>dp</mi><mrow><mo>(</mo><msub><mi>n</mi><mn>1</mn></msub><mo>)</mo></mrow></mrow></mfrac><mo>]</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000042.GIF" wi="1412" he="155" /></maths>其中,f为成像方程,p(1)…p(n<sub>1</sub>)是需要调整的摄像机参数向量元素;光束平差法以其待调整参数的雅克比矩阵组成迭代因子,对模型参数进行调整,得到使成像误差S(θ)达到最小的模型参数最优解;(2)待测区域特征点的提取本发明针对待测区域特征点的提取采用的是Harris角点检测算法;Harris提出了用高斯函数代替方形区域计算梯度的算法,若目标点的像素坐标为x,y,其在x方向和y方向的位移分别为u和v,ω(u,v)为窗口函数,则该目标点的灰度变化由如下公式表示:<maths num="0015" id="cmaths0015"><math><![CDATA[<mrow><mfenced open='' close=''><mtable><mtr><mtd><mi>E</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><munder><mi>&Sigma;</mi><mrow><mi>u</mi><mo>,</mo><mi>v</mi></mrow></munder><mi>&omega;</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo></mrow><msup><mrow><mo>[</mo><mi>I</mi><mrow><mo>(</mo><mi>x</mi><mo>+</mo><mi>u</mi><mo>,</mo><mi>y</mi><mo>+</mo><mi>v</mi><mo>)</mo></mrow><mo>-</mo><mi>I</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>]</mo></mrow><mn>2</mn></msup></mtd></mtr><mtr><mtd><mo>&ap;</mo><munder><mi>&Sigma;</mi><mrow><mi>u</mi><mo>,</mo><mi>v</mi></mrow></munder><mi>&omega;</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo></mrow><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo></mrow><mfenced open='[' close=']'><mtable><mtr><mtd><msubsup><mi>I</mi><mi>x</mi><mn>2</mn></msubsup></mtd><mtd><msub><mi>I</mi><mi>x</mi></msub><msub><mi>I</mi><mi>y</mi></msub></mtd></mtr><mtr><mtd><msub><mi>I</mi><mi>x</mi></msub><msub><mi>I</mi><mi>y</mi></msub></mtd><mtd><msubsup><mi>I</mi><mi>y</mi><mn>2</mn></msubsup></mtd></mtr></mtable></mfenced><mfenced open='(' close=')'><mtable><mtr><mtd><mi>u</mi></mtd></mtr><mtr><mtd><mi>v</mi></mtd></mtr></mtable></mfenced></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>14</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000051.GIF" wi="1412" he="347" /></maths>令<maths num="0016" id="cmaths0016"><math><![CDATA[<mrow><mi>M</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>=</mo><munder><mi>&Sigma;</mi><mrow><mi>u</mi><mo>,</mo><mi>v</mi></mrow></munder><mi>&omega;</mi><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo></mrow><mfenced open='[' close=']'><mtable><mtr><mtd><msubsup><mi>I</mi><mi>x</mi><mn>2</mn></msubsup></mtd><mtd><msub><mi>I</mi><mi>x</mi></msub><msub><mi>I</mi><mi>y</mi></msub></mtd></mtr><mtr><mtd><msub><mi>I</mi><mi>x</mi></msub><msub><mi>I</mi><mi>y</mi></msub></mtd><mtd><msubsup><mi>I</mi><mi>y</mi><mn>2</mn></msubsup></mtd></mtr></mtable></mfenced></mrow>]]></math><img file="FDA0000741774530000052.GIF" wi="766" he="200" /></maths>则:<maths num="0017" id="cmaths0017"><math><![CDATA[<mrow><mi>E</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mo>&ap;</mo><mrow><mo>(</mo><mi>u</mi><mo>,</mo><mi>v</mi><mo>)</mo></mrow><mi>M</mi><mrow><mo>(</mo><mi>x</mi><mo>,</mo><mi>y</mi><mo>)</mo></mrow><mfenced open='(' close=')'><mtable><mtr><mtd><mi>u</mi></mtd></mtr><mtr><mtd><mi>v</mi></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>15</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000053.GIF" wi="1413" he="172" /></maths>计算矩阵M(x,y)的特征值λ<sub>1</sub>,λ<sub>2</sub>,如果两个特征值都比较大,说明目标点的图像灰度自相关函数的两个正交方向上的值均较大,则该点即为特征点;通过Harris算法对特征点进行初提取后;采用亚像素角点提取算法提取出精度更高的角点坐标;对于理想角点,它附近的像素点的灰度梯度方向均垂直于该点与理想角点的连线;这一特点可以用公式表达为:<maths num="0018" id="cmaths0018"><math><![CDATA[<mrow><mo>&dtri;</mo><mover><mi>H</mi><mo>&RightArrow;</mo></mover><mrow><mo>(</mo><mover><mi>&alpha;</mi><mo>&RightArrow;</mo></mover><mo>-</mo><mover><mi>&beta;</mi><mo>&RightArrow;</mo></mover><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000054.GIF" wi="1324" he="126" /></maths>其中<img file="FDA0000741774530000055.GIF" wi="90" he="70" />是理想角点的灰度梯度方向,<img file="FDA0000741774530000056.GIF" wi="50" he="61" />向量为图像原点指向理想角点的坐标,<img file="FDA0000741774530000057.GIF" wi="55" he="83" />向量为图像原点指向理想角点附近任一边缘点的坐标;实际上,由于受到图像噪声的影响,公式(10)通常都不为零,可以将其视为误差e,即<maths num="0019" id="cmaths0019"><math><![CDATA[<mrow><mi>e</mi><mo>=</mo><mo>&dtri;</mo><mover><mi>H</mi><mo>&RightArrow;</mo></mover><mrow><mo>(</mo><mover><mi>&alpha;</mi><mo>&RightArrow;</mo></mover><mo>-</mo><mover><mi>&beta;</mi><mo>&RightArrow;</mo></mover><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>17</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000058.GIF" wi="1359" he="141" /></maths>在以角点为中心的邻域内,对所有点按上式计算,误差和为E,则<maths num="0020" id="cmaths0020"><math><![CDATA[<mrow><mi>E</mi><mo>=</mo><munder><mi>&Sigma;</mi><mi>i</mi></munder><mi>e</mi><mo>=</mo><munder><mi>&Sigma;</mi><mi>i</mi></munder><mo>&dtri;</mo><mover><mi>H</mi><mo>&RightArrow;</mo></mover><mrow><mo>(</mo><mover><mi>&alpha;</mi><mo>&RightArrow;</mo></mover><mo>-</mo><mover><mi>&beta;</mi><mo>&RightArrow;</mo></mover><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>18</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000059.GIF" wi="1391" he="139" /></maths>其中,i是邻域内第i个点;这样求使误差和E最小的点即是亚像素级别的角点坐标;(3)特征点的立体匹配采用归一化交叉算法(NCC)实现特征点的初匹配,以CCD摄像机2拍得的图片作为基准图片,并以特征点p<sub>l</sub>为中心,构造一个图像像素N*N,的局部图像块作为模版图T,并令模版图T在待匹配图像S的包含有特征点的一定范围内进行遍历,搜索模版所覆盖的子图记做S<sup>i,j</sup>,i,j为子图中心点在匹配图像S中的像素坐标,通过使用归一化交叉算法(NCC)来计算模版图T和搜索子图S<sup>i,j</sup>之间的归一化互相关系数,算法如下:<maths num="0021" id="cmaths0021"><math><![CDATA[<mrow><mi>R</mi><mrow><mo>(</mo><mi>i</mi><mo>,</mo><mi>j</mi><mo>)</mo></mrow><mo>=</mo><mfrac><mrow><munder><mi>&Sigma;</mi><mrow><mi>m</mi><mo>,</mo><mi>n</mi></mrow></munder><mo>|</mo><msup><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msup><mrow><mo>(</mo><mi>m</mi><mo>,</mo><mi>n</mi><mo>)</mo></mrow><mo>-</mo><mi>E</mi><mrow><mo>(</mo><msup><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msup><mo>)</mo></mrow><mo>|</mo><mo>|</mo><mi>T</mi><mrow><mo>(</mo><mi>m</mi><mo>,</mo><mi>n</mi><mo>)</mo></mrow><mo>-</mo><mi>E</mi><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>|</mo></mrow><msqrt><munder><mi>&Sigma;</mi><mrow><mi>m</mi><mo>,</mo><mi>n</mi></mrow></munder><msup><mrow><mo>(</mo><msup><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msup><mrow><mo>(</mo><mi>m</mi><mo>,</mo><mi>n</mi><mo>)</mo></mrow><mo>-</mo><mi>E</mi><mrow><mo>(</mo><msup><mi>S</mi><mrow><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msup><mo>)</mo></mrow><mo>)</mo></mrow><mn>2</mn></msup><munder><mi>&Sigma;</mi><mrow><mi>m</mi><mo>,</mo><mi>n</mi></mrow></munder><msup><mrow><mo>(</mo><mi>T</mi><mrow><mo>(</mo><mi>m</mi><mo>,</mo><mi>n</mi><mo>)</mo></mrow><mo>-</mo><mi>E</mi><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>)</mo></mrow><mn>2</mn></msup></msqrt></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>19</mn><mo>)</mo></mrow></mrow>]]></math><img file="FDA0000741774530000061.GIF" wi="1650" he="311" /></maths>其中:m,n代表各个像素的坐标,T(m,n)为点m,n处的灰度值,E(S<sup>i,j</sup>)和E(T)分别是搜索子图S<sup>i,j</sup>与模版图T的灰度平均值;互相关系数R(i,j)的值越大,则搜索子图与模板图的匹配程度越高;设定互相关系数R(i,j)的阈值后,选出特征点p<sub>l</sub>的候选匹配点;使用NCC算法进行初匹配后还需要利用外极限约束和距离约束来对匹配的特征点对进行校正;通过使用Longguet‑Higgins提出的归一化8点算法在摄像机标定时计算左右两个CCD摄像机的基本矩阵F,根据外极线理论,对基准图像特征点p<sub>l</sub>,在待匹配图像S中对应的外极线l<sub>r</sub>可以表示为:l<sub>r</sub>=F·p<sub>l</sub>              (20)初匹配的候选特征点若在外极线l<sub>r</sub>附近则可进一步确定该候选特征点为匹配点;最后统计每一候选特征点到其他特征点的距离和基准图像中特征点到其他特征点点距离相等的个数num;num最大值时对应的候选特征点则为基准图像中特征点的正确匹配点;(4)三维坐标的求取首先计算出双目图像中特征点经亚像素算法提取出的图像坐标值,并大致估算出特征点在世界坐标系中的估计值;在优化程序中将左右两个CCD摄像机参数对应雅克比矩阵设为零,同时在优化程序计算迭代因子部分中嵌入关于空间坐标x、y、z的雅克比矩阵方程;将已经优化好的左右两个CCD摄像机参数以及特征点点图像坐标值和世界坐标估计值输入光束平差法优化程序,就可以对空间点的三维坐标进行优化,得到特征点的三维空间测量值,从而完成三维坐标点的求取并最终实现特征点的定位。
地址 116100 辽宁省大连市凌工路2号
您可能感兴趣的专利