CN103700090A - 一种基于各向异性热核分析的三维图像多尺度特征提取方法 - Google Patents

Abstract

本发明提出一种基于各向异性热核分析的三维图像多尺度特征提取方法,包括如下步骤：（1）借助于黑塞矩阵的特征值和特征向量，设计一种三维图像局部几何结构的向量化表示方法；（2）基于定义的向量场，在三维图像上数学严格地定义了一种各向异性热传递描述模型；（3）基于定义的各向异性热传递模型，提出用黑塞矩阵的第一特征值来定义热量场；（4）通过引入特征保持的下采样方法，设计了一种由粗到细的层次化的热核加速策略；本发明的算法最核心的思想是：通过将局部和全局特征结构融入到各向异性热传递的理论框架来实现三维图像的多尺度特征提取。

Description

一种基于各向异性热核分析的三维图像多尺度特征提取方法

技术领域

本发明涉及基于各向异性热核分析的三维图像多尺度特征提取方法，具体涉及一种基于热核理论的三维医学图像分析模型，通过将图像的局部和全局谱分析融入到热传递的理论框架中，实现鲁棒的多尺度点特征提取。

背景技术

三维医学影像数据分析的数学建模方法涉及到热核理论应用、图像特征提取和基于热扩散的图像处理三个方面。

在热核理论应用方面：近年来，物理上的热核和热扩散理论在计算机图形图像处理领域受到了众多研究者的青睐，已在二维流形的多尺度特征提取、平滑处理以及形状检索等方面得到了广泛应用。其中，最重要的原因是：热核天然地蕴含了多尺度的性质，并具有等距变换不变性，这对图形或图像的鲁棒特征提取具有重要意义。如：一种尺度不变的热核信号计算方法被提出，并被应用于形状检索，取得了较好的效果。同时，为了提高在大规模网格模型上进行热核计算的时间效率，提出了一种多分辨率的计算策略来加速热核信号的计算。

在图像的特征提取方面：图像的局部特征可被看作是与周围领域具有显著不同、并表征特定物理意义的一种模态，他们一般蕴含丰富的信息并可通过一些微分量来进行验证。现有的以Harris、Hessian、LOG以及DOG等算子为核心的图像特征提取算法的主要思想是：通过对图像的微分属性分析，提取某种物理量的局部极值，并将其作为特征。因此，这些算法一般对尺度变化、噪声和图像形变等因素都比较敏感。将Harris算子和Hessian算子相结合并辅以尺度选择，是提高图像特征提取鲁棒性的一个可行的策略，而且该方法已被推广到了三维情形。应该说，目前2D SIFT算法仍是最为稳定的二维图像特征提取算法，而且它的高维形式的扩展已经开始出现。

在基于热扩散理论的图像处理方面：在过去的20年中，借助高斯核函数，热扩散理论在图像的多尺度分析领域被广为使用。但是，高斯核一般只能一个很小的邻域内刻画热扩散过程，因为其本质上等价于用各向同性的热传递系数作为权重对原始图像进行卷积处理。自从双边滤波方法被提出来用于边界保持的图像平滑和去噪处理以后，基于各向异性热传递理论的滤波器才开始在二维图像处理领域逐渐流行起来。特别值得指出的是，在二维流形上首先提出的热核信号已被推广用于三维流形，但是这种扩展尚没有考虑三维流形内部材质的差异性。这也是启发我们直接在三维医学图像上研究各向异性热核的一个主要因素。

发明内容

热核理论在流形网格模型建模分析领域已取得引人瞩目的新进展，将其扩展，并与谱图分析理论相结合，在三维医学影像数据建模分析方面具有重要的研究价值。本部分研究内容将热核理论和谱图分析理论进行优势结合，探索将三维医学数据从几何空间非线性映射到特征空间的数学模型，为进而在特征空间对原始医学影像数据进行分析奠定基础。

本发明提出了一种基于各向异性热核分析的三维图像多尺度特征提取方法，包括如下步骤：

（1）借助于黑塞矩阵的特征值和特征向量，设计一种三维图像局部几何结构的向量化表示方法；

（2）通过将特征向量、特征值和上述局部几何结构的比率值进行组合，来定义一个高维的向量场；借助三维图像的全局拉普拉斯矩阵，该向量场的向量化表示将被用于各向异性热传递模型的构建；基于定义的向量场，在三维图像上数学严格地定义了一种各向异性热传递描述模型；

（3）基于定义的各向异性热传递模型，提出用黑塞矩阵的第一特征值来定义热量场；

（4）通过引入特征保持的下采样方法，设计了一种由粗到细的层次化的热核加速策略。

附图说明

图1为算法整体处理流程示意图。

图2不同特征值的不同取值范围所对应的图像。

图3为特征保持的下采样结果。

图4为人头MRI图像特征提取结果及其与3D SIFT方法的效果对比。

图5为人体盆腔CT图像的特征提取结果及其与3D SIFT算法的效果对比。

图6为人体腹部CT图像的特征提取结果及其与3D SIFT算法的效果对比。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。具体实施如下：

1.三维图像的局部和全局结构分析

1.1三维图像局部结构的向量化

黑塞矩阵H是刻画三维图像局部几何结构的一种有效方式，它能够很好的描述表面法向量沿某一等值面的变化情况。作为一种实值的对称矩阵，H具有实数特征值，并且满足|λ1|≥|λ2|≥|λ3|，其相应的特征向量可被记为：e1,e2,e3。从某种意义上说，最大特征值所对应的特征向量表征了从一种材质过度到另一种材质时最为迅速的变化方向，而最小特征值所对应的方向则指明了某种材质在图像内部流动的方向。并且，最大和最小特征值恰恰蕴含了图像的主曲率信息。因此，从各向异性热传递的角度来说，同时考虑黑塞矩阵的特征值和特征向量无疑可以克各向同性滤波器的一些局限性。同时，提取三维图像的局部结构并对三维图像内部的材质属性变化进行定量刻画也正是我们的目的所在。Alfiansyah等人通过将黑塞矩阵的特征值进行不同形式的组合运算，分析了它们与图像局部结构（平面、块和噪声）之间的关系：

R_sheet＝(|λ₃|+|λ₂|)/|λ₁|

R_blob＝(2|λ₁|-|λ₂|-|λ₃|)/2|λ₁|

$R_{\mathrm{noise}}=\sqrt{{\mathrm{\λ}}_3^2+{\mathrm{\λ}}_2^2+{\mathrm{\λ}}_1^2}$

公式中的λ₁、λ₂、λ₃都是管状曲面上的体素的特征值。对位于管状曲面上的体素，它的特征值具有如下规律：|λ₁|=0and|λ₂|=|λ₃|>>0，因此，我们也可以类似地定义一种特征值的比率来描述这类几何结构：

R_tube＝(|λ₃|+|λ₂|-|λ₁|)/|λ₁|

公式中的λ₁、λ₂、λ₃都是管状曲面上的体素的特征值。当它们所代表的结构出现时，上述4种特征值的比率值会相应的变为0。据此，我们可以通过将特征向量、特征值和上述描述局部几何结构的比率值进行组合，来定义一个高维的向量场，其向量化表示形式可被定义为：ELSV R=(e₁ ^T,e₂ ^T,e₃ ^T,λ₁,λ₂,λ₃,R_sheet,R_blob,R_tube,R_noise)。借助三维图像的全局拉普拉斯矩阵，ELSVR将被用于各向异性热传递模型的构建。

1.2基于局部结构变化加权的全局拉普拉斯算子

在ELSVR所定义的流形上，热的各向异性传播可由图的全局拉普拉斯算子来决定。若将三维图像定义为一个无向图G=(V,E)，其中V代表所有的体素，E满足E∈V×V，并且当且仅当相应体素满足一环领域要求时，边E_ij∈E的定义才成立（这里我们将某个体素的邻域体素记为N_i）。为了诱使热沿相同的材质传播，而抑制其跨材质传播，以ELSVR空间的距离度量作为边的权重来对图的全局拉普拉斯算子进行计算。因此，首先需要度量两个相邻体素V_i和V_j在其所属材质类型方面的相似度：

$S_{i,j}=\left\{\begin{array}{cc}0,& \mathrm{if}{R}_{\mathrm{noise}}^i=0\mathrm{or}{R}_{\mathrm{noise}}^j=0\\ e^{-d_r},& \mathrm{otherwise}\end{array}\right.$

这里d_r=ΣR_k∈R(Rⁱ _k-R^j _k)²，而R={R_sheet,R_blob,R_tube}，权重W(I,j)可被定义为：

$W_{i,j}=\left\{\begin{array}{cc}{e}^{\frac{-(1-S_{i,j})d_{i,j}{d}_{i,j}^{\mathrm{\α}}}{{\mathrm{\σ}}^2}},& V_j\∈N_i\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中d_i,j代表V_i和V_j间的欧氏距离d_x,d_y,d_z代表图像采集时的采样间距，d^σ _i；j=(rcos(eⁱ ₁·e^j ₁)+rcos(eⁱ ₂·e^j ₂)+rcos(eⁱ ₃·e^j ₃))，而σ则是一个控制参数，它正比于两个像素间结构距离的方差。因此，图的全局拉普拉斯矩阵可被定义为：

$L_{i,j}=\left\{\begin{array}{cc}{W}_{i,j},& V_j\∈N_i\\ -{\mathrm{\Σ}}_{j\∈N_i}{W}_{i,j},& i=j\\ 0,& \mathrm{otherwise}\end{array}\right.$

三维图像的各向异性热核算子h^v _t(V_i,V_j)(VAHK)可使用拉普拉斯矩阵的有限个特征向量来进行近似。由于ELSVR加权的图的拉普拉斯算子自然的将图像局部结构的差异性纳入其中，因此，相应的热核算子更具材质和结构敏感性。

2.三维图像的各向异性热核定义

2.1热核算子

在高维ELSVR流形上的热传递过程可用一个热控制方程来描述。假设(M；μ)是一个具有度量μ的紧支黎曼流形，在(M；μ)上面的热扩散过程由如下的方程控制：

$\frac{\∂u(x,t)}{\∂t}=-{\mathrm{\Δ}}_{\mathrm{\μ}}u(x,t),$

其中，△μ表示定义在流形(M；μ)上的拉普拉斯贝尔特拉米算子,而h_t(x,y):R⁺×M×M-→R⁺则被称作热核。并且对所有的f∈L²(M；μ),t>0,x∈M，满足：

H_tf(x)＝∫_Mh_t(x，y)f(y)dμ(y).

实质上，热核函数h_t是对算子H_t的一种积分。假设M的拉普拉斯贝尔特拉米算子△μ能够被特征分解，并具有特征系统

那么热核函数则可被定义为：

$h_t(x,y)=\underset{k=1}{\overset{\∞}{\mathrm{\Σ}}}{e}^{-t{\mathrm{\λ}}_k}{\mathrm{\φ}}_k\left(x\right){\mathrm{\φ}}_k\left(y\right).$

从物理意义的角度来说，若在x位置给定一个单位脉冲热源，热核描述了在t时刻热从x到y的传播量。热核算子具有对称性、等距不变性、多尺度、信息丰富和稳定性，并且它同时满足半群特性：h_t+s(x,y)=∫_Mh_t(x,z)h_s(y,z)dz。

2.2初始热量场定义

最简单的方式是直接使用图像的颜色信息来定义初始热量场，但是这种方法会忽略掉ELSVR中所蕴含的结构信息。如图2所示，ELSVR中的特征值可用来进行特征敏感的初始热量场定义，因为：（1）第一特征向量代表了最大曲率的方向；（2）特征值是旋转不变的；（3）第一特征值和第三特征值在图像的边界区域往往呈现出一种负对称性，而通过第二特征值所提取的图像对比度较差。因此，我们选择第一特征值来定义初始热量场。

然后，以各向异性热核函数作为卷积核对初始热量场进行多尺度卷积，进而通过定位跨尺度的极值点来进行特征提取。而初始热量场f₀(λⁱ ₁)在经过一次卷积后，将变为：

$f_t\left({\mathrm{\λ}}_1^i\right)=\underset{V_j\∈{\mathrm{\Ω}}_i}{\mathrm{\Σ}}{h}_t^v(V_j,V_i)f_0\left({\mathrm{\λ}}_1^j\right),$

其中，Ωi是体素Vi所对应的卷积窗口。将公式以矩阵的形式进行重写，我们可以得到F_t(λ₁)=A_tF₀(λ₁)。通过应用热核函数的半群性质，经过n次卷积后的热量场可用数学公式表达为F_nt(λ₁)=A_t0F_(n-1)t(λ₁)=…=A^n-1 _t0F_t(λ₁)=A_nt0F₀(λ₁)。将t=t₀时刻的热量场看做初始热量场经一步扩散后的结果，那么随着n的增加，热量的扩散结果可通过将其和一个稀疏矩阵的乘法结果来近似，因此，为了提高计算效率，这个卷积过程可用基于CUDA的矩阵乘法来加速。

2.3层次化的特征提取策略

对于三维图像来说，最大挑战源于其巨大的数据量给热核计算所带来的极高的计算复杂度。实际上，多分辨率分解方法非常适合用来对此进行加速。其核心思想是：首先将原始图像进行下采样并进行类似金子塔形式的层次化组织，先在数据规模较小的低分辨率图像上进行热核计算和特征点位置的粗略定位，然后在原始图像相应的局部区域内，进行精细的热核计算，进而最终找出特征点的精确位置。然而，由于下采样会不可避免的导致部分信息的的损失，因此我们提出了一种特征保持的下采样方法，它可被定义为：

其中，

是下采样后的图像灰度值，Ω_i是

所对应的原始图像的下采样区域。图3给出了特征保持的下采样方法与平均下采样方法之间的效果对比。从中可以看出，我们的方法对局部曲率信息更为敏感，并且在计算过程中可以重用ELSVR中业已计算好的曲率信息。在下采样后，对于各向异性热核的计算，我们使用拉普拉斯矩阵最小的300个特征值及其相应的特征向量来进行近似计算。热核信号对大尺度的特征提取非常鲁棒，而下采样的三维图像恰恰可以较好的保留原始图像中较大尺度的特征信息。因此，下采样后三维图像的热核信号(VAHKS)的计算公式为：

$h_t^v(V_x,V_x)=\underset{k=1}{\overset{n}{\mathrm{\Σ}}}{e}^{-t{\mathrm{\λ}}_k}{\mathrm{\φ}}_k\left(V_x\right){\mathrm{\φ}}_k\left(V_x\right).$

一个比较重要的参数是时间t，它直接与热传递区域的范围直接相关。在下采样的图像中，最小的时间t的取值被设置为D_size×min(d_x,d_y,d_z)，而D_size代表下采样间距。通过给定不同的t值，我们可以得到多个尺度下的VAHKS，因此，可在该尺度空间通过查找跨尺度的极值点来定位特征点的大致位置。然后，再在原始图像中对这些特征点的位置进行细化。在实际实现过程中，我们采用Matlab的eigs函数来执行拉普拉斯矩阵的特征分解，根据我们的经验，当矩阵的规模超过40,000×40,000时，这个分解过程将变得异常耗时。此外，由于我们所采用的最大的下采样间距为16倍的体素间距，因此，对于在下采样图像中所定位的每一个特征点，在原始图像中相应的选取它的一个Ω_refine=32×32×32邻域来重新进行热核计算，足以满足特征点求精计算的需要。然而，还有一个需要说明的问题是局部热扩散区域的大小与扩散时间之间的关系。由于热核h(V,·)可被看作是布朗运动的概率传递函数，而一个粒子发生布朗运动的影响范围可用一个测地球

来限定，因此，当原始图像中的某个体素Vi作为热源时，我们可以定义它的热传递区域为：

${\mathrm{\Ω}}_t^i=\left\{\begin{array}{ccc}{V}_j|V_j\∈{\mathrm{\Ω}}_{\mathrm{refine}}& \mathrm{and}& h_t^v(V_j,V_i)>C\left(t\right)\end{array}\right\},$

其中C(t)是一个阈值常数，可被设定为

所以，鉴于上述分析，在边界附近所提取的特征点是不稳定的，应当被舍弃掉。在我们的试验中，对于常规规模的三维图像来说，在舍弃不稳定的点后，我们一般可以提取100-1500个稳定的特征点，这足以满足配准之类的基于特征的图像处理方面的应用。

Claims (5)

1.一种基于各向异性热核分析的三维图像多尺度特征提取方法，其特征是包括以下步骤:

（1）借助于黑塞矩阵的特征值和特征向量，设计一种三维图像局部几何结构的向量化表示方法；

（2）通过将特征向量、特征值和上述局部几何结构的比率值进行组合，来定义一个高维的向量场；借助三维图像的全局拉普拉斯矩阵，该向量场的向量化表示将被用于各向异性热传递模型的构建；基于定义的向量场，在三维图像上数学严格地定义了一种各向异性热传递描述模型；

（3）基于定义的各向异性热传递模型，提出用黑塞矩阵的第一特征值来定义热量场；

（4）通过引入特征保持的下采样方法，设计了一种由粗到细的层次化的热核加速策略。

2.根据权利要求1所述的基于各向异性热核分析的三维图像多尺度特征提取方法，其特征是：所述的借助于黑塞矩阵的特征值和特征向量中，设计了一种三维图像局部几何结构的向量化表示方法，并且向量空间的距离度量用于表征三维图像材质连续性；黑塞矩阵H是刻画三维图像局部几何结构的一种有效方式，它能够很好的描述表面法向量沿某一等值面的变化情况；作为一种实值的对称矩阵，H具有实数特征值，并且满足：

|λ1|≥|λ2|≥|λ3|，λ_i分别为矩阵特征值，

其相应的特征向量可被记为：e1,e2,e3；从某种意义上说，最大特征值所对应的特征向量表征了从一种材质过度到另一种材质时最为迅速的变化方向，而最小特征值所对应的方向则指明了某种材质在图像内部流动的方向，并且，最大和最小特征值恰恰蕴含了图像的主曲率信息，通过将黑塞矩阵的特征值进行不同形式的组合运算，分析了它们与图像局部结构之间的关系，所述局部结构是指平面、块、噪声；

R_sheet＝(|λ₃|+|λ₂|)/|λ₁|

R_blob＝(2|λ₁|-|λ₂|-|λ₃|)/2|λ₁|

$R_{\mathrm{noise}}=\sqrt{{\mathrm{\λ}}_3^2+{\mathrm{\λ}}_2^2+{\mathrm{\λ}}_1^2}$

对位于管状曲面上的体素，它的特征值具有如下规律：|λ1|=0and|λ2|=|λ3|>>0，因此，在此定义一种特征值的比率来描述这类几何结构：

R_tube＝(|λ₃|+|λ₂|-|λ₁|)/|λ₁|

当它们所代表的结构出现时，上述4种特征值的比率值会相应的变为0，据此，通过将特征向量、特征值和上述描述局部几何结构的比率值进行组合，来定义一个高维的向量场，其向量化表示形式可被定义为：

$\mathrm{ELSVR}=(e_1^T,e_2^T,e_3^T,{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,{\mathrm{\λ}}_3,R_{\mathrm{sheet}},R_{\mathrm{blob}},R_{\mathrm{tube}},R_{\mathrm{noise}}),$

借助三维图像的全局拉普拉斯矩阵，ELSVR将被用于各向异性热传递模型的构建。

3.根据权利要求1所述的基于各向异性热核分析的三维图像多尺度特征提取方法，其特征是：所述的基于定义的向量场中，在三维图像上数学严格地定义了一种各向异性热传递描述模型，该模型实现了三维图像局部结构的向量化，引入了基于局部结构变化加权的全局拉普拉斯算子，在ELSVR所定义的流形上，热的各向异性传播由图的全局拉普拉斯算子来决定，若将三维图像定义为一个无向图G=(V,E)，其中V代表所有的体素，E代表所有的边且满足E∈V×V，并且当且仅当相应体素满足一环领域要求时，边E_ij∈E的定义才成立，这里将某个体素的邻域体素记为N_i，为了诱使热沿相同的材质传播，而抑制其跨材质传播，以ELSVR空间的距离度量作为边的权重来对图的全局拉普拉斯算子进行计算，因此，首先需要度量两个相邻体素V_i和V_j在其所属材质类型方面的相似度：

$S_{i,j}=\left\{\begin{array}{cc}0,& \mathrm{if}{R}_{\mathrm{noise}}^i=0\mathrm{or}{R}_{\mathrm{noise}}^j=0\\ e^{-d_r},& \mathrm{otherwise}\end{array}\right.$

这里

而R={R_sheet,R_blob,R_tube}，权重W(i,j)可被定义为：

$W_{i,j}=\left\{\begin{array}{cc}{e}^{\frac{-(1-S_{i,j})d_{i,j}{d}_{i,j}^{\mathrm{\α}}}{{\mathrm{\σ}}^2}},& V_j\∈N_i\\ 0,& \mathrm{otherwise}\end{array}\right.$

其中d_ij代表V_i和V_j间的欧氏距离d_x,d_y,d_z代表图像采集时的采样间距，

$d_{i,j}^a=(\arccos (e_1^i\·e_1^j)+\arccos (e_2^i\·e_2^j)+\arccos (e_3^i\·e_3^j)),$ 而σ则是一个控制参数，它正比于两个像素间结构距离的方差，因此，图的全局拉普拉斯矩阵可被定义为：

$L_{i,j}=\left\{\begin{array}{cc}{W}_{i,j},& V_j\∈N_i\\ -{\mathrm{\Σ}}_{j\∈N_i}{W}_{i,j},& i=j\\ 0,& \mathrm{otherwise}\end{array}\right.$

W_i,j为相邻体素V_i和V_j的权重三维图像的各向异性热核算子

可使用拉普拉斯矩阵的有限个特征向量来进行近似，基于ELSVR加权的图的拉普拉斯算子自然的将图像局部结构的差异性纳入其中，相应的热核算子更具材质和结构敏感性，使得热的传递能受到局部和全局几何结构的控制。

4.根据权利要求1所述的基于各向异性热核分析的三维图像多尺度特征提取方法，其特征是：所述的基于定义的各向异性热传递模型中，用黑塞矩阵的第一特征值来定义热量场，并通过将图像颜色转换为热核信号，来辅助三维图像多尺度特征的提取，并定义尺度空间热核信号的极值点为特征点。

5.根据权利要求1所述的基于各向异性热核分析的三维图像多尺度特征提取方法，其特征是：所述的通过引入特征保持的下采样方法中，设计了一种由粗到细的层次化的热核加速策略，首先将原始图像进行下采样并进行类似金子塔形式的层次化组织，先在数据规模较小的低分辨率图像上进行热核计算和特征点位置的粗略定位，然后在原始图像相应的局部区域内，进行精细的热核计算，进而最终找出特征点的精确位置。

Images