数学建模全国一等奖论文-供参考

关于数码相机定位问题的数学模型
摘要
本文是关于系统标定中的双目定位问题，对靶标上面的特征点在像平面内进行标定。

我们首先对像图进行了灰度处理，然后通过sobel算子进行边缘检测，将像平面转化为边界处理，计算区域由二维降低为一维，显著的降低了计算量。

对于问题一，我们根据成像原理，得到像图中某些特征点和特征关系用以确定圆心，并由此建立了割线逼近模型，并且利用遍历搜索的算法进行求解。

同时，我们根据图像的边界点特征，建立了霍夫变换模型和基于K-平均算法的数学模型。

割线逼近模型主要依据物与像之间的不变的几何关系，通过寻找靶标的像图中的切线和切点来确定圆心，在求解过程中使用了灭点理论。

这种算法虽然在边界识别以及切线确定中存在误差，但它可以有效的避免对像平面上点集形状的认定。

霍夫变换以及K-平均算法都对像图中的像进行了假设，并认为物的圆心的像在像图中是具有几何特征的点（分别是圆的圆心和椭圆的中心），从而将问题转化为求像图的几何特征点。

霍夫变换把图形从坐标空间映射到参数空间，在参数空间搜索极值点来确定图形的实际参数，只要图形可以通过较少的参数来描述，该方法就可以在可接受的时间里获得较为精确的结果；而K-平均算法是模式识别中常用的聚类算法，利用该算法我们把图像中五个圆的点分为五类，对每一类求取重心点来获得其圆心所在位置。

K平均算法可以在极短时间内获得精确的圆心坐标，有着良好的性能，但对图像质量要求较高，需要对较为严格的预处理过程。

而霍夫变换抗干扰能力相对要强一些，但是运算时间较长，且只能用于图形近似为圆的场合，对椭圆的场合运算速度则会过慢。

我们用问题一中的三个模型和算法对问题二进行了求解计算。

结果显示，基于K-平均算法的数学模型最合适本题的求解，此时，我们到靶标在像平面内的圆心坐标为：（-187，193）、（-86，187）、（130，171）、（-225，-119）、（72，-119）。

问题三中，我们通过三维空间中的不同的摄像角度，给出了不同的像图。

通过靶标圆心标定求得靶标圆心在像平面内的实际坐标；同时利用问题一中基于K-平均算法的数学模型给出计算坐标。

我们将其做距离差比较，得到两者圆心的像素坐标距离误差为[4.4721;4.1231;0;1.0000;0]像素。

这说明我们的模型准确度比较高；在不同角度都能给出相应的结果，说明模型具有一定的稳定性。

问题四中，我们建立了一个确定两相机相对位置的数学模型。

利用像平面中像点的坐标，通过坐标变换逆过程处理，通过确定物坐标系与像坐标系关系向量，间接确定两个像坐标系的空间位置关系。

论文在最后给出了模型的评价。

本文的特色是将数学模型与计算机算法，计算机图形学相结合，并较灵活的对切线等有关概念在离散的点空间中进行了重新定义，使问题得到较好的解决。

关键词：双目定位系统标定霍夫变换K-平均算法
一问题的提出
1.1背景
数码相机定位在交通监管（电子警察）等方面有广泛的应用。

所谓数码相机定位是指用数码相机摄制物体的相片确定物体表面某些特征点的位置。

最常用的定位方法是双目定位，即用两部相机来定位。

对物体上一个特征点，用两部固定于不同位置的相机摄得物体的像，分别获得该点在两部相机像平面上的坐标。

只要知道两部相机精确的相对位置，就可用几何的方法得到该特征点在固定一部相机的坐标系中的坐标，即确定了特征点的位置。

于是对双目定位，精确地确定两部相机的相对位置就是关键，这一过程称为系统标定。

标定的一种做法是：在一块平板上画若干个点，
同时用这两部相机照相，分别得到这些点在它们像平
面上的像点，利用这两组像点的几何关系就可以得到
这两部相机的相对位置。

然而，无论在物平面或像平
面上我们都无法直接得到没有几何尺寸的“点”。

实
际的做法是在物平面上画若干个圆（称为靶标），它
们的圆心就是几何的点了。

而它们的像一般会变形，
如图1所示，所以必须从靶标上的这些圆的像中把圆
心的像精确地找到，标定就可实现。

图1靶标上圆的像有人设计靶标如下，取1个边长为100mm的正方形，分别以四个顶点（对应为A、C、D、E）为圆心，12mm为半径作圆。

以AC边上距离A点30mm处的B为圆心，12mm 为半径作圆，如图2所示。

用一位置固定的数码相机摄得其像，如图3所示。

图2靶标示意图图3靶标的像
1.2问题
（1）建立数学模型和算法以确定靶标上圆的圆心在该相机像平面的像坐标,这里坐标系原点取在该相机的光学中心，x-y平面平行于像平面；
（2）对由图2、图3分别给出的靶标及其像，计算靶标上圆的圆心在像平面上的像坐标,该相机的像距（即光学中心到像平面的距离）是1577个像素单位(1毫米约
为3.78个像素单位)，相机分辨率为1024×768；
（3）设计一种方法检验你们的模型，并对方法的精度和稳定性进行讨论；
（4）建立用此靶标给出两部固定相机相对位置的数学模型和方法。

二基本假设
1、不考虑光的衍射，色散，只考虑反射和折射；
2、假设所有场景中感兴趣的点都在镜头前面；
3、将相机简化为一个小孔成像机构，对于产生的相差，以及相机对像平面成像的调整
作用不予考虑；
4、像平面是由一个个有大小的像素点组成，是一个不连续的点空间，而几何定义大都
是在连续空间内定义的，这里假定，在几何推理中，像平面是连续的面，也即每个象素点除了表示坐标外，不再具有实际的大小；
5、为了在离散空间中求解，在解答过程中对相切、相割的含义做了一些调整，具体论
述见正文；
6、对于题中随给相机中的几何关系有如下描述：焦点，透镜中心，像平面中点三点共
线，此直线即为相机主轴；透镜与像平面平行；相机主轴垂直于像平面，且垂足为像平面中心。

三符号说明
为了便于描述问题和表达，我们对文中的变量符号进行解释说明，如表1所示：表1符号说明一览表
符号意义
K ，i K 像空间中点的集合（经过简化，是边界点的集合）
r b a ,,像空间中参数
R B A ,,分别为r b a ,,的集合
P ,i P 聚类算法中参数空间的点集
'r 霍夫变换中给定i i b a ,，曲线系方程得到的r 的理论值
J 霍夫变换中'r r −的集合
f 相机的像距
α坐标系间X 轴间的夹角
β坐标系间Y 轴间的夹角
γ坐标系间Z 轴间的夹角
z y D D ,,D x 为像坐标系原点在物坐标系中的坐标
T 坐标系变换算法中平移变换的变换矩阵
321R R R ，，坐标系变换算法中绕X Y Z 三个方向的变换矩阵
R
坐标系变换算法中旋转变换的变换矩阵
四问题分析
本问题涉及的是数码相机的定位问题，问题出现在双目定位的背景下，要解决的问题是如何实现物体表面某些特征点在物平面和像平面之间的坐标转换，其中如何在像平面中找到发生变形的靶标（圆）像的圆心是需要解决的核心问题。

4.1基本思路
题目中给出了已经得到的靶标的像和一些可测量的参数，如相机的像距，分辨率及靶标圆的半径和各圆之间的相对位置等。

精确确定两部相机的相对位置是最终目的；实现物体表面某些特征点在物平面和像平面之间的坐标转换，是解决问题的基础；而像的变形决定了找到圆心在像上的投影是问题的核心也是难点。

我们考虑先找出实现物体表面某些特征点在物平面和像平面之间的坐标转换的方法，在此基础上找到靶标圆心在像平面上的圆心的具体位置，然后通过坐标转换，把空间点从图像中恢复出来，最后根据两次恢复中求出的空间点位置的不一致关系，通过坐标系转换推出两部相机的相对位置。

4.2具体分析
问题1，需要确定靶标上圆的圆心在相机像平面的像坐标。

有两种方法可以考虑实现，第一是利用物与像之间某些不变的相对位置关系找到对应与圆心的像点；另一种方法就是借助计算机图形处理中的霍夫变换和聚类算法在一定模型假设基础上求解。

问题2主要是模型的求解。

根据问题一中得到的模型和算法，具体的编制程序进行实现，具体用到的有++
C，lingo,matlab帮助求解。

问题3是对问题1中的模型以及模型2中的计算结果的检验。

我们考虑通过在三维空间内寻找另外一个摄影角度，数码相机在此角度可以得到图2的另外一个像，同时可以给出在像平面内圆心的靶标的像坐标。

利用问题1的模型和算法求得的圆心坐标与之进行比较，如果圆心坐标可以重合的很好，说明模型的精度比较好；同时在本文中我们还对模型和算法的精度进行了初步的讨论。

问题4是需要找出可以确定两部固定相机相对位置的方法。

我们借助该靶标在两部相机的像平面中所成的像，求出其在像平面中的坐标，并通过逆变换，把空间点从像平面的二维图像中恢复到三维坐标系中，通过两次恢复得到的坐标之间的不一致关系确定出两部相机相对位置。

五模型的建立与求解
5.1模型的准备
5.1.1几个相关坐标系统的说明
由于本问题设计到三维空间和二维像平面，同时有长度坐标和像素的麻烦，因此，在建立模型之前，先要了解几个相关的坐标系统，以及它们之间的关系：
1)世界坐标系统
也称物坐标系统XYZ ，一般的D 3−场景都是用这个坐标系统来表示的，本文中的靶标和像平面之间的关系必须通过此坐标来联系起来；
2)相机坐标系统
又称为观察坐标系统或光心坐标系。

以相机的光心为中心设定的坐标系统xyz ，一般常取相机的光轴为z 轴。

3)图像平面坐标系统
又称像坐标系。

在相机内所形成的图像平面坐标系统),(O '''y x 。

取像平面与相机坐标系统的),(O y x 平面平行，且x 轴与'x 轴，y 轴与'y 轴分别重合，这样像平面原点就在相机光轴上。

5.1.2初期像图的转换处理
在本问题中，像平面是一个二维空间，相应的像素值即为空间内各点的横纵坐标。

在图像平面坐标系统中。

首先，我们通过边缘检测的方法确定靶标在像平面中各个圆像的边界。

即通过检测每个像素和其直接邻域点灰度是否连续，以决定该像素是否确实处于一个物体的边界上。

因此，我们需要将像平面图转化为二维灰度图，当图像中各个像素的灰度发生很大的跳跃产生不连续时，这时出现像的边界，同时记下边界在相图平面中的坐标。

两个具有不同灰度值的相邻区域之间总存在边缘，边缘是灰度值不连续的结果。

边缘检测的基本思想是利用边缘增强算子，突出图像中的局部边缘，然后定义像素的“边缘强度”，通过设置门限的方法提取边缘点集。

具体操作步骤如下：
Step1：把二维像图读入MATLAB的变量中；
Step2：提取靶标图像的边界——确定边缘的方法
Step3：靶标图像的边界是由若干个点组成的，各点坐标可以通过读取图像的过程获得。

常用的边缘检测算子有很多种，我们使用Sobel算子、Robert算子、Prewitt算子、LOG算子、Canny算子和Zerocross算子对图3中的靶标图像进行边缘检测和确定，结果如图4所说。

通过不同算子边缘提取效果比较，sobel算子的效果最好，zerocross效果最差，我们最终确定选用sobel算子进行靶标像图的边缘检测，然后提取出图像的边界。

通过边缘检测和边界提取，我们将图像从二维转化为边界线离散点，实质上是一个降维操作过程，有利于模型的简化，同时减小计算量，提高运行速度。

图4各种算子的靶标的图像边界
5.1.3有关切线的说明
切线定义：
在平面几何中，和圆只有一个公共交点的直线叫做圆的切线。

这种定义不适用于一般的曲线；在图5中，PT是曲线C在点P的切线，但它和曲线C还有另外一个交点；相反，直线l尽管和曲线C只有一个交点，但它却不是曲线C的切线。

另一种对曲线切线的定义如下：
P和Q是曲线C上邻近的两点，P是定点，当Q点沿着曲线C无限地接近P点时，割线PQ的极限位置PT叫做曲线C在点P的切线，P点叫做切点（无限逼近的思想）。

图5曲线的切线示意图
为了能够读取照片数据便于计算机处理，我们利用无限逼近的思想，在实际的像平面中作切线做如下的说明：首先，我们认为，在物平面内存在切线的图形的像在像平面中也具有切线。

若设一个精度系数q（q即为被直线所截点集形成的类劣弧点集中包含点的个数，q为正整数，如下图6），则我们对像平面中的切线做如下定义：
经过像（被证明存在切线的图形的像）点集的直线族中与直线有<=q个交点的直线，即近似认为是该图像在q精度下的切线。

由此可推出，在某一精确度下，像平面中像图的公切线并不一定唯一，而是一组符合要求的直线，但在每次取切线时则根据具体情况唯一的确定一条切线。

精度系数的确定：
在像图中，我们精度系数q的取值如图6所示：
在提取出的像图边界图形上任做一条直线，则直线将图像（点集）分为两部分，点少的一侧的点的个数，即为我们精度系数q的取值。

显然，精度系数k越小，得到的切线越精确，当q=0条件下的切线是切线的极限状态。

图6精度系数示意图
5.2问题1的模型建立和算法设计
问题1的分析
在问题1中需要对像平面和世界坐标系中求图像圆心坐标。

如图7所示，像平面的中心O’是像平面的像素坐标的原点位置，而像距是焦点到像平面距离，焦点在OO’这条线上，同时x-y平面平行于像平面。

因此，问题转化为求在像平面坐标系中靶标圆心的像坐标。

图7空间靶标及其在像平面内坐标关系示意图
5.2.1割线逼近模型和算法设计
割线逼近模型
相机成像应用的是小孔成像的原理,根据基本假设,我们认为相机是一个小孔成像机构。

小孔成像有关的性质如下:
1）在小孔成像中直线的像是直线；
2）物平面上的元素（特征点、圆和直线）都能反映在象平面上；
3）在原像是平面图形的情况下像空间较之物空间不会发生扭曲,即虽然可能发生形变,但点之间的相对位置关系保持不变。

我们的目标是在像平面内像图中存在的固定不变的几何关系中找到稳定关系从而确定圆心在像平面的位置。

由如图8所示，靶标示意图，我们可以得出如下关系：
1、靶标是由五个圆和四条止于相应圆心的线段构成；
2、对于顶点上的圆相邻两条线段交于圆心；
3、做靶标图形的外切正方形（外切正方形如图，其存在性显然），这样存在四条切线，和四个由切线两两相交产生的顶点；
4、像图中存在明显的灰度突变，可以用来确定边界。

在成像的过程中,如图8（b），几何图形的形状会发生改变，但是图形间的对应关系不会发生改变。

因此在分析过程中，我们力图抓住这些关系确定像平面内的“圆心”。

（a）（b）
图8通过切点定圆心示意图
如图8所示，在图8（a）中构造圆A、圆B、圆C的公切线。

由于正方形的性质和圆的直径的判定方法可知A1-E1、A2-C2与圆A直径重合，则根据“两直径交点为圆心”这一性质可知A1-E1、A2-C2交点OA为圆A的圆心。

在图8(b)中，由小孔成像保持图像物平面中各元素的相对关系这一性质，我们可以想象，如果照相时靶标上存在切线A1-B1-C1（A1、B1、C1为切点，以下同理）、C2-D2、D1-E1、A2-E2，由相切的性质，在像平面上必存在相应的A1’-B1’-C1’、C2’-D2’、D1’-E1’和E2’-A2’并与相应圆有且只有一个交点。

在前面的原理中我们还说到直线在小孔成像中像平面的对应图形依然是直线，因此，A1’-B1’-C1’、C2’-D2’、D1’-E1’和E2’-A2’就是靶标上圆在像空间中像的切线。

如图所示分别连接相对的两个切点A1’，E1’和A2’，C2’。

由于两点确定一条直线，我们在像平面中可以得到A1’-E1’和A2’-C2’两条直线。

仍然由于相应图形的相对关系的不变性，对比（a）图，我们可知A1’-E1’和A2’-C2’的交点即对应于原靶标的圆心。

圆心A，C，D，E 对应的像点同理可得。

圆心B点对应的像点坐标：
根据灭点的定义：立体图形各点延伸线向消失延伸的相交点，如果是平行透视，只有一个灭点，在对象中间的后方。

方法是延长物体左右纵深的两条有会聚趋势的线，向后方会聚与一点，平行透视能产生纵深感。

连接A圆心、E圆圆心的延伸线l ae与连接C 圆圆心和D圆圆心的延伸线l cd向消失延伸的相交点为G，G就是平行透视的灭点，有透视学理论和几个关系可知，灭点G点、B的切点B1点与圆B的圆心O2必定共线，于是，又构造出一条过B圆圆心O2的直线，连同直线L ac一起构造直线方程组，求解两直线方程的交点，此交点即为B圆圆心O2。

割线逼近模型的算法设计
Step1：做边界探测，找到各个像点区域的边界点集K ；
Step 2：在区域点集j -i K （j -i K 表示这两个点集间做公切线）中，各取一个点配对
连线，对这样的连线进行穷搜索，知道找到符合精度q 要求的直线作为这
两个点集的公切线；
Step 3：重复2的操作，直到找到全部四条公切线，并记录下各个切点；
Step 4：连接相应的切点，这些切点连线的交点即为所求圆心的像。

5.2.2霍夫变换模型和算法设计
关于霍夫变换
霍夫变换（Hough Transform）是图像处理中对图像进行边缘检测和几何形状识别的基本方法，应用很广泛。

其中，在图像的倾斜角度检测方面有着较为重要的应用。

最基本的霍夫变换是从黑白图像中检测直线(线段)。

它的基本原理是利用函数空间和参数空间的转换，来确定参数向量。

霍夫变换的算法说明
求靶标中的圆心在像平面中对应的点的坐标，这里我们忽略由于摄像机角度（即像平面与物平面成一定角度）而产生的圆的像的形变。

我们认为像平面中的图形近似的为圆形。

根据投影原理，在二维投影中物平面与像平面中点的相对关系不发生改变，即像不发生扭曲。

因此，如果假设像为圆形，那么靶标上圆的圆心也应该对应于其像的圆心。

我们这时只要求出像平面中各个像的圆心即可。

现在，我们利用霍夫变换来求解圆心在像平面的像（以下一段简称为像圆）。

我们依然利用前文方法对像圆进行边界识别。

这样我们得到了像圆边界上的有限点集}{
为像边界的点),(),(y x k y x K ==。

由于在像平面的二维坐标系中圆的方程
总可以表示为：222)(r b y a x =−+−）（我们只要确定参数a,b,r 即可。

(题中所给像
平面图的像素尺寸为实际尺寸的1/3，因此首先把坐标系缩小三倍)。

此时，显然有2810≤≤a ，2810≤≤b ，500≤≤r 。

我们取步长为1个像素，即以1个像素为单位，将a,b,r 在取值范围内离散化，这样就对应有限个在三维空间ABR 中的点。

即有：
{}{}{}
N
r r r N
b b b N a a a ∈≤≤=∈≤≤=∈≤≤=,500R ,2810B ,2810A 从而，得到点集：{}R
r B b A a r b a P P i i i i i i i ∈∈∈==,,),,(现在将前面在边界识别中获得的点带入方程，得到一个ABR 中的曲线系}{K y x B b A a r y b x a C C i i i i i i ∈∈∈=−−+−=),(,,,0)()(:C 2
'22在该曲线系中，取B b A a i i ∈∈,，带入上述曲线系的方程，则可唯一确定'r ，有定
义数量集合
}
{0,,,''>∈∈−==r P r b a R r r r j J i i i ）且（即取最接近于i r 的'r ，则认为i P 为符合条件的点，则计i P 的数量加一。

对所有曲线系C 中的曲线做相同处理，取}{0Z n N i ≥∈=i i i i n n P n 且的个数，为为ABR 坐标系中各点符合条件的次数，取i n max ，此时),,(i i i i r b a P =对应的i a ，i b ，i r 即为所求像圆的参数。

此时，i a ，i b ，i r 即为所求的像平面中对应于圆的参数（i a ，i b ）即为所求像图中对应圆心的点坐标。

5.2.3K-平均算法的数学模型和算法设计
K-平均算法的数学模型
聚类算法和霍夫算法都是利用图形的特点,进行量化确定数量关系的方法。

其中动态聚类法是先行选择若干样本点作为聚类中心,再按某种聚类准则(通常采用最小距离原则)使各样本点向各个中心聚集,从而获得初始分类,再计算其重心,作为新的中心,然后进行第二次分类,一直到所有样本不再调整为止。

动态聚类法有K-平均算法、ISODATA 算法、模拟退火算法等。

本文采用K-平均算法,初始聚类中心点的选择很大程度地影响聚类效果。

K-平均算法以欧式距离作为相似性测度,它是求对应某一初始聚类中心向量
()
T k 1U ,U V �����，……=的最优分类;使得评价指标v J 最小。

2K 1i x i t v i i Z X minJ ∑
∑=Ω∈−=���s.t (1)d j 1,K i 1,UX u LX j ij j ≤≤≤≤≤≤�(2)
{}n tj j X min LX =d j 1≤≤t=1，……，N (3)
{}
tj j X max UX =d j 1≤≤t=1，……，N (4)∑Ω∈Ω=i t X t i
i X 1Z ��,其中i i n =ΩK i 1≤≤，N
t 1≤≤(5)其中K ：预期的聚类中心数目；
i Ω：—类i 的样本集合；
i i n =Ω：为类i 的样本数目；
d ：为样本的属性数目（维数）；
N ：全部样本数目；
i Z �：第i 类的聚类中心向量；
j LX ：所有样本的第j 维下限；
j UX ：所有样本的第j 维上限；
算法设计
首先，我们依然对形变进行假设，这里，我们认为由于相机位置以及光的角度，像比起原像的在形状上产生了一定的形变，并且，形变的结果是靶标上圆在像平面中的像为规则的椭圆形。

比起霍夫变换中把像认为是圆，这一点改进是很重要的，它使情况更接近实际，误差更小。

椭圆依然是对称图形，认为椭圆的中心即对应于物平面中圆的圆心。

对于椭圆（包括圆），其重心与几何中心是重合的，可以利用这一性质来求取像中圆（椭圆）的圆心，对每个圆或椭圆，求得其重心后即得到圆心的位置。

因此首先应该对属于不同的圆上的点进行分类再求重心，避免不同的圆上的点之间的干扰。

在这里我们采用了聚类的方法，采用K-平均聚类算法，对预处理后的图像进行分类处理。

k-平均聚类分析是一种广为人知的方法，它通过尽量缩小一个分类中的项之间的差异，同时尽量拉大分类之间的距离，来分配分类成员身份。

k-平均中的"平均"指的是分类的“中点”，它是任意选定的一个数据点，之后反复优化，直到真正代表该分类中的所有数据点的平均值。

“k ”指的是用于为聚类分析过程设种子的任意数目的点。

k-平均算法计算一个分类中的数据记录之间的欧几里得距离的平方，以及表示分类平均值的矢量，并在和达到最小值时在最后一组k 分类上收敛。

k-平均算法仅仅将每个数据点分配给一个分类，并且不允许成员身份存在不确定性。

分类中的成员身份表示为与中点的距离。

k-平均算法大致流程如下：
每个簇用该簇中对象的平均值来表示。

Step 1：首先将所有对象随机分配到k 个非空的簇中。

Step 2：计算每个簇的平均值，并用该平均值代表相应的簇。

Step 3：根据每个对象与各个簇中心的距离，分配给最近的簇。

Step 4：
然后转第二步，重新计算每个簇的平均值。

这个过程不断重复直到满足某
个准则函数才停止。

为了实际的计算方法，我们对问题作如下合理假设：
1由于我们已经将像点集的图形近似认为是椭圆，根据椭圆的性质，椭圆的中心就是它的重心，因此求椭圆中心实际上转换为求重心。

2根据题中所给照片，我们发现像平面上点的分布是除边界外处处稠密的，因此可以认为题中像平面上的像点集都是表面均匀的图形。

而均匀图形的边界图形与其所围成的平面图形具有相同的性质。

为了减少算法中的计算量，我们这里将求椭圆面重心转化为求椭圆形的重心。

3图像中不存在相交的圆。