图像平移缩放旋转匹配

图像平移缩放旋转匹配

图像相位匹配，或者是傅里叶梅林变换匹配吧，能解决两幅图像之间的平移，缩放，旋转的匹配问题。研究了很久，弄清楚了里头的原理，才发现这些都是别人早就做好的东西。为什么我之前怎么查都查不到。

反正都研究了，不如认真点写成论文。老师说我的论文没有创新点，所以直接贴到这儿算了

csdn不支持粘贴图片我也是醉了。所以我把我的文档转成图片贴上来了，省得我再一个个保存图像复制图像 
 
 
 
 
 
 
 
 
 
 
 
 

下面我把我用于验证的代码贴上来：

这是只有平移的匹配函数： 
function [xoffset,yoffset,phase]=phase_match(img1,img2)%计算后一幅图相对于前一幅图的位移,原理为相位相关,输出为【行偏移，列偏移】 
%相位匹配好像不抗目标形变

if(size(img1,3)>1&&strcmp(class(img1),’uint8’))%输入的两幅图应该是一样大的 
img1=rgb2gray(img1); 
end 
if(size(img2,3)>1&&strcmp(class(img1),’uint8’)) 
img2=rgb2gray(img2); 
end 
img1=double(img1(:,:,1)); 
img2=double(img2(:,:,1)); 
img1f=fft2(img1); 
img2f=fft2(img2); 
%phasef=img2f./img1f; 
phasef=exp(log(img2f./abs(img2f))-log(img1f./abs(img1f)));%此处的提取相位有多种方法，只要原理对怎么都行，不过记得避免频谱图的有0值的地方容易发生除0错。。。 
phase=ifft2(phasef); 
[xoffset,yoffset]=find(abs(phase)==max(max(abs(phase)))); 
[a,b]=size(img1);

if(xoffset>a/2) 
xoffset=a+1-xoffset; 
xoffset=-xoffset; 
else 
xoffset=xoffset-1; 
end

if(yoffset>b/2) 
yoffset=b+1-yoffset; 
yoffset=-yoffset; 
else 
yoffset=yoffset-1; 
end 
end 
%本脚本已完

下面贴只有缩放的匹配，这里我直接给了个例子： 
img=imread(‘lena.jpg’); %自己搞个256*256的lena图，就跟我用的一模一样了 
img=rgb2gray(img); 
img=double(img); 
[a,b]=size(img);

img2=imresize(img(1:220,1:220),size(img)); 
disp(‘实际缩放比例’) 
256/220 
imgf=fft2(img); 
img2f=fft2(img2);%取绝对值能去掉平移因子 
window=hann(256)*hann(256)’; 
imgf=imgf.*window; 
img2f=img2f.*window;%抑制特别高的低频区 
imgf=imresize(imgf(1:end/2,1:end/2),[256,256]); 
img2f=imresize(img2f(1:end/2,1:end/2),[256,256]); 
[imgc,imgc_scalea,imgc_scaleb]=Cartesian2log_coordinate(abs(imgf));%直角坐标到对数坐标的转换函数 
[img2c,img2c_scalea,img2c_scaleb]=Cartesian2log_coordinate(abs(img2f));

window=hann(81)*hann(81)’; 
imshow([img,img2],[]) 
title(‘前一幅图 后一幅图’); 
[x,y,phase]=phase_match(img2c,imgc); 
disp(‘计算得到的缩放比例’) 
2^(x/imgc_scalea) 
%本脚本已完

要注意这里转对数坐标是直接对图像转对数坐标，后面还有个是对频谱图转对数坐标。因为图像时域缩放，频域放缩。。此处只有缩放没有平移和旋转的干扰，没必要去计算频谱图的缩放值。这些我在文章里应该说了

function [out,scale_a,scale_b]=Cartesian2log_coordinate(img)%输入一幅图，输出对数坐标图和行的缩放倍数和列的缩放倍数 
img=double(img); 
[a,b,c]=size(img); 
outa=floor(log2(a)); 
outb=floor(log2(b));%对数坐标图的实际长宽 
scale_a=a/outa; 
scale_b=b/outb;%将对数坐标图放大到和输入图一样大小 
out=zeros(size(img)); 
for j1=1:c 
for j2=1:a 
for j3=1:b %循环变换后的图 
tmp_outa=j2; 
tmp_outb=j3;%获取对数坐标图的坐标 
tmp_outa=tmp_outa/scale_a; 
tmp_outb=tmp_outb/scale_b; 
tmp_a=2^tmp_outa; 
tmp_b=2^tmp_outb;%将对数坐标图的坐标按照缩放比例找到放大后的坐标 
out(j2,j3,j1)=get_decimal_pixel(img,tmp_a,tmp_b,j1); 
end 
end 
end 
end 
%本脚本已完

function pixel=get_decimal_pixel(img,a,b,c)%使用双线性插值原理找到原图小数坐标的像素值 
pixel=0;

pixel=pixel+img(floor(a),floor(b),c) * (1-(a-floor(a))) * (1-(b-floor(b)));%对左上点的采样

if(a>floor(a)) 
pixel=pixel+img(floor(a)+1,floor(b),c) * (a-floor(a)) * (1-(b-floor(b)));%对左下点的采样 
end

if(b>floor(b)) 
pixel=pixel+img(floor(a),floor(b)+1,c) * (1-(a-floor(a))) * (b-floor(b));%对右上点的采样 
end

if(a>floor(a) && b>floor(b)) 
pixel=pixel+img(floor(a)+1,floor(b)+1,c) * (a-floor(a)) * (b-floor(b));%对右下点的采样 
end

end 
%本脚本已完

下面贴个缩放和平移组合后的例子： 
img=imread(‘lena.jpg’); 
img=rgb2gray(img); 
img=double(img); 
[a,b]=size(img);

window=hann(256)*hann(256)’; 
img2=imresize(img(1:220,1:220),size(img));%手动设置放大 
%img2=circshift(img2,[60,60]); 
img2=circshift(img2,[20,20]);%手动设置平移 
% 
% img2(1:20,:)=100*rand(size(img2(1:20,:)));%各种加噪声 
% img2(:,1:20)=100*rand(size(img2(:,1:20))); 
% img=img+10*rand(size(img)); 
% img=img+290; 
img=img.*window; 
img2=img2.*window; 
% 
imgf=fft2(img); 
img2f=fft2(img2); 
imgf=imgf.*window; 
img2f=img2f.*window; 
imgf=imresize(imgf(1:floor(end/2),1:floor(end/2)),[256,256]); 
img2f=imresize(img2f(1:floor(end/2),1:floor(end/2)),[256,256]); 
% imgf=imresize(rot90(rot90(imgf(floor(end/2):end,floor(end/2):end))),[256,256]); 
% img2f=imresize(rot90(rot90(img2f(floor(end/2):end,floor(end/2):end))),[256,256]);%这两句是验证后1/4的频谱也有缩放 
[imgc,imgc_scalea,imgc_scaleb]=Cartesian2log_coordinate(abs(imgf)); 
[img2c,img2c_scalea,img2c_scaleb]=Cartesian2log_coordinate(abs(img2f));

imshow([img,img2],[]) 
title(‘前一幅图 后一幅图’); 
[x,y,phase]=phase_match(img2c,imgc); 
figure; 
mesh(flipdim(abs(phase),1))

disp(‘实际缩放比 256/220’) 
256/220

disp(‘算得的缩放比’) 
2^(7/32) 
%这个代码只算了缩放值，因为缩放解决了，平移很好弄 
%本脚本已完

我的x和y我也记不清哪个代表行哪个代表列了，大家差不多看，嘿嘿~ 
下面贴个只有旋转的，，，，函数： 
function [x,y,phase]=rotate_match(img1,img2) 
if(size(img1,3)>1&&strcmp(class(img1),’uint8’))%输入的两幅图应该是一样大的 
img1=rgb2gray(img1); 
end 
if(size(img2,3)>1&&strcmp(class(img1),’uint8’)) 
img2=rgb2gray(img2); 
end 
img1=double(img1(:,:,1)); 
img2=double(img2(:,:,1)); 
img1=Cartesian2polar_coordinate(img1); 
img2=Cartesian2polar_coordinate(img2); 
[x,y,phase]=phase_match(img1,img2); 
end 
%本脚本已完

function [out,scale_a,scale_b]=Cartesian2polar_coordinate(img)%输入一幅图，输出对数坐标图和行的缩放倍数和列的缩放倍数 
img=double(img); 
[a,b,c]=size(img); 
%要保证转换前后的图像大小一致 
outa=sqrt(a*b);%极坐标下的行数代表原图的距原点最长的像素 
outb=b;%极坐标下的列数一共长90度，代表以坐标原点为基点的极坐标，图像只在旋转90度的范围内有。 
outc=c; 
out=zeros(outa,outb,outc); 
single_col=90/outb; %一列代表多少度 
for j1=1:outc 
for j2=1:outa 
for j3=1:outb %循环变换后的图 
tmp_outa=j2*cos(single_col*j3*pi/180); 
tmp_outb=j2*sin(single_col*j3*pi/180); 
out(j2,j3,j1)=get_decimal_pixel(img,tmp_outa,tmp_outb,j1); 
end 
end 
end 
end 
% 图： |———————| 
% | \ 角度 | 
% | \ | 
% | \ | 
% |———————| 
% 长度 
function pixel=get_decimal_pixel(img,a,b,c)%使用双线性插值原理找到原图小数坐标的像素值 
pixel=0; 
if(a<1) 
a=1; 
end 
if(b<1) 
b=1; 
end 
if(a>size(img,1)) 
a=size(img,1); 
end 
if(b>size(img,2)) 
b=size(img,2); 
end 
pixel=pixel+img(floor(a),floor(b),c) * (1-(a-floor(a))) * (1-(b-floor(b)));%对左上点的采样

if(a>floor(a)) 
pixel=pixel+img(floor(a)+1,floor(b),c) * (a-floor(a)) * (1-(b-floor(b)));%对左下点的采样 
end

if(b>floor(b)) 
pixel=pixel+img(floor(a),floor(b)+1,c) * (1-(a-floor(a))) * (b-floor(b));%对右上点的采样 
end

if(a>floor(a) && b>floor(b)) 
pixel=pixel+img(floor(a)+1,floor(b)+1,c) * (a-floor(a)) * (b-floor(b));%对右下点的采样 
end

end

%本脚本已完

下面是个只有旋转的例子： 
img=imread(‘lena.jpg’); 
img=rgb2gray(img); 
img=double(img); 
[a,b]=size(img); 
img2=zeros(a*2,b*2); 
img2(end-a+1:end,end-b+1:end)=img; 
img3=imrotate(img2,10); %旋转了10度 
[a2,b2]=size(img3); 
astart=floor(a2/2)+1; 
bstart=floor(b2/2)+1; 
img4=img3(astart:astart+255,bstart:bstart+255); %得到绕坐标原点旋转后的图像 
[x,y,phase]=rotate_match(img,img4); 
disp(‘实际的旋转度数’) 
10 
disp(‘计算得到的旋转度数’) 
(90/256)*y

%本脚本已完

function [x,y,phase]=rotate_match(img1,img2) 
if(size(img1,3)>1&&strcmp(class(img1),’uint8’))%输入的两幅图应该是一样大的 
img1=rgb2gray(img1); 
end 
if(size(img2,3)>1&&strcmp(class(img1),’uint8’)) 
img2=rgb2gray(img2); 
end 
img1=double(img1(:,:,1)); 
img2=double(img2(:,:,1)); 
img1=Cartesian2polar_coordinate(img1); 
img2=Cartesian2polar_coordinate(img2); 
[x,y,phase]=phase_match(img1,img2); 
end

function [out]=Cartesian2polar_coordinate(img)%输入一幅图，输出对数坐标图 
img=double(img); 
[a,b,c]=size(img); 
%要保证转换前后的图像大小一致 
outa=sqrt(a*b);%极坐标下的行数代表原图的距原点最长的像素 
outb=b;%极坐标下的列数一共长90度，代表以坐标原点为基点的极坐标，图像只在旋转90度的范围内有。 
outc=c; 
out=zeros(outa,outb,outc); 
single_col=90/outb; %一列代表多少度 
for j1=1:outc 
for j2=1:outa 
for j3=1:outb %循环变换后的图 
tmp_outa=j2*cos(single_col*j3*pi/180); 
tmp_outb=j2*sin(single_col*j3*pi/180); 
out(j2,j3,j1)=get_decimal_pixel(img,tmp_outa,tmp_outb,j1); 
end 
end 
end 
end 
% 图： |———————| 
% | \ 角度 | 
% | \ | 
% | \ | 
% |———————| 
% 长度 
function pixel=get_decimal_pixel(img,a,b,c)%使用双线性插值原理找到原图小数坐标的像素值 
pixel=0; 
if(a<1) 
a=1; 
end 
if(b<1) 
b=1; 
end 
if(a>size(img,1)) 
a=size(img,1); 
end 
if(b>size(img,2)) 
b=size(img,2); 
end 
pixel=pixel+img(floor(a),floor(b),c) * (1-(a-floor(a))) * (1-(b-floor(b)));%对左上点的采样

if(a>floor(a)) 
pixel=pixel+img(floor(a)+1,floor(b),c) * (a-floor(a)) * (1-(b-floor(b)));%对左下点的采样 
end

if(b>floor(b)) 
pixel=pixel+img(floor(a),floor(b)+1,c) * (1-(a-floor(a))) * (b-floor(b));%对右上点的采样 
end

if(a>floor(a) && b>floor(b)) 
pixel=pixel+img(floor(a)+1,floor(b)+1,c) * (a-floor(a)) * (b-floor(b));%对右下点的采样 
end

end

%本脚本已完

%下面是平移缩放旋转都有的例子，我只是算出来了缩放值和旋转值，剩下的工作就好解决了

img=imread(‘lena.jpg’); 
img=rgb2gray(img); 
img=double(img); 
img=imresize(img,1.1); 
img1=img(1:256,1:256); 
img2=img(21:276,21:276); %制造平移效果 平移10像素

[a,b]=size(img2); 
img3=100*rand(a*2,b*2); 
img3(end-a+1:end,end-b+1:end)=img2; 
img4=imrotate(img3,10); %旋转了 度 
[a2,b2]=size(img4); 
astart=floor(a2/2)+1; 
bstart=floor(b2/2)+1; 
img5=img4(astart:astart+255,bstart:bstart+255); %得到绕坐标原点旋转后的图像

img6=img5; 
img6=imresize(img6,1.1); %放大了1.1倍 
img6=img6(1:256,1:256);

% img1 和 img6 是具有平移缩放旋转关系的两幅图

img1f=fft2(img1); 
img2f=fft2(img6);

img1f=abs(img1f); 
img2f=abs(img2f); %取绝对值消除平移

window=hann(256)*hann(256)’; 
img1f=img1f.*window; 
img2f=img2f.*window; 
% img1f=log(img1f); 
% img2f=log(img2f);

img1f=imresize(img1f(1:floor(size(img1f,1)/2),1:floor(size(img1f,2)/2)),[256,256]); 
img2f=imresize(img2f(1:floor(size(img2f,1)/2),1:floor(size(img2f,2)/2)),[256,256]); %取1/4

img1c=Cartesian2polar_coordinate(img1f); 
img2c=Cartesian2polar_coordinate(img2f); %转为极坐标

img1c_2=polar2polarlog_coordinate(img1c); 
img2c_2=polar2polarlog_coordinate(img2c); %由极坐标转为对数极坐标

[xoffset,yoffset,phase]=phase_match(img1c_2,img2c_2);

%本脚本已完

function [out,scalea]=polar2polarlog_coordinate(img)%相当于只在行方向上做log，因为列方向是极坐标的角度参数 
[a,b,c]=size(img); 
la=log2(a); 
scalea=a/la;%计算出ρ上的缩放比 
out=zeros(size(img)); 
for j1=1:c 
for j2=1:a 
for j3=1:b 
tmp_outa=j2/scalea; 
tmp_outa=2^tmp_outa; 
out(j2,j3,j1)=get_decimal_pixel(img,tmp_outa,j3,j1); 
end 
end 
end 
end

function pixel=get_decimal_pixel(img,a,b,c)%使用双线性插值原理找到原图小数坐标的像素值 
pixel=0; 
if(a<1) 
a=1; 
end 
if(b<1) 
b=1; 
end 
if(a>size(img,1)) 
a=size(img,1); 
end 
if(b>size(img,2)) 
b=size(img,2); 
end 
pixel=pixel+img(floor(a),floor(b),c) * (1-(a-floor(a))) * (1-(b-floor(b)));%对左上点的采样

if(a>floor(a)) 
pixel=pixel+img(floor(a)+1,floor(b),c) * (a-floor(a)) * (1-(b-floor(b)));%对左下点的采样 
end

if(b>floor(b)) 
pixel=pixel+img(floor(a),floor(b)+1,c) * (1-(a-floor(a))) * (b-floor(b));%对右上点的采样 
end

if(a>floor(a) && b>floor(b)) 
pixel=pixel+img(floor(a)+1,floor(b)+1,c) * (a-floor(a)) * (b-floor(b));%对右下点的采样 
end

end

哦对，中间有一个点需要注意，转坐标的时候，从直角坐标转到对数坐标再转到极坐标是不可以等价到对数极坐标的~~，不懂的可以查一查对数极坐标的定义

我也把我的代码打包上传到我的CSDN资源里头了，代码有可能会有点乱，，已经开始忘了不想整理了，不过大部分都是能跑的。

需要注意的一点是，希望大家还是懂了原理之后再跑代码，不然代码结果可能看不懂哦，毕竟我这么懒，没写多少说明