(lunwen)數(shù)碼相機定位的數(shù)學(xué)模型.doc

數(shù)碼相機定位的數(shù)學(xué)模型摘 要隨著數(shù)碼相機定位在各領(lǐng)域的廣泛應(yīng)用，對相關(guān)問題機器視覺的研究也成為熱點。因此建立一個精度較高，穩(wěn)定性好的數(shù)碼相機定位的數(shù)學(xué)模型，具有很好的現(xiàn)實意義。問題1要求給出確定靶標(biāo)上圓的圓心在給定相機像平面的像坐標(biāo)的算法，問題2利用問題1的模型對給定數(shù)據(jù)求解。為此，首先建立了四個空間直角坐標(biāo)系，在MATLAB中把圖3的數(shù)字信息提取出來，主要是五個橢圓的邊緣點的信息；同時為了便于運算，通過坐標(biāo)變換將計算機圖像坐標(biāo)變換為圖像坐標(biāo)；并用提取的圖像邊界坐標(biāo)擬合出5個橢圓的方程，利用“曲線切線的投影仍與曲線的投影相切，而且切點的投影仍為投影的切點”這一引理，提取出靶標(biāo)上圓及其像上的公切點的坐標(biāo)作為特征點，利用RAC兩步法標(biāo)定過程和最小二乘法建立了計算世界坐標(biāo)系到相機坐標(biāo)系的旋轉(zhuǎn)變換矩陣和平移向量及徑向畸變系數(shù)的算法。利用16個公切點作為特征點，通過Matalb編程求得靶標(biāo)上圓的圓心在文中給定相機像平面的五個坐標(biāo)(單位：mm)：A(-49.7132, 51.1289 417.1958)，B(-23.3475, 49.1539 417.1958)，C(33.8194, 44.8716, 417.1958)， D(18.8173,-31.5798, 417.1958)，E(-59.7830, -31.1754, 417.1958)。問題3的解決分為兩步：一是通過對模型計算出的焦距及畸變系數(shù)及上面五個坐標(biāo)值的分析得出模型的精度較高的結(jié)論；二是采用改變特征點數(shù)的方法或利用“,三個標(biāo)靶的中心的像應(yīng)在一條直線上”驗證模型的穩(wěn)定性。問題4采用二目立體視覺模型確定了給出兩部固定相機相對位置的數(shù)學(xué)模型和方法。本文建立的算法可操作性強，精度較高，穩(wěn)定性好，對解決類似問題的計算有一定的推廣價值。關(guān)鍵詞：擬合橢圓 特征點提取 RAC兩步法 坐標(biāo)旋轉(zhuǎn)矩陣公切點- 23 -數(shù)碼相機定位的數(shù)學(xué)模型一問題的提出數(shù)碼相機定位的數(shù)學(xué)模型來源于2008年全國大學(xué)生數(shù)學(xué)建模競賽的A題。一般地在物平面上畫若干個圓（稱為靶標(biāo)），它們的圓心就是幾何的點了，它們的像一般會變形為橢圓，從靶標(biāo)上的這些圓的像中把圓心的像精確地找到，標(biāo)定就可實現(xiàn)?，F(xiàn)設(shè)計靶標(biāo)如下，取1個邊長為100mm的正方形，分別以四個頂點（對應(yīng)為A、C、D、E）為圓心，12mm為半徑作圓。以AC邊上距離A點30mm處的B為圓心，12mm為半徑作圓，如圖1所示，用一位置固定的數(shù)碼相機攝得其像，如圖2所示。圖1 靶標(biāo)示意圖圖2 靶標(biāo)的像(1) 建立數(shù)學(xué)模型和算法以確定靶標(biāo)上圓的圓心在該相機像平面的像坐標(biāo), 這里坐標(biāo)系原點取在該相機的光學(xué)中心，x-y平面平行于像平面；(2) 對由圖2、圖3分別給出的靶標(biāo)及其像，計算靶標(biāo)上圓的圓心在像平面上的像坐標(biāo), 該相機的像距（即光學(xué)中心到像平面的距離）是1577個像素單位(1毫米約為3.78個像素單位)，相機分辨率為1024768；(3) 設(shè)計一種方法檢驗?zāi)銈兊哪Ｐ?，并對方法的精度和穩(wěn)定性進行討論；(4) 建立用此靶標(biāo)給出兩部固定相機相對位置的數(shù)學(xué)模型和方法。二模型的假設(shè)及符號說明1模型的假設(shè)(1) 假設(shè)題目中給出的圖3的尺寸是實際物理尺寸；(2) 圖像坐標(biāo)系的坐標(biāo)原點是圖3的中心，稱為主心坐標(biāo)；(3) 相機不需做任何的運動；(4) 標(biāo)定物為共面點，將靶標(biāo)所在的面作為世界坐標(biāo)系的坐標(biāo)平面；(5) 相機的有效焦距不變。2. 符號說明(1) ：表示世界坐標(biāo)系到相機坐標(biāo)系旋轉(zhuǎn)變換矩陣，；(2) ：表示世界坐標(biāo)系到相機坐標(biāo)系的平移向量，；(3) ：表示相機的有效焦距，；(4) ：表示世界坐標(biāo)系下物點的坐標(biāo)；(5) ：表示圖像坐標(biāo)系下物點有徑向畸變的實際像坐標(biāo)；(6) ：表示圖像坐標(biāo)系下物點的針孔成像的理想坐標(biāo)；(7) ：表示相機坐標(biāo)系下物點的坐標(biāo)；(8) ：主點坐標(biāo)，即圖像坐標(biāo)系下坐標(biāo)原點的像素坐標(biāo)；(9) ：表示計算機坐標(biāo)系下的像的像素坐標(biāo)。(10) ：表示徑向畸變系數(shù)。(11) 、：每個像素在軸與軸方向上的物理尺寸。三問題的分析為了確定靶標(biāo)上圓的圓心在文中給定相機像平面的像坐標(biāo)，要把圖3中的數(shù)據(jù)信息提取出來。由于圖3中的數(shù)據(jù)信息是以像素為單位的，為了與圖2的毫米單位一致，同時便于運算，通過坐標(biāo)變換，將計算機圖像坐標(biāo)變換為圖像坐標(biāo)，并用圖像邊界坐標(biāo)的信息擬合出5個橢圓的方程，從而可以求出橢圓的幾何中心。1坐標(biāo)系的建立計算過程中需要建立如下的四個坐標(biāo)系，如圖4所示。(1) 相機坐標(biāo)系：原點定義在相機的光學(xué)中心，軸與光軸重合；(2) 圖像坐標(biāo)系：原點（主點）定義為相機光軸與圖像平面的交點，軸和軸與，軸平行，為相機的有效焦距；(3) 計算機圖像坐標(biāo)系：原點位于CCD圖像平面的左上角，軸和軸分別與軸和軸平行。、分別表示該像素在數(shù)組中的列數(shù)和行數(shù)，且以像素為單位。在計算機圖像坐標(biāo)系中，主點坐標(biāo)記為，假設(shè)主點坐標(biāo)是圖像坐標(biāo)的中點，本文中可表示為。(4) 世界坐標(biāo)系：原點位于靶標(biāo)所在正方形的中心，和軸分別平行、（軸的正向為向量的正方向），因此物點的坐標(biāo)可表示為。OO2圖4 考慮畸變的徑向攝像機模型圖設(shè)為物點在小孔成像的理想圖像坐標(biāo)中的像點坐標(biāo)，是由透鏡徑向畸變引起的物點實際圖像點，考慮畸變的徑向攝像機模型圖，如圖4所示。一般情況下，CCD鏡頭畸變主要為一階徑向畸變，在這里我們只考慮一階徑向畸變。用一個二階多項式近似： 式中為畸變系數(shù)，為畸變坐標(biāo)。2提取圖像數(shù)據(jù)(1) 提取圖3并以.bmp格式保存。(2) 在MATLAB中用imread命令讀入該圖片并轉(zhuǎn)換成RGB圖像矩陣。(3) 使用rgb2gray和im2bw命令將上述圖像矩陣轉(zhuǎn)化成二值圖像矩陣。(4) 使用bwlabel命令對二值圖像矩陣進行分塊標(biāo)識（每個橢圓域標(biāo)記為一塊，共5塊）。(5) 對每一塊橢圓域用edge命令提取邊緣點，并求出這些邊緣點在計算機圖像坐標(biāo)系下的坐標(biāo)。3計算機圖像坐標(biāo)與圖像坐標(biāo)的轉(zhuǎn)換如圖4所示，計算機圖像坐標(biāo)系以像素為單位，為便于計算，需將其轉(zhuǎn)換為圖像坐標(biāo)系。設(shè)坐標(biāo)系原點在坐標(biāo)系中坐標(biāo)為，每個像素在軸與軸方向上的物理尺寸為、，則圖像中任意一個像素在兩個坐標(biāo)系下的坐標(biāo)有如下關(guān)系：(1)用齊次坐標(biāo)與矩陣形成，將式(1)表示為逆關(guān)系可寫成(2)已知相機分辨率為，若題中給定圖3尺寸為實際物理尺寸，則 ， 考慮到題中給出，為方便計算，取，即1毫米均表示3.78個像素單位。設(shè)取在靶標(biāo)像平面的中心位置，根據(jù)式(2)，可將像平面五個橢圓邊緣點坐標(biāo)變換到圖像坐標(biāo)系下的坐標(biāo)。4橢圓擬合1橢圓擬合是對已提取的圖3的五個橢圓邊緣點的數(shù)據(jù)，進行曲線擬合。設(shè)橢圓曲線的一般表達式為其橢圓的中心點可由式(3)求得：，(3)為了用數(shù)據(jù)作曲線擬合，首先設(shè)，利用線性最小二乘法擬合，用MATLAB編程得到五個擬合橢圓曲線方程系數(shù)如表1所示，橢圓A的擬合曲線如圖5所示。表1 擬合橢圓曲線方程系數(shù)abcdegA1.873110-4-1.701610-51.935610-40.0196-0.02071.0000B3.276210-4-4.568310-53.285910-40.0177-0.03361.0000C3.886810-4-9.845710-53.631410-4-0.0219-0.02951.0000D6.273510-4-2.388710-56.926410-4-0.03110.04821.0000E2.341610-4-5.376410-52.875410-40.02650.01471.0000圖5橢圓A的擬合曲線將橢圓方程標(biāo)準(zhǔn)化，由式(3)得到A、B、C、D、E五個靶標(biāo)的像的幾何中心點的坐標(biāo)如表2所示。表2 五個靶標(biāo)的像的幾何中心點的坐標(biāo)ABCDE-49.9705-23.493533.886818.7601-60.092651.402449.430345.1492-31.5326-31.2101Z=1577/3.78 單位：mm5用橢圓的公切線提取特征點(1) 引理：曲線切線的投影仍與曲線的投影相切，而且切點的投影仍為投影的切點。3(2)兩橢圓公切線的計算設(shè)像平面上任意兩個橢圓（基于圖像坐標(biāo)系）的曲線方程為：(4)(5)其公切線方程為：(6)將式(6)代入式(4)并整理得：(7)其中，根據(jù)切線的判別法則：，得：(8)其中， ，同理，將式(6)代入式(5)，可得與式(8)形式相同的結(jié)果：(9)其中， ，錯誤求解(8)、(9)聯(lián)立的方程組，可得出四組，即確立了任意兩個橢圓的四條公切線。對于題中給出的靶標(biāo)，我們僅提取出如圖6所示的四條外圍公切線的投影（如圖7所示），并計算出八個公切點的坐標(biāo)。圖7對應(yīng)的MATLAB求解如圖8所示。圖7EDCBA(8)(16)(7)(15)(14)(13)(6)(5)(11)(3)(12)(4)(1)(2)(9)(10)ABCCYwXwDEOw圖6圖8圖6對應(yīng)的MATLAB求解四條外圍公切線的投影方程分別為：: y=8.8393x+590.9984: y=-0.0840x+58.1670: y=4.8613x-166.4789: y=0.0040x-40.3090四條內(nèi)圍公切線的投影方程分別為：: y=7.5705x+345.7467: y=-0.0650x+37.2174: y=5.2918x-83.1902: y=-0.0121x-22.577716個公切點及其對應(yīng)像的坐標(biāo)如表3所示。表316個公切點及其對應(yīng)像的坐標(biāo)序號(1)(2)(3)(4)原坐標(biāo)(-50，-62)(-62，-50)(-62，50)(-50，62)像坐標(biāo)(-61.1197，-40.5540)(-70.3928，-31.2242)(-60.9633，52.1263)(-48.5424，62.2443)序號(5)(6)(7)(8)原坐標(biāo)(50，62)(62，50)(62，-50)(50，-62)像坐標(biāo)(35.9167，55.1501)(43.3636，44.3244)(27.7327，-31.6618)(17.1397，-40.2403)序號(9)(10)(11)(12)原坐標(biāo)(-38，-50)(-50，-38)(-50，38)(-38，50)像坐標(biāo)(-49.8120，-31.3566)(-58.8831，-21.8669)(-51.1894，40.5449)(-39.0025，50.4772)序號(13)(14)(15)(16)原坐標(biāo)(38，50)(50，38)(50，-38)(38，-50)像坐標(biāo)(24.3757，45.8002)(32.0257，35.1357)(20.5297，-22.8256)(9.7596，-31.5445)四、模型的建立與求解1問題(1)的求解算法為了數(shù)學(xué)模型和算法以確定靶標(biāo)上圓的圓心在該相機像平面的像坐標(biāo)，采用zhang提出的RAC兩步法標(biāo)定過程，由于此方法考慮了鏡頭畸變，使精度有了一定的提高。此方法采用徑向排列約束關(guān)系求出世界的坐標(biāo)系到相機坐標(biāo)系得旋轉(zhuǎn)矩陣和平移向量，然后求得利用從而求出世界的坐標(biāo)系與相機坐標(biāo)系之間的對應(yīng)關(guān)系，同時又求得鏡頭畸變系數(shù)。設(shè)表示在已建立的世界坐標(biāo)系下物點的坐標(biāo)，表示物點在相機坐標(biāo)系下的像坐標(biāo)。則世界坐標(biāo)系到相機坐標(biāo)系之間的變換關(guān)系為：(10)其中，表示旋轉(zhuǎn)矩陣，表示平移向量。下面采用徑向約束的兩步法（Two-stage）進行計算：第一步：計算和1)過渡參數(shù)的求解利用徑向排列約束原理及公式(10)得到：(11)由于，式(11)可表示為(12)即:(13)若對每一個特征點來說，知道了其空間坐標(biāo)和相應(yīng)于圖像坐標(biāo)系下坐標(biāo)，就有一個方程(13)與之對應(yīng)，取5個以上的點，式(12)即成為超定方程組，利用最小二乘法求解得到和。利用最小二乘法求解式(13)的具體方法：取個物點，和個對應(yīng)的圖像坐標(biāo)點，得超定方程組為：其中，則的最小二乘估計為。記：。2) 的求解求出上述過度參數(shù)后，就可以根據(jù)R正交性求出的大小。當(dāng)參數(shù)，不兩兩同時為0時：B=否則為其余兩個參數(shù)的平方和的倒數(shù)。3）判斷的符號在運用徑向平行約束時，包括兩平行向量同向和反向的兩種情況，故的符號有兩種可能。可采用下列方法判斷：先設(shè)為正，由過渡參數(shù)求出，由此將標(biāo)定點再投影到圖像平面上。計算出對應(yīng)得圖像坐標(biāo)，比較兩者的符號，如果與同號，那么為正；否則為負(fù)，則上面求出的其余外部參數(shù)都作相應(yīng)的改變。4）根據(jù)R的正交性，計算R的其余參數(shù)（如果的符號為正，則在前加負(fù)號）如果下面計算的為負(fù)號，則的符號與上面相反。4）過渡焦距，過渡徑向畸變系數(shù)和過渡平移分量的求解由于經(jīng)整理可得如下含未知數(shù)，得線性方程組：(14)式中，。用最小二乘法聯(lián)合解(14)這個線形方程組，即可得到，的解。進一步判斷的符號，如果為負(fù)，則將反號。進一步解得，其余參數(shù)符號按上面3)所述作相應(yīng)得改變。至此，與，已經(jīng)求出。第二步：計算平移分量和徑向畸變系數(shù)由于，與已知，根據(jù)(14)式可以得到：(15)式中，。用最小二乘法聯(lián)合解(15)這個線形方程組，即可得到，的解。以上即為兩步法的完整的求解過程。這樣通過式(10)得到物點在相機坐標(biāo)系下坐標(biāo)，再利用關(guān)系式(16)得到物點的像在相機坐標(biāo)系下坐標(biāo)，問題(1)得到解決。2問題（2）的求解(1) 利用上述方法計算靶標(biāo)上圓的圓心在像平面上的像坐標(biāo)在前面利用橢圓的公切線，已經(jīng)找出的8個公切點的世界坐標(biāo)和圖像坐標(biāo)如下：序號(1)(2)(3)(4)原坐標(biāo)(-50，-62)(-62，-50)(-62，50)(-50，62)像坐標(biāo)(-61.1197，-40.5540)(-70.3928，-31.2242)(-60.9633，52.1263)(-48.5424，62.2443)序號(5)(6)(7)(8)原坐標(biāo)(50，62)(62，50)(62，-50)(50，-62)像坐標(biāo)(35.9167，55.1501)(43.3636，44.3244)(27.7327，-31.6618)(17.1397，-40.2403)把它們作為特征點代入公式(13)、(15)得到超定方程組，利用MATLAB編程求解超定方程組得到 和 最后利用公式 和，得到五個靶標(biāo)的中心點的像在相機坐標(biāo)系下坐標(biāo)如表4所示。 表4 利用8個特征點求得五個靶標(biāo)的中心點的像在相機坐標(biāo)系下坐標(biāo)ABCDEX-49.8498-23.400333.940518.9014-59.9297Y51.286149.299744.9934-31.6823-31.2713Z=1577/3.78 單位：mm(2) 計算結(jié)果的分析 的值：通過模型(14)計算的過渡值為，而實際的給出的，值的相對誤差：,說明了計算方法的可靠性。 的值：通過模型(14)計算徑向畸變系數(shù)的值為,說明在此問題的求解過程中可以忽略徑向畸變帶來的誤差。 靶標(biāo)中心點的像的坐標(biāo)：通過對表2和表4的比較可以看出，靶標(biāo)的像的幾何中心和靶標(biāo)的像的中心坐標(biāo)并不一致，但誤差并不是很大，說明如果要求精度不高時，也可以用靶標(biāo)的像的幾何中心代替靶標(biāo)的像的中心。3問題(3)：模型的檢驗 通過上面的計算結(jié)果的分析可以看出算法的精度較高。下面用以下兩種方法檢驗上述模型的穩(wěn)定性(1) 利用改變特征點數(shù)的方法驗證模型的穩(wěn)定性在計算靶標(biāo)上圓的圓心在像平面上的像坐標(biāo)時，前面利用8個公切點作為特征點計算了結(jié)果?，F(xiàn)在加大和減少特征點的個數(shù)計算模型，如果加大和減少特征點對計算結(jié)果的影響不大說明，我們建立的模型是穩(wěn)定的。通過對靶標(biāo)上的的16個外公切點作為特征點，計算靶標(biāo)上圓的圓心在像平面上的像坐標(biāo)及相應(yīng)的參數(shù)，。 和 表5 利用16個特征點求得五個靶標(biāo)的中心點的像在相機坐標(biāo)系下坐標(biāo)ABCDEX-49.7132-23.347533.819418.8173-59.7830Y51.128949.153944.8716-31.5798-31.1754Z=1577/3.78 單位：mm計算結(jié)果與利用8個特征點求得的結(jié)果比較更精確，但差別并不大，說明我們的算法的穩(wěn)定性。(2) 利用成像原理檢驗引理：一條空間直線的投影到任意一個不與之垂直的平面上像仍是一條直線。已知,三個標(biāo)靶的中心在一條直線上，從而其像的三個中心坐標(biāo)也應(yīng)該在一條直線上。若利用我們的模型求得的結(jié)果，的中心在一條直線或近似一條直線，說明模型較穩(wěn)定。檢驗方法：取中心的三個坐標(biāo)，擬合直線， 很小，說明算法穩(wěn)定。算例：利用16個特征點，根據(jù)模型(13)求得的,E中心的像坐標(biāo)如表5所示，即有： (-49.7132, 51.1289)， (-23.3475, 49.1539)， (33.8194, 44.8716)的擬合方程為：，殘差。從圖9中可以看出三點幾乎構(gòu)成一條直線。同理對特征點數(shù)為8的情形，擬合直線如圖10所示。對應(yīng)的的擬合方程為：，殘差。圖9對16個特征點的情形三點擬合的直線圖10對8個特征點的情形三點擬合的直線4問題（4）的模型(1) 二目立體視覺模型二目立體測試系統(tǒng)利用兩臺不同角度和位置安放的相機來模擬人眼的功能，獲取像點與相機的三維空間中物體的幾何信息。未來利用二目立體測量系統(tǒng)獲得三維幾何信息，需要是兩個相機在一個統(tǒng)一的坐標(biāo)系下工作，因此不僅要得到單個相機的相機模型，還需要精確確定兩個相機的相互位置關(guān)系，即建立精確的立體視覺模型。 立體視覺模型以單相機模型為基礎(chǔ)來表示兩個相機的位置幾何關(guān)系。如果左相機的外部參數(shù)為，右相機的外部參數(shù)為，則兩個相機之間的幾何關(guān)系表示如下（以左相機為基準(zhǔn)）：(2) 雙目視覺特征點三維重建算法利用基于視差原理的雙目視覺測量原理。對一般情況進行考慮，如圖11所示，假定左相機位于世界坐標(biāo)系的原點處且無旋轉(zhuǎn)，圖像坐標(biāo)系為，焦距為，右相機坐標(biāo)系為，圖像坐標(biāo)系為，焦距為，根據(jù)相機透視變換模型有與之間可通過空間轉(zhuǎn)換矩陣表示為：，P圖11與之間的旋轉(zhuǎn)矩陣和原點之間的平移變量矢量。對于坐標(biāo)系中的空間點，兩相機像平面點的對應(yīng)關(guān)系為：由此，空間點三維坐標(biāo)可表示為： 當(dāng)已知兩相機焦距、及空間低昂在左右相機中的圖像坐標(biāo)時，通過標(biāo)定出的旋轉(zhuǎn)矩陣和平移矢量就可以得到被測物體點的三維空間坐標(biāo)。五、模型的評價及推廣優(yōu)點：算法簡單，適合用計算機編寫程序計算，誤差積累少，精度較高，穩(wěn)定性好，算法可操作性強。缺點：對算法精度的分析只是利用數(shù)值分析得直觀方法，缺少理論依據(jù)，有待進一步分析。鑒于以上特點，本文建立的模型對解決類似問題的計算有一定的推廣價值。參考文獻1 張廣軍，機器視覺，北京：科學(xué)出版社，2005；2 譚曉軍、余志、李軍，一種改進的立體攝像機標(biāo)定方法，測繪學(xué)報，第35卷第2期，2006.5；3 /zhuanti/jpkc/jxzt/huafajihe/6-2.htm；4 蘇金明，王永利，MATLAB圖形圖像，北京：電子工業(yè)出版社，2005.11；5 黃宣國，空間解析幾何與微分幾何，上海：復(fù)旦大學(xué)出版社，2003.96 胡培成等，一種改進的基于圓環(huán)點的攝像機自標(biāo)定方法，光電工程，第34卷12期；7 鄒鳳嬌，蘇顯渝，李美菊，基于共面點的攝像機線性標(biāo)定法，光電工程，第32卷第4期，2005.48 何文章等，數(shù)學(xué)建模與實驗，哈爾濱：哈爾濱工程大學(xué)出版社，2002附錄主程序，文件名：QiuYuanXin_Over.mclc;clear;tic %啟動計時器%=I=imread(Xiang2.bmp);J=rgb2gray(I);K=im2bw(J,100/255);K=K;I1,num1=bwlabel(K);clear I J%*for num=1:num1 %依次對塊進行處理%=r,c = find(I1=num); %尋找塊內(nèi)點的坐標(biāo)n1,m1=size(r);n,m=size(K);M=zeros(n,m);for i=1:n1 M(r(i),c(i)=1; %構(gòu)建單個橢圓的圖像end%=M1=edge(M,sobel); %尋找邊界線%imshow(M1)y,x=find(M1=1);%=y=(abs(y-n)-n/2)/3.78; %坐標(biāo)變換x=(x-m/2)/3.78;%=xishu,bi,r,rint,s=NiHeTuoYuan(x,y); %橢圓擬合s_tuoyuan(num,:)=s;%figure,rcoplot(r,rint)%= %存儲擬合方程的系數(shù)%系數(shù)順序：ax2+bxy+cy2+dx+ey+f(=1)=0XS(num,:)=xishu;%=xc,yc=TuoYuanZhonxin(xishu); %求橢圓中心點%=if num=1 %存儲擬合橢圓的中心點坐標(biāo) TZhongXinx=xc; TZhongXiny=yc;else len1=length(TZhongXinx); TZhongXinx(len1+1)=xc; TZhongXiny(len1+1)=yc;end%=rp=abs(max(x)-min(x); %繪制擬合橢圓x1=(min(x)-rp/3):0.5:(max(x)+rp/3);%確定橫軸的范圍和步長xn,ynu,ynd=NHTYBian(xishu,x1);%=if num=1 %存儲各橢圓邊的擬合點坐標(biāo) NiheBianx=xn;NiheBianyu=ynu;NiheBianyd=ynd;else len1=length(NiheBianx); len2=length(xn); for g=(len1+1):(len2+len1) NiheBianx(g)=xn(g-len1); NiheBianyu(g)=ynu(g-len1); NiheBianyd(g)=ynd(g-len1); endend%=figure,plot(x,y,+,xc,yc,+,xn,ynu,o,xn,ynd,o);axis(equal)%=clear rp xc yc r c n1 m2 n m M M1 xn ynu ynd len1 len2endclear I J K I1 rint s %清理內(nèi)存%=disp(擬合橢圓方程的系數(shù))XSdisp(擬合橢圓的中心點坐標(biāo))TZhongXin=TZhongXinx,TZhongXiny%= %求橢圓A,B,C幾何中心的擬合直線 xcc=TZhongXinx;ycc=TZhongXiny;x=xcc(2) xcc(3) xcc(5); y=ycc(2) ycc(3) ycc(5);n=length(y);X=x,ones(n,1);Y=y;nhxs,bi,r,rint,sz=regress(Y,X);disp(橢圓A,B,C幾何中心的擬合直線系數(shù))nhxssz%figure,rcoplot(r,rint)%y1=nhxs(1).*x+nhxs(2);%figure,plot(x,y,+,x,y1)clear xcc ycc x y n X Y bi r rin%=%*%=mn1=QeiXian(XS(1,:),XS(2,:);mn2=QeiXian(XS(2,:),XS(5,:);mn3=QeiXian(XS(5,:),XS(4,:);mn4=QeiXian(XS(1,:),XS(4,:);WaiGongQeiXian=mn1(1,:);mn2(3,:);mn3(2,:);mn4(2,:)X1,Y1=Qeidian(XS(1,:),mn4(2,:);mn1(1,:);X2,Y2=Qeidian(XS(2,:),mn1(1,:);mn2(3,:);X5,Y5=Qeidian(XS(5,:),mn2(3,:);mn3(2,:);X4,Y4=Qeidian(XS(4,:),mn3(2,:);mn4(2,:);Q1=X1(:,1),Y1(:,1);Q2=X2(:,1),Y2(:,1);Q3=X5(:,1),Y5(:,1);Q4=X4(:,1),Y4(:,1);WaiQeiDian=Q1;Q2;Q3;Q4;WaiQeiDian=real(WaiQeiDian)%=NeiGongQeiXian=mn1(2,:);mn2(2,:);mn3(1,:);mn4(3,:)X1,Y1=Qeidian(XS(1,:),mn1(2,:);mn4(3,:);X2,Y2=Qeidian(XS(2,:),mn2(2,:);mn1(2,:);X5,Y5=Qeidian(XS(5,:),mn3(1,:);mn2(2,:);X4,Y4=Qeidian(XS(4,:),mn4(3,:);mn3(1,:);Q1=X1(:,1),Y1(:,1);Q2=X2(:,1),Y2(:,1);Q3=X5(:,1),Y5(:,1);Q4=X4(:,1),Y4(:,1);NeiQeiDian=Q1;Q2;Q3;Q4;NeiQeiDian=real(NeiQeiDian)clear Q1 Q2 Q3 Q4 X1 X2 X3 X4 X5 Y1 Y2 Y3 Y4 Y5%=X=WaiQeiDian(:,1);Y=WaiQeiDian(:,2); %求直線作圖邊界xmin=min(X); xmax=max(X);ymin=min(Y);ymax=max(Y);lix=xmin:0.1:xmax;liy=ymin:0.1:ymax;clear X Y%=%外切線。XJ=WaiGongQeiXian;lw1s=line(XJ(1,1),XJ(1,2),liy,2);lw2h=line(XJ(2,1),XJ(2,2),lix,1);lw3s=line(XJ(3,1),XJ(3,2),liy,2);lw4h=line(XJ(4,1),XJ(4,2),lix,1);%=%內(nèi)切線。XJ=NeiGongQeiXian;ln1s=line(XJ(1,1),XJ(1,2),liy,2);ln2h=line(XJ(2,1),XJ(2,2),lix,1);ln3s=line(XJ(3,1),XJ(3,2),liy,2);ln4h=line(XJ(4,1),XJ(4,2),lix,1);clear XJ%=figure,plot(NiheBianx,NiheBianyu,+,NiheBianx,. NiheBianyd,+,lw1s,liy,lw3s,liy,ln1s,liy,. ln3s,liy,lix,lw2h,lix,lw4h,lix,ln2h,lix,ln4h,. NeiQeiDian(:,1),NeiQeiDian(:,2),o,. WaiQeiDian(:,1),WaiQeiDian(:,2),o)axis(equal);clear lw1s lw2h lw3s lw4h ln1s ln2h ln3s ln4h%=%*%=disp(基于8個特征點的相機內(nèi)外部參數(shù)與圓中心像點求解)%順序為 %求相機內(nèi)外參數(shù)% 2(A) 3(B) 5(C)% 1(E) 4(D)Xf1=WaiQeiDian(:,1);Yf1=WaiQeiDian(:,2);%Word為世界坐標(biāo)系下的坐標(biāo)，Z坐標(biāo)取0Word1=dlmread(Yuanxindata8.m);xw=Word1(:,1);yw=Word1(:,2);fb=1577/3.78R1 T1 f1 k11 k1=TwoStage(Xf1,Yf1,xw,yw,fb)fc1=abs(fb-f1)*100/fbclear xw yw%=W1=-50 -50 0;-50 50 0;-20 50 0;50 -50 0;50 50 0;image1=WordToImage(R1,T1,fb,W1)%=ximage=image1(:,1);yimage=image1(:,2);x=ximage(2) ximage(3) ximage(5);y=yimage(2) yimage(3) yimage(5);n=length(y);X=x,ones(n,1);Y=y;pl1,bi,r,rint,s1=regress(Y,X);pl1,s1%figure,rcoplot(r,rint)y1=pl1(1).*x+pl1(2)figure,plot(x,y,+,x,y1,LineWidth,2),axis(equal)%=disp(基于16個特征點的相機內(nèi)外部參數(shù)與圓中心像點求解)Xf2=WaiQeiDian(:,1);NeiQeiDian(:,1);Yf2=WaiQeiDian(:,2);NeiQeiDian(:,2);%Word為世界坐標(biāo)系下的坐標(biāo)，Z坐標(biāo)取0Word2=dlmread(Yuanxindata16.m);xw=Word2(:,1);yw=Word2(:,2);R2 T2 f2 k21 k2=TwoStage(Xf2,Yf2,xw,yw,fb)fc2=abs(fb-f2)*100/fb%=image2=WordToImage(R2,T2,fb,W1)%=ximage=image2(:,1);yimage=image2(:,2);x=ximage(2) ximage(3) ximage(5);y=yimage(2) yimage(3) yimage(5);n=length(y);X=x,ones(n,1);Y=y;pl2,bi,r,rint,s2=regress(Y,X);pl2,s2%figure,rcoplot(r,rint)y1=pl2(1).*x+pl2(2)figure,plot(x,y,+,x,y1,LineWidth,2),axis(equal)clear bi r rint X Y x y y1 ximage yimage i lix liy m1 n . num num1 x1 xh xH yH xw yw%=%*toc %程序運行完畢，關(guān)閉計時器，并顯示程序運行時間函數(shù)程序，文件名：NiHeTuoYuan.mfunction xishu,bi,r,rint,s=NiHeTuoYuan(X,Y)%橢圓擬合n=length(Y); Xnew=X.2 X.*Y Y.2 X Y;Ynew=-ones(n,1);B,bi,r,rint,s=regress(Ynew,Xnew);%系數(shù)順序：ax2+bxy+cy2+dx+ey+1=0xishu=B,1;函數(shù)程序，文件名：TuoYuanZhonxin.mfunction X,Y=TuoYuanZhonxin(xishu)%求橢圓中心點坐標(biāo)a=xishu(1);b=xishu(2);c=xishu(3);d=xishu(4);e=xishu(5);temp=b.*b-4*a.*c;X=(2*c*d-b*e)/temp; %橢圓中心點橫坐標(biāo)Y=(2*a*e-b*d)/temp; %橢圓中心點縱坐標(biāo)函數(shù)程序，文件名：NHTYBian.mfunction xn,ynu,ynd=NHTYBian(xishu,x)a=xishu(1);b=xishu(2);c=xishu(3);d=xishu(4);e=xishu(5);longx=length(x);j=1;for i=1:longx pd=(b*x(i)+e)2-4*c*(a*x(i)2+d*x(i)+1); if pd=0 ynu(j)=(-(b*x(i)+e)+sqrt(pd)/(2*c); ynd(j)=(-(b*x(i)+e)-sqrt(pd)/(2*c); xn(j)=x(i);j=j+1; endend函數(shù)程序，文件名：QeiXian r.mfunction mn=QeiXian(k,l)s1=k(5)2-4*k(3)*k(6);s2=4*k(3)*k(4)-2*k(2)*k(5);s3=k(2)2-4*k(1)*k(3);s4=2*k(4)*k(5)-4*k(2)*k(6);s5=2*k(2)*k(4)-4*k(1)*k(5);s6=k(4)2-4*k(1)*k(6);q1=l(5)2-4*l(3)*l(6);q2=4*l(3)*l(4)-2*l(2)*l(5);q3=l(2)2-4*l(1)*l(3);q4=2*l(4)*l(5)-4*l(2)*l(6);q5=2*l(2)*l(4)-4*l(1)*l(5);q6=l(4)2-4*l(1)*l(6);syms m0 n0fun1=s1*m02+s2*m0*n0+s3*n02+s4*m0+s5*n0+s6;fun2=q1*