rigid-registration.rar

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %几何变换及互信息计算 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [mi]=PV(x,y,ang,I,J) a=double(I); b=double(J); [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; coordinate_y(i,j)=j; end end [w z]=tforminv(tform,coordinate_x,coordinate_y); f_New=uint8(zeros(M,N)); % for i2=1:M % for j2=1:N % source_x=w(i2,j2); % source_y=z(i2,j2); % if (source_x>=width-1||source_y>=height-1||double(uint16(source_x))<=0||double(uint16(source_y))<=0) % f_New(i2,j2)=0; % else % if ((source_x/double(uint16(source_x))==1.0)&(source_y/double(uint16(source_y))==1.0)) % f_New(i2,j2)=b(int16(source_x),int16(source_y)); % else %双线性插值法 % p=double(uint16(source_x)); % q=double(uint16(source_y)); % x11=double(b(p,q)); % x12=double(b(p,q)); % x21=double(b(p+1,q)); % x22=double(b(p+1,q+1)); % f_New(i2,j2)=uint8((q+1-source_y)*((source_x-p)*x21+(p+1-source_x)*x11)+(source_y-q)*((source_x-p)*x22+(p+1-source_x)*x12)); % end % end % end % end % f_New=double(f_New); % a=(a-min(min(a)))/(max(max(a))-min(min(a))); % f_New=(f_New-min(min(f_New)))/(max(max(f_New))-min(min(f_New))); % a=double(int16(a*255))+1; % f_New=double(int16(f_New*255))+1; % for i1=1:M % for j1=1:N % index_x=a(i1,j1); % index_y=f_New(i1,j1); % hab(index_x,index_y)=hab(index_x,index_y)+1; % end % end 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 habsum=sum(sum(hab)); index=find(hab~=0); pab=hab/habsum; Hab=sum(sum(-pab(index).*log2(pab(index)))); pa=sum(pab'); index=find(pa~=0); Ha=sum(sum(-pa(index).*log2(pa(index)))); pb=sum(pab); index=find(pb~=0); Hb=sum(sum(-pb(index).*log2(pb(index)))); mi=Ha+Hb-Hab end