迟到一年的建模：再探碎纸复原

前言：一年前国赛的时候，很初级地做了一下B题，做完之后还写了个《碎纸复原：一个人的数学建模》。当时就是对题目很有兴趣，然后通过一天的学习，基本完成了附件一二的代码，对附件三也只是有个概念。而今年我们上的数学建模课，老师把这道题作为大作业让我们做，于是我便再拾起了一年前的那份激情，继续那未完成的一个人的数学建模...

与去年不同的是，这次将所有代码用Python实现了，更简洁，更清晰，甚至可能更高效~~以下是论文全文。

研究背景 #

2011年10月29日，美国国防部高级研究计划局（DARPA）宣布了一场碎纸复原挑战赛（Shredder Challenge），旨在寻找到高效有效的算法，对碎纸机处理后的碎纸屑进行复原。[1]该竞赛吸引了全美9000支参赛队伍参与角逐，经过一个多月的时间，有一支队伍成功完成了官方的题目。

近年来，碎纸复原技术日益受到重视，它显示了在碎片中“还原真相”的可能性，表明我们可以从一些破碎的片段中“解密”出原始信息来。另一方面，该技术也和照片处理领域中的“全景图拼接技术”有一定联系，该技术是指通过若干张不同侧面的照片，合成一张完整的全景图。因此，分析研究碎纸复原技术，有着重要的意义。

本文以2013年全国大学生数学建模竞赛的B题为契机，初步研究了一些理想化的碎纸片的复原，基本完成了官方所提供的碎纸片的复原。我们的研究表明，完成理想化的碎纸片复原，完全是有可能的，并且能够在较快的时间内完成。

理想化假设 #

我们假设碎纸片满足以下的理想化假设：

1、每张碎纸片都是大小一样的矩形；
2、纸片的切割方式为平行切割，无倾斜；
3、碎纸片都是正面向上的；
4、纸片无缺失和混杂，即所有碎片能够拼合成一张完成的原图

模型的建立 #

在上述的理想化假设之下，我们建立了碎纸复原的数学模型。其中，我们利用并改进了贪心算法，并对匹配指标和人工干预进行了分析。

理想化的图像碎片是一个矩形，将其像素化之后就是一个像素矩阵。为了比较两张碎片是否相邻，我们主要是比较两者的边缘是否吻合，也就是比较像素矩阵的边缘列向量的相似性。而在匹配度指标这一问题上，我们对两个不同的指标——距离和相关系数——进行了比较，最终选取了距离为指标。最后我们针对附件三和附件四的拼接，加入了两种人工干预模式，结果表明这两种干预模式工作情况良好，但是有待进一步优化。

需要说明的话，为了保证算法的适用范围，我们并没有将图像进行二值化。虽然就本题而言，二值化碎片图像能获得更有效的处理，但是二值化并不适合大多数类型图片（这里的类型指的是碎片的内容，碎片同样要满足理想化的假设。）的复原处理。因此为了让我们的算法有较大的适用范围，我们直接处理原图像而不是先对图像进行二值化。

贪心算法

在碎纸的匹配中，我们主要用到了贪心算法。主要的步骤是人工选择一个起始碎片作为起点，然后把这个碎片像素矩阵的最右边缘向量跟其他碎片的最左边缘向量进行比较，选择匹配度最高的碎片进行匹配，然后以新匹配到的碎片为开始，再对余下碎片进行最优匹配，以此类推。

贪心算法是一种寻求局域最优以求达到全局最优的算法，它具有易于编程、效率高的特点。但是由于它只是每一步寻求局部最优，我们无法预料结果是否全局最优的。该算法在碎纸复原中，可能会导致两个异常：

1、不同的起点拼接的结果不同。
2、在边缘信息较少时，频繁出现匹配错误。

为了（部分地）避免贪心算法的缺点，我们在人工干预方面对贪心算法进行了改进。当然，测试表明附件一和二不需要额外的人工干预，因此我们人工干预只针对附件三和附件四设计。改进的内容将在$\ref{sec:rengongganyu}$详述。

匹配度

相邻的碎纸片

如上图所示，一般来讲，从图像的连通性考虑，相邻的两个碎纸片的边缘像素应当是相似、甚至是相同的。这样很容易联想到两个边缘的匹配指标为两个向量的距离：
$$\begin{equation}d^2=|\boldsymbol{x}-\boldsymbol{y}|^2=\sum_{i=1}^{n}(x_i-y_i)^2\end{equation}$$

如果两个边缘向量的距离为0，我们有把握说这两个碎片是匹配的。但是当两个向量距离很大时，我们有多大把握说这两个碎片不匹配呢？我们知道，一般的图像不仅仅具有连通性，还具有渐变性。对于一张一般的图片来说，相邻两个列向量可以认为是线性渐变的，这样一来，匹配度似乎用“线性相关性”衡量更为妥当。线性相关系数[3]的计算公式为：
$$\begin{equation}\rho_{X,Y}={\mathrm{cov}(X,Y) \over \sigma_X \sigma_Y} ={E[(X-\mu_X)(Y-\mu_Y)] \over \sigma_X\sigma_Y}\end{equation}$$

线性相关系数的计算公式比距离公式复杂很多。我们对同一组碎片采用两个指标进行拼接测试，经过对不同内容碎片的测试，发现总体来讲采用线性相关系数的效果优于采用距离公式的效果。然而，线性相关系数的优势体现在照片内容的碎纸（该碎片为我们队伍自行生成用来测试算法。）复原上，而对于官方提供的文字内容的碎纸复原，两者效果并无明显区别。最终，从简单性考虑，我们选择了向量的距离作为了我们的匹配度指标。

人工干预

本节为附件三和附件四而设计。为了方便人工干预的进行，我们在拼接碎纸片的时候，在每张碎片上标出了文件名，但没有修改原来的碎片文件。

首先，如果直接把附件一、二的代码用于附件三，会出现下面的效果（局部）：

附件三的初步拼接

所有的碎纸片被拼成了一个长条，长条中有拼接正确的地方，也有不少拼接错误。换个角度看，目前碎纸片被拼成了一堆（大概20~30 个）稍大的碎片块，也就是说，我们可以把已经拼接好的确定下来，然后再对这些碎块进行拼接。这样一来，从200多个碎片块变成了20 多个碎片块，自然大大减少了工作量。这里需要人工从图片上辨别出哪些已经拼接好的，然后记录到代码之中。事实上，这就相当于一个人工聚类的过程。（在这里，我们并没有采用算法自动聚类，具体原因将在后面的困难分析中解释。）

除此之外，我们引入了“人工排除”。在有两个距离很近的时候，贪心算法可能会出现匹配错误的现象。但是我们必须认识到这种情况不会频繁出现，当边缘信息足够时，或者已经对碎片进行了聚类之后，匹配错误是一个低概率的事情。因此，假如我们通过观察发现出现匹配的错误时，就可以人工地排除这种匹配方式，使得程序在剩下的匹配方案中选择最优的。我们有理由相信，在有限次（很少次）的排除之后，会使得程序自动匹配到正确的碎片。

当有了进一步的正确拼接信息之后，可以继续把正确的拼接补充到代码中，然后继续排除、人工拼接。重复这个过程，就能够在“可以接受的”时间内完成碎纸的复原。

问题求解 #

这里我们汇总我们所得到的求解结果。

附件一和附件二

这里只以附件一为例子进行说明。附件一的拼接代码位于附录中。首先，随意指定一个起始图片，即代码中的$m$值，然后根据拼接结果判断出正确的起始图片（008.bmp），然后修改$m$值，重新再拼一次即可。测试表明不需要额外的人工干预即可完成。最终得到附件一的拼接顺序为

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

附件二的碎片复原过程跟附件一基本一样，得到附件二的拼接顺序为：

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

附件三和附件四

附件三的代码由附件一的代码修改而成，主要是加入了文件名标注、人工拼接和人工排除三个方面的功能。直接运行改代码，得到拼接序列（部分）

[14, 128, 3, 159, 82, 199, 7, 208, 29, 64, 111, 201, 5, 92, 180, 48, 37, 75, 38, 148, ...

观察得知，[14, 128, 3, 159, 82, 199]这个序列是拼接正确的，[7, 208]是拼接正确的，[29, 64, 111, 201, 5, 92, 180, 48, 37, 75]也是拼接正确的，等等。于是依次在代码的“known=[[] for i in range(0,num)]”下一行加入：

known[14] = [128, 3, 159, 82, 199]
known[7] = [208]
known[29] = [64, 111, 201, 5, 92, 180, 48, 37, 75]
...

并且观察得知199后面匹配7是错误的，208后面匹配29也是错误的，75后面也不能匹配38，因此，在“impossible=[[] for i in range(0,num)]”下一行加入：

impossible[199] = [7]
impossible[208] = [29]
impossible[75] = [38]
...

修改完成后，重新运行脚本，得到新的拼接序列。把拼接序列中新增加的正确方案，以及人工可以容易找到的正确方案都加到known里边去，把明显的不正确拼接都加入到impossible里边。重复该操作。在两个小时左右的时间里，可以完成附件三每一行的拼接。也就是说，先拼好每一行，生成每一行的图片，然后进行纵向的拼接。

同样的做法可以完成附件四的拼接。

困难分析 #

对于附件三和附件四的复原，我们各用上了两个半小时的时间，通过简单的自动聚类，可以缩短用时，但经计算只能是常数级的优化。这表明碎纸的自动复原是一个相当困难的问题。下面就我们队的研究过程，分析该问题的难度所在。

聚类难以奏效

虽然聚类可以降低部分工作量，但实际上聚类能够起到的作用是非常小的。在我们的代码中，即使我们人工找出了在同一行的19个碎片，然后用算法进行自动拼接，拼接效果也不理想，这是由于边缘信息过少所引起的后果，难以避免。

其次，聚类局限性太大。赛题的碎片是文章切割而成，文章具有比较明显的规律（字号一样，行距一样等），因此有办法进行初步的聚类。然而，对于一般的碎片，比如照片的碎片，就没有类似的聚类方式。因此聚类是难以奏效的。聚类在本题中只能起到常数级优化作用。

纵向信息难以利用

我们队伍在完成附件三的每一行的拼接后，得到了11个横向的切条。但是即使以这11个横向切条为原料，自动拼接还是出了很多错误，无法自动拼接完成。也就是说，即使已经完成了每一行的拼接，纵向的边缘所能够提供的信息还是相当少，如果一开始在200多个小碎片时就想利用纵向信息，几乎是不可能的。因此，几乎可以说，在完成拼接的过程中，只有横向的信息可以使用，这使得难度变大。

前车之鉴

碎纸复原本身就是一个相当艰难的问题，本文开头所提到的DARPA的碎纸复原挑战赛，胜利的队伍也用了一个多月的时间才完成任务。[4]虽然两个问题的难度无法相比并论，但这侧面说明了碎纸复原的难度之大。因此，本文复原附件三所花的两个多小时，应该属于在合理的时间之内。

改进方向 #

我们构思了一个改进算法的方向，但由于时间限制，无法进行详细验证，简单论述如下。

人的肉眼判断两个碎片是否相邻是比较容易的，而转化为计算机算法则比较困难。我们注意到，当我们在观察一幅碎片的边缘的时候，并非只是看边缘的一个像素的（我们肉眼分辨率并没有这么精细），而是看到边缘的几列像素的平均效果。因此，需要改进的是匹配度不能仅以边缘的一列像素进行计算，而应当综合考虑边缘的多列像素。

根据肉眼识别的思路，我们可以用图像的连通度作为匹配度的指标。然而，图片的连通度还没有明确的定义。如何根据多列向量判断图片是否连通，还需要深入分析。仅以此抛砖引玉，望读者不吝赐教。

参考文献 #

[1] 果壳网：http://www.guokr.com/article/78259/
[2] 《离散数学结构》，Bernard Kolman等 著，罗平 译
[3] 皮尔逊积矩相关系数：http://zh.wikipedia.org/zh/皮尔逊积矩相关系数
[4] Solidot：http://www.solidot.org/story?sid=27531

代码列表 #

以下代码的运行环境为Python 3.4，Win32版本，需要安装对应的Numpy和Pillow库。

附件1和2的代码

import os
import numpy
from PIL import Image

#获取当前目录列表
images=os.listdir(os.getcwd())
images.remove('c.py') #脚本的文件名
num=len(images) #图片数目

hang=Image.open(images[0]).size[1] #图片高度
lie=Image.open(images[0]).size[0] #图片宽度
bianyuan=numpy.zeros((hang,2*num)) #矩阵用于储存边缘信息。

#打开图片，获取边缘值。
for i in range(0,num):
    img=Image.open(images[i])
    bianyuan[:,2*i+1]=numpy.array(img)[:,0] #!!! 左边缘信息，放到奇数列
    bianyuan[:,2*i]=numpy.array(img)[:,lie-1] #!!! 右边缘信息，放到偶数列
#边缘值获取完毕

i=0;m=31 #m是起始点
xulie=[m] #储存拼接顺序
temp=[k for k in range(0,num)]
while i<num-1:
    m1=temp[m]
    temp.remove(temp[m])
    bijiao=numpy.zeros(len(temp)) #用以比较矩阵
    for j in range(0,len(temp)):
        bijiao[j]=numpy.sum((bianyuan[:,2*m1]-bianyuan[:,2*temp[j]+1])**2) #每个右边缘跟所有的左边缘比较
    m=numpy.argmin(bijiao) #差别最小意味着吻合度最高
    xulie.append(temp[m])
    i=i+1
print(xulie) #比较完成，输出序列

comb_img=Image.new ("RGBA", (lie*num, hang), (255, 0, 0)) #新图片，用来合成
j=0
for i in xulie:
    img=Image.open(images[i])
    region=img.crop((0,0,lie,hang))
    comb_img.paste(region,(lie*j,0,lie*(j+1),hang))
    j=j+1
comb_img.show()

附件3和4的代码

import os
import numpy
from PIL import Image
from PIL import ImageDraw
from PIL import ImageFont

#获取当前目录列表
images=os.listdir(os.getcwd())
images.remove('c.py')
num=len(images) #图片数目

hang=Image.open(images[0]).size[1] #图片高度
lie=Image.open(images[0]).size[0] #图片宽度
bianyuan=numpy.zeros((hang,2*num)) #矩阵用于储存边缘信息。

#打开图片，获取边缘值。
for i in range(0,num):
    img=Image.open(images[i])
    bianyuan[:,2*i+1]=numpy.array(img)[:,0] #!!! 左边缘信息，放到奇数列
    bianyuan[:,2*i]=numpy.array(img)[:,lie-1] #!!! 右边缘信息，放到偶数列
#边缘值获取完毕

known=[[] for i in range(0,num)] #已知的拼接

impossible=[[] for i in range(0,num)] #否定的拼接

i=0;m=14 #m是起始点
m1=m
xulie=[m] #储存拼接顺序
temp=[k for k in range(0,num)]
temp.remove(m1)

for kn in known:
    for knn in kn:
        try:
            temp.remove(knn) #删除已知情况
        except:pass

while i<num-1:
    if len(known[m1]) != 0:
        m2=0
        for kn in known[m1]:
            xulie.append(kn)
            m2=kn
            i=i+1
        m1=m2
    else:
        temp1=temp[:]
        for imp in impossible[m1]:
            try:
                temp.remove(imp)
            except:pass
        bijiao=numpy.zeros(len(temp)) #用以比较矩阵
        for j in range(0,len(temp)):
            bijiao[j]=numpy.sum((bianyuan[:,2*m1]-bianyuan[:,2*temp[j]+1])**2) #每个右边缘跟所有的左边缘比较
        m=numpy.argmin(bijiao) #差别最小意味着吻合度最高
        xulie.append(temp[m])
        i=i+1
        m1=temp[m]
        temp=temp1[:]
        temp.remove(m1) #过河拆桥
print(xulie) #比较完成，输出序列

#加入水印
myfont = ImageFont.truetype("C:\Windows\Fonts\simsun.ttc",36)    #水印字体与大小

comb_img=Image.new ("RGBA", (lie*num, hang), (255, 0, 0)) #新图片，用来合成
j=0
for i in xulie:
    img=Image.open(images[i])
    d=ImageDraw.Draw(img) #加入水印
    d.ink = 150 #水印颜色
    d.text((0,0),str(i),font=myfont) #加入水印
    region=img.crop((0,0,lie,hang))
    comb_img.paste(region,(lie*j,0,lie*(j+1),hang))
    j=j+1
comb_img.show()

最后修改结果（附件三）

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

impossible=[[] for i in range(0,num)] #否定的拼接
impossible[199]=[7,29,38,49,61,62,67,71,80,89,94,125]
impossible[1]=[146,129,102,134]
impossible[79]=[71,146,89,94]
impossible[123]=[94,125]
impossible[176]=[7]
impossible[13]=[135,94]
impossible[160]=[143,168,146,25,94,7,29,38,49]
impossible[95]=[168,129]
impossible[124]=[136,8,46,138,129]
impossible[58]=[161,46,182,83]
impossible[46]=[9]
impossible[97]=[144]
impossible[25]=[168]

最后修改结果（附件四）

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

impossible=[[] for i in range(0,num)] #否定的拼接
impossible[16]=[30,178,195,150,186,2,19,70,81,132]
impossible[122]=[109,197,32]
impossible[137]=[32]
impossible[204]=[4,54]
impossible[22]=[82]
impossible[12]=[120,149,159]
impossible[175]=[24,1]
impossible[158]=[30,195]
impossible[190]=[54,4]
impossible[144]=[197,32,178,195]
impossible[121]=[167,169]
impossible[150]=[130]
impossible[176]=[88]
impossible[141]=[182,70,81]
impossible[124]=[3]
impossible[138]=[102,0]
impossible[65]=[169,184,2]
impossible[165]=[9]
impossible[119]=[32,195,70]
impossible[27]=[177,195]
impossible[169]=[4]
impossible[62]=[126]
impossible[193]=[85,177]
impossible[162]=[70]
impossible[31]=[149]
impossible[184]=[202]
impossible[85]=[102]
impossible[57]=[88]

转载到请包括本文地址：https://kexue.fm/archives/3134

更详细的转载事宜请参考：《科学空间FAQ》

打赏

微信打赏

支付宝打赏

因为网站后台对打赏并无记录，因此欢迎在打赏时候备注留言。
你还可以点击这里或在下方评论区留言来告知你的建议或需求。