基于互信息和归一化互信息的医学图像配准算法matlab仿真-摩杜云开发者社区

1.算法仿真效果 matlab2022a仿真结果如下：

2.算法涉及理论知识概要信息论中将互信息定义为信息之间的关系，可以表示为两个随机变量之间统计相关性的度量，由此可以得出图像互信息的计算方法。作为图像多模态配准中的度量，图像互信息利用对图像灰度值的统计数据形成单个图像的灰度值概率函数和两个图像相似部分对应的灰度值联合概率函数，以此来衡量两幅图像的相关程度。在图像配准的过程中，认为两幅图像的相关性最大时对应的互信息也最大。基于互信息的图像配准最早在1995年被用于医学图像配准中，此后研究人员对于互信息的配准方法做了大量研究。

2.1互信息互信息(mutual information)来自于信息论，是用于图像配准的一种经典相似性测度。它最初是由Collignon等和Viola等于1995年引入图像处理领域，用于配准多模图像。目前由于其无需预处理、自动化程度高以及鲁棒性强等特点，利用最大互信息法进行多模图像配准成为了图像处理领域的研究热点。因此，互信息的计算成为多模图像配准的一个关键研究问题，在工业、生物医学、计算机视觉、遥感等应用领域中将可能发挥更大的作用。目前已有许多基于互信息的图像配准算法的研究文献，所采用的互信息计算方法也各有不同。互信息的计算方法是关系到配准精度和效率的关键因素，尚未有文献对此进行系统全面的讨论；由于算法模型的多样性，使人们在实际应用常常很难确定采用哪一种计算方法。

一.熵简介

   熵，热力学中表征物质状态的参量之一，用符号S表示，其物理意义是体系混乱程度的度量。任何一种能量在空间中分布得越均匀，熵就越大，一个体系的能量完全均匀分布时，这个系统的熵就达到最大值。它在控制论、概率论、数论、天体物理、信息论等领域都有重要应用，在不同的学科中也有引申出的更为具体的定义，是各领域十分重要的参量。

1.状态函数

  熵S是状态函数，具有加和（容量）性质，是广度量非守恒量，因为其定义式中的热量与物质的量成正比，但确定的状态有确定量。其变化量ΔS只决定于体系的始终态而与过程可逆与否无关。由于体系熵的变化值等于可逆过程热温商δQ/T之和，所以只能通过可逆过程求的体系的熵变。孤立体系的可逆变化或绝热可逆变化过程ΔS=0。

2.宏观量

  熵是宏观量，是构成体系的大量微观离子集体表现出来的性质。它包括分子的平动、振动、转动、电子运动及核自旋运动所贡献的熵，谈论个别微观粒子的熵无意义。

3.绝对值

  熵的绝对值不能由热力学第二定律确定。可根据量热数据由第三定律确定熵的绝对值，叫规定熵或量热法。还可由分子的微观结构数据用统计热力学的方法计算出熵的绝对值，叫统计熵或光谱熵。

4.信息熵

  任何一个消息的自信息量都代表不了信源所包含的平均自信息量。不能作为整个信源的信息测度，因此定义自信息量的数学期望为离散信源的平均自信息量：

   称之为信源的信息熵。H是从整个信源的统计特性来考虑的，它是从平均意义上来表征信源的总体特性的。对于某特定的信源，其信息熵只有一个;不同的信源因统计特性不同，其熵也不同。 

   图像处理中常常需要比较两幅图像的相似度，例如在图像配准中，将互信息值作为配准结果的评价指标。互信息值是其中一种较为常用的方法，其核心思想是熵，即图像所包含的信息。假设有图像A,B，则它们互信息值计算公式为：I(A,B)=H(A)+H(B)-H(A,B)公式的意义很明显。H(A,B)为A,B的联合熵，是使用A,B的联合直方图计算出的结果，可以理解为A,B共同包含的信息。联合熵H(A,B)是检测随机变量A和B相关性的统计量。对于两个随机变量A、B，它们的概率分布分别为pA(a)和pB(b)，联合分布为pAB(a,b).则他们的联合熵为：

若A、B共同包含信息越少(A、B信息重复越多)，则H(A,B)越小，因此互信息值I越大。

2.2归一化互信息尽管互信息测度成功地应用于图像配准中，由于两幅图像重叠都分的大小对互信息的量度有很大影响，重叠部分减小，参与统计互信息的像素个数减小导致互信息值减小，互信息与两个图像重叠部分多少成正比，误配数量增加可能导致互信息值增大。因此互信息值达到最大并不能保证得到正确的配准结果。为了解决这个问题，使目标函数能更加准确反映互信息量和配准参数之间的关系。Studholme等提出了一个归一化互信息测度，归一化互信息使配准函数更平滑，它能减少对图像重叠部分敏感性，配准精度更高。归一化互信息的定义如下：

3.MATLAB核心程序

    l = 1;
    if x <= mr
        for y = wzy1-rfield+wzy2:step:wzy1+rfield+wzy2
            if y <= nr
                Im_F = image2(mf-x+1:mf,nf-y+1:nf);
                Im_R = image1(1:x,1:y);
            else
                Im_F = image2(mf-x+1:mf,1:nf+nr+1-y);
                Im_R = image1(1:x,y-nf+1:nr);
            end
            if method == 'MI'
                MeV(k,l) = mi(Im_R,Im_F,Mi_method);
                l = l+1;
            end
        end
    else
        for y = wzy1-rfield+wzy2:step:wzy1+rfield+wzy2
            if y <= nr
                Im_F = image2(1:mf+mr+1-x,nf-y+1:nf);
                Im_R = image1(x-mf+1:mr,1:y);
            else
                Im_F = image2(1:mf+mr+1-x,1:nf+nr+1-y);
                Im_R = image1(x-mf+1:mr,y-nf+1:nr);
            end
            if method == 'MI'
                MeV(k,l) = mi(Im_R,Im_F,Mi_method);
                l = l+1;
            end
        end
    end
    k = k+1;
end
....................................................
figure
mesh(x,y,MeV);
title('测度曲线');
[Maxm,ind] = max(MeV(:));
[X,Y] = ind2sub(size(MeV),ind);
X = (X-1)*step;Y = (Y-1)*step; % 平移像素个数
 
% 模板图像分别沿x，y轴平移(X-1)*step，(Y-1)*step像素点后与参考图像配准
% 此时，浮动图像中心位于参考图像坐标系(wzx1-rfield+(X-1)*step,wzy1-rfield+(Y-1)*step)处
H.Position = [132 258 260 260];
figure(H)
imagesc(image2)
title('原始浮动图像');
colormap(gray)
H.Position = [402 258 260 260];
figure(H)
imagesc(image1)
title('原始参考图像')
colormap(gray)
h1 = get(gca,'Position');
x = h1(1);y = h1(2);w = h1(3);h  = h1(4);
X1 = wzx1-rfield+X;Y1 = wzy1-rfield+Y;
H.Position = [672 258 260 260];
figure(H)
axes('Position',[w/mr*(X1-(mf-1)/2)+x,h/nr*(Y1-(nf-1)/2)+y,w/mr*mf,h/nr*nf]);
imagesc(image2)
title('配准图像')
colormap (gray)