07数码相机定位问题全国一等奖

1、数码相机模型摘要本文依据靶标的像图与原图的点线关系，建立了公切线模型；通过拟合切点连线求交点的方法,得到圆心的坐标。在模型实现上，我们首先对图像数据进行了除噪点、二值化和边缘检测的处理，可以将畸变看成误差，然后将割线在区域内单方向移动，其与图像各点的关系，最后得出公切线。我们引入了 DLT 模型来对问题一的模型进行检验。问题一：在进行有关图像处理和相机标定时，大部分模型都是采用了坐标转换的方法，即由世界坐标、相机坐标向图像坐标转换或逆转换。而在这个问题中，我们发现从原图到像图的变化中，点和线的位置关系（点是否上）不变，连接圆上两条平行切线切点的直线必经过圆心，因此我们建立了公切线模型。这个模型

2、是通过做平行于连结圆心的正方形的 8 条切线，并在靶标图像上寻找与之对应的这 8 条切线而建立。考虑到畸变对图像的影响，将得到的各个切点，线性拟合，并找出其交点，即是圆心。问题二：我们首先对图像进行了预处理，即除噪点、二值化和边缘检测，然后将 5个图形划分在 5 个不相交的矩形内，在两个矩形内各取一点，连结这两点的直线必是两圆的公割线，再让公割线沿一个方向在靶标点上移动，当不能移动时，公割线即是公切线，同时得到各个切点的坐标。拟合各个切点对应靶标原图正方形的连线，连线的交点就是圆心。最后通过编程，可得 A,B,C,D,E 对应靶标像图各点的坐标分别是：(-186.383,197.253) (-

3、91.250,189.178) (132.353,170.198)(78.000,-118.478) (-224.73,-116.581)(注：是像素)问题三：我们采用了基于 DLT 模型的一种数码相机检验方法，首先将问题二中得到的切点作为已知点，然后利用二维 DLT 变换求解出内外元素初值，再得到像中心的坐标，最后与第二问中得到的像中心坐标对比，发现相对误差都在 5%以内，这说明问题一中的模型精度很高。问题四：首先利用传统方法对单个相机的内参数进行标定，然后将靶标当作双目匹配点来用,表示出两相机与匹配点的矩阵转化关系，列出矩阵方程，从而线性地求解“本质矩阵”,最后根据“本质矩阵”的原理找出相

4、机的相对位置。1公切线畸变DLT 标定双目标定本质矩阵一、 问题重述数码相机在交通监管（电子）等方面有广泛的应用。所谓数码相机是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的方法是双目，即用两部相机来。对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。于是对双目为系统标定。，精确地确定两部相机的相对位置就是关键，这一过程称标定的一种做法是：在一块平板上画若干个点， 同时用这两部相机照相，分别得到这些点在它们像平面上的

5、像点，利用这两组像点的几何关系就可以得到这两部相机的相对位置。然而，无论在物平面或像平面上我们都无法直接得到没有几何的“点”。实际的做法是在物平面上画若干个圆（称为靶标），它们的圆心就是几何的点了。而它们的像一般会变形，所以必须从靶标上的这些圆的像中把圆心的像精确地找到，标定就可实现要求解决下列问题：（1）建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标, 这里坐标系原点取在该相机的光学中心，x-y 平面平行于像平面；对图中分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标, 该相（2）机的像距（即光学中心到像平面的距离）是 1577 个像素(1 毫米约为 3.78 个)，

6、相机分辨率为 1024×768；像素设计（3）（4）法检验的模型，并对方法的精度和稳定性进行讨论；建立用此靶标给出两部固定相机相对位置的数学模型和方法 。二、问题的分析问题立模型确定靶标的圆的圆心在该相机像平面的像坐标。首先用 PHOTOSHOP 图像处理软件将靶标的像图为数据矩阵，然后求出靶标像的公切线，再通过对畸变因素的分析，进而利用所求的切点的数据求出靶标中心点在像素坐标系中的坐标，最后利用上的坐标。几何中的坐标平移，得出靶标的圆心在相机像平面问题二：在已经给定的靶标原像图和靶标的像图的条件下，确定靶标的圆的圆心在2该相机像平面的像坐标。将上述问题一的模型在中实现，得出靶标的圆

7、心在相机像平面上的坐标。问题三：寻找另法得到圆心坐标，并与之对比得到两种方法的相对精确度。同入多个点，分析其稳定性。问题四：题目要求给出两部固定相机相对位置的数学模型和方法。首先建立单个摄像机模型，再建立双目模型，最后通过本质矩阵可以求出两部相机的相对位置。三、模型的假设.假设靶标原图在世界坐标 X-Y 平面内。假设靶标是精确制造的假设相机是固定的光照等客观条件不变四、模型的建立与求解4.1 凸透镜的成像原理如图一图一其中 u 为物距，f 为透镜焦距，V 为像距，三者满足下列关系:1 + 1 = 1uvf（1）f » v 。这里 u>>v,由（1）式可取34

8、.2 靶标的像的公切线模型的建立4.2.1 图像预处理(1)(2)(3)(4)(5)用 PHOTOSHOP 图像处理软件将 RGB 图像转化为亮度图像； 对图像进行去噪处理；将图像二值化；对图像进行边缘检测，做出检测后图像； 建立图像坐标系，如图二。图二4.2.2 公切线算法的建立首先用五个尽量小的矩形将五个圆的图像框起来,具体数据见表一。然后设点(x1, y1) 和点(x2 , y2 ) 分别为步骤三所求的其中两个图的行和列的平均值. 最后以求两个图形上方的公切线为例来说明求公切线的方法（如图三）：4表一：五个矩形边缘的值图三具体实现求公切线的算法如下：（1）建立通过点(x1, y1) 和点

9、(x2 , y2 ) 的直线 L1 方程2 - y1 (x - x )f (x12 - x1（2）利用点与直线位置关系：点(x,y)在直线上方的条件为，f(x)>0。在圆周上进行搜索，若(x,y)满足 f(x)>0，则将直线 L1 向上移动到经过(x,y)，另一个经过的点不变。一个圆周搜索完毕，则在另一个圆周上搜索。如此一直重复下去。（3）当平移至不满足（2）中的条件时，也就是平移至 o1”和 o2”时，就会得到直线 L2,此直线即为两个图形的公切线。其他各切点，切线均可按该搜索算法得到。求两个图形下方，左方，右方的公切线的方法同上。这样，我们便求出了所有的切点的坐标和公切线，如图

10、四所示。其中 H1、H2、H3、H4、H5、H6、H7、H8、H9、H10 为公切线；a, a1,a2,a3, a4, a5, a6, a7, a8, a9, b,b1, b2, b3, b4, b5, b6, b7,b8,b9 为切点，切点的坐标见表二。5上（行）下（行）平均（行）左（列）右（列）平均（列）A538620579282365323.5B532612572383464423.5C517593555604677640.5D233298265.5549618583.5E231302266.5246324285表二：靶标像切点的坐标图四4.3 畸变因素在近轴光线和近轴物的条件下，垂直于

11、主轴的物所成的像仍然是垂直于主轴的，如y1 。y图五。定义像的横向大小与物的大小之比值为横向放大率b ，即b =6上切点下切点左切点右切点Aa（361,618）a7 (309,538)a3 (284,609)a4 (365,572)Ba1 (433,612)a8 (400,532)Ca2 (670,592)a9 (628,517)a5 (609,588)a6 (677,554)Db2 (618,298)b6 (618,233)b2 (554,297)b3 (618,233)Ea1 (324,302)a7 (324,231)a4 (248,291)a5 (324,231)OyQP1Py1SS1O

12、1图五当离轴远近不同的物点成像时，横向放大率 b 不一定相同，这就形成了像差，称为畸变。畸变并不影响成像的清晰程度，只改变像的几何形状。例如物面是垂直于主轴且具有正方形网眼的平面如图六（a），它的像仍是平面的，但成为由曲线组成的网。当物点与主轴间的距离增大时，如果放大率随之增大，则称为枕形畸变如图六（c）；如果随之减小，则称为桶形畸变如图六（c）。图六从图四可以看出，图 A、C 的切线为直线。如图七，上、下两条直线分别为 A、C的上公切线和下公切线,应用点到直线的距离公式：设上公切线直线方程为 Ax + By + C = 0 ,B 的切点的坐标为 a (a1, a2) ，则点 a 到上切线的距

13、离 d1 为： d1 =.同理可求 d2 的长度。A2 + B2通过编程，得到 d1=0.158319448273571，d2=0.319058614126851。和原图中的 A、B、C 的公切线比较，可以看出其误差非常小。所以，可以将畸变看成”误差”。7Aa1 + Ba2 + C图七4.4 靶标像的中心坐标模型如图四，首先通过中的 regress 命令，对，(b2,b3,b4,b5), (a1,a8),(a2,a9,b1,b7), （ a6,a4,a5,a6 ） ,( a,a7,b,b6), 分别进行拟合，得出直线K1,K2,K3,K4,K5 的方程，方程如下。然后通过中的 solve 命令

14、求出上述直线的交点。也就是 A、B、C、D、E 五点在像素坐标系中的坐标，见表三。y1 = 8.1841x1 - 2083.62y2 = 4.7059x2 -1406.82y3 = 5.3112x3 - 2868.06y4 = -0.0849x4 + 608.8909y5 = -0.0063x5 + 269.2195置信区间取 0.5，概率值p 见下表p0.00950.03370.01470.0019从上表数据可以看出：通过拟合得到的直线方程无拆很小，所以得到的靶标像的中心在像素坐标系中的坐标也比较准确。表三：靶标像的中心在像素坐标系中的坐标 (以像素为)8图形ABCDE坐标(325.617,

15、581.253)(420.750,573.178)(644.353,554.198)(590.000,265.522)(287.270,267.419)4.5平移转换坐标模型（问题 1、问题 2）将像素坐标系设为 oxy,像平面的坐标系设为 oxy,如图八。YyOXox图八为方便称oxy 为旧系，oxy为新系，设o在旧系下的坐标为(x0,y0),于相机分变率为 1024×768，所以 x0=1024/2,y0=768/2.则可以得到由ìx ' = x - 1024 ,ï2í768ï y ' = y -.ïî2

16、.将 4.4 中所计算出的坐标带入上式，则可以得到靶标上圆的圆心在该相机像平面的像坐标。也就是 A、B、C、D、E 五点在像平面的像坐标，见表四。表四：靶标像的中心在像平面坐标系中的坐标(以像素为)4.6 DLT 模型（问题 3）对数码相机的检校，我们采用的方法是直接线性变换(DLT)模型法，其基本思想是在空间设置足够数量的点，事先测定每个点的空间坐标位置，用相机对点进行拍摄，将像点坐标和空间坐标代入构像方程，就可计算出该相机的内外方位元素。9图形ABCDE坐标(-186.383,197.253)(-91.250,189.178)(132.353,170.198)(78.000,-118.47

17、8)(-224.73,-116.581)DLT 方法最大的优点是解算过程不需要迭代，标定速度快。对一般精度要求来说，DLT 方法是一种比较好的标定方法。二维 DLT 与光束法平差法是首先利用二维 DLT 变换求解出内外元素初值，然后用摄影测量中常用的光束法进行相机标定，即一步法。二维直接线性变换解法的基本关系式为：ìx =a1 ´ X + a2 ´Y + a3ïa ´ X + a ´Y +1ï78ía ´ X + a ´Y + aï y = 456ïî为直接线性变换

18、系数a7 ´ X + a8 ´Y +1式中,对上式进行适当变换，求解出内外方位元素初值。这里，问题一中求得的切点作为已知点，将点的世界坐标和像素坐标坐标带入上式中，可以得到，若再将世界坐标中 A、C、D、E 圆心的坐标带入，可得到相机各参数的数值和像素平面内各圆心的坐标分别如表五、六。表五：相机各参数的数值表六：像素平面内各圆心的坐标(以像素为)数码相机在 x, y 方向上像素的大小 dX 和 dY 是一样大的，但也存在过程的误差造成在 x、y 方向上像素的大小 dX 和 dY,不一致的现象。如果没有系统误差改正数时，d、和 d5,的不等，并不影响像点的正确位置。为lyx

19、向为lx ，Y 向, x 向分辨率为 nx ， Y 向分辨假设 CCD率为 ny ,，有l， ly = nyd y10图形ACDE新像素坐标(X1,Y1)(331.312, 577.957）(638.435,544.880)(587.151,264.931)(294.857,266.620)像素坐标(X2,Y2)(325.617,581.253)(644.353,554.198)(590.000,265.522)(287.270,267.419)a1a2a3a4a5a6a7a83.3540.1518497.79750.16632.7749447.66840.0007-0.0006a1a8a1a

20、8则误差改正数 a 值表达式为:a = dy / dx = nxly / (nylx )那么，a 值的相对中误差是:ma / a = (mlx / lx ) + (mly2/ ly )2 1/2/ lx » mly / ly = ml近似地取 mlx/ l ，则有ml / l = (1/2)(ma / a)ml / l = (1/2)(d Dx '/ Dx ')可近似认为得到ml / l 必须测定精度为ml / l = (1/ 22)(d Dx '/ Dx)经检验，相对误差都在 5%以内，这说明方法的精度很高。4.7 确定两个相机相对位置模型（问题 4）4.7

21、.1 单个相机模型：图九所示的是单个相机的成像模型，其中使用了 4 个坐标系;1) 世界坐标，即图中的Ow - xw ywzw 。在世界坐标系中 3表示为 R0 .中的一点可以2) 相机坐标，即图中的O - xc yczc 。该坐标系以镜头的光心为原点，光轴为 z轴。相机坐标系上的一点可以可以表示为(xc , yc , zc ) 。3) 相机成像平面坐标，即图中的Oi XY .这是基于成像平面建立的 2 维坐标系，成像平面与光轴垂直，与光心的距离为相机的焦距 f.4) 计算机图像坐标系。这是以像素为长度的一个坐标系，图像上的某一点= ( X f ,Yf ) 。记为 Pf11图九找出3 维世界坐

22、标的点 Pw =（xw, yw, zw) 与2 维图像坐标的点 Pf=（X f ,Yf ) 的对应关系，就是单个相机标定所要解决的问题。其对应关系可以分为 4 个步骤。1.世界坐标到相机坐标的变换。由于两个坐标系都是 3 维直角坐标系，变换的形式满足Pc = R0Pw + T 0 , Pc = R0Pw + T 0（1）其中， R0 是一个3´ 3 的矩阵，反映旋转关系，T0 是一个3´1的矩阵形式，反映平移关系。2.相机坐标到成像平面坐标的变换。利用针孔模型可以将 3 维的相机坐标变换为理想的（不考虑畸变的）图像坐标，其转换关系为xvc zccX= f,Y = f（2）u

23、uzc其中 f 是相机的等效焦距。3. 真实的图像坐标到计算机图像坐标转换。其转换关系为X f = hx Xd + Cx ,Yf = hyYd + Cy(4)其中Cx , Cy 是主点在计算机图像坐标系中的坐标， hx , hy 反映了两个坐标系转换时横、纵坐标分别满足的线性关系。以上 R0 ,T0 , f , k1, k2 , hx , hy , Cx , Cy 都是单个相机标定所要求解的参数。其中 R0 ,T0 称为12外部参数，反映的是相机的空间位置， f , k1, k2 , hx , hy , Cx , Cy 称为内部参数。4.7.2 双目模型在双目相机系统中, 每个相机都分别符合

24、4.7.1 所描述的成像模型。而通过单个相机模型，可以分别知道左、右相机的图像点与 3 维世界坐标的点的对应关系。即Pw = (xw, yw, zw ) ® Pf 1 = (X f 1,Yf 1) ，Pw = (xw, yw, zw ) ® Pf 2 = (X f 2,Yf 2 ) 的对应关系已分别知道。而双目标定所需要解决的问题是两个相机坐标系之间的对应关系，即(xc1, yc1, zc1) 与yc2 , zc2 ) 之间转换关系。而由于它们分别属于 3 维直角坐2 +T '标系, 因此它们之间的转换包含有旋转以及平移的关系, 形如其中, 旋转矩阵 R以及平移向量

25、 T 就是双目定标所要求解的内容4.7.3 本质矩阵的基本原理4.7.3-1 符号说明pu1 = ( Xu1,Yu1)1. 对左相机成像平面的 2 维坐标， 对应与相机坐标系的p= ( X ,Y ) 点，将其改写为其次向量的形式u = ( X/ f ,Y/ f ,1)T ；同理，对右c1u1 u11u11 u11相机u = ( X/ f ,Y/ f ,1)Tu2u22222. 把 3 维的世界坐标亦记作齐次向量的形式 x = (Xw ,Yw, Zw,1)T3. 依照式 (1)、（2） , 世界坐标到成像平面的 2 维坐标的变换都是线性的, 因此可以把左右相机的针孔成像过程写成以下的矩阵形式:

26、u1 = A1x, u2 = A2 x ，其中，A1, A2 都是3´ 4 矩阵；4. 为下面讨论方便起见, 把 Ai 看作是一个方阵与一个向量的组合, 即A1 = (M1 | -M1T1) 和 A2 = (M 2 | -M 2T2 ) ，其中Mi 是3´ 3 的方阵，Ti 是3´1的向量；5. 运算符号：向量t = (tx ,ty ,tz ) ，定义运算T-tz0txæö0tyt´ c = ç÷-ttçx ÷zç÷-t0èøy13定义W1 »

27、W2 表示两个矩阵W1, W2 只相差一个比例因子，即W1 = kW2 ，k 为实数不难看出，这样定义的符号“ » ”满换性和传递性。4.7.3-2 本质矩阵及其性质使用 4.6.3-1 的记号方法，把左右相机成像平面的坐标分别记为向量形式，u1 = ( Xu1 / f ,Yu1 / f1,1）T , P1H » P1 '且P2H » P2 ' 。必存在一个 3´ 3 的方阵 Q ，使得 u TQu = 0 ，矩阵 Q 被称质矩阵。其有如下的性质;12定理一：若左右相机矩阵分别是 A1 = (M1- M1T1)和A2 = (M 2- M

28、2T2) ，那么相当于这两台相机的本质矩阵为Q » (M T )-1 M éM (T - T )ù´ c1 ë1û221由这个定理知道，本质矩阵 Q 可以分解成一个满秩的方阵 R 与一个称矩阵 S的乘积，即Q » RS .以下定理表明，这样的分解几乎是惟一的(最多相差一个尺度因子)。定理二：若一个3´ 3 的方阵 Q 能够有两种的分解方法Q » R1S1 » R2S2 ，其中Si 是非称矩阵， Ri 是满秩的方阵，那么必有 S1 » S2 ，并且，如果 Si = t ´ c

29、，那么» R1 + atT ,a 是某三维的向量。零R2定理三：若左右相机的相对位置固定， A , A和 A , A '对应于相同的本质矩阵 Q，'1212当且仅当存在一个满秩的4´ 4 的方阵 H，使得 P1H » P1 '且P2H » P2 ' .4.7.3-3 利用本质矩阵标定双目相机上一节的定理反映了一个事实:两台相机的相对位置关系与本质矩阵几乎是等价的。即给定了两台相机的相对位置关系 R 和 T，就给定了系统的本质矩阵 Q；反之，给定了本质矩阵 Q，两台相机的旋转平移关系 R 和 T 也几乎确定了, 最多相差一个

30、比例常数因子。“几乎”的意思。就 R和 T 的确定还需要标定一个比例因子 k 。这里给出一个方便的通过Q 求解 R 和T 的方法, 即先把 Q 进行特征值分解（SVD），得Q = UDVT ,-1 00æ 00 ö, 其中 G = ç 10 ÷t´= VZVT根据 4.7.3-2 定理 一 ，得到 R = UGVT ，ç÷ç 01 ÷èø14-1æ 0Z = ç 10 ö00 ÷,而 t 是把平移矩阵 T 的 3 个元素看成向量。÷&#

31、247;ø如此, 我们确定了旋转和平移关系的基本形式 R , T 。而 R,T 与 R , T相差一个比例因子 k ,即 R=kR,T=kT, 可以借助世界坐标中已知距离的两个点' p = ( X ,Y , Z ) 和 p = ( X '',Y '', Z '') '''''wwwwwwww4.7.4 标定过程标定的具体步骤如下：1.利用传统的方法分别对左右两台相机的内参数进行标定。点通过靶标图案提取，该图案沿平行于光轴方向多次放置 ,精确下各次放置的位置。2. 把已完成单目标定的两台相机

32、安装到实际系统中,固定。3. 分别在左右像中确立匹配像素对 ,列出方程。4. 根据 4.7.3-3 的原理求出两个相机的相对位置 R和 T, 待定比例因子 k 。5. 最后, 标定比例因子 k 。具体做法为, 首先, 设定两个标志点, 测出它们的距离 d0 ; 其次, 通过 R和 T求出其在相机坐标系下的距离 d, 因为 R和 T是 k的线性函数, 所以 d也是 k 的线形函数，即 d ' = f (k) 。解一个一次的线性方程f (k) = d0 ，即可得到比例因子 k 。五、模型的评价与改进问题一的模型采用做切线求交点得圆心的方法，计算比较准确，但是对图像要求比较高，要求图像边缘比

33、较平滑。本题中图像比较小，图像边缘有锯齿，这增大了求切点的误差。对模型进行改进，应先考虑使图像边缘平滑，尽可能消除锯齿，减小误差。问题三模型 DLT 模型没有迭代，计算快，测量效率高，但精度并不是很高。对模型的改进，应采用精确的参考点或增加参考点的数目，对标定进行处理，这样就会提高精度。问题四本质矩阵标定模型不仅可以充分考虑到立体视觉系统中每个相机存在的非线性失真，具有与传统方法相当的精度，同时双目标定不限于实验环境，具有简化标定的灵活性。考文献:宋卫东，几何，北京：高等教育参12，2003,12:208209。Rafael C. Gonzalez,Richard E. Woods 著，Ste

34、ven L. Eddins,阮秋琦等译，数字图像处理（版），北京：电子工业，2005,9:17、289290。3 姚启钧 原著，华东师大光学2005,279280、186187。4 郑阿奇 主编，曹弋 编著. 2004,108109。编写组改编,光学，北京：高等教育，实用（第 2 版），北京：电子工业，5 谭晓军,余只,李军，一种改进的立体摄像机标定方法，测绘图报，第 35 卷，第二期：13814216附录：附录 1：问题 2 中所用到的程序clear clcA=imread('clip002.jpg'); for i=1:768for j=1:1024if (A(i,j)&g

35、t;=180)A(i,j)=0;elseA(i,j)=255;endend end c=edge(A);for i=1:768cc(769-i,:)=c(i,:);endload cc pd=zeros(5,8);%* p1=323,579;p2=640,555;pp1=0,0;pp2=0,0;AC*while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=282:365 for j=538:620if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(

36、1)+p1(2)-j)<0)p1=i,j;endendendendfor i=604:677 for j=517:593if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)<0)17p2=i,j;endendendendend pd(1,1:2)=p1;pd(3,1:2)=p2;%* p1=323,579;p2=640,555;pp1=0,0;pp2=0,0;while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=

37、282:365 for j=538:620if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)>0)p1=i,j;endendendendfor i=604:677 for j=517:593if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)>0)p2=i,j;endendendendend pd(1,3:4)=p1;pd(3,3:4)=p2;%* p1=583,265;p2=285,266;pp1=0,0;pp2=0,0;DE*while(p1(1)=

38、pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;18for i=549:618 for j=233:298if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)<0)p1=i,j;endendendendfor i=246:324 for j=231:302if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)<0)p2=i,j;endendendendend pd(4,1:2)=p

39、1;pd(5,1:2)=p2;%* p1=583,265;p2=285,266;pp1=0,0;pp2=0,0;while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=549:618 for j=233:298if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p1(2)/(p2(1)-p1(1)+p1(2)-j)>0)p1=i,j;endendendendfor i=246:324 for j=231:302if(cc(j,i)=0)if(i-p1(1)*(p2(2)-p

40、1(2)/(p2(1)-p1(1)+p1(2)-j)>0)p2=i,j;endend19endendend pd(4,3:4)=p1;pd(5,3:4)=p2;%* p1=323,579;p2=285,266;pp1=0,0;pp2=0,0;AE*while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=282:365 for j=538:620if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)>0)p1=i,

41、j;endendendendfor i=246:324 for j=231:302if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)>0)p2=i,j;endendendendend pd(1,5:6)=p1;pd(5,5:6)=p2;%* p1=323,579;p2=285,266;pp1=0,0;pp2=0,0;while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=282:365 for j=538:620i

42、f(cc(j,i)=0)20if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)<0)p1=i,j;endendendendfor i=246:324 for j=231:302if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)<0)p2=i,j;endendendendend pd(1,7:8)=p1;pd(5,7:8)=p2;%* p1=640,555;p2=583,265;pp1=0,0;pp2=0,0;CD*while(p1(1)=pp1(1)|(p1(2)=pp

43、1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=604:677 for j=517:593if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)>0)p1=i,j;endendendendfor i=549:618 for j=233:298if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)>0)p2=i,j;endendendendend21pd(3,5:6)=p1;pd(4,5:6)=p

44、2;%* p1=640,555;p2=583,265;pp1=0,0;pp2=0,0;while(p1(1)=pp1(1)|(p1(2)=pp1(2)|(p2(1)=pp2(1)|(p2(2)=pp2(2) pp1=p1;pp2=p2;for i=604:677 for j=517:593if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)-p1(2)+p1(1,1)-i)<0)p1=i,j;endendendendfor i=549:618 for j=233:298if(cc(j,i)=0)if(j-p1(2)*(p2(1)-p1(1)/(p2(2)

45、-p1(2)+p1(1,1)-i)<0)p2=i,j;endendendendend pd(3,7:8)=p1;pd(4,7:8)=p2;%* p1=423,572;for i=383:464 for j=532:612if(cc(j,i)>0)B*if(pd(3,2)-pd(1,2)*p1(1)-(pd(3,1)-pd(1,1)*p1(2)>(pd(3,2)-pd(1,2)*i-( pd(3,1)-pd(1,1)*j)p1=i,j;endendendend pd(2,1:2)=p1;22%* p2=423,572;for i=383:464 for j=532:612if(

46、cc(j,i)>0)if(pd(3,4)-pd(1,4)*p2(1)-(pd(3,3)-pd(1,3)*p2(2)<(pd(3,4)-pd(1,4)*i-( pd(3,3)-pd(1,3)*j)p2=i,j;endendendend pd(2,3:4)=p2;w_up=(pd(3,2)-pd(1,2)*(pd(2,1)-pd(1,1)+(pd(3,1)-pd(1,1)*(pd(1,2)-pd(2,2)/sqrt(pd(3,2)-pd(1,2)2+(pd(3,1)-pd(1,1)2);w_down=(pd(3,4)-pd(1,4)*(pd(2,3)-pd(1,3)+(pd(3,3)-

47、pd(1,3)*(pd(1,4)-pd(2,4)/sqrt(pd(3,4)-pd(1,4)2+(pd(3,3)-pd(1,3)2);clear;clc; load pd%*cir=1,3;1x=pd(cir(1),5),pd(cir(1),7),pd(cir(2),5),pd(cir(2),7)'y=pd(cir(1),6),pd(cir(1),8),pd(cir(2),6),pd(cir(2),8)' X=ones(4,1),x;b,bint,r,rint,stats=regress(y,X,0.5); par(:,1)=b;p_stats3(1)=stats(3); x=2

48、00:700;y=b(1)+x.*b(2); figure(2) plot(x,y);hold on; axis(1,1024,1,768); x=;y=;%*cir=4,5;2x=pd(cir(1),5),pd(cir(1),7),pd(cir(2),5),pd(cir(2),7)'y=pd(cir(1),6),pd(cir(1),8),pd(cir(2),6),pd(cir(2),8)' X=ones(4,1),x;b,bint,r,rint,stats=regress(y,X,0.5); par(:,2)=b;p_stats3(2)=stats(3);23x=200:70

49、0;y=b(1)+x.*b(2); plot(x,y); axis(1,1024,1,768); x=;y=;%* cir=1,5;3x=pd(cir(1),1),pd(cir(1),3),pd(cir(2),1),pd(cir(2),3)'y=pd(cir(1),2),pd(cir(1),4),pd(cir(2),2),pd(cir(2),4)' X=ones(4,1),x;b,bint,r,rint,stats=regress(y,X,0.5); par(:,3)=b;p_stats3(3)=stats(3); y=200:700;x=(y-b(1)./b(2);plot(x,y); axis(1,1024,1,768); x=;y=;%*cir=3,4;4x=pd(cir(1),1),pd(cir(1),3),pd(cir(2),1),pd(cir(2),3)'y=pd(cir(1),2),pd(cir(1),4),pd(cir(2),2),pd(cir(2),4)' X=ones(4,1),x;b,bint,r,rint,stats=regress(y,X,0.5); par(:,4)=b;p_sta