#113 基于图像信号分析的碎纸片的拼接复原

Created by 欧 长坤 & 郭 瞾阳 & 黄 沐 on 13-9-16. Copyright (c) 2013年 欧 长坤 & 郭 瞾阳 & 黄 沐. All rights reserved.

摘 要

碎片化程度不同的碎纸片,其复原难度各有不同。好的复原筛选算法有效的帮助人们减少复原的时间,提高复原的效率。

针对三个问题中复原碎纸片的难度递进关系,本文同样建立了具备递进关系的三个数学模型,并给出其算法。分别为信号相似性算法、边界积分匹配算法和行距匹配算法。针对复原文件的特性,本文将碎片的复原过程分划为三个独立的代码模块,分别为整行筛选模块、整行拼接模块和各行拼接模块。每个模块适当组合使用三个算法,达到精准筛选与匹配。

在问题一中,由于附件1和附件2所提供的待复原图片数量较少,因此直接采用信号相似性分析,来拼接待复原的图片,拼接结果显示,拼接的匹配率达到100%,无需人工干预。 在问题二中,随着碎片数的增加,为了减小拼接问题的复杂度,先筛选出整行碎片,则将该问题转化为问题一进行研究。而整行筛选的方式,中文碎片和英文碎片的处理方式上有细微差异。对于中文碎片而言,采用了行距分析算法等,半自动化的完成了整块碎片文件的拼接,其间人工干预3次。对于英文碎片而言,实际拼接过程中直接使用信号相似性算法来筛选整行更为有效,匹配率为84.21%,完成整块209张图像碎片拼接只需人工干预33次。 对于更复杂的双面碎片问题三而言,综合运用了三个模块进行全面筛选匹配,将匹配率优化为63.31%,完成整个518张图像碎片的拼接复原需要人工干预95次。

关键词: 一维信号, 相似性分析, 边际积分, 行距提取, OpenCV

限于篇幅,下面谈谈摘要中提到的核心模块和核心算法。

一维信号相似分析模型

文[1]对于测量信号\(x(t),y(t)\),使用了Hilbert空间上的内积理论获得了对两个信号的相似性判定的衡量方法。设\(x=[\xi_{1},\xi_{2},...]\)\(x=[\eta_{1},\eta_{2},...]\)为两个离散信号,为使能量误差

\[ \begin{equation} Q=\Vert x-\lambda_{0}y\Vert ^{2}=\sum_{n}(\xi_{n}-\lambda_{0}\eta_{n})^{2} \end{equation} \]

只需要令\(\frac{dQ}{d\lambda_{0}}=0\),得到

\[ \begin{equation} \lambda_{0}=\frac{\sum_{n}\xi_{n}\eta_{n}}{\sum_{n}\eta_{n}^{2}} \end{equation} \]

此时对应的最小能量误差为

\[ \begin{equation} Q_{min}=\sum_{n}(\xi_{n}-\lambda_{0}\eta_{n})^{2} = \sum_{n}\xi_{n}^{2} - \frac{\sum_{n}\xi_{n}\eta_{n}}{\sum_{n}\eta_{n}^{2}} \end{equation} \]

而最小相对能量误差为:

\[ \begin{equation} Q_{min} = 1 - \frac{\left(\sum_{n}\xi_{n}\eta_{n}\right)^{2}}{\sum_{n}\xi_{n}^{2}\sum_{n}\eta_{n}^{2}} = 1-R_{xy}^{2} \end{equation} \]

根据Cauchy-Schwarz不等式得结论:规范相关系数

\[ \begin{equation} R_{xy} = \frac{\sum_{n}\xi_{n}\eta_{n}}{\sqrt{\sum_{n}\xi_{n}^{2}\sum_{n}\eta_{n}^{2}}} \end{equation} \]

它表征了\(x,y\)的相似性,且规范相关系数的范数越大信号越相似,两信号的能量误差越趋近于0。

信号相似性算法:

这里对序列的相关性系数计算,直接编写代码,具体实现见略,判定算法描述如下:

  1. 图像\(i\)为需要查找后(或上、下、前)继图像\(j\)位于指定的待选空间(集合)中,提取图像\(i\)与图像\(j\)的边际信号序列(若是寻找后继,则图像\(i\)提取右侧信号,图像\(j\)提取左侧信号,若是寻找前继,则图像\(i\)提取左侧信号,图像\(j\)提取右侧信号,寻找上继和下继可类推),为了防止计算时导致的数据溢出,将信号序列归一化(将\([0,255]\)上的整型值转化为\([0,1]\)上的双精度浮点值);
  2. 计算两者的规范相关系数\(R_{xy}\)
  3. 初始化筛选下界\(temp = 0.5\)和标记变量\(flag = 0\)
  4. \(R_{xy}>temp\),则\(temp = R_{xy}\)\(flag = j\)\(j++\)
  5. 若遍历完所有图片,则第\(flag\)张图片即为与图像\(i\)最匹配的图像。否则,返回(1)。

边际积分模型

对图像边界上图像信号的相似性分析理论上是可行的,但是在实际操作中由于筛选量较大、时间复杂度较高,因此我们通过引入边际积分的概念来辅助分析图片的前驱与后继。

图像\(G\)在边际\(left\)上的边际积分\(f\)

\[ \begin{equation} f=count(G,D,n) \end{equation} \]

其中,\(count\)是一个从\((G,D,n)\xrightarrow{f} N\)的映射,\(G\)为图像空间,\(D\)为积分方向,\(n\)为对在该方向上的积分灰度值,\(D\)的取值仅有\(\{left,right,top,bottom\}\)\(n\)的取值为\([0,255]\)上的整数,积分值为在积分方向上的灰度值为\(n\)的像素个数。 当\(G\)为某图片时,若令\(D=left\)\(n=255\),则返回的是该图片左边第一列的全部像素中白色像素的个数。

边际积分匹配算法: 开源计算机视觉库OpenCV提供了非常方便的访问图像数据的结构,更容易统计该列某颜色值像素个数。其算法描述为:

  1. 对于阈值\(\Phi_{1}\)(可根据实际结果进行调整)将图像二值化,保证图像中仅有黑白两色。
  2. 对于图像\(i\)(不妨设为当前行的首张图片),计算它在边际\(right\)上的边际积分,记为\(p\)
  3. 在剩余图像中取一图像\(j\),计算它在\(left\)边际上的边际积分,记为\(q\)
  4. 计算的值\(\Vert p-q \Vert\),若小于给定阈值\(\Phi_{2}\)并保存\(j\)的值在数组\(space\)中,作为一个可能后继。
  5. 若未遍历完所有图片,则返回(3),否则进入(6)。
  6. 数组\(space\)保存的值即为筛选出的可能匹配对象的集合。

部分关键代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
// Created by 欧 长坤 on 13-9-16.
// Copyright (c) 2013年 欧 长坤 . All rights reserved.
// 接口:arr_1,信号1所在数组;arr_2,信号2所在数组;lenth_1,信号1所在数组长度;lenth_2,信号2所在数组长度
// 使用:函数返回信号1与信号2的规范相关系数。
double R(double *arr_1, int lenth_1, double *arr_2, int lenth_2)
{
double sum_1 = 0;
for (int i = 0; i < lenth_1; ++i)
{
sum_1 += arr_1[i]*arr_2[i];
}
double sum_2 = 0;
double sum_3 = 0;
for (int i = 0; i < lenth_2; ++i)
{
sum_2 += arr_1[i]*arr_1[i];
sum_3 += arr_2[i]*arr_2[i];
}
double sum_4 = sqrt(sum_2*sum_3); // need to #include
return (sum_1/sum_4);
}
// Created by 郭 瞾阳 on 13-9-16.
// Copyright (c) 2013年 郭 瞾阳 . All rights reserved.
// 信号提取部分
int shang_1(IplImage *img,double* black_arr)
{
int i = 0;
double x1,x2,x3;
uchar*par = (uchar *)(img-&imageData);
for(int x=0;xwidth;x++,i++)
{
x1 = (double)par[3*x]/255;
x2 = (double)par[3*x+1]/255;
x3 = (double)par[3*x+2]/255;
black_arr[i] = (x1+x2+x3)/3;
}
return i;
}
int xia_1(IplImage * img,double* black_arr)
{
int i = 0;
double x1,x2,x3;
uchar* par = (uchar*)(img-&imageData+(img-&height-1)*(img-&widthStep));
for(int x=0 ;xwidth;x++,i++)
{
x1 = (double)par[3*x]/255;
x2 = (double)par[3*x+1]/255;
x3 = (double)par[3*x+2]/255;
black_arr[i] = (x1+x2+x3)/3;
}
return i;
}
int zuo_1(IplImage *img,double* black_arr)
{
int i = 0;
double x1,x2,x3;
for(int y = 0;yheight;y++,i++)
{
uchar * par =(uchar *)(img-&imageData+y*img-&widthStep);
x1 = (double)par[0]/255;
x2 = (double)par[1]/255;
x3 = (double)par[2]/255;
black_arr[i] = (x1+x2+x3)/3;
}
return i;
}
int you_1(IplImage *img,double* black_arr)
{
int i = 0;
int x = img-&width-1
double x1,x2,x3;
for(int y = 0;yheight;y++,i++)
{
uchar * par = (uchar *)(img-&imageData+y*img-&widthStep);
x1 = (double)par[3*x]/255;
x2 = (double)par[3*x+1]/255;
x3 = (double)par[3*x+2]/255;
black_arr[i] = (x1 + x2 + x3)/3;
}
return i;
}
// Created by 欧 长坤 on 13-9-16.
// Copyright (c) 2013年 欧 长坤 . All rights reserved.
// main()函数中具体操作的部分
char filepath[100] = &quot;C:\\B\\5\\000a.bmp&quot;;
char filepath_2[100] = &quot;C:\\B\\5\\001a.bmp&quot;;
char filepath_3[100] = &quot;C:\\B\\5\\001b.bmp&quot;;
std::ofstream fout(&quot;a.txt&quot;);
for (int dota = 0; dota <= 208; dota++)
{
numChangeFilePath(filepath,dota);
IplImage* img = cvLoadImage(filepath);
IplImage* img_do = cvLoadImage(filepath_2);
IplImage* img_do_2 = cvLoadImage(filepath_3);
double arr_0[500];
int arr_0_lenth = you_1(img,arr_0);
double arr_1[500];
int arr_1_lenth = zuo_1(img_do,arr_1);
double arr_2[500];
int arr_2_lenth = zuo_1(img_do,arr_2);
double temp = 0.5;
double temp_2 = 0.5;
int flag_1 = 0;
int flag_2 = 0;
for (int i = 1; i <= 208; ++i)
{
numChangeFilePath(filepath_2,i);
cvReleaseImage(&img_do);
img_do = cvLoadImage(filepath_2);
arr_1_lenth = zuo_1(img_do,arr_1);
if (R(arr_0,arr_0_lenth,arr_1,arr_1_lenth)&temp)
{
temp = R(arr_0,arr_0_lenth,arr_1,arr_1_lenth);
flag_1 = i;
}
}
for (int i = 1; i <= 208; ++i)
{
numChangeFilePath(filepath_3,i);
cvReleaseImage(&img_do_2);
img_do_2 = cvLoadImage(filepath_3);
arr_2_lenth = zuo_1(img_do_2,arr_2);
if (R(arr_0,arr_0_lenth,arr_2,arr_2_lenth)&temp_2)
{
temp_2 = R(arr_0,arr_0_lenth,arr_2,arr_2_lenth);
flag_2 = i;
}
}
//printf(&quot;(img:%d,next_a:%d,next_b_%d)&quot;,dota,flag_1,flag_2);
fout << &quot;(&quot; << dota << &quot;, &quot; << flag_1 << &quot;, &quot; << flag_2 << &quot;)&quot; << std::endl;
}
// Created by 郭 瞾阳 on 13-9-16.
// Copyright (c) 2013年 郭 瞾阳 . All rights reserved.
// 接口:img,进行边界积分的图片
// 使用:对图片左侧的边界进行积分,返回积分值
int countLeftEdgeIntegral(IplImage* img)
{
int black = 0;
for(int y = 0; yheight; y++)
{
unsigned char * par =(unsigned char *)(img-&imageData+y*img-&widthStep);
if(par[0]==0 && par[1]==0 && par[2]==0)
black++;
}
return black;
}
// 接口:img,进行边界积分的图片
// 使用:对图片右侧的边界进行积分,返回积分值
int countRightEdgeIntegral(IplImage* img)
{
int black = 0;
int x = img-&width-1;
for(int y = 0;yheight;y++)
{
unsigned char * par = (unsigned char *)(img-&imageData+y*img-&widthStep);
if(par[3*x+0]==0 && par[3*x+1]==0 && par[3*x+2]==0)
black++;
}
return black;
}
// 接口:img,进行边界积分的图片
// 使用:对图片上侧的边界进行积分,返回积分值
int countTopEdgeIntegral(IplImage* img)
{
int black = 0;
for(int x=0; xwidth; x++)
{
unsigned char * par =(unsigned char *)(img-&imageData);
if(par[3*x] == 0 && par[3*x+1] == 0 && par[3*x+2] == 0)
black++;
}
return black;
}
// 接口:img,进行边界积分的图片
// 使用:对图片下侧的边界进行积分,返回积分值
int countBottomEdgeIntegral(IplImage* img)
{
int black = 0;
uchar * par = (uchar *)(img-&imageData+(img-&height-1)*(img-&widthStep));
for(int x=0 ;xwidth;x++)
{
if(par[3*x]==0 && par[3*x+1]==0 && par[3*x+2]==0)
black++;
}
return black;
}
// Created by 欧 长坤 on 13-9-16.
// Copyright (c) 2013年 欧 长坤 . All rights reserved.
// 接口:img_left,左图片
// img_right,右图片
// 使用:将右图片拼接到左图片的右侧,保存在一个新的IplImage中,并返回该地址。
IplImage * stitchingImage(IplImage* zuo,IplImage * you)
{
CvSize tt;
tt.width = zuo-&width+you-&width;
tt.height = zuo-&height;
IplImage* heti = cvCreateImage(tt,zuo-&depth,zuo-&nChannels);
for(int y = 0;yheight;y++)
{
uchar* par = (uchar*)(heti-&imageData+y*heti-&widthStep);
uchar* par2 = (uchar*)(zuo-&imageData+y*zuo-&widthStep);
uchar* par3 = (uchar*)(you-&imageData+y*you-&widthStep);
for(int x = 0;xwidth;x++)
{
par[3*x+0] = par2[3*x+0];
par[3*x+1] = par2[3*x+1];
par[3*x+2] = par2[3*x+2];
}
for(int x = 0;xwidth;x++)
{
par[3*zuo-&width+3*x+0] = par3[3*x+0];
par[3*zuo-&width+3*x+1] = par3[3*x+1];
par[3*zuo-&width+3*x+2] = par3[3*x+2];
}
}
return heti;
}
IplImage * stitchingImage_shangxia(IplImage* shang,IplImage * xia)
{
CvSize tt;
tt.height = shang-&height+xia-&height;
tt.width = shang-&width;
IplImage* heti = cvCreateImage(tt,shang-&depth,shang-&nChannels);
for(int y = 0;yheight;y++)
{
uchar* par = (uchar*)(heti-&imageData+y*heti-&widthStep);
uchar* parr = (uchar*)(heti-&imageData+(y+shang-&height)*heti-&widthStep);
uchar* par2 = (uchar*)(shang-&imageData+y*shang-&widthStep);
uchar* par3 = (uchar*)(xia-&imageData+y*xia-&widthStep);
for(int x = 0;xwidth;x++)
{
par[3*x+0] = par2[3*x+0];
par[3*x+1] = par2[3*x+1];
par[3*x+2] = par2[3*x+2];
}
}
for(int y = 0;yheight;y++)
{
uchar* par = (uchar*)(heti-&imageData+y*heti-&widthStep);
uchar* parr = (uchar*)(heti-&imageData+(y+shang-&height)*heti-&widthStep);
uchar* par2 = (uchar*)(shang-&imageData+y*shang-&widthStep);
uchar* par3 = (uchar*)(xia-&imageData+y*xia-&widthStep);
for(int x = 0;xwidth;x++)
{
parr[3*x+0] = par3[3*x+0];
parr[3*x+1] = par3[3*x+1];
parr[3*x+2] = par3[3*x+2];
}
}
return heti;
}
// Created by 黄 沐 on 13-9-16.
// Copyright (c) 2013年 黄 沐 . All rights reserved.
// 接口:filepath,图片文件的系统路径
// 使用:利用Add来访问到下一张图片路径
char* nextImageFilePath(char *filepath){
int i;
char *p;
char* save = filepath;
p=filepath;
for(i=0;filepath[i]!='\0';i++,p++){
if(filepath[i]<='9' && filepath[i+1]<='9' && filepath[i+2]<'9' && filepath[i+3]=='.'){
p++;
p++;
*p=*p+1;
return save;
}
if(filepath[i]<'9' && filepath[i+1]=='9' && filepath[i+2]=='.'){
*p=*p+1;
p++;
*p=*p-9;
return save;
}
if(filepath[i]=='9' && filepath[i+1]=='9' && filepath[i+2]=='.'){
p--;
*p=*p+1;
p++;
*p=*p-9;
p++;
*p=*p-9;
return save;
}
}
}
// 接口:filepath,图片文件的系统路径,num,修改文件名的数字
// 使用:传入num将filepath中最后的文件xxx.bmp改为num.bmp
char* numChangeFilePath(char *filepath, int num) {
int i,j;
int n_bai,n_shi,n_ge;
char* p;
char* save = filepath;
p=filepath;
n_ge=num%10;
n_bai=num/100;
if(num&=100)
{
n_shi=num/10;
n_shi=n_shi%10;
}
else n_shi=num/10;
for(i=0;filepath[i]!='\0';i++,p++){
if(filepath[i]<='9' && filepath[i+1]<='9' && filepath[i+3]=='.'){
*p='0';
for(j=0;j<n_bai;j++){
*p=*p+1;
}
p++;
*p='0';
for(j=0;j<n_shi;j++){
*p=*p+1;
}
p++;
*p='0';
for(j=0;j<n_ge;j++){
*p=*p+1;
}
return save;
}
}
}
// 接口:filepath,图片文件的系统路径,num,修改文件名的字符串
// 使用:仅适用于附件5
char* strChangeFilePath(char *filepath, char num[])
{
int i;
char* p;
char* save;
p=filepath;
save=filepath;
for(i=0;filepath[i]!='\0';i++,p++)
{
if(filepath[i]<='9' && filepath[i+1]<='9' && filepath[i+4]=='.')
{
*p='0';
*p=num[0];
p++;
*p='0';
*p=num[1];
p++;
*p='0';
*p=num[2];
p++;
*p='0';
*p=num[3];
return save;
}
}
return save;
}
[1]辛玉忠, 刘常凯, 判断两个信号相似性的内积表达式, 潍坊学院学报, Vol.2, No.4, 11-12, 2002
如果我的文章对你起到了帮助,你可以选择金额不限的捐助,帮助我写出更多的文章。