【图像配准】多模态非刚性图像配准算法附matlab代码
✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。
🍎个人主页：Matlab科研工作室
🍊个人信条：格物致知。
更多Matlab完整代码及仿真定制内容点击👇
智能优化算法 神经网络预测 雷达通信 无线传感器 电力系统
信号处理 图像处理 路径规划 元胞自动机 无人机
🔥 内容介绍
图像配准是计算机视觉领域中一个重要的问题，它涉及将多个图像对齐在同一个坐标系统中。在许多应用中，我们需要将来自不同传感器或不同时间点的图像进行配准，以便进行比较、分析或融合。而多模态非刚性图像配准则是一种特殊的图像配准问题，它需要解决不同模态（如磁共振成像和计算机断层扫描）之间的配准问题，并且还要考虑到图像中的非刚性变形。
在本文中，我们将介绍一种常用的多模态非刚性图像配准算法，并详细阐述其步骤和原理。
第一步是特征提取。在图像配准中，我们需要找到一些具有代表性的特征点来进行匹配。对于多模态图像配准，我们可以使用一些特定的特征提取方法，如局部特征描述子（如SIFT、SURF）或基于深度学习的方法（如CNN）。这些方法可以从图像中提取出鲁棒性较强的特征点，以便进行后续的匹配。
第二步是特征匹配。在这一步中，我们需要将来自不同图像的特征点进行匹配。这可以通过计算特征点之间的距离或相似性来实现。一种常用的方法是使用最近邻算法（如K最近邻算法）来找到两个图像中最相似的特征点对。匹配的结果将会是一组对应的特征点对。
第三步是变换模型的选择。在多模态非刚性图像配准中，我们需要考虑到图像中的非刚性变形，因此我们需要选择一个合适的变换模型来描述图像之间的变换关系。常用的变换模型包括仿射变换、弹性变形模型和非线性变形模型等。选择合适的变换模型可以更好地适应图像之间的变形情况。
第四步是优化。在这一步中，我们需要通过优化算法来估计变换模型的参数，使得匹配的特征点对之间的误差最小化。常用的优化算法包括最小二乘法、梯度下降法和Levenberg-Marquardt算法等。通过迭代优化算法，我们可以得到最佳的变换模型参数。
第五步是图像插值和重采样。在配准完成后，我们需要将图像进行插值和重采样，以便得到最终的配准结果。插值算法可以根据配准结果的像素坐标在原始图像中找到对应的像素值，并进行插值计算。重采样算法可以将配准后的图像调整为特定的大小和分辨率。
最后，我们需要评估配准的质量。在多模态非刚性图像配准中，我们可以使用一些评估指标来评估配准结果的质量，如均方差（MSE）、互信息（MI）和结构相似性指数（SSIM）等。这些指标可以帮助我们了解配准结果与原始图像之间的差异和相似性。
综上所述，多模态非刚性图像配准是一个复杂而重要的问题，在许多医学图像处理和计算机视觉应用中起着关键作用。通过特征提取、特征匹配、变换模型选择、优化、图像插值和重采样以及质量评估等步骤，我们可以实现准确而鲁棒的多模态非刚性图像配准。这将为我们提供更多的图像信息，帮助我们进行更深入的分析和研究。
📣 部分代码
function [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,Options)% This function register_images is the most easy way to register two% images both affine and nonrigidly.%% Features:% - It can be used with images from different type of scans or modalities.% - It uses both a rigid transform and a nonrigid registration.% - It uses multilevel refinement% - It can be used with images of different sizes.% - The function will automaticaly detect single modality or multiple% modalities, and choose the right registration method.%% [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,Options);%% Inputs,% Imoving : The image which will be registerd% Istatic : The image on which Imoving will be registered% Options : Registration options, see help below%% Outputs,% Ireg : The registered moving image% Bx, By : The backwards transformation fields of the pixels in % x and y direction seen from the static image to the moving image.% Fx, Fy : The (approximated) forward transformation fields of the pixels in % x and y direction seen from the moving image to the static image.% (See the function backwards2forwards)% Options,% Options.SigmaFluid : The sigma of the gaussian smoothing kernel of the pixel% velocity field / update field, this is a form of fluid% regularization, (default 4)% Options.SigmaDiff : The sigma for smoothing the transformation field% is not part of the orignal demon registration, this is % a form of diffusion regularization, (default 1)% Options.Interpolation : Linear (default) or Cubic.% Options.Alpha : Constant which reduces the influence of edges (and noise)% and limits the update speed (default 4). % Options.Similarity : Choose 'p' for single modality and 'm' for% images of different modalities. (default autodetect)% Options.Registration: Rigid, Affine, NonRigid % Options.MaxRef : Maximum number of grid refinements steps.% Options.Verbose: Display Debug information 0,1 or 2 % % Notes,% In case of Multiple Modalities affine registration is done with mutual% information. The non-rigid registration is done by first doing a% modality transformation (paints regions in image 1 with the intensity% pallette of those regions in image2 and visa versa), and than % using "normal" pixel based demon registration. See MutualTransform.m% %% Example,% % Read two greyscale images of Lena% Imoving=imread('images/lenag1.png'); % Istatic=imread('images/lenag3.png');% % % Register the images% [Ireg,Bx,By,Fx,Fy] = register_images(Imoving,Istatic,struct('Similarity','p'));%% % Show the registration result% figure,% subplot(2,2,1), imshow(Imoving); title('moving image');% subplot(2,2,2), imshow(Istatic); title('static image');% subplot(2,2,3), imshow(Ireg); title('registerd moving image');% % Show also the static image transformed to the moving image% Ireg2=movepixels(Istatic,Fx,Fy);% subplot(2,2,4), imshow(Ireg2); title('registerd static image');%% % Show the transformation fields% figure,% subplot(2,2,1), imshow(Bx,[]); title('Backward Transf. in x direction');% subplot(2,2,2), imshow(Fx,[]); title('Forward Transf. in x direction');% subplot(2,2,3), imshow(By,[]); title('Backward Transf. in y direction');% subplot(2,2,4), imshow(Fy,[]); title('Forward Transf. in y direction');%% % Calculate strain tensors% E = strain(Fx,Fy);% % Show the strain tensors% figure,% subplot(2,2,1), imshow(E(:,:,1,1),[]); title('Strain Tensors Exx');% subplot(2,2,2), imshow(E(:,:,1,2),[]); title('Strain Tensors Exy');% subplot(2,2,3), imshow(E(:,:,2,1),[]); title('Strain Tensors Eyx');% subplot(2,2,4), imshow(E(:,:,2,2),[]); title('Strain Tensors Eyy');%% Example Multi-Modalities% % Read two brain images % Imoving=im2double(imread('images/brain_T1_wave.png')); % Istatic=im2double(imread('images/brain_T2.png'));%% % Register the images% [Ireg,Bx,By] = register_images(Imoving,Istatic,struct('SigmaFluid',4));%% figure,% subplot(1,3,1), imshow(Imoving); title('moving image');% subplot(1,3,2), imshow(Istatic); title('static image');% subplot(1,3,3), imshow(Ireg); title('registerd moving image');%% % Read normal T1 image and transformation field% Inormal=im2double(imread('images/brain_T1.png'));% load('images/wave_field.mat');%% % Show the difference with ideal image% figure, imshow(Imoving-Inormal,[-0.5 0.5]); title('unregistered')% figure, imshow(Ireg-Inormal,[-0.5 0.5]); title('registered');% disp(['pixel abs difference : ' num2str(sum(abs(Imoving(:)-Inormal(:))))])% disp(['pixel abs difference : ' num2str(sum(abs(Imoving(:)-Ireg(:))))])%% % Show Warp field% figure,% subplot(2,2,1), imshow(BxNormal,[-20 20]); title('Bx Normal');% subplot(2,2,2), imshow(Bx,[-20 20]); title('Bx');% subplot(2,2,3), imshow(ByNormal,[-20 20]); title('By Normal');% subplot(2,2,4), imshow(By,[-20 20]); title('By');% Function is written by D.Kroon University of Twente (March 2009)% add all needed function pathstry functionname='register_images.m'; functiondir=which(functionname); functiondir=functiondir(1:end-length(functionname)); addpath([functiondir '/functions']) addpath([functiondir '/functions_affine']) addpath([functiondir '/functions_nonrigid'])catch me disp(me.message);end% Disable warningwarning('off', 'MATLAB:maxNumCompThreads:Deprecated')% Process inputsdefaultoptions=struct('Similarity',[],'Registration','NonRigid','MaxRef',[],'Verbose',2,'SigmaFluid',4,'Alpha',4,'SigmaDiff',1,'Interpolation','Linear');if(~exist('Options','var')), Options=defaultoptions; else tags = fieldnames(defaultoptions); for i=1:length(tags) if(~isfield(Options,tags{i})), Options.(tags{i})=defaultoptions.(tags{i}); end end if(length(tags)~=length(fieldnames(Options))), warning('register_images:unknownoption','unknown options found'); endend% Set parametersMaxRef=Options.MaxRef;% Start time measurementif(Options.Verbose>0), tic; end% Store the class of the inputsIclass=class(Imoving);% Convert the inputs to doubleImoving=im2double(Imoving);Istatic=im2double(Istatic);% Resize the moving image to fit the static imageif(sum(size(Istatic)-size(Imoving))~=0) Imoving = imresize(Imoving,size(Istatic),'bicubic');end% Make smooth images for histogram and fast affine registrationISmoving=imgaussian(Imoving,2.5,[10 10]);ISstatic=imgaussian(Istatic,2.5,[10 10]);% Detect if the mutual information or pixel distance can be used as % similarity measure. By comparing the histograms.if(isempty(Options.Similarity)) Hmoving= hist(ISmoving(:),60)./numel(Imoving); Hstatic = hist(ISstatic(:),60)./numel(Istatic); Hmoving(1)=0; Hstatic(1)=0; if(sum(log(abs(Hmoving-Hstatic)+1))>0.3), Options.Similarity='m'; if(Options.Verbose>0), disp('Multi Modalities, Mutual information is used'); drawnow; end else Options.Similarity='p'; if(Options.Verbose>0), disp('Same Modalities, Pixel Distance is used'); drawnow; end endendif(Options.Similarity(1)=='p'), type_affine='sd'; else type_affine='mi'; end% Register the moving image affine to the static image% Affine register the smoothed images to get the registration parametersif(strcmpi(Options.Registration(1),'R')) if(Options.Verbose>0), disp('Start Rigid registration'); drawnow; end % Parameter scaling of the Translation and Rotation scale=[1 1 1]; % Set initial affine parameters x=[0 0 0];elseif(strcmpi(Options.Registration(1),'A')) if(Options.Verbose>0), disp('Start Affine registration'); drawnow; end % Parameter scaling of the Translation, Rotation, Resize and Shear scale=[1 1 1 0.01 0.01 1e-4 1e-4]; % Set initial affine parameters x=[0 0 0 100 100 0 0];elseif(strcmpi(Options.Registration(1),'N')) if(Options.Verbose>0), disp('Start Affine part of Non-Rigid registration'); drawnow; end % Parameter scaling of the Translation, Rotation, Resize and Shear scale=[1 1 1 0.01 0.01 1e-4 1e-4]; % Set initial affine parameters x=[0 0 0 100 100 0 0];else warning('register_images:unknownoption','unknown registration method'); endfor refine_itt=1:2 if(refine_itt==2) ISmoving=Imoving; ISstatic=Istatic; end % Use struct because expanded optimset is part of the Optimization Toolbox. optim=struct('GradObj','off','GoalsExactAchieve',1,'Display','off','MaxIter',100,'MaxFunEvals',1000,'TolFun',1e-14,'DiffMinChange',1e-6); if(Options.Verbose>0), optim.Display='iter'; end x=fminlbfgs(@(x)affine_registration_error(x,scale,ISmoving,ISstatic,type_affine),x,optim); end% Scale the translation, resize and rotation parameters to the real valuesx=x.*scale;if(strcmpi(Options.Registration(1),'R')) % Make the rigid transformation matrix M=make_transformation_matrix(x(1:2),x(3));else % Make the affine transformation matrix M=make_transformation_matrix(x(1:2),x(3),x(4:5)); end% Make center of the image transformation coordinates 0,0[x,y]=ndgrid(0:(size(Imoving,1)-1),0:(size(Imoving,2)-1));xd=x-(size(Imoving,1)/2); yd=y-(size(Imoving,2)/2);% Calculate the backwards transformation fieldsBx = ((size(Imoving,1)/2) + M(1,1) * xd + M(1,2) *yd + M(1,3) * 1)-x;By = ((size(Imoving,2)/2) + M(2,1) * xd + M(2,2) *yd + M(2,3) * 1)-y;% Initialize the modality transformed image variablesM_TF=[]; F_TF=[]; % The nonrigid part of the registrationif(strcmpi(Options.Registration(1),'N')) % Demon registration parameters refinements=floor(log2(min(size(Imoving))/16)); if(refinements>MaxRef), refinements=MaxRef; end parameters.sigma_diff=Options.SigmaDiff; % Non-rigid registration if(Options.Verbose>0), disp('Start non-rigid demon registration'); drawnow; end % Do every refinements step twice if modality transformation enabled if(Options.Similarity(1)=='m'), loop=2; else loop=1; end % Loop trough all refinements steps. for j=0:refinements for l=1:loop % Set scaling parameters.resizepercentageentage resizepercentage=1/2^(refinements-j); if(resizepercentage>1), resizepercentage=1; end parameters.alpha=Options.Alpha*sqrt(resizepercentage); parameters.sigma_fluid=Options.SigmaFluid; if(Options.Verbose>0), disp(['Scaling resizepercentageentage : ' num2str(resizepercentage)]), end % Incase of multiple modalities, transform both images to their % opposite modalities. if(Options.Similarity(1)=='m') if(Options.Verbose>0), disp('Start modality transformation'); drawnow; end Bx_large=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1)); By_large=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2)); [Imoving_TF,Istatic_TF]=MutualTransform(Imoving,Istatic,15*sqrt(1/resizepercentage),4,Bx_large,By_large); if(Options.Verbose>0), disp('Finished modality transformation'); drawnow; end end sigma = 0.3/resizepercentage; % Set and resize the moving image and static image M=imresize(imgaussian(Imoving,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic'); F=imresize(imgaussian(Istatic,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic'); % Resize the modality transformed images if(Options.Similarity(1)=='m') M_TF=imresize(imgaussian(Imoving_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic'); F_TF=imresize(imgaussian(Istatic_TF,sigma,[sigma*6 sigma*6]),resizepercentage,'bicubic'); end % Resize the transformation field to current image size Bx=imresize(Bx,size(M),'bicubic')*(size(M,1)/size(Bx,1)); By=imresize(By,size(M),'bicubic')*(size(M,2)/size(By,2)); % Put transformation fields in x and y direction in one variable B=zeros([size(M) 2]); B(:,:,1)=Bx; B(:,:,2)=By; % Store the dimensions of transformation field, and make a long vector from T sizes=size(B); B=B(:); % Parameters options.sigma_fluid=parameters.sigma_fluid; options.sigma_diff=parameters.sigma_diff; options.alpha=parameters.alpha; options.interpolation=Options.Interpolation; % Optimizer parameters optim=struct('Display','off','StoreN',10,'GoalsExactAchieve',0,'HessUpdate','lbfgs','GradObj','on','OutputFcn', @store_transf,'MaxIter',200,'TolFun',1e-14,'DiffMinChange',1e-5); if(l==loop), optim.TolX = 0.02; else optim.TolX = 0.1; end if(Options.Verbose>1), optim.Display='iter'; end % Start the demon energy registration optimizer B=fminlbfgs(@(x)demons_energy(M,F,M_TF,F_TF,x,sizes,options),B,optim); % Reshape B from a vector to an x and y transformation field B=reshape(B,sizes); Bx=B(:,:,1); By=B(:,:,2); end end % Scale everything back if not already if(resizepercentage~=1) Bx=imresize(Bx,size(Imoving),'bicubic')*(size(Imoving,1)/size(Bx,1)); By=imresize(By,size(Imoving),'bicubic')*(size(Imoving,2)/size(By,2)); endend % Transform the input imageIreg=movepixels(Imoving,Bx,By,[], 3);if ( nargout>3 ) % Make the forward transformation fields from the backwards [Fx,Fy]=backwards2forwards(Bx,By);end% Set the class of output to input classif(strcmpi(Iclass,'uint8')), Ireg=im2uint8(Ireg); endif(strcmpi(Iclass,'uint16')), Ireg=im2uint16(Ireg); endif(strcmpi(Iclass,'uint32')), Ireg=im2uint32(Ireg); endif(strcmpi(Iclass,'int8')), Ireg=im2int8(Ireg); endif(strcmpi(Iclass,'int16')), Ireg=im2int16(Ireg); endif(strcmpi(Iclass,'int32')), Ireg=im2int32(Ireg); endif(strcmpi(Iclass,'single')), Ireg=im2single(Ireg); end% End time measurementif(Options.Verbose>0), toc, end
⛳️ 运行结果
🔗 参考文献
[1] 石跃祥,陈才.基于最优Atlas多模态图像的非刚性配准分割算法[J].光学学报, 2019, 39(4):11.DOI:10.3788/AOS201939.0410002.
[2] 苏孟超.基于SURF和光流场的多模态医学图像配准技术研究[D].南昌航空大学,2019.
[3] 王丽芳,王雁丽,史超宇,等.基于ZMLD与GC离散优化的非刚性多模态医学图像配准方法:CN201810563047.7[P].CN108711168A[2023-10-26].
