數(shù)學(xué)建模：帶你再探碎紙片復(fù)原

數(shù)學(xué)建模：帶你再探碎紙片復(fù)原

研究背景

2011年10月29日，美國(guó)國(guó)防部高級(jí)研究計(jì)劃局（DARPA）宣布了一場(chǎng)碎紙復(fù)原挑戰(zhàn)賽（Shredder Challenge），旨在尋找到高效有效的算法，對(duì)碎紙機(jī)處理后的碎紙屑進(jìn)行復(fù)原。[1]該學(xué)術(shù)活動(dòng)吸引了全美9000支參賽隊(duì)伍參與角逐，經(jīng)過(guò)一個(gè)多月的時(shí)間，有一支隊(duì)伍成功完成了官方的題目。

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

本文以2013年全國(guó)大學(xué)生數(shù)學(xué)建模學(xué)術(shù)活動(dòng)的B題為契機(jī)，初步研究了一些理想化的碎紙片的復(fù)原，基本完成了官方所提供的碎紙片的復(fù)原。我們的研究表明，完成理想化的碎紙片復(fù)原，完全是有可能的，并且能夠在較快的時(shí)間內(nèi)完成。

2

理想化假設(shè)

我們假設(shè)碎紙片滿(mǎn)足以下的理想化假設(shè)：

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

3

模型的建立

在上述的理想化假設(shè)之下，我們建立了碎紙復(fù)原的數(shù)學(xué)模型。其中，我們利用并改進(jìn)了貪心算法，并對(duì)匹配指標(biāo)和人工干預(yù)進(jìn)行了分析。

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

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

4

貪心算法

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

貪心算法是一種尋求局域最優(yōu)以求達(dá)到全局最優(yōu)的算法，它具有易于編程、效率高的特點(diǎn)。但是由于它只是每一步尋求局部最優(yōu)，我們無(wú)法預(yù)料結(jié)果是否全局最優(yōu)的。該算法在碎紙復(fù)原中，可能會(huì)導(dǎo)致兩個(gè)異常：

1、不同的起點(diǎn)拼接的結(jié)果不同。
2、在邊緣信息較少時(shí)，頻繁出現(xiàn)匹配錯(cuò)誤。

為了（部分地）避免貪心算法的缺點(diǎn)，我們?cè)谌斯じ深A(yù)方面對(duì)貪心算法進(jìn)行了改進(jìn)。當(dāng)然，測(cè)試表明附件一和二不需要額外的人工干預(yù)，因此我們?nèi)斯じ深A(yù)只針對(duì)附件三和附件四設(shè)計(jì)。

5

匹配度

如上圖所示，一般來(lái)講，從圖像的連通性考慮，相鄰的兩個(gè)碎紙片的邊緣像素應(yīng)當(dāng)是相似、甚至是相同的。這樣很容易聯(lián)想到兩個(gè)邊緣的匹配指標(biāo)為兩個(gè)向量的距離：

如果兩個(gè)邊緣向量的距離為0，我們有把握說(shuō)這兩個(gè)碎片是匹配的。但是當(dāng)兩個(gè)向量距離很大時(shí)，我們有多大把握說(shuō)這兩個(gè)碎片不匹配呢？我們知道，一般的圖像不僅僅具有連通性，還具有漸變性。對(duì)于一張一般的圖片來(lái)說(shuō)，相鄰兩個(gè)列向量可以認(rèn)為是線性漸變的，這樣一來(lái)，匹配度似乎用“線性相關(guān)性”衡量更為妥當(dāng)。線性相關(guān)系數(shù)[3]的計(jì)算公式為：

線性相關(guān)系數(shù)的計(jì)算公式比距離公式復(fù)雜很多。我們對(duì)同一組碎片采用兩個(gè)指標(biāo)進(jìn)行拼接測(cè)試，經(jīng)過(guò)對(duì)不同內(nèi)容碎片的測(cè)試，發(fā)現(xiàn)總體來(lái)講采用線性相關(guān)系數(shù)的效果優(yōu)于采用距離公式的效果。然而，線性相關(guān)系數(shù)的優(yōu)勢(shì)體現(xiàn)在照片內(nèi)容的碎紙（該碎片為我們隊(duì)伍自行生成用來(lái)測(cè)試算法。）復(fù)原上，而對(duì)于官方提供的文字內(nèi)容的碎紙復(fù)原，兩者效果并無(wú)明顯區(qū)別。最終，從簡(jiǎn)單性考慮，我們選擇了向量的距離作為了我們的匹配度指標(biāo)。

6

人工干預(yù)

本節(jié)為附件三和附件四而設(shè)計(jì)。為了方便人工干預(yù)的進(jìn)行，我們?cè)谄唇铀榧埰臅r(shí)候，在每張碎片上標(biāo)出了文件名，但沒(méi)有修改原來(lái)的碎片文件。

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

所有的碎紙片被拼成了一個(gè)長(zhǎng)條，長(zhǎng)條中有拼接正確的地方，也有不少拼接錯(cuò)誤。換個(gè)角度看，目前碎紙片被拼成了一堆（大概20~30 個(gè)）稍大的碎片塊，也就是說(shuō)，我們可以把已經(jīng)拼接好的確定下來(lái)，然后再對(duì)這些碎塊進(jìn)行拼接。

這樣一來(lái)，從200多個(gè)碎片塊變成了20 多個(gè)碎片塊，自然大大減少了工作量。這里需要人工從圖片上辨別出哪些已經(jīng)拼接好的，然后記錄到代碼之中。事實(shí)上，這就相當(dāng)于一個(gè)人工聚類(lèi)的過(guò)程。（在這里，我們并沒(méi)有采用算法自動(dòng)聚類(lèi)，具體原因?qū)⒃诤竺娴睦щy分析中解釋。）

除此之外，我們引入了“人工排除”。在有兩個(gè)距離很近的時(shí)候，貪心算法可能會(huì)出現(xiàn)匹配錯(cuò)誤的現(xiàn)象。但是我們必須認(rèn)識(shí)到這種情況不會(huì)頻繁出現(xiàn)，當(dāng)邊緣信息足夠時(shí)，或者已經(jīng)對(duì)碎片進(jìn)行了聚類(lèi)之后，匹配錯(cuò)誤是一個(gè)低概率的事情。

因此，假如我們通過(guò)觀察發(fā)現(xiàn)出現(xiàn)匹配的錯(cuò)誤時(shí)，就可以人工地排除這種匹配方式，使得程序在剩下的匹配方案中選擇最優(yōu)的。我們有理由相信，在有限次（很少次）的排除之后，會(huì)使得程序自動(dòng)匹配到正確的碎片。

當(dāng)有了進(jìn)一步的正確拼接信息之后，可以繼續(xù)把正確的拼接補(bǔ)充到代碼中，然后繼續(xù)排除、人工拼接。重復(fù)這個(gè)過(guò)程，就能夠在“可以接受的”時(shí)間內(nèi)完成碎紙的復(fù)原。

7

問(wèn)題求解

這里我們匯總我們所得到的求解結(jié)果。

附件一和附件二

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

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

附件二的碎片復(fù)原過(guò)程跟附件一基本一樣，得到附件二的拼接順序?yàn)椋?/span>

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

附件三和附件四

附件三的代碼由附件一的代碼修改而成，主要是加入了文件名標(biāo)注、人工拼接和人工排除三個(gè)方面的功能。直接運(yùn)行改代碼，得到拼接序列（部分）

[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]這個(gè)序列是拼接正確的，[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是錯(cuò)誤的，208后面匹配29也是錯(cuò)誤的，75后面也不能匹配38，因此，在“impossible=[[] for i in range(0,num)]”下一行加入：

impossible[199]=[7]

impossible[208]=[29]

impossible[75]=[38]...

修改完成后，重新運(yùn)行腳本，得到新的拼接序列。把拼接序列中新增加的正確方案，以及人工可以容易找到的正確方案都加到known里邊去，把明顯的不正確拼接都加入到impossible里邊。

重復(fù)該操作。在兩個(gè)小時(shí)左右的時(shí)間里，可以完成附件三每一行的拼接。也就是說(shuō)，先拼好每一行，生成每一行的圖片，然后進(jìn)行縱向的拼接。

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

困難分析

對(duì)于附件三和附件四的復(fù)原，我們各用上了兩個(gè)半小時(shí)的時(shí)間，通過(guò)簡(jiǎn)單的自動(dòng)聚類(lèi)，可以縮短用時(shí)，但經(jīng)計(jì)算只能是常數(shù)級(jí)的優(yōu)化。這表明碎紙的自動(dòng)復(fù)原是一個(gè)相當(dāng)困難的問(wèn)題。下面就我們隊(duì)的研究過(guò)程，分析該問(wèn)題的難度所在。

聚類(lèi)難以奏效

雖然聚類(lèi)可以降低部分工作量，但實(shí)際上聚類(lèi)能夠起到的作用是非常小的。在我們的代碼中，即使我們?nèi)斯ふ页隽嗽谕恍械?9個(gè)碎片，然后用算法進(jìn)行自動(dòng)拼接，拼接效果也不理想，這是由于邊緣信息過(guò)少所引起的后果，難以避免。

其次，聚類(lèi)局限性太大。賽題的碎片是文章切割而成，文章具有比較明顯的規(guī)律（字號(hào)一樣，行距一樣等），因此有辦法進(jìn)行初步的聚類(lèi)。然而，對(duì)于一般的碎片，比如照片的碎片，就沒(méi)有類(lèi)似的聚類(lèi)方式。因此聚類(lèi)是難以奏效的。聚類(lèi)在本題中只能起到常數(shù)級(jí)優(yōu)化作用。

縱向信息難以利用

我們隊(duì)伍在完成附件三的每一行的拼接后，得到了11個(gè)橫向的切條。但是即使以這11個(gè)橫向切條為原料，自動(dòng)拼接還是出了很多錯(cuò)誤，無(wú)法自動(dòng)拼接完成。也就是說(shuō)，即使已經(jīng)完成了每一行的拼接，縱向的邊緣所能夠提供的信息還是相當(dāng)少，如果一開(kāi)始在200多個(gè)小碎片時(shí)就想利用縱向信息，幾乎是不可能的。因此，幾乎可以說(shuō)，在完成拼接的過(guò)程中，只有橫向的信息可以使用，這使得難度變大。

前車(chē)之鑒

碎紙復(fù)原本身就是一個(gè)相當(dāng)艱難的問(wèn)題，本文開(kāi)頭所提到的DARPA的碎紙復(fù)原挑戰(zhàn)賽，勝利的隊(duì)伍也用了一個(gè)多月的時(shí)間才完成任務(wù)。[4]雖然兩個(gè)問(wèn)題的難度無(wú)法相比并論，但這側(cè)面說(shuō)明了碎紙復(fù)原的難度之大。因此，本文復(fù)原附件三所花的兩個(gè)多小時(shí)，應(yīng)該屬于在合理的時(shí)間之內(nèi)。

改進(jìn)方向

我們構(gòu)思了一個(gè)改進(jìn)算法的方向，但由于時(shí)間限制，無(wú)法進(jìn)行詳細(xì)驗(yàn)證，簡(jiǎn)單論述如下。

人的肉眼判斷兩個(gè)碎片是否相鄰是比較容易的，而轉(zhuǎn)化為計(jì)算機(jī)算法則比較困難。我們注意到，當(dāng)我們?cè)谟^察一幅碎片的邊緣的時(shí)候，并非只是看邊緣的一個(gè)像素的（我們?nèi)庋鄯直媛什](méi)有這么精細(xì)），而是看到邊緣的幾列像素的平均效果。因此，需要改進(jìn)的是匹配度不能僅以邊緣的一列像素進(jìn)行計(jì)算，而應(yīng)當(dāng)綜合考慮邊緣的多列像素。

根據(jù)肉眼識(shí)別的思路，我們可以用圖像的連通度作為匹配度的指標(biāo)。然而，圖片的連通度還沒(méi)有明確的定義。如何根據(jù)多列向量判斷圖片是否連通，還需要深入分析。僅以此拋磚引玉，望讀者不吝賜教。

參考文獻(xiàn)

[1] 果殼網(wǎng)：http://www.guokr.com/article/78259/
[2] 《離散數(shù)學(xué)結(jié)構(gòu)》，Bernard Kolman等 著，羅平 譯
[3] 皮爾遜積矩相關(guān)系數(shù)：http://zh.wikipedia.org/zh/皮爾遜積矩相關(guān)系數(shù)
[4] Solidot：http://www.solidot.org/story?sid=27531

代碼列表

以下代碼的運(yùn)行環(huán)境為Python 3.4，Win32版本，需要安裝對(duì)應(yīng)的Numpy和Pillow庫(kù)。

附件1和2的代碼

importosimportnumpyfromPIL importImage

#獲取當(dāng)前目錄列表

images=os.listdir(os.getcwd())

images.remove('c.py')#腳本的文件名

num=len(images)#圖片數(shù)目

hang=Image.open(images[0]).size[1]#圖片高度

lie=Image.open(images[0]).size[0]#圖片寬度

bianyuan=numpy.zeros((hang,2*num))#矩陣用于儲(chǔ)存邊緣信息。

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

img=Image.open(images[i])

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

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

i=0;m=31#m是起始點(diǎn)

xulie=[m]#儲(chǔ)存拼接順序

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

附件3和4的代碼

importosimportnumpyfromPIL importImagefromPIL importImageDrawfromPIL importImageFont

#獲取當(dāng)前目錄列表

images=os.listdir(os.getcwd())

images.remove('c.py')

num=len(images)#圖片數(shù)目

hang=Image.open(images[0]).size[1]#圖片高度

lie=Image.open(images[0]).size[0]#圖片寬度

bianyuan=numpy.zeros((hang,2*num))#矩陣用于儲(chǔ)存邊緣信息。

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

img=Image.open(images[i])

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

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

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

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

i=0;m=14#m是起始點(diǎn)

m1=m

xulie=[m]#儲(chǔ)存拼接順序

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

temp.remove(m1)

forkn inknown:

forknn inkn:

try:

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

except:pass

whilei

最后修改結(jié)果（附件三）

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]

最后修改結(jié)果（附件四）

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]

?