基于互信息的图像配准方法

- 6 -
http://www.nwpu.edu.cn
Image Registration Based on Mutual Information
Wang Guoxing
College of Automation, Northwestern Polytechnical University , Xi 'an, PRC(710129)
Hy=0; for i=1:c; if y_marg(i)~=0 Hy = Hy + -(y_marg(i)*(log2(y_marg(i)))); end end Hx=0; for i=1:r; if x_marg(i)~=0 Hx = Hx + -(x_marg(i)*(log2(x_marg(i)))); end end hxy = -sum(sum(b.*(log2(b+(b==0))))); f = -(Hx+Hy-hxy); end
- 5 -
http://www.nwpu.edu.cn for ii=1:rows; for jj=1:cols; h(IM1(ii,jj)+1,IM 2(ii,jj)+1)= h(IM1(ii,jj)+1,IM 2(ii,jj)+1)+1; end end [r,c] = size(h); b= h./(r*c); %正常化的联合直方图 y_marg=sum(b); x_marg=sum(b'); %正常化联合直方图的行的和 %正常化联合直方图的列的和
西北工业大学自动化学院，陕西西安(710129) E-mail: wangguoxing@mail.nwpu.edu.cn
摘 要： 本文介绍了图像配准的基本原理与方法，然后提出了一种基于互信息 的图像配准方法，并利用 MATLAB 对两幅图像进行了配准实验，通过实验，证 明该方法非常有效，达到了图像配准的目的。 关键词：MATLAB; 信息熵；互信息；图像配准； 中图分类号：TN91 文献标识码：A
http://www.nwpu.edu.cn
基于互信息的图像配准方法
学院： 专业： 班级： 学号： 姓名：
自动化学院 信息工程专业 09030701 2007302147 王国星
时间：2009 年 11 月 1 日
- 1 -
http://www.nwpu.edu.cn
基于互信息的图像配准方法
王国星
%图像1的边际熵
%图像2的边际熵
%联合熵
%互信息的相反数，因为主函数要求f的最小值，所以f最小时，互信息最大。
4. 结论
本文介绍了图像配准的基本原理与方法， 然后提出了一种基于互信息的图像 配准方法，并利用 MATLAB 对两幅图像进行了图像配准实验，通过实验，证明 该方法非常有效，达到了图像配准的目的。
1. 引 言
随着医学、计算机技术及生物工程技术的发展，医学影像学为临床诊断 提供了多种模态的医学图像，不同的医学图像提供了相关脏器的不同信息： CT(Computed Tomography，电子计算机 X 射线断层扫描 )和 MRI(Magnetic resona nce ima ging， 核磁共振成像 ) 以较高的空间分辨率提供了脏器的解剖结 构信息。在实际临床应用中，单一模态的图像往往不能提供医生所需要的 足够的信息，通常需要将不同模态的图像融合在一起，得到更丰富的信息， 以便了解病变组织或器官的综合信息，从而做出准确的诊断或制订出合适 的治疗方案。而图像配准是图像融合的重要前提，图像配准是指对一幅图 像进行一定的几何变换而映射到另一幅图像中，使得两幅图像中的相关点 达到空间上的一致。 图像配准主要有两大类方法，基于灰度的方法和基于特征的方法[1] 。基 于灰度的配准方法直接利用图像的灰度数据进行配准， 从而避免了因分割而 带来的误差，因而具有精度较高、鲁棒性强、不需要预处理而能实现自动配 准的特点。在基于灰度的配准方法中，基于互信息的方法包括互信息和归一 化互信息方法，它们已经被广泛使用并具有最高的精度。本文使用的是基于 互信息的配准方法。
图3
配准结果
从配准结果可以看出，对图像 2进行顺时针旋转 13.5506 °角 ， 压缩0.5379 倍后，两幅图像的互信息达到最大值1.3596 ，此时图像 配准已经完成，整个配准所用时间为46.176秒，所以本文所用的算 法的非常有效，可以满足像素级的配准精度要求。 3.2.所用到的 M 文件及其源代码 .m 3.2.1．ImageRegistration ImageRegistration.m
Abstract
This paper describes the basic principles and methods of the image registration, and proposes a method of image registration based on mutual information, and do an image registration experiment of two images using MATLAB, through experiments, we can see that the method is very effective, reaching the purposes of image registration..
Keywords : MATLAB; Entropy; Mutual Information; Image Registration
3.2.2 MI.m
function f= MI(x) %这是ImageRegistration.m子程序, 基于互信息的图像配准方法. load image IM1=double(IM 1); IM2=double(IM 2); IM2=imresize(IM 2, x(4), 'bilinear'); %x(4) 倍IM2大小的图像 J=imrotate(double(IM 2), x(3),'bilinear'); %逆时针旋转x(3) 度 J=abs(J)*255/max(max(J)); [n1 n2]=size(IM 1); [n3 n4]=size(J); if n1>n3-x(1)/2 f=1000; message=strvcat('The scaling factor is too small.', 'Press Ctrl+C to stop.','Increase x0(4) and restart.'); disp('Press Ctrl+C to stop.') Errordlg(message) pause; else if x(1)>n3-n1 x(1)=n3-n1-1; IM1(1:n1, 1:n2)=255; end if x(2)>n4-n2 x(2)=n4-n2-1; IM1(1:n1, 1:n2)=255; end if x(1)<0 x(1)=0; IM1(1:n1, 1:n2)=255; end if x(2)<0 x(2)=0; IM1(1:n1, 1:n2)=255; end xt=1:n1; yt=1:n2; xx=round(xt+x(1)); yy=round(yt+x(2)); IM2=round(J(xx, yy)); %选择图像2的部分去匹配图像1的大小。 rows=size(IM 1,1); cols=size(IM 2,2); N=256; h=zeros(N,N);
%互信息配准优化工具箱 % x0(1): 为裁剪旋转x0(3) 角的IM2行的第一个索引 % x0(2): 为裁剪旋转x0(3) 角的IM2列的第一个索引 % x0(3): 旋转的角度 % x0(4): 比例因子 clear all tic x0=[50; 50; -15; 0.5]; [x, fval]=fminsearch(@M I,x0) load image IM1=double(IM 1); IM2=double(IM 2); IM2=imresize(IM 2, x(4), 'bilinear'); J=imrotate(double(IM 2), x(3),'bilinear');
参考文献
[1] 高智勇 .钱勇先,林家瑞等 《多模态医学图像配准技术的分类与研究进展》 . 国外医学生物医学工程分 册,2000,23~206 [2] 赵永强,《认知系统中的信息论》.西北工业大学出版社,2009,127~148 [3] 陈显毅,《基于互信息的多模医学图像配准》,琼州学院学报, 第 l5 卷第 5 期，2008,37~40 [4] 刘正君,《MATLAB 科学计算与可视化仿真宝典》电子工业出版社,2009 ,405~409. [5] 董霖,《MATLAB 使用详解》科学出版社,2008 ,116~119.
图1
参考的MRI 图像
图2
待配准的CT图像
- 3 -
http://www.nwpu.edu.cn
然后利用 [x, fval]= fminsearch(@MI,x0) 进行配准 , fminsearch [5] 在 x0 附近寻 找函数MI的最小值，MI 函数值为图像1和变换后图像2的互信息的相反数。运 行该程序后，得出： x = [53.7948 49.6512 -13.5506 0.5379]T fval = -1.3596 Elapsed time is 46.176000 seconds. 并得出了配准后的图像，如图3所示。
3. 图像配准的实现
3.1. 实现过程：
配准所用的图像1为230 *230的核磁共振（MRI ）成像，即参考图像，如 图1所示，图像 2为512*512计算机 X射线断层扫描（ CT）成像，即待配准图 像，如图2所示。将两幅图像读入 MATLAB[4] ,利用rgb2gray转为灰度图像，并 把数据保存在ima ge.mat 中。
%启动秒表计时器 %初始点，它们选择尽可能接近的匹配点
%x(4) 倍im2大小的图像 %逆时针旋转x(3) 度,
- 4 -
http://www.nwpu.edu.cn [n1 n2]=size(IM 1); [n3 n4]=size(J); position=1:n1; xx=round(position+x(1)); yy=round(position+x(2)); IM2=round(J(xx, yy)); subplot(1,2,1), imshow(IM1, [ ]), title(' 参考图像1') subplot(1,2,2), imshow(IM2, [ ]), title(' 配准后的图像2') toc %阅读秒表计时器
Biblioteka Baidu
其中，H(A) 和 H(B) 分别是系统 A 和 B 的熵，
H ( A) = 1∑ pA( a) log pA (a)
a
H ( B ) = −∑ pB (b) log pB (b)
b
H(A，B) 是它们的联合熵,
H ( A , B) = − ∑ pAB (a, b) log pAB (a , b )
a ,b
2.2. 配准方法
基于互信息的配准方法 [3] 直接利用图像的灰度值实现两幅图像间的配 准，其算法过程如下： 首先根据两幅图像的基本情况预设一个初始参数 x0 ，其中 x0(1) 为裁剪 旋转 x0(3) 角的图像 2 行的第一个索引。x0(2) 为裁剪旋转 x0(3) 角的图像 2 列 的第一个索引，x0(3) 为旋转的角度，x0(4) 为比例因子。然后按照给定的初始 参数对图像 2 进行变换，并计算图像 1 和图像 2 的互信息，然后利用最优化 工具箱中的 fminsearch 函数在 x0 附近寻找使图像 1 和图像 2 互信息最大的 点，直至搜索到满足精度要求的参数；最后输出配准参数。
2. 基于互信息的配准方法
2.1.互信息的计算
互信息[2] 是信息论中的一个基本概念，通常用于描述两个系统间的统计 相关性，或者是一个系统中所包含的另一个系统中信息的多少，它可以用信 息熵来描述 ： I(A,B)= H(A)+H(B)-H(A,B)
- 2 -
http://www.nwpu.edu.cn