fouriermellin代码.rar_victoryg71_傅里叶变换_傅里叶梅林_图像处理_梅林变换

function RegnisterFourierMellin() I1 = imread('待配准图像需自己设置相应图像.bmp'); I2 = imread('lena_cropped_rotated_shifted.bmp'); SizeX = size(I1, 1); SizeY = size(I1, 2); FA = fftshift(fft2(I1)); FB = fftshift(fft2(I2)); IA = hipass_filter(size(I1, 1),size(I1,2)).*abs(FA); IB = hipass_filter(size(I2, 1),size(I2,2)).*abs(FB); L1 = transformImage(IA, SizeX, SizeY, SizeX, SizeY, 'nearest', size(IA) / 2, 'valid'); L2 = transformImage(IB, SizeX, SizeY, SizeX, SizeY, 'nearest', size(IB) / 2, 'valid'); THETA_F1 = fft2(L1); THETA_F2 = fft2(L2); a1 = angle(THETA_F1); a2 = angle(THETA_F2); THETA_CROSS = exp(i * (a1 - a2)); THETA_PHASE = real(ifft2(THETA_CROSS)); THETA_SORTED = sort(THETA_PHASE(:)); SI = length(THETA_SORTED):-1:(length(THETA_SORTED)); [THETA_X, THETA_Y] = find(THETA_PHASE == THETA_SORTED(SI)); DPP = 360 / size(THETA_PHASE, 2); Theta = DPP * (THETA_Y - 1); R1 = imrotate(I2, -Theta, 'nearest', 'crop'); R2 = imrotate(I2,-(Theta + 180), 'nearest', 'crop'); R1_F2 = fftshift(fft2(R1)); a1 = angle(FA); a2 = angle(R1_F2); R1_F2_CROSS = exp(i * (a1 - a2)); R1_F2_PHASE = real(ifft2(R1_F2_CROSS)); R2_F2 = fftshift(fft2(R2)); a1 = angle(FA); a2 = angle(R2_F2); R2_F2_CROSS = exp(i * (a1 - a2)); R2_F2_PHASE = real(ifft2(R2_F2_CROSS)); MAX_R1_F2 = max(max(R1_F2_PHASE)); MAX_R2_F2 = max(max(R2_F2_PHASE)); if (MAX_R1_F2 > MAX_R2_F2) [y, x] = find(R1_F2_PHASE == max(max(R1_F2_PHASE))); R = R1; else [y, x] = find(R2_F2_PHASE == max(max(R2_F2_PHASE))); if (Theta < 180) Theta = Theta + 180; else Theta = Theta - 180; end R = R2; end Tx = x - 1; Ty = y - 1; if (x > (size(I1, 1) / 2)) Tx = Tx - size(I1, 1); end if (y > (size(I1, 2) / 2)) Ty = Ty - size(I1, 2); end input2_rectified = R; move_ht = Ty; move_wd = Tx; total_height = max(size(I1,1),(abs(move_ht)+size(input2_rectified,1))); total_width = max(size(I1,2),(abs(move_wd)+size(input2_rectified,2))); combImage = zeros(total_height,total_width); registered1 = zeros(total_height,total_width); registered2 = zeros(total_height,total_width); if((move_ht>=0)&&(move_wd>=0)) registered1(1:size(I1,1),1:size(I1,2)) = I1; registered2((1+move_ht):(move_ht+size(input2_rectified,1)),(1+move_wd):(move_wd+size(input2_rectified,2))) = input2_rectified; elseif ((move_ht<0)&&(move_wd<0)) registered2(1:size(input2_rectified,1),1:size(input2_rectified,2)) = input2_rectified; registered1((1+abs(move_ht)):(abs(move_ht)+size(I1,1)),(1+abs(move_wd)):(abs(move_wd)+size(I1,2))) = I1; elseif ((move_ht>=0)&&(move_wd<0)) registered2((move_ht+1):(move_ht+size(input2_rectified,1)),1:size(input2_rectified,2)) = input2_rectified; registered1(1:size(I1,1),(abs(move_wd)+1):(abs(move_wd)+size(I1,2))) = I1; elseif ((move_ht<0)&&(move_wd>=0)) registered1((abs(move_ht)+1):(abs(move_ht)+size(I1,1)),1:size(I1,2)) = I1; registered2(1:size(input2_rectified,1),(move_wd+1):(move_wd+size(input2_rectified,2))) = input2_rectified; end if sum(sum(registered1==0)) > sum(sum(registered2==0)) plant = registered1; bleed = registered2; else plant = registered2; bleed = registered1; end combImage = plant; for p=1:total_height for q=1:total_width if (combImage(p,q)==0) combImage(p,q) = bleed(p,q); end end end imshow(combImage, [0 255]); function [r,g,b] = transformImage(A, Ar, Ac, Nrho, Ntheta, Method, Center, Shape) global rho; theta = linspace(0,2*pi,Ntheta+1); theta(end) = []; switch Shape case 'full' corners = [1 1;Ar 1;Ar Ac;1 Ac]; d = max(sqrt(sum((repmat(Center(:)',4,1)-corners).^2,2))); case 'valid' d = min([Ac-Center(1) Center(1)-1 Ar-Center(2) Center(2)-1]); end minScale = 1; rho = logspace(log10(minScale),log10(d),Nrho)'; xx = rho*cos(theta) + Center(1); yy = rho*sin(theta) + Center(2); if nargout==3 if strcmp(Method,'nearest'), r=interp2(A(:,:,1),xx,yy,'nearest'); g=interp2(A(:,:,2),xx,yy,'nearest'); b=interp2(A(:,:,3),xx,yy,'nearest'); elseif strcmp(Method,'bilinear'), r=interp2(A(:,:,1),xx,yy,'linear'); g=interp2(A(:,:,2),xx,yy,'linear'); b=interp2(A(:,:,3),xx,yy,'linear'); elseif strcmp(Method,'bicubic'), r=interp2(A(:,:,1),xx,yy,'cubic'); g=interp2(A(:,:,2),xx,yy,'cubic'); b=interp2(A(:,:,3),xx,yy,'cubic'); else error(['Unknown interpolation method: ',method]); end mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1); r(mask)=0; g(mask)=0; b(mask)=0; else if strcmp(Method,'nearest'), r=interp2(A,xx,yy,'nearest'); elseif strcmp(Method,'bilinear'), r=interp2(A,xx,yy,'linear'); elseif strcmp(Method,'bicubic'), r=interp2(A,xx,yy,'cubic'); else error(['Unknown interpolation method: ',method]); end mask= (xx>Ac) | (xx<1) | (yy>Ar) | (yy<1); r(mask)=0; end function H = hipass_filter(ht,wd) res_ht = 1 / (ht-1); res_wd = 1 / (wd-1); eta = cos(pi*(-0.5:res_ht:0.5)); neta = cos(pi*(-0.5:res_wd:0.5)); X = eta'*neta; H=(1.0-X).*(2.0-X);