图像配准的方法

转自：http://blog.sina.com.cn/s/blog_4b9b714a0100d5k5.html

 图像配准的方法

1  基于特征的图像配准

基于特征的图像配准首先提取图像信息的特征，然后以这些特征为模型进行配准。特征提取的结果是一含有特征的表和对图像的描述，每个特征由一组属性表示，对属性的进一步描述包括边缘的定向和弧度、区域的大小等。局部特征之间存在着相互关系，如几何关系、辐射度量关系、拓扑关系等。可以用这些局部特征之间的关系描述全局特征。通常基于局部特征配准大多都是基于点、线或边缘的，而全局特征的配准则是利用局部特征之间的关系进行配准的方法。

由于图像的特征点比图像的像素点要少很多，因此大大减少了配准过程的计算量，但特征提取方法的计算代价通常较大，不便于实时应用。特征点的配准度量值对位置的变化比较敏感，可以大大提高配准的精确程度。对于纹理较少的图像区域提取的特征的密度通常比较稀少，局部特征的提取就比较困难。特征点的提取过程可以减少噪声的影响，对灰度变化、图像形变和遮挡等都有较好的适应能力。因此，在图像配准领域得到了广泛应用。基于特征的图像配准方法有两个重要环节：特征提取和特征配准。

2  基于互信息的图像配准

医学图像配准技术从基于特征的配准方法发展到基于统计的配准方法有其突破性的意义。与基于特征的配准方法相比，基于统计的配准方法的突出优点为鲁棒性好、配准精度高、人工干预少。基于统计的配准方法通常是指最大互信息的图像配准方法。

基于互信息的图像配准是用两幅图像的联合概率分布与完全独立时的概率分布的广义距离来估计互信息，并作为多模态医学图像配准的测度。当两幅基于共同的解剖结构的图像达到最佳配准时，它们的对应像素的灰度互信息应为最大。由于基于互信息的配准对噪声比较敏感，首先，通过滤波和分割等方法对图像进行预处理。然后进行采样、变换、插值、优化从而达到配准的目的。

基于互信息的配准技术属于基于像素相似性的方法。它基于图像中所有的像素进行配准，基于互信息的图像配准引入了信息论中的概念，如熵、边缘熵、联合熵和互信息等，可使配准精度达到亚像素级的高精度。    

基于互信息只依赖于图像本身的信息，不需要对图像进行特征点提取和组织分类等预处理，是一种自动而有效的配准算法。该算法可靠，对图像中的几何失真、灰度不均匀和数据的缺失等不敏感。不依赖于任何成像设备，可应用于多模态医学图像配准。基于互信息的图像配准也有其缺点，它运算量大，对噪声敏感，要求待配准图像间联合概率分布函数必须是严格正性的。

【例1】计算两幅图像的互信息。

function mi = MI(a,b)

%Caculate MI of a and b in the region of the overlap part

%计算重叠部分

[Ma,Na] = size(a);

[Mb,Nb] = size(b);

M=min(Ma,Mb);

N=min(Na,Nb);

%初始化直方图数组

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;

%统计直方图

for i=1:M

for j=1:N

indexx = a(i,j);

indexy = b(i,j) ;

hab(indexx,indexy) = hab(indexx,indexy)+1;%联合直方图

ha(indexx) = ha(indexx)+1;%a图直方图

hb(indexy) = hb(indexy)+1;%b图直方图

end

end

%计算联合信息熵

hsum = sum(sum(hab));

index = find(hab~=0);

p = hab/hsum;

Hab = sum(sum(-p(index).*log(p(index))));

%计算a图信息熵

hsum = sum(sum(ha));

index = find(ha~=0);

p = ha/hsum;

Ha = sum(sum(-p(index).*log(p(index))));

%计算b图信息熵

hsum = sum(sum(hb));

index = find(hb~=0);

p = hb/hsum;

Hb = sum(sum(-p(index).*log(p(index))));

%计算a和b的互信息

mi = Ha+Hb-Hab;

%计算a和b的归一化互信息

%mi = hab/(Ha+Hb);

3  医学图像配准

目前，图像配准技术在医学领域应用得相当广泛。医学图像配准具有很重要的临床应用价值。对各种使用相同或者不同的成像手段所获得的医学图像进行配准不仅可以应用于医疗诊断，也可以应用于手术计划的制定、治疗计划的制定、病理情况的跟踪和治疗效果的评价等各个方面。

下面给出一个具体的图像配准的例子来形象说明配准的应用。

(1) 读取图像数据。因为源图像以矩阵形式存在一个二进制的文件里，用fread可将其读取到变量矩阵中。将读取文件编制成一个子函数(RTIread.m)。源代码如下：

function imMatrix=RTIread(FILENAME,SIZE)

%RTIread Read the image matrix from binary "Registration Test Image" file.

% imMatrix=RTIread(FILENAME,SIZE) opens the file FILENAME, and reads the

% number of elements specified by SIZE.

%

% FILENAME is a string containing the name of the file to be opened.

% Valid entries for SIZE are:

% N read N elements into a column vector.

% inf read to the end of the file.

% [M,N] read elements to fill an M-by-N matrix, in column order.

% N can be inf, but M can't.

%

% It returns the image matrix.

fid=fopen(FILENAME,'r');

imMatrix=fread(fid,SIZE,'uint8=>uint8');

fclose(fid);

这里选取了两张256×256的图片，文件名为sag1.tif和sag2.tif。

运行以下代码读取图像矩阵：

base=RTIread('sag1.tif',[256,256]);

input=RTIread(sag2.tif,[256,256]);

(2) 选取匹配点。根据预定的配准方法，选定足够的匹配点对。运行下列代码：

cpselect(input,base); %please select 15 points for test.

注意：

这里允许选择15个点进行测试。

操作很简单，只需注意选点要均匀布开，以增加其代表性。选定完毕，再选择File→Save Points to Workspace命令将数据保存到工作区中。Workspace立刻多出两个N×2的数组(其中N为选定的匹配点对数)，分别为input_points和base_points。

(3) 利用十字相关法调整选定了的匹配点。如果对上面一个步骤的选点工作不满意，可以进行这一步。如果感觉尚可，则跳过这一步，以减小运算量。

运行代码：

input_points_corr = cpcorr(input_points,base_points,input,base);

%优化选择点的作用

input_points_corr为优化后在输入图片的对应匹配点。

(4) 计算变换公式的参数。利用cp2tform选定变换类型(配准方法)，计算变换参数。以下只需选定一种即可。

% (1) not Fine-tune points

Tlinear = cp2tform(input_points,base_points,'linear conformal');

Taffine = cp2tform(input_points,base_points,'affine');

Tprojective = cp2tform(input_points,base_points,'projective');

Tpolynomial2 = cp2tform(input_points,base_points,'polynomial',2);

Tpolynomial3 = cp2tform(input_points,base_points,'polynomial',3);

Tpolynomial4 = cp2tform(input_points,base_points,'polynomial',4);

Tpiecewise = cp2tform(input_points,base_points,'piecewise linear');

Tlwm = cp2tform(input_points,base_points,'lwm');

% (2)Fine-tune points

fTlinear = cp2tform(input_points_corr,base_points,'linear conformal');

fTaffine = cp2tform(input_points_corr,base_points,'affine');

fTprojective = cp2tform(input_points_corr,base_points,'projective');

fTpolynomial2 = cp2tform(input_points_corr,base_points,'polynomial',2);

fTpolynomial3 = cp2tform(input_points_corr,base_points,'polynomial',3);

fTpolynomial4 = cp2tform(input_points_corr,base_points,'polynomial',4);

fTpiecewise = cp2tform(input_points_corr,base_points,'piecewise linear');

fTlwm = cp2tform(input_points_corr,base_points,'lwm');

(5) 变换图像。

% Transform the unregistered image to bring it into alignment.

title('image registration polynomial method');

subplot(2,2,1);

imshow(base);

title('Base image');

subplot(2,2,2);

imshow(input);

title('Input image');

subplot(2,2,3);

imshow(imtransform(input,Tpolynomial2));

title('registered image');

subplot(2,2,4);

imshow(imtransform(input,fTpolynomial2));

title('registered image(fine-tune points)');

小结：

(1) 选择适当的方法来建立转换参数，并非算法越复杂越好，应参考成像因素(退化因素)。

(2) 利用cpselect选择匹配点，cpselect可以返回一个GUI句柄。在使用中，应当尽量选择恰当的匹配区域，这样能够得到较好的效果。

上面的这个例子典型地说明了图像配准的步骤：提取图像特征，根据特征对进行几何变换，执行整个变换。