2008全国数学建模竞赛A题一等奖论文.pdf

1 基于切线特征的数码相机定位摘要本文依据成像原理，通过应用图像处理和物、像的不变性质来确定坐标。采用了几何线性定标模型来标定相机的相对位置。对于问题一，本文建立了基于图像处理的圆心搜索模型，模型在对公切线定点可行性的论证基础上，建立了图像圆心搜索算法，从而确定靶标上圆的圆心在该相机像平面的像坐标。对于问题二，运用了Matlab技术对具体图像进行去噪处理，按照圆心搜索算法对任一个圆相对其余各圆求解圆心在像坐标上的多个估计值，计值求均值消除误差后作为圆心在该相机像平面的像坐标。具体结果见下表（以像素坐标系的结果见正文表3）：圆 A 圆 B 圆 C 圆 D 圆 E 图像坐标系x 坐标值-54.2857-49.51-44.67 31 31.35 y 坐标值-50.1455-24.46 34.34 19.74-59.25 对于问题三，基于三个在同一条直线上的圆心其像在像平面坐标仍在一条直线上，通过相机拍摄试验分别找到不同相机分辨力不同角度下的相片对其处理，检验三个圆心两两连线的夹角，得到误差偏角。误差偏角的均值反应了模型的精度，方差变化影响稳定性。这里求得误差偏角的均值在0.02，方差在 0.002 以内。对于问题四，建立几何线性定位模型，得到两个相机的外部参数。通过两个相机外部参数求解得到两个相机间的旋转矩阵，平移向量及距离见正文表6。关键词：几何模型算法切线几何线性定标圆心坐标2 一、问题重述数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。最常用的定位方法是双目定位，即用两部相机来定位。该方法的基本思想是首先用两部固定于不同位置的相机摄得物体的像，分别获得物体上一个特征点在两部相机像平面上的坐标；然后根据两部相机精确的相对位置，利用几何的知识确定该点的位置。对于双目定位，精确地标定两部相机的相对位置就是关键。标定的一种方法是设计一个靶标，靶标可以由若干个圆组成，同时用这两部相机照相，分别得到各圆圆心在他们像平面的像点，根据像点的几何关系就可以得到这两部相机的相对位置。现在设计靶标如下：取1 个边长为 100mm 的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。以 AC 边上距离 A 点 30mm 处的 B 为圆心，半径作圆。同时给定了一部固定相机所摄的像。现就相机定位理论以及给定的靶标解决以下几个问题：1.建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,这里坐标系原点取在该相机的光学中心，yx平面平行于像平面；2.根据问题给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标,该相机的像距（即光学中心到像平面的距离）是 1577个像素单位(1 毫米约为 3.78 个像素单位)，相机分辨率为 1024768；3.设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；4.建立用此靶标给出两部固定相机相对位置的数学模型和方法。二、模型假设与符号说明1.模型的假设：1）.问题中所给的图片真实可靠。2）.假设相机的摄像过程中，靶标同一条直线上的所有点的投影像仍然在同一条直线上。不考虑非线性畸变的影响。2.符号说明：:A图像信息储存矩阵.iB 第i个圆在图像中的信息矩阵.:wwwwzyxo世界坐标系.:cccczyxo照相机坐标系.:oxy图像坐标系.:ouv 像素坐标系.:固定照相机偏航角.:固定照相机俯仰角.:固定照相机旋转角.:R旋转矩阵，属于相机的外部参数.:T平移向量，属于相机的外部参数.3 三、问题分析所谓照相机的定位就是确定照相机内部和外部参数。而相机参数的确定关键在于点在像平面的坐标值的计算。点在像平面的坐标值越精确，相机参数的计算越理想，标定结果就越优。对于问题一，考虑到任意两个远离的圆的两条平行的公切线的切点在同一个圆上的连线必过圆心这一性质，利用三个圆间的关系在每个圆内作两条直线，就可以对圆心进行标定，又考虑到照相机成像的过程可以用针孔摄像机模型进行模拟，所以只要找寻同一个圆在像坐标平面对应的标定点通过平面两条不平行直线有且只有一个交点即可求到靶标上圆的圆心在该相机像平面的像坐标。对于问题二，对图 2和图 3 按问题一的算法思路进行图像识别，得到任意两个圆的公切线确定的切点所在直线，有 4 或 3 条。对其进行两两组合，若有 4 条则可得到624C种直线相交的方式。即有6 个交点。这 6 个点都是靶标圆心在像平面的像坐标。对其进行平均处理即可求解。对于问题三，考虑三个在同一条直线上的圆心在像平面的投影点仍在一条直线上，可以通过相机拍摄试验分别找到不同相机分辨率不同角度下的相片对其处理求取圆心在像平面的坐标。检验三个圆心两两连线的夹角，通过角度的大小及变化关系讨论模型的精度及稳定性。对于问题四，数码相机的定位过程就是相机内外参数的确定的过程。模型可以从摄像坐标系转换考虑，建立照相机成像过程中最主要的三个坐标系：世界坐标系、照相机坐标系和图像坐标系。找到各个坐标系的转换关系即可以建立靶标圆心坐标和图像圆心坐标的线性关系。求解线性方程，即可得到相机的各种参数。四、模型的建立与求解1.问题一基于图像处理的圆心搜索模型的建立从问题分析的思路，靶标上圆的圆心在相机像平面的像坐标的确定过程应该包括模型思路可行性的讨论以及具体算法的研究两个步骤。1.1.模型思路可行性讨论：分析靶标上圆心在相机像平面上的像坐标，可以从靶标圆心的确定出发，通过圆的一般性质确定圆心的具体位置。然后考虑相机摄像的线性畸变特性（即任一条直线通过相机摄像后仍然为直线），确定圆心的两条直线在像坐标仍然为两条相交的直线。1）.靶标上圆心的几何标定：在平面几何中，在同一个平面内，两条相互不平行的直线有且只有一个交点。所以为了标定圆心的位置，至少要在圆内找到两条以上的过圆心的直线。考虑两点确定一条直线，所以至少要三个圆心不在同一条直线上的圆两两作平行公切线才能确定某个圆的圆心。在靶标中任意取三个圆心不在同一条直线上的圆，分别对其两两求平行的公切线，如下图所示：4 图 1.靶标上圆的圆心标定示意图分析图形可得，任一个圆上平行切线对应的切点连线在一条直线上并且经过圆心。（可以先证明四边形BBAA11和33AABB分别为矩形，再证明 A1AA3为平角。）故当考虑三个以上圆两两作平行公切线时，任意考虑一个圆A，其与圆 B和圆 C所作两条平行公切线对应的切点连线A1A3和 A2A4必然都过圆 A的圆心 A。根据两点确定一条直线，可以准确确定 A圆的圆心。通过上述方式可以对各个圆的圆心进行标定。2）.靶标上圆与切线关系的投影不变性讨论：照相机的摄像过程可以通过针孔摄像机模型进行模拟1，针孔摄像机摄像过程满足光线的直线传播的性质。任意一条直线的经针孔摄像后在像平面仍然是直线。所以对于任意一个圆在物平面的切线，投影到像平面的像必然是靶标圆的像的切线。这里可以用反证法进行证明。对于实物平面任意一个圆m，直线 l 为它的切线。假如它们经过照相机成像后不再是相切的关系。只有两种情况，直线的像与圆的像远离或者相割，如果两者远离时，则物平面的切点必然对应像平面的两个不同点，一个在圆的像上，另一个在直线的像上。而根据针孔成像过程中光线的直线传播过程，这是不可能的。同理可以证明相割也是不合理的。排除以上两种情况即得证结论。根据以上证明可以对标定实平面的两条过切点的直线的像进行处理，其交点即为所求的圆心的像平面坐标。1.2.圆心搜索算法的研究确定图像中的各圆的圆心的像坐标首先可以读取图片得到像素平面上各个圆的信息，然后找寻靶标平面切线的像的关系，最后将像素平面转换到图像坐标系，具体算法思想如下：:1Step创建二维矩阵 A，存储灰度化后的照相机照摄到的像平面信息。:2Step对图像信息矩阵 A进行边缘检测，得到像平面中各圆的边缘线。:3Step创建圆信息矩阵nBBB21,，储存各个圆的投影在图像中的边缘信息。n 为靶标上圆的个数。:4Step取两个圆的信息矩阵jiBB,，通过全局搜索找寻圆上平行公切线对应公切点所确定的直线，圆心必在这条直线上。当求得两条以上此直线时即可以确定圆心。:5Step将像素坐标系转换到以光心为原点的像平面坐标系。4Step的像平面公切线具体算法步骤：1.取两个圆的信息矩阵jiBB,，设iB 中有1n 个像素，jB 中有2n 个像素，分别从中各C4 B2 C3 C A1 B1 A B A2 A3 A4 B3 B4 C1 C2 5 取 一 个 点),(rriyxB和),(kkjyxB，作 出 过 这 两 点 的 直 线l：0cbyax。.1,.121nknr2.遍历iB 上的点),(rrrriyxB，令cybxarrKrrrr*)(1。1.1 nrr3.遍历jB 上的点),(kkkkjyxB，令cybxakkKkkkk*)(2。2.1 nkk4.令)min(),min(min(*)max(),max(max(2121KKKKH，若,0H则直线l为圆jiBB 与的外公切线，点),(rriyxB和),(kkjyxB为所求的切点。5.由同一个圆上的两个公切点作直线。则圆心必在这条直线上。算法具体流程图见附录图。2.问题二的求解：2.1.图像的预处理由于图 3 在拍摄过程会受外界噪声等因素的影响，所以在图像读取的过程中，个别点像素值出现奇异值，对问题的精确求解产生比较大的干扰。图像处理中减少噪声的常用方法主要是图像的平滑，而图像的平滑处理一般由噪声本身的特性而定。这里考虑噪声点的随机性比较大，并且两个噪声信号间的关系脉宽比较大，调用Matlab命令2m edfilt利用十字丝窗口对图像中的点进行中值滤波进行去除噪声处理2。2.2.对图 2 进行物平面的标定图 2.圆 B与其它圆的公切线关系以圆 B为例，分别画出圆B与其它各圆的公切线如上图所示。一般情况下，与圆B相切的平行的外公切线有4 对，由于圆 A,B,C 的圆心在同一条直线上，所以这里只有3对。由这 3 对公切线在圆 B上的切点可以确定三条直线。由问题一的证明这三条直线必然相交于圆心。在像平面内，在圆的像内的三条直线从理论上也是相交于一点的。此点就是圆B的圆心在像坐标系的像点。2.3.基于问题一算法的求解：对预处理后的图像进一步处理，按照问题一算法的思路，得到5 个圆信息矩阵，按照4Step的像平面公切线具体算法步骤求解5 个圆两两间的公切线与圆的切点确定的直线。根据各直线间的相交关系，得到每个圆圆心在像平面的像素坐标估计值。程序见附录。这里以圆 A为例得到各种情况下的坐标值：AB CED1A3AA4 4B2B1B5B6B2A6 表 1：圆 A在像素坐标系的坐标值作用圆BE-A BD-A CE-A CD-A ED-A 平均值u坐标188.52 186.83 188.52 186.83 188.29 187.8 v坐标322.6 322.95 322.6 322.95 321.14 322.45 注：表格“作用圆”行中，横杆以前的两个字母表示确定公切线关系的两个圆，横杆后面的点表示所求圆的圆心。如BE-A表示圆 B 和 E与圆 A作平行公切线得到切点的确定的坐标值，这里为（188.52，322.6）。从表格可以看出，不同情况下求得的圆心的像素坐标值是有波动的。为了消除这种波动误差，这里对不同情况下求得的圆心坐标进行平均，即：nuuii/，nvvii/，即可得到 A 圆圆心在像坐标系的坐标值（见上表格中的平均值）。像素坐标系的坐标值是建立在像素单元的基础上，坐标原点是图像的左上角，按照问题的要求，靶标圆的圆心在像平面的坐标值，应该是原点在相机的光心，x-y 平面平行于像平面的坐标系上的值。根据题设给出的已知条件可以得到：单个像素单元的宽度是mmdydx78.3/1像坐标系的原点为dxvu512,384,00；坐标转换公式为：00vvdyyuudxx根据以上坐标系转换公式得到像平面的坐标值为表 2：表 2：圆 A在像平面坐标系的坐标值作用圆BE-A BD-A CE-A CD-A DE-A 平均值x 坐标-54.10-54.54-54.10-54.54-54.16-54.29 y 坐标-50.11-50.01-50.11-50.01-50.50-50.15 表格中“作用圆”行的解释与表格1 相同。这里单位为 mm。用同样方法可以得到各个圆的圆心坐标如下表：表 3：其余各圆在像平面坐标系的均值坐标值圆 B 圆 C 圆 D 圆 E u 坐标196.87215.15501.17502.52v坐标419.53641.79586.62288.02x 坐标-49.51-44.673131.35y 坐标-24.4634.3419.74-59.25这里yx,坐标的单位为 mm。综合考虑上面两个表格，由于像素坐标系原点从图像的左上角开始，靶标平面圆心在一条线上的A，B，C三个圆心坐标在像素平面上依次增大，呈现在图像上是依次下降的过程，这与照相机图像中圆位置变化是吻合的。其余相关数值见文章附录表1 到表 4。3.问题三检验方法的确定与讨论本文采用试验的方法进行检验。建立靶标如下图 3 圆心坐标7 图3.试验靶标其中，A,B,C三个圆的圆心在同一直线上，圆D 是求解 A,B,C三个圆的圆心在相机像平面的像素坐标的参考圆。利用模型一求出A,B,C三个圆的圆心在相机像平面的像素坐标),(yxaa，),(yxbb，),(yxcc，模型算法的误差导致三个像坐标不在同一直线上，CA 与CB 之间存在一个夹角，定义该夹角为误差偏角。用三个不同分辨率的数码相机分别从不同的角度拍摄三张靶标的像，部分图像（不同分辨率图像各附一张）见下图4 分辨率为 640 480 分辨率为 1600 1200 分辨率为 2048 1536 图4.不同分辨率图像通过问题一模型算法求得),(),(),(yxyxyxccbbaa进而求得误差偏角，具体数据如下：表4：误差偏角的数据31 对试验数据的处理：分别对三种不同分辨率条件下的三个图像所对应的三个误差偏角求平均值和方差得到如下结果：表 5：误差偏角的平均值与方差分辨率误差偏角的平均值误差偏角的方差640 480 0.011 0.001242 1600 1200 0.006 0.001327 2048 1536 0.004 0.001196 32 精度及稳定性分析：模型的精度是指模型的测量值与实际值之间的偏差，本文问题一的模型的精度可由误差偏角的大小来反应；而稳定性是指在某些参数改变的情况下，而不影响模型的适用性，该模型可由误差偏角的方差来反映。分辨率640 480 1600 1200 2048 1536 图像1 2 3 1 2 3 1 2 3 误差偏角0.012 0.011 0.011 0.007 0.006 0.006 0.004 0.004 0.005 A B C D 8 分析上述试验数据，可以得到模型算法的误差偏角控制在0.02 弧度以内，充分说明模型算法的精度是比较高的。针对模型的精度随分辨率增大而提高的分析，我们所建立的模型的关键是对拍摄的图像的处理，当图像的分辨率增大时，其必然提高我们图象处理的精度，误差偏角必然会下降，因此该试验数据对模型准确性的验证提供了一个有利的依据。针对稳定性的分析，从试验数据中我们可以发现在任何一种分辨率的条件下误差偏角的方差均控制在0.002 以内并且方差相对比较接近，这说明在任何一种分辨率的条件下模型算法是相对稳定的。4.问题四模型的建立与求解：问题四是对两个相机相对位置的标定过程。相机位置标定充分展现了二维图像提取三维空间信息的过程。根据光线的直线传播原理，即透镜成像的基本原理，空间点通过相机在像平面上所成像点的位置与其本身位置之间存在着一定的几何关系。具体分析这种关系既可得到相机成像的几何模型。1.几何线性定标模型的建立1，3空间点通过相机的成像的几何模型其实就是空间三个坐标系的转换过程。该几何模型的参数分为由相机内在属性决定的内参数和相机坐标系与世界坐标系的位置关系决定的外部参数。而这些参数一般都是可以确定的。1.1 各坐标系的确定。在照相机成像的几何模型中，主要包括三个主要坐标系：世界坐标系、照相机坐标系和图像坐标系。世界坐标系(wwwwzyxo)是为了描述任意位置的相机和实物而假象的一个坐标系，一般用作参考，可以根据具体情况来选取。照相机坐标系(cccczyxo)是安置固定相机的坐标系，与相机的位置有关。原点定义在相机的光心上，cz 轴为光轴，cx 和cy 与像平面的水平轴和垂直轴相平行。图像坐标系（oxy）是经过数模转换器转换后的数字图像（这里考虑灰色图像，若是彩色图像有RGB三个颜色通道，则为三维坐标系）。图像坐标系的坐标原点为光轴与像平面的交点，x和y轴即为像平面的水平轴和垂直轴。1.2.从世界坐标系到照相机坐标系的变换过程：从世界坐标系到照相机坐标系的变换过程考虑固定相机的角度变化。主要包括偏航角，俯仰角和旋转角。由这三个角度可以确定这两个坐标系的旋转矩阵R。综合考虑两个坐标系间的水平距离关系，可以得到两个坐标系的关系如下：TzyxRzyxwwwccc1101wwwTccczyxTRzyx（1）其中ccczyx为照相机坐标系的三个坐标值，wwwzyx为世界坐标系的三个坐标值R 为两个坐标系间的旋转矩阵，可得)()()(xyzRRRR9 coscoscossinsincossinsinsincoscoscossinsinsinsincossinsincossincossinsinsincoscos，而321tttT为它们之间的平移向量。1.3.从照相机坐标系到图像坐标系的变换过程照相机的摄像过程，即从照相机坐标系到图像坐标系的变换可以用针孔摄像机来模拟。针孔摄像机摄像的过程其实就是实物各个方向上的大小同比例缩小的过程。对于一个圆通过针孔摄像机的过程可以从下图看出：图 3.针孔摄像的原理从图 1 中任取一点 P，其在世界坐标系中的坐标为)(wwwzyxP,对应在照相机坐标系中的坐标值为)(ccczyxP。由于光线的直线传播，P在像平面的投影点P各方向的尺寸是对应成比例的。不妨将图形变换如下：图 4.对 P点变化后的图形从上图可以看出：cccyyxxzf（2）其中f为光心到像平面的距离。在以上基础上可以得到P在图像上的像素坐标),(vu为：00vdyyvudxxu（3）cxcocyO czx y f ABB A 10 其中，dx，dy是每一个像素在图像平面x轴和y轴上的物理尺寸，),(00vu为摄像机光轴与图像平面的交点。由式)3)(2)(1(同时考虑图像中u轴和v轴的倾斜，可得到点P在世界坐标系下的坐标与其投影点P的图像坐标),(yx的关系，如下式所示：)4(1100100000012100MXXMMZYXTRvuvuZwwwc其中，001,/,/vudydxfMdyfdxf由，决定，由于这些参数只与摄像机内部结构有关，故这些参数称为相机内部参数;2M 则只与旋转矩阵R和平移向量 t有关，它是完全由相机相对于世界坐标系的位置关系决定的，因此被称为相机外部参数。而M为33的矩阵，称为投影矩阵。2.模型的求解分析由前面的几个问可以求出靶在像平面的坐标，而00,vudydxf这些变量在题设中已经给出，所以本问题的参数的求解只需求解外部参数R和 T。对公式（4）进行线性展开，即可求解出2M 矩阵的各个元素。由于这里最多只有6 个变量，所以最多找到三个点就可以求解 R，T 的值。当对两个照相机进行定位时，相同相机的内部参数都是一样的，由于相机的位置不同其外部参数也是不同的。当求出两个相机的外部参数时，两个相机间的位置就可以通过 平移 向 量 进 行 求 得。假 设两 种 情 况 下 的 旋 转 矩 阵和 平 移 向量 已 知，分 别为10111TRM和10222TRM则两个相机照相机坐标系的坐标值都可以表示为111wwwccczyxMzyx（5）112wwwccczyxMzyx（6）对（6）式进行变换得：)1()1(12cccwwwzyxMzyx（7）综合式)7)(5(可得：)1()1(121cccccczyxMMzyx令21MMM,且10TRM，1021TRM则1011110101011TTRRRTRTRTRM11 相机间的距离即为向量T 的模。综上所述，两个相机间的旋转矩阵RRR1；平移量为TTRT1；其中TR,；相机间的空间距离即为平移量T 的模可由 M2求逆求解。表 6：相机定位参数分析表格定标的参数公式备注两个相机间的旋转矩阵RRRR11R 为一个相机的旋转矩阵；R由另一个相机外部参数矩阵逆矩阵求到平移量TTTRT11T 为一个相机的旋转矩阵；T由另一个相机外部参数矩阵逆矩阵求到空间距离|T即为平移量的模五、模型评价1.模型优点：1）.模型通过公切线定圆心的思路很好的解决了像坐标系中圆心确定问题，为以后的模型求解奠定了基础2）.模型试验充分考虑问题假设，与前面衔接较严密。采用了相机拍摄试验的方式对模型算法进行分析，比较合理。2.模型的不足：通过相机摄像试验是选取的样本数量有限，试验的精度有一定的局限。参考文献1 胡凌山，移动机器人双目立体视觉技术研究J,2005.1。2 夏良正，李久贤，数字图像处理（第2 版）M，东南大学出版社，2005.8。3 唐志豪，基于双目立体视觉的测量技术研究J，2006.12。4 姜启源,谢金星数学模型第三版 M 北京：高等教育出版社，20035 苏金明,阮沈勇,王永利，MATLAB工程数学，北京：电子工业出版社，2005。12 附录Start;0,0;0,0kkrrkr;1rr;1kk过点),(rriyxB),(kkjyxB做直线l;1rrrr:1kkkk计算)(1rrK计算)(2kkK?2nkk)max(1K，)min(1K)max(2K，)min(2K?0H?2nkl是外公切线,储存切点?1nrEnd?1nrrN N Y N N N 13 附图：圆心搜索算法流程图附表：表 1 作用圆AE-B AD-B CE-B CD-B DE-B 平均值u 坐标195.1 198.82 195.1 198.82 196.5 196.87 v 坐标419.2 418.36 419.2 418.36 422.5 419.53 x 坐标-49.9735-48.9894-49.9735-48.9894-49.6032-49.5053 y 坐标-24.5503-24.7725-24.5503-24.7725-23.6772-24.463 表 2：作用圆AE-C AD-C BE-C BD-C DE-C 平均值u 坐标215.61 214.88 215.61 214.88 214.78 215.15 v 坐标641.85 642.03 641.85 642.03 641.22 641.79 x 坐标-44.5476-44.7407-44.5476-44.7407-44.7672-44.6693 y 坐标34.351 34.39947 34.35185 34.39947 34.18519 34.33598 表 3：作用圆AB-D AC-D AE-D BC-D BE-D CE-D 平均值u 坐标503 503.49 495.7 503.48 497.25 504.1 501.17 v 坐标582.5 581.87 591.86 581.74 591.55 590.18 586.62 x 坐标31.48 31.61 29.550 31.60 29.96 31.77 30.99 y 坐标18.650 18.484 21.126 18.449 21.044 20.68254 19.74 表 4：作用圆AB-E AC-E AD-E BC-E BD-E CD-E 平均值u 坐标503.11 502.63 502.77 500.31 501.64 504.64 502.52 v 坐标293.59 286.04 288.27 283.85 288.47 287.93 288.02 x 坐标31.510 31.38 31.42 30.76 31.12 31.91 31.35 y 坐标-57.78-59.77-59.18-60.35-59.13-59.27-59.25 附录程序：clear;clc 唐俨 数学建模论文2008b.jpg);%imshow(A);A=rgb2gray(A);%A1=edge(A,roberts);a1 b1=size(A);for i=1:a1 for j=1:b1 if A(i,j)=120 A1(i,j)=1;else A1(i,j)=0;14 end end end A2=edge(A1,roberts);imshow(A2);a,b=find(A2);n1=0;n2=0;n3=0;n4=0;n5=0;for i=1:length(a)if a(i)140&a(i)260&b(i)140&a(i)370&b(i)140&a(i)580&b(i)450&a(i)220&b(i)450&a(i)530&b(i)=0;m=m+1;b12(m,:)=B1(i,:);b21(m,:)=B2(j,:);end end end m m=0;x1=0;x2=0;x3=0;x4=0;f1=0;f2=0;for i=1:length(B1)for j=1:length(B3)x1=B1(i,1);y1=B1(i,2);x2=B3(j,1);y2=B3(j,2);for ii=1:length(B1)x3=B1(ii,1);y3=B1(ii,2);f1(ii)=(y3-y1)*(x2-x1)-(y2-y1)*(x3-x1);end ma1=max(f1);mi1=min(f1);for jj=1:length(B3)x4=B3(jj,1);y4=B3(jj,2);f2(jj)=(y4-y1)*(x2-x1)-(y2-y1)*(x4-x1);16 end ma2=max(f2);mi2=min(f2);MA=max(ma1,ma2);MI=min(mi1,mi2);if MA*MI=0;m=m+1;b13(m,:)=B1(i,:);b31(m,:)=B3(j,:);end end end m m=0;x1=0;x2=0;x3=0;x4=0;f1=0;f2=0;for i=1:length(B1)for j=1:length(B4)x1=B1(i,1);y1=B1(i,2);x2=B4(j,1);y2=B4(j,2);for ii=1:length(B1)x3=B1(ii,1);y3=B1(ii,2);f1(ii)=(y3-y1)*(x2-x1)-(y2-y1)*(x3-x1);end ma1=max(f1);mi1=min(f1);for jj=1:length(B4)x4=B4(jj,1);y4=B4(jj,2);f2(jj)=(y4-y1)*(x2-x1)-(y2-y1)*(x4-x1);end ma2=max(f2);mi2=min(f2);MA=max(ma1,ma2);MI=min(mi1,mi2);if MA*MI=0;m=m+1;b14(m,:)=B1(i,:);b41(m,:)=B4(j,:);end end end m=0;x1=0;x2=0;x3=0;x4=0;f1=0;f2=0;for i=1:length(B1)for j=1:length(B5)x1=B1(i,1);y1=B1(i,2);x2=B5(j,1);y2=B5(j,2);17 for ii=1:length(B1)x3=B1(ii,1);y3=B1(ii,2);f1(ii)=(y3-y1)*(x2-x1)-(y2-y1)*(x3-x1);end ma1=max(f1);mi1=min(f1);for jj=1:length(B5)x4=B5(jj,1);y4=B5(jj,2);f2(jj)=(y4-y1)*(x2-x1)-(y2-y1)*(x4-x1);end ma2=max(f2);mi2=min(f2);MA=max(ma1,ma2);MI=min(mi1,mi2);if MA*MI=0;m=m+1;b15(m,:)=B1(i,:);b51(m,:)=B5(j,:);end end end x=b12;b14;O124=zuobiao(x);x=b12;b15;O125=zuobiao(x);x=b13;b14;O134=zuobiao(x);x=b13;b15;O135=zuobiao(x);x=b14;b15;O145=zuobiao(x);O1=(O124+O125+O134+O135+O145)/5 function c=zuobiao(x)k1=(x(2,2)-x(1,2)/(x(2,1)-x(1,1);k2=(x(4,2)-x(3,2)/(x(4,1)-x(3,1);aa=k1,-1;k2,-1;bb=-x(1,2)+k1*x(1,1);-x(3,2)+k2*x(3,1);syms u v;c=u,v;c=aabb;