Phase Correlation_matlab_dirtyacc_phasecorrelation_

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