2013全国数学建模竞赛 ——碎纸片拼接复原

2013高教社杯全国大学生数学建模竞赛

重庆工商大学姜木北小组作品

编号专用页

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

赛区评阅记录（可供赛区评阅时使用）：

评

阅

人

评

分

备

注

全国统一编号（由赛区组委会送交全国前编号）：

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

碎纸片的拼接复原

摘要

目前，“碎片拼接复原”技术在司法物证复原、历史文物修复及社会生活各项领域扮演着重要角色，对于碎片数量特别巨大而人工又难以在短时间内完成碎片拼接时，要找到一种高效快捷的自动拼接方法已变得尤为重要。本文针对只有中英文的碎片拼接问题，综合分析了从单一的纵切到纵横切以及纵横切双面碎片这三个不同的情况，提出了碎片拼接复原的解决方案.

在问题一中，对于仅有“纵切”且数量相对较少的碎纸片，我们基于边缘去噪和采用构建碎纸图片的左右边缘二值矩阵提取相似度分析的方法，再通过两张图片左右相似度匹配排序，得到附件1和附件2中的碎纸排序（见表2和表3），并运用Matlab的图像处理工具箱，按排列顺序导入碎纸片得到相应拼接结果（见附录附件一）.

在问题二中，由于碎纸片数量相对较多，同时存在横切和纵切的情况，在问题一的基础上增加了碎纸片的上下边缘相似度匹配。在进行人工干预，找到第一张起始碎纸片作为匹配起点后，我们基于索贝尔算子的原理，对碎纸片灰度值进行边缘相似度的旋转检测和比较匹配，最后进行二叉树搜索排序（见表4和表5）。对附件3和4的碎纸图片拼接出的结果详见附录中的附件二.

在问题三中，由于碎纸片是两面的并且碎纸片数量更多，若采用第二问的求解方案则加大了求解难度同时也存在较大误差。因此，我们基于蚁群算法（ACA）的SIFT特征点匹配原理来求解。先提取碎纸图片特征点，然后基于蚁群算法的最优化快速比对匹配，最后基于ACA的搜索排序对碎纸片拼接。Matlab编程所求得的排序结果详见表6和表7，附件5中的碎纸片拼接复原结果见附录中的附件三.

在问题的解决中，我们得出结论:碎纸图片导入量越小，图片匹配出的效果越佳，在相似度的匹配上，出现的误差减小,最后拼出的图像效果好，人工干预量也相对小。本建模考虑到了图片噪声对图片拼接时的影响，选择了去噪效果较好的邻域平均法对图片进行处理。但是，为了解题方便，我们忽略了碎纸机切纸时碎纸片可能产生的边界遗失破损。在这里我们的改进是，对于边界遗失图像碎片的修复，根据复原后的整体形状，可以根据线连续性来拟合此类线段,从而得到较为完整的图像.

关键词：碎纸拼接、腐蚀算法、蚁群算法、图像特征匹配、邻域平均法

一、问题的重述

破碎文件的拼接在司法物证复原、历史文献修复以及军事情报获取等领域都有着重要的应用。传统上，拼接复原工作需由人工完成，准确率较高，但效率很低。特别是当碎片数量巨大，人工拼接很难在短时间内完成任务。随着计算机技术的发展，人们试图开发碎纸片的自动拼接技术，以提高拼接复原效率。现解决如下问题：

1. 对于给定的来自同一页印刷文字文件的碎纸机破碎纸片（仅纵切），建立碎纸片拼接复原模型和算法，并针对附件1、附件2给出的中、英文各一页文件的碎片数据进行拼接复原。如果复原过程需要人工干预，请写出干预方式及干预的时间节点。复原结果以图片形式及表格形式表达.

2. 对于碎纸机既纵切又横切的情形，请设计碎纸片拼接复原模型和算法，并针对附件3、附件4给出的中、英文各一页文件的碎片数据进行拼接复原。如果复原过程需要人工干预，请写出干预方式及干预的时间节点。复原结果表达要求同上.

3. 上述所给碎片数据均为单面打印文件，从现实情形出发，还可能有双面打印文件的碎纸片拼接复原问题需要解决。附件5给出的是一页英文印刷文字双面打印文件的碎片数据。请尝试设计相应的碎纸片拼接复原模型与算法，并就附件5的碎片数据给出拼接复原结果，结果表达要求同上.

二、问题的分析

本题是一个关于碎纸片拼接复原的问题，针对题目的要求，我们对于“碎纸片拼接复原”等一系列问题的处理特点和处理方案，就每个问题做除了以下的具体分析：2．1针对问题1的分析

结合给出的碎纸机碎片情况，分别是一页中文和英文文件的碎片数据，也就是说，在每张碎片中都会出现相当多的完整或不完整的汉字（英文）或者结构，要解决把纵切的碎纸片拼接复原的问题：

首先，利用MATLAB编程建立图片的二值矩阵。鉴于纸片的大小相同，并且只有黑白两色，通过观察附件一的碎片，每张碎片最左端或者最右端不是全空，这些碎片就可以全部转化为二值图，0代表白色，1代表黑色。其次，采用邻域平均法对图像进行边缘去噪，进一步利用分割最佳阈值的迭带算法来达到消除或减少噪声影响的效果。最后，比较图像左右边缘的相似情况,从而得出两个边缘图像的相似度,以此来度量两个图像

能否匹配。复原过程中需要人工干预，将依顺序显示的未拼接碎片，采用半自动拼接，即人工选择碎片并拼接到计算机屏幕上。将序号排列好后；再通过基于MATLAB编程将19张图片分别排好序通过调用imshow()函数合成图片，这样就能得到最终求解结果. 2．2针对问题2的分析

求解要求中我们可以知道该问题在第一问的基础上增加了“横切”碎片，出现的纸张既有横切又有纵切。解决方案与第一问有所不同；增加了“横切”过后，在第一问的左右相似度比较之上增加了上下边缘的相似度比较和匹配。但考虑到求解精度我们采用灰度值（0~255）求解的图片处理方法；同样基于MATLAB编程对附件3和附件4中的碎片分别进行灰度处理得到图片矩阵的灰度值；并采用基于Sobel 的算法对个碎片进行四个边的边缘化处理，提取边缘特征，并对其进行边缘相似度的计算。同时对边缘化处理后的灰度值进行聚类分析确定一个阀值，也就是确定一个误差值。只考虑切缝处并对其相似度一致的缝合。又由于该问中需要拼接的图片数量众多，我们采用基于二叉树的搜索排序算法对相似度矩阵排序并得到该相似度矩阵所对应的图片序列号数组，该排列的

序列号数组即为碎片的组合方式；最后采用MATLAB 循环调用序列号数组对应的碎片名拼接出所求结果达到解决该问题的目的. 2.3针对问题3的分析

对于问题3其实是在问题二的基础上在增加了正反面的拼接。另外通过分析题目中告诉了的类容可以发现附件中同一个序号表示是一张碎纸图片，a 和b 表示同一张碎片的了两面，但并不确定a 和b 的正反。这样同第二问中的解决思路相似，但由于问题三中对应的附件5的碎图片数据量过大如若采用第二问的求解方式也并不理想所以我们采用基于蚁群算法（ACA ）的SIFT 特征点图像碎片拼接方法。针对附件五中的图片基于MATLAB 采用所设计的算法提取其特征点并建立特征点匹配矩阵。再利用这些特征点匹配矩阵根据理论建模设计基于MATLAB 的算法求解出顺序排列的部分序列号和误差节点。将其中求得的序列号对应的碎纸图片导入MATLAB 中求解出拼接图像，最后通过人工干预的方式补全各个残缺节点的碎片得到完整的纸片图像已达到拼接出碎纸片的正反两面的目的.

三、模型假设

1、 假设未碎纸张的文字行方向沿水平方向，字与字之间有间隔且字宽度与高度比值

1/3；

2、假设碎片模型为理想模型，碎片表面光滑平整无磨损且厚度为零；

3、假设在切割过程中除边界外，其他切割线都切割文字；

4、假设每张碎纸片的大小一致；

5、假设碎纸机切的纸片无损坏；

6、假设碎图片中的文字符字体格式相同，英文字符字体格式相同.

四、符号的定义

(,)g x y 图像平面上的一系列平均噪声点，,0,1,2,,1x y N =???-

i X 第i 个二值矩阵 i g

图像的分割阀值，即灰度值

k T

阀值，0,1,2,,k n =??? i Z

灰度均值，1,2,,i n =???

bw

待测边缘图像 (),x i j

边缘点

()x y G G

经横（纵向）边缘检测的图像

Θ

梯度方向

,i k C

第i 像素和它的K 连通像素

()EdSim x

边缘相似度

j i α

不同中心像素点的协议变量，,1,2,,i j n =???

α

为启发式因子的相对重要程度 β

为信息素的相对重要程度 ρ

为信息素蒸发系数 1ρ-

表示信息素的持久性系数

τ

为窗口信息素含量

五、模型的建立及求解

5.1 问题一模型的建立及求解 算法的流程图：

图1

建立图片的二值矩阵

采用邻域平均法进行边缘去噪

比较图像左右边缘的相似情况

运用腐蚀算法轮廓提取方法来消除图像的左右边界的作用

载入图像

5.1.1、把图像转化为二值矩阵

二值图像是指每个像素不是黑就是白，其灰度值没有中间过渡的图像。二值图像一般用来描述文字或者图形，其优点是占用空间少，人们经常用黑白、B&W 、单色图像表示二值图像，但是也可以用来表示每个像素只有一个采样值的任何图像，例如灰度图像等.

因此针对问题一，我们先将对应的碎纸图片用Matlab 调用im2bw 函数对导入的碎纸图片求解其二值矩阵.

5.1.2、图像边缘去噪——邻域平均法

邻域平均法是一种利用Box 模板对图像进行模板操作（卷积操作）的图像平滑方法最简单的平滑滤波是将原图中一个像素的灰度值和它周围邻近像素的灰度值相加，然后将求得的平均值作为新图中该像素的灰度值。它采用模板计算的思想，模板操作实现了一种邻域运算，即某个像素点的结果不仅与本像素灰度有关，而且与其邻域点的像素值有关.

将附件1或2中已经进行二值化的碎纸片矩阵惊醒边缘去噪处理，消弱噪声点起到平滑作用，减小相似度误差，增大匹配准确度.

设),(j i f 为给定的含有噪声的图像，经过邻域平均处理后的图像为),(y x g ，则邻域平均法也可以用数学公式表达：

(1)

式中：,0,1,2,,1x y N =???-，；s 是以(,)x y 为中心的邻域的集合，M 是S 内的点数.

邻域平均法的思想是通过一点和邻域内的像素点求平均来去除突变的像素点，从而滤掉一定的噪声，其主要在实际应用中，也可以根据不同的需要选择使用不同的模板尺寸，如3×3、5×5、7×7、9×9等。邻域平均处理方法是以图像模糊为代价来减小噪声的，且模板尺寸越大，噪声减小的效果越显著。如果

),(j i f 是噪声点，其邻近像素灰

度与之相差很大，采用邻域平均法就是用邻近像素的平均值来代替它，这样能明显消弱噪声点，使邻域中灰度接近均匀，起到平滑灰度的作用. 从理论上看:

a 、要改变均值滤波的滤波能力，就要改变模板的大小，模板越大去噪能力越强，反之越弱；

b 、由于是在整幅图上以取局部灰度平均的方法来淡化噪声，所以不可避免的会使原有图像的清晰度、对比度下降，会损失原有图像的细节，模板越大损失的越多；

c 、对整幅图像用同一个模板，不能适应局部具体情况，会出现有滤波不均匀的情况.

5.1.3、分割阀值的迭代分析法

阀值分割方法是把图像的灰度分成不同的等级，然后用设置灰度门限的方法确定欲分割物体的边界。当用阀值来分割口标与背景时，如果某一灰度值g 是某图像的分割阀值，即小于g 的灰度点将构成口标(不妨如此假设)，而大于g 的灰度点就构成背景一般而言，如果灰度值g 可以作为图像的一个阀值，那么它应该使按这个阀值划分口标和背景的错误分割的图像像素点数为最小.

∑∈=s

j i j i f M y x g ),(),(1

),(

如果前景物体的内部具有均匀一致的灰度值，并分布在另一个灰度值的均匀背景上，那么图像的灰度直方图应具有明显的双峰。可是在许多情况下，噪声的干扰使峰谷的位置难以判定或者结果不稳定本文采用迭代算法，有效地消除或减少噪声对灰度门限值g 的影响.

设有一幅混入噪声的图像g (x),力是由原始图像f (x,y)和e(x,y)叠加而成的。即:

()()()

,,,g x y f x y e x y =+ （2）

这里假设各点的噪声是互不相关的，且具有零均值，标准差为ε。通过阀值分割将图像分割为两部分，由于噪声是随机作用于图像的像素点上，则可以认为在分割出口标

1

g 和背景

2

g 图像上噪声干扰仍为()

,e x y ，即:

()()()11,,,g x y f x y e x y =+ （3）

()()()22,,,g x y f x y e x y =+ （4）

在迭代算法中，需要对分割出的图像分别求其灰度均值，则:

(){}()(){}(){}111,,,,E g x y E f x y e x y E f x y =+= （5） (){}()(){}(){}222,,,,E g x y E f x y e x y E f x y =+= （6）

上式说明，随着迭代次数的增加，平均灰度值将趋向于真值。因此，用迭代算法求得的最佳阀值不受噪声干扰的影响根据上述分析，针对该问题的迭代算法描述如下.

（1）首先选择一个近似阀值作为估训一值的初始值，然后进行分割，产生子图像，并根据子图像的特性来选择新的阀值，再用新的阀值分割图像，经过几次循环，使错误分割的图像像素点降到最少。这样做的效果好于用初始阀值直接分割图像的效果，阀值的改进策略是迭代算法的关键。算法步骤如下:

选择一个初始阀值的估算值{}

00k T T k ==

min max

2

Z Z T +=

（7） 式中，min Z ，max Z 分别表示图像中的最小和最大灰度值. （2）利用阀值k T 把图像分割成两组，1R 和2R ，其中：

()(){}1,,k R f x y f x y T =

≥ （8）

()(){

}

2,0,k R f x y f x y T =<< （9）

(3) 计算区域1R 和2R 的灰度值均值1Z ，2Z ，其中：

()()

()()

(),1,,,,k

k

f

i j T

f

i j T

f

i j N i j Z N i j <>?=

∑

∑

（10）

()()

()()

(),2

,,,,k

k

f

i j T

f

i j T

f

i j N i j Z N

i j >>?=

∑

∑

（11）

式中，f (x,y)是图像上(i ，j)点的灰度值，N(i,j)是(i,j)点的权重系数，一般N(i,j)=1.0。 （4）选择新的阈值1

k T

+

1

12

2

k Z Z T ++=

（12） （5）如果1

k k T T

+=，则结束,否则1k k =+,转步骤(2)。

经过图像分割处理后,有效减少了噪声对图像的干扰范围,并使图像的边缘邻域像素点离噪声干扰的敏感区,从而提高了图像边缘检测抗干扰能力. 5.1.4、腐蚀算法轮廓提取方法

经过上述三求解阶段过后，图像轮廓提取算法就变得非常简单。基于二值图像轮廓提取的算法就是掏空内部点,如果原图中有一点为黑,且它的8个相领点都是黑色时,判定该点是原图像的内部点,则将该点删除。最后,经过这样的算法处缘检测[4, 5]。在数学形态学运算中,腐蚀具有消除物体边界点的作用。结构元素取3×3的黑点块,腐蚀将使物体的边界沿周边减少一个像素。那么边缘检测实际上相当于用3×3块的9个点结构元素对原图进行腐蚀,再用原图像减去腐蚀的图像。令X 为图像,B 为结构元素,BZ 表示结构元素B 平移Z 后的结果, s B 代表结构元素关于原点的对称集合。其数学表示如下:

{};s

B b b B =-∈ （13）

则腐蚀的运算定义:

{};s

z X B Z B X Θ=? (14)

式中, :33E ?的结构元素; ED(X):代表图像X 的边界.

该方法检测到的物体边缘宽度仅为一个像素,因此具有较高的定位精度。同时,二值化处理后的图像具有完整的轮廓,所以轮廓提取检测到的边缘具有连续性. 5.1.5、边缘相似度

边缘相似度的基本思想是用标准的无噪图像产生出标准的边缘图像，设为bw ，将含噪图像用噪声抑制算法滤波处理后产生待测边缘图像，设为bw ,比较两个边缘图像上各个边缘点的相似情况，从而得出两个边缘图像的相似度，以此来度量滤波算法的边缘

保持能力.

对于bw 上的一个边缘点(),x i j ，其边缘相似度定义为：

()()()'''max ,'','''

,min 1

det int i i j j x BW e

x i j BW x x x x EdSim x x erroneous ection po s ---∈?∈-=-?=?

?-∈? （15）

由上式可见，单点边缘相似度反映的是检测边缘点与对应的原标准边缘点之间的位置偏差:如果完全重合，则为最大相似度1;如果相差k 个像素点，则其下降为k e -，如果是一个误检点包括错检点和漏检点，即原来不是边缘点，但检测成了一个边缘点；或

者原来有边缘，但是没检测出来，则为-1.

整幅图像的边缘相似度则定义为所有边缘点相似度的平均公式：

()1x bw

EdSim EdSim x M

∈=

∑ （16）

式中M 为bw 中的边缘点数.

5.1.6、基于计算机MATLAB 的模型求解 （1）计算机编程求解得到图片的二值矩阵

该问题所见模型基于计算机仿真的MATLAB 求解碎纸片的二值矩阵可以得到一个1272?的0、1矩阵列.

以对附件一的纵切碎纸片拼接为例,罗列出前6张碎纸图片求得的二值矩阵如下：

、

表1: 二值矩阵

000.bmp 001.bmp 002.bmp

0 0 0 0

1 1 0 0 0 1 1 0 0 1 1 1 0 0 1 1 1 1 0 1 1 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 1 1 0 1 0 1 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 0 0 0 1 0 1

1 0

1 1 1

003.bmp 004.bmp 005.bmp

0 1 0010

1 1 0100

1 1 1101

0 1 1001

0 0 0100

1 0 1011

0 1 0100

1 0 1000

1 0 1000

0 0 1010

0 1 1111

0 1 0101（2）相似度匹配结果

通过相似度匹配分析可以得到每张碎纸图片右边缘与另一张碎纸图片左边缘的相似度匹配分析。若高于90%则匹配吻合。以附件一为例：先人工干预确定第一张图片“008.bmp”为起始点，它的二值矩阵右边与“014.bmp”相似度匹配结果最高.

（3）求得的排序结果

表2：对于附件一基于相似度分析匹配求得的碎纸片自动排序

8 14 12 15 3 10 2 16 1 4 5 9 13 18 11 7 17 0 6

表3：对于附件二基于相似度分析匹配求得的碎纸片自动排序

4 7 3 8 16 19 12 1 6 2 10 14 11 9 13 1

5 18 17 5 （说明：“8”表示第九张图片“008.bmp”排在第一列；附件1和2碎纸图片的拼接结果详见附件）5.2 问题二模型建立及求解

算法流程图：

导入图像

Sobel 的算法对个碎片进行四个边的边缘化处理

对边缘化处理后的图像进行聚类分析确定阀值

二元相似性度量算法对图像进行分析

二叉树的搜索排序算法对图片排序

图2

在实际图像边缘检测问题中，图像的边缘作为图像的一种基本特征，经常被应用到较高层次的图像应用中去。图像边缘是图像最基本的特征之一，往往携带着一幅图像的大部分信息。而边缘存在于图像的不规则结构和不平稳现象中，也即存在于信号的突变点处，这些点给出了图像轮廓的位置，这些轮廓常常是我们在图像边缘检测时所需要的非常重要的一些特征条件，这就需要我们对一幅图像检测并提取出它的边缘. 5.2.1、Sobel operator

索贝尔算子（Sobel operator ）是图像处理中的算子之一，主要用作边缘检测。在技术上，它是一离散性差分算子，用来运算图像亮度函数的梯度之近似值。在图像的任何一点使用此算子，将会产生对应的梯度矢量或是其法矢量.

该算子包含两组3x3的矩阵，分别为横向及纵向，将之与图像作平面卷积，即可分别得出横向及纵向的亮度差分近似值。如果以A 代表原始图像，Gx 及Gy 分别代表经横向及纵向边缘检测的图像，其公式如下:

101202101x G A -+??

?

=-+* ? ?

-+??

121000121y G A +++?? ?=* ? ?---??

图像的每一个像素的横向及纵向梯度近似值可用以下的公式结合，来计算梯度的大

小。

22

x y

G G G =+ (17) 然后可用以下公式计算梯度方向

arctan()y

x G

G Θ= (18)

在以上例子中，如果以上的角度Θ等于零，即代表图像该处拥有纵向边缘，左方较右方暗。在边沿检测中，常用的一种模板是Sobel 算子。Sobel 算子有两个，一个是检测水平边沿的；另一个是检测垂直平边沿的。Sobel 算子对于象素的位置的影响做了加权，因此效果更好。Sobel 算子另一种形式是各向同性Sobel(Isotropic Sobel)算子，也有两个，一个是检测水平边沿的，另一个是检测垂直平边沿的。各向同性Sobel 算子和普通Sobel 算子相比，它的位置加权系数更为准确，在检测不同方向的边沿时梯度的幅度一致.

5.2.2、调整簇阀值的加速聚类方法

调整簇阀值的加速聚类方法的基本思想是:设置簇调整阀值,在每次迭代完成后,只有簇调整量大于阀值的簇才参与下一次迭代调整。簇调整阀值可以设置为簇中心迁移的调整阀值ε和簇中样本数目调整阀值θ二种情况:

(1)记某一簇在第1次迭代后聚类中心1c 与初始的簇中心0c 之间的距离为10(,)d c c ,在第i 次迭代后新的簇中心c ，与原簇中心c 之间的距离()',d c c ,如果

()

()

10,',d

c c

d c c <ε,则

以c ，为代表的簇不再参与下轮迭代调整,否则,继续参与迭代调整. (2)对于任一个簇i V ,假设样本数目i n

N k

>

,在每次迭代完成后,该簇中样本调整个数

为i N ?(增加或减少的数目绝对值),如果i N ?满足

i

i

N N ?<θ,那么,该簇将不再参与下一轮迭代调整。针对航运信息数据,通过多次Monte-Carlo 仿真实验测试,在该算法中,通常设置阀值ε在0～0.1之间,阀值θ在0～0.05之间。阀值可以只设置任意一个,也可以两个都设置。但是,值得注意的是,如果阀值设置过大,会严重影响聚类精度,因此,要合理地设置阀值的大小。这样,有了调整阀值条件后,参与迭代的样本数目大大减少,从而可以节约大量的运算时间. 算法的基本流程如下:

输入:聚类个数k,数据集X,阀值θ和ε 输出:k 个簇集V 和簇中心C

1S ：初始化簇中心{}1k

oi i c =;可以使用传统的k-means 随机初始化法,也可以采用

k-means 初始化法。

2S ：设最终簇集和中心都为空集V=Φ和C=Φ,迭代次数t=1。计算每个对象与中心

的距离,按最小距离原则将数据划分到相应的簇,形成k 个簇{}1k

oi i V =,并记录各簇中

样本个数{}1k

oi i N =,设置初始的调整标号{}{}0111k k

i i f ==.

3S ：重新计算每个簇中心11

i

i x V i

c X

V ∈=

∑和各个簇中心移动的初始位移()10,i i d c c ,

其中,i V 表示簇i V 中样本的个数.

4S :计算标号1tj f =的簇中每个样本对象与参与调整簇的中心的距离,按最小距离原

则将数据划分到相应的簇,形成k 个新簇{}tj V ,并记录各簇中样本个数{}tj N ,求得调整的样本个数tj N ?.

5S ：重新计算参与调整簇中心1tj

tj x V tj

c x V ∈=

∑和簇中心移动的位移(),tj

tj

d c c .

6S ：如果簇tj V 的中心tj c 满足

()

()

112,,tj j j j d c c d c c <ε或者

tj tj

N N ?<θ,那么,将簇tj V 中的调整标

号设置为0,即tj f =0,并将tj V 和tj c 添加到最终簇及簇中心集合中,即{}tj V V V =?和

{}tj c c c =?.

7S ：如果形成了包含k 个簇的集合,即V=k,那么,终止迭代,算法结束,否则,迭代次

数t=t+1,跳转到步骤4.

聚类结果的评价主要两个方面:聚类精度和算法时间复杂度. 算法性能分析实验如下: (1)聚类精度分析

传统的k-means 算法,所求的解往往是局部最优解,使用k-means 算法在一定程度上减少了局部最优可能性,但从聚类过程看,求的解也往往是局部最优解,本文提出的调整簇阀值的加速算法,由于采用了近似的k-means 算法,导致该算法只能找到局部近似解。适当地设置阀值,也能得到和k-means 算法相当的聚类精度. (2)时间复杂度分析

传统k-means 算法的时间复杂度为()O nkdt ,其中,n 为样本个数,k 为簇数目,d 为数据维数,t 为迭代次数. 5.2.3、相似度计算

二元相似性度量算法分析

传统的度量是采用基于两幅图像对应像素位置的位与位匹配的方法。但是在被动篡改检测中,是没有可以用来比对的原始图像的。所以,在基于二元相似性度量的图像被动篡改检测中,我们考虑比较图像位平面的二元纹理统计特征. 首先,我们用

{}

,1,...,i i k x x k =表示第i 像素和它的K 连通像素(或K 邻居像素,K=4或

K=8)的顺序关系。在本文中,我们选择K=4,表示选择N,W,S,E(上!下!左!右)四个连通像素,i 表示图像所有像素点,假定图像的大小是M 和N 。我们定义模板函数.r s x 如式(6)所示：

,,,,,1,00,2,01,3,10,4,1 1.

r r s r r s r s

r r s r r s if x and x if x and x x if x and x if x and x ==??==?=?==??==?

（19）

这样，对于一个中间像素点i ，它与其四连通像素的关系可以表示为式：

{},1,2,3,4,,,i i i i i x x x x x = （20）

对于中心像素i x ，我们定义它的协议变量(agreement)为下式：

(),1

,,1,...,4,4K

j

i i k

k a x

j j K ==

δ==∑ （21）

其中，(),a b δ是克罗内克函数(Kronecker delta function)，其定义如式：

()1,,0,a b

a b a b =?δ=?

≠?

（22）

因此对于,1, (4)

i a j =,有：

1,1,2,3,42,1,2,3,43

,1,2,3,44,1,2,3,4(,1)(,1)(,1)(,1),

(,2)(,2)(,2)(,2),(,3)(,3)(,3)(,3),(,4)(,4)(,4)(,4).

i i i i i i i i i i i

i i i i i i i i i x x x x x x x x x x x x x x x x αδδδδαδδδδαδδδδαδδδδ=+++=+++=+++=+++ （23）

很显然，j

i a 函数表示中心像素点i x 和其j 连通像素点之间的转换关系。对于中心像

素点i x ，它与其四连通像素的关系可以重新表示为式：

,1,1,1,1,11,00,

2,01,3,10,4,1 1.i i i i i i i i i if

x and x if

x and x x if x and x if x and x ==??==?=?

==??==?

,2,2,2,2,21,00,2,01,3,10,4,1 1.

i i i i i i i i i if x and x if

x and x x if

x and x if x and x ==??

==?=?

==??==?

,3,3,3,3,31,00,

2,01,3,10,4,1 1.i i i i i i i i i if

x and x if

x and x x if x and x if x and x ==??==?=?

==??==?

4,44,4,44,44,41,00,2,01,3,10,4,1 1.

i i i i i if x and x if

x and x x if x and x if x and x ==??

==?=?

==??==?

累积一协议量(accumulated agreements)的定义如下所示:

11

i

i

a MN

=

α

∑ ， 21

i

i

b MN

=α∑

31i

i

c MN

=

α

∑ ， 41

i

i

d MN

=α∑

这四个变量{},,,a b c d 是一幅二进制图像的一阶伴随变量.

5.2.4、二叉树搜索排列

二叉树搜索排序定义：二叉树的定义：空或有一个根，根有左子树、右子树；而左右子树本身又是二叉树。如图所示：

图3：

1、二叉树的性质：

性质1：在二叉树的第i 层上至多有2i-1个结点； 性质2：高度为k 的二叉树至多有2k-1个结点；

性质3：二叉树的叶子结点数n0等于度为2的结点数n2+1； 性质4：具有n 个结点的完全二叉树高度为log2n+1；

性质5：对一棵有n 个结点的完全二叉树按照从第一层（根所在的层次〕到最后一层，并且每一层都按照从左到右的次序进行编号。根结点的编号为1，最后一个结点的编号为 n.

1：对任何一个编号为i 的结点而言，它的左儿子的编号为2i( 若2i<= n) ,而右儿子的编号为2i+1(若 2i +1 <= n). 2：对任何一个编号为j 的结点而言，它的父亲结点的的编号为j/2。根结点无父结点. 2、二叉树的遍历:

设 N 代表根节点,L 代表左子树,R 代表右子树。 a. 前序（或先序）：如果二叉树为空，则操作为空：否则访问根结点；前序遍历左子树；前序遍历右子树。记为：NLR

b. 中序：如果二叉树为空，则操作为空：否则中序遍历左子树；访问根结点；中序遍历右子树。记为：LNR.

c. 后序 ：如果二叉树为空，则操作为空：否则后序遍历左子树；后序遍历右子 树；访问根结点。记为：LRN.

图4：

3、二叉树的作用:

1查找当前第K 大元素；

2查询某元素在当前所有元素中的排名.

B

C

L E

D

F B

C

L

E

D

F B X

W

(),1

,,1,...,4,4K

j

i i k

k a x

j j K ==

δ==∑ （24）

其中，

()

,a b δ是克罗内克函数(Kronecker delta function)，其定义如式：

()1,,0,a b

a b a b =?δ=?≠?

（25）

因此对于

,1,...,4

j i a j =,有：

1,1,2,3,42,1,2,3

,43,1,2,3,41,1,2,3

,4

(,1)(,1)(,1)(,1),(,2)(,2)(,2)(,2),

(,3)(,3)(,3)(,3),(,4)(,4)(,4)(,4).

i i i i i i i i i i i i i i i i i i i i x x x x x x x x x x x x x x x x αδδδδαδδδδαδδδδαδδδδ=+++=+++=+++=+++

很显然，

j

i a 函数表示中心像素点

i

x 和其j 连通像素点之间的转换关系.

对于中心像素点i

x ，它与其四连通像素的关系可以重新表示为式：

,1,1,1,1,11,00,2,01,3,10,4,1 1.i i i i i i i i i if x and x if x and x x if x and x if x and x ==??==?=?

==??==?

,2,2,2

,2,21,00,2,01,3,10,4,1 1.

i i i i i i i i i if x and x if x and x x if x and x if x and x ==??==?=?==??==?

,3,3,3

,3,31,00,2,01,3,10,4,1 1.

i i i i i i i i i if x and x if x and x x if x and x if x and x ==??==?=?==??==?

4,44,4,4

4,44,41,00,2,01,3,10,4,1 1.

i i i i i if x and x if x and x x if x and x if x and x ==??==?=?==??==?

累积一协议量(accumulated agreements)的定义如式所示

11i

i

a MN

=

α

∑，21i

i

b MN

=

α

∑

31i i

c MN

=

α∑，41i

i

d MN

=

α

∑

这四个变量

{},,,a b c d 是一幅二进制图像的一阶伴随变量.

5.2.5、构建模型的计算机仿真求解

通过基于MATLAB 编程的算法求解分析我们得到了如下表格的碎纸图片序号排列矩阵：

表4: 附件三的碎纸片排列矩阵

0 1 2 3 4 5 6 1 049 054 065 143 186 002 2 061 019 078 067 069 099 3 168 100 076 062 142 030 4 038 148 046 161 024 035 5 071 156 083 132 200 017 6 014 128 003 159 082 199 7 094 034 084 183 090 047 8 125 013 182 109 197 016 9 029 064 111 201 005 092 10 007 208 138 158 126 068 11 089 146 162 154 114 040 7 8 9 10 11 12 057 192 178 118 190 095 162 096 131 079 063 116 041 023 147 191 050 179 081 189 122 103 130 193 080 033 022 198 015 133 135 012 073 160 203 169 121 042 124 144 077 112 184 110 187 066 106 150 180 048 037 075 055 044 175 045 147 000 137 053 151 207 155 140 185 108

13 14 15 16 17 18 19

011 022 129 028 091 188 141 163 072 006 177 020 052 036 120 086 195 026 001 087 018 088 167 025 008 009 105 074 170 205 085 152 165 027 060 134 039 031 051 017 115 176 149 097 136 164 127 058 043 021 170 157 181 204 139 145 206 010 104 098 172 171 059 056 093 153 070 166 032 196 117 004 101 113 194 119 123

表5：对于附件四的碎纸片排列矩阵

0 1 2 3 4 5 6

1 191 075 011 154 190 184

2 201 148 170 196 198 094

3 086 051 107 029 040 158

4 019 194 093 141 088 121

5 159 139 001 129 063 238

6 020 041 108 116 136 073

7 208 021 007 049 061 119

8 070 084 060 014 068 174

9 132 181 095 069 167 163

10 171 042 066 205 010 157

11 081 077 128 200 131 052

7 8 9 10 11 12 13 002 104 180 064 106 004 149 113 164 078 103 091 080 101 186 098 024 117 150 005 059 126 105 155 114 176 182 151 153 053 038 123 120 175 085 036 207 135 015 076 043 199 033 142 168 062 169 054 192 137 195 008 047 172 156 096 166 188 111 144 206 003 130 074 145 083 134 055 018 056 125 140 193 087 089 048 072

14 15 16 17 18 19 032 204 065 039 067 147 026 100 006 017 028 146 058 092 030 037 046 127 022 057 202 071 165 082 050 160 187 097 203 031 045 173 079 161 179 143 133 118 189 162 197 112 023 099 122 090 185 109 034 013 110 025 027 178 035 016 009 183 152 044 012

177 124

000 102 115

（附件3和4的碎纸片拼接结果详见附件）

5.3、问题三模型的建立及求解

5.3.1、蚁群算法实现图像最优化比对搜索

蚁群算法（ACA)是受到蚂蚁群体寻找食物行为的启发而提出的一种基于蚁群的模拟进化算法。一般来讲群体随机搜索算法常用于解决特定的组合优化问题。自然界中蚂蚁搜索食物过程是一个不断聚类的过程，食物就是聚类中心，其主要思想是:将每个数据i X 看作一个蚂蚁，蚂蚁分别聚集到j 个聚类中心()

1,2,...,j C j k =，i X 到

j

C 的距离为

ij

d ，采

用欧式距离计算:

()2

1

m

ij k ik

jk k d P X

C ==

-∑ （26）

其中:m 表示每个数据特征的维数;P 加权因子，根据像素分量对信息测度影响程度设定。 假设算法中的蚂蚁具有一定的记忆能力，能根据2幅待拼接图像上互信息浓度选择转移方向，从而引导蚂蚁向互信息测度最大值的方向移动.

蚁群算法处理的问题一般具有如下特点:①搜索空间是离散的;②有一组有限的约束条件;③有一个代价函数，为搜索算法生成的解计算对应的代价，解的每一部分都会对解的代价产生影响;④有一个有限的节点集合，用于构建解;⑤有一个有限的节点间的可能转移的集合;⑥有一个节点序列的有限集合，用于表示所有的有效组合，来定义完整的搜索空间，组合的有效性、可行性由约束条件决定.

可以将图像拼接过程中图像A 中的窗口在B 搜索迭代的过程看成是“i 只蚂蚁的城市模型”，从而使得蚂蚁可以在这个模型图上进行爬行，并且保证了蚂蚁爬行所获得的每个哈密顿路都对应一个互信量最大的组合. 蚁群算法的参数进行如下设置: α为信息素的相对重要程度; β为启发式因子的相对重要程度;

ρ为信息素蒸发系数，()1-ρ表示信息素的持久性系数，且

()()()(),0

AB AB AB AB N N N N αβ

α

β?τη??ρ=τη?????????

???∑满足搜索条件，其他 （27） 其中:τ为窗口信息素含量;η为启发函数;N 为循环遍历次数.

5.3.2、基于ACA 搜索的碎纸片图像SIFT 特征点匹配

基于图像特征的方法，首先要对待配准的两幅图像进行处理，提取满足特定应用要求的特征集，然后将这两组特征集进行匹配对应，生成一组对应特征对集，最后利用这组特征对之间的对应关系估计出全局变换参数。基于图像特征的方法，在特征提取后得到的特征点的数量将会大大减少，因此可以提高配准的速度，但其配准的效果很大程度上还取决于特征点的提取精度以及特征点匹配的准确度。基于图像特征的配准方法主要困难在于如何提取和选择鲁棒的特征，以及如何对特征进行匹配，其中要克服，由于图像噪声和场景中出现遮挡现象所引起的误匹配的问题。常用的图像匹配特征有点、直线、

曲线等.

目前大多数文献都是采用点特征进行图像之间的配准，从而实现图像的拼接。基于图像特征配准的方法的主要优点是它提取了图像的显著特征，大大压缩了图像的信息量，计算量小，速度较快，而且它对图像灰度的变化具有鲁棒性。但另一方面，正是由于只有一小部分的图像灰度信息被使用了，所以这种方法对特征提取和特征匹配的错误更敏感，需要可靠的特征提取和特征一致性.

图5：基于特征点的图像拼接流程图

该问题总结的算法的搜索过程可以理解为:在2幅或多幅图上利用蚁群作为搜索窗

口，通过互信息测度的方式代替搜索出的2幅图像之间在目标、背景、边界和噪声等内容中的特征向量的过程。依据互信息量最大作为搜索并且不断改变的方向，最终朝着使互信息量最大的方向搜索，从而确定出2幅图像之间存在的匹配区域. 碎片全局拼接

如果有两个碎片证明是匹配的，那么把匹配之后的碎片组成一个新的人工碎片，重