利用最小二乘法估算仿射变换参数

对opencv的理解还比较肤浅，所以会问一些比较弱智的问题。今天在做基于仿射变换模型图像配准的时候，H矩阵正确的计算出来，但是将目标图像进行H变换的时候，发现变换之后的图像出现了被裁减的情况。想来想去可能是目标图像乘以映射矩阵H之后出现负值，按0处理，于是想到了图像的平移….但是opencv中根本没有平移的函数，猛然间想到计算映射矩阵的时候其实已经用到了tx和ty这两个平移量，于是乎，对H矩阵的tx和ty进行操作之后，问题解决！

下图中原图像（左），转换后被裁减图像（中），调整后图像（右）

假设匹配点已经存储在matchpt矩阵中。

void CMTDoc::OnImageRegistration() 
{
	// TODO: Add your command handler code here
	if(RegistrateSet.DoModal()==IDOK)
	{
		int i,j;
		int m=matchpt->rows;
    //*********************计算H矩阵****************************	
		CvMat* A=cvCreateMat(m*2,1,CV_64FC1);
		CvMat* B=cvCreateMat(m*2,6,CV_64FC1);
		CvMat* tempH=cvCreateMat(6,1,CV_64FC1);
		CvMat* H=cvCreateMat(2,3,CV_64FC1);
		CvMat* temp=cvCreateMat(6,2*m,CV_64FC1);
		cvSetZero(B);
		for(i=0;i<m;i++)
			for(j=0;j<2;j++)
			{
				CV_MAT_ELEM(*B,double,i,j)=CV_MAT_ELEM(*matchpt,double,i,j+2);
				CV_MAT_ELEM(*B,double,i+m,j+2)=CV_MAT_ELEM(*matchpt,double,i,j+2);
			}
		for(i=0;i<m;i++)
		{
			CV_MAT_ELEM(*B,double,i,4)=1;
			CV_MAT_ELEM(*B,double,i,5)=0;
			CV_MAT_ELEM(*B,double,i+m,4)=0;
			CV_MAT_ELEM(*B,double,i+m,5)=1;
		}
		for(i=0;i<m;i++)
		{
			CV_MAT_ELEM(*A,double,i,0)=CV_MAT_ELEM(*matchpt,double,i,0);
			CV_MAT_ELEM(*A,double,i+m,0)=CV_MAT_ELEM(*matchpt,double,i,1);
		}
		cvInvert(B,temp,CV_SVD);
		cvMatMul(temp,A,tempH);

		CV_MAT_ELEM(*H,double,0,0)=CV_MAT_ELEM(*tempH,double,0,0);
		CV_MAT_ELEM(*H,double,0,1)=CV_MAT_ELEM(*tempH,double,1,0);
		CV_MAT_ELEM(*H,double,0,2)=CV_MAT_ELEM(*tempH,double,4,0);
		CV_MAT_ELEM(*H,double,1,0)=CV_MAT_ELEM(*tempH,double,2,0);
		CV_MAT_ELEM(*H,double,1,1)=CV_MAT_ELEM(*tempH,double,3,0);
//如出现图像的上方被裁减的情形，则可以对ty进行修正
		CV_MAT_ELEM(*H,double,1,2)=CV_MAT_ELEM(*tempH,double,5,0)+100;
		
 		cvSetImageROI(showImg1,rect2);
 		IplImage* tempimg=cvCreateImage(cvSize(rect2.width,rect2.height),IPL_DEPTH_8U,3);
 		cvCopy(showImg1,tempimg,NULL);

		cvSetImageROI(showImg1,cvRect(100,100,srcImg1->width+200,srcImg1->height+200));

		cvWarpAffine(tempimg,showImg1,H,CV_WARP_FILL_OUTLIERS);
		cvResetImageROI(showImg1);
		
		UpdateAllViews(NULL);
	}
}

 

功能：利用最小二乘法估算仿射变换参数
% [ ui vi ]' = [ a1 a2; a3 a4 ]*[ xi yi ]' + [ tx ty ]'
% A * x = B;
% A = [ x1  y1  0  0  1  0;
%          x2  y2  0  0  1  0;
%          :  :  :  ......…………;
%          ..................………;
%          xn  yn  0  0  1  0;
%          0  0  x1  y1  0  1;
%          0  0  x2  y2  0  1;
%          :  :  :  ......…………;
%          ..................………;
%          0  0  xn  yn  0  1;  ];
% B = [ x1'; x2';...; xn';y1';...; yn' ];
% H = [a1; a2; a3; a4; tx; ty];
% H = pinv(A)*B;
%%
clc;
clear;
% addpath E:\Master\Etone\SURFmex\examples\panorama  %添加你的数据所在路径
po=importdata('point_O.txt');%读取数据，引号内为文件名
pt=importdata('point_T.txt');

[m,n]=size(po);
A=zeros(2*m,2*n);
A(1:m,1:n)=po;
A(1+m:end,1+n:end)=po;
m2=size(A,1);
n2=size(A,2);
A(1:m2/2,(n2+1):(n2+2))=repmat([1 0],m2/2,1);
A(m2/2+1:end,(n2+1):(n2+2))=repmat([0 1],m2/2,1);

Bx=pt(:,1);
By=pt(:,2);
Bn=[Bx;By];

H=pinv(A)*Bn;
H2=reshape(H,2,3);//这句话有错…
disp(H);
disp(H2);

转载于:https://www.cnblogs.com/doucontorl/archive/2010/12/12/1903573.html