17、锥形束螺旋计算机断层扫描（CBCT）中的图像重建技术

锥形束螺旋计算机断层扫描（CBCT）中的图像重建技术

1. 引言

在锥形束螺旋计算机断层扫描（CBCT）领域，图像重建是一个关键环节。为了实现高质量的图像重建，需要解决诸多技术难题，如确定重建平面的最佳参数、进行投影值的近似和校正等。本文将详细介绍CBCT中相关的重建方法和技术，包括确定重建平面的最佳倾斜角度、纵向近似方法、投影值的校正以及离散实现的ASSR算法等内容。

2. 确定重建平面的最佳倾斜角度

首先，我们要找到函数 (R_f \tan t \sin a_h) 和 (\frac{k a_h}{2\pi}) 之间的最小平均距离，对于给定参数 (a)。为此，我们将 (D_{z_{mean}}) 的导数设为零，即从方程 (\frac{dD_{z_{mean}}}{da_h}=0) 计算 (a_h) 的值。

通过一系列计算，我们得到 (a_h) 的值为：
(a_h = \frac{1}{2}(1 + \cos(a p)))

结合相关方程，我们可以得出重建平面的最佳倾斜角度 (t) 的值：
(t = \arctan\left(\frac{k - \arccos\left(\frac{1}{2}(1 + \cos a p)\right)}{2\pi R_f \sin\left(\arccos\left(\frac{1}{2}(1 + \cos a p)\right)\right)}\right))

3. 纵向近似方法

3.1 平面特性

纵向近似利用了从投影物理获得的真实射线和近似投影的射线都在与 (z) 轴平行的同一平面这一事实。若该平面用符号 (A) 表示，则有 (A \parallel z = \begin{bmatrix}0\0\1\end{bmatrix})。

平面 (A) 由向量 (A \parallel u_0(a_0p)) 和 (O + s_0 - s_0(a_0p) \in A) 唯一确定，平面内所有点满足方程：
(u_0(a_0p) \times z \cdot \begin{bmatrix}x\y\z\end{bmatrix} - (O + s_0 - s_0(a_0p)) = 0)

3.2 确定交点

平面 (A) 与移动 X 射线管的螺旋路径的交点可通过矩阵方程 (u_0(a_0p) \times z \cdot (\text{ellipse}(a_h) - (O + s_0 - s_0(a_0p))) = 0) 确定。

经过一系列推导，我们得到以下方程：
(R_f \cos a_0p \sin \hat{a}_h - \sin a_0p \cos \hat{a}_h \cos t + s_0 \cos t = 0)

进一步简化为全局坐标系下的变量形式：
(R_f \sin (\hat{a}_h - a_p) - s = 0)

由此可确定投影系统的相对旋转角度：
(\hat{a}_h(a_p, s) = a_p - \arcsin\frac{s}{R_f})

3.3 投影角度范围

对于重建过程所需的投影角度 (a_h) 的范围，如果 (a_p) 在 ([-a_p, a_p]) 范围内变化，且每个投影值所需的参数 (s) 的范围是 ([-R, R])，则需要在角度 (a_h \in [a_{h_p} - a_p + \arcsin\frac{R}{R_f}, a_{h_p} + a_p - \arcsin\frac{R}{R_f}]) 进行投影。

3.4 确定最佳射线

为了使投影值的近似误差尽可能小，我们需要选择一条与计算投影的虚拟射线路径偏差尽可能小的射线。为此，我们可以进行两个几何构造：
- 构造一 ：建立一个平面，该平面与重建切片的表面在点 (O)、(\text{ellipse}(a_h - \frac{\pi}{2})) 和 (\text{ellipse}(a_h + \frac{\pi}{2})) 相交，并包含管的位置 (\text{focus}(a_h))。该平面与重建切片的交线为：
(\text{line} 1 = O + a_1 (\text{ellipse}(a_h + \frac{\pi}{2}) - \text{ellipse}(a_h - \frac{\pi}{2})) = O - 2a_1R_f \begin{bmatrix}\cos a_h\\sin a_h\\tan t \cdot \cos (a_h - a {h_p})\end{bmatrix})
其中 (a_1) 是任意实系数。
- 构造二 ：包含点 (\text{focus}(a_h)) 和与坐标系中心 (O) 距离为 (s_0 - s(a_0p)) 的直线（由向量 (u_0(a_0p)) 表示）的平面。该平面与重建切片椭圆的交点所在直线方程为：
(\text{line}_2 = O + s_0 - s_0(a_0p) + a_2u_0(a_0p))
其中 (a_2) 是任意实系数。

两条直线的交点满足向量方程：
(O + a_1 (\text{ellipse}(a_h + \frac{\pi}{2}) - \text{ellipse}(a_h - \frac{\pi}{2})) = O + s_0 - s_0(a_0p) + a_2u_0(a_0p))

通过一系列计算，我们可以得到交点的坐标：
(\text{point} = \begin{bmatrix}x\y\z\end{bmatrix} = O + \frac{s_0 \cos t}{\cos \hat{a}_h \cos a_0p + \sin \hat{a}_h \sin a_0p \cos t} \begin{bmatrix}\cos a_h\\sin a_h\\tan t \cdot \cos \hat{a}_h\end{bmatrix})

将该点投影到屏幕上的坐标为：
(w_{\text{point}} = \frac{R_f + R_d}{R_f} \cdot \frac{-s_0 \cdot \cos t}{\cos \hat{a} h \cos a_0p + \sin \hat{a}_h \sin a_0p \cos t})
(v {\text{point}} = \frac{R_f + R_d}{R_f} \cdot \frac{s_0 \cos \hat{a}_h \sin t}{\cos \hat{a}_h \cos a_0p + \sin \hat{a}_h \sin a_0p \cos t} - \frac{k \hat{a}_h}{2\pi})

转换为全局坐标系下的形式为：
(w_{\text{point}} = -\frac{R_f + R_d}{R_f} \cdot \frac{-s}{\cos (a_h - a_{h_p} - a_p)})
(v_{\text{point}} = \frac{R_f + R_d}{R_f} \cdot \frac{s \cdot \cos (a_h - a_{h_p}) \cdot \tan t}{\cos (a_h - a_{h_p} - a_p)} - \frac{k (\hat{a} h - a {h_p})}{2\pi})

3.5 投影值的计算

平行 X 射线穿过具有衰减函数 (\lambda(x, y, z)) 的测试对象的重建平面时测量的投影值可由公式 (p_0(s_0, a_0p) = \int \lambda(x, y, z) \cdot \delta(x_0 \cos a_0p + y_0 \sin a_0p - s_0)dx_0dy_0) 确定。

为避免将重建结果转换到全局坐标系，我们需要对投影值进行一系列转换。首先，转换狄拉克δ函数的参数，然后计算积分变量的雅可比行列式。最终，积分方程变为：
(p_{0p}(s_0, a_0p) = \int \lambda(x, y, z) \cdot \frac{\sqrt{\cos^2 a_p + \cos^2 t \sin^2 a_p}}{\cos t} \cdot \delta(x \cos (a_{h_p} + a_p) + y \sin (a_{h_p} + a_p) - s)dxdy)

由于虚拟屏幕上假设的探测器密度随平行光束的入射角而变化，我们需要使用几何因子将重建切片的椭圆变形为从投影系统角度观察的圆形，得到：
(\tilde{p} p(s, a_p) = \frac{\cos t}{\sqrt{\cos^2 a_p + \cos^2 t \sin^2 a_p}} \cdot p {0p}(s_0, a_0p) = \int \lambda(x, y, z) \cdot \delta(x \cos (a_{h_p} + a_p) + y \sin (a_{h_p} + a_p) - s)dxdy)

此外，由于实际射线通过组织的路径比插值射线长，我们需要进行校正：
(p_p(s, a_p) = \text{CORRECTION} \cdot \tilde{p}_p(s, a_p))

其中校正因子 (\text{CORRECTION} = \frac{p_p(s, a_p)}{\tilde{p}_p(s, a_p)} = \frac{u}{\tilde{u}} = \cos \epsilon)，(\epsilon) 是真实射线和虚拟射线之间的角度。

校正因子的计算公式为：
(\cos \epsilon = \frac{(\text{screen}(w, v, a_h) - \text{focus}(a_h)) \cdot u_0(a_0p)}{|\text{screen}(w, v, a_h) - \text{focus}(a_h)| \cdot |u_0(a_0p)|})

经过整理，得到可计算的校正因子公式：
(\text{CORRECTION} = \cos \epsilon = \frac{w \sin (\hat{a}_h - a_p) \cos t + (R + R_d) \cos (\hat{a}_h - a_p) \cos t - v \sin a_p \sin t}{\sqrt{w^2 + v^2 + (R + R_d)^2} \cdot \sqrt{\sin^2 a_p + \cos^2 t \sin^2 a_p}})

4. 卷积/滤波和反投影重建

现在我们有了 (p_p(s, a_p))，可以进行 ASSR 重建方法的最终信号处理阶段，即使用为平行投影系统设计的任何重建方法。

在最终的反投影操作中，我们需要考虑参数 (a) 可能不等于 1 的情况，将反投影操作的关系转换为：
(\lambda(x, y) = \frac{\int_{-a_p}^{a_p} \tilde{p}_p(s, a_p)da_p}{2a})

如果扫描仪的机架倾斜，重建算法需要考虑这一点，相关实际解决方案可参考相关文章。

5. ASSR 算法的离散实现

5.1 投影系统参数

在螺旋锥形束投影系统中，ASSR 算法只能利用在特定角度获得的投影，并且仅在部分圆柱形屏幕上的特定点进行测量。

辐射束到达屏幕上的各个探测器行 (k = 1, 2, \cdots, K)（(K) 为偶数），每行中的选定射线撞击探测器，每个探测器由变量 (g = -\frac{H - 1}{2}, \cdots, 0, \cdots, \frac{H - 1}{2}) 索引（(H) 为每行探测器的奇数数量）。

在螺旋断层扫描检查期间，辐射束（实际上是其中心）仅通过 (z) 轴上的任何点一次。实际中，只进行有限数量的投影，每个投影由变量 (h = 0, \cdots, H - 1) 索引。

重建算法可用的投影值为 (\hat{p}_h(g, h, k))，范围为 (g = -\frac{H - 1}{2}, \cdots, 0, \cdots, \frac{H - 1}{2})；(h = 0, \cdots, H - 1)；(k = 1, 2, \cdots, K)。

5.2 确定初始参数

在使用算法之前，我们需要确定重建切片在 (z) 轴上的中间位置，该位置由 (z_p) 或投影系统在锥形束对称中心与之相交时的旋转角度 (a_{h_p}) 表示。可通过相关方程将 (z_p) 转换为 (a_{h_p})。

然后，使用方程 (t = \arctan\left(\frac{k - \arccos\left(\frac{1}{2}(1 + \cos a p)\right)}{2\pi R_f \sin\left(\arccos\left(\frac{1}{2}(1 + \cos a p)\right)\right)}\right)) 计算投影系统在角度 (a_{h_p}) 时重建切片绕 (y) 轴的旋转角度 (t)。

5.3 重排策略

确定初始参数 (a_{h_p}) 和 (t) 后，我们使用重排策略进行重建。将锥形束系统的投影结果与平行束系统的光栅相关联，平行束系统的探测器索引为 (l = -\frac{L - 1}{2}, \cdots, 0, \cdots, \frac{L - 1}{2})，投影索引为 (w = 0, \cdots, W - 1)。

其中 (s_l = l \cdot \Delta s)，(a_{p_w} = w \cdot \Delta p_a)，且 (L = 2 \lfloor\frac{R}{\Delta p_s}\rfloor + 1)，(W = \text{Trunc}(a \cdot 2\pi, \Delta p_a))。

5.4 获取平行投影值的步骤

1. 计算最佳角度位置 ：对于每个 (l) 和 (w) 的值，计算进行纵向近似的投影系统的最佳角度位置 (a_{h_{lw}} = w\Delta p_a - \arcsin\frac{l\Delta p_s}{R_f} + a_{h_p})。
2. 计算虚拟探测器位置 ：使用方程 (w_{lw} = \frac{R_f + R_d}{R_f} \cdot \frac{-l\Delta s}{\cos (a_{h_{lw}} - a_{h_p} - w\Delta a)}) 和 (v_{lw} = \frac{R_f + R_d}{R_f} \cdot \frac{l\Delta p_s \cdot \cos (a_{h_{lw}} - a_{h_p}) \cdot \tan t}{\cos (a_{h_{lw}} - a_{h_p} - w\Delta p_a)} - \frac{k (\hat{a} {h {lw}} - a_{h_p})}{2\pi}) 计算平面屏幕上虚拟探测器的位置。
3. 转换为实际屏幕坐标 ：将 (w_{lw}) 和 (v_{lw}) 转换为实际部分圆柱形屏幕的坐标，(b_{lw} = -\arctan\frac{w_{lw}}{R_f + R_d})，(\hat{z} {lw} = v {lw} \cdot \cos b \cdot \frac{R_f}{R_f + R_d})。
4. 定义相对值 ：定义相对值 (g_{lw} = \frac{b_{lw}}{\Delta h_b})，(h_{lw} = \frac{a_{h_{lw}}}{\Delta h_a})，(k_{lw} = \frac{K + 1}{2} + \frac{\hat{z}_{lw}}{\Delta h_z} \cdot \frac{R_f + R_d}{R_f})。
5. 定义邻域投影值 ：定义用于三线性插值的邻域投影值，如 (\hat{p} h(g {lw}^+, h_{lw}^+, k_{lw}^+)) 等八个投影值，其中 (g_{lw}^- = \text{Trunc}(g_{lw}, 1))，(g_{lw}^+ = g_{lw}^- + 1) 等。
6. 进行三线性插值 ：使用三线性插值计算投影值 (\hat{\tilde{p}} p(l, w) = \tilde{p}_p(s_l, a {p_w}))。

以下是获取平行投影值步骤的 mermaid 流程图：

graph TD;
    A[开始] --> B[计算最佳角度位置 a_{h_{lw}}];
    B --> C[计算虚拟探测器位置 w_{lw} 和 v_{lw}];
    C --> D[转换为实际屏幕坐标 b_{lw} 和 \hat{z}_{lw}];
    D --> E[定义相对值 g_{lw}, h_{lw}, k_{lw}];
    E --> F[定义邻域投影值];
    F --> G[进行三线性插值];
    G --> H[结束];

通过以上步骤，我们可以在螺旋锥形束投影系统中实现 ASSR 算法的离散实现，从而完成图像的重建。这种方法综合考虑了多种因素，如射线的几何关系、投影值的校正等，有助于提高重建图像的质量。

6. 核心公式总结

为了方便大家理解和使用上述介绍的各种方法和计算过程，下面对核心公式进行总结：

公式用途	公式内容
计算 (a_h)	(a_h = \frac{1}{2}(1 + \cos(a p)))
计算重建平面最佳倾斜角度 (t)	(t = \arctan\left(\frac{k - \arccos\left(\frac{1}{2}(1 + \cos a p)\right)}{2\pi R_f \sin\left(\arccos\left(\frac{1}{2}(1 + \cos a p)\right)\right)}\right))
确定投影系统相对旋转角度 (\hat{a}_h)	(\hat{a}_h(a_p, s) = a_p - \arcsin\frac{s}{R_f})
计算投影值 (p_0)	(p_0(s_0, a_0p) = \int \lambda(x, y, z) \cdot \delta(x_0 \cos a_0p + y_0 \sin a_0p - s_0)dx_0dy_0)
转换后的投影值 (p_{0p})	(p_{0p}(s_0, a_0p) = \int \lambda(x, y, z) \cdot \frac{\sqrt{\cos^2 a_p + \cos^2 t \sin^2 a_p}}{\cos t} \cdot \delta(x \cos (a_{h_p} + a_p) + y \sin (a_{h_p} + a_p) - s)dxdy)
变形后的投影值 (\tilde{p}_p)	(\tilde{p} p(s, a_p) = \frac{\cos t}{\sqrt{\cos^2 a_p + \cos^2 t \sin^2 a_p}} \cdot p {0p}(s_0, a_0p) = \int \lambda(x, y, z) \cdot \delta(x \cos (a_{h_p} + a_p) + y \sin (a_{h_p} + a_p) - s)dxdy)
校正后的投影值 (p_p)	(p_p(s, a_p) = \text{CORRECTION} \cdot \tilde{p}_p(s, a_p))
校正因子 (\text{CORRECTION})	(\text{CORRECTION} = \cos \epsilon = \frac{w \sin (\hat{a}_h - a_p) \cos t + (R + R_d) \cos (\hat{a}_h - a_p) \cos t - v \sin a_p \sin t}{\sqrt{w^2 + v^2 + (R + R_d)^2} \cdot \sqrt{\sin^2 a_p + \cos^2 t \sin^2 a_p}})
反投影操作计算 (\lambda(x, y))	(\lambda(x, y) = \frac{\int_{-a_p}^{a_p} \tilde{p}_p(s, a_p)da_p}{2a})
计算 (a_{h_{lw}})	(a_{h_{lw}} = w\Delta p_a - \arcsin\frac{l\Delta p_s}{R_f} + a_{h_p})
计算 (w_{lw})	(w_{lw} = \frac{R_f + R_d}{R_f} \cdot \frac{-l\Delta s}{\cos (a_{h_{lw}} - a_{h_p} - w\Delta a)})
计算 (v_{lw})	(v_{lw} = \frac{R_f + R_d}{R_f} \cdot \frac{l\Delta p_s \cdot \cos (a_{h_{lw}} - a_{h_p}) \cdot \tan t}{\cos (a_{h_{lw}} - a_{h_p} - w\Delta p_a)} - \frac{k (\hat{a} {h {lw}} - a_{h_p})}{2\pi})

7. 关键技术点分析

7.1 纵向近似的意义

纵向近似利用射线与 (z) 轴平行平面的特性，通过确定平面与螺旋路径的交点，为投影角度范围的确定提供了依据。这种近似方法有助于在有限的投影角度下，更准确地获取重建所需的信息，减少不必要的投影操作，提高重建效率。

7.2 投影值校正的必要性

由于实际射线和插值射线在穿过组织时路径长度不同，以及虚拟屏幕探测器密度随入射角变化等因素，会导致投影值存在误差。通过对投影值进行校正，如使用几何因子变形椭圆和引入余弦校正因子，可以有效减小这些误差，提高重建图像的质量。

7.3 三线性插值的作用

在离散实现 ASSR 算法时，由于实际射线参数很难与计算值完全匹配，使用三线性插值可以从邻域投影值中获取更准确的投影值。这种插值方法考虑了多个方向的信息，能够在一定程度上弥补数据的不足，使重建结果更加精确。

8. 实际应用建议

8.1 初始参数的确定

在实际应用中，准确确定重建切片在 (z) 轴上的中间位置 (z_p) 或投影系统的旋转角度 (a_{h_p}) 是关键。可以通过相关的测量设备和计算方法，确保这些初始参数的准确性，为后续的重建操作奠定基础。

8.2 投影角度的选择

根据投影角度范围的计算结果，合理选择投影角度进行投影操作。避免投影角度范围过大或过小，以保证获取足够的有效信息，同时减少不必要的辐射剂量。

8.3 校正因子的计算

在计算校正因子时，要准确确定实际射线和虚拟射线的方向和角度。可以通过向量运算等方法，确保校正因子的计算准确，从而有效减小投影值的误差。

9. 总结

本文详细介绍了锥形束螺旋计算机断层扫描（CBCT）中的图像重建技术，包括确定重建平面的最佳倾斜角度、纵向近似方法、投影值的校正以及 ASSR 算法的离散实现等内容。通过对这些方法和技术的应用，可以在螺旋锥形束投影系统中实现高质量的图像重建。

在实际应用中，需要注意初始参数的确定、投影角度的选择和校正因子的计算等关键环节，以确保重建结果的准确性和可靠性。同时，随着技术的不断发展，这些方法和技术也可能会得到进一步的优化和改进，为医学影像等领域的发展提供更有力的支持。

以下是整个重建过程的 mermaid 流程图：

graph TD;
    A[开始] --> B[确定初始参数 a_{h_p} 和 t];
    B --> C[使用重排策略关联投影结果];
    C --> D[获取平行投影值];
    D --> E[进行卷积/滤波和反投影重建];
    E --> F[输出重建图像];
    F --> G[结束];

希望本文的介绍能够帮助读者更好地理解和应用 CBCT 中的图像重建技术，为相关领域的研究和实践提供有益的参考。