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

Published at发布于:： 2013-09-21 | PV/UV: / | Reading阅读:： 7 min

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。

信号相似性算法：

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

图像$i$为需要查找后（或上、下、前）继图像$j$位于指定的待选空间（集合）中，提取图像$i$与图像$j$的边际信号序列（若是寻找后继，则图像$i$提取右侧信号，图像$j$提取左侧信号，若是寻找前继，则图像$i$提取左侧信号，图像$j$提取右侧信号，寻找上继和下继可类推），为了防止计算时导致的数据溢出，将信号序列归一化（将$[0,255]$上的整型值转化为$[0,1]$上的双精度浮点值）；

计算两者的规范相关系数$R_{xy}$；

初始化筛选下界$temp = 0.5$和标记变量$flag = 0$；

若$R_{xy}>temp$，则$temp = R_{xy}$，$flag = j$，$j++$；

若遍历完所有图片，则第$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提供了非常方便的访问图像数据的结构，更容易统计该列某颜色值像素个数。其算法描述为：

对于阈值$\Phi_{1}$（可根据实际结果进行调整）将图像二值化，保证图像中仅有黑白两色。

对于图像$i$（不妨设为当前行的首张图片），计算它在边际$right$上的边际积分，记为$p$。

在剩余图像中取一图像$j$，计算它在$left$边际上的边际积分，记为$q$。

计算的值$\Vert p-q \Vert$，若小于给定阈值$\Phi_{2}$并保存$j$的值在数组$space$中，作为一个可能后继。

若未遍历完所有图片，则返回（3），否则进入（6）。

数组$space$保存的值即为筛选出的可能匹配对象的集合。

部分关键代码

// 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

推荐订阅源

Posts on Changkun's Blog