“十字线”交叉点定位法求相交线焦点位置.docx

2014南京理工大学大学生数学建模竞赛

承 诺 书

我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括电话、电子邮件、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。如有违反竞赛规则的行为，我们将受到严肃处理。
我们授权全国大学生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。

我们参赛选择的题号是（从A/B中选择一项填写）： A
我们的参赛报名号为（报名编号）： 84
所属学院（请填写完整的全名）： 理学院
参赛队员 (打印并签名) ：1. 徐 昊 18362963697
2. 李业良 13270801839
3. 方昱暐 18362963699

日期： 2014 年 8 月 14 日

评阅编号（由组委会评阅前进行编号）：

13

“十字线”交叉点定位法求相交线焦点位置
摘要
本文是一个有关求“十字线”相交线焦点位置的解决方案。我们针对相交线的不同情况提出了不同的解决方案，通过查阅资料，综合运用多项式函数拟合，分析高斯模糊的原理及其逆的形式等方法建立起数学模型成功准确地确定“十字线”的交叉点位置，并给出了相应的置信水平。
针对问题一，我们重点考虑了曲线相交的情形。我们以一个像素作为基本单位建立起了平面直角坐标系。在这个坐标系框架内分别对横竖两条曲线取得二十个参考点。依据这些参考点，我们通过多项式函数拟合的方式用多项式函数近似替代了图中的曲线。然后联立两个多项式函数并求解，得到的即是要求的交点坐标。求解过程中我们为了避免在多项式次数过高导致的产生振荡现象，影响精度，没有采用直接插值模拟函数图像，而是选择大量参考点并使用低次方程进行拟合以提高精度，防止出现较大的震荡。最后我们给出了交点的置信水平与误差，证明该方法准确可靠。
针对问题二，我们认为在亚像素条件下的图形只有线条宽度上的改变，其他条件与问题一的条件无异。此时我们选取亚像素线条中间的点作为参考点，于是本问题即回到问题一。通过建立相应的平面直角坐标系，使用多项式函数对曲线进行函数模拟，再将两个函数联立找出近似数值解的坐标即为交点坐标。
针对问题三，我们的重点放在了寻找高斯模糊算法的逆的具体形式，通过对图像进行模糊算法的逆运算来恢复到亚像素级别的情形，进而退化到问题一从而得以解决问题。我们将素材图片转化为显示每个像素点灰度值的矩阵，将高斯模糊看作是对灰度值矩阵元素的线性操作，进而将高斯模糊算法的逆看成是一个多元非齐次线性方程组问题。求解这个方程组我们即可得到模糊操作之前的灰度值矩阵，也就得到了原图片。最后我们选取一张素材图作为例子，运用高斯模糊算法的逆算法消除原图像的模糊效果，成功将其还原到了亚像素条件下的曲线相交问题，并求解出交点坐标。
最后，本模型将实际操作图片的每一步细致呈现在论文中，并附上相应的程序源代码。从最终结论来看，模型较好解决了“十字线”相交线焦点位置问题，模型结论与和实际符合地较好。

关键词：“十字线”交叉点定位 函数拟合 高斯模糊 灰度值矩阵

一、 问题重述
在实际工作中，为了确定目标的位置，通常采用“十字线”交叉点定位法确定目标的距离和方向，它的关键问题为：如何准确地确定“十字线”的交叉点位置？解决这个问题的难点是因为有以下2种情况的出现：
(1)从目标返回的图像，由原来的水平和竖直的直线变成了不水平或者不垂直的状态，甚至变成弯曲状态。
(2)由于噪声、图像设备缺陷、光线影响等很多因素导致实际图像存在模糊、光强不均匀、线条变形等等。
针对上述思考，本文着重解决和回答了以下三个问题：
(1)建立“十字线”交叉点求解模型，确定水平、倾斜、弯曲这三种类型的直线或者曲线交叉点位置坐标。
(2)建立亚像素级的“十字线”交叉点定位模型，确定水平、倾斜、弯曲这三种类型的直线或者曲线交叉点位置。
(3)针对实际捕捉到的图像建立相应的“十字线”交叉点位置求解模型，并求得具体交点坐标位置。

二、 模型假设
（1） 我们得到的图像素材真实可靠；
（2） 我们得到的是以单个像素为基本单位的位图；
（3） 单像素不考虑RGB值，只有灰度值，灰度值取值从0到255闭区间内的整数；
（4） 亚像素情况下认为图像只有线条宽度上的改变，各像素点灰度值取值依然从0到255闭区间内的整数；
（5） 实际捕捉到的图像的对焦模糊看作是高斯模糊；
（6） 高斯模糊的模糊半径远小于图像的长或宽，模糊权重按照距离被操作像素点的远近而正态分布。

三、 模型符号说明
f(x)：横向曲线拟合函数
g(x)：纵向曲线拟合函数
A：f(x)和g(x)的交点
B：分析图上的交点
d：交点坐标绝对误差
P：原像素灰度值矩阵
P'：模糊处理后的像素灰度值矩阵
G：高斯模糊映射
G-1：高斯模糊逆映射
r：模糊范围
w：权重矩阵，是一个(2r+1)*(2r+1)矩阵，且矩阵所有元素非负，其和为1
D：线性方程系数矩阵

四、 模型建立与求解
4.1问题一
4.1.1问题分析
问题一是有关在单像素线条情况下的交点坐标求解问题。实际情形中，我们“十字线”交叉点定位法求相交线焦点位置时会遇到图像中的线条出现倾斜、弯曲等变形的情况，我们需要在各种变形情况中找到交点的准确坐标。考虑到一般情况是弯曲变形，其他各种情形都可以看做是弯曲变形的特例，我们在这里着重考虑弯曲变形情形下的交点坐标求解问题。
我们通过建立相应的平面直角坐标系，使用多项式函数对曲线进行函数模拟，再将两个函数联立找出近似数值解的坐标即为交点坐标。
4.1.2模型建立与求解
以原题中的素材为例，如下图所示：

图一：素材全图

我们以图一中右上角的交点为例，试图找出其交点的具体坐标值。首先我们以一个像素为最小单位，取于左上角水平方向距离475个像素，铅直方向距离200像素的点为原点，水平向右为x轴正向，铅直向上为y轴正向，建立整个图片的直角坐标系。如图二所示：

图二：素材图建立平面直角坐标系

由于曲线清晰，我们认为构成此曲线的点中，每一个点只占一个像素。从原点起每隔25个像素取曲线上一点，分别在横竖两条曲线上取二十个参考点，得到以下坐标数据：

表一：f（x）与g（y）上部分参考点坐标

单位：像素
由以上坐标数据我们可以用多项式函数对两条曲线分别进行拟合。在拟合的过程中为了避免多项式拟合点数多的情况下，次数过高导致的产生振荡现象，影响精度，我们选取了五次函数进行拟合，使用Matlab 软件得到以下两个拟合方程：
①、对设第一条曲线拟合到5次幂：

拟合曲线f(x)给出了大致趋势，结果很好，通过运行相关代码，在置信水平为95%的情况下,我们得到了各参数的估计值以及其置信区间（括号内即为置信区间）。
a5 = 3.849e-11, (2.786e-11, 4.911e-11)
a4 = -8.924e-9, (-1.114e-8, -6.711e-9)
a3 = -3.242e-6, (-3.726e-6, -2.758e-6)
a2 = 0.001764, (0.001692, 0.001836)
a1 = 0.1686, (0.1623, 0.1748)
a0 = 0.5527, (0.1228, 0.9827)
由此得到拟合曲线f(x)图像：

图三：拟合曲线f(x)图像

②、设第二条曲线g(x)，其为一多值函数，即
g (x) =
, x ∈ [-71, -56]
, x ∈ [-56, 43]

对于，拟合到5次幂，即：

拟合曲线给出了大致趋势，结果很好，通过运行相关代码，在置信水平为95%的情况下,我们得到了各参数的估计值以及其置信区间（括号内即为置信区间）。
b5 = -3.305e-007, (-0.02764, 0.02764)
b4 = -7.829e-005, (-8.85, 8.85)
b3 = -0.008266, (-1131, 1131)
b2 = -0.336, (-7.22e+004, 7.22e+004)
b1 = 0.1905, (-2.3e+006, 2.3e+006)
b0 = 0.9706, (-2.924e+007, 2.924e+007)
由此得到拟合曲线g1(x)图像：

图四：拟合曲线g1(x)图像

对于 ，拟合到五次项，即：

拟合曲线给出了大致趋势，结果很好，通过运行相关代码，在置信水平为95%的情况下,我们得到了各参数的估计值以及其置信区间（括号内即为置信区间）。
c5 = 8.326e-008, (3.51e-008, 1.314e-007)
c4 = 5.225e-006, (1.49e-006, 8.96e-006)
c3 = 3.623e-005, (-9.182e-005, 0.0001643)
c2 = -0.002043, (-0.01027, 0.006188)
c1 = 1.317, (1.176, 1.459)
c0 = 0.0275, (-3.51, 3.565)
由此得到拟合曲线g2(x)图像：

图五：拟合曲线g2(x)图像

通过以上的拟合方程我们可以得到如下交点附近的函数拟合图像：

图六：交点附近的函数拟合图像

图像显示，交点在（0,0）之间，所以就是解方程在（0,0）附近的根。通过运行相关的代码可知，在90.25%的置信水平下，交点为A(0.4580,0.6303)由此我们可以认为两条曲线的交点在我们所建立的平面直角坐标系体系下的坐标为(0.4580,0.6303)。

4.1.3模型检验
我们在建立坐标系时，通过肉眼对素材的直接操作，将坐标原点选取在了离真实交点差距极小的地方，故可以认为以上方法得到的交点坐标和原点的距离即为该算法的绝对误差。由此我们可以计算出其绝对误差：
d=|AO|=0.7794<1（像素）
可见使用函数拟合的方式可以在较高的置信度下得到相对精确的结果。而其他各种情形都可以看做是弯曲变形的特例，均可用函数拟合的方法求得交点处的坐标值。

4.2问题二
4.2.1问题分析
经查阅相关资料我们知道，亚像素的形成原因相对比较复杂，若单独处理其图形则十分麻烦，于是我们可以考虑其和问题一只是像素取值上面的区别。在亚像素条件下的图形我们认为只有线条宽度上的改变，其他条件与问题一的条件无异。这里我们不妨假设线条宽度由问题一的一个像素增加到了三个像素，考虑到一般情况是弯曲变形，其他各种情形都可以看做是弯曲变形的特例，我们在这里着重考虑弯曲变形情形下的交点坐标求解问题。
此时我们选取亚像素线条中间的点作为参考点，于是本问题即回到问题一。通过建立相应的平面直角坐标系，使用多项式函数对曲线进行函数模拟，再将两个函数联立找出近似数值解的坐标即为交点坐标。
4.2.2模型建立与求解
同问题一，我们选取原题所提供的素材四曲线图中右上角的交点为例，建立相应的平面直角坐标系，如下图所示：

图七：亚像素素材图建立平面直角坐标系

类似的，我们在亚像素图形的中心选择参考点，该问题即退化为问题一。模仿问题一，我们通过相应的函数模拟得到拟合图像，再求得拟合函数交点坐标，即为实际交点坐标。
当线条宽度由问题一的一个像素增加到了三个像素时，亚像素素材图右上角的交点坐标为(1.374,1.891)。

4.3问题三
4.3.1问题分析
由于噪声、图像设备缺陷、光线影响等很多因素导致实际图像存在模糊，导致最终成像上的像素点有了灰度值上的差异。一般情况下我们认为像素点的灰度值取值介于从0到255闭区间内的整数，灰度值为零表示全白，灰度值为255表示全黑。由于最终成像上的图形总是呈现出其几何中心的像素灰度值最大，四周像素灰度值逐渐减小，故我们认为最终成像接近于经过高斯模糊的图片。
本题的解答着重于寻找高斯模糊算法的逆，通过对图像进行模糊算法的逆运算来恢复到亚像素级别的情形，进而退化到问题一从而得以解决问题。
4.3.2模型建立与求解
由于我们得到的是以单个像素为基本单位的位图，且单像素不考虑RGB值，只有灰度值，灰度值取值从0到255闭区间内的整数。故我们可以将一张分辨率为a*b的图像等价于一个a*b的像素灰度值矩阵P，矩阵内的每个元的值是图像相应位置像素的灰度值。
我们定义一个关于矩阵P的映射G，G的作用是在给定模糊范围r和2r+1)*(2r+1)权重矩阵w的情况下将原矩阵P映射为模糊矩阵，即有：
G(P,r,w)=P' （1）
在问题三中，我们已知经模糊处理后的像素灰度值矩阵P'，选取适当的r和w之后，（1）式可以看成是一个a*b元非齐次线性方程组，其系数矩阵D与r、w有关，且|D|≠0。则（1）式存在唯一解P，即G为一一映射。
此时存在G的逆映射G-1，使得：
G-1(P',r,w)=P （2）
（2）式表明：已知经过模糊处理后的图像，若选取适当的模糊半径r和权重矩阵w，经过高斯模糊的逆映射可以得到原图像。再将灰度值高于127的像素替换成255灰度值，不足127（含）的像素替换成0灰度值，即恢复到亚像素级别的情形，进而退化到问题一从而得以解决问题。
下面给出实例：
我们选取本题素材中的一张图片如下：

图八：实际中采集的图像样张

考虑上方的交点，我们将其周边150*150像素邻域单独剪裁出来作为分析图，取灰度，如下图所示：

图九：交点邻域（分析图）

此时若取r=5，w取均值为0方差为1的正态分布作为权重分布，经过模糊逆处理及黑白化处理之后的图像如下图所示：

图十：处理后的分析图

再将边缘平滑化，图像即可以被退化到亚像素级别情况（问题二），如下图所示：

图十一：退化为亚像素级别焦点问题的分析图

自此，素材中采集自实际照片的焦点问题已经退化至问题二。经过相关计算可知：若以距离分析图左边缘像素数为横坐标，距离分析图下边缘像素数为纵坐标，可以求得这里的交点坐标为B(80.0,64.7)。

五、 模型思考与改进
5.1多项式插值的震荡现象
在解决问题一和问题二的过程中，我们在函数模拟时为了避免在多项式次数过高导致的产生振荡现象，影响精度，没有采用直接插值模拟函数图像，而是选择大量参考点并使用低次方程进行拟合以提高精度，防止出现较大的震荡。

5.2高斯模糊图像的复原处理
对于图像处理，本题中主要涉及的是消除图像的模糊。事实上我们还可以消除图像中出现的噪声。通过将图像中不平整的像素点用其周边像素点的平均效果替代可以实现一定程度的降噪效果。综合使用去模糊和降噪两种方法会使得本问题的解决更加精确。
六、 参考文献
[01]吴宪君；高斯模糊算法的改进及图像处理应用[D],广东，广东石油化工学院，2013,129-131
[02]咸兆勇；图像模糊检测与模糊区域分割研究[D]，广西，广西大学，2013
[03]赵静，但琦主编，严尚安，杨秀文副主编，数学建模与数学实验[M]，北京：高等教育出版社；海德堡：施普林格出版社，2000。
[04]张德荣，王新民，高安民，计算方法与算法语言[M]，北京：高等教育出版社，1981。
[05]易大义，陈道琦，数值分析引论[M]，杭州：浙江大学出版社，1998。

七、 附录
此代码均在matlab下运行。

% 拟合曲线f(x)
x=[ fliplr(-10:-10:-250) 10:10:180 ];
x=x';
y=[ 48 46 44 42 39 37 33 30 26 23 19 16 12 9 5 2 1 -1 -3 -2 -4 -2 -2 -3 -1 3 5 7 10 13 16 20 24 27 31 35 39 44 49 53 58 62 67 ];
y=y';
f1=fittype('a5*x^5+a4*x^4+a3*x^3+a2*x^2+a1*x+a0','independent','x','coefficients',{'a5','a4','a3','a2','a1','a0'});
f=fit(x,y,f1) %显示拟合曲线f(x)
x1=[ fliplr(-10:-10:-250) 10:10:180 ];
y1=f(x1);
plot(x1,y1,x,y,'r*') %得到拟合曲线f(x)与原始数据的对照图

% 拟合曲线g1(x)

m=[-56 -60 -64 -66 -68 -70 -71];
m=m';
n=[ -200 -190 -180 -170 -160 -150 140];
n=n';
q1=fittype('b5*x^5+b4*x^4+b3*x^3+b2*x^2+b1*x+b0','independent','x','coefficients',{'b5','b4','b3','b2','b1','b0'});
g1=fit(m,n,q1) %显示拟合曲线g1(x)
m1=[-56 -60 -64 -66 -68 -70 -71];
n1=g1(m1);
plot(m1,n1,m,n,'r*') %得到拟合曲线g1(x)与原始数据的对照图

% 拟合曲线g2(x)

z=[ -71 -70 -69 -67 -63 -59 -55 -50 -44 -37 -30 -23 -15 -7 8 15 21 26 32 36 40 42 43];
z=z';
t=[ -140 -130 -120 -110 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 10 20 30 40 50 60 70 80 90];
t=t';
q2=fittype('c5*x^5+c4*x^4+c3*x^3+c2*x^2+c1*x+c0','independent','x','coefficients',{'c5','c4','c3','c2','c1','c0'});
g2=fit(z,t,q2) %显示拟合曲线g2(x)
z1=[ -71 -70 -69 -67 -63 -59 -55 -50 -44 -37 -30 -23 -15 -7 8 15 21 26 32 36 40 42 43];
t1=g2(z1);
plot(z1,t1,z,t,'r*') %得到拟合曲线g2(x)与原始数据的对照图
%求出交点坐标。
%建立M文件F.m,源代码见附件F.m。

% 调用fzero函数求根。
h=fzero('F',0) %求f(x)=g2(x)在（0,0）附近的根。
d=f(h);
D=[h,d] %显示交点坐标。

%求拟合后原始数据和曲线的整体比较图。
plot(x,y,x1,y1,'r*',m,n,m1,n1,'r*',z,t,z1,t1,'r*')