數(shù)碼相機數(shù)學(xué)建模

1、2008高教社杯全國大學(xué)生數(shù)學(xué)建模競賽承諾書我們仔細閱讀了中國大學(xué)生數(shù)學(xué)建模競賽的競賽規(guī)則.我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的，如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻的表述方式在正文引用處和參考文獻中明確列出。我們鄭重承諾，嚴格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴肅處理。我們參賽選擇的題號是（從A/B/C/D中選擇一項填寫）：A我們的參賽報名號為（如果賽區(qū)設(shè)置報名號的話

2、）：1324615所屬學(xué)校（請?zhí)顚懲暾娜褐袊刭|(zhì)大學(xué)（武漢）參賽隊員（打印并簽名）：1.H2.3. 代永力指導(dǎo)教師或指導(dǎo)教師組負責(zé)人（打印并簽名）：付麗華日期：2008年9月22日賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：2008高教社杯全國大學(xué)生數(shù)學(xué)建模競賽編號專用頁賽區(qū)評閱編號（由賽區(qū)組委會評閱前進行編號）：賽區(qū)評閱記錄（可供賽區(qū)評閱時使用）：評閱人評分備注全國統(tǒng)一編號（由賽區(qū)組委會送交全國前編號）：全國評閱編號（由全國組委會評閱前進行編號）：數(shù)碼相機定位摘要相機定位是計算機視覺領(lǐng)域里從二維圖像獲取三維信息的基本要求，是完成許多視覺工作必不可少的步驟。題目利用由五個圓組成的靶標模

3、板（圖2）對相機進行標定。對于給出的相機拍攝的靶標的像（圖3）,本文對其進行了預(yù)處理：首先應(yīng)用MATLAB軟件對圖3進行去噪，然后利用自行設(shè)計的算法（掃描法）提取去噪后圖形的邊緣（圖8）。事實上，標定就是找到空間坐標系和像平面坐標系之間的相互關(guān)系。而這是由相機成像的幾何模型決定的，具成像參數(shù)包括相機內(nèi)部幾何和光學(xué)特性（內(nèi)部參數(shù)），以及相機世界坐標系的三維位置和方向（外部參數(shù)）。求解相機內(nèi)外參數(shù)是解決問題的關(guān)鍵。本文考慮了理想的情況，即直線的投影仍是直線，以及交點的投影仍然是像的交點。利用圖2中A、B、C、D、E兩兩圓的外公切線的切點與圖3中像的外公切線的切點相對應(yīng)，通過最小二乘法求出該相機的內(nèi)

4、、外參數(shù)。然后，利用相機成像原理，將圖2中圓心的坐標Ac、Bc、Cc、Dc、Ec和相機內(nèi)外參數(shù)代入，即可得到靶標上圓的圓心在該相機像平面的像坐標Ac、Bc、Cc、Dc、Ec（見表2）。為了驗證該模型，本文利用留一法，用四個圓（共有三組：A、B、D、E;B、C、D、E;A、C、D、E）的切點與對應(yīng)像點求相機內(nèi)外參數(shù)，剩下的一個圓用于做檢驗。以第一組為例，具體過程是：根據(jù)A、B、D、E算出內(nèi)外參數(shù)，然后求出C圓圓心的像的坐標CJ,并與Ac、Bc所確定的直線的距離來檢驗方法的精度和穩(wěn)定性。結(jié)果顯示（見本文13-14頁），本文所用方法精度較高，穩(wěn)定性較好。雙目標定時，將雙目系統(tǒng)的坐標系建立在左相機上，

5、把右相機的坐標系相對于其的旋轉(zhuǎn)和平移參數(shù)利用最小二乘法求出，從而可以得到兩部相機的相對位置。關(guān)鍵詞：標定；內(nèi)、外參數(shù)；深度優(yōu)先搜索；最小二乘法一、問題的提出1 .背景目前，交通管理中使用人工監(jiān)控的方式已經(jīng)不適合日益壯大的交通系統(tǒng)，而有效的電子監(jiān)管（例如數(shù)碼相機定位）能更有效的管理交通系統(tǒng)運行，于是電子監(jiān)管逐漸被人們所認識。隨著數(shù)碼相機的普及，數(shù)碼相機定位在交通監(jiān)管（電子警察）等方面應(yīng)用越來越廣泛。所謂數(shù)碼相機定位是指用數(shù)碼相機攝制物體的相片確定物體表面某些特征點的位置。當前，最常用的定位方法是雙目定位，即用兩部相機來進行定位。對物體上一個特征點，用兩部固定于不同位置的相機攝得物體的像，分別獲得

6、該點在兩部相機像平面上的坐標。如果知道兩部相機精確的相對位置，就可使用幾何的方法得到該特征點在固定一部相機的坐標系中的坐標，即確定了特征點的位置。由此可見，對于雙目定位，精確地確定兩部相機的相對位置就是關(guān)鍵，這一過程稱為系統(tǒng)標定。系統(tǒng)標定方法很多，常用的有張正友平面標定方法、孟胡平面圓標定方法、平行圓標定方法1。目前有人提出這樣一種方法：在一塊平板上畫若干個點，同時用這兩部相機照相，分別得到這些點在它們像平面上的像點，利用這兩組像點的幾何關(guān)系就可以得到這兩部相機的相對位置。然而，無論在物平面或像平面上我們都無法直接得到?jīng)]有幾何尺寸的“點”。實際的做法是在物平面上畫若干個圓（稱為靶標），它們的圓

7、心就是幾何的點。而它們的像一般會變形，如圖1所示。問題的關(guān)鍵是，能否從靶標上的這些圓的像中把圓心的像精確地找到，如果能找到，標定就可實現(xiàn)。圖1像的變形示意圖2 .問題有人設(shè)計靶標如下，取1個邊長為100mm的正方形，分別以四個頂點（對應(yīng)為A、C、D、E）為圓心，12mm為半徑作圓。以AC邊上距離A點30mm處的B為圓心，12mm為半徑作圓，如圖2所示。圖2靶標示意圖用一個位置固定的數(shù)碼相機攝得其像，如圖3所示。圖3靶標的像研究的問題如下：（1）建立數(shù)學(xué)模型和算法以確定靶標上圓的圓心在該相機像平面的像坐標，這里坐標系原點取在該相機的光學(xué)中心，x-y平面平行于像平面；（2）對由圖2、圖3分別給出的

8、靶標及其像，計算靶標上圓的圓心在像平面上的像坐標，該相機的像距（即光學(xué)中心到像平面的距離）是1577個像素單位（1毫米約為3.78個像素單位），相機分辨率為1024X768;（3）設(shè)計一種方法檢驗?zāi)Ｐ?，并對方法的精度和穩(wěn)定性進行討論；（4）建立用此靶標給出兩部固定相機相對位置的數(shù)學(xué)模型和方法。二、符號定義OwXwywZw世界坐標系Oc數(shù)碼相機透視中心OcXcycZc數(shù)碼相機坐標系Ou主光軸與像平面的交點OuXuyu坐標軸分別平行于像素坐標系的坐標軸建立的物理坐標系fx0U0F數(shù)碼相機內(nèi)參矩陣，F(xiàn)0fyV0,fx，fy分別表示在X001方向和y方向上像點的物理坐標到圖像坐標的比例系數(shù)，U0,V0

9、是主光軸與像平面交點的圖像像素坐標RT數(shù)碼相機的外參矩陣，R和T分別為模板相對數(shù)碼相機的r1,1r1,2r1,3t1旋轉(zhuǎn)分量和平移分量，R3收22,3,Tt2r3,1r3,2r3,3t3Ac、Bc、Cc、Dc、Ec圓A、B、C、D、E在靶標平面上的圓心A、Bc、Cc、Dc、Ec圓心Ac、Bc、Cc、Dc、Ec在像平面上對應(yīng)的點三、問題的分析根據(jù)題意，要找到靶標圓的圓心在像平面中的對應(yīng)點，就應(yīng)該確定數(shù)碼相機的內(nèi)、外參數(shù)，然后利用所得到的內(nèi)、外參數(shù)來確定靶標平面上圓心在像平面上對應(yīng)點的坐標。只要尋找靶標上的點在像平面上的對應(yīng)點，就可以利用成像原理求出相機的內(nèi)、外參數(shù)。利用圓及其切線對投影面的特點，

10、即兩圓及其公切線在投影面上的圖形仍保持直線的投影公切于兩圓的投影，且它們在靶標平面上的切點對應(yīng)于投影面上相應(yīng)的切點，利用這一性質(zhì)就可以找到靶標平面與投影面上部分相對應(yīng)的點，然后就可以用所建立的模型來求出數(shù)碼相機的內(nèi)、外參數(shù)，從而就找到了靶標平面上點與投影面上點的對應(yīng)關(guān)系，再代入數(shù)據(jù)就可以求出靶標上圓的圓心在像平面上的像坐標。模型穩(wěn)定性和精度的驗證，主要是考查該模型對數(shù)據(jù)變化這一影響下所得結(jié)果的考量。這里，我們將利用留一法，同時用3組數(shù)據(jù)來進行驗證。問題4中兩部相機相對位置的確定，主要是兩相機像平面坐標系之間的轉(zhuǎn)化，我們將以左相機為基準，將右相機坐標系統(tǒng)一到左相機上，從而得到兩者相對位置。四、模

11、型假設(shè)1 .假設(shè)數(shù)碼相機成像的原理為針孔相機模型；2 .假設(shè)相機本身的配置沒有誤差；3 .假設(shè)數(shù)碼相機成像產(chǎn)生較小的畸變，對結(jié)果影響不大，忽略不計；4 .假設(shè)像平面中心與光學(xué)中心連線垂直于像平面。五、模型的建立與問題的解（一）模型的建立建立物點（世界坐標）和像點（像坐標）位置對應(yīng)關(guān)系方程。如圖4所示：在模板平面上，分別以兩條正交的直線為xy軸，以靶標上正方形的中心為原點，并按右手準則建立世界坐標系OwXwyw4,然后以數(shù)碼相機透視中心Oc為坐標原點，并使其x軸和y軸分別平行于像面坐標系的x軸和y軸,且方向相同，z軸為數(shù)碼相機鏡頭的主光軸，按右手準則建立數(shù)碼相機坐標系OcXcycZc0以像平面上

12、像素坐標值0,0的點為坐標原點O,以水平像元陣列方向為X方向，以垂直像元陣列方向為y方向，建立像素平面坐標系OUV。最后以主光軸與像平面的交點Ou為坐標原點，并使坐標軸分別平行于像素坐標系的坐標軸來建立物理坐標系Ou凡九，模板平面上Pw的理想成像點為R。圖4透視投影模型示意圖假設(shè)模型為理想透視模型根據(jù)針孔成像原理，由世界坐標點到理想像素點的齊次變換2,3如下：u s v1F R T ywzw 1其中，s為照相機坐標向像平面坐標轉(zhuǎn)化的比例系數(shù)，F(xiàn)為數(shù)碼相機內(nèi)參矩陣,RT為數(shù)碼相機的外參矩陣，R和T分別為模板相對數(shù)碼相機的旋轉(zhuǎn)分量和平1,11,21,32,12,22,3,Tr3,13,23,3t1

13、t2t3移分量，分別定義如下3,4fx0UoF0fyVo,R001其中，fx,fy分別表示在x方向和y方向上像點的相機坐標到物理坐標的比例系數(shù)，U0,V0表示主光軸與像平面交點的圖像坐標。由于世界坐標系的xy軸在模板平面上，因此對于模板上的點zw00這樣式（1）可重寫為:fx000fy0Uori,iV01,1,1消去s即得到以下理想透視數(shù)學(xué)模型：fxri,iXwri,2ywr1,2r2,2,2tit2t3xwyw13,1xwr3,2ywt3fy2,1xw2,2Ywt23,1Xw3,2ywt3t1U0V0（二）模型的求解1.圖像處理題目中給定一幅使用待標定數(shù)碼相機拍攝的圖片（圖為1024*768

14、大小的位圖。圖3中有5個黑色非正規(guī)圓斑，分別標號為ABC、黑色非正規(guī)圓斑周圍為暗灰色，與圓斑顏色相差很大。使用3靶標的像），規(guī)格D、E,如圖5所示。MATLAB數(shù)學(xué)軟件讀入圖片后，將圖片轉(zhuǎn)化為1024*768的0-1矩陣，圖中黑色部分標記為0,非黑色部分標記為1。這樣，將圖片轉(zhuǎn)化為一個離散化的矩陣，并用Morlet小波進行去噪處理，去噪結(jié)果如圖5所示。圖5黑斑標號圖2,邊緣提取邊緣提取就是提取圖片中特征區(qū)域的邊界，該處指提取圖5中黑斑的邊界。邊緣提取算法很多，MATLAB軟件中有邊緣提取功能，但效果并不理想。由于黑斑的特點，本文設(shè)計新的算法進行邊緣提取，稱為掃描法（算法實現(xiàn)C程序碼見附錄一）。

15、掃描法的主要原理是：分兩次掃描，第一次按行掃描，第二次按列掃描。按行掃描時，對每一行從左往右遍歷，如果第I個點到第J個點全為黑色點，且第I-1、J+1點不為黑色，則標記第I點和第J點，其余點不標記。按列掃描與按行掃描方法一樣。按行按列掃描完后，被標記的點即為邊緣點。圖6單行像素掃描前示意圖圖6表小一行像素點，黑球表小黑色點，白球表小不為黑色的點。用按行掃描掃描圖6,從左往右遍歷，可得到兩段黑色點段，分別為36、911,標記點3、6、9、11,其余點不標記。保留下來的邊緣點如圖7所示：圖7單行像素掃描后示意圖利用掃描法對圖5進行處理，得到效果圖如圖8。圖8對圖5進行邊緣提取的效果圖通過上面算法，

16、得到黑斑邊緣矩陣后，再對該矩陣使用深度優(yōu)先搜索將5個黑斑邊緣分離出來，存儲各自的邊緣點，并將5個黑斑邊緣標號（見圖8）。3.計算機模擬求切點對圖8中任意兩邊緣曲線的外公切線，其他點都在該直線的同一側(cè)，利用該特點，通過計算機模擬來求兩邊緣曲線的公切線及其切點。具體的模擬算法：對兩個邊緣曲線，以A、B為例，在曲線A上任取一點a,曲線B上任取一點b,作直線L,然后窮舉A和B上的其他點，看是否在該直線的同側(cè)，如果在同側(cè)，既可認為該直線為A、B的公切線，這兩點即為切點。例如，根據(jù)A邊緣點集和B邊緣點集，求兩邊緣曲線的兩條外公切線，并求得公切點1,3,5,6。根據(jù)上述算法，編寫部分C+程序（見附錄二），求

17、得的A和B外公切線效果圖如圖9:圖9切線效果圖根據(jù)已建立的坐標系，將此算法應(yīng)用于AB、C、D、E圓，可求得像上18個點坐標，這18個點分布如圖10所示。圖10像坐標系及切點對應(yīng)點靶標中各圓之間外公切線與圓切點標號，如圖11所示。圖11靶標坐標系及切點編號圓上切點及其像圖10和圖11上切點在相應(yīng)坐標系內(nèi)的坐標見表1。表1圓上切點及對應(yīng)像的在相應(yīng)坐標系中的坐標、坐標切）圓切點在世界坐標系中的坐標（mm）圓上切點的像在像坐標系中的坐標（mm）xyxy1-506248.14815-62.433862-6250-61.11111-53.439153-5038-52.38095-41.0052914-38

18、50-39.41799-50.264555-2062-22.48677-60.317466-20P38-26.98413-39.4179917385024.07407-46.031758506236.77249-55.291019625042.85714-43.9153410503832.01058-35.4497411-50-38-57.6719621.6931212-62-50-70.6349230.6878313-50P-62-61.1111140.21164114-38-50-50.2645531.746031538-509.5238130.952381650-3822.2222222

19、.751321762-5027.2486832.010581850-6218.7830739.682544.數(shù)碼相機的內(nèi)、外參數(shù)在實際求解的過程中最終的像坐標采用物理坐標形式（x,y）,單位為mm,且fxfxd,fyfyd,其中d為單位像素的物理尺寸（1/3.78）。這樣靶標上點與像平面上的點對應(yīng)關(guān)系如下：Xsy1fX0X0fy000011,1r2,13,1r1,2t1r2,2t23,2t3Xwyw1(5)由上式可得：Xsy1一fX1,1fy2,1r3,1一fX1,2fyr2,2r3,2一fXt1fyt：t3IXw2yw1(6)令單應(yīng)矩陣45A3,4a1a4a2a3a5a6fX1,1fy2,1

20、fX1,2fy2,2fXt1fyt2(7)a7a813,13,2t3則可把參數(shù)求解問題轉(zhuǎn)化為先求單應(yīng)矩陣A,然后從A矩陣分離出待求參數(shù)的問題。為求取A陣，變換式(6),消去s/3,4,并分離待求參數(shù)，則可得aiXwyw1000XwXywXa2x(8)000Xwyw1XwyywyMya8本文獲得了18個切點，將其帶入該方程組，共36個方程。已知點的個數(shù)遠遠超過解該線性方程組的需要，這是一個超定方程組。在此使用最小二乘法，解這個超定方程組，求出最優(yōu)的參數(shù)(程序見附錄三)。求出單應(yīng)矩陣A后，即可提取出fx，fy以及外部參數(shù)矩陣：,42222a1a2a4a5a1a2a4a51/1/f21 Xf2y22

21、a8a7122a2a5fX2fy22a82a1fX212a42一,2a7fy1,12,13,11,22,23,2t1t2t3fX 0 00 fy 0100 A1(10)(11)旋轉(zhuǎn)矩陣R的第三列可由其正交性求得：這里為向量叉積。根據(jù)表1中點的坐標求出相機的內(nèi)、外參數(shù):5.求靶標圓心像在像坐標系中的坐標將靶標圓心在世界坐標系中的坐標代入方程，使用物理坐標公式，求靶標圓心像在像坐標系中的坐標:fxrm例L(12)3,lXwr3,2ywt3.fyr2,lXwr2,2ywt2r3,lXwr3,2ywt3求出的圓心的像坐標晨bC、cC、dC、eC見表2(根據(jù)題目，若以相機光學(xué)中6.兩相機相對位置的確定假

22、設(shè)兩部相機的內(nèi)部參數(shù)相同，兩部相機的坐標系建立同已提出的模型。為 了確定兩部相機的相對位置，將左、右兩部相機中左相機固定，即以左相機坐標 系為基準，則右相機相對于左相機的坐標有如下轉(zhuǎn)換關(guān)系：XxrV R yr tzzrxryrzr(13)心為原點，則像點的z坐標均為像距417.19576mm):表2靶標圓心及其像點坐標坐標圓編號、靶標圓心坐標(mm)靶標圓心的像坐標(mm)XYXYA-5050-50.460651.7391B-2050-23.801549.6625C505033.897445.1681D-50-50-59.9735-30.9809E50-5019.3974-31.5215其中x

23、l,y,zT,xr,yr,zrT分別是左相機和右相機的像點坐標，r1,1r1,2r1,3Rr2,1r2,2r2,3為變換的旋轉(zhuǎn)分量,r3,1r3,2r3,3t1Tt2為變換的位移分量t3由于在建立坐標系時,相機的像平面上的點位于xy平面上，故z坐標為零，因此轉(zhuǎn)換關(guān)系簡化為：xlr1,1r1,2Lxrylr2,1r2,2t2yr(14)1r3,1r3,2t31通過模板上圓切點在兩相機上的像，可用最小二乘法把九個參數(shù)標定。旋轉(zhuǎn)矩陣的第三列可根據(jù)其正交性由前兩列算出。兩坐標系的轉(zhuǎn)換關(guān)系求出，即可確定它們的相對位置關(guān)系。六、模型的精度和穩(wěn)定性分析在上一節(jié)模型的求解過程中，我們利用A、B、C、D、E5個

24、圓上18個切點與像坐標求出相機內(nèi)、外參數(shù)，為了驗證該模型，我們利用留一法，只利用4個圓(A、B、D、E;B、C、D、E;A、C、D、E)來計算另一個圓心像坐標，以A、B、D、E這一組為例，代入圓A、B、D、E上切點坐標及其相對應(yīng)的投影點坐標，求出數(shù)碼相機的一組內(nèi)、外參數(shù)，利用這組參數(shù)求出圓C的圓心在投影面上對應(yīng)點的坐標C；,該點應(yīng)該在圓A和圓B的圓心在投影面上像點，即aC,bC所成直線上面，然而由于不同的參數(shù)導(dǎo)致了誤差，通過求C；到該直線的距離dc應(yīng)反映所求參數(shù)誤差的大小。同理，用B、C、D、E和A、C、D、E求出另外兩組參數(shù)后，進行類似的處理。三組數(shù)據(jù)如下：1 .用AB、D、E上切點求相關(guān)參

25、數(shù)得：Cc的坐標為(32.5312,45.2120),與AB。所成直線距離為0.062355mm2 .用RC、D、E上切點求相關(guān)參數(shù)得：Ac的坐標為(-51.5430,51.8816),與BC。所成直線距離為0.058033mm3 .用AC、D、E上切點求相關(guān)參數(shù)：B；的坐標為(-25.3578,52.7422),與AC；所成直線距離為2.949541mm分析數(shù)據(jù)可知本文提出的模型精確度比較高，前兩組誤差均不到一個像素，即基本重合，完全可以滿足要求。第三組沒有考慮B圓相應(yīng)的切點，誤差較大，可見B圓的設(shè)置對計算精度有較大影響。這可能與B圓離中心較近，產(chǎn)生畸變較小所致5,6。七、模型的優(yōu)缺點優(yōu)點：

26、本模型在提取圖像邊緣時，采用掃描法，對于噪音極小的簡單的圖案，提取完全精確，比一般的通用邊緣檢測算法效果要好；該模型采用最小二乘法求取變換參數(shù)，容易實施，且不失其精確性和穩(wěn)定性；本文采用留一法進行檢驗，檢驗結(jié)果精度高穩(wěn)定性好。缺點：模型沒有考慮相機物理因素及制作工藝缺陷而導(dǎo)致的成像畸變，所以和實際的成像結(jié)果會有一定的偏差；在使用最小二乘法計算相機內(nèi)、外參數(shù)時，由線性方程組所得的線性方程沒有經(jīng)過規(guī)格化處理，所求結(jié)果會有一定的影響；計算旋轉(zhuǎn)矩陣第三列是直接使用了其正交性，使用叉積計算，而由于計算誤差，已得參數(shù)可能并不嚴格滿足正交性，這樣計算得出的結(jié)果并不十分精確(但對以后計算圓心的像坐標沒有影響)

27、。參考文獻1李偉、呂曉旭、錢曉凡、鐘麗云?；谄矫婺０宓臄z像機標定方法比較J,激光雜志，第27卷，第2期，54-55,2006年。2毛劍飛、鄒細勇、諸靜。改進的平面模板兩步法標定攝像機J,中國圖象圖形學(xué)報，第9卷，第7期，846852,2004年。3雷彥章、趙慧潔、姜宏志。一種單雙目視覺系統(tǒng)結(jié)合的三維測量方法J,光學(xué)學(xué)報，第28卷，第7期，1338-1342,2008年。4吳謹，王俊。基于投影矩陣的攝像機參數(shù)標定J,福建電腦，第8期，1527,2008年。5孟曉橋，胡占義。一種新的基于圓環(huán)點的攝像機自標定方法J,軟件學(xué)報，第13卷，第5期，957965,2002年。6J.Heikkila.Ge

28、ometriccameracalibrationusingcircularcontrolpointsJ.IEEETransactionsonPatternAnalysisandMachineIntelligence,October,22(10),10661077,2000.附錄一：掃描法算法C代碼由圖像處理已求出圖片的0-1矩陣，設(shè)矩陣為a口口，aIJ表示位置為第I行第J列的值，0為黑色點，1為非黑色點。處理完后，圖像矩陣保存在s口口中，值為0表示為邊緣點，值為-1的點位非邊緣點。include#inclueinta1025769;ints1025769;intmain#()inti,j；fr

29、eopen(data.txt,r,stdin);freopen(point.txt,w,stdout);memset(s,-1,sizeof(s);for(i=0;i768;i+)for(j=0;j1024;j+)intt;scanf(%d,&t);a皿i=t;flagji=false;for(intii=0;ii1024;ii+)for(intjj=0;jj768;jj+)if(aiijj=0)sii皿=0;intkk=jj;while(aiikk=0&kk768)kk+;kk-;jj=kk+1;siikk=0;for(intjss=0;jss768;jss+)for(intiss=0;is

30、s1024;iss+)if(aissjss=0)sissjss=0;intkss=iss;while(akssjss=0&kss1024)kss+;kss-;iss=kss+1;skssjss=0;for(inti=0;i1024;i+)for(intj=0;j768;j+)printf(%d,sij);printf(n);return0;附錄二：求公切線C+弋碼structnodeintx,y;r2500;/r0保存A的所有邊緣點，r1保存B所有邊緣點intt2;/t0表示A邊緣點個數(shù)，t1表示B邊緣點個數(shù)boolCheck(doublea,doubleb,doublec,intt1,int

31、t2,doublew)inti;for(i=0;itt1;i+)doublew1=a*rt1i.y+b*rt1i.x+c;if(w1*w0)returnfalse;for(i=0;itt2;i+)doublew1=a*rt2i.y+b*rt2i.x+c;if(w1*w0)returnfalse;returntrue;voidGetQieDian()boolok=false;for(inti=0;it0;i+)for(intj=0;jt1;j+)if(r0i.x=r1j.x|r0i.y=r1j.y)continue;doublek=(double)(r0i.y-r1j.y)/(r0i.x-r1j.x);doublea=1;doubleb=-k;doublec=k*r0i.x-r0i.y;doublew=c;ok=Check(a,b,c,0,1,-w);/abc為直線方程參數(shù)if(ok)point1=r0i;/保存切點坐標point5=r1j;break;if(ok)break;for(inti=0;it0;i+)for(intj=0;jt1;j+)if(r0i.x=r1j.x|r0i.y=r1j.y)contin