數學建模：帶你再探碎紙片復原

數學建模：帶你再探碎紙片復原

研究背景

2011年10月29日，美國國防部高級研究計劃局（DARPA）宣布了一場碎紙復原挑戰賽（Shredder Challenge），旨在尋找到高效有效的算法，對碎紙機處理后的碎紙屑進行復原。[1]該學術活動吸引了全美9000支參賽隊伍參與角逐，經過一個多月的時間，有一支隊伍成功完成了官方的題目。

近年來，碎紙復原技術日益受到重視，它顯示了在碎片中“還原真相”的可能性，表明我們可以從一些破碎的片段中“解密”出原始信息來。另一方面，該技術也和照片處理領域中的“全景圖拼接技術”有一定聯系，該技術是指通過若干張不同側面的照片，合成一張完整的全景圖。因此，分析研究碎紙復原技術，有著重要的意義。

本文以2013年全國大學生數學建模學術活動的B題為契機，初步研究了一些理想化的碎紙片的復原，基本完成了官方所提供的碎紙片的復原。我們的研究表明，完成理想化的碎紙片復原，完全是有可能的，并且能夠在較快的時間內完成。

2

理想化假設

我們假設碎紙片滿足以下的理想化假設：

1、每張碎紙片都是大小一樣的矩形；
2、紙片的切割方式為平行切割，無傾斜；
3、碎紙片都是正面向上的；
4、紙片無缺失和混雜，即所有碎片能夠拼合成一張完成的原圖

3

模型的建立

在上述的理想化假設之下，我們建立了碎紙復原的數學模型。其中，我們利用并改進了貪心算法，并對匹配指標和人工干預進行了分析。

理想化的圖像碎片是一個矩形，將其像素化之后就是一個像素矩陣。為了比較兩張碎片是否相鄰，我們主要是比較兩者的邊緣是否吻合，也就是比較像素矩陣的邊緣列向量的相似性。而在匹配度指標這一問題上，我們對兩個不同的指標——距離和相關系數——進行了比較，最終選取了距離為指標。最后我們針對附件三和附件四的拼接，加入了兩種人工干預模式，結果表明這兩種干預模式工作情況良好，但是有待進一步優化。

需要說明的話，為了保證算法的適用范圍，我們并沒有將圖像進行二值化。雖然就本題而言，二值化碎片圖像能獲得更有效的處理，但是二值化并不適合大多數類型圖片（這里的類型指的是碎片的內容，碎片同樣要滿足理想化的假設。）的復原處理。因此為了讓我們的算法有較大的適用范圍，我們直接處理原圖像而不是先對圖像進行二值化。

4

貪心算法

在碎紙的匹配中，我們主要用到了貪心算法。主要的步驟是人工選擇一個起始碎片作為起點，然后把這個碎片像素矩陣的最右邊緣向量跟其他碎片的最左邊緣向量進行比較，選擇匹配度最高的碎片進行匹配，然后以新匹配到的碎片為開始，再對余下碎片進行最優匹配，以此類推。

貪心算法是一種尋求局域最優以求達到全局最優的算法，它具有易于編程、效率高的特點。但是由于它只是每一步尋求局部最優，我們無法預料結果是否全局最優的。該算法在碎紙復原中，可能會導致兩個異常：

1、不同的起點拼接的結果不同。
2、在邊緣信息較少時，頻繁出現匹配錯誤。

為了（部分地）避免貪心算法的缺點，我們在人工干預方面對貪心算法進行了改進。當然，測試表明附件一和二不需要額外的人工干預，因此我們人工干預只針對附件三和附件四設計。

5

匹配度

如上圖所示，一般來講，從圖像的連通性考慮，相鄰的兩個碎紙片的邊緣像素應當是相似、甚至是相同的。這樣很容易聯想到兩個邊緣的匹配指標為兩個向量的距離：

如果兩個邊緣向量的距離為0，我們有把握說這兩個碎片是匹配的。但是當兩個向量距離很大時，我們有多大把握說這兩個碎片不匹配呢？我們知道，一般的圖像不僅僅具有連通性，還具有漸變性。對于一張一般的圖片來說，相鄰兩個列向量可以認為是線性漸變的，這樣一來，匹配度似乎用“線性相關性”衡量更為妥當。線性相關系數[3]的計算公式為：

線性相關系數的計算公式比距離公式復雜很多。我們對同一組碎片采用兩個指標進行拼接測試，經過對不同內容碎片的測試，發現總體來講采用線性相關系數的效果優于采用距離公式的效果。然而，線性相關系數的優勢體現在照片內容的碎紙（該碎片為我們隊伍自行生成用來測試算法。）復原上，而對于官方提供的文字內容的碎紙復原，兩者效果并無明顯區別。最終，從簡單性考慮，我們選擇了向量的距離作為了我們的匹配度指標。

6

人工干預

本節為附件三和附件四而設計。為了方便人工干預的進行，我們在拼接碎紙片的時候，在每張碎片上標出了文件名，但沒有修改原來的碎片文件。

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

所有的碎紙片被拼成了一個長條，長條中有拼接正確的地方，也有不少拼接錯誤。換個角度看，目前碎紙片被拼成了一堆（大概20~30 個）稍大的碎片塊，也就是說，我們可以把已經拼接好的確定下來，然后再對這些碎塊進行拼接。

這樣一來，從200多個碎片塊變成了20 多個碎片塊，自然大大減少了工作量。這里需要人工從圖片上辨別出哪些已經拼接好的，然后記錄到代碼之中。事實上，這就相當于一個人工聚類的過程。（在這里，我們并沒有采用算法自動聚類，具體原因將在后面的困難分析中解釋。）

除此之外，我們引入了“人工排除”。在有兩個距離很近的時候，貪心算法可能會出現匹配錯誤的現象。但是我們必須認識到這種情況不會頻繁出現，當邊緣信息足夠時，或者已經對碎片進行了聚類之后，匹配錯誤是一個低概率的事情。

因此，假如我們通過觀察發現出現匹配的錯誤時，就可以人工地排除這種匹配方式，使得程序在剩下的匹配方案中選擇最優的。我們有理由相信，在有限次（很少次）的排除之后，會使得程序自動匹配到正確的碎片。

當有了進一步的正確拼接信息之后，可以繼續把正確的拼接補充到代碼中，然后繼續排除、人工拼接。重復這個過程，就能夠在“可以接受的”時間內完成碎紙的復原。

7

問題求解

這里我們匯總我們所得到的求解結果。

附件一和附件二

這里只以附件一為例子進行說明。附件一的拼接代碼位于附錄中。首先，隨意指定一個起始圖片，即代碼中的mm值，然后根據拼接結果判斷出正確的起始圖片（008.bmp），然后修改mm值，重新再拼一次即可。測試表明不需要額外的人工干預即可完成。最終得到附件一的拼接順序為

[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的代碼

importosimportnumpyfromPIL importImage

#獲取當前目錄列表

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))#矩陣用于儲存邊緣信息。

#打開圖片，獲取邊緣值。fori inrange(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 fork inrange(0,num)]whilei

附件3和4的代碼

importosimportnumpyfromPIL importImagefromPIL importImageDrawfromPIL importImageFont

#獲取當前目錄列表

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))#矩陣用于儲存邊緣信息。

#打開圖片，獲取邊緣值。fori inrange(0,num):

img=Image.open(images[i])

bianyuan[:,2*i+1]=numpy.array(img)[:,0]#!!! 左邊緣信息，放到奇數列

bianyuan[:,2*i]=numpy.array(img)[:,lie-1]#!!! 右邊緣信息，放到偶數列#邊緣值獲取完畢

known=[[]fori inrange(0,num)]#已知的拼接

impossible=[[]fori inrange(0,num)]#否定的拼接

i=0;m=14#m是起始點

m1=m

xulie=[m]#儲存拼接順序

temp=[k fork inrange(0,num)]

temp.remove(m1)

forkn inknown:

forknn inkn:

try:

temp.remove(knn)#刪除已知情況

except:pass

whilei

最后修改結果（附件三）

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]

?