✅博主简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，Matlab项目合作可私信。
🍎个人主页：海神之光
🏆代码获取方式：
海神之光Matlab王者学习之路—代码获取方式
⛳️座右铭：行百里者，半于九十。
更多Matlab仿真内容点击👇
Matlab图像处理（进阶版）
路径规划（Matlab）
神经网络预测与分类（Matlab）
优化求解（Matlab）
语音处理（Matlab）
信号处理（Matlab）
车间调度（Matlab）
⛄一、简介
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);
⛄三、运行结果
⛄四、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1]杨帆,张汗灵.蚁群算法和Powell法结合的多分辨率三维图像配准[J].电子与信息学报. 2007,(03)
3 备注
简介此部分摘自互联网，仅供参考，若侵权，联系删除
🍅 仿真咨询
1 各类智能优化算法改进及应用
生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化
2 机器学习和深度学习方面
卷积神经网络（CNN）、LSTM、支持向量机（SVM）、最小二乘支持向量机（LSSVM）、极限学习机（ELM）、核极限学习机（KELM）、BP、RBF、宽度学习、DBN、RF、RBF、DELM、XGBOOST、TCN实现风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断
3 图像处理方面
图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知
4 路径规划方面
旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、车辆协同无人机路径规划、天线线性阵列分布优化、车间布局优化
5 无人机应用方面
无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配
6 无线传感器定位及布局方面
传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化
7 信号处理方面
信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化
8 电力系统方面
微电网优化、无功优化、配电网重构、储能配置
9 元胞自动机方面
交通流 人群疏散 病毒扩散 晶体生长
10 雷达方面
卡尔曼滤波跟踪、航迹关联、航迹融合
4758
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
