图像相似性度量——互信息计算Matlab实现

约定不等于承诺〃 2022-06-12 05:27 275阅读 0赞

## 前言 ##

之前写了一篇关于计算互信息的[*文章*][Link 1]，基本的原理已经在那篇文章里面写出来了，这里就不在赘述了。在这片文章中将介绍采用的是直方图计算的方式，比之前的方法快了很多，软件的代码是Matlab实现的，这里贴出来与大家分享。

## 1. 代码 ##

%将图像B经过几何变换之后再计算互信息
    %x：水平方向上的位移
    %y：垂直方向上的位移
    %ang：图像的旋转角度
    %A：参考图像
    %B：浮动图像
    function [mi]=PV(x, y, ang, A, B)
    a=A;
    a=double(a);
    b=B;
    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
            coordinate_y(i,j)=j;
        end
    end
    [w z]=tforminv(tform,coordinate_x,coordinate_y);
    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_x-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
    hadsum=sum(sum(hab));
    index=find(hab~=0);
    pab=hab/hadsum;
    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))));
    fprintf('图片A/B的互信息为：%f, x=%f,y=%f,angle=%f \n', (Ha+Hb-Hab),x,y,ang);
    mi = Ha+Hb-Hab;

[Link 1]: http://blog.csdn.net/m_buddy/article/details/53579886