CN106296633A - 一种基于多尺度图像域双能分解算法 - Google Patents

Abstract

本发明公开了一种基于多尺度图像域双能分解的算法。采用滤波反投影(FBP)算法根据实测高能与低能投影分别重建高、低能CT图像；针对高、低能CT图像，分别进行相同参数的图像域多尺度分解和方差图估计，得到各尺度下高、低能CT图像及方差图；针对多尺度分解中相同尺度高、低能CT图像及方差图，进一步做双能分解，得到对应尺度下各基材料的CT图像；把所有尺度下双能分解得到的各基材料CT图像分别进行累加，得到各基材料的CT图像。本发明的多尺度图像域双能分解的方法得到骨和软组织基材料图像与单尺度相比，在相同的噪声水平下，保持更好的对比度及空间分辨率，并且尺度因子越大，基材料图像的对比度和空间分辨率越高。

Description

一种基于多尺度图像域双能分解算法

技术领域

本发明涉及医学工程技术领域，具体涉及一种多尺度图像域双能分解算法。

背景技术

目前，现代X射线CT成像广泛应用于医学临床诊断与治疗应用中，具有巨大的社会价值与意义。近些年，随着各种CT成像新模态的不断涌现，双能CT及其图像域分解表现出巨大的临床潜力。双能CT及其图像域分解首先可采用任意高、低能切换技术，达到实测双能投影的目的；然后采用滤波反投影(FBP)算法根据实测高能与低能投影分别重建高、低能CT图像；最后根据针对高、低能CT图像，进行图像域分解。图像域分解存在多种方法，包括直接求逆法、极小化目标函数(如单尺度PWLS)法等。

直接求逆法，优点是原理简单、实现简便，缺点是双能分解得到的基材料图像噪声非常大，对临床诊断与治疗质量有着明显的负面影响。

单尺度PWLS法，在保持直接求逆法优点的同时，可以有力地克服直接求逆法基材料图像信噪比差的缺陷，极大地提升分解图像的质量。但是，需要注意的是单尺度PWLS法在信噪比与空间解析度之间的平衡方面，仍然未取得最优分解性能。

该缺点的可能原因是，CT图像边缘与噪声在尺度空间中的分布特点具有显著差异——非相关噪声几乎均匀地分布于整个尺度空间，而台阶式图像边缘的傅里叶变换随着频率增大近似呈现出指数衰减特点。单尺度PWLS法，对信号与噪声的尺度不加区分，因此可能造成二者互相干扰。

发明内容

针对现有技术不足，本发明提出了一种多尺度图像域双能分解算法，其中包括：

1.高低能投影数据重建CT图像

本发明将实测投影包括高、低能投影数据采用滤波反投影算法重建出高、低能CT图像。

2.图像域多尺度分解

本发明采用如下各向同性热扩散方程进行高低能CT图像在图像域的多尺度分解：

$\frac{\∂}{\∂t}f(x,y,t)=\mathrm{d\Δf}(x,y,t),t\≥0,---\left(1\right)$

f_s(x，y)＝f(x，y，t_s)-f(x，y，t_s-1)，t_s＞t_s-1，t₀＝0， (2)

其中(x，y)为二维CT图像中任意一点；f(x，y)为二维高或低能CT图像；t为时间；d为热扩散常数；Δ为拉普拉斯算子；f_s(x，y)为多尺度分解得到的尺度s下的CT图像，t_s为对应尺度下的热扩散时间。

在多尺度图像分解中使用高斯卷积核函数，在t＝0时，定义f(x，y，0)≡f(x，y)，则：

$f(x,y,t)=g(x,y,{\mathrm{\σ}}^2)\⊗f(x,y,0),{\mathrm{\σ}}^2=2\mathrm{dt},---\left(3\right)$

符号表示关于(x，y)的二维卷积操作，g(x，y，σ²)是高斯卷积核函数，定义为：

$g(x,y,{\mathrm{\σ}}^2)=\frac{1}{2\mathrm{\π}{\mathrm{\σ}}^2}\exp (-{\left|\right|(x,y)\left|\right|}_2^2/\left(2{\mathrm{\σ}}^2\right)),---\left(4\right)$

符号||·||₂表示L₂正则；对高斯卷积核函数做二维傅里叶变换可得：

$G(v,w,{\mathrm{\σ}}^2)=\exp (-{\left|\right|2\mathrm{\π}(v,w)\left|\right|}_2^2{\mathrm{\σ}}^2/2),---\left(5\right)$

(v，w)为空域变量(x，y)的频域对偶变量。

3.方差图估计

本发明为了估计每个尺度下的方差图，首先在高或低能CT图像中选择一个均匀区域并估计数值方差常数C，然后高或低能CT图像方差图可按如下来假设：

V(x，y)≡C， (6)

最后我们仿照图像域多尺度分解的方式，按如下公式估计每个尺度下的方差图：

V_s(x，y)＝V(x，y，t_s)+V(x，y，t_s-1)，V(x，y，0)＝V(x，y)， (7)

其中忽略协方差影响，V(x，y，t_s)根据公式(1)的离散形式和公式(6)的结果可很容易地估计出来。值得注意的是，公式(7)使用了两个随机变量(RVs)和的方差公式来设计一个算法，自动地估计每个尺度的方差图，这比手动估计方差更有效。

4.基于PWLS的双能分解

本发明采用惩罚加权最小均方(PWLS)算法对高、低能图像双能分解：

(3-1)图像域多尺度分解的尺度总数S为大于1的整数如4，即s＝1，2，…，S；图像域双能分解的基材料总数为2，分别对应如骨和软组织；

(3-2)采用惩罚加权最小均方(PWLS)算法对高、低能图像双能分解，可实现基材料边分解边降噪，其中极小化目标函数表示如下：

$\mathrm{\Φ}\left({\stackrel{\→}{\mathrm{\μ}}}_s\right)={(A{\stackrel{\→}{\mathrm{\μ}}}_s-{\stackrel{\→}{f}}_s)}^{\′}{V}_s^{-1}(A{\stackrel{\→}{\mathrm{\μ}}}_s-{\stackrel{\→}{f}}_s)+{\stackrel{\→}{\mathrm{\β}}}_s^{0}R\left({\stackrel{\→}{\mathrm{\μ}}}_s\right),---\left(8\right)$

其中A为材料成分矩阵，不依赖于尺度s，维数为2N×2N，N为一个二维CT图像的像素总数； ${\stackrel{\→}{\mathrm{\μ}}}_s={\left[\begin{array}{cc}{\stackrel{\→}{\mathrm{\μ}}}_{1,s}^{\text{'}}& {\stackrel{\→}{\mathrm{\μ}}}_{2,s}^{\text{'}}\end{array}\right]}^{\text{'}}$ 为2N×1维列向量，与分别是将两种基材料图像拉直后的列向量； ${\stackrel{\→}{f}}_s={\left[\begin{array}{cc}{\stackrel{\→}{f}}_{H,s}^{\text{'}}& {\stackrel{\→}{f}}_{L,s}^{\text{'}}\end{array}\right]}^{\text{'}}$ 为2N×1维列向量，与分别是多尺度分解得到的高、低能CT图像拉 直后的列向量；V_s为对角矩阵，维数为2N×2N，对角元素分别是高、低能CT图像中每像素的噪声方差；为正则项，决定双能分解图像的变动强度；为尺度s下的初设正则系数，决定极小化目标函数中的保真项与正则项的平衡关系；

(3-3)采用极小化目标函数对高、低能图像进行双能分解，对初设正则系数可通过如下方法调整：

${\stackrel{\→}{\mathrm{\β}}}_s(x,y)={\stackrel{\→}{\mathrm{\β}}}_s^0/(1+k_s|\▿{\stackrel{\‾}{\mathrm{\μ}}}_s(x,y)|/\max |\▿{\stackrel{\‾}{\mathrm{\μ}}}_s\left|\right),---\left(9\right)$

其中k_s为伸缩因子，影响双能分解得到图像的空间分辨率；为采用直接求逆得到的双能分解基材料图像；为梯度算子；||表示求梯度矢量幅度。

所述步骤(3-2)具体如下：

(3-2-1)其中是基材料i＝1，2在能量j＝H，L下的等效线性衰减系数，I是N×N维单位阵；

(3-2-2)V_s为对角矩阵，维数为2N×2N，对角元素分别是高、低能CT图像在尺度s下每像素的噪声方差。

$V_s=D\left(\begin{array}{cc}{\stackrel{\→}{V}}_{H,s}^{\′}& {\stackrel{\→}{V}}_{L,s}^{\′}\end{array}\right),---\left(10\right)$

其中D表示对行向量进行对角线化操作，和为高低能CT图像在尺度s下估计的表示为列向量的方差。

5.各尺度的分解材料对应相累加

本发明把所有尺度(s＝1，2，...，S)下双能分解得到的骨图像和软组织图像进行累加，即得骨分解图像和软组织分解图像

与现有技术相比，本发明的基于多尺度图像域双能分解的方法得到骨和软组织基材料图像与单尺度相比，在相同的噪声水平下，保持更好的对比度及空间分辨率。

附图说明

图1为本实施例的多尺度图像域双能分解方法的流程图；

图2为原始高能(125kV)与低能(75kV)二维CT图像，显示窗为[0.015 0.045]mm^-1；

图3为高、低能二维CT图像直接矩阵求逆双能分解的骨和软组织基材料图像，显示窗[0.2 1.2]；

图4为高、低能二维CT图像单尺度PWLS双能分解的骨和软组织基材料图像，显示窗[0.2 1.2]；

图5为高、低能二维CT图像多尺度PWLS(k＝20)双能分解的骨和软组织基材料图像，显示窗[0.2 1.2]；

图6为高、低能二维CT图像多尺度PWLS(k＝80)双能分解的骨和软组织基材料图像，显示窗[0.2 1.2]；

图7为在不同方法分解的骨图像里分别从最大的线对开始在半个圆周内采样1800个点的一维曲线比较图；

图8为在不同方法分解的软组织图像里分别从最大的线对开始在半个圆周内采样1800个点的一维曲线比较图。

具体实施方式

下面将结合附图和具体实施例对本发明进行详细描述。

1.在一个实施例中，将实测投影包括高、低能投影数据采用滤波反投影算法重建出高、低能CT图像，验证本发明相关算法对双能分解基材料图像噪声、对比度、空间解析度的有效程度；

2.上述技术方案中，所述的采用滤波反投影算法重建出高、低能CT图像如图(2)，骨和软组织基材料图像在同幅图中；

3.上述技术方案中，所述的高、低能CT图像直接矩阵求逆分解的骨和软组织基材料图像如图(3)，可发现分解的基材料图像中含有很大的噪声；

4.上述技术方案中，所述的高、低能CT图像单尺度PWLS基材料分解图像如图(4)，骨和软组织的正则系数分别为1.9×10^-5，7.6×10^-5，可发现分解的基材料图像中噪声有很大程度的降低；

5.上述技术方案中，所述的高、低能CT图像多尺度PWLS基材料分解图像如图(5)，其中各尺度下尺度因子相同为k＝20并且最大尺度下正则系数初始值与图4一样，其它尺度下的骨和软组织的正则系数初始值分别为2×10^-5，8×10^-5，为了在分解的基材料图像中有相同的噪声水平(选一个均匀区域如图3中矩形框所示，计算其方差值作为噪声水平)，其它尺度(除最大尺度)下的骨和软组织的正则系数初始值会有轻微的调整；

6.上述技术方案中，所述的高、低能CT图像多尺度PWLS基材料分解图像如图(6)，其中各尺度下尺度因子相同为k＝80并且最大尺度下正则系数初始值与图4一样，其它尺度下的骨和软组织的正则系数初始值分别为2.5×10^-5，1.5×10^-4，为了在分解的基材料图像中有相同的噪声水平(选一个均匀区域如图3中矩形框所示，计算其方差值作为噪声水 平)，其它尺度(除最大尺度)下的骨和软组织的正则系数初始值会有轻微的调整；

7.上述技术方案中，所述的不同方法双能分解基材料图像一维曲线对比图如图(7)、(8)，可发现，所提多尺度图像域双能分解方法对分解的基材料图像降噪的同时可保持很好的对比度和空间分辨率，并且尺度因子k越大，基材料图像的对比度和空间分辨率越高。

8.以上所述的具体实施方式对本发明的技术方案和有益效果进行了详细说明，应理解的是以上所述仅为本发明的最优选实施例，并不用于限制本发明，凡在本发明的原则范围内所做的任何修改、补充和等同替换等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种多尺度图像域双能分解算法，其特征在于，包括：

(1)采用滤波反投影(FBP)算法根据实测高能与低能投影分别重建高、低能CT图像；

(2)针对高、低能CT图像，分别进行相同参数的图像域多尺度分解及方差图估计，得到各尺度下高、低能CT图像及方差图；

(3)针对多尺度分解中相同尺度高、低能CT图像及方差图，进一步做双能分解，得到对应尺度下各基材料的CT图像；

(4)把所有尺度下双能分解得到的各基材料CT图像分别进行累加，得到各基材料的CT图像。

2.如权利要求1所述的多尺度图像域双能分解算法，其特征在于，所述步骤(1)可采用任意高、低能切换技术，达到实测双能投影的目的。

3.如权利要求1所述的多尺度图像域双能分解算法，其特征在于，所述步骤(2)图像域多尺度分解采用如下各向同性热扩散方程进行：

f_s(x，y)＝f(x，y，t_s)-f(x，y，t_s-1)，t_s＞t_s-1t₀＝0， (2)

其中(x，y)为二维CT图像中任意一点；f(x，y)为二维高或低能CT图像；t为时间；d为热扩散常数；Δ为拉普拉斯算子；f_s(x，y)为多尺度分解得到的尺度s下的CT图像，t_s为对应尺度下的热扩散时间。

4.如权利要求1所述的多尺度图像域双能分解算法，其特征在于，所述步骤(2)方差图估计首先分别在高或低能CT图像中选择一个均匀区域估计数值方差常数C，然后高或低能CT图像方差图可按如下来假设：

V(x，y)≡C， (3)

最后仿照图像域多尺度分解，按下式估计每个尺度下的方差图：

V_s(x，y)＝V(x，y，t_s)+V(x，y，t_s-1)，V(x，y，0)＝V(x，y)， (4)

其中忽略协方差影响，V(x，y，t_s)根据式(1)的离散形式和式(3)的结果可很容易地估计出来。

5.如权利要求1所述的多尺度图像域双能分解算法，其特征在于，所述步骤(3)通过如下方法进行多尺度高、低能图像双能分解：

(3-1)图像域多尺度分解的尺度总数S为大于1的整数如4，即s＝1，2，...，S；图像域双能分解的基材料总数为2，分别对应如骨和软组织；

(3-2)采用惩罚加权最小均方(PWLS)算法对高、低能图像双能分解，可实现基材料边分解边降噪，其中极小化目标函数如下表示：

其中A为材料成分矩阵，不依赖于尺度s，维数为2N×2N，N为一个二维CT图像的像素总数；为2N×1维列向量，与分别是将两种基材料图像拉直后的列向量；为2N×1维列向量，与分别是多尺度分解得到的高、低能CT图像拉直后的列向量；V_s为对角矩阵，维数为2N×2N，对角元素分别是高、低能CT图像中每像素的噪声方差；为正则项，决定双能分解图像的变动强度；为尺度s下的初设正则系数，决定极小化目标函数中的保真项与正则项的平衡关系；

(3-3)采用极小化目标函数对高、低能图像进行双能分解，对初设正则系数可通过如下方法调整：

其中k_s为伸缩因子，影响双能分解得到图像的空间分辨率；为采用直接求逆得到的双能分解基材料图像；为梯度算子；||表示求梯度矢量幅度。

6.如权利要求4所述的多尺度图像域双能分解算法，其特征在于，所述步骤(3-2)具体如下：

(3-2-1)其中是基材料i＝1，2在能量j＝H，L下的等效线性衰减系数，I是N×N维单位阵；

(3-2-2)V_s为对角矩阵，维数为2N×2N，对角元素分别是高、低能CT图像在尺度s下每像素的噪声方差。

其中D表示对行向量进行对角线化操作，和为高低能CT图像在尺度s下估计的表示为列向量的方差。

7.如权利要求1所述的多尺度图像域双能分解算法，其特征在于，所述步骤(4)把所有尺度(s＝1，2，...，S)下双能分解得到的骨图像和软组织图像进行累加，即得骨分解图像和软组织分解图像