【图像配准】基于互信息的图像配准算法：MI、EMI、ECC算法

原文：http://blog.csdn.net/hujingshuang/article/details/47910949

简介：

        基于互信息的图像配准算法以其较高的配准精度和广泛的适用性而成为图像配准领域研究的热点之一，而基于互信息的医学图像配准方法被认为是最好的配准方法之一。基于此，本文将介绍简单的基于互信息的图像配准算法。

预备知识

熵

        熵（entropy）是信息论中的重要概念，用来描述系统不确定性的测度，反映一个系统本身所能提供的信息总量。除去枯燥的概念，信息熵的数学表达式为：



当然，对于一幅图像来说，其熵的计算表达式如下：



       hi表示图像Y中灰度值为i的像素点总数，N表示图像Y的灰度级数。显然Pi表示灰度i出现的概率，于是很自然的就会想到用直方图来计算。（对于单幅图像MATLAB中可由entropy()函数求得）

联合熵

联合熵反映了随机变量X、Y的相关性，设两个随机变量X和Y，则X和Y的联合信息熵表示为：



对于两幅图像X、Y来说，利用联合直方图，显然可以计算出二者的联合熵。

互信息

         互信息（Mutual Information，MI）是信息论中的重要概念，描述了两个系统之间的相关性，或互相包含信息的多少。在图像配准中，两幅图的互信息是通过它们的熵以及联合熵，来反映它们之间信息的相互包含成都程度（02-21-2017更正）。对于图像R、F来说，其互信息表示为：



当两幅图像相似度越高或重合部分越大时，其想关性也越大，联合熵越小，也即互信息越大。

配准及改进：

        通过上述预备知识的了解，可以得到一个匹配准则：寻找模板与各子图之间互信息（MI）的最大者，即为配准图像。

        当模板与源图中的多个子图具有相似的灰度分布式，互信息容易出现误匹配的情况；互信息对两幅图像之间的重叠区域比较敏感，如果两幅图像的重叠区太小，互信息就会很小，配准精度随之降低。基于以上问题，Studholme和Maes分别对互信息测度进行了改进，分别提出了归一化互信息（Normalization
Mutual Information，NMI）、熵相关系数（Entropy Corrleation Coefficient，ECC）改进后的表达式如下：



        改进后的NMI、ECC算法，依然是寻找模板与子图之间最大NMI、ECC值作为配准图像。根据上面的三种匹配准则，我们编写代码，实验如下。

实例验证：

MATLAB代码：

[cpp]
view plain
copy

print?

%//*****************************************       
%//Copyright (c) 2015 Jingshuang Hu       
       
%//@filename:demo.m       
%//@datetime:2015.08.23       
%//@author:HJS       
%//@e-mail:eleftheria@163.com       
%//@blog:http://blog.csdn.net/hujingshuang       
%//*****************************************      
%%   
%//--------------------------------------------------------------------------  
%//图像互信息配准  
clear all;  
close all;  
img=rgb2gray(imread('lena.jpg'));%//源图  
[M N]=siz
c55a
e(img);  
%//--------------------------------------------------------------------------  
x0=30;y0=40;  
dx=64;dy=64;%//模板、子图尺寸  
img1=img(x0:x0+dx-1,y0:y0+dy-1);%//模板  
ET=entropy(img1);%//模板熵  
%//--------------------------------------------------------------------------  
%//联合熵  
[m n]=size(img1);%//模板尺寸  
MI=zeros(M-dx,N-dy);%//互信息  
NMI=zeros(M-dx,N-dy);%//归一化互信息  
ECC=zeros(M-dx,N-dy);%//熵相关系数  
for i=1:M-dx  
    for j=1:N-dy  
        img2=img(i:i+dx-1,j:j+dy-1);%//子图  
        ES=entropy(img2);%//模板熵  
        histq=zeros(256,256);%//联合直方图，清空  
        %//联合直方图  
        for s=1:m  
            for t=1:n  
                x=img1(s,t)+1;y=img2(s,t)+1;%//灰度<—>坐标  
                histq(x,y)=histq(x,y)+1;  
            end  
        end  
        p=histq./sum(sum(histq));%//联合概率密度  
        EST=-sum(sum(p.*log(p+eps)));%//联合熵（越小说明相似度越高）  
        MI(i,j)=ES+ET-EST;%//MI互信息越大，说明相互包含的信息多，即越匹配  
        NMI(i,j)=(ES+ET)/EST;%//NMI,越大越匹配  
        ECC(i,j)=2*MI(i,j)/(ES+ET);%//ECC，越大越匹配  
    end  
end  
%//--------------------------------------------------------------------------  
mi_max=max(max(MI));  
nmi_max=max(max(NMI));  
ncc_max=max(max(ECC));  
[xt1,yt1]=find(MI==mi_max);  
[xt2,yt2]=find(NMI==nmi_max);  
[xt3,yt3]=find(ECC==ncc_max);  
src=img1;  
dst1=img(xt1:xt1+dx-1,yt1:yt1+dx-1);  
dst2=img(xt2:xt2+dx-1,yt2:yt2+dx-1);  
dst3=img(xt3:xt3+dx-1,yt3:yt3+dx-1);  
figure;imshow(src);title('模板');  
figure;imshow(img);hold on;rectangle('position',[yt1,xt1,n-1,m-1],'edgecolor','r');title('MI配准图');hold off;  
figure;imshow(img);hold on;rectangle('position',[yt2,xt2,n-1,m-1],'edgecolor','r');title('NMI配准图');hold off;  
figure;imshow(img);hold on;rectangle('position',[yt3,xt3,n-1,m-1],'edgecolor','r');title('NCC配准图');hold off;  
%//over  

%//*****************************************
%//Copyright (c) 2015 Jingshuang Hu

%//@filename:demo.m
%//@datetime:2015.08.23
%//@author:HJS
%//@e-mail:eleftheria@163.com
%//@blog:http://blog.csdn.net/hujingshuang
%//*****************************************
%%
%//--------------------------------------------------------------------------
%//图像互信息配准
clear all;
close all;
img=rgb2gray(imread('lena.jpg'));%//源图
[M N]=size(img);
%//--------------------------------------------------------------------------
x0=30;y0=40;
dx=64;dy=64;%//模板、子图尺寸
img1=img(x0:x0+dx-1,y0:y0+dy-1);%//模板
ET=entropy(img1);%//模板熵
%//--------------------------------------------------------------------------
%//联合熵
[m n]=size(img1);%//模板尺寸
MI=zeros(M-dx,N-dy);%//互信息
NMI=zeros(M-dx,N-dy);%//归一化互信息
ECC=zeros(M-dx,N-dy);%//熵相关系数
for i=1:M-dx
for j=1:N-dy
img2=img(i:i+dx-1,j:j+dy-1);%//子图
ES=entropy(img2);%//模板熵
histq=zeros(256,256);%//联合直方图，清空
%//联合直方图
for s=1:m
for t=1:n
x=img1(s,t)+1;y=img2(s,t)+1;%//灰度<—>坐标
histq(x,y)=histq(x,y)+1;
end
end
p=histq./sum(sum(histq));%//联合概率密度
EST=-sum(sum(p.*log(p+eps)));%//联合熵（越小说明相似度越高）
MI(i,j)=ES+ET-EST;%//MI互信息越大，说明相互包含的信息多，即越匹配
NMI(i,j)=(ES+ET)/EST;%//NMI,越大越匹配
ECC(i,j)=2*MI(i,j)/(ES+ET);%//ECC，越大越匹配
end
end
%//--------------------------------------------------------------------------
mi_max=max(max(MI));
nmi_max=max(max(NMI));
ncc_max=max(max(ECC));
[xt1,yt1]=find(MI==mi_max);
[xt2,yt2]=find(NMI==nmi_max);
[xt3,yt3]=find(ECC==ncc_max);
src=img1;
dst1=img(xt1:xt1+dx-1,yt1:yt1+dx-1);
dst2=img(xt2:xt2+dx-1,yt2:yt2+dx-1);
dst3=img(xt3:xt3+dx-1,yt3:yt3+dx-1);
figure;imshow(src);title('模板');
figure;imshow(img);hold on;rectangle('position',[yt1,xt1,n-1,m-1],'edgecolor','r');title('MI配准图');hold off;
figure;imshow(img);hold on;rectangle('position',[yt2,xt2,n-1,m-1],'edgecolor','r');title('NMI配准图');hold off;
figure;imshow(img);hold on;rectangle('position',[yt3,xt3,n-1,m-1],'edgecolor','r');title('NCC配准图');hold off;
%//over

输出结果：

 

 

 

参考文献：

1、别术林,基于互信息的医学图像配准算法研究[D],2014.

2、陈伟卿,基于互信息的医学图像刚性配准研究[D],2010.