图像的简单相位配准 with MATLAB

相位配准概述
相位配准是一种非常经典的理论,理论论述互联网上有很多,我直接给出一个流程图,简述重点。

算法流程图

  • 相位配准的原始理论互联网上已经很多,我不过多赘述,使用完全相同的两个图像在空间内平移进行推导。只强调一下,在实际使用时,可以证明,只要两个图像有完全相同的部分,那么就可以使用相位配准
    结论图像

  • 带有旋转/缩放的相位配准不要看普通教材,普通教材愚蠢的试图先绕图像中心旋转缩放,再进行平移。事实上,只要有空间变换(或刚体运动学)知识,就可以知道如果刚体可以仅通过缩放和旋转从A变为B,那么必然可以确定一个空间点O,使刚体A仅绕O进行一次缩放和一次旋转到达B,也就是说,只要我们确定了点O,那么进行对数极换元就可以直接确定旋转和缩放参数。

  • 如果只是绕中心旋转的图片,那么做对数极换元,再做平移的相位配准即可,我必须说,这毫无意义,这里面的工作最大头为对数极换元,which和图像处理没有什么关系。
  • MATLAB 2014a后自带了一个 imregcor 函数,它对重合区域较大(目测至少20%)的两个图像进行相位配准非常高效,平移,旋转 缩放都OK。本文末尾附带的是一个我的实现,它的重难点在于附带了一个图像拼接程序,要比相位配准中求出坐标更复杂。

function runrun()
PhCorrelationRotation('cameraman.png','cameraman_part_rotate_scaling.png',0)
end

function PhCorrelationRotation(baseimgpath,sampleimgpath,deltaH)
if nargin<3
    deltaH=0;
end
fixed=imread(baseimgpath);
moving=imread(sampleimgpath);
if(ndims(fixed)~=ndims(moving))
    if(ndims(fixed)==3)
        fixed=rgb2gray(fixed);
    end
    if(ndims(moving)==3)
        moving=rgb2gray(moving);
    end
end
if deltaH==0
    tformEstimate = imregcorr(moving,fixed);
    movingReg = imwarp(moving,tformEstimate);
else
    movingReg=moving;
end
PhCorrelation(fixed,movingReg,moving)
end

function PhCorrelation(baseImg,sampleImg,Origin)
%calc
[deltaDim1,deltaDim2,R]=phaseCalc(baseImg,sampleImg);
combineImg=stichImg(baseImg,sampleImg,deltaDim1,deltaDim2);
imwrite(combineImg,'imans.bmp')
%show ans
figure
subplot(2,2,1);
imshow(baseImg);
title('图像1')
subplot(2,2,2);
imshow(Origin);
title('图像2')
subplot(2,2,3);
imshow(combineImg);
title('拼接图像')
subplot(2,2,4);
mesh(flipud(abs(R)));
title('脉冲图像')
end

function [indI,indJ,R]=phaseCalc(im1,im2)

if(ndims(im1)==3)
    im1=rgb2gray(im1);
    fprintf('n非8bit图像,可能会丢失信息n')
end
if(ndims(im2)==3)
    im2=rgb2gray(im2);
    fprintf('n非8bit图像,可能会丢失信息n')
end
maxsize=max([size(im1);size(im2)]);
imm1=zeros(maxsize);
imm2=zeros(maxsize);

imm1(1:size(im1,1),1:size(im1,2),:)=im2double(im1);
imm2(1:size(im2,1),1:size(im2,2),:)=im2double(im2);
F1=fft2(imm1);
F2=fft2(imm2);

div=F1.*conj(F2)./abs(F1.*conj(F2));
R=ifft2(div);
maxvalue=max(max(abs(R)));
[indI,indJ]=find(abs(R)==maxvalue);
if(max(max(R))==1) %完全相关,即表明两幅图像完全一样
    indJ=0;
    indI=0;
else
    if (indI==size(im1,1))
        indI=0;
    end
    if(indJ==size(im1,2))
        indJ=0;
    end
end
end

function combineImg=stichImg(baseImg,sampleImg,deltaDim1,deltaDim2)
maxsize=max([size(baseImg);size(sampleImg)]);
szbase=size(baseImg);szsample=size(sampleImg);
windosize=min([szbase(1:2)-[deltaDim1,deltaDim2];szsample(1:2)]);%edgeofbase=1 downedge 2 rightedge
if any(windosize<0)
    combineImg=stichLandR(sampleImg,baseImg,maxsize(1)-deltaDim1,maxsize(2)-deltaDim2);
else
    samplearea=sampleImg(1:windosize(1),1:windosize(2),:);
    [d1,d2,~]=phaseCalc(baseImg,samplearea);
    if norm([d1,d2]-[deltaDim1,deltaDim2])>10;
        combineImg=stichLandR(sampleImg,baseImg,maxsize(1)-deltaDim1,maxsize(2)-deltaDim2);
    else
        combineImg=stichLandR(baseImg,sampleImg,deltaDim1,deltaDim2);
    end
end
end

function combineImg=stichLandR(leftImg,rightImg,deltaDim1,deltaDim2)
if ndims(leftImg)==3
    maxsize=max([size(leftImg);size(rightImg)+[deltaDim1,deltaDim2,0]]);
else
    maxsize=max([size(leftImg);size(rightImg)+[deltaDim1,deltaDim2]]);
end

combineImg=zeros(maxsize);
combineImg(1:(size(leftImg,1)),1:(size(leftImg,2)),:)=leftImg;
mid=combineImg(deltaDim1:((size(rightImg,1)+deltaDim1-1)),deltaDim2:((size(rightImg,2)+deltaDim2)-1),:)+double(rightImg);
combineImg(deltaDim1:((size(rightImg,1)+deltaDim1-1)),deltaDim2:((size(rightImg,2)+deltaDim2)-1),:)=mid;
combineImg=uint8(combineImg);
end