一种变换域HMT模型的脑磁共振体数据自适应增强方法与流程

本发明属于医学图像处理领域，涉及一种脑磁共振体数据自适应增强方法，特别涉及一种变换域hmt(hiddenmarkovtree)模型的脑磁共振体数据自适应增强方法。

背景技术：

磁共振成像(magneticresonanceimaging,mri)是一种十分重要的先进医学成像技术，是一种革命性的医学诊断工具。该成像技术在不会对被检测对象造成危害的情况下，可以对人进行全方位扫描，利用人体不同组织含氢核数量的不同对应原子核运动产生的能级和相位变化也不同的特点，重构含有人体组织病理诊断信息的高精度断层立体图像。但是，在成像过程中，许多因素会导致图像质量下降，如噪声干扰、低对比度、成像伪影和几何变形等。这给基于磁共振图像的临床诊断和科学研究带来了挑战。为了克服磁共振图像的噪声，各种图像增强方法便应运而生。然而，图像增强存在一对矛盾，即增强边缘的同时也会增强噪声，而滤除噪声的同时又会使图像边缘在一定程度上模糊化。因此，进一步探究更好的图像增强方法，从而有效抑制磁共振图像的噪声是一项非常有意义和价值的研究课题。

在众多图像增强方法中，基于小波分析的方法取得了巨大成功。目前的小波方法有很多子方法且各有特点。最早由weaver等人提出的小波阈值噪声抑制方法是其中的子方法之一，该方法需要选定一个阈值，若小波系数比阈值大则置零。紧接着，wood和johnson对磁共振图像原始数据的实部和虚部分别进行基于小波包阈值噪声抑制方法，提高了方法的性能。随后，gregg等人提出一种幅值平方的小波阈值噪声抑制方法，该方法的性能优于上述小波阈值方法。同时，alexander等人提出既对幅值又对相位进行噪声抑制的小波方法，此方法有效提高了磁共振图像的信噪比和对比度。然后，yang等人提出小波多尺度噪声抑制方法，该方法考虑了小波系数随尺度的变化，方法的性能变得更好。然而，基于小波的阈值噪声抑制方法存在一些问题和不足：对阈值的选择比较敏感，阈值过大则不能有效抑制图像噪声，过小则容易丢失图像中的有用信息。

统计模型小波噪声抑制方法则可以一定程度上克服小波阈值方法的不足。由于小波系数具有服从非高斯分布的特性。同时，小波系数还具有其他两个特性，即聚集性和持续性：聚集性表现为如果一个小波系数是大或小的值，那么其相邻小波系数也极有可能为大或为小的值；持续性表现为在不同尺度上小波系数的大、小值沿着尺度传递。根据小波系数的这些特性，bouman等人在1991年提出建立独立的非高斯小波系数统计模型，该模型忽略了小波系数间的相关性，而lee等人提出建立联合高斯模型，则没有考虑小波系数分布的非高斯性。为了同时考虑小波系数的这3个特性，course等人提出小波域的隐马尔科夫模型，来对图像进行增强。

本发明为了克服上述方法中存在的问题，弥补现有小波阈值噪声抑制方法不足，通过充分利用隐马尔科夫模型能够对小波系数特性有效建模的优势，提出了基于小波域的隐马尔科夫树模型对脑磁共振体数据进行自适应增强的方法。

技术实现要素：

本发明的目的在于提出一种变换域hmt模型的脑磁共振体数据自适应增强方法，该方法将小波变换和隐马尔科夫链结合起来。具体来说，根据单个小波系数的概率密度函数呈高峰值、长拖尾的非高斯分布的特性，对单个小波系数的随机性建立高斯混合模型。同时，小波系数在尺度间传递的持续性采用隐马尔科夫树来描述。通过上述小波域隐马尔科夫树模型，实现脑磁共振体数据的自适应增强，以期进一步有效抑制数据的噪声。

本发明采用的技术方案由八个步骤组成，即如图1所示的流程图。各个步骤具体描述如下。

第一步，三维小波变换。

小波系数具有聚集性、持续性、局部性、多分辨性、边缘检测性、去相关性、以及能量紧束性等良好特性，是一种极其重要的时频分析工具。因此，选择三维小波变换对体数据进行增强。

第二步，建立双状态隐马尔科夫树模型。

单个小波系数的概率密度函数呈高峰值、长拖尾的非高斯分布，故使用两个高斯分布拟合小波系数的非高斯分布。其中小方差、零均值的高斯分布对应大部分幅值小的小波系数(小状态)，另外大方差、零均值的高斯分布对应少数幅值大的小波系数(大状态)。进一步结合小波系数的其他两个特性，即簇性聚集性和持续性，建立小波系数的双状态隐马尔科夫树模型。

定义描述单个小波系数非高斯分布的概率密度函数为

式中m为状态数(m＝2)；s为状态变量，即s＝1代表小状态，s＝2代表大状态；ps(m)表示小波系数在m状态的概率；f(w|s＝m)表示在状态m下小波系数对应的高斯概率密度函数。

确定双状态隐马尔科夫树模型的参数集θ＝(π,a,b)。π为概率分布函数，即根节点的概率分布函数其中m＝1,2,…,m；a为状态转移概率矩阵。所有父子节点的状态转移概率组成状态转移概率矩阵，即表示父节点小波系数状态sρ(i)＝n时，子节点系数状态si＝m的概率，其中m,n＝1,2,…,m；b为观测值概率分布矩阵。矩阵中的元素为小波系数对应不同状态时的均值μi,m和方差其中m＝1,2,…,m。故双状态隐马尔科夫树模型的参数集可具体表示为

式中s为隐状态随机变量，w为小波系数随机变量，s1∈s为小波树的根节点的隐状态，w1∈w为小波树的根节点的小波系数，si∈s为小波树第i个节点处小波系数的隐状态，wi∈w为小波树第i个节点处小波系数，ρi为节点i的父节点。

第三步，捆绑具有相同统计特性的小波系数树。

对三维数据进行小波变换，在七个高频方向得到各自的模型参数集，记为θ＝{θ^hhh,θ^hhl,θ^hlh,θ^lhh,θ^hll,θ^lhl,θ^llh}。由于每棵小波树上每个节点对应的模型参数不同，这样会导致模型参数量过大而难于估计。故假设在同一尺度上，不同位置处的小波系数具有相同的统计特性，则位于同一尺度上的节点具有相同的模型参数集，对这些具有相同统计特性的小波系数树进行捆绑，从而降低增强方法的复杂度。

第四步，初始化模型参数集。

初始化状态数m＝2，循环计数l＝0。初始化模型参数集θ⁰为：初始化概率分布初始化状态转移概率初始化均值μi,m＝0，初始化方差其中表示节点i所在尺度j上小波系数的方差。

第五步，确定小波系数的条件概率。

为了在模型参数θ^l下，计算条件概率p(si＝m|w,θ^l)和定义如下变量。以节点i为根，子树ti包括小波系数wi及其所有后代。tp是ti的一棵子树，定义ti\p表示在ti中除去子树tp后其余的所有小波系数。特别的，t1表示整棵小波树，小波系数w1位于整棵小波树的根部。

定义条件似然函数为bi(m)＝f(ti|si＝m,θ)，定义联合概率函数为ai(m)＝p(si＝m,t1\i|θ)。根据贝叶斯公式得到小波系数的条件概率为

第六步，模型的em算法求解。

e步计算

(1)定义高斯概率密度函数表示为

(2)为了计算条件概率p(si＝m|w,θ^l)和首先需要求得条件似然函数值bi(m)、bi,ρ(i)(m)、bρ(i)\i(m)和联合概率函数值ai(m)。而bi(m)、bi,ρ(i)(m)和bρ(i)\i(m)值是通过传播粗尺度小波系数的信息来影响细尺度小波系数的状态来求得，这称之为向上算法；ai(m)值是通过转移细尺度小波系数的信息来影响粗尺度小波系数的状态来求得，这称之为向下算法。

对单棵小波树进行e步计算—向上算法

①初始化j＝1(最细尺度)上节点i处的似然函数为

②针对尺度j上的节点，计算

③令j＝j+1(由细尺度向粗尺度迭代，即向树的上一层迭代)；

④当j＝j(最粗尺度)则停止，否则回到第②步循环执行。

对单棵小波树进行e步计算—向下算法

①初始化最粗尺度(j＝j)上根节点处的联合概率密度函数为

②令j＝j-1(由粗尺度向细尺度迭代，即向树的下一层迭代)；

③针对尺度j上的节点，计算

④如果j＝1则停止，否则回到第②步循环执行。

(3)对k棵捆绑的小波树进行e步计算。假设第k棵小波树有p个节点，其小波系数记为对应隐状态记为则所有小波系数记为w＝(w¹,w²,…,w^k)，对应隐状态为s＝(s¹,s²,…,s^k)，使用向上算法和向下算法分别对每一棵小波树计算条件似然函数值bi(m)、bi,ρ(i)(m)、(ρ(i)\i(m)和联合概率函数值ai(m)，则可以得到第k棵小波树对应的条件概率和

m步计算

(4)更新模型参数θ^l+1为

(5)判断迭代是否收敛。若‖θ^l+1-θ^l‖足够小，循环停止；否则，l＝l+1从“e步计算”继续进行模型更新。

第七步，估计经噪声抑制后的小波系数。

首先，根据求解的hmt模型，使用最小均方误差准则，第k棵小波树的第i个节点的小波系数yi^k在隐状态si下的条件期望为

式中am是增益因子，对于噪声系数，am＝0；对于边缘系数，am＞1；是n层分解噪声小波系数方差的估计。对体数据进行第一次小波分解后得到高频子带hhh1，它主要由噪声小波系数组成，而该子带中也含有极少数小幅值的真实信号小波系数。为了消除真实信号对小波系数噪声方差的估计，需要对高频子带hhh1小波系数进行中值滤波。于是使用donoho等人提出的小波域噪声方差估计策略，得到各层分解的噪声小波系数方差的估计第一层分解的噪声小波系数方差估计为

由于噪声小波系数在尺度间具有持续性，其他各层的噪声小波系数方差估计为

然后，使用通过em算法求解得到的第k棵小波系数树对应条件概率计算在双状态隐马尔科夫树模型下小波系数yi^k的条件期望为

上式为经噪声抑制后的小波系数。

第八步，三维小波逆变换。

使用第k棵小波系数树节点i处小波系数的条件期望替代对这些经过噪声抑制后的小波系数进行三维小波逆变换，得到增强的脑磁共振体数据。

为了评价本发明所提出技术方案的有效性，采取主观和客观相结合的评价方式。

(一)主观评价方式

所谓主观评价方式，指的是通过肉眼直接对增强的图像进行观察，同时评判图像增强的效果。

(二)客观评价方式

在一般情况下，无法采集原始无噪声图像，传统以原始无噪声图像为参考的图像质量评价指标则不能使用。为克服这个问题，本发明采用视觉信息保真度(visualinformationfidelity,vif)来客观地评价增强图像的质量。视觉信息保真度的定义具体描述如下。

在视觉信息保真度的评价中，图像增强的图像源作为ir参考图像，增强图像it作为测试图像。首先，分别对参考图像及增强图像进行j层小波分解，并将小波子带分成3×3不重叠的系数块，从小波系数块的相同位置中提取小波系数分别构成向量cri和dti，这里i＝1,…,n,n＝9×(7j+1)。

假设图像源ir符合高斯尺度混合(gaussianscalemixtures,gsms)模型，则从参考图像的小波子带不重叠块中提取的小波系数cri为该模型中的随机向量

式中uri是一个0均值的高斯向量，是一个正标量。来自相同层小波分解子带的系数块向量cri之间具有相同的协方差结构cu。将图像增强运算作为图像信号失真通道。增强图像it可视为图像源ir在小波域的信号cri在增益gri作用下输出与附加噪声v叠加的信号

dti＝gricri+v

式中gri是一个确定的标量增益场，而v是一个独立平稳的加性0均值高斯白噪声场，其协方差为i是单位矩阵。

其次，人眼视觉失真模型在小波变换域中建模为平稳0均值信号附加白高斯噪声

ei＝cri+n

fi＝dti+n′

其中ei和fi分别为人眼视觉感知参考图像和增强图像的相同小波子带系数块的系数向量，同时n和n′为0均值高斯噪声随机场，他们的协方差矩阵均为i是单位矩阵。

再者，令表示随机向量场cr的n个实现集合。同理，对于dti，ei和fi，类似地定义d，e和f，且令同时，假设模型参数gri,i＝1,…,n，和均已知，各小波子带互相独立的，于是可以将特定子带互信息扩展到k＝(7j+1)个子带求和。对于参考图像ir，给定z，则各子带各系数块c和e之间的总互信息为

这里，h(c)为连续随机向量的熵微分。同理，对于增强图像，有

。对于第j层小波分解子带的系数块，由于其协方差矩阵是对称的，它可写成对角矩阵，其对角线上的一组特征值分别为于是以上两式可写为

i(c；e|z)和i(c；f|z)分别为人类视觉系统(humanvisualsystem,hvs)从参考图像和增强图像k个子带中提取的信息，则增强图像质量的客观评价指标vif定义为

为了计算vif的值，还需要对模型参数zri、gri、j＝1,…,j进行估计。

附图说明

图1是小波自适应增强方法的流程图。

图2是含有噪声的原始脑磁共振图像以及分别进行二、三、四、五层分解的小波自适应增强结果。图2(a)为原始脑磁共振噪声图像；图2(b)为二层分解的小波自适应增强结果；图2(c)为三层分解的小波自适应增强结果；图2(d)为四层分解的小波自适应增强结果；图2(e)为五层分解的小波自适应增强结果。

图3是不同分解层数小波自适应增强方法的vif评价。

图4是含有噪声的原始脑磁共振图像、小波阈值增强结果、以及小波自适应增强结果。图4(a)为原始脑磁共振噪声图像；图4(b)为小波阈值增强结果；图4(c)为小波自适应增强结果。

图5是3层分解小波自适应增强方法和小波阈值增强方法的vif评价。

具体实施方式

本发明提出了一种变换域hmt模型的脑磁共振体数据自适应增强方法，该方法的技术方案流程图如图1所示。具体的实施方式说明如下。

第一步，采集脑磁共振体数据。采集的3t脑磁共振图像分辨率为0.6×0.6×0.6mm³。由于全脑磁共振图像比较大，不易于分析增强方法的性能，截取全脑磁共振图像的一部分，其大小为64×64×64。

第二步，三维离散小波变换。对于脑磁共振体数据，三维离散小波变换通过分别沿着三个坐标轴方向的离散小波变换来实现。离散小波变换分解的层数为2至5层。

第三步，建立双状态隐马尔科夫树模型。由于单个小波系数的概率密度函数呈高峰值、长拖尾的非高斯分布。本发明采用两个高斯分布来拟合小波系数的非高斯分布：一个是小方差、零均值的高斯分布对应大部分幅值小的小波系数，即小状态；另一个是大方差、零均值的高斯分布对应少数幅值大的小波系数，即大状态。

第四步，求解隐马尔科夫树模型。小波域隐马尔科夫树模型的求解使用em算法。首先初始化模型参数集θ⁰，然后通过em算法对模型参数进行更新，并根据最终的模型参数θ计算小波系数所属大小状态的条件概率。

第五步，估计经噪声抑制后的小波系数。根据小波系数所属大小状态的条件概率，估计在双状态隐马尔科夫树模型下小波系数yi^k的条件期望

第六步，离散小波逆变换。在对所有高频子带小波系数进行了经噪声抑制的估计后，使用离散小波逆变换得到增强的脑磁共振体数据。

本发明提出了一种变换域hmt模型的脑磁共振体数据自适应增强方法。通过主观和客观相结合的评价方式可知，提出方法在三层小波分解情况下的脑磁共振图像增强效果最好(图2和3)。另外，提出方法的图像增强效果比小波阈值增强方法更好(图4)。同时，增强图像的vif值也明显要大(图5)，这进一步表明提出方法得到的增强图像有更好的效果，即具有更优的视觉信息保真度。