图像处理基本算法--仿射变换

淩亂°似流年 2022-06-16 10:17 283阅读 0赞

几何空间变换和图像配准

几何空间变换又称为橡皮膜变换，因为他可以看做是在一幅橡皮膜上印制图像，然后根据一定规则拉伸橡皮膜。由两个基本操作组成：

1）坐标的空间变换

2）灰度内插

最常用的是仿射变换一般形式如下：

\[x,y,1\] = \[v,w,1\]\*T

\[t11   t12   0\]

= \[v,w,1\]\*\[t21   t22  0\]

\[t31   t32   0\]

这一变换可以根据矩阵T中元素的值，对一组坐标点做尺度变换，旋转，平移，偏移。

![0_1325689236ysZ9.gif][]

我们有两种方法使用该式子，第一种是前向映射，它由扫描输入图像的像素，并在每一个位置（v,w）直接计算输出图像中像素的空间位置（x,y）组成。

第二种方法是反向映射，扫描输出图像的位置，并在每一个位置（x,y）使用(v,w) = T-1(x,y)反向计算输入图像相应的位置，然后进行灰度内插（一般使用最邻近，双线性，双三次内插）。

前向映射一个问题是，输入图像的的两个或多个像素会映射到输出图像的同一个位置，另一种可能是某些输出位置完全没有像素，因此反向映射比前向映射更加有效，被很多的商业公司选用，包括mathworks的 matlab使用的也是反向映射。

一、图像平移

x = v + tx

y = w + ty

原图：

![0_1325689345716a.gif][]

平移后的图像：

![0_1325689393KiAK.gif][]

二、图像旋转

x = v\*cosA  - w \* sinA

y = v \*sinA + w \* cosA

进行图像旋转时，需要注意的是当图像的大小发生改变时此时若以左上角点为原点，那么原点会发生改变，一般会选定图像的中心作为原点，这样在图像进行坐标变换时中心不变。但在在变换时要记得加减一个常数，使原点相对应。

旋转后的图像：

![0_1325689495TJuM.gif][]

观察旋转后的图像，会有很明显的锯齿现象。

若采用双线性插值或者三次立方插值那么锯齿依然难以消除。

三、图像的尺度变换

x  = cx \* v

y = cy \* w

这个就不附上图了，比较简单

四、图像的偏移

垂直偏移变换：

x =v\*sv + w

y = w

水平偏移变换：

x = v

y = sh\*v+w

水平偏移后的图像：

![0_1325690871sFTe.gif][]

仿射变换的源代码：

#include<cv.h>
    #include<highgui.h>
    #include <cmath>
    
    #include <stdlib.h>
    const double PI = 3.141592653; 
    int main(){
    	IplImage * image,*image2;
    	image = cvLoadImage("E:\\image\\base.jpg",0);
    	cvNamedWindow("image",CV_WINDOW_AUTOSIZE);
    	//cvSaveImage("E:\\image\\pollen.jpg",image,0);
    	cvShowImage("image",image);
    	//cvWaitKey(0);
    	//平移 dx dy  
    	int dx = 10;
    	int dy = 20;
    	image2 = cvCreateImage(cvSize(image->width + dx, image->height + dy),image->depth,1);
    	unsigned char * ptr,*dst;
    	int i,j ,m,n;//i j为新图中坐标，m n 为原图中坐标
    	
    		for(i = 0 ; i < image2->height;i++){
    			for(j = 0; j< image2->width;j++){
    				dst = (unsigned char *)image2->imageData + i*image2->widthStep + j; 
    				n = j - dx;
    				m = i - dy;
    				if(m < 0 || n <  0 || m >= image->height  || n >= image->width){
    					*dst = 0;
    				}
    				else {
    					ptr = (unsigned char * )(image->imageData + m * image->width + n);
    					*dst =  *ptr; 
    				}
    
    
    			}
    		}		
    	cvNamedWindow("image2",1);
    	cvShowImage("image2",image2);
    	//cvWaitKey(0);
    	cvSaveImage("E:\\image\\basepingyi.jpg",image,0);
    
    	/**************************************************************************************/
    	//旋转图像 这里以逆时针为正  插值选用双线性插值 围绕左上角点进行旋转  双线性插值
    	IplImage* image3;
    	double  angle = 15;
    	double Hangle = angle*PI/180*(-1);
    	int Width,Height;
    	//double r,s;
    	double  r,s;
    
    	int temp[4];
    	double z1,z2;
    	//原点会发生移动
    	double  fx,fy;
    
    	double Hcos = cos(Hangle);
    	double Hsin = sin(Hangle);
    
    	//x y 存放原图中四个点的位置，以中心为原点 左上角开始 顺时针数
    	int x[4];
    	int y[4];
    	int x1[4];
    	int y1[4];
    	x[0] = -(image->width-1)/2;
    	x[1] = -x[0];
    	x[2] = -x[0];
    	x[3] = x[0];
    	y[0] = -(image->height -1)/2;
    	y[1] = y[0];
    	y[2] = -y[0];
    	y[3] = -y[0];
    	//x1 y1 分别存放新图中图像的四个角点的位置 以中心为原点 左上角点开始 顺时针数
    	for (i = 0; i < 4 ; i++)
    	{
    		x1[i] = (int)(x[i]*Hcos - y[i]*Hsin + 0.5);
    		y1[i] = (int)(x[i]*Hsin + y[i]*Hcos + 0.5);
    		printf("x:   %d \n",x[i] );
    		printf("y:   %d\n",y[i]);
    		printf("x1:   %d \n",x1[i] );
    		printf("y1:   %d\n",y1[i]);
    		
    	}
    	//确定新图像的长宽
    	if (abs(y1[2] - y1[0]) > abs(y1[3] - y1[1]))
    	{
    		Height = abs(y1[2] - y1[0]);
    		Width  = abs(x1[3] - x1[1]);
    	}
    	else{
    		Height = abs(y1[3] - y1[1]);
    		Width =  abs( x1[2] - x1[0]);
    	}
    	image3 = cvCreateImage(cvSize(Width,Height),image->depth,1);
    	//两个偏移常量
    	fx = -1*(Width -1)* Hcos*0.5 -(Height -1)*Hsin*0.5 + (image->width-1)/2;
    	fy =  (Width -1)*Hsin*0.5 -(Height -1)*Hcos*0.5 + (image->height -1)/2; 
    	
    	for (i = 0 ; i< Height ; i++)
    	{
    		for (j =0 ; j< Width; j++)
    		{
    			dst = (unsigned char *)(image3->imageData + i*image3->widthStep + j);
    			
    			r =   (int)(i*Hcos- j * Hsin + fy + 0.5);
    			s =   (int)(i*Hsin  + j*Hcos + fx + 0.5 );
    			 
    			if (r <0 || s < 0 || r >= image->height || s >= image->width)
    			{
    				*dst = 0;
    			}
    			else{
    				ptr = (unsigned char *)(image->imageData + image->widthStep*(int)r + (int)s );
    				temp[0] = *ptr;
    				temp[1] = *(ptr+1);
    				ptr = (unsigned char*)(image->imageData + image->widthStep *(int)(r+1) + int (s));
    				temp[2] = *ptr;
    				temp[3] = *(ptr+1);
    
    				z1 = (temp[1] - temp[0])*(s - int(s)) + temp[0];
    				z2 = (temp[3] -temp[2])*(s- int(s)) + temp[2];
    				*dst = (int)(z1 + (z2 -z1)*(r- (int)r));
    			}
    			
    		}
    	}
    	cvSaveImage("bldpingyi.bmp",image2);
    	/**************************************************************************************/
    
    	/**************************************************************************************/
    	/* //这里是采用最简单的取离映射点最近的那个点 另外 的方法是采用双线性插值的方式 
    	//旋转图像 这里以逆时针为正  插值选用双线性插值 围绕左上角点进行旋转 
    	IplImage* image3;
    	double  angle = 15;
    	double Hangle = angle*PI/180*(-1);
    	int Width,Height;
    	//double r,s;
    	int r,s;
    	//原点会发生移动
    	double  fx,fy;
    
    	double Hcos = cos(Hangle);
    	double Hsin = sin(Hangle);
    
    	//x y 存放原图中四个点的位置，以中心为原点 左上角开始 顺时针数
    	int x[4];
    	int y[4];
    	int x1[4];
    	int y1[4];
    	x[0] = -(image->width-1)/2;
    	x[1] = -x[0];
    	x[2] = -x[0];
    	x[3] = x[0];
    	y[0] = -(image->height -1)/2;
    	y[1] = y[0];
    	y[2] = -y[0];
    	y[3] = -y[0];
    	//x1 y1 分别存放新图中图像的四个角点的位置 以中心为原点 左上角点开始 顺时针数
    	for (i = 0; i < 4 ; i++)
    	{
    		x1[i] = (int)(x[i]*Hcos - y[i]*Hsin + 0.5);
    		y1[i] = (int)(x[i]*Hsin + y[i]*Hcos + 0.5);
    		printf("x:   %d \n",x[i] );
    		printf("y:   %d\n",y[i]);
    		printf("x1:   %d \n",x1[i] );
    		printf("y1:   %d\n",y1[i]);
    
    	}
    	//确定新图像的长宽
    	if (abs(y1[2] - y1[0]) > abs(y1[3] - y1[1]))
    	{
    		Height = abs(y1[2] - y1[0]);
    		Width  = abs(x1[3] - x1[1]);
    	}
    	else{
    		Height = abs(y1[3] - y1[1]);
    		Width =  abs( x1[2] - x1[0]);
    	}
    	image3 = cvCreateImage(cvSize(Width,Height),image->depth,1);
    	//两个偏移常量
    	fx = -1*(Width -1)* Hcos*0.5 -(Height -1)*Hsin*0.5 + (image->width-1)/2;
    	fy =  (Width -1)*Hsin*0.5 -(Height -1)*Hcos*0.5 + (image->height -1)/2; 
    
    	for (i = 0 ; i< Height ; i++)
    	{
    		for (j =0 ; j< Width; j++)
    		{
    			dst = (unsigned char *)(image3->imageData + i*image3->widthStep + j);
    
    			r =   (int)(i*Hcos- j * Hsin + fy + 0.5);
    			s =   (int)(i*Hsin  + j*Hcos + fx + 0.5 );
    
    			if (r <0 || s < 0 || r >= image->height || s >= image->width)
    			{
    				*dst = 0;
    			}
    			else{
    				ptr = (unsigned char *)(image->imageData + image->widthStep*r + s );
    				*dst = *ptr;
    			}
    
    		}
    	}
    	*/
    	/**************************************************************************************/
    	cvNamedWindow("image3",1);
    	cvShowImage("image3",image3);
    	cvSaveImage("bldxuanzhuan2.bmp",image3);
    	//cvWaitKey(0);
    
    	//垂直偏移  x = v*Sv + w    y = w
    	//水平偏移  x = v    y = Sh*v + w   
    	//这里实现水平偏移 
    	double sh = 0.8;
    	//确定新图的长和宽
    	Width = (int)(image->width + sh* image->height);
    	Height = image->height;
    	int t;
    	IplImage * image4 = cvCreateImage(cvSize(Width , Height),image->depth,1);
    	for (i = 0 ; i< Height ; i++)
    	{
    		for (j =0 ; j< Width; j++)
    		{
    			dst = (unsigned char *)(image4->imageData + i*image4->widthStep + j);
    
    			//r =   (int)(i*Hcos- j * Hsin + fy + 0.5);
    			//s =   (int)(i*Hsin  + j*Hcos + fx + 0.5 );
    			t = i;
    			s =  j - sh*i;
    
    			if (t <0 || s < 0 || t >= image->height || s >= image->width)
    			{
    				*dst = 0;
    			}
    			else{
    				ptr = (unsigned char *)(image->imageData + image->widthStep*t + (int)s );
    				z1 = *ptr;
    				z2 = *(ptr+1);
    				*dst = (int)(z1 + (z2 -z1)*(s- (int)s));
    			}
    
    		}
    	}
    	cvNamedWindow("image4",1);
    	cvShowImage("image4",image4);
    	cvSaveImage("bldpianyi.bmp",image4);
    	cvWaitKey(0);
    
    	cvReleaseImage(&image);
    	cvReleaseImage(&image2);
    	cvReleaseImage(&image3);
    	cvDestroyAllWindows();
    	return 0;
    }

//图像的仿射变换和图像配准
    //平移、缩放、旋转、偏移 
    //若想变换之后的图像时可以恢复的，那么在进行变换的时候图像的大小应该可以变化，是为了保留更多的信息
    //使用反向映射的方式

五、图片配准

图像配准，一般用于对齐两幅或多幅相同场景的图像。前面已经讲了所有的仿射变换的函数形式，更为复杂的是多种仿射变化的连续作用，显然可以得出，多种仿射变换的连续作用可以用如下式子概括：

x  =c1\*v + c2\*w + c3\*v\*w + c4

y  =c5\*v + c6\*w + c7\*v\*w + c8

方程组含有8个未知量，假设我们已知一对图像中相对应的4个点的坐标，那个既可以解出8个未知量，从而实现两张图片的配准。

使用photoshop打开一种图片，放上至少4个小方块，能够肉眼识别，但又不要太大，太大会使误差增大，对图像进行一定的仿射变换，并记下四对匹配点的坐标。

这里

image1中取四个点

47 44

433 39

47 430

462 444

image2中取四个点

140 40

461 106

55 305

395 390

采用matlab解矩阵，求出8个参数，然后根据该参数进行图像配准：

matlab程序：\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*\*

% 图像配准 首先是解方程组获取参数 
    d = ones(4,4);
    d(1:2,:) = [47 433 47 462;
                44, 39 ,430 ,444];
    d(3,:) = d(1,:).*d(2,:)
    
    x= [140 461 55 395]
    C1 = x/d
    %得 C1 =  0.8290   -0.2200   -0.0000  110.7284
    %事实上这个图像当时因为只是进行了缩放和旋转，因此不存在偏移，C（3） = 0是合理的
    y = [40 ;106 ;305 ;390];
    C2 = y/d
    %C2 = 0.1797    0.6863    0.0000    1.3466
    %理由同上，事实上这个图像当时因为只是进行了缩放和旋转，因此不存在偏移，C（3） = 0是合理的

图像配准程序，也是基于反向映射

/*******************************************************************
    已知两幅图像上的四个对应点，很据变换式，采用matlab求出二者之间的变换方程
    根据变换式，假如需要将图片2变换到图片1，那么需要求出图片1到2 的变换式，反向映射，
    假如映射区域在2区域之内，那么赋予该处的值，否者需要赋予0 
    ************************************************************************/
    //中值滤波和均值滤波
    #include<cv.h>
    #include<highgui.h>
    
    int main(){
    	IplImage * image,*image2,*image3;
    	image = cvLoadImage("E:\\image\\base.jpg",0);//以灰度图像的形式读入图片
    	image2 = cvLoadImage("E:\\image\\base2.jpg",0);
    	cvNamedWindow("image",CV_WINDOW_AUTOSIZE);
    	cvNamedWindow("image2",CV_WINDOW_AUTOSIZE);
    	cvNamedWindow("image3",CV_WINDOW_AUTOSIZE);
    	cvShowImage("image",image);
    	image3 = cvCreateImage(cvGetSize(image),image->depth,1);
    	//cvWaitKey(0);
    	unsigned char * ptr,*dst;
    	int i,j;
    	double m, n;
    	int temp[4];
    	double C[8] = {0.8290, -0.2200, -0.0000,  110.7284,0.1797,0.6863,0.0000,1.3466};//图片1到2 的变换矩阵
    	double Z1,Z2;
    	for( i = 0 ; i < image3->height;i++){
    		for( j = 0; j< image3->width;j++){
    			dst = (unsigned char *)image3->imageData+ i*image3->widthStep+ j;
    			n = C[0]*j + C[1]*i + C[2]*i*j + C[3];
    			m = C[4]*j + C[5]*i + C[6]*i*j + C[7];
    			if (m < 0.0 || n < 0.0 || m >= image2->height || n >= image2->width)
    			{
    				*dst = 0;
    			}
    			else{
    				ptr = (unsigned char *)(image2->imageData + image2->widthStep * (int)m + (int)n);
    				temp[0] = *ptr;
    				temp[1] = *(ptr+1);
    				ptr = (unsigned char *)(image2->imageData + image2->widthStep * (int)(m +1) + (int)n);
    				temp[2] = *ptr;
    				temp[3] = *(ptr+1);
    				Z1 = (temp[1] - temp[0])*(m - int(m)) + temp[0];
    				Z2 = (temp[3] -temp[2])*(m- int(m)) + temp[2];
    				*dst = (int)(Z1 + (Z2 -Z1)*(n- (int)n));
    			}
    			
    		}
    	}
    	
    
    	cvShowImage("image2",image2);
    	cvShowImage("image3",image3);
    	cvSaveImage("bldpeizhun.bmp",image3);
    
    
    	//比较配准效果 
    	IplImage* image4;
    	image4 = cvCreateImage(cvSize(image->width , image->height),image->depth , 1);
    	cvSub(image , image3 , image4,0);
    	cvNamedWindow("sub");
    	cvShowImage("sub", image4);
    	cvSaveImage("bldsub.bmp",image4);
    
    
    	cvWaitKey(0);
    	
    	return 0;
    }

配准的两个源图像：

图1：

![0_1325690246nM1z.gif][]

图2

![0_1325690265HzSY.gif][]

配准之后的图像：

![0_1325690365Lmnu.gif][]

配准图像和原图像的差值：

![0_13256903962u7u.gif][]

[0_1325689236ysZ9.gif]: /images/20220616/c0c2472ecbb64e3dbbbf81ec3d48c0c8.png
[0_1325689345716a.gif]: /images/20220616/062ad7e617e94405a183823786a34cd1.png
[0_1325689393KiAK.gif]: /images/20220616/155be96c23ce4cb6b94d4c15148c4697.png
[0_1325689495TJuM.gif]: /images/20220616/6f6b7ce7eb8645849c12fb93b5eb55d7.png
[0_1325690871sFTe.gif]: /images/20220616/d53d01bca9944849b45a157fdee1f6b7.png
[0_1325690246nM1z.gif]: /images/20220616/3cd1876313ea4a349f3c65c3459c8260.png
[0_1325690265HzSY.gif]: /images/20220616/7f4e983e403a40cf8bdd8eebfde8b158.png
[0_1325690365Lmnu.gif]: /images/20220616/6d0bdf7c91fe41c1ae42767136a951b0.png
[0_13256903962u7u.gif]: /images/20220616/d1b9300874494f688d952486c27df17b.png