function [Km, rho, theta0, y,x, maxr, maxtranr] = phascorr(Im, Jm)
% This function performs phase correlation for rotation and then
% translation using Reddy and Chatterji method according to their article
% "An FFT-Based Technique for Translation, Rotation, and Scale-Invariant
% Image Registration", IEEE TRANSACTIONS ON IMAGE PROCESSING, VOL. 5,
% NO. 8, AUGUST 1996
% Input:
% Im - Base image. Must be square!
% Jm - Transformed image. Similiar to Im except for translation and
% rotation (scaling not implemented. Must be square!
% Im and Jm must have the same size.
% Output:
% Kmm - image of Jm rotated and translated to fit Im. Kmm and Im
% should be perfectly matched.
% offset_trans - cartesian coordinates of translation by pixels
% offset_rot - polar coordinates [radius, angle] of scaling and
% rotation in degrees
if size(Im,1) ~= size(Im,2)
warning('Base image is not square! Computation invalid.')
end
if size(Jm,1) ~= size(Jm,2)
warning('Transformed image is not square! Computation invalid.')
end
if size(Im) ~= size(Jm)
error('Images must be of the same size.')
end
FII = fft2(Im);
FI = log(abs(fftshift(FII)));
FJ = log(abs(fftshift(fft2(Jm))));
% Highpass Filter
% x=-0.5:1/croz:(0.5-1/croz);
% y=x;
% [X,Y]=meshgrid(x,y);
% HPFX = cos(pi*X).*cos(pi*Y);
% HPF = (1-HPFX).*(2-HPFX);
% FI = HPF.*FI;
% FJ = HPF.*FJ;
% H = fspecial('gaussian',15,1);
% FIfilt = imfilter(FI, H );
% FJfilt = imfilter(FJ, H);
% FI = FI - FIfilt;
% FJ = FJ - FJfilt;
% FI = FI(4:end-3,4:end-3);
% FJ = FJ(4:end-3,4:end-3);
% Building the polar coordinates
[m,n] = size(FI); c = round([m/2,n/2]);
rs = 1:m;
thetas = 0:2*pi/n:(2*pi-pi/n);
[thetas_img,rs_img] = meshgrid(thetas,rs);
[xs,ys] = pol2cart(thetas_img,rs_img);
zoff = complex(xs,ys) + complex(c(1),c(2));
% Converting the images to polar coordinates
FIpol = interp2(FI,real(zoff),imag(zoff),'linear', 0);
FJpol = interp2(FJ,real(zoff),imag(zoff),'linear', 0);
% polFI2 = interp2(FI,c(2)+rs_img.*cos(thetas_img),c(1)+rs_img.*sin(thetas_img),'nearest',0); %mean(FI(:)));
% Phase correlation for log-magnidute-polar image
rotR = fft2(FIpol).*conj(fft2(FJpol));
rotR = rotR./(abs(rotR)+eps);
rotr = abs(ifftshift(ifft2(rotR))); % 'symmetric' ?
[m,n] = size(FIpol); c = round([m/2,n/2]);
rotr(c(1)+1,c(2)+1)=0;
[rho, theta0] = find(rotr==max(rotr(:)));
rho = (rho-round(size(FIpol,1)/2)-1)/length(rs); % normalization of rs is unclear. relates to scaling.
theta0 = 360*(theta0-round(size(FIpol,2)/2)-1)/length(thetas);
maxr=max(rotr(:));
% Finding translation
Km0 = imrotate(Jm,-theta0, 'bilinear', 'crop');
[mK,nK] = size(Km0);
[mI,nI] = size(Im);
Km0 = imcrop(Km0, [ round((mK-mI)/2) round((nK-nI)/2) nI mI]);
FK = fft2(Km0);
tranR = FII.*conj(FK);
tranR = tranR./(abs(tranR)+eps);
tranr0 = abs(ifftshift(ifft2(tranR)));
maxtranr0 = max(tranr0(:));
% Checking ambiguity in theta0
Km1 = imrotate(Jm,-(180+theta0), 'bilinear', 'crop');
[mK,nK] = size(Km1);
Km1 = imcrop(Km1, [ round((mK-mI)/2) round((nK-nI)/2) nI mI]);
FK = fft2(Km1);
tranR = FII.*conj(FK);
tranR = tranR./(abs(tranR)+eps);
tranr1 = abs(ifftshift(ifft2(tranR)));
maxtranr1 = max(tranr1(:));
if maxtranr1 > maxtranr0
tranr = tranr1;
theta0 = theta0+180;
Km = Km1;
else
tranr = tranr0;
Km = Km0;
end
if theta0<0
theta0 = theta0+360;
end
% Finding Translation
[x,y] = find(tranr==max(tranr(:)));
x = -(x-round(mI/2)-1);
y = -(y-round(nI/2)-1);
if length(x)>1
x=0; y=0;
warning('Predication for translation is ambigious. Defaulting to no-translation.')
end
offset_trans = -[y,x];
Km = imtranslate(Km,-[y,x]);
maxtranr = max(tranr(:));
offset_rot = [rho, theta0];
end
评论1