关于交叉累计剩余熵

经过好几天的修改，我自己编写的交叉累计剩余熵基本上算是完成了，现在上传与大家分享此程序是在pv_mi基础上改的，欢迎大家批评指出function out=pv_mi(tx,ty,ang,R,F)%pv插值求互信息[m,n]=size(R);hist=zeros(256,256);ha = zeros(1,256);hb = zeros(1,256);%归一化到256级灰度% if max(max(r))~=min(min(r)) %max(max(a))结果是A中最大的元素，max(A)结果一个行向量，元素分别是A的每个列向量的最大的元素% r = (r-min(min(r)))/(max(max(r))-min(min(r)));% else% r = zeros(M,N);% end% % if max(max(f))-min(min(f))% f = (f-min(min(f)))/(max(max(f))-min(min(f)));% else% f = zeros(M,N);% end% % r = double(int16(r*255))+1;% f = double(int16(f*255))+1;%MATLAB数组从1开始R=R+1;F=F+1;theta=ang*pi/180; %旋转角度转弧度cx=floor(n/2); %旋转的中心点cy=floor(m/2);% cx=125;% cy=149;%先求边缘概率密度,通过一副图像的直方图算,pv插值不能这样算% for i=1:m% for j=1:n% indexx = R(i,j);% indexy = F(i,j) ; % ha(indexx) = ha(indexx)+1;% hb(indexy) = hb(indexy)+1;% end% end%然后求联合概率密度for j=1:m for i=1:n %参考图像在浮动图像平移后的对应点 % x=i-tx; %列 % y=j-ty; % x=(i-cx)*cos(theta)-(j-cy)*sin(theta)+cx-tx; y=(i-cx)*sin(theta)+(j-cy)*cos(theta)+cy-ty; x1=floor(x); y1=floor(y); rval=R(j,i); %图像重合部分 if(x1>=1&&x1=1&&y1<m)% x1=floor(x); % y1=floor(y); a=x-x1; b=y-y1; w1=(1-a)*(1-b); w2=a*(1-b); w3=(1-a)*b; w4=a*b; fval1=F(y1,x1); fval2=F(y1,x1+1); fval3=F(y1+1,x1); fval4=F(y1+1,x1+1); hist(rval,fval1)=hist(rval,fval1)+w1; hist(rval,fval2)=hist(rval,fval2)+w2; hist(rval,fval3)=hist(rval,fval3)+w3; hist(rval,fval4)=hist(rval,fval4)+w4; % elseif (x1==n&&y1<m) %右边界% b=y-y1;% w1=1-b;% w2=b;% rval1=r(y1,x1);% rval3=r(y1+1,x1);% hist(fval,rval1)=hist(fval,rval1)+w1;% hist(fval,rval3)=hist(fval,rval3)+w2; % elseif (y1==m&&x1<n) %下边界% a=x-x1;% w1=1-a;% w2=a;% rval1=r(y1,x1);% rval2=r(y1,x1+1);% hist(fval,rval1)=hist(fval,rval1)+w1;% hist(fval,rval2)=hist(fval,rval2)+w2;% elseif(x1==n&&y1==m) %右下角点% rval=r(y1,x1);% hist(fval,rval)=hist(fval,rval)+1; end endendgab=zeros(256,256);for i=1:256 gab(i,:)=sum(hist(i:1:256,:));end%下面求边缘概率密度for i=1:256 ha(i)=sum(hist(i,:)); hb(i)=sum(hist(:,i));endhabsum = sum(sum(hist));hasum=sum(sum(ha));Gab=gab/habsum; %%联合分布函数概率Ha=ha/hasum; %%概率SP=sum(Gab,2);index1=find(Ha~=0);index2=find(SP~=0);GUV1=-sum(sum(SP(index2).*log(SP(index2))));len1=length(index1);Gab_N=zeros(256,len1);Gab_Ha=zeros(256,len1);for i=1:len1 Gab_N(:,i)=Gab(:,index1(i)); Gab_Ha(:,i)=Gab_N(:,i)/Ha(index1(i));endindex3=find(Gab_N~=0);GUV2=sum(sum(Gab_N(index3).*log(Gab_Ha(index3))));out=GUV1+GUV2;