跳转至

矫正方案：水束硬化

水束硬化 (Water beam-hardening)

方案0. 硬件矫正

由 “光束硬化” 的定义得知，X射线经过物体时，整体能谱会向高能方向移动。一般使用较高密度 (如铜片) 物理滤波，使入射物体的 X 光具有较高的能量占比，从而减轻水束硬化带来的杯状伪影。

方案1. 投影域 postlog 拟合

根据上式 (2) 和图2. ，很自然的，我们想让实际的投影值曲线拟合到理想直线上，这样就能矫正水束硬化带来的杯状伪影。

Q：那么怎么求得理想直线的斜率呢？

A：\(I_p\) 在原点位置的切线满足：

\[ \begin{equation}\begin{split} \lim_{L\to 0}\frac{P(L)}{L}&=\lim_{L\to 0}-\frac{ln\frac{\int\Omega(E)e^{-\mu(E)L}dE}{\int\Omega(E)dE}}{L}\\ &=\lim_{L\to 0}-(ln{\int\Omega(E)e^{-\mu(E)L}dE})'\\ &=\lim_{L\to 0}-\frac{\int\Omega(E)e^{-\mu(E)L}(-\mu(E))dE}{\int\Omega(E)e^{-\mu(E)L}dE}\\ &=\lim_{L\to 0}\frac{\int\Omega(E)e^{-\mu(E)L}\mu(E)dE}{\int\Omega(E)e^{-\mu(E)L}dE}\\ &=\frac{\int\Omega(E)\mu(E)dE}{\int\Omega(E)dE}\\ &=\mu_{water,ref} \end{split}\end{equation}\tag{3} \]

所以，图2. 中物质水的理想直线的斜率为 \(\mu_{water,ref}=\frac{\int \Omega (E)\mu_{water}(E)dE}{\int \Omega (E) dE }\)。

细节

根据不同的探测器种类，\(\mu_{water,ref}\) 的计算并不相同。对于 能量积分型探测器 (EID)，\(\mu_{water,ref}=\frac{\int \Omega (E)\cdot E \cdot\mu_{water}(E)dE}{\int \Omega (E)\cdot E dE }\)；对于 光子计数探测器 (PCD) 或模拟投影，\(\mu_{water,ref}=\frac{\int \Omega (E)\mu_{water}(E)dE}{\int \Omega (E) dE }\)。

那么，我们只需要用多项式拟合这两者关系即可，即：

\[ \mu_{water,ref}\times L = \alpha_1P(L) + \alpha_2P^2(L) + O(P^2(L))\tag{4} \]

拟合曲线获得的系数 \(\alpha_i\) 对实际测量值 \(I_p\) 矫正，即可矫正水束硬化引起的杯状伪影。

注意

这种方案依赖能谱，当能谱未知时 \(\mu_{water,ref}\) 无法得知，拟合更无从谈起。同样，这也意味着当能谱发生变化时，拟合参数 \(\alpha_i\) 需要重新计算。（如何估计能谱？）

如图3. 左所示，为去掉头盖骨的头部模体重建图，整体值分布不均匀，杯状伪影严重。进行上述矫正后，得到如图3. 右所示，整体分布清晰且均匀。

Fig 3. 软组织硬化伪影矫正

方案2. 图像域拟合 (ECC)

图像域拟合的经典处理手段参考 Empirical cupping correction (ECC)1，这里简述处理过程：

统一符号

为与原文表述一致，该方案使用原文符号，必要部分提供解释。

ECC 方案不需要能谱，但需要额外的一组真实水模数据来做拟合，具体过程如下：

  • 仍然假设理想的投影由真实投影的多项式构成：
\[ P(q)=c_0+c_1q+c_2q^2+\cdots=\sum_nc_nP_n(q)\tag{5} \]

其中，\(q\) 代指上文postlog值 \(I_p\)，\(P_i(q)=q^i\)。

  • 不同于方案一，我们把每个幂下的真实投影值分别重建，得到：
\[ f(r)=\Re^{-1}P(q)=\sum_{n=0}^Nc_nf_n(r)\tag{6} \]

其中 \(\Re^{-1}\) 为 iradon 变换，即重建。N=4 时精度已足够。

  • 优化目标：
\[ arg\min_{c} E^2=\int w(r)[f(r)-t(r)]^2d^2r\tag{7} \]

其中 \(w(r)\) 为权重 weight 图，是一张二值图，指代需要在哪些位置做拟合。\(t(r)\) 为模板 template，在 \(w(r)\neq0\) 的位置填上该处材料的衰减系数。最终让真实重建图像尽可能贴近模板，这样的一组 \(c\) 为最终结果，作用于以后该机器的重建图像拟合参数。

Q：那么，如何求 \(c\) 呢？

A：只要对目标函数 (7) 求导即可：

\[ \nabla_c\int w(r)[f(r)-t(r)]^2d^2r=0\tag{8} \]

即求解 \(a=B\cdot c\)，其中：

\[ a_i=\int w(r)f_i(r)t(r)d^2r,\ \ B_{ij}=\int w(r)f_i(r)f_j(r)d^2r\tag{9} \]

那么，\(c=B^{-1}\cdot a\)。

示例

假设我们有水模的真实投影，只需顺序执行如下几步：

  • 根据式 (5) (6)，将水模投影 \(I_p\) 分别做多次幂后重建，得到 \(f_i(r)\)，如图4. 左，为 \(f_1(r)\)。

  • 转为二值图，使用圆形结构元 腐蚀 (如何自己实现？) 床板和水模壁，得到横断面一定是水的位置，如下图4. 中。该区域在原图的最大值作为模板 \(t_{water}(r)\) 取值。

  • 阈值分割出重建视野中 FOV 空气位置，取一个很小的数作为模板 \(t_{air}(r)\) 取值，如下图4. 右。

  • 结合式 (9)，计算 \(c\)，之后把 \(c\) 用于该 CT 机需要进行水束硬化矫正的重建图像中。

Fig 4. 杯状伪影矫正

提示

水束矫正的方案大体分为以上 投影域 与 图像域 两类，实际使用过程中可以将这些方法组合使用。例如，先做物理滤波，再结合软件方法矫正。或如果得知具体能谱，ECC 方案的 \(t(r)=\mu_{water}\) 位置，采用 \(\mu_{water,ref}\) 以取代腐蚀后，取区域最大值带来的偏差。

4. 代码实现

方案1. 投影域 postlog 拟合

杯状伪影矫正关键代码
1
pass

方案2. 图像域拟合 (ECC)

杯状伪影矫正关键代码
1
pass

完整实现 在这里

[注]：以上图片来自 ICRP 110 数字体模投影，模拟投影以及小动物 CT，仅供学习参考

  1. Kachelriess M, Sourbelle K, Kalender WA. Empirical cupping correction: a first-order raw data precorrection for cone-beam computed tomography. Med Phys. 2006 May;33(5):1269-74. doi: 10.1118/1.2188076. PMID: 16752561. ↩