優(yōu)秀1數(shù)碼相機(jī)定位算法研究

高教社杯大學(xué)生數(shù)學(xué)建模競賽區(qū)評閱編號（由賽區(qū)評閱前進(jìn)行編號賽區(qū)評閱記錄（可供賽區(qū)評閱時使用中w統(tǒng)一編號（由賽區(qū)送交前編號）：評閱編號（由 評閱前進(jìn)行編號數(shù)學(xué)數(shù)學(xué)ww針對模型的求解，本文使用模擬退火算法計算出了像平面上四條公切線交點的坐標(biāo)，并使用基于最小二乘法的優(yōu)化工具箱的工具求解出靶標(biāo)的位置，進(jìn)一步求出了圓心在像平面上的坐標(biāo)，五個坐標(biāo)分別為：A0(-190.26,-196.77),B0(-88.88,-189.15),C0(129.74,-172.72),D0(72.85,119.30),E0(-229.13,119.21)在1像素以內(nèi)，說明以上模型的精度很高。偏移距離也只有0.2682，不到一個像素，說明模型具有很好的穩(wěn)定性。一、問題重數(shù)學(xué)標(biāo)定的一種做法是：在一塊平板上畫若干個點，同時用這兩部相機(jī)照相，分別得實際的做法是在物平面上畫若干個圓（稱為靶標(biāo)數(shù)學(xué)uxux像素平面坐標(biāo)系縱軸xyMRtKAB（粗體向量AB，粗體表三、問題分?jǐn)?shù)學(xué)四、模型數(shù)學(xué)1w2w34五、模型準(zhǔn)(一)靶圖像矩陣表首先將題目中的保存出來，得到的圖像可以很方便的放到里面進(jìn)行處將文件讀入，使用imread()函對于用以上方式得到的矩陣，有兩個值：0、150代表像素為白色的點，15(二)圖像輪廓的提FBGCADHEI1素，只要其四鄰域有一個值為0，則認(rèn)為這個像素為邊界上的像素。個角上的像素）值都為0，則此點為孤立點，予以剔除。數(shù)學(xué)數(shù)學(xué)2(三坐標(biāo)定O Op 23456…m數(shù)學(xué)mvuv4數(shù)圖中Op-uv坐標(biāo)系為像素平面坐標(biāo)系，Or-xryr為物理平面坐標(biāo)系。兩者的轉(zhuǎn)換關(guān)數(shù)為 d

中 wvwv vy 1 0r1 1 dxdy為像素點xy軸方向上的物理尺寸(在這里使用像素為單位)，u0和v0為點Or在像素坐標(biāo)系中的坐標(biāo)，s為圖像傾斜度，對于數(shù)碼相機(jī)，可以認(rèn)為圖像傾斜度為一個極小量，一般為0。[1]再次，我們選擇dx和dy分別為1，s=0，由于行、列像素數(shù)均為偶數(shù)，因此選擇u0=512.5，v0=384.5，此時，要計算u、v時，需要代入像素的中心的坐標(biāo)來計算。相機(jī)坐標(biāo)即整個物理世界的坐標(biāo)，現(xiàn)根據(jù)假設(shè)1、假設(shè)2做如下坐標(biāo)：以相機(jī)光學(xué)中心為原O，O與像的中Or連成的直線z軸，以O(shè)r平行于像平Or-xryr且平行Orxr的直線x軸，以O(shè)r平行于Or-xryr且平行Oryr的直y軸，建立zOxy數(shù)學(xué)5w(四)分別提取五個圓輪廓的像素坐w11。因此我們可以使用如下方i=i+1,j=j+1，跳至第ii步。若i=1024j=768，終將(i,j)記錄至Bk，并判斷(i-1,j)、(i+1,j)、(i,j+1)、(i,j-1)、(i-1,j-1)、(i-1,j+1)、(i+1,j-1)、(i+1,j+1)是否等于1，并令(i,j)為以上檢驗出等于1的數(shù)值，跳至下一步第v步通過以上步驟即可得到五個圓的輪廓的坐標(biāo)（源程序見附錄三六、模型建(一)針孔相機(jī)模型的建對于三維坐標(biāo)(X,Y,Z)，其齊次坐標(biāo)表示為(X1,X2,X3,X4)，其中X=X1/X4，Y=X2/X4，針孔相機(jī)模數(shù)

p 0R0MKTGK 0

R為旋轉(zhuǎn)矩陣，t fu K

fv/

v0

1fu,fv：相機(jī)焦距，分別以像素為單位沿水平方向和垂直方向測量。根據(jù)假設(shè)3、4，這里取fu=fv=f=1577。u0,v0：相機(jī)主點，即Or點，根據(jù)假設(shè)1，這里均為0 00K 0 因此，相機(jī)坐標(biāo)(x,y,z)映射到像物理平面坐標(biāo)(xr,yr)就可以寫成如下形式x 1 R

x2 x 0 3 xxx1,yxx x (二)確定靶標(biāo)平面相對相機(jī)位置數(shù)通過以上分析，物理世界的點可以通過次線性變換轉(zhuǎn)換成像物理平面坐標(biāo)上的點，然而這個過程并不是可逆的。三上的點不能與二上的點建立起一一對應(yīng)的關(guān)系，所以單純從像平面上的點映射回三上的點是不可能的。數(shù)盡管如此，如果我們控制一些參數(shù)（比如靶標(biāo)大小從像 在靶像上找到四個定點圓的外切線，如下圖所示，確定四個交點分別為A’、B’、C’、 6 - 01xy - 0 2 -1得到的四個點A’’、B’’、C’’、D’’既為靶標(biāo)在z=-1CO數(shù)7z=-1數(shù)四個點，它們是w由此我們可以設(shè)OA’’’=λ1OA’’,OB’’’=λ2OB’’,OC’’’=λ3OC’’,OD’’’=λ4OD’’，w而得到A’’’B’’’、A’’’C’’、B’’’D’’、C’’’D’’’的表達(dá)式。這樣我們就可以建立以下方程求解 其中l(wèi)為靶標(biāo)上線段A’’’B’’’的長個偏差最小的解，所以這8個方程都不能舍棄。利用以上方程組求出4個參數(shù)之后，就可以確定靶標(biāo)平面上四個點A’’’、B’’’、C’’’D’’’的坐標(biāo)，靶標(biāo)平面在相機(jī)坐標(biāo)中的位置也隨(三)尋找圓心的方x 1 Rx2 x 0 3 由于已經(jīng)在以上的模型中求出了靶標(biāo)所在的相對位置，因此其中的R、t分別變?yōu)閱挝痪仃嚭土阆蛄?，f已知，利用這些條件就可以求出像物理平面上圓心的坐標(biāo)。假設(shè)靶標(biāo)上的圓心坐標(biāo)為Oo(xo,yo,zo)，則圓心的坐標(biāo)計算公式為：x1 0xoyx 0o 0o數(shù)2 z數(shù)x 3 w

11 XY 七、模型的求(一)像物理坐標(biāo)系中切線的求數(shù)數(shù)w8S2

abcd 數(shù)數(shù)w9則繼續(xù)產(chǎn)生可行域，直到加入擾動后在可行域上，加入擾動后的解為r1計算ΔE=E1-一開始給定了一個最大最小范圍Lmax,Lmin，在每次迭代時，L都做如下計算：T1T數(shù)LLmaxLminT1T

T1T0w我們選T0=100,T1=0.01,Lmax=10,Lmin=1,k=10,a=0.9進(jìn)行計算（源程序見附錄四，w10求解切線結(jié)果圖（放大后(二)圓心坐標(biāo)的計具體求法我們使用里面的fsolve()函數(shù)。fsolve函數(shù)是軟件優(yōu)化工具箱零點。如果一個系統(tǒng)能表示成一個非線性方程組，便可用fsolve函數(shù)求零解。[4]fsolve函數(shù)時，初始解的選取非常重要。若初始解不合理，則fsolve函數(shù)很在像中，四邊形的長大約為400像素，即大約有105mm。而靶標(biāo)的真實長度有數(shù)數(shù)同時計算所得的四條邊平均長度誤差為3.19像素，相對誤差0.844%，誤差比較小，w因此靶標(biāo)上四條外公切線交 A’’’、B’’’、C’’’、D’’’的坐標(biāo)分別為w1657.63),(313.64,-177.82,-1815.08),(-200.12,240.73,-(-126.48,-191.92,-1984.68)如果以A’’’B’’’和A’’’C’’’作為如果以A’’’B’’’和A’’’C’’’作為平面的基向量，則五個圓的圓心分別可以表A'''ww

A'''C'''A'''C'''

A'''B''A'''B''A'''C112A'''C'''

12A''' A'''D112A'''C'''112A'''B'A'''E

A'''C'''112真空相機(jī)模型的轉(zhuǎn)換矩陣就可以很容易的找到像物理平面上對應(yīng)的圓心的坐標(biāo)。（像物理坐標(biāo)系）A0(-190.26,-196.77),B0(-88.88,-189.15),C0(129.74,-D0(72.85,119.30),E0(-（像素坐標(biāo)系）A0(322,188),B0(424,195),C0(642,212),D0(585,504),E0數(shù)12數(shù)w八、模型的精度分析與檢(一)誤差產(chǎn)生的原因分被拍攝的物體一般是連續(xù)的，然而物體拍攝到數(shù)碼相機(jī)后，原來連續(xù)的點被離散化，Photoshop中用“閾值”工具將圖像處理為黑白顏色。然而這往往是不精確的。比如我們選擇閾值為205，當(dāng)選擇206時，又有一些像素數(shù)下，這些假設(shè)總是近似符合（否則模型就會出現(xiàn)不符合實際的錯誤，而不是誤差），這數(shù)(二)精度檢驗的基本思ww(三)檢驗算法及求OA13數(shù)現(xiàn)假設(shè)像平面上任意一點A坐標(biāo)為(a,b,c)，則AO=(-a,-b,-c)，若假設(shè)OA’=λAO，數(shù)現(xiàn)在我們?nèi)“袠?biāo)平面上的點A’’’(m,n,p)，同 A’’’B’’’=(ux,uy,uz),A’’’C’’’=(vx,vy,vz)則OA’還可以表示為OA’=OA’’’+λ1A’’’B’’’+λ2A’’’C’’’=(m,m,p)+λ1(ux,uy,uz)+λ(vx,vy,vz)，用矩陣可以表示為 cz 2 m z

c u

m p

z假設(shè)映射到靶標(biāo)平面上的點位A’，預(yù)期對應(yīng)的圓心為O’，則下一步計算兩個點之間的距離|O’A’|。對于原邊緣的整個序列上的點，計算所有的距離，并計算出所有這點距離相應(yīng)圓心的距離的均值μ與標(biāo)準(zhǔn)差σμσ標(biāo)準(zhǔn)差系數(shù)ABCDE從以上表可以看出，以|O’A’|2像九、模型的穩(wěn)數(shù)數(shù)ww值只會比減少錢更小此不妨令模擬退火算法的初始點為以上我們模型確定的4個 p為輪廓上點的保留概率，則輪廓上的點被刪除的概率1-p，我們3個不同的p值下利用以上模型計算圓心的坐p值做了10次模（模擬程序見附錄七），模擬p1151離即10在一個p下10次模擬每個點偏移距離的綜合；每個p下做10次模擬，每次模擬有5個坐標(biāo)，因此每個p下有50個點，平均偏移距離既為總偏移距離除以50。表中從以上數(shù)據(jù)可以看出，在輪廓發(fā)生改變后，即使在輪廓損失了近30%的信息量時，圓十、兩部固定相機(jī)相對位置的A為原點，ACx軸正方向，AEy軸正方向，ACAE的叉積的方向為z軸正方向建立直角坐標(biāo)系。建立了世界坐標(biāo)系后，就可以確定兩個相機(jī)的光學(xué)中心O1、O2以及像平面中心Or1、Or2的坐標(biāo)。利用這四個點的坐標(biāo)，就可以確定兩臺相機(jī)的相對位置。數(shù)將給出O1的世界坐標(biāo)系的求法，其他三個點方法類似。數(shù)如下圖所示，回到針孔相機(jī)模型的世界坐標(biāo)系中，O1仍為坐標(biāo)原點，A,B,C,D,E各DD OCA14O1度為h，根據(jù)體積守恒，可以根據(jù)下式計算h：根據(jù)O1O1’的長度h以及O1O1’與AE,AC的垂直關(guān)系，求O1’的坐標(biāo)。得到坐標(biāo)后，計算a和b，使得：AO1’=aAC+bAE。A為原點的世界坐標(biāo)系中，O1(a,b,h)利用上述方法，得到四個點O1、O2、Or1、Or2在世界坐標(biāo)系中的坐標(biāo)，就可以得到兩十一 模型擴(kuò)數(shù)rn

12nnnp(X,Y,Z12nnn p(X, , 式中n和cn分別是標(biāo)準(zhǔn)化的像素坐標(biāo)，Xn，Yn和Zn分別是標(biāo)準(zhǔn)化的地面坐標(biāo)。所謂標(biāo)準(zhǔn)化是指2個像點坐標(biāo)和3個地面點坐標(biāo)分別平移和縮放使它們的值落在paijk*Xi*Yj* i0j0pp pp Xn，Yn的反解形母多項式的不同組合而變化，共有9種情況，如下表：3p2333p2222p2716184 p(X,Y,Z)

n2nnn2nnn

p(X,Y,ZXnXn,n)

XY 創(chuàng)建像格創(chuàng)建像格15十二 模型評十三 參考文陸宏偉，基于編隊遙感圖像的對地定位算法研究，國防科技大學(xué)，2004年數(shù) 日數(shù)，基于多像滅點進(jìn)行相機(jī)標(biāo)定的方法研究，大學(xué)，2004年11月1日w513w十四 附附錄一原圖像數(shù)據(jù)在中的處forifa(i,j)==15

cleara;cleari;clearfori=1:length(dat(:,1))ifdat(i,j)==1ifdat(i-1,j)==1&dat(i+1,j)==1&dat(i,j-1)==1&dat(i,j+1)==1

數(shù)%檢測孤立數(shù)fori=1:length(dat(wifdat(i,j)==1wifdat(i-1,j)==0&dat(i+1,j)==0&dat(i,j-1)==0&dat(i,j+1)==0&dat(i-1,j+1)==0&dat(i-1,j-1)==0&dat(i+1,j+1)==0&dat(i+1,j-1)==0

cleari;clearwhileforforifififi==768&j==1024

whileif數(shù)數(shù)ww

forifa(i,j)==1ifa(i-1,j)==0&a(i+1,j)==0&a(i,j-1)==0&a(i,j+1)==0&a(i-1,j+1)==0&a(i-1,j-1)==0&a(i+1,j+1)==0&a(i+1,j-1)==0

cleara;clearh;cleari;clearj;clearhh;clear數(shù)functionf=isMatch(r,zb_y)數(shù)wfori=1:length(b)wif~all(a-%目標(biāo)函functionf=mbFunction(r)%生成隨機(jī)擾%n為向量維%mfunctionfunctionf=zb_qx(r)if

forifa(i,1)>=min([r(3),r(1)])&a(i,1)<=max([r(3),r(1)])

endwwif

forifa(i,2)>=min([r(2),r(6)])&a(i,2)<=max([r(2),r(6)])if

forifa(i,2)>=min([r(8),r(4)])&a(i,2)<=max([r(8),r(4)])if

數(shù)數(shù)wforifa(i,1)>=min([r(7),r(5)])&a(i,1)<=max([r(7),r(5)])dE能量差whileT>=TnwhilewhileifisMatch(r1,zb_y)ifdE<=0

數(shù)ifp>=rand(1)數(shù)wwif

clearE;clearT;clearTn;cleardE;clearh;clearhh;cleark;clearl;clearp;clearr1;cd'..'%cd zb0=[best(2)best(1);best(4)best(3);best(6)best(5);best(8)best(7)];fori