www.gusucode.com > 《图像配准技术及其Matlab编程实现》--源码程序 > 《图像配准技术及其Matlab编程实现》/match/PartialVolumeInterpolation.m

    function [hrf]=PartialVolumeInterpolation(x,y,ang,imgR,imgF)
%参数：imgR：参考图像
%      imgF：浮动图像
%      x：水平平移（像素）
%      y：垂址平移（像素）
%      ang：旋转角度（角度）
%      hrf：联合直方图
%功能：利用PV插值法计算两幅图像的联合直方图的值，结果返回联合直方图数组的值

a=imgR;
a=double(a);
b=imgF;
b=double(b);
[M,N]=size(a);
hab=zeros(256,256);
ha=zeros(1,256);
hb=zeros(1,256);

if max(max(a))~=min(min(a))
    a=(a-min(min(a)))/(max(max(a))-min(min(a)));
else
    a=zeros(M,N);
end

if max(max(b))-min(min(b))
    b=(b-min(min(b)))/(max(max(b))-min(min(b)));
else
    b=zeros(M,N);
end
a=double(int16(a*255))+1;
b=double(int16(b*255))+1;

%仿射空间几何变换
[width,height]=size(b);
u=(width-1)/2;
v=(height-1)/2;
rad=pi/180*ang;
t1=[1 0 0;0 1 0;x y 1];
t2=[1 0 0;0 1 0;-u -v 1];
t3=[cos(rad) -sin(rad) 0;sin(rad) cos(rad) 0;0 0 1];
t4=[1 0 0;0 1 0;u v 1];
T=t2*t3*t4*t1;
tform=maketform('affine',T);
coordinate_x=zeros(width,height);
coordinate_y=zeros(width,height);
for i=1:width
    for j=1:height
        coordinate_x(i,j)=i;
    end
end
for i=1:width
    for j=1:height
        corrdinate_y(i,j)=j;
    end
end
[w z]=tforminv(tform,coordinate_x,coordinate_y);

%利用PV插值法统计联合直方图
for i=1:width
    for j=1:height
        source_x=w(i,j);
        source_y=z(i,j);
        if(source_x>width-1 || source_y>height-1 || double(uint16(source_x))<=1 || double(uint16(source_y))<=1)
            hab(a(1,1),a(1,1))=hab(a(1,1),a(1,1))+1;
        else
            m=fix(source_x);
            n=fix(source_y);
            index_b=b(i,j);
            index_a0=a(m,n);
            index_a1=a(m+1,n);
            index_a2=a(m,n+1);
            index_a3=a(m+1,n+1);
            dx=source_x-m;
            dy=source_y-n;
            hab(index_a0,index_b)=hab(index_a0,index_b)+(1-dx)*(1-dy);
            hab(index_a1,index_b)=hab(index_a1,index_b)+dx*(1-dy);
            hab(index_a2,index_b)=hab(index_a2,index_b)+(1-dx)*dy;
            hab(index_a3,index_b)=hab(index_a3,index_b)+dx*dy;
        end
    end
end
hrf=hab;
figure,imshow(hrf);