图像配准中的互信息

声明

本文主体部分参考 胡永祥 基于互信息的多模态医学图像非刚性配准研究[D] 中南大学 2012

对其中部分不清晰及符号使用问题进行了调整。

旨在对图像配准过程中使用的互信息及其梯度计算进行大致了解。

疏漏之处敬请指出。

熵和互信息的基本概念

嫡和互信息都是信息论中的重要概念。熵在信息系统屮作为事物不确定性的表征，反映了系统所包含的信息总量。五信息用于表示信息之间的关系，是两个随机变量统计相关性的测度。

香农熵

消息中所含的信息量是与消息所描述的事件的概率或不确定性有关。消息所表达的事件的概率越小，其包含的信息量就越大，反之，其信息量就越小。如果某事件发生的概率为1，则其包含的信息量为0。如果某事件发生概率为0，则消息含有无限的信息量。设某事件En的发生概率为Pn，则其包含的信息量为:

H(En)=−logaPn

其中，a为对数的底。通常取a的值为2，这时H的单位为比特(bit)。当要描述事件集合的平均信息量时，就需引入用香农熵。设事件集合E=(E1,E2,...,EN)的发生

概率分别为(p1,p2,...,pn)，则E所能提供的平均信息量即为香农熵，定义为：H(E)=−∑i=1npilogapi

香农熵的一些数学性质：

非负性，H(E)≤0，当且仅当描述单个概率为1的事件时其值为0；

当所有事件的发生概率相等时，即当pi=1/N,i=1,2,...,N时熵达到最大值logaN；

熵是多个事件自信息量的加权和，因此概率分量的位置互换，并不影响熵的值；

熵是概率分布(p1,p2,...,pn)的凸函数；

联合熵与条件熵

联合熵是检测两个随机变量之间相关性的一种统计量。设p(x,y)为随机变量X、Y的联合概率密度，则联合熵定义为：

H(X,Y)=−∑xyp(x,y)logap(x,y)

两个随机变量越相关，它们的联合熵的值就越小；反之，联合熵的值就越大。联合熵大于等于任一个变量的独立熵，即H(X,Y)≥max(H(X),H(Y))，联合熵小于独立熵的和，即H(X,Y)<H(X)+H(Y)

条件熵用来衡量在已知随机变量X的前提下，随机变量Y的信息量。其定义为：

H(Y|X)=H(X,Y)−H(X)

当根据X就能完全确定Y时，条件熵为0，当X，Y相互独立时H(Y|X)=H(Y)

互信息

对于两个变量X，Y来说，互信息可以定义为：

MI(X,Y)=H(X)+H(Y)−H(X,Y)

根据条件熵的定义，可以得出：

MI(X,Y)=H(X)+H(Y)−H(X,Y)=H(X,Y)−H(X|Y)−H(Y|X)=H(X)−H(X|Y)=H(Y)−H(Y|X)

将熵的定义带入公式，对数以2为底，可得：

MI(X,Y)=∑xyp(x,y)log2p(x,y)−∑xp(x)log2p(x)−∑yp(y)log2p(y)=∑x,yp(x,y)log2p(x,y)p(x)p(y)

式中p(x,y)为随机变量X、Y的联合概率密度函数，p(x),p(y)分别为变量X、Y的边缘概率密度函数。

互信息有如下几个性质：

对称性：MI(X,Y)=MI(Y,X);

非负性：MI(X,Y)>0;

极值性：MI(X,Y)≤H(X);

凸函数性；

图像互信息

根据互信息的定义，可以得出在图像上的互信息的定义：

给定两幅图像F和M，它们的像素灰度分别为fij,mij,i=1,...,m,j=1,...,n，其中m、n表示图像的行和列的大小。计算互信息时，将两个图像的对应坐标位置处的像素构成一个二维向量(fij,mij),i=1,...,m,j=1,...,n。利用这个二维向量组估计出两幅图像的联合概率密度pF,M(f,m)和边缘概率pF(f)=∑mpF,M(f,m)，pM(m)=∑fpF,M(f,m)，最后计算互信息。

估计两个图像的互信息的关键在于两个图像的联合概率密度和边缘概率密度。概率密度估计的方法可以分为参数估计法和非参数估计法两类。参数估计法是指概率密度函数的形式已知，而函数的参数未知，通过训练数据来估计参数的方法，如最大似然估计，贝叶斯(Bayes)估计。非参数估计法对密度函数的形式不作任何假设，而是直接利用数据对概率密度进行估计，如核密度估计法(Parzen法)，k-近邻法。由于不同图像的灰度数据差别很大，无法预知其概率密度形式，因此常用非参数估计法。最简单的估计法是直方图估计法，利用该方法可得到离散的概率密度。由于在刚性配准中通常使用Powell最优化算法，该最优化方法不需求目标函数的梯度，因此该方法在刚性配准中得到广泛应用。在非刚性图像配准中，最常用的概率密度估计法为核密度估计法。

图像概率密度估计

直方图法和核密度估计法是最常用的两种概率密度估计法。直方图估计法原理简单、计算速度快，但估计误差较大且只能得到离散的概率密度。核密度估计法比较准确，能得到连续的概率密度，但计算量大。

直方图概率密度估计法是最简便的密度估计法,

直方图估计法

直方图概率密度估计法是最简便的密度估计法，它用灰度出现的频率来近似表示灰度的概率。

核密度估计法

定义1 设x1,x2,...,xn为取值于R的独立同分布随机变量，其所服从的分布密度函数为f(x),x∈R。定义函数：

f^h(x)=1nh∑i=1nK(xi−xh)

为概率密度函数f(x)的核密度估计，K(⋅)称为核函数。为方便，记Kh(u)=K(u/h)/h，则上式可表示为：

f^h(x)=1n∑i=1nKh(xi−x)

核密度估计不仅与样本点集有关，还与核函数的选择和箱宽h的选择有关。常用的核函数有以下几种：

高斯函数

K(x)=12π−−√exp(−x22

Epanechnikov核函数

K(x)={3(1−x2)20|x|≤1|x|>1

Biweight核函数

K(x)=⎧⎩⎨⎪⎪15(1−x2)2160|x|≤1|x|>1

基于样条的核函数

K(x)=⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪(x+2)3−x3−2(x+1)3+6(x+1)x3+2(x−1)3−6(x−1)(2−x)30−2≤x<−1−1≤x<00≤x<11≤x<2else

定义 设X1,X2,...,Xn为取值独立于Rd的独立分布随机变量，器所服从的分布密度函数为f(x),x∈R。定义函数：

f^x=1n∑i=1nKd,h(xi−x)

其中

Kd,h(X)=∏i=1dKhi(xi)

特别的，采用d维高斯函数构造的Rd上的核密度函数具有如下形式：

Kd,σ(X)=(2π−−√σ)−d∏i=1de−(x−xi)22σ2

基于核密度估计法的图像概率密度估计

设图像I的像素为rij,i=1,...,N,j=1,...,M，其灰度概率密度函数的估计可表示为：

p(r)=1MN∑i=1M∑j=1NKh(rij−r)

如果选用高斯函数作为核函数，则可表示为：

p(r)=1MN2π−−√σ

方差的选取对基于高斯核函数的概率密度函数估计有较大的影响。(此处尚未完全理解，按下不表，日后填坑)

根据核密度估计法,很容易的将单个图像概率密度估计法推广到两个图像的联合概率密度估计。设参考图像F和浮动图像M构成的二维联合向量为(fij,mij),i=1,...,M,j=1,...,N则联合概率密度的核估计法可表示为：

p^(f,m)=1MN∑i=1M∑j=1NKhf(f−fij)Khm(m−mij)

其中hf,hm分别表示核函数的箱宽（带宽？不清楚）采用高斯函数作为核函数时可以表示为：

p^(f,m)=1MN∗2πσ1σ2∑i=1M∑j=1Ne−((f−fij)22σ21+(m−mij)22σ22)

基于快速高斯变换的图像互信息估计(先挖坑，以后填)

快速高斯变换

快速高斯变换最早由Grecngard等提出,用于从N个点的数据空间中快速估计M个点的高斯函数的值。

设给定一个N点的数据集{sj}j=1,...,N，用这个数据集来估计M个点的高斯函数可以表示为：

G(ti)=∑j=1Nfje−(ti−sj)2σ2,i=1,...,M

其中fi是加权系数，σ为高斯函数的带宽（箱宽？方差？），{ti}i=1,...,M为目标点。

互信息梯度

非刚性配准通常是一个最优化过程。为找到最优的配准参数，最优化算法的每一次迭代都需计算目标函数的梯度或Hessian矩阵。

数学定义

为便于计算梯度,我们假设图像为连续图像并引入了几何形变参数,因此在这里先对图像互信息进行重新定义。令fF(x),fM(x)分别表示参考图像和浮动图像，x∈Ω∈Rd，d表示图像维数，Ω表示图像空间。g(x;μ1,μ2,...)表示几何变换，μ=(μ1,μ2,...)表示形变参数。LF,LM分别表示参考图像和浮动图像的灰度集，l∈LF,k∈LM。则使用核函数估计法的两个图像的联合概率密度可表示为：

p(l,k;μ)=α∑x∈ΩKhF(l−fF(x))⋅KhM(k−fM(g(x;μ)))

其中α为归一化因子，确保∑f∈lF∑m∈LM(f,m;μ)=1。

边缘概率密度为：

pF(l;μ)=∑kp(l,k;μ)

pM(k;μ)=∑lp(l,k;μ)

其中，f∈LF,m∈LM。

图像互信息可记为：

MI(F,M;μ)=−∑l∈LF∑k∈LMp(l,k;μ)log2p(l,k;μ)pF(l;μ)pM(k;μ)

当互信息作为相似性测度时,我们将其记为

SMI(F,M;μ)=MI(F,M;μ)

几种互信息相似性测度的梯度计算方法

有限差分法(finite difference method)

有限差分法定义为：

∂SMI(T,F;μ)∂[μ]i≈SMI(F,M;μ+Δei)−SMI(F,M;μ−Δei)2Δ

其中，∂[μ]i表示对参数μ的第i个分量求偏导，Δ为小的标量值，ei是除了第i个分量为1外其他位置都为0的单位向量。

基于参数的解析法(analytic approach)

互信息的梯度可以表示为：

∇SMI(F,M;μ)=[∂SMI∂μ1,∂SMI∂μ2,...,∂SMI∂μi,...]

其中∂SMI∂μi,i=1,2...表示梯度向量的第i个分量，可表示为：

∂SMI∂μi=∂∂μi(−∑l∈LF∑k∈LMp(l,k;μ)log2p(l,k;μ)pF(l;μ)pM(k;μ))=(−∑l∈LF∑k∈LM∂∂μip(l,k;μ)log2p(l,k;μ)pF(l;μ)pM(k;μ))=−∑l∈LF∑k∈LM∂∂μi(p(l,k;μ)(log2p(l,k;μ)−log2pF(l;μ)−log2pM(k;μ)))=−∑l∈LF∑k∈LM(∂p(l,k;μ)∂μilog2p(l,k;μ)pF(l;μ)pM(k;μ)−ln2∗p(l,k;μ)∗(1p(l,k;μ)∂p(l,k;μ)∂μi−1pF(l;μ)∂pF(l;μ)∂μi−1pM(k;μ)∂pM(k;μ)∂μi))=−∑l∈LF∑k∈LM(∂p(l,k;μ)∂μi(log2p(l,k;μ)pF(l;μ)pM(k;μ)+ln2)+1pF(l;μ)∂pF(l;μ)∂μi+1pM(k;μ)∂pM(k;μ)∂μi)

在上式中，可以注意到：

pF(l;μ) 事实上与μ无关，所以1pF(l;μ)∂pF(l;μ)∂μi=0

∑l∈LF∑k∈LMp(l,k;μ)pM(k;μ)∂pM(k;μ)∂μi=∑k∈LM∂pM(k;μ)∂μi∑l∈LFp(l,k;μ)pM(k;μ)=∂∂μi(∑k∈LMpM(k;μ))=0

∑l∈LF∑k∈LMp(l,k;μ)∂μi=∂∂μi(∑l∈LF∑k∈LMp(l,k;μ))=0

所以有：

∂SMI∂μi=−∑l∈LF∑k∈LM∂p(l,k;μ)∂μilog2p(l,k;μ)pM(k;μ)

根据核变换下联合概率密度计算公式有：

∂p(l,k;μ)∂μi=α∑x∈ΩKhF(l−fF(x))(dKhMdξ∣∣∣ξ=k−fM(g(x;μ)))(−dfMdζ∣∣∣ζ=g(x;μ))T∂g(x;μ)∂μi

其中dfMdζ表示浮动图像在形变后位置ζ=g(x;μ)的梯度。当空间形变采用度为n的B-spline时，有：

dfMdζ∣∣∣ζ=g(x;μ)=∑xi∈Ωc(xi)dβn(μ)dμ∣∣∣u=t−xi=∑xi∈Ωc(xi)⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢dβn(μ)dμ∣∣∣u=(t−xi)1βn((t−xi)2)...βn((t−xi)1)∂βn(μ)∂μ∣∣∣u=(t−xi)2...⋮⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥

参考Optimization of mutual information for multiresolution image registration一文所写，∂g(x;μ)∂μi描述某位置关于形变参数的变化，是一个紧依赖于几何学的项（问题是到底怎么求啊！！！）。

参考文献

1 胡永祥 基于互信息的多模态医学图像非刚性配准研究[D] 中南大学 2012

2 Thévenaz P, Unser M. Optimization of mutual information for multiresolution image registration[J]. Image Processing, IEEE Transactions on, 2000, 9(12): 2083-2099.