CN105512487A - 图像域双能ct多材料统计分解方法 - Google Patents

Abstract

本发明公开了一种图像域双能CT多材料统计分解方法。采用滤波反投影(FBP)算法根据模拟和实测高能与低能投影分别重建高、低能CT图像；针对高、低能CT图像，进行图像域直接多材料分解，得到直接分解的基材料图像；针对高、低能CT图像，分别选一个均匀区域估计其方差，得到高、低能CT图像的数值方差；针对直接分解的基材料图像和高、低能CT图像的数值方差，采用基于逐像素点可分离的目标函数进行双能多材料统计分解，得到基于统计方式分解的基材料图像。本发明的图像域双能CT多材料统计分解方法得到的基材料图像与直接分解的基材料图像相比，可以大幅度地降低基材料图像的噪声水平，并且保持很好的空间分辨率，且对分解精确度有极小的影响。

Description

图像域双能CT多材料统计分解方法

技术领域

本发明涉及医学工程技术领域，具体涉及一种图像域双能CT多材料统计分解方法。

背景技术

目前，现代X射线CT成像广泛应用于医学临床诊断与治疗应用中，具有巨大的社会价值与意义。近些年，随着各种CT成像新模态的不断涌现，双能CT及其材料分解的特性表现出巨大的临床潜力。理论上，两种基材料能从双能CT测量中精确地分解出来，然而很多临床应用中需要知道三种或更多的材料成分。从双能CT测量中分解出L₀≥3种基材料需要外加的约束和假设，面临很大的挑战。双能CT多材料分解方法包括直接图像域分解法，极小化目标函数(如统计方法)法等。

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

极小化目标函数分解法可在投影域中进行也可在图像域中进行，文献中提出的投影域双能CT多材料分解方法分解精度高，但实施复杂，计算量大。文献中提出的图像域分解的方法目前分解的基材料图像最多为3种，极大地限制了多材料分解的应用。

发明内容

本发明的目的在于针对现有技术的不足，提供一种图像域双能CT多材料统计分解方法。

本发明的目的是通过以下技术方案来实现的：一种图像域双能CT多材料统计分解方法，该方法包括以下步骤：

(1)采用滤波反投影(FBP)算法根据模拟和实测高能与低能投影分别重建高、低能CT图像；所述高能为120-150kVp，低能为60-80kVp；

(2)针对高、低能CT图像，进行图像域直接多材料分解，得到直接分解的基材料图像；

(3)针对高、低能CT图像，分别选一个均匀区域估计其方差，得到高、低能CT图像的数值方差；

(4)针对直接分解的基材料图像和高、低能CT图像的数值方差，采用基于逐像素点可分离的目标函数进行双能多材料统计分解，得到基于统计方式分解的基材料图像。

进一步地，所述步骤(1)可采用任意高、低能切换技术，达到实测双能投影的目的。

进一步地，所述步骤(2)图像域直接多材料分解的方程如下所示：

${\stackrel{\→}{\mathrm{\μ}}}_E=\underset{l=1}{\overset{L_0}{\Σ}}{\mathrm{\μ}}_{lE}{\stackrel{\→}{x}}_l,---\left(1\right)$

其中为l材料的体积分数图像；L₀为分解的基材料种类数目；为E能量水平下的CT图像；μ_lE为E能量水平下的l材料的有效线性衰减系数；并且应该满足下面的约束和假设：

$\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1,\∀p,---\left(2\right)$

$0\≤x_{lp}\≤1,\∀l,p---\left(3\right)$

$\underset{l=1}{\overset{L_0}{\Σ}}{1}_{\{x_{lp}\≠0\}}\≤3,\∀p.---\left(4\right)$

其中p为体积分数图像的第p像素点；公式(2)为体积守恒约束；公式(3)为箱约束；公式(4)为每像素点p最多包含3种基材料的假设。

进一步地，所述步骤(3)数值方差估计为分别在高、低能CT图像中选择一个均匀区域估计其数值方差作为CT图像每像素的噪声方差。

进一步地，所述步骤(4)采用基于逐像素点可分离的目标函数进行双能多材料统计分解通过如下方法进行：

(4-1)图像域双能CT多材料统计分解的目标函数如下所示：

$\mathrm{\Φ}\left(\stackrel{\→}{x}\right)={(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})}^T{V}^{-1}(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})+R\left(\stackrel{\→}{x}\right),---\left(5\right)$

其中A为材料成分矩阵，表示为维数为2N_p×L₀N_p，N_p为1个CT图像的像素个数，I是N_p×N_p维单位矩阵，H，L表示为高、低能谱；为2N_p×1维列向量，与分别是高、低能CT图像拉直后的列向量；T表示转置操作；为L₀N_p×1维列向量，到均为基材料拉直后的列向量；维数为2N_p×2N_p，对角元素分别是高、低能CT图像每像素的噪声方差，其中D表示对行向量进行对角线化操作，和为高、低能CT图像估计的表示为列向量的噪声方差；为正则项，公式表示如下：

$R\left(\stackrel{\→}{x}\right)=\underset{l=1}{\overset{k}{\Σ}}{\mathrm{\β}}_l{R}_l\left({\stackrel{\→}{x}}_l\right),---\left(6\right)$

$R_l\left({\stackrel{\→}{x}}_l\right)=\underset{p=1}{\overset{N_p}{\Σ}}\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ψ}}_l(x_{lp}-x_{lk})---\left(7\right)$

${\mathrm{\ψ}}_l\left(t\right)=\frac{{\mathrm{\δ}}_l^2}{3}(\sqrt{1+3{\left(\frac{t}{{\mathrm{\δ}}_l}\right)}^2}-1)---\left(8\right)$

其中x_lk是x_lp像素点的邻域，N_lp是x_lp像素点邻域的集合，k_lp和k_lk是激励均匀空间分辨率的参数；正则系数β_l和边界保持参数δ_l可以根据不同的材料来选择以实现不同的平滑正则水平；

(4-2)在公式(2)-(4)的约束下最小化目标函数可以获得分解的基材料图像公式如下：

$\stackrel{\→}{x}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}subjectto\left(2\right)-\left(4\right)}{\arg \min }\mathrm{\Φ}\left(\stackrel{\→}{x}\right).---\left(9\right)$

但此目标函数很难直接最小化，本发明提出应用最优转换原理设计出一个目标函数的逐像素点可分离的二次替代函数在每次迭代过程中来单调地降低目标函数直至达到最大预设迭代次数；这个替代函数在第n次迭代过程中很容易最小化，得出的即为此替代函数的最小值，表示如下：

${\stackrel{\→}{x}}^{(n+1)}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}}{\arg }{\mathrm{min\φ}}^{\left(n\right)}\left(\stackrel{\→}{x}\right).---\left(10\right)$

分为保真项的二次替代函数和正则项的二次替代函数

(4-3)对p像素点图像域双能CT多材料统计分解方法的目标函数进行公式(2)-(4)约束下的求解。

进一步地，所述步骤(4-2)保真项的二次替代函数如下：

$\begin{array}{c}{L}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\\ +\frac{1}{2}{(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})}^T{D}_L^{\left(n\right)}(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\end{array}---\left(11\right)$

其中表示梯度操作；是关于p＝1,...,N_p的分块矩阵，公式如下：

$D_L^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{D_{L,p}^{\left(n\right)}\right\},---\left(12\right)$

$D_{L,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}2A_0^T{V}_p^{-1}{A}_0,---\left(13\right)$

其中V_p＝diag([var_H,pvar_L,p])，var_H,p，var_L,p分别为高、低能图像中p点的噪声方差；

保真项的梯度表示为：

${\stackrel{\·}{L}}_p=2A_0^T{V}_p^{-1}{A}_0{\stackrel{\→}{x}}_p{|}_{{\stackrel{\→}{x}}_p={\stackrel{\→}{x}}_p^{\left(n\right)}}-2A_0^T{V}_p^{-1}{\stackrel{\→}{\mathrm{\μ}}}_p.---\left(14\right)$

其中 $\begin{array}{cc}{\stackrel{\→}{x}}_p^{\left(n\right)}={\[x_{1p}^{\left(n\right)},\mathrm{...},x_{L_{0}p}^{\left(n\right)}\]}^T,& {\stackrel{\→}{\mathrm{\μ}}}_p={\left[\begin{array}{cc}{\mathrm{\μ}}_{H,p}& {\mathrm{\μ}}_{L,p}\end{array}\right]}^T.\end{array}$

所述正则项的二次替代函数如下：

$\begin{array}{c}{R}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\\ +\frac{1}{2}{(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})}^T{D}_R^{\left(n\right)}(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\end{array}---\left(15\right)$

$D_R^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{D_{R,p}^{\left(n\right)}\right\},---\left(16\right)$

$D_{R,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{4{\mathrm{\β}}_l\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ω}}_{{\mathrm{\ψ}}_l}(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)})\right\},---\left(17\right)$

其中 ${\mathrm{\ω}}_{{\mathrm{\ψ}}_l}\left(t\right)\stackrel{\mathrm{\Δ}}{=}{\stackrel{\·}{\mathrm{\ψ}}}_l\left(t\right)/t.$

正则项的梯度表示为：

${\stackrel{\·}{R}}_p=\left({\mathrm{\β}}_1\frac{\∂}{\∂x_{1p}}{R}_1\right({\stackrel{\→}{x}}_1){|}_{{\stackrel{\→}{x}}_1={\stackrel{\→}{x}}_1^{\left(n\right)}};\mathrm{...};{\mathrm{\β}}_{L_0}\frac{\∂}{\∂x_{L_{0}p}}{R}_{L_0}({\stackrel{\→}{x}}_{L_0}\left){|}_{{\stackrel{\→}{x}}_{L_0}={\stackrel{\→}{x}}_{L_0}^{\left(n\right)}}\right)---\left(18\right)$

$\frac{\∂}{\∂x_{lp}}{R}_l\left({\stackrel{\→}{x}}_l\right){|}_{{\stackrel{\→}{x}}_l={\stackrel{\→}{x}}_l^{\left(n\right)}}=\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\stackrel{\·}{\mathrm{\ψ}}}_l(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)}),l=1,\mathrm{...},L_0.---\left(19\right)$

进一步地，所述步骤(4-3)的最大预设迭代次数的设置具体为：使得分解的基材料图像的降噪程度和空间分辨率保持达到最优。

进一步地，所述步骤(4-3)p像素点图像域双能CT多材料统计分解方法的目标函数表示如下：

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\≡{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right){\stackrel{\→}{x}}_p+\frac{1}{2}\left|\right|{\stackrel{\→}{x}}_p-{\stackrel{\→}{x}}_p^{\left(n\right)}|{|}_{D_p^{\left(n\right)}}^2---\left(20\right)$

其中 $D_p^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}{D}_{L,p}^{\left(n\right)}+D_{R,p}^{\left(n\right)}.$

进一步地，所述步骤(4-3)对p像素点的目标函数按如下方式求解：

$({\widehat{\stackrel{\→}{x}}}_p,\widehat{\mathrm{\τ}})=\underset{{\stackrel{\→}{x}}_p,\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)$

$s.t.\left\{\begin{array}{c}\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1\\ 0\≤x_{lp}\≤1,l\∈\mathrm{\τ}\\ x_{lp}=0,l\∉\mathrm{\τ}\end{array}\right.---\left(21\right)$

其中

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^{T}H{\stackrel{\→}{x}}_p+{\stackrel{\→}{q}}^T{\stackrel{\→}{x}}_p$

$H\stackrel{\mathrm{\Δ}}{=}{D}_p^{\left(n\right)}$

$\begin{array}{c}\stackrel{\→}{q}\stackrel{\mathrm{\Δ}}{=}{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right)-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}\\ ={\stackrel{\·}{L}}_p+{\stackrel{\·}{R}}_p-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}\end{array}.---\left(22\right)$

τ为材料库Ω里的一个三材料组合，即τ∈Ω，Ω为包含所有预选材料三材料组合的一个三材料组合库；这个最优的目的是估计出最优的和材料类型

进一步地，所述p像素点目标函数的求解具体如下：

对每个τ∈Ω，找到最优的和相关的函数值对一个给定的τ＝(i,j,k)，其中i，j，k表示L₀材料类型里的3种材料，最优公式如下：

${\widehat{\stackrel{\→}{x}}}_p\left(\mathrm{\τ}\right)=\underset{{\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right(\mathrm{\τ}\left)\right)$

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right(\mathrm{\τ}\left)\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^T\left(\mathrm{\τ}\right)H\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)+{\stackrel{\→}{q}}^T\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)---\left(23\right)$

$s.t.\{\begin{array}{c}\underset{l}{\Σ}{x}_{lp}=1\\ a_l\≤x_{lp}\≤b_l\end{array},l=i,j,k$

其中 $\begin{array}{cc}{\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip},x_{jp},x_{kp}\]}^T,& {\stackrel{\→}{x}}_p^{\left(n\right)}\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip}^{\left(n\right)},x_{jp}^{\left(n\right)},x_{kp}^{\left(n\right)}\]}^T,\end{array}$ H(τ)和是H和中与τ＝(i,j,k)相关的值；通过比较所有的来决定最优的公式如下：

$\widehat{\mathrm{\τ}}=\underset{\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\widehat{\stackrel{\→}{x}}}_p\right(\mathrm{\τ}\left)\right).---\left(24\right)$

对l∈τ， $\begin{array}{ccc}{\widehat{\stackrel{\→}{x}}}_{lp}\≡{\widehat{\stackrel{\→}{x}}}_{lp}\left(\widehat{\mathrm{\τ}}\right);& l\∉\mathrm{\τ},& {\widehat{\stackrel{\→}{x}}}_{lp}\≡0;\end{array}$

由于目标函数是逐像素点可分离的，本发明采用并行的方式同时更新所有像素点的基材料值，即得到基于统计方式分解的基材料图像。

与现有技术相比，本发明的图像域双能CT多材料统计分解方法得到的基材料图像与直接分解的基材料图像相比，可以大幅度地降低基材料图像的噪声水平，并且保持很好的空间分辨率，且对分解精确度有极小的影响。

附图说明

图1为本发明实施例的图像域双能CT多材料统计分解方法的流程图；

图2为数字模体的二维CT图像，(a)为低能(70kVp)，(b)为高能(125kVp)，显示窗为[0.010.05]mm^‐1；

图3为数字模体高、低能二维CT图像直接多材料分解的基材料图像，本实施例中取L₀＝4，(a)为骨图像，(b)为碘图像，(c)为软组织图像，(d)为空气图像，显示窗[0.21]；

图4为数字模体高、低能二维CT图像利用本发明提出方法分解的基材料图像，本实施例中取L₀＝4，(a)为骨图像，(b)为碘图像，(c)为软组织图像，(d)为空气图像，显示窗[0.21]；

图5为模体二维CT图像，(a)为低能(75kVp)，(b)为高能(125kVp)，显示窗为[0.010.04]mm^‐1；

图6为模体高、低能二维CT图像直接多材料分解的基材料图像，本实施例中取L₀＝4，(a)为骨图像，显示窗为[0.221]；(b)为碘图像，显示窗为[0.40.7]；(c)为软组织图像，显示窗为[0.21]；(d)为空气图像，显示窗[0.21]；

图7为模体高、低能二维CT图像利用本发明提出方法分解的基材料图像，本实施例中取L₀＝4，(a)为骨图像，显示窗为[0.221]；(b)为碘图像，显示窗为[0.40.7]；(c)为软组织图像，显示窗为[0.21]；(d)为空气图像，显示窗[0.21]。

具体实施方式

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

如图1所示，本发明提供的一种图像域双能CT多材料统计分解方法，包括以下步骤：

(1)采用滤波反投影(FBP)算法根据模拟和实测高能与低能投影分别重建高、低能CT图像；所述高能为120-150kVp，低能为60-80kVp；

(2)针对高、低能CT图像，进行图像域直接多材料分解，得到直接分解的基材料图像；图像域直接多材料分解的方程如下所示：

${\stackrel{\→}{\mathrm{\μ}}}_E=\underset{l=1}{\overset{L_0}{\Σ}}{\mathrm{\μ}}_{lE}{\stackrel{\→}{x}}_l,---\left(1\right)$

其中为l材料的体积分数图像；L₀为分解的基材料种类数目；为E能量水平下的CT图像；μ_lE为E能量水平下的l材料的有效线性衰减系数；并且应该满足下面的约束和假设：

$\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1,\∀p,---\left(2\right)$

$0\≤x_{lp}\≤1,\∀l,p---\left(3\right)$

$\underset{l=1}{\overset{L_0}{\Σ}}{1}_{\{x_{lp}\≠0\}}\≤3,\∀p.---\left(4\right)$

其中p为体积分数图像的第p像素点；公式(2)为体积守恒约束；公式(3)为箱约束；公式(4)为每像素点p最多包含3种基材料的假设。

(3)针对高、低能CT图像，分别选一个均匀区域估计其方差，得到高、低能CT图像的数值方差；

(4)针对直接分解的基材料图像和高、低能CT图像的数值方差，采用基于逐像素点可分离的目标函数进行双能多材料统计分解，得到基于统计方式分解的基材料图像，具体如下：

(4-1)图像域双能CT多材料统计分解的目标函数如下所示：

$\mathrm{\Φ}\left(\stackrel{\→}{x}\right)={(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})}^T{V}^{-1}(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})+R\left(\stackrel{\→}{x}\right),---\left(5\right)$

其中A为材料成分矩阵，表示为维数为2N_p×L₀N_p，N_p为1个CT图像的像素个数，I是N_p×N_p维单位矩阵，H，L表示为高、低能谱；为2N_p×1维列向量，与分别是高、低能CT图像拉直后的列向量；T表示转置操作；为L₀N_p×1维列向量，到均为基材料拉直后的列向量；维数为2N_p×2N_p，对角元素分别是高、低能CT图像每像素的噪声方差，其中D表示对行向量进行对角线化操作，和为高、低能CT图像估计的表示为列向量的噪声方差；为正则项，公式表示如下：

$R\left(\stackrel{\→}{x}\right)=\underset{l=1}{\overset{k}{\Σ}}{\mathrm{\β}}_l{R}_l\left({\stackrel{\→}{x}}_l\right),---\left(6\right)$

$R_l\left({\stackrel{\→}{x}}_l\right)=\underset{p=1}{\overset{N_p}{\Σ}}\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ψ}}_l(x_{lp}-x_{lk})---\left(7\right)$

${\mathrm{\ψ}}_l\left(t\right)=\frac{{\mathrm{\δ}}_l^2}{3}(\sqrt{1+3{\left(\frac{t}{{\mathrm{\δ}}_l}\right)}^2}-1)---\left(8\right)$

其中x_lk是x_lp像素点的邻域，N_lp是x_lp像素点邻域的集合，k_lp和k_lk是激励均匀空间分辨率的参数；正则系数β_l和边界保持参数δ_l可以根据不同的材料来选择以实现不同的平滑正则水平；(4-2)在公式(2)-(4)的约束下最小化目标函数可以获得分解的基材料图像公式如下：

$\stackrel{\→}{x}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}subjectto\left(2\right)-\left(4\right)}{\arg \min }\mathrm{\Φ}\left(\stackrel{\→}{x}\right).---\left(9\right)$

但此目标函数很难直接最小化，本发明提出应用最优转换原理设计出一个目标函数的逐像素点可分离的二次替代函数在每次迭代过程中来单调地降低目标函数直至达到最大预设迭代次数；这个替代函数在第n次迭代过程中很容易最小化，得出的即为此替代函数的最小值，表示如下：

${\stackrel{\→}{x}}^{(n+1)}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}}{\arg }{\mathrm{min\φ}}^{\left(n\right)}\left(\stackrel{\→}{x}\right).---\left(10\right)$

分为保真项的二次替代函数和正则项的二次替代函数

所述保真项的二次替代函数如下：

$\begin{array}{c}{L}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\\ +\frac{1}{2}{(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})}^T{D}_L^{\left(n\right)}(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\end{array}---\left(11\right)$

其中表示梯度操作；是关于p＝1,...,N_p的分块矩阵，公式如下：

$D_L^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}dian\left\{D_{L,p}^{\left(n\right)}\right\},---\left(12\right)$

$D_{L,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}2A_0^T{V}_p^{-1}{A}_0,---\left(13\right)$

其中V_p＝diag([var_H,pvar_L,p])，var_H,p，var_L,p分别为高、低能图像中p点的噪声方差；

保真项的梯度表示为：

${\stackrel{\·}{L}}_p=2A_0^T{V}_p^{-1}{A}_0{\stackrel{\→}{x}}_p{|}_{{\stackrel{\→}{x}}_p={\stackrel{\→}{x}}_p^{\left(n\right)}}-2A_0^T{V}_p^{-1}{\stackrel{\→}{\mathrm{\μ}}}_p.---\left(14\right)$

其中 $\begin{array}{cc}{\stackrel{\→}{x}}_p^{\left(n\right)}={\[x_{1p}^{\left(n\right)},\mathrm{...},x_{L_{0}p}^{\left(n\right)}\]}^T,& {\stackrel{\→}{\mathrm{\μ}}}_p={\left[\begin{array}{cc}{\mathrm{\μ}}_{H,p}& {\mathrm{\μ}}_{L,p}\end{array}\right]}^T.\end{array}$

所述正则项的二次替代函数如下：

$\begin{array}{c}{R}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\\ +\frac{1}{2}{(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})}^T{D}_R^{\left(n\right)}(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)})\end{array}---\left(15\right)$

$D_R^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{D_{R,p}^{\left(n\right)}\right\},---\left(16\right)$

$D_{R,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{4{\mathrm{\β}}_l\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ω}}_{{\mathrm{\ψ}}_l}(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)})\right\},---\left(17\right)$

其中 ${\mathrm{\ω}}_{{\mathrm{\ψ}}_l}\left(t\right)\stackrel{\mathrm{\Δ}}{=}{\stackrel{\·}{\mathrm{\ψ}}}_l\left(t\right)/t.$

正则项的梯度表示为：

${\stackrel{\·}{R}}_p=\left({\mathrm{\β}}_1\frac{\∂}{\∂x_{1p}}{R}_1\right({\stackrel{\→}{x}}_1){|}_{{\stackrel{\→}{x}}_1={\stackrel{\→}{x}}_1^{\left(n\right)}};\mathrm{...};{\mathrm{\β}}_{L_0}\frac{\∂}{\∂x_{L_{0}p}}{R}_{L_0}({\stackrel{\→}{x}}_{L_0}\left){|}_{{\stackrel{\→}{x}}_{L_0}={\stackrel{\→}{x}}_{L_0}^{\left(n\right)}}\right)---\left(18\right)$

$\frac{\∂}{\∂x_{lp}}{R}_l\left({\stackrel{\→}{x}}_l\right){|}_{{\stackrel{\→}{x}}_l={\stackrel{\→}{x}}_l^{\left(n\right)}}=\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\stackrel{\·}{\mathrm{\ψ}}}_l(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)}),l=1,\mathrm{...},L_0.---\left(19\right)$

所述最大预设迭代次数的设置具体为：使得分解的基材料图像的降噪程度和空间分辨率保持达到最优。

(4-3)对p像素点图像域双能CT多材料统计分解方法的目标函数进行公式(2)-(4)约束下的求解。具体为：

p像素点图像域双能CT多材料统计分解方法的目标函数表示如下：

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\≡{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right){\stackrel{\→}{x}}_p+\frac{1}{2}\left|\right|{\stackrel{\→}{x}}_p-{\stackrel{\→}{x}}_p^{\left(n\right)}|{|}_{D_p^{\left(n\right)}}^2---\left(20\right)$

其中 $D_p^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}{D}_{L,p}^{\left(n\right)}+D_{R,p}^{\left(n\right)}.$

对p像素点的目标函数按如下方式求解：

$({\widehat{\stackrel{\→}{x}}}_p,\widehat{\mathrm{\τ}})=\underset{{\stackrel{\→}{x}}_p,\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)$

$s.t.\left\{\begin{array}{c}\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1\\ 0\≤x_{lp}\≤1,l\∈\mathrm{\τ}\\ x_{lp}=0,l\∉\mathrm{\τ}\end{array}\right.---\left(21\right)$

其中

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^{T}H{\stackrel{\→}{x}}_p+{\stackrel{\→}{q}}^T{\stackrel{\→}{x}}_p$

$H\stackrel{\mathrm{\Δ}}{=}{D}_p^{\left(n\right)}$

$\begin{array}{c}\stackrel{\→}{q}\stackrel{\mathrm{\Δ}}{=}{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right)-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}\\ ={\stackrel{\·}{L}}_p+{\stackrel{\·}{R}}_p-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}\end{array}.---\left(22\right)$

τ为材料库Ω里的一个三材料组合，即τ∈Ω，Ω为包含所有预选材料三材料组合的一个三材料组合库；这个最优的目的是估计出最优的和材料类型

p像素点目标函数的求解具体如下：

对每个τ∈Ω，找到最优的和相关的函数值对一个给定的τ＝(i,j,k)，其中i，j，k表示L₀材料类型里的3种材料，最优公式如下：

${\widehat{\stackrel{\→}{x}}}_p\left(\mathrm{\τ}\right)=\underset{{\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right(\mathrm{\τ}\left)\right)$

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right(\mathrm{\τ}\left)\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^T\left(\mathrm{\τ}\right)H\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)+{\stackrel{\→}{q}}^T\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)---\left(23\right)$

$s.t.\{\begin{array}{c}\underset{l}{\Σ}{x}_{lp}=1\\ a_l\≤x_{lp}\≤b_l\end{array},l=i,j,k$

其中 $\begin{array}{cc}{\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip},x_{jp},x_{kp}\]}^T,& {\stackrel{\→}{x}}_p^{\left(n\right)}\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip}^{\left(n\right)},x_{jp}^{\left(n\right)},x_{kp}^{\left(n\right)}\]}^T,\end{array}$ H(τ)和是H和中与τ＝(i,j,k)相关的值；通过比较所有的来决定最优的公式如下：

$\widehat{\mathrm{\τ}}=\underset{\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\widehat{\stackrel{\→}{x}}}_p\right(\mathrm{\τ}\left)\right).---\left(24\right)$

对l∈τ， $\begin{array}{ccc}{\widehat{\stackrel{\→}{x}}}_{lp}\≡{\widehat{\stackrel{\→}{x}}}_{lp}\left(\widehat{\mathrm{\τ}}\right);& l\∉\mathrm{\τ},& {\widehat{\stackrel{\→}{x}}}_{lp}\≡0\end{array};$

由于目标函数是逐像素点可分离的，本发明采用并行的方式同时更新所有像素点的基材料值，即得到基于统计方式分解的基材料图像。

实施例

1、在两个实施例中，将模拟和实测投影包括高、低能投影数据采用滤波反投影算法重建出高、低能CT图像，验证本发明相关算法的有效程度；

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

3、上述技术方案中，所述两个实施例高、低能CT图像直接多材料分解的基材料图像如图3、6，可发现分解的基材料图像中含有很大的噪声；

4、上述技术方案中，所述的数字模体高、低能CT图像利用本发明提出方法分解的图像如图4，其中骨、碘、软组织和空气的正则系数分别为10⁰，10¹，10¹，10^‐1，边界保持参数分别为0.01，0.005，0.01，0.1，权衡分解基材料图像的降噪程度和空间分辨率保持，迭代次数为200，和直接多材料分解的基材料图像的噪声水平相比(选一个均匀区域如图3(c)中矩形框所示，计算其方差值作为噪声水平)，分解的骨、碘、软组织和空气图像噪声分别降低了100％，97％，96％和95％，如图3(c)中黑色箭头所示，本发明的方法可以有效的抑制射束硬化伪影；

5、上述技术方案中，所述的模体高、低能CT图像利用本发明提出方法分解的图像如图7，其中骨、碘、软组织和空气的正则系数分别为8×10^‐1，2×10^‐1，10⁰，10¹，边界保持参数分别为0.005，0.02，0.012，0.01，权衡分解基材料图像的降噪程度和空间分辨率保持，迭代次数为200，和直接多材料分解的基材料图像的噪声水平相比(选一个均匀区域如图7(c)中矩形框所示，计算其方差值作为噪声水平)，分解的骨、碘、软组织和空气图像噪声分别降低了49％，76％，92％和96％，本发明提出的方法和直接多材料分解的方法的电子密度均方根误差分别为14.91％和14.98％，可发现本发明可以大幅度地降低基材料图像的噪声，同时对分解精确度有极小的影响；

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

Claims (10)

1.一种图像域双能CT多材料统计分解方法，其特征在于，包括以下步骤：

(1)采用滤波反投影(FBP)算法根据模拟和实测高能与低能投影分别重建高、低能CT图像；所述高能为120-150kVp，低能为60-80kVp；

(2)针对高、低能CT图像，进行图像域直接多材料分解，得到直接分解的基材料图像；

(3)针对高、低能CT图像，分别选一个均匀区域估计其方差，得到高、低能CT图像的数值方差；

(4)针对直接分解的基材料图像和高、低能CT图像的数值方差，采用基于逐像素点可分离的目标函数进行双能多材料统计分解，得到基于统计方式分解的基材料图像。

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

3.如权利要求1所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(2)图像域直接多材料分解的方程如下所示：

${\stackrel{\→}{\mathrm{\μ}}}_E=\underset{l=1}{\overset{L_0}{\Σ}}{\mathrm{\μ}}_{lE}{\stackrel{\→}{x}}_l,---\left(1\right)$

其中为l材料的体积分数图像；L₀为分解的基材料种类数目；为E能量水平下的CT图像；μ_lE为E能量水平下的l材料的有效线性衰减系数；并且应该满足下面的约束和假设：

$\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1,\∀p,---\left(2\right)$

$0\≤x_{lp}\≤1,\∀l,p---\left(3\right)$

$\underset{l=1}{\overset{L_0}{\Σ}}{1}_{\left\{x_{lp}\≠0\right\}}\≤3,\∀p.---\left(4\right)$

其中p为体积分数图像的第p像素点；公式(2)为体积守恒约束；公式(3)为箱约束；公式(4)为每像素点p最多包含3种基材料的假设。

4.如权利要求1所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(3)数值方差估计为分别在高、低能CT图像中选择一个均匀区域估计其数值方差作为CT图像每像素的噪声方差。

5.如权利要求1所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(4)采用基于逐像素点可分离的目标函数进行双能多材料统计分解通过如下方法进行：

(4-1)图像域双能CT多材料统计分解的目标函数如下所示：

$\mathrm{\Φ}\left(\stackrel{\→}{x}\right)={(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})}^T{V}^{-1}(A\stackrel{\→}{x}-\stackrel{\→}{\mathrm{\μ}})+R\left(\stackrel{\→}{x}\right),---\left(5\right)$

其中A为材料成分矩阵，表示为维数为2N_p×L₀N_p，N_p为1个CT图像的像素个数，I是N_p×N_p维单位矩阵，H，L表示为高、低能谱；为2N_p×1维列向量，与分别是高、低能CT图像拉直后的列向量；T表示转置操作；维列向量，到均为基材料拉直后的列向量；维数为2N_p×2N_p，对角元素分别是高、低能CT图像每像素的噪声方差，其中D表示对行向量进行对角线化操作，和为高、低能CT图像估计的表示为列向量的噪声方差；为正则项，公式表示如下：

$R\left(\stackrel{\→}{x}\right)=\underset{l=1}{\overset{L_0}{\Σ}}{\mathrm{\β}}_l{R}_l\left({\stackrel{\→}{x}}_l\right),---\left(6\right)$

$R_l\left({\stackrel{\→}{x}}_l\right)=\underset{p=1}{\overset{N_p}{\Σ}}\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ψ}}_l(x_{lp}-x_{lk})---\left(7\right)$

${\mathrm{\ψ}}_l\left(t\right)=\frac{{\mathrm{\δ}}_l^2}{3}(\sqrt{1+3{\left(\frac{t}{{\mathrm{\δ}}_l}\right)}^2}-1)---\left(8\right)$

其中x_lk是x_lp像素点的邻域，N_lp是x_lp像素点邻域的集合，k_lp和k_lk是激励均匀空间分辨率的参数；正则系数β_l和边界保持参数δ_l可以根据不同的材料来选择以实现不同的平滑正则水平；(4-2)在公式(2)-(4)的约束下最小化目标函数可以获得分解的基材料图像公式如下：

$\stackrel{\→}{x}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}subjectto\left(2\right)-\left(4\right)}{\arg min}\mathrm{\Φ}\left(\stackrel{\→}{x}\right).---\left(9\right)$

但此目标函数很难直接最小化，本发明提出应用最优转换原理设计出一个目标函数的逐像素点可分离的二次替代函数在每次迭代过程中来单调地降低目标函数直至达到最大预设迭代次数；这个替代函数在第n次迭代过程中很容易最小化，得出的即为此替代函数的最小值，表示如下：

${\stackrel{\→}{x}}^{(n+1)}\stackrel{\mathrm{\Δ}}{=}\underset{\stackrel{\→}{x}}{\arg }{\mathrm{min\φ}}^{\left(n\right)}\left(\stackrel{\→}{x}\right).---\left(10\right)$

分为保真项的二次替代函数和正则项的二次替代函数

(4-3)对p像素点图像域双能CT多材料统计分解方法的目标函数进行公式(2)-(4)约束下的求解。

6.如权利要求5所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(4-2)保真项的二次替代函数如下：

$\begin{array}{c}{L}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿L\left({\stackrel{\→}{x}}^{\left(n\right)}\right)\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)\\ +\frac{1}{2}{\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)}^T{D}_L^{\left(n\right)}\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)\end{array}---\left(11\right)$

其中表示梯度操作；是关于p＝1,...,N_p的分块矩阵，公式如下：

$D_L^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{D_{L,p}^{\left(n\right)}\right\},---\left(12\right)$

$D_{L,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}2A_0^T{V}_p^{-1}{A}_0,---\left(13\right)$

其中A₀＝[μ_1H,...,μ_L0H；μ_1L,...,μ_L0L]，V_p＝diag([var_H,pvar_L,p])，var_H,p，var_L,p分别为高、低能图像中p点的噪声方差；

保真项的梯度表示为：

${\stackrel{\·}{L}}_p=2A_0^T{V}_p^{-1}{A}_0{\stackrel{\→}{x}}_p{|}_{{\stackrel{\→}{x}}_p={\stackrel{\→}{x}}_p^{\left(n\right)}}-2A_0^T{V}_p^{-1}{\stackrel{\→}{\mathrm{\μ}}}_p.---\left(14\right)$

其中 ${\stackrel{\→}{x}}_p^{\left(n\right)}={\[x_{1p}^{\left(n\right)},...,x_{L_{0}p}^{\left(n\right)}\]}^T,{\stackrel{\→}{\mathrm{\μ}}}_p={\left[\begin{array}{cc}{\mathrm{\μ}}_{H,p}& {\mathrm{\μ}}_{L,p}\end{array}\right]}^T.$

所述正则项的二次替代函数如下：

$\begin{array}{c}{R}^{\left(n\right)}\left(\stackrel{\→}{x}\right)=R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)+\▿R\left({\stackrel{\→}{x}}^{\left(n\right)}\right)\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)\\ +\frac{1}{2}{\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)}^T{D}_R^{\left(n\right)}\left(\stackrel{\→}{x}-{\stackrel{\→}{x}}^{\left(n\right)}\right)\end{array}---\left(15\right)$

$D_R^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{D_{R,p}^{\left(n\right)}\right\},---\left(16\right)$

$D_{R,p}^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}diag\left\{4{\mathrm{\β}}_l\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\mathrm{\ω}}_{{\mathrm{\ψ}}_l}(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)})\right\},---\left(17\right)$

其中 ${\mathrm{\ω}}_{{\mathrm{\ψ}}_l}\left(t\right)\stackrel{\mathrm{\Δ}}{=}{\stackrel{\·}{\mathrm{\ψ}}}_l\left(t\right)/t.$

正则项的梯度表示为：

${\stackrel{\·}{R}}_p=\left({\mathrm{\β}}_1\frac{\∂}{\∂x_{1p}}{R}_1\right({\stackrel{\→}{x}}_1){|}_{{\stackrel{\→}{x}}_1={\stackrel{\→}{x}}_1^{\left(n\right)}};...;{\mathrm{\β}}_{L_0}\frac{\∂}{\∂x_{L_{0}p}}{R}_{L_0}\left({\stackrel{\→}{x}}_{L_0}\right){|}_{{\stackrel{\→}{x}}_{L_0}={\stackrel{\→}{x}}_{L_0}^{\left(n\right)}})---\left(18\right)$

$\frac{\∂}{\∂x_{lp}}{R}_l\left({\stackrel{\→}{x}}_l\right){|}_{{\stackrel{\→}{x}}_l={\stackrel{\→}{x}}_l^{\left(n\right)}}=\underset{k\∈N_{lp}}{\Σ}{k}_{lp}{k}_{lk}{\stackrel{\·}{\mathrm{\ψ}}}_l(x_{lp}^{\left(n\right)}-x_{lk}^{\left(n\right)}),l=1,...,L_0.---\left(19\right)$

7.如权利要求5所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(4-2)的最大预设迭代次数的设置具体为：使得分解的基材料图像的降噪程度和空间分辨率保持达到最优。

8.如权利要求5所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(4-3)p像素点图像域双能CT多材料统计分解方法的目标函数表示如下：

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\≡{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right){\stackrel{\→}{x}}_p+\frac{1}{2}\left|\right|{\stackrel{\→}{x}}_p-{\stackrel{\→}{x}}_p^{\left(n\right)}|{|}_{D_p^{\left(n\right)}}^2---\left(20\right)$

其中 $D_p^{\left(n\right)}\stackrel{\mathrm{\Δ}}{=}{D}_{L,p}^{\left(n\right)}+D_{R,p}^{\left(n\right)}.$

9.如权利要求5所述的图像域双能CT多材料统计分解方法，其特征在于，所述步骤(4-3)对p像素点的目标函数按如下方式求解：

$({\widehat{\stackrel{\→}{x}}}_p,\widehat{\mathrm{\τ}})=\underset{{\stackrel{\→}{x}}_p,\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)$

$s.t.\left\{\begin{array}{c}\underset{l=1}{\overset{L_0}{\Σ}}{x}_{lp}=1\\ 0\≤x_{lp}\≤1,l\∈\mathrm{\τ}\\ x_{lp}=0,l\∉\mathrm{\τ}\end{array}\right.---\left(21\right)$

其中

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^{T}H{\stackrel{\→}{x}}_p+{\stackrel{\→}{q}}^T{\stackrel{\→}{x}}_p$

$H\stackrel{\mathrm{\Δ}}{=}{D}_p^{\left(n\right)}$

$\begin{array}{c}\stackrel{\→}{q}\stackrel{\mathrm{\Δ}}{=}{\▿}_{{\stackrel{\→}{x}}_p}\mathrm{\Φ}\left({\stackrel{\→}{x}}^{\left(n\right)}\right)-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}\\ ={\stackrel{\·}{L}}_p+{\stackrel{\·}{R}}_p-{\left(D_p^{\left(n\right)}\right)}^T{\stackrel{\→}{x}}_p^{\left(n\right)}.\end{array}---\left(22\right)$

τ为材料库Ω里的一个三材料组合，即τ∈Ω，Ω为包含所有预选材料三材料组合的一个三材料组合库；这个最优的目的是估计出最优的和材料类型

10.如权利要求9所述的图像域双能CT多材料统计分解方法，其特征在于，所述p像素点目标函数的求解具体如下：

对每个τ∈Ω，找到最优的和相关的函数值对一个给定的τ＝(i,j,k)，其中i，j，k表示L₀材料类型里的3种材料，最优公式如下：

${\widehat{\stackrel{\→}{x}}}_p\left(\mathrm{\τ}\right)=\underset{{\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)\right)$

${\mathrm{\φ}}_p^{\left(n\right)}\left({\stackrel{\→}{x}}_p\right(\mathrm{\τ}\left)\right)\stackrel{\mathrm{\Δ}}{=}\frac{1}{2}{\stackrel{\→}{x}}_p^T\left(\mathrm{\τ}\right)H\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)+{\stackrel{\→}{q}}^T\left(\mathrm{\τ}\right){\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)---\left(23\right)$

$s.t.\left\{\begin{array}{c}\underset{l}{\Σ}{x}_{lp}=1\\ a_l\≤x_{lp}\≤b_l\end{array}\right.,l=i,j,k$

其中 ${\stackrel{\→}{x}}_p\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip},x_{jp},x_{kp}\]}^T,{\stackrel{\→}{x}}_p^{\left(n\right)}\left(\mathrm{\τ}\right)\stackrel{\mathrm{\Δ}}{=}{\[x_{ip}^{\left(n\right)},x_{jp}^{\left(n\right)},x_{kp}^{\left(n\right)}\]}^T,$ H(τ)和是H和中与τ＝(i,j,k)相关的值；通过比较所有的来决定最优的公式如下：

$\widehat{\mathrm{\τ}}=\underset{\mathrm{\τ}\∈\mathrm{\Ω}}{\arg }{\mathrm{min\φ}}_p^{\left(n\right)}\left({\widehat{\stackrel{\→}{x}}}_p\right(\mathrm{\τ}\left)\right).---\left(24\right)$

对l∈τ， ${\widehat{\stackrel{\→}{x}}}_{lp}\≡{\widehat{\stackrel{\→}{x}}}_{lp}\left(\widehat{\mathrm{\τ}}\right);l\∉\mathrm{\τ},{\widehat{\stackrel{\→}{x}}}_{lp}\≡0;$

由于目标函数是逐像素点可分离的，本发明采用并行的方式同时更新所有像素点的基材料值，即得到基于统计方式分解的基材料图像。