在Matlab中实现图像的几何变换
图像的几何变换，包括平移、旋转、切变、缩放等规则的变换，还包括一些不规则的变换。主要的区别就体现在变换矩阵上。一般来说，当使用Homogeneous coordinates时，任何一个几何变换都可以用一个三阶矩阵h来表示。该矩阵有两类，一类对应于平移、旋转、切变、缩放等规则的变换，它的特点是第三行的第一二个元素是0，即h=
另一类称为Homogeneous矩阵，即h=
则对应于一些无法预测的变换。
可以看出，Homogeneous矩阵是最一般化的情形。
下面的问题是如何用Matlab来实现Homogeneous矩阵变换。在很多参考书上，讲到几何变换的时候都是用Matlab中的现成函数，如平移函数、旋转函数等，极少有直接运用矩阵的，本文则试图做这方面的尝试。
下面是一段Matlab代码：
function r=translate(H,I) %H是homography matrix,I是输入图像
A=[];
[m,n,z]=size(I);
for i=1:m
for j=1:n
N=H*[i;j;1];
x=round(N(1)/N(3));
y=round(N(2)/N(3));
if(x>0 & y>0) %N是变换后的坐标,N(1)是横坐标,N(2)是纵坐标
A(x,y,1)=I(i,j,1); %变换前后对应像素的值保持不变，red
A(x,y,2)=I(i,j,2); %green
A(x,y,3)=I(i,j,3); %blue
end
end
end
A=uint8(A);
imshow(A)
但是上述代码有个问题，那就是对变换后的坐标值用了取整，这样可能会导致图像中出现一些黑点，影响图像质量，一种较好的方法是利用Matlab中提供的interp2函数进行插值。部分代码如下：
X = H \[ Xp(:) Yp(:) ones(n,1) ]\';
xI = reshape( X(1,:)./X(3,:),wp,hp);
yI = reshape( X(2,:)./X(3,:),wp,hp);
Ip(:,:,1) = interp2(double(I(:,:,1)), xI, yI, \'bilinear\'); % red
Ip(:,:,2) = interp2(double(I(:,:,2)), xI, yI, \'bilinear\'); % green
Ip(:,:,3) = interp2(double(I(:,:,3)), xI, yI, \'bilinear\'); % blue
下面将展示一个具体例子：图像拼接。即把两张有部分相似的图片拼接成一张完整的图像。要做到这一点，首先要找出两张图片中相同点的对应坐标，可以利用Matlab提供的ginput函数来手动寻找（8组点）；随后，根据这8组坐标值来计算出Homogeneous矩阵H，利用的方法是最小二乘法；然后利用Homogeneous矩阵把一张图片进行变换后叠加到另一张图片上，拼接即告完成。下面是完整的代码：
%Getting correspondences
I=imread(\'e:\uttower1.jpg\');
I1=imread(\'e:\uttower2.jpg\');
figure(1);
imshow(I);
[x,y]=ginput(8); %先在I上取8个点
figure(2);
imshow(I1);
[xa,ya]=ginput(8); %再在I1上取8个对应的点
%Computing the homography parameters
x1=x(1);
x2=x(2);
x3=x(3);
x4=x(4);
x5=x(5);
x6=x(6);
x7=x(7);
x8=x(8);
x1a=xa(1);
x2a=xa(2);
x3a=xa(3);
x4a=xa(4);
x5a=xa(5);
x6a=xa(6);
x7a=xa(7);
x8a=xa(8);
y1=y(1);
y2=y(2);
y3=y(3);
y4=y(4);
y5=y(5);
y6=y(6);
y7=y(7);
y8=y(8);
y1a=ya(1);
y2a=ya(2);
y3a=ya(3);
y4a=ya(4);
y5a=ya(5);
y6a=ya(6);
y7a=ya(7);
y8a=ya(8);
A=[x1,y1,1,0,0,0,-x1a*x1,-x1a*y1; %利用最小二乘法求H矩阵
0,0,0,x1,y1,1,-y1a*x1,-y1a*y1;
x2,y2,1,0,0,0,-x2a*x2,-x2a*y2;
0,0,0,x2,y2,1,-y2a*x2,-y2a*y2;
x3,y3,1,0,0,0,-x3a*x3,-x3a*y3;
0,0,0,x3,y3,1,-y3a*x3,-y3a*y3;
x4,y4,1,0,0,0,-x4a*x4,-x4a*y4;
0,0,0,x4,y4,1,-y4a*x4,-y4a*y4;
x5,y5,1,0,0,0,-x5a*x5,-x5a*y5;
0,0,0,x5,y5,1,-y5a*x5,-y5a*y5;
x6,y6,1,0,0,0,-x6a*x6,-x6a*y6;
0,0,0,x6,y6,1,-y6a*x6,-y6a*y6;
x7,y7,1,0,0,0,-x7a*x7,-x7a*y7;
0,0,0,x7,y7,1,-y7a*x7,-y7a*y7;
x8,y8,1,0,0,0,-x8a*x8,-x8a*y8;
0,0,0,x8,y8,1,-y8a*x8,-y8a*y8];
b=[x1a;y1a;x2a;y2a;x3a;y3a;x4a;y4a;x5a;y5a;x6a;y6a;x7a;y7a;x8a;y8a];
h=A\b;
h=[h;1]; %H矩阵的最后一个元素是1
H=reshape(h,3,3)\';
%Warping between image planes
[h w d] = size(I);
cp = H*[ 1 1 w w ; 1 h 1 h ; 1 1 1 1 ];
Xpr = min([cp(1,:) 0]) : max([cp(1,:) w]); % min x : max x
Ypr = min([cp(2,:) 0]) : max([cp(2,:) h]); % min y : max y
[Xp,Yp] = meshgrid(Xpr,Ypr);
[wp hp] = size(Xp);
n = wp*hp;
X = H \[ Xp(:) Yp(:) ones(n,1) ]\'; %计算X=H\X\'
clear Ip;
xI = reshape( X(1,:)./X(3,:),wp,hp);
yI = reshape( X(2,:)./X(3,:),wp,hp);
Ip(:,:,1) = interp2(double(I(:,:,1)), xI, yI, \'bilinear\'); % red
Ip(:,:,2) = interp2(double(I(:,:,2)), xI, yI, \'bilinear\'); % green
Ip(:,:,3) = interp2(double(I(:,:,3)), xI, yI, \'bilinear\'); % blue
Ip=uint8(Ip);
figure(3);
imshow(Ip); %warping后的结果
%Create the output mosaic
hold on
imshow(I1); %最终结果