一種基于Delaunay剖分三角形內插值的攝像機標定方法

徐 侃，李文龍，尹周平

（華中科技大學 數字制造裝備與技術國家重點實驗室，湖北 武漢 430074）

視覺定位因其具有非接觸、高精度、高速度等特點，在RFID標簽封裝機、裝片機、劃片機、引線鍵合機電子封裝設備上得到廣泛應用[1~2]。攝像機標定是視覺定位的重要環節之一，用于校正圖像采集、傳輸及處理過程中引入的變形，并建立由圖像坐標到世界坐標的轉換關系[3]。目前攝像機標定方法大多基于三維測量原理，從不同位置拍攝的多幅標定圖像計算攝像機內、外參數[4~5]。電子封裝設備標定對象為平面距離和角度，要求盡可能減少人工干預，提升在線運行效率。因此，開發操作方便、計算復雜度小且精度滿足電子封裝設備運行需要的標定算法，具有重要的理論意義和實用價值。

針對在線攝像機標定中減少人工干預的需求，國內外研究者提出一種基于多項式擬合的攝像機標定方法[6~7]。該方法基于小孔成像模型，將常見圖像變形（平移、旋轉、非線性畸變）的數學模型，以一個二元五次多項式表示。通過帶入變形前后特征點的像素坐標擬合計算得到上述多項式的系數，從而完成標定計算。該方法只用一幅圖像即可完成電子封裝設備所需的標定計算且精度滿足要求，因此，得到了普遍應用。但該方法的計算過程涉及到了高階（21階）矩陣的求逆運算，計算復雜度較高且易出現奇異情況（矩陣的逆不存在）。

本文提出一種基于Delaunay剖分三角形內插值的攝像機標定方法，簡稱CIDT(Calibration based on the Interpolation ofDelaunay Triangles)方法。該方法利用了特征點網格Delaunay剖分三角形的性質以及有效剖分三角形區域的內插值線性關系不變性，使用計算復雜度更低的三角形內插值算法取代原有的多項式擬合算法。在不降低標定精度的前提下，簡化了標定計算，降低了標定算法實現的難度。

1 標定特征點信息提取

1.1 特征點圖像坐標的提取

目前，常用的標定板靶標有棋盤格型與圓點陣列型。棋盤格型靶標采用棋盤格四周的角點作為特征點，這類靶標制作簡單、適用性廣泛。對于這類靶標特征點的提取，可以采用基于灰度值相關的Harris角點檢測算法，通過求取特定區域的梯度相對于相鄰區域的梯度均值的變化規律以及相鄰區域之間梯度分布的相關性來求取CIM (Cornersof Image Marking)[8]，最后根據CIM中大于指定閾值的局部極大值來確定角點[9]。

圓點陣列型靶標采用圓形區域中心作為特征點，具有更高的標定精度，在IC封裝設備上被廣泛應用。對于這類靶標特征點的提取，可以采用區域輪廓橢圓擬合算法，通過提取特征點區域的輪廓點以橢圓方程作為目標式進行最小二乘擬合得到其對應的橢圓方程，方程對應的橢圓中心即為該特征點區域對應的中心坐標[10]。

1.2 特征點網格坐標的計算

CIDT方法對標定靶標的成像有如下三個要求：第一，特征點盡可能地充滿整個圖像且特征點網格中每一行的點個數相等；第二，特征點網格中每一行特征點的圖像行坐標范圍與其它行沒有交集；第三，用戶應輸入特征點網格的行數R 以及列數C。

在滿足上述要求的前提下，假設特征點集P 中每個特征點Pi（i 為特征點坐標提取時的順序序號）的像素坐標為(xi，yi)，按照如下步驟完成P 中每個特征點像素坐標與記錄了其在特征點網格上所處位置的網格坐標之間的對應：

步驟1：根據像素坐標(x，y)中的行坐標y 對P 中的所有特征點按照從小到大的順序排序，得到的排序后特征點集記作P'。

步驟2：將P'中的特征點分為R個小組，每個小組C個點。在每個小組內根據像素坐標(x，y)中的列坐標x 對組內的特征點按照從小到大的順序排序，得到的再次排序后特征點集記作P"。

步驟3：對特征點集P"中的每個特征點Pj（j 為經過兩次排序后的特征點順序序號），按照式（1）計算得到其網格坐標(rj，cj)。其中，符號int表示向下取整運算，符號%表示取余運算。

2 標定計算

2.1 有效剖分三角形區域的內插值線性關系不變性

圖像采集系統的光學畸變會給標定圖像帶來非線性變形，其中，攝像機鏡頭的一階徑向畸變對圖像中的非線性變形起到了主要作用，可以近似地認為圖像非線性變形全部來自于攝像機鏡頭的一階徑向畸變[12~14]，其數學模型，如式（2）所示。

式中，(xu，yu)、(xd，yd)分別對應畸變前、后的像素坐標，(x0，y0)是畸變主點的像素坐標，z 為圖像采集系統的一階徑向畸變系數。在完成像素坐標的提取后，按照圖像坐標對特征點集P 進行Delaunay三角剖分，可以得到剖分三角形集合τ。由于圖像非線性畸變的影響，τ 可能會包含對于本文所述標定計算起到干擾作用奇異剖分三角形，剔除這些三角形后得到的剩余有效剖分三角形集合存在一些有用性質。

2.1.1 有效剖分三角形的定義

對于有特征點集合按照圖像坐標進行Delaunay三角剖分得到的剖分三角形集合τ，定義其中“頂點由一個特征點及其在特征點網格上的兩個四鄰域點組成的三角形”為有效剖分三角形。例如，圖1（a）顯示了一個3×3網格點集的Delaunay剖分三角形集合，圖1（b）顯示了其中的有效剖分三角形集合。

圖1 有效剖分三角形定義示例

2.1.2 有效剖分三角形的建立

按照式（3）定義一個剖分三角形的形狀偏差角η，其中，ω1、ω2、ω3分別為該三角形三個內角。

根據實際圖像非線性變形程度設定合適的形狀偏差角閾值η0（非線性變形越大，該閾值取得越高，取值范圍一般在5°~25°），從三角形集合τ 中剔除其中η>η0的三角形可以得到有效剖分三角形集合τe。

2.1.3 有效剖分三角形的內插值精度分析

在大小為512×512的圖像上，以像素點(64，64)為起點，以32為像素間距，均勻陣列分布著的13×13標準陣列點集合，記作ξ。對其進行Delaunay三角剖分并按照上述方法建立有效剖分三角形集合τ，如圖2（a）所示。對τ 中的每個三角形τk（k 表示有效剖分三角形的序號）隨機生成一組三角形內插值系數mk、nk（mk、nk應滿足條件：0

圖2 CDIT方法精度分析實驗用仿真點集

按照式（2）對點集合ξ 和點集合ε 加入一階徑向畸變系數z=1×10-6的非線性畸變，得到變形后的點集合ξ'和ε'，由ξ'亦可得到變形后的有效剖分三角形集合τ'，結果如圖2（c）所示。利用三角形集合τ'中的每個三角形τk'的頂點圖像坐標按照與之前相同的插值系數mk、nk帶入式（4）計算可以得到點集合ε"，結果如圖2（d）所示。對點集合ε'和點集合ε"中每對位于同一有效剖分三角形內的點εk'與εk"，計算它們在圖像上的像素距離，結果如表1所示（點個數太多，本文僅顯示其中的一部分）。

表1 點集合ε'和ε"中對應點的像素距離

從表1可以發現，在畸變條件下，利用變形前像素點與其所在的有效剖分三角形頂點的內插值關系，由有效剖分三角形頂點插值估算得到的像素點位置與實際的像素點位置相差不大。對在大小為512×512的圖像上以像素點(64，64)為起點的不同規格的標準陣列網格點集，分別加入不同程度的一階徑向畸變并按照上述方法計算實際像素點位置與內插值估算位置的像素距離，結果如圖3所示。

圖3 不同情況下實際位置與內插值估算位置的像素距離

可以看出，畸變程度越小，有效剖分三角形的面積越小（網格越密集），利用內插值估算出的像素點位置越精確。對于一副512×512圖像而言，當陣列網格點的密度達到使得有效剖分三角形面積不大于250個像素時，即便畸變系數已經高達1.5×106（遠遠超過測量用圖像采集系統最大畸變不超過1%的要求），這種估算的誤差已經只有不到0.1個像素。

在實際應用過程中，只要陣列網格點足夠密集，就可以認為由變形后的有效剖分三角形頂點按照變形前像素點與其所在的有效剖分三角形頂點的內插值關系估算得到的變形后像素點位置就是該像素點變形后實際的位置。這種在一定條件下，變形前后圖像素點與其所在的有效剖分三角形頂點之間的內插值關系保持不變的性質，稱作有效剖分三角形區域的內插值線性關系不變性。

2.2 標定文件的制作

根據2.1節中討論的有效剖分三角形區域的內插值線性關系不變性，可以利用變形前每個像素點與其所在有效剖分三角形頂點的內插值關系，計算變形前每個圖像像素點對應在變形后圖像上的位置。將該對應關系以及變形前特征點的像素間距與特征點的物理間距存儲起來，完成標定文件的制作。

2.2.1 變形前圖像上特征點位置的估算

d=max（（H+2）/（R-1），（W+2）/（C-1））（5）

設標定圖像的高、寬分別為H、W，按照式（5）估算變形前特征點之間的像素間距d。并根據由用戶輸入的標定特征點物理間距D，計算得到坐標單位轉換系數T=D/d。根據像素間距d，按照式（6）估算特征點網格最左上角點Plu在變形前對應點（記作Qlu）的圖像坐標（XXlu，YYlu）。

根據點Qlu的圖像坐標（XXlu，YYlu）與每個特征點Pi的網格坐標(ri，ci)，按照式（7）計算變形前每個特征點Qi（此處i 的含義與1.2節提到的一致）的圖像坐標（xxi，yyi）。按照這種方法估算得到的變形前特征點集，可以保證變形前圖像每個像素點都被一個有效剖分三角形所包含。

2.2.2 內插值系數的計算

將有效剖分三角形集合τe中每個三角形△k（此處k 的含義2.1.3節一致）的頂點Pk1、Pk2、Pk3根據其對應的特征點序號i 按照式（1）計算其網格坐標，帶入式（7）計算得到其對應的變形前點Qk1、Qk2、Qk3的圖像坐標。從而得到△k對應的變形前三角形△k'。由△k'組成的三角形集合稱為變形前的有效剖分三角形集合，記作τe'。

設變形前的標定圖像為I'，遍歷I'的每一個像素點Bl（l 為圖像像素點的順序序號），在τe'中找到并記錄包含Bl的三角形△k'（τe'中有且僅有一個三角形包含Bl）。設△k'頂點Qk1、Qk2、Qk3的圖像坐標分別為（xxk1，yyk1）、（xxk2，yyk2）、（xxk3，yyk3）,解 式（8）所 示 的二元一次方程組，得到三角形內插值系數ml、nl并分別存入WH×1矩陣M、N 的第l行中，其中（xxl，yyl）為像素點Bl的像素坐標。

2.2.3 變形前后像素點對應關系計算

遍歷I'的每一個像素點Bl，根據Bl所在的三角形△k'的序號k 在三角形集合τe中找到三角形△k，設其頂點Pk1、Pk2、Pk3的圖像坐標分別為（xk1，yk1）、（xk2，yk2）、（xk3，xk3）；根據像素點序號l在矩陣M、N 中三角形內插值系數ml、nl。按照式（9）計算Bl在變形后圖像I 中對應位置的圖像坐標（xl，yl）。從而構建I'中每個像素點Bl對應在I 中對應位置的圖像坐標映射關系“（xxl，yyl）→（xl，yl）”。

最后，將上述坐標映射關系與坐標單位轉換系數T 一起存儲到硬盤上，便完成了標定文件，便完成了標定文件的制作。綜上所述標定文件的制作方法可以用如圖4所示流程圖來表示。

2.3 圖像校正與坐標轉換的實現

圖4 標定文件的制作方法流程圖

載入待處理圖像II（II 與標定圖像I 由相同的攝像機在相同條件下拍攝得到，其高度、寬度也分別為H、W）以及標定文件，創建一副高度H、寬度W 的空圖像II'，用于存放校正后的圖像。遍歷II'中每個像素點Fl'，按照標定文件中存儲的對應關系找到II'中每個像素點Ql'（xl'，yl'）變形后對應于圖像II 的位置圖像坐標（xl，yl）。用（xl，yl）的臨近像素點Qllu（int(xl，)，int(yl)）、Qlld（int(xl,)+1，int(yl)）、Qlru（int(xl,),int(yl)+1）、Qlrd（int(xl,)+1,int(yl)+1）的像素值pix(Qllu)、pix(Qlld)、pix(Qlru)、pix(Qlrd)（函數int()表示向下取整操作），雙線性插值法得到像素點Ql'的像素值pix(Ql')。遍歷結束后，II'中儲存的就是圖像II 通過校正后的理想圖像。

讀取標定文件中儲存的坐標單位轉換系數T，對于在II'計算得到的特征尺寸像素值，都可以通過乘以該系數的方式轉換為世界坐標值。

3 實驗測試

圖5（a）顯示了一副存在較大非線性畸變的標定板圖像，圖5（b）顯示了一副與圖5（a）在相同條件（相同的攝像機、鏡頭、焦距、物距等）下拍攝得到的電路板圖像。

圖5 實驗測試圖像

已知上述標定圖像中特征點間距的世界坐標值D=1mm，提取其中的特征點像素坐標并對特征點集合進行Delaunay三角剖分，以形狀偏差角閾值η0=15°篩選有效剖分三角形，結果如圖6所示。

圖6 特征點集的有效剖分三角形集合

按照CIDT方法完成標定文件的制作并載入該標定文件對圖5（a）、圖5（b）所示的圖像進行去畸變校正，結果分別如圖7（a）、圖7（b）所示。另外，利用圖5（a）所示的標定圖像采用成熟商用軟件包MIL制作標定文件對圖5（a）、圖5（b）所示的圖像進行去畸變校正，結果分別如圖7（c）、圖7（d）所示。可以看出，CIDT方法能夠對存在較大程度畸變的圖像達到良好的去畸變校正效果。由于在估計變形前的特征點分布時需要保證所有像素點都至少在一個有效剖分三角形內部，故CIDT方法校正圖像會帶來部分圖像邊緣區域信息損失的問題。一般而言，需要測量、定位的特征常位于圖像中心，上述圖像邊緣區域信息損失可以忽略不計。

圖7 MIL與CIDT方法分別得到的去畸變校正圖像

采用CIDT方法與MIL軟件包校正所得到的圖7（b）與圖7（d）中選擇的3個特征尺寸，分別如圖8（a）、8（b）所示。

圖8 測試用特征尺寸

對這些特征尺寸都采用MIL軟件包求取像素距離，并分別經過CIDT方法與MIL制作的標定文件轉換為物理尺寸，結果見表2。可以發現，使用CIDT方法轉換得到的特征尺寸物理值與通過MIL軟件包轉換得到結果差距可以保持在0.2%內。

表2 特征尺寸測試結果

4 結束語

（1）本文通過分析特征點集Delaunay剖分三角形的性質，定義了有效剖分三角形并提出了其篩選方法。同時，分析了有效剖分三角形區域的內插值線性關系不變性，利用這一性質提出一種新的平面標定方法，該方法采用分塊線性插值代替了傳統方法需要采用二元高次多項式擬合才能完成的變形前像素點與其在變形后圖像上對應位置的映射運算，從而簡化標定計算。

（2）測試結果表明，本文所述方法在應用于圖像去畸變校正時，對于超出一般測量要求的非線性圖像畸變仍然能有效去除。同時，該方法在校正圖像的過程中也帶來圖像邊緣信息損失的問題，當需要測量、定位的特征位于圖像中心區域時，這一問題帶來的影響不大。

（3）測試結果表明，本文所述方法在應用于圖像坐標向世界坐標的轉換時，精度與成熟商用軟件包MIL相當，滿足電子封裝設備的精度要求。

[1]Li Junlan，Zhang Dawei，Zhao Xingyu,etal.Micro-vision positioning systems for IC packaging[J].Optics and Precision Engineering,2010,18(4):995-972.

[2]Zhang Kai,Wang Yuhui,Yin Zhoupin.A fast vision alignment in high performance chipmounter[J].Optical Technique，2005，31(4):604-607.

[3]Qiu Maolin,Ma Desong,Li Yi.Overview of Camera calibration forcomputervision[J].Acta Automatica Sinica,2000,26(1):43-55.

[4]Zhang Guangjun.Machine vision[M].Beijing：Science Press,2005.

[5]Zhao Xiaosong,Zhang Hongwei,Zhang Guoxiong,et al.Study on camera calibration technology[J].Chinese Journal of MechanicalEngineering,2002,38(3):149-151.

[6]Liao Shizhong,Gao Peihuan,Su Yi,etal.A geometric rectificationmethod for lens camera[J].Journal of Image and Graphics,2000,5(7):593-596.

[7]He Xiaolan,Jiang Guoquan,Du Shangfeng.Camera calibration arithmetic based on integration ofmultinomial[J].J.Cent.South Univ.(Science and Technology),2007,38(1):1117-1122.

[8]LOWE D G.Distinctive image features from scale-invariant key-points[J].Pattern Recognition and Artificial Intelligence,2009,22(4):503-513.

[9]Liu Ning,Lu Shengrong,Xia Ruixue,et al.Sub-pixel Harris corner detection algorithm based on Gaussian surfacemodel[J].Electronic MeasurementTechnology,2011,34(12):49-53.

[10]Xia Ruixue,Lu Shengrong,Liu Ning.A Method of Automatic Extracting Feature Point Coordinates Based on Circle Array Target[J].China Mechanical Engineering,2010,21(16):1906-1910.