高教社杯全国大学生数学建模竞赛B题论文

碎纸片的拼接复原
摘要
本文利用Manhattan距离，聚类分析，图像处理等方法解决了碎纸片的拼接复原问题。

由于碎纸机产生的碎纸片是边缘规则且等大的矩形，此时碎纸片拼接方法就不能利用碎片边缘的尖角特征等基于边界几何特征的拼接方法，而要利用碎片内的字迹断线或碎片内的文字位置搜索与之匹配的相邻碎纸片。

拼接碎片前利用数学软件MATLAB软件对碎片图像进行数据化处理，得到对应的像素矩阵，后设置阈值对像素矩阵进行二值化处理，得到相应的0-1矩阵。

下面分别对三个问题的解决方法和算法实现做简单的阐述：
问题一，分别对附件1和附件2的碎片数据进行处理得到相应的0-1矩阵，依次计算某个0-1矩阵最右边一列组成向量与其他所有0-1矩阵的最左边向量的Manhattan距离，可以得到某个最小距离值、说明最小距离值对应的碎片是可与基准碎片拼接的，最终得到碎片拼接完整的图像。

问题二，同样对于附件3和附件4中的碎片数据进行处理得到相应的数值矩阵，并计算得到每个碎片顶部空白高度和文字高度，即指每行像素点都为255
的行数、一行中存在像素点为非255的行数，根据空白高度和文字高度对碎片进行聚类分类，聚类阀值取3像素，得到11组像素矩阵，进而得到11类可能在同一行的碎片类。

其中对附件4中的英文的处理中，我们还采用水平像素投影累积的方法，进一步分类出可能在同一行的碎片类。

用问题一的方法，计算Manhattan 距离可以对每一类碎片按次序排列好，得到11行已经排列好的碎片，再应用曼哈顿距离在竖直方向上进行聚合得到完整的图像。

问题三，首先，对于附件5中的碎片数据我们采用正反相接，本文将b面最左边的一列像素拼接到a面最右边的一列像素的下面，构成360×1的向量，再把其他的碎片采用相同的办法得到360×1的向量，再用问题一的方法，计算出各碎片之间的Manhattan距离。

其次，根据每个碎片顶部的空白高度或者文字高度对碎片进行区间分类，得到22组矩阵，然后应用曼哈顿距离将得到的22组矩阵聚成两类，每类各包含两面的11组矩阵，最后利用Manhattan距离在竖直方向上进行聚合得到完整的图像。

本文最后，我们根据算法的效率实现进行了改进和优化，实现算法的移植性、灵活性、运行效率等得以提升。

关键词：曼哈顿距离，聚类分析，二值化处理
一、问题重述
破碎文件的拼接在司法物证复原、历史文献修复以及军事情报获取等领域都有着重要的应用。

传统上，拼接复原工作需由人工完成，准确率较高，但效率很低。

特别是当碎片数量巨大，人工拼接很难在短时间内完成任务。

随着计算机技术的发展，人们试图开发碎纸片的自动拼接技术，以提高拼接复原效率。

请讨论以下问题：
1. 对于给定的来自同一页印刷文字文件的碎纸机破碎纸片（仅纵切），建立碎纸片拼接复原模型和算法，并针对附件1、附件2给出的中、英文各一页文件
的碎片数据进行拼接复原。

如果复原过程需要人工干预，请写出干预方式及干预的时间节点。

复原结果以图片形式及表格形式表达。

2. 对于碎纸机既纵切又横切的情形，请设计碎纸片拼接复原模型和算法，并针对附件3、附件4给出的中、英文各一页文件的碎片数据进行拼接复原。

如果复原过程需要人工干预，请写出干预方式及干预的时间节点。

复原结果表达要求同上。

3. 上述所给碎片数据均为单面打印文件，从现实情形出发，还可能有双面打印文件的碎纸片拼接复原问题需要解决。

附件5给出的是一页英文印刷文字双面打印文件的碎片数据。

请尝试设计相应的碎纸片拼接复原模型与算法，并就附件5的碎片数据给出拼接复原结果，结果表达要求同上。

二、问题分析
我们从附件中的碎片数据可知由于碎纸机产生的碎纸片边缘是规则的，此时碎纸片计算机拼接方法就不能利用碎片边缘的尖点特征、尖角特征、面积特征等基于边界几何特征的拼接方法，而要利用碎片内的字迹断线或碎片内的文字内容是否匹配搜索与之匹配的相邻碎纸片并进行拼接。

首先，我们对碎片内图像进行数据化处理，得到对应的像素值矩阵；然后，我们设置阈值对像素值矩阵进行二值化处理得到相应的数值矩阵；最后，由于曼哈顿距离公式计算快、数值小，数值矩阵与数值矩阵之间应用最小曼哈顿距离对碎纸片进行拼接复原。

问题一中碎纸机破碎纸片只有纵切，每页纸被切为19条碎片，经过处理可以得到19个数值矩阵。

对于每个数值矩阵，我们依次取出最左边一列从上至下各格的值组成一个向量，同样我们依次取出最右边一列从上至下各格的值组成一个向量。

计算出每一数值矩阵的左边向量与所有非同源数值矩阵的右边向量的曼哈顿距离，再将得到的距离值进行排序，当某个距离值最小时、说明相应的左边向量与右边向量的匹配率最大，则该距离对应的左、右边认为是可拼接的。

若得到的最小距离值不止一个，则此时需要进行人工干预。

问题二是对碎纸机既纵切又横切的情形进行讨论，比问题一多了横切条件，此时每页纸被切为209个碎片。

首先，我们利用文件最左边碎片与最上面碎片的特殊性对这209个碎片进行聚类，得到两类特殊的碎片，分别是文件最左边一列碎片和最上面一行碎片，然后类似于问题一的处理方法，应用最小曼哈顿距离对每一类碎片按正确顺序拼接，此后对其余碎片再应用最小曼哈顿距离逐一进行拼接，直至剩余所有的碎片都拼接上。

问题三中，题目要求考虑双面打印文件的碎纸拼接复原问题的解决方案，此时每页纸虽然也是被切为209个碎片，但每个碎片却有正反两面，因此经过处理
得到418个数值矩阵，，此时我们分别对每一面各自进行类似问题一的处理，然后综合每一面的聚类情况再应用最小曼哈顿距离对双面碎纸片进行拼接复原。

三、模型假设
1. 假设碎纸机破碎纸片（纵切或横切）得到的碎纸片是规则且边缘是整齐的等大的矩形；
2.假设我们对文档碎纸片拼接复原不考虑碎片边缘的尖点特征、尖角特征、面积特征等基于边界几何特征；
3.假设附件中给出的所有中、英文文件中的文字排版是按标准格式排版的。

4.假设附件中给出的所有中、英文字符都是统一格式，且内容为普通文章。

5.1 问题一（曼哈顿距离）
➢ 模型一的建立
题目要求对于给定的来自同一页印刷文字文件的碎纸机破碎纸片（仅纵切）建立碎纸片拼接复原模型和算法，并且要对中、英文各一页文件的碎片数据分别进行拼接复原。

首先，我们利用数学软件MATLAB 软件将19条碎片数据化，得到19个像素值矩阵，像素值的变化范围是从0变化到255，此时我们设置127τ=为阈值对像素值矩阵进行二值化处理，当矩阵某位置像素值小于等于τ时，则将对应位置的数值设为0；当矩阵某位置像素值大于τ时，则将对应位置的数值设为127。

这样我们就得到19个二值化了的数值矩阵
i
A ，对于每个数值矩阵
i
A ，我
们依次取出最左边一列从上至下各格的值组成一个向量，记为
i
X ，同样的我们
依次取出最右边一列从上至下各格的值组成一个向量，记为i Y 。

计算出每一数值矩阵的左边向量与所有非同源数值矩阵的右边向量的曼哈顿距离
(,)
i j d X Y 。

➢ 模型一的求解 对于得到的向量
12(,,...,)(1,2,...,)
T
i i i ik X x x x k m ==和向量
12(,,...,)(1,2,...,)
T
i i i ik Y y y y k n ==，两向量的曼哈顿距离为
1
(,)||
(,1,2,...,)
n
i j ik jk k d X Y x y i j m i j ==-=≠∑且。

可求出附件1碎片与碎片
之间的曼哈顿距离，如下表所示。

表1 附件1碎片与碎片间的曼哈顿距离
从而可得到附件1碎片序号按复原后顺序如下表所示。

表2 附件1碎片序号复原后顺序
附件1碎片复原图片如附录中图8.1所示。

同法可求出附件2碎片与碎片之间的曼哈顿距离，如下表所示。

表3 附件2碎片与碎片间的曼哈顿距离
从而可得到附件2碎片序号按复原后顺序如下表所示。

附件2碎片复原图片如附录中图8.2所示。

问题一人工干预情况如下表所示。

5.2 问题二（Manhattan 距离）
➢ 模型二的建立
在中文文件中，两个连续的汉字中间的空白间隔所占像素宽度与其左边或者
右边的汉字所占像素宽度的比值最大的约为
2
13，则对于每一行文字，碎纸机纵
切未切到文字的概率为
2
13，对于每两行文字碎纸机纵切未切到文字的概率为
4
169，而对于每三行文字碎纸机纵切未切到文字的概率更小，可以忽略不计，所
以对于总共209个碎片，每个碎片上面的文字至少有两行（碎片上不完整的一行也算一行），所以出现某个碎片上面的文字完全没被碎纸机切割到（即文字完整
无缺）的概率至多为
4
169，我们把这样的碎片称之为干扰碎片。

我们知道，整篇文件的最上面一行字的上边缘是空白的，我们可以利用此特殊性对209个碎纸片进行聚类，可以得到一个特殊的类，即碎纸片上边缘为空白的类，此类碎纸片个数大于等于11；出现个数大于11的情形即为混入上面提到
的干扰碎片，此概率最大不超过
4
169，可知此类碎纸片应该拼接在文件最上面一
行，应用最小曼哈顿距离对此类碎片按正确顺序拼接。

同理可聚类出另一个特殊的类，即碎纸片左边缘为空白、拼接在文件最左边一列的类，并且也应用最小曼哈顿距离对此类碎片按正确顺序拼接。

然后以此拼接好的第一行和第一列碎片为基准，再应用最小曼哈顿距离拼接其余剩下的碎片，最后拼接复原出原中文文件。

在英文文件中，一个英文单词中两个连续的英文字母中间的空白间隔所占像
素宽度与其左边或者右边的英文字母所占像素宽度的比值最大的约为
1
11，则对
于每一行英文单词，碎纸机纵切未切到英文单词的概率为
1
11，对于每两行英文
单词碎纸机纵切未切到英文单词的概率为
1
121，而对于每三行英文单词碎纸机纵
切未切到英文单词的概率为，然后同上述中文文件的分析过程可知，此时对拼接在文件最左边一列归类时混入上面提到的干扰碎片的概率最大不超过
，最后拼接复原出原英文文件。

➢模型二的求解
我们利用SPSS软件根据每个碎片顶部空白高度或者文字高度的不同，应用聚类分析方法将碎片聚成11类，结果如下图所示。

图1 根据碎片顶部文字高度聚类
图2 根据碎片顶部空白高度聚类
结合上面的聚类图，可得出附件3的乱序矩阵，如下表所示。

表6 附件3的乱序矩阵
49 22 129 178 118 143 188 192 57 141 91 190 28 186 2 54 11 95 65 61 79 116 78 72 20 69 52 163 177 36 99 96 19 67 63 162 131 6 168 179 1 30 23 142 191 87 147 62 76 86 195 18 26 120 100 41 50 38 167 74 46 103 148 88 35 9 8 24 193 161 105 189 25 130 122 81 71 205 27 200 60 85 15 33 156 170 198 132 17 202 152 83 165 133 80 14 115 159 128 199 12 107 176 82 160 73 31 51 203 169 3 135 39 134 94 58 90 149 77 42 34 112 144 136 124 84 164 97 47 127 121 183 43 125 13 187 173 139 66 150 197 182 16 106 181 145 109 21 110 184 157 204 29 10 104 172 55 48 171 5 98 37 206 59 92 201 64 44 180 111 75
7 0 93 32 56 175 153 166 196 137 45 208 174 68 158 138 53 70 126 89 151 114 140 102 207 155 101 146 194 119 4 117 40 123 108 154 185 113
同样的方法可得出附件4的乱序矩阵，如下表所示。

表7 附件4的乱序矩阵
然后我们先求出附件3碎片与碎片之间的曼哈顿距离，从而得到附件3碎片序号按复原后顺序如下表所示。

表8 附件3碎片序号复原后顺序
同法我们再求出附件4碎片与碎片之间的曼哈顿距离，从而得到附件4碎片序号按复原后顺序如下表所示。

表9 附件4碎片序号复原后顺序
问题二人工干预情况如下表所示。

5.3 问题三（曼哈顿距离）
➢ 模型三的建立
问题三在问题二的基础上继续加大碎片拼接复原难度，此时我们对双面碎纸片进行类似问题一的处理，得到418个数值矩阵，我们根据每个碎片顶部的空白高度或者文字高度对碎片进行区间分类，得到22组矩阵，再根据曼哈顿距离将得到的22组矩阵聚成两类，每类各包含某一面的11组矩阵，然后综合每一面的聚类情况再应用最小曼哈顿距离对双面碎纸片进行拼接复原。

然后再利用曼哈顿距离对碎纸片在竖直方向上进行聚合得到最终图像。

➢ 模型三的求解
问题三的解决方法与问题二的类似，不过我们分两步进行聚类分析。

第一步，我们根据每个碎片顶部空白高度的不同进行聚类，第二步，我们根据每个碎片底部空白高度的不同进行聚类。

然后我们选取第一、二步聚类产生的公共类，若得到的公共类数量小于22类，则再从单独由第一步聚类产生的类中选取，直到数量达到22类。

对于这22类碎纸片，我们再利用问题二的方法聚成两组，
每组数量都为11类。

后面类似模型二的处理过程，结果顺序如下表所示。

问题三人工干预情况如下表所示。

1.模型的评价
对于问题一，由于题目中给的样本较为简单，所以模型一能很好的解决附件1、附件2给出的中、英文文件碎纸片拼接复原问题。

对于问题二，模型二也能较好的解决问题，但模型二也有不足之处。

比如模型二只考虑根据每个碎片顶部的空白高度和文字高度对碎片进行区间分类，分为11组矩阵。

而没有综合考虑每个碎片顶部与底部的空白高度和文字高度对碎片进行区间分类，因此分类准确率降低。

对于问题三，由于每个碎片都有正反两面，而且模型三对碎片聚类时只对每一面单独进行分析，所以模型三解决问题时错误匹配的数量明显多于前面两题。

因此，我们除了挖掘算法潜力，还能对模型三做出进一步的改进，比如对碎片聚类时综合对两面进行分析以及对某面分析时综合考虑已经排好顺序的第一行与第一列碎片等都能进一步优化模型，减少错误匹配的数量，提高效率，增加模型的适应性。

2.模型的推广
我们建立的模型在处理碎纸片较大且碎纸片数量不是很多的时候，模型可以较好的解决问题，但在实际应用中，通常会涉及碎纸片被切割得很细很小，并且要对大量碎纸片数据进行管理和处理工作。

所以我们要进一步优化算法和程序结构，改善模型，真正建立起快速有效的计算机辅助碎纸片自动拼接复原模型，从而才能将此模型广泛地应用到我们的实际生活中。

七、参考文献
[1]. 张翠. 基于点线的文档图片数字水印与碎片拼接[D]. 青岛:中国海洋大学, 2011. [26-34]。

[2]. 张艳. 图像拼接技术在文档图像扭曲识别中的应用与研究[D]. 北京:北方工业大学, 2011. [23-29]。

[3]. 贾海燕, 朱良家, 周宗潭, 胡德文. 一种碎纸自动拼接中的形状匹配方法[J]. 计算机仿真, 2006, 23(11): [180-183]。

[4]. 汪晓银, 邹庭荣. 数学软件与数学实验[M]. 北京: 科学出版社,2008. [1-27]。

[5]. 张铮,杨文平,石博强,李海鹏.?Matlab程序设计与实例应用[M].?北京:?中国铁道出版社,2003.[276-306]。

[6]. 姚文敏. 数字图像处理[M]. 北京:机械工业出版社, 2006. [24-45]。

八、附录
图像拼接结果
图8.1 附件1碎片复原图片
图8.2 附件2碎片复原图片
图8.3 附件3碎片初次拼接结果图
图8.4 附件3碎片复原图片
图8.5 附件4碎片初次拼接结果图
图8.6 附件4碎片复原图片
图8.7 附件5碎片复原图片一面
图8.8 附件5碎片复原图片另一面
程序源代码
问题一
function Q1() %处理问题1的相关程序
clear all;
close all；
clc; %清理显示屏幕
var=0;
N=18;
m=0;
n=0;
pi=0;
pj=0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %导入图像部分
for i=0:N
if(i<10)
filename=sprintf('00%d.bmp',i); %构造需要导入的文件名字
else %而相应的文件存放在工作目录中
filename=sprintf('0%d.bmp',i); %构造需要导入的文件 end
I=imread(filename); %导入图片文件
Im{i+1}=I; %将变量Im用cell格式存放数据，后续调用这个变量
end
[m n]=size(I); %提取矩阵的行与列 for i=0:N
IM{i+1}=im2bw(Im{i+1}); %Im变量中所有的图像数据进行二值化
i=i+1;
end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=0:N
subplot(1,N+1,i+1)
imshow(Im{i+1}); %显示没有经过处理时候的情况
end
%}
for i=0:N
var=-1;
for j=0:N
if(i~=j)
if(var==-1)
var=sum(abs(IM{i+1}(:,72)-IM{j+1}(:,1)));
pi=i;
pj=j;
else
if(sum(abs(double(IM{i+1}(:,72))-double(IM{j+1}(:,1))))<var)
var=sum(abs(double(IM{i+1}(:,72))-double(IM{j+1}(:,1))));
%通过残差绝对值再求和判断
pi=i; %左边位置的图片编号 pj=j; %右边位置的图片编号 end
end
end
end
res(i+1,:)=[pi pj var];
end
res %显示的结果为：左边位置图片编号，右边位置图片编号，曼哈顿距离值
k=0;
for i=0:N
if(res(i+1,3)==0) %将误差为0的编号放在最前面
sq(k+1)=res(i+1,2);
end
end
for i=0:N-1
for j=0:N
if(sq(i+1)==res(j+1,1))
sq(i+2)=res(j+1,2); %对序列进行排序
end
end
end
sq %显示排列的顺序ResIm=Im{sq(1)+1};
for i=1:N
ResIm=[ResIm,Im{sq(i+1)+1}]; %将图片比编号按顺序进行拼接
end
hold off
imshow(ResIm) %显示图片拼接结果hold on
%title('附件1拼接结果图')
问题二—函数1
function Q21()
clear all;
close all;
clc;
N=208;
var=0;
m=0;
n=0;
pi=0;
pj=0;
R=11;
C=19;
for i=0:N
if(i<10)
filename=sprintf('00%d.bmp',i); %构造文件的名字
else if(i<100)
filename=sprintf('0%d.bmp',i);
else
filename=sprintf('%d.bmp',i);
end;
end
I=imread(filename); %将像素文件读入程序
[m n]=size(I);
Im{i+1}=I;
end
for i=0:N
IM{i+1}=im2bw(Im{i+1}); %像素矩阵二值化
end
for i=1:N+1
[H(i) B(i)]=LocateHeight(IM{i}); %调用LocateHeight函数求一张图片中最上方白色部分的宽度
L(i)=i-1;
end
RE(:,1)=L; %L表示图像标识位置
RE(:,2)=H; %H表示图片中顶部的空白宽度或是文字宽度
RE(:,3)=B;
OUT=sortrows(RE,3);
temp=sum(OUT(:,3));
O1=OUT(1:(N+1-temp),:);
O2=OUT((N+2-temp):N,:);
sortrows(O1,2) %按照图像中顶部空白的图像高度进行排序
sortrows(O2,2) %按照图像中顶部黑色的图像高度进行排序
function [res1 res2]=LocateHeight(II) %LocateHeight函数，计算图像顶部文字高度或者空白位置，为.m内部函数
t1=-1;
t2=-1;
for i=0:179
if(i==0&&sum(II(i+1,:))==72)
t1=i;
end
if(i==0&&sum(II(i+1,:))~=72)
t2=i;
end
if(t1~=-1&&i>t1&&sum(II(i+1,:))~=72)
res1=i-t1;
res2=1;
break;
end
if(t2~=-1&&i>t1&&sum(II(i+1,:))==72)
res1=i-t2;
res2=0;
break;
end
end
问题二—函数2
function Q22()
clear all;
close all;
clc;
N=208; %图像数量
var=0; %方差
m=0;
n=0;
pi=0;
pj=0;
R=11;
C=19;
%Z=[91,95,118,190,178,186,143,129,49,188,57,11,2,141,192,28,65,22,54; 163,177,63,72,162,116,78,61,52,79,6,67,99,20,69,96,131,36,19;23,142,1 68,26,100,76,147,30,195,86,191,62,179,1,50,120,41,87,18;25,130,148,46 ,74,24,9,193,103,81,8,88,38,122,35,161,189,105,167;152,202,71,165,133 ,60,198,200,17,33,205,85,132,83,170,156,15,80,27;128,135,31,73,115,16 9,3,203,82,199,107,160,12,159,14,134,176,39,51;183,77,121,127,90,58,1 44,42,136,43,124,34,94,164,84,47,97,112,149;204,184,16,66,187,125,110 ,157,139,106,150,197,173,21,109,181,13,145,182;201,5,111,206,48,98,10 ,37,55,75,104,64,171,180,172,92,29,59,44;137,126,7,174,68,56,158,153, 70,138,0,208,45,175,53,166,93,32,196;123,151,140,108,114,4,113,207,15 4,117,194,146,101,40,89,102,185,155,119;];
Z=[191 147 11 67 204 106 104 2 184 190 154 65 39 180 149 75 4 64 32 201 80 91 101 198 100 94 6 26 196 103 113 17 28 148 78 146 164 170 86 59 51 117 24 29 92 58 186 107 46 158 127 40 98 37 5 30 150 19 121 114 57 88 176 82 194 151 22 155 182 126 141 105 93 202 71 165 159 203 187 53 1 120 160 153 85 31 97 138 129 50 139 123 63 38 175 20 136 76 135 36 43 143 41 173 79 199 179 161 45 73 207 116 108 15 208 189 168 49 112 118 169 33 142 119 54 197 61 62 7 133 21 192 162 70 23 68 109 195 60 84 99 174 90 137 8 96 156 47 172 14 122 185 132 163 181 110 25 188 206 27 95 166 69 178 3 111 130 34 167 13 144 171 16 74 134 152 35 83 55 9 157 205 42 145 44 66 56 18 183 10
81 177 52 72 48 12 89 102 140 87 128 125 0 115 193 77 200 131 124
] %利用聚类得到的矩阵表for i=0:N
if(i<10)
filename=sprintf('00%d.bmp',i);
else if(i<100)
filename=sprintf('0%d.bmp',i);
else
filename=sprintf('%d.bmp',i);
end;
end
I=imread(filename);
[m n]=size(I);
Im{i+1}=I;
end
for i=0:N
IM{i+1}=im2bw(Im{i+1}); %将导入的像素矩阵进行二极值化变换，得到0-1矩阵
end
a=1;
S=Z;
[m,n]=size(S);
for i=1:R
for j=1:C
IMG{i,j}=IM{S(i,j)+1}; %将同一行的像素矩阵放置在一个4维数组中
end
end
for r=1:R %在水平方向上根据两个矩阵的曼哈顿距离进行拼接
for j=2:C-1
var=-1;
for i=j:C
if(var==-1)
var=sum(abs(IMG{r,j-1}(:,72)-IMG{r,i}(:,1)));
p=i;
else if(sum(abs(IMG{r,j-1}(:,72)-IMG{r,i}(:,1)))<var) var=sum(abs(IMG{r,j-1}(:,72)-IMG{r,i}(:,1)));
p=i;
end
end
end
temp1=IMG{r,p};
IMG{r,p}=IMG{r,j};
IMG{r,j}=temp1;
end
end
for i=1:R
ResIm{i}=IMG{i,1};
for j=2:C
ResIm{i}=[ResIm{i},IMG{i,j}]; %在水平方向对矩阵进行合并
end
end
RES=ResIm{1};
for i=2:R
RES=[RES;ResIm{i}]; %在竖直方向将矩阵合并成大矩阵
end
%hold off
imshow(RES); %显示结果图
%hold on
title('附件4的初次拼接效果')。