【图像配准】基于matlab GUI Powell+蚁群算法图像配准【含Matlab源码 928期】

一、简介

1 引言
医生在对病人疾病进行诊断和治疗时, 经常需要采用X线, CT, MRI, PET, SPECT等多种成像模式对病人进行成像, 以提供病人解剖形态方面和功能代谢方面的互补信息。但不同模式下的图像会出现分辨率不同、灰度属性不同、病灶的位置不同等情况, 为了能将不同模式的图像中的信息融合在一起, 就必须首先进行图像的配准。多模医学图像配准就是将一幅图像中的像素空间位置进行变换, 以使其与另一幅图像的像素在空间位置上对齐, 从而将多种成像模式的图像数据融合起来, 利用各自图像信息的特点, 在一幅图像上表达来自人体多方面的信息。

基于最大互信息的配准方法由于直接使用图像像素灰度信息的统计特性即互信息作为配准的依据, 不需要提取图像的解剖特征, 因此它是一种精度高、稳健性强的方法, 在医学图像配准领域得到了普遍关注和广泛应用。在基于互信息的医学图像配准中目前使用得最多的优化算法是Powell法。Powell法不需要计算导数, 在每一维中使用Brent算法迭代搜索, 搜索速度比较快, 局部寻优能力极强, 在局部搜索中精度要高于其它的优化算法;但由于互信息函数存在很多局部极值, 并且Powell法存在初始点依赖问题, 即在配准过程中很容易落入局部最优中, 最终的优化结果很大程度上依赖于初始点所在的位置, 这就使得Powell法的优化过程常常收敛到局部极值而得到错误的配准参数, 给互信息函数优化带来很多不便。为了克服Powell法的缺点, Plattard等人在配准优化过程中将Powell法和单纯形法结合起来, Jenkinson和Smith把多起始点策略引入到Powell法优化中, Wachowiak等人则提出了基于Powell法、禁忌搜索、粒子群优化算法等优化方法的多种混合优化策略并应用到图像配准中。我们采用基于小波变换的多分辨率策略, 将蚁群算法和Powell法结合起来对目标函数进行寻优, 提高了配准的准确度和速度。该方法很好地解决了Powell法的初始点依赖问题, 并且其Powell优化过程对参数优化顺序不敏感, 不像现有的很多方法需要根据成像特点来设定Powell法的参数优化顺序, 因此具有很好的实用性。实验结果表明, 本文方法能够有效地克服局部极值问题, 大大地提高了配准精度。

2 基于互信息的配准方法
2.1 基于互信息的图像配准原理
互信息是信息论中的一个基本概念, 用来描述两个随机变量间的统计相关性, 是一个变量包含另一个变量的信息量的多少的度量。它可用熵来描述:

其中H (A) 和H (B) 分别为图像A和B的熵, H (A, B) 为二者的联合熵。在多模图像配准中, 当两幅图像的空间位置完全一致时, 其中一幅图像中表达的关于另一幅图像的信息, 也就是互信息I (A, B) 为最大。

由于互信息对重叠区域的变化比较敏感, Studholme和Maes分别提出了两种归一化互信息的表现形式:

归一化互信息能更好地反映配准函数的变化。

2.2 配准变换模型
本文中变换模型采用的是刚体变换。对待配准的两幅图像, 首先建立统一的立体坐标系统, 坐标原点定义在图像的灰度重心, 将两幅图像进行粗配准。选择一幅图像作为参考图像R, 另一幅图像作为浮动图像F, 从浮动图像的空间坐标PF到参考图像的空间坐标PR的刚体变换可以用下式描述:

其中VR和VF为3×3的对角阵, 分别表示图像F和R的像素大小;CF和CR分别是两幅图像的中心;Rx (φx) , Ry (φy) , Rz (φz) 都是3×3的旋转矩阵, φx, φy, φz分别是绕x轴, y轴和z轴的旋转角度;t是平移向量, tx, ty, tz分别是在x轴, y轴和z轴上的平移距离。

基于互信息的配准过程是一个多参数的优化过程, 即搜索使两幅图像间的互信息最大的6个空间变换参数tx, ty, tz, φx, φy, φz。

2.3 插值与多分辨率配准
浮动图像上的点通过空间变换后得到的点的坐标不一定是整数, 需要通过插值方法来得到变换点的灰度值。常用的插值方法有最近邻法、线性插值、三次样条等, 在基于互信息的三维图像配准中经常采用三线性PV插值算法 (trilinear Partial Volume distribution interpolation) 。PV插值算法不会引入新的灰度值, 浮动图像中的一点p的灰度f § 对联合直方图的贡献是由参考图像中的点q的周围最近邻的8个点取与三线性插值算法相同的权重加权而得。

基于互信息的配准方法需要计算联合直方图和进行图像插值, 配准时间比较长, 尤其是数据量很大的三维图像配准。为了提高配准速度, 我们利用多分辨率的配准方法。小波变换能够产生图像数据的多分辨率表示。即使是在低分辨率时, 小波分解后的近似分量子图像也保留了原始数据的大多数重要特征。使用这种多分辨率数据, 可以一开始在低分辨率图像进行搜索, 然后上升到高分辨率图像上来进一步细化配准结果, 这样可以大大地减少搜索的数据量。对基于互信息的三维多模医学图像配准来说, 二级多分辨率策略即可满足配准精度要求。

3 优化算法
3.1 蚁群算法
蚁群算法 (ant colony algorithm) 是意大利学者Dorigo受自然界中真实蚁群集体行为的启发而提出来的一种基于种群的模拟进化算法, 属于带构造性特征的随机搜索算法。蚁群算法易于与其它的优化方法结合, 具有并行性和很强的鲁棒性。
假设优化的问题为

其中f (x) :Rs→R为已知的多维函数, [x0, xf]为已知的s维解空间。设m只人工蚂蚁, 每只蚂蚁刚开始随机地位于解空间的 (n1×n2×…×ns) 个等分区域的某处, 并按下式定义的概率实现状态转移:

其中pij表示蚂蚁从解空间区域i转移到区域j的概率, τj为区域j的吸引强度, 期望值ηij定义为ηij=fjmaxfimax, 即蚁群在区域j与区域i已经搜索到的空间位置上目标函数最大值的差值, 给定参数α, β>0为启发式因子, 分别表示蚂蚁在状态转移过程中区域吸引强度τj和期望值ηij对蚂蚁转移概率所起的不同作用。
区域j吸引强度的更新方程为

式中∆τjk反映本次循环中第k只蚂蚁在区域j的局部搜索中吸引强度的增加, Ljk表示本次循环中第k只蚂蚁在区域j的局部搜索中目标函数值的变化量;参数ρ∈ (0, 1) , 体现解空间中各个等分区域中吸引强度的持久性;参数Q是一常数, 为蚂蚁释放的信息素密度;算法中有关的初始值取为τj (0) =C, ∆τj (0) =0。区域i中的蚂蚁的转移及搜索规则定义为

其中neighbor i表示区域i的相邻区域。每只蚂蚁要么以上述规则从当前区域转移到其他相邻区域中作局部随机搜索, 要么在当前区域内进行局部随机搜索。于是多维函数f (x) 的寻优就借助于m只蚂蚁在解空间的等分区域间的不断移动以及某些区域内的局部随机搜索来进行, 这种寻优方式相当于一群蚂蚁对定义域进行有穷尽的且在先验知识引导下的随机搜索, 最终收敛到问题的近似全局最优解。

函数优化问题的蚁群算法步骤如下:
(1) count←0 (count为迭代次数) , 设置算法参数及解空间的分区数, 将m只蚂蚁随机放置于初始区域上;
(2) 每只蚂蚁以式 (6) , 式 (9) 给出的概率规则转移或作局部搜索;
(3) 存储各区域当前已搜索到的目标函数最大值;
(4) 记录当前最好解xmax及最优值fmax;
(5) 按式 (7) , 式 (8) 更新各区域的吸引强度∆τj, count←count+1;
(6) 若count小于预定的迭代次数, 则转到 (2) ;
(7) 输出最佳结果。

蚁群算法中的主要参数一般设置如下:1≤α, β≤5;0.5≤ρ≤1, 取0.7左右为最佳;1≤Q≤10000, Q的取值对算法影响不大;蚂蚁数目m和解空间的分区数有关, 分区数越大则蚂蚁数目越多。

3.2 蚁群算法与Powell法结合的配准优化算法
蚁群算法在一般情况下均能找到比较满意的结果, 但得到的解不一定是解空间的最优解, 通常是全局最优解附近的一个解。而Powell算法有极强的局部寻优能力, 所以在本文的配准优化过程中我们将蚁群算法的全局搜索能力和Powell算法的局部寻优能力有机结合起来。由于蚁群算法中目标函数的计算次数比较多, 从而优化时间比较长, 为此我们采用了基于小波变换的多分辨率策略, 优化过程分为两步, 第1步首先在较低分辨率的图像上采用蚁群算法进行配准, 此时图像比较小, 互信息计算速度快, 优化过程能较快完成;第2步采用Powell法在高分辨率图像上进行寻优, 算法的初始点为上一步中蚁群算法得到的最好解。该算法的具体步骤描述如下:

(1) 首先对待配准的浮动图像和参考图像进行小波分解, 得到比较小的子图像;对子图像进行配准, 采用蚁群算法进行寻优。优化过程中插值方法采用最近邻法, 联合灰度直方图的灰度级数设置为64。在图像分辨率比较低时, 采用最近邻法插值速度快, 效果和PV插值差不多, 将图像灰度归一化到较少的灰度级别能够提高配准精度和速度。
(2) 将步骤 (1) 中蚁群算法得到的配准参数作为Powell法寻优的起始点, 在高分辨率的图像上进行配准, 得到配准的最优解。为提高配准精度, 优化过程中插值方法采用PV插值, 图像灰度级数设置为256。
由于步骤 (1) 中蚁群算法的配准结果为步骤 (2) 中的Powell法优化提供了一个非常有效的初始点, 使得Powell法的参数优化顺序对其优化结果的影响很小, 因此不需要根据成像特点来设定参数优化顺序, 同时也使得步骤 (2) 的优化时间要比从默认初始点或者随机初始点开始搜索的普通的Powell法花费的时间短很多。步骤 (1) 中的蚁群算法是在分辨率低的图像上进行寻优, 并且插值方法简单, 所以互信息的计算比较快, 步骤 (1) 的优化过程能在比较短的时间内完成, 整个算法总的运行时间与普通的Powell法相当。

二、部分源代码

function varargout = ImageRegistration(varargin)
% IMAGEREGISTRATION M-file for ImageRegistration.fig
%      IMAGEREGISTRATION, by itself, creates a new IMAGEREGISTRATION or raises the existing
%      singleton*.
%
%      H = IMAGEREGISTRATION returns the handle to a new IMAGEREGISTRATION or the handle to
%      the existing singleton*.
%
%      IMAGEREGISTRATION('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in IMAGEREGISTRATION.M with the given input arguments.
%
%      IMAGEREGISTRATION('Property','Value',...) creates a new IMAGEREGISTRATION or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before ImageRegistration_OpeningFcn gets called.  An
%      unrecognized property name or invalid value makes property
%      application
%      stop.  All inputs are passed to ImageRegistration_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help ImageRegistration

% Last Modified by GUIDE v2.5 24-May-2021 17:01:50

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @ImageRegistration_OpeningFcn, ...
                   'gui_OutputFcn',  @ImageRegistration_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin && ischar(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT

addpath(pwd);
% --- Executes just before ImageRegistration is made visible.
function ImageRegistration_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to ImageRegistration (see VARARGIN)

% Choose default command line output for ImageRegistration
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes ImageRegistration wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = ImageRegistration_OutputFcn(hObject, eventdata, handles) 
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

% Get default command line output from handles structure
varargout{1} = handles.output;


% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
clc
%%%%%%%%%%%%%%%调用OpenImage.m读入参考图像并获取文件名、图像大小%%%%%%%%%%%%
[fname,pname,index]=uigetfile({'*.jpg;*.bmp;*.png;*.tif;*.gif;*.jpeg'},'选择参考图像');
if index
    name=[pname fname];
    temp=imread(name);
    [m,n]=size(temp);
    m=num2str(m);
    n=num2str(n);
    imsize=['  ',m,'*',n];
    handles.ImsizeI=imsize;
    handles.filenameI=name;
    handles.names_dispI=[fname,imsize];
end
set(handles.text2,'String',handles.names_dispI);
guidata(hObject,handles);

%%%%%%%%%%%%%%%%%%显示参考图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes(handles.axes1)
I=imread(handles.filenameI);
save I;
imshow(I)
figure(5)
imshow(I)


% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton2 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
clc
%%%%%%%%%%%调用OpenImage.m读入浮动图像并获取并获取文件名、图像大小%%%%%%%%%%
[fname,pname,index]=uigetfile({'*.jpg;*.bmp;*.png;*.tif;*.gif;*.jpeg'},'选择浮动图像');
if index
    name=[pname fname];
    temp=imread(name);
    [m,n]=size(temp);
    m=num2str(m);
    n=num2str(n);
    imsize=['  ',m,'*',n];
    handles.ImsizeJ=imsize;
    handles.filenameJ=name;
    handles.names_dispJ=[fname,imsize];
end
set(handles.text3,'String',handles.names_dispJ);
guidata(hObject,handles);

%%%%%%%%%%%%%%%%%%显示浮动图像%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes(handles.axes2)
J=imread(handles.filenameJ);
imshow(J)
figure(6)
imshow(J)


% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject    handle to pushbutton3 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
clc
%%%%%%%%%%%%%%%%检测是否已输入参考图像与浮动图像%%%%%%%%%%%%%%%%%%%%%%%%%%
axesIbox=get(handles.axes1,'box');
axesJbox=get(handles.axes2,'box');
if strcmp(axesIbox,'off')||strcmp(axesJbox,'off')
    errordlg('Please select Image for Registration','Error')
    error('NO Image!')
end

%%%%%%%%%%%%%%%检测参考图像与浮动图像大小是否相同%%%%%%%%%%%%%%%%%%%%%%%%%
handles.isSameSizeIJ=strcmp(handles.ImsizeI,handles.ImsizeJ);
if handles.isSameSizeIJ~=1
    errordlg('Please select the Same Size Image','Error')
    error('Image Size doesn"t match!')
end

%%%%%%%%%%%%%读入并复制图像，一幅用于配准过程，另一幅用于配准后输出%%%%%%%%%
I1=imread(handles.filenameI);
J1=imread(handles.filenameJ);
handles.Old_I=I1;
handles.Old_J=J1;
I=I1;
J=J1;
handles.I=I;
handles.J=J;

guidata(hObject,handles);

%%%%%%%%%%%%%显示配准前参考图像与浮动图像的“融合"效果图%%%%%%%%%%%%%%%%%%%%
axes(handles.axes4)
imshow(uint8(I+J))

%%%%%%%%%%%%%%%调用函数GLPF.m(高斯低通滤波器)，完成高斯低通预处理%%%%%%%%%%%%%
[I,J]=GLPF(handles);
handles.I=I;
handles.J=J;
guidata(hObject,handles);

%%%%%%%%%%%%%%调用函数Powell.m，实现图像配准%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
RegistrationParameters=Powell(handles);
%RegistrationParameters=AntColony(handles);
%RegistrationParameters=AntColonyimprove(handles);
toc
ElapsedTime=toc;

%%%%%%%%%%%%%%显示配准结果%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
handles.RegistrationParameters=RegistrationParameters;
y=RegistrationParameters(1);
x=RegistrationParameters(2);
ang=RegistrationParameters(3);
RegistrationResult=sprintf('X,Y,Angle=[%.5f][%.5f][%.5f]',x,y,ang);
MI_Value=RegistrationParameters(4);
MI_Value=sprintf('MI_Value=[%.4f]',MI_Value);
ElapsedTime=sprintf('Elapsed Time=[%.3f]',ElapsedTime);
axes(handles.axes5)
[FusionImage,RegistrationImage]=Fusion(handles);
imshow(FusionImage)
axes(handles.axes3)
imshow(RegistrationImage)
figure(7)
imshow(RegistrationImage)
set(handles.text8,'string',RegistrationResult);
set(handles.text9,'string',MI_Value);
set(handles.text10,'string',ElapsedTime);

  

 

  
1
  
2
  
3
  
4
  
5
  
6
  
7
  
8
  
9
  
10
  
11
  
12
  
13
  
14
  
15
  
16
  
17
  
18
  
19
  
20
  
21
  
22
  
23
  
24
  
25
  
26
  
27
  
28
  
29
  
30
  
31
  
32
  
33
  
34
  
35
  
36
  
37
  
38
  
39
  
40
  
41
  
42
  
43
  
44
  
45
  
46
  
47
  
48
  
49
  
50
  
51
  
52
  
53
  
54
  
55
  
56
  
57
  
58
  
59
  
60
  
61
  
62
  
63
  
64
  
65
  
66
  
67
  
68
  
69
  
70
  
71
  
72
  
73
  
74
  
75
  
76
  
77
  
78
  
79
  
80
  
81
  
82
  
83
  
84
  
85
  
86
  
87
  
88
  
89
  
90
  
91
  
92
  
93
  
94
  
95
  
96
  
97
  
98
  
99
  
100
  
101
  
102
  
103
  
104
  
105
  
106
  
107
  
108
  
109
  
110
  
111
  
112
  
113
  
114
  
115
  
116
  
117
  
118
  
119
  
120
  
121
  
122
  
123
  
124
  
125
  
126
  
127
  
128
  
129
  
130
  
131
  
132
  
133
  
134
  
135
  
136
  
137
  
138
  
139
  
140
  
141
  
142
  
143
  
144
  
145
  
146
  
147
  
148
  
149
  
150
  
151
  
152
  
153
  
154
  
155
  
156
  
157
  
158
  
159
  
160
  
161
  
162
  
163
  
164
  
165
  
166
  
167
  
168
  
169
  
170
  
171
  
172
  
173
  
174
  
175
  
176
  
177
  
178
  
179
  
180
  
181
  
182
  
183
  
184
  
185
  
186
  
187
  
188
  
189
  
190
  
191
  
192
  
193
  
194
  
195
  
196
  
197
  
198
  
199
  
200
  
201
  
202
  
203
  
204
  
205
  
206
  
207
  
208
 

三、运行结果

四、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]杨帆,张汗灵.蚁群算法和Powell法结合的多分辨率三维图像配准[J].电子与信息学报. 2007,(03)

文章来源: qq912100926.blog.csdn.net，作者：海神之光，版权归原作者所有，如需转载，请联系作者。

原文链接：qq912100926.blog.csdn.net/article/details/117226429