【图像配准】基于matlab互信息图像配准【含Matlab源码 1210期】

一、获取代码方式

获取代码方式1：
完整代码已上传我的资源：【图像配准】基于matlab互信息图像配准【含Matlab源码 1210期】

获取代码方式2：
通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。

备注：
订阅紫极神光博客付费专栏，可免费获得1份代码（有效期为订阅日起，三天内有效）；

二、互信息配准简介

1 互信息配准流程
信息论中将互信息定义为信息之间的关系，可以表示为两个随机变量之间统计相关性的度量，由此可以得出图像互信息的计算方法。作为图像多模态配准中的度量，图像互信息利用对图像灰度值的统计数据形成单个图像的灰度值概率函数和两个图像相似部分对应的灰度值联合概率函数，以此来衡量两幅图像的相关程度。在图像配准的过程中，认为两幅图像的相关性最大时对应的互信息也最大。基于互信息的图像配准最早在1995年被用于医学图像配准中，此后研究人员对于互信息的配准方法做了大量研究。

对于以互信息作为相似性测度的图像配准来说，目标就是寻找一个空间变换模型，使得两幅图像中的对应点相互对齐，以便进行比较或融合。基于互信息的图像配准流程如图1所示。

图1 基于互信息的图像配准流程图
它的流程可以概括为：首先提取两幅图像的特征空间信息，计算紫外图像与可见光图像之间的互信息，然后通过优化搜索算法对空间变换矩阵进行更新迭代，从中找出使两幅图像相似性测度互信息最大的最优变换矩阵T，最后根据变换矩阵实现对紫外图像的空间变换，使其在空间上与可见光图像配准，达到图像融合的要求。

2 图像互信息相关理论
基于图像灰度信息的相似性测度是一种通过统计和概率估计的方法对图像进行匹配的表示相似程度的参数。常见的图像灰度信息相似性测度包括互相关法、均方差和法及互信息法等。在多模态图像中最常用的为互信息法，互信息法是一种不需要对源图像做预处理的方法，被广泛应用于医学图像配准中。互信息法对于各种模式的图像配准都具有有效性，鲁棒性及可靠性较高。本文采用基于互信息的相似性测度来提高紫外与可见光图像的匹配精度，下面介绍一些图像互信息的基本概念。
1）熵
在信息论中，熵表示随机变量的不确定性。假设随机事件xi出现的概率为pi，则信息量为：I(xi)=-log 2pi，单位为bit。而信息量的数学期望为信源的平均信息量，即信源的信息熵。
假设一个随机变量X内存在N个随机事件X={x1,x2，…，xn}，随机事件出现的概率集合为P={p1,p2，…，pn}，也就是随机变量X的概率密度函数为p(x），则随机变量X的熵可以定义为

此外联合熵表示两个随机变量的相关性。联合熵越小表示这两个随机变量越相关，若两个随机变量越不相关，它们的联合熵会越大。假设pXY(x,y）是两个随机变量X,Y的联合概率密度函数，则联合熵的数学定义为

2）图像互信息
图像互信息可以度量两幅图像的统计独立程度，基于前面熵和联合熵的定义，给定紫外与可见光图像UV和V，他们的互信息可以定义为

或

式（3）中，MI(V,UV）为可见光与紫外光图像的互信息，H(V),H(UV）和H(V,UV）分别表示可见光图像的熵、紫外光图像的熵、可见光与紫外光图像的联合熵。式（4）中，p(v,uv）表示可见光与紫外光图像的联合概率密度函数，p(v）和p(uv）分别表示可见光图像的概率密度函数和紫外光图像的概率密度函数。根据互信息的定义，对于可见光与紫外光的互信息计算流程如图2所示。

图2 紫外光与可见光图像互信息计算流程
可见光图像和紫外光图像的像素灰度值分别表示为LV(i,j),LUV(i,j):

其中，i=1,2，…，m,j=1,2，…，n，其中i,j分别表示像素的横纵坐标，m,n分别表示图像的宽高像素数。先将可见光图像和紫外图像对应坐标的像素值整合成一个二维向量（LV(i,j),LUV(i,j)),i=1,2，…，m,j=1,2，…，n。然后通过直方图估计法估计出两个图像的联合概率密度p(v,uv）和边缘概率密度p(v）=∑p(v,uv),p(uv）=∑p(v,uv）。最后，根据式（4）计算出两幅图像的互信息。

三、部分源代码

close all
clear allticx0=[50; 50; -15; 0.5]; % Initial points, [X, Y, angle, scaling] % Select them as close to the matching points as possible% by guessing; otherwize, it will fail.[x, fval]=fminsearch(@image_registr_MI,x0) % Optimization using 'fminsearch'%
% Display
%
load imageIM1=double(IM1);
IM2=double(IM2);
IM2=imresize(IM2, x(4), 'bilinear');
J=imrotate(double(IM2), x(3),'bilinear'); %rotated cropped IMAGE2
function f=image_registr_MI(x)load imageIM1=double(IM1);
IM2=double(IM2);
IM2=imresize(IM2, x(4), 'bilinear');
J=imrotate(double(IM2), x(3),'bilinear'); %rotated cropped IMAGE2if n1>n3-x(1)/2f=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;
elseif x(1)>n3-n1x(1)=n3-n1-1;IM1(1:n1, 1:n2)=255;endif x(2)>n4-n2x(2)=n4-n2-1;IM1(1:n1, 1:n2)=255;endif x(1)<0x(1)=0;IM1(1:n1, 1:n2)=255;endif x(2)<0x(2)=0;IM1(1:n1, 1:n2)=255;endxt=1:n1;yt=1:n2;xx=round(xt+x(1));yy=round(yt+x(2));IM2=round(J(xx, yy)); % selecting part of IMAGE2 matching the size of IMAHE1rows=size(IM1,1);cols=size(IM2,2);N=256;h=zeros(N,N);for ii=1:rows;    %  col for jj=1:cols;   %   rowsh(IM1(ii,jj)+1,IM2(ii,jj)+1)= h(IM1(ii,jj)+1,IM2(ii,jj)+1)+1;endend[r,c] = size(h);b= h./(r*c); % normalized joint histogramy_marg=sum(b); %sum of the rows of normalized joint histogramx_marg=sum(b');%sum of columns of normalized joint histogranHy=0;for i=1:c;    %  colif( y_marg(i)==0 )%do nothingelseHy = Hy + -(y_marg(i)*(log2(y_marg(i)))); %marginal entropy for image 1endendHx=0;for i=1:r;    %rowsif( x_marg(i)==0 )%do nothingelseHx = Hx + -(x_marg(i)*(log2(x_marg(i)))); %marginal entropy for image 2end   endh_xy = -sum(sum(b.*(log2(b+(b==0))))); % joint entropyf=-(Hx+Hy-h_xy);% Mutual information%x
end

四、运行结果

五、matlab版本及参考文献

1 matlab版本
2014a

2 参考文献
[1] 蔡利梅.MATLAB图像处理——理论、算法与实例分析[M].清华大学出版社，2020.
[2]杨丹,赵海滨,龙哲.MATLAB图像处理实例详解[M].清华大学出版社，2013.
[3]周品.MATLAB图像处理与图形用户界面设计[M].清华大学出版社，2013.
[4]刘成龙.精通MATLAB图像处理[M].清华大学出版社，2015.
[5]侯思祖,陈宇,刘雅婷.基于互信息的紫外成像仪中图像配准研究[J].半导体光电. 2020,41(04)