数字图像处理(4版)——第 5 章——图像恢复和重建(下)(Rafael C.Gonzalez&Richard E. Woods)

原创于 2026-04-22 17:38:21 发布 · 466 阅读

5 ·

本内容遵循CC 4.0 BY-SA版权协议

GEO检测

标签

#图像处理 #计算机视觉 #图像恢复 #图像重建

数字图像处理专栏收录该内容

55 篇文章

订阅专栏

5.7 逆滤波(Inverse filtering)

5.8 最小化均方误差(Wiener)滤波(Minimum mean square error (Wiener) filtering)

5.9 约束最小平方滤波(Constrained least squares filtering)

5.10 几何均值滤波(Geometric mean filter)

5.11 基于投影的图像重建(Image reconstruction from projections)

5.11.1 简介(Introduction)

5.11.2 X射线计算机分片扫描(CT)原理(Principles of X-ray computed tomography (CT))

5.11.3 投影和Radon变换(Projections and the Radon transform)

5.11.4 反向投影(Backprojections)

5.11.5 Fourier分片定理(The Fourier-slice theorem)

5.11.6 使用平行束滤波反投影进行图像重建(Reconstruction using parallel-beam filtered backprojections)

5.11.7 使用扇形束滤波反投影进行图像重建(Reconstruction using fan-beam filtered backprojections)

5.7 逆滤波(Inverse filtering)

本节内容是我们研究图像恢复的第一步，该图像恢复由给定的退化函数 ℋ 实现，或者可以通过上一节讨论的方法获得。最简单的恢复方法是直接逆滤波，即通过将退化图像的变换 G( u , v) 除以退化转移函数，计算出原图像变换的估计值 $\hat{F}(u,v)$ ：

(5-78) $\displaystyle \hat{F}(u,v)=\frac{G(u,v)}{H(u,v)}$

除法是逐元素进行的，如 2.6 节所定义，且与公式 (5-65) 相关。在公式 (5-78) 中用公式 (5-2) 的右边替代 G( u , v ) 得到

(5-79) $\displaystyle \hat{F}(u , v) = F(u , v) + \frac{N( u ,v)}{H( u ,v)}$

这是一个有趣的表达式。它告诉我们，即使我们知道退化函数，我们也无法恢复未退化的图像( F(u , v) 的逆 Fourier 变换），因为 N(u , v) 是未知的。还有更坏的消息。如果退化函数为零或非常小，那么 N(u , v) 与 H(u , v) 的比值很容易超过 F (u , v) 项。事实上，这种情况经常发生，你很快就会看到。

解决零值或小值问题的一种方法是将滤波器频率限制在原点附近的值。根据公式 (4-92) 的讨论，我们知道 H (0,0) 通常是 H(u , v) 在频域中的最高值。因此，通过将分析限制在原点附近的频率，可以降低遇到零值的可能性。以下示例说明了这种方法。

例子 5.9：通过逆滤波对图像进行去模糊。

图 5.25(b) 中的图像通过公式 (5-78) 进行逆滤波，其滤波结果是生成该图像的退化函数的逆函数。也就是说，所使用的退化函数是

$\displaystyle H(u , v) = \displaystyle e^{\displaystyle -k[(u+M/2)^{2}+(v-N/2)^{2}]^{5/6} }$

其中，k = 0.0025 。M/2 和 N/2 常数是偏移值；它们使函数居中，使其与中心 Fourier变换相对应，如上一章所述。(记住，我们不对这些函数使用填充。)在本例中，M = N = 480。我们知道 Gauss 函数没有零点，所以这里不用担心。然而，尽管如此，退化值仍然非常小，以至于完全逆滤波的结果[图 5.27(a)]毫无用处。造成这种糟糕结果的原因已在公式 (5-79) 中讨论过。

图 5.27(b) 至 (d) 分别展示了在半径 40、70 和 85 之外截取 G(u , v) H(u , v) 比值的结果。截取是通过对比率应用 10 阶 Butterworth 低通滤波器来实现的。这在所需半径处提供了尖锐(但平滑)的过渡。半径接近 70 时，视觉效果最佳 [图 5.27(c)]。半径低于 70 时，图像会变得模糊，如图 5.27(b) 所示，该图是使用半径 40 时获得的。半径高于 70 时，图像质量开始下降，如图 5.27(d) 所示，该图是使用半径 85 时获得的。在这幅图像中，图像内容几乎可以看见，虽然被一层噪声“幕布”遮挡，但噪声仍然是最终结果的主要因素。进一步增加半径值，图像看起来越来越像图 5.27(a)。

-----------------------------图 5.27：使用公式 (5-78) 恢复图 5.25(b)。(a) 使用完整滤波器的结果。(b) H 值在 40 半径之外截止的结果。(c) H 值在70 半径之外截止的结果。(d) H 值在85 半径之外截止的结果。-----------------------------------

上述示例的结果说明了直接逆滤波通常性能不佳。接下来三个部分的基本主题是如何改进直接逆滤波。

5.8 最小化均方误差(Wiener)滤波(Minimum mean square error (Wiener) filtering)

上一节讨论的逆滤波方法没有明确处理噪声。本节将讨论一种将退化函数和噪声统计特征结合到恢复过程中的方法。该方法基于将图像和噪声视为随机变量，目标是求得未损坏图像 f 的估计值 $\hat{f}$ ，使得它们之间的均方误差最小化。该误差度量定义为

(5-80) $e^{2} = E \{ (f-\hat{f})^{2} \}$

其中 E { • } 是自变量的期望值。我们假设噪声和图像不相关，其中一方的均值为零，并且估计值的强度水平是退化图像强度水平的线性函数。基于这些假设，公式 (5-80) 中误差函数的最小值在频域中由以下表达式给出：

(5-81) $\begin{array}{rl} \displaystyle \hat{F}(u,v)&\displaystyle=\left [ \frac{H^{*}(u,v)S_{f}(u,v)}{S_{f}(u,v)|H(u,v)|^{2}+S_{{\eta}(u,v)}} \right ]G(u,v) \\ \\ &\displaystyle=\left [ \frac{H^{*}(u,v)}{|H(u,v)|^{2}+S_{{\eta}(u,v)}/S_{f}(u,v)} \right ]G(u,v) \\ \\ &\displaystyle=\left [\frac{1}{H(u,v)} \frac{|H(u,v)|^{2}}{|H(u,v)|^{2}+S_{{\eta}(u,v)}/S_{f}(u,v)} \right ]G(u,v) \end{array}$

其中，我们利用了这样一个事实：复数与其共轭的乘积等于复数的平方。这个结果称为 Wiener 滤波器，以 N.Wiener [1942] 的名字命名，他在所示年份首次提出了这一概念。由括号内的项组成的滤波器通常也称为最小均方误差滤波器或最小二乘(或平方)误差滤波器(least square error filter)(译注：“二乘法”在古汉语中指的是“两个相同的数相乘，即平方”)。我们在本章末尾提供了包含 Wiener 滤波器详细推导的资料来源。从公式 (5-81) 的第一行可以看出，Wiener 滤波器不会遇到与退化函数中为零的逆滤波器相同的问题，除非对于相同的 u 和 v 值，整个分母都为零。

公式 (5-81) 中的项解释如下：

(1) $\hat{F}(u , v)=$ 退化图像的估计之 Fourier 变换。

(2) $G(u , v) =$ 退化图像的 Fourier 变换。

(3) $H(u , v) =$ 退化转移函数(空间退化的 Fourier 变换)。

(4) $H^{*}(u , v) = H(u , v)$ 的复共轭。

(5) $|H(u,v)|^{2} = H^{*}(u , v) H(u , v)$ 。

(6) $S_{\eta}(u , v) = |N(u,v)|^{2} =$ 噪声的能量谱[见公式 (4-89)] (注：项 $|N(u,v)|^{2}$ 又称为噪声的自相关。该术语来自相关性定理(表 4.4 第 7 条第一行)。当这两个函数相同时，相关性就成了自相关，而这个条件的右边就成了 $H^{*}(u , v) H(u , v)$ , 其等于 $|H(u,v)|^{2}$ 。类似的评注适用于 $|F(u,v)|^{2}$ ，它是图像的自相关，我们将在第 12 章详细讨论相关性 )。

(7) $S_{f}(u , v) = |F(u,v)|^{2} =$ 退化函数的能量谱。

空域中的恢复图像由频域估计值 $\hat{F}(u,v)$ 的逆Fourier变换给出。注意，如果噪声为零，则噪声功率谱消失，Wiener滤波器退化为逆滤波器。另外，请记住第 5.5 节末尾关于本章所有变换工作均未进行填充的讨论。

许多有用的测量方法都基于噪声和未退化图像的功率谱。其中最重要的一个是信噪比，它使用频域量来近似，例如

(5-82) $\displaystyle SNR = \sum_{u=0}^{M-1}\sum_{v=0}^{N-1}\bigg |F(u,v) \bigg |^{2} \bigg{/}\sum_{u=0}^{M-1}\sum_{v=0}^{N-1}\bigg |N(u,v) \bigg |^{2}$

该比率衡量了信息承载信号功率(即原未降质图像的功率)与噪声功率的比值。低噪声图像往往具有较高的 SNR，反之，噪声水平较高的同一图像则具有较低的SNR。该比率是表征恢复算法性能的重要指标。

公式 (5-80) 中以统计形式给出的均方误差也可以通过原图像和恢复图像的总和来近似表示：

(5-83) $\displaystyle SNR = \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\left [f(x,y)-\hat{f}(x,y) \right ]^{2}$

事实上，如果将恢复的图像视为“信号”，并将该图像与原图像之间的差异视为“噪声”，我们可以将空间域中的信噪比定义为

(5-84) $\displaystyle$ $\displaystyle SNR = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\left [\hat{f}(x,y) \right ]^{2}\bigg{/} \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\left [f(x,y)-\hat{f}(x,y) \right ]^{2}$

f 和 $\hat{f}$ 越接近，信噪比越大。有时也会使用前两个指标的平方根，在这种情况下，它们分别称为均方根误差比(the root-mean-square-error ratio)和均方根信噪比(the root-mean-square-signal-to-noise ratio)。正如我们之前提到的，请记住，定量指标并不一定与感知图像质量有很好的相关性。处理白噪声时，频谱是一个常数，这大大简化了计算过程。然而，未退化图像的功率谱很少为人所知。当这些量未知或无法估计时，一种常用的方法是将公式 (5-81) 近似为

(5-85) $\displaystyle F (u , v) = \left [\frac{1}{H(u ,v)} \frac{|H(u,v)|^{2}}{|H(u,v)|^{2}{}+K}\right ]G(u , v)$

其中，K 是一个指定的常量，其加入到 $|H(u,v)|^{2}$ 的所有其中。下面的例子说明了这个表达式的用法。

例子 5.10：逆滤波和 Wiener 滤波去模糊的比较。

图 5.28 展示了 Wiener 滤波相对于直接逆滤波的优势。图 5.28(a) 是图 5.27(a) 的完整逆滤波结果。类似地，图 5.28(b) 是图 5.27(c) 的径向限制逆滤波结果。为了便于比较，此处复制了这些图像。图 5.28(c) 显示了使用公式 (5-85) 和示例 5.9 中使用的退化函数获得的结果。K 值是交互式选择的，以获得最佳视觉效果。Wiener 滤波相对于直接逆方法的优势在本例中显而易见。通过比较图 5.25(a) 和 5.28(c)，我们发现Wiener 滤波得到的结果在外观上非常接近原始的、未退化的图像。

-------------------图 5.28：逆滤波与 Wiener 滤波的比较。(a)图 5.25(b) 的完全逆滤波结果。(b) 径向限制逆滤波结果。(c) Wiener 滤波结果。--------------------

例子 5.11：更多使用 Wiener 滤波的去模糊示例。

图 5.29 的第一行从左到右依次显示了图 5.26(b) 的模糊图像(该图像受到均值为零、方差为 650 的加性 Gauss 噪声严重破坏)；直接逆滤波的结果；以及 Wiener 滤波的结果。文中使用了公式 (5-85) 中的 Wiener 滤波器，H(u ,v) 来自示例 5.8，K 值以交互方式选择，以获得最佳的视觉效果。正如预期的那样，直接逆滤波产生了无法使用的图像。请注意，逆滤波图像中的噪声非常强，以至于完全掩盖了图像的内容。Wiener滤波的结果并非完美，但它确实为我们提供了关于图像内容的提示。只需稍加努力即可阅读文本。

图 5.29 的第二行显示了刚才讨论的相同序列，但噪声方差降低了一个数量级。这种降低对逆滤波器的影响很小，但 Wiener 滤波结果得到了显著改善。例如，文本现在更容易阅读了。在图 5.29 的第三行中，噪声方差比第一行降低了五个数量级以上。事实上，图 5.29(g) 中的图像没有可见的噪声。在这种情况下，逆滤波的结果很有趣。噪声仍然很明显，但文本可以透过一层噪声“幕布”看到(参见问题 5.30)。图 5.29 (i) 中的 Wiener 滤波结果非常出色，在视觉上与图 5.26(a) 中的原始图像非常接近。在实践中，恢复滤波的结果很少如此接近原始图像。本例以及下一节中的示例 5.12 略微理想化，以重点关注噪声对恢复算法的影响。

---------------------------图 5.29：(a) 受运动模糊和加性噪声影响的 8 位图像。( b ) 逆滤波结果。(c) 维纳滤波结果。( d )–( f ) 相同序列，但噪声方差减小一个数量级。( g )–( i ) 相同序列，但噪声方差比 (a) 降低了五个数量级。请注意 (h) 中去模糊后的图像如何透过噪声“幕布”清晰可见。----------------------------------

5.9 约束最小平方滤波(Constrained least squares filtering)

本章讨论的所有方法都存在一个共同的问题：需要了解退化函数 H 的一些信息。然而，Wiener 滤波器还存在一个额外的难题：还必须知道未退化图像和噪声的功率谱。我们在上一节中表明，在某些情况下，使用公式 (5-85) 中的近似值可以获得可接受的结果，但功率谱比的常数值并不总是合适的解决方案。

本节讨论的方法只需要知道噪声的均值和方差即可。如第 5.2 节所述，这些参数通常可以根据给定的退化图像计算得出，因此这是一个重要的优势。另一个区别是，Wiener滤波器基于最小化统计准则，因此它在平均意义上是最优的。本节介绍的算法有一个显著的特点，即它能为应用它的每幅图像提供最优结果。当然，重要的是要记住，这些最优准则虽然从理论角度来看令人欣慰，但与视觉感知的动态无关。因此，选择哪种算法几乎总是取决于对最终图像的感知视觉质量。

通过使用公式 (4-94) 中给出的卷积定义，以及如第 2.6 节中所述，我们可以将公式 (5-64) 表示为向量矩阵形式：

(5-86) $g = Hf + \eta$

例如，假设 g(x, y) 的大小为 M × N 。我们可以利用 g(x, y) 第一行的图像元素，以及第二行的后 N 个元素，以此类推，构造向量 g 的前 N 个元素。所得向量的维数为 MN × 1。由于 f 和 H 的构造方式相同，因此维数也为 MN × 1。矩阵 H 的维数为 MN × MN 。其元素由公式 (4-94) 中卷积的元素给出。

得出这样的结论是合理的：恢复问题现在可以简化为简单的矩阵运算。遗憾的是，事实并非如此。例如，假设我们处理的是中等大小的图像，例如 M = N = 512。那么公式 (5-86) 中的向量维数为 262,144 × 1，矩阵 H 维数为 262,144 × 262,144。处理如此大小的向量和矩阵并非易事。由于 H 对噪声高度敏感(鉴于我们在前两节中对噪声影响的经验，这一点并不令人意外)，问题变得更加复杂。以矩阵形式表示恢复问题的关键优势在于它有助于推导恢复算法。

虽然我们尚未完全推导即将介绍的约束最小二乘法(最小平方法)，但该方法的根源在于矩阵公式。我们在本章末尾提供了详细介绍推导过程的参考资料。该方法的核心是 H 对噪声的敏感性问题。降低噪声敏感性影响的一种方法是基于平滑度度量(例如图像的二阶导数，我们的老朋友，Laplace 算子)来实现复原的最优性。为了使复原有意义，必须受所研究问题的参数约束。因此，我们需要找到一个准则函数 C 的最小值，其定义为

(5-87) $\displaystyle C = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}[\nabla^{2}f(x,y)]^{2}$

其遵从约束

(5-88) $\displaystyle \Vert g-H\hat{f}{\Vert}^{2}=\Vert \eta \Vert^{2}$

其中， ${\Vert}a\Vert^{2} \triangleq a^{T}a$ 是一个 Euclid 范数(见 2.6 节 )，而 $\hat{f}$ 是这个未退化图像的估算值。Laplace 算子 $\nabla^{2}$ 的定义见公式 (3-50) 。

该优化问题的频域解为表达式

(5-89) $\begin{array}{rl} \displaystyle \hat{F}(u,v)&\displaystyle=\left [ \frac{H^{*}(u,v)}{|H(u,v)|^{2}+{\gamma}|P(u,v)|^{2} } \right ]G(u,v) \end{array}$

其中，参数 γ 必须调整以满足公式 (5-88) 中的约束。P (u ,v ) 是函数

(5-90) $\displaystyle p(x,y)=\begin{bmatrix} 0&-1&0\\ -1&4&-1\\ 0&-1&0 \end{bmatrix}$

的 Fourier 变换。我们将这个函数视为来自图 3.45 的一个 Laplace 核。注意，若 γ = 0 ，则公式 (5-89) 简化为逆滤波。

函数 P(u ,v ) 和 H(u ,v ) 的大小必须相同。如果 H 的大小为 M × N，则意味着 p (x ,y ) 必须嵌入到 M × N 零数组的中心。为了保持 p (x ,y ) 的偶对称性，M 和 N 必须是偶数，如示例 4.10 和 4.15 所述。如果用于计算 H 的给定退化图像的维度不是偶数，则在计算 H 以用于公式 (5-89) 之前，必须根据需要删除相应的行和/或列。

例子 5.12：Wiener 去模糊和约束最小二乘滤波的比较。

图 5.30 展示了使用约束最小二乘滤波器处理图 5.29(a)、(d) 和 (g) 的结果，其中 g 的值是手动选择的，以获得最佳视觉效果。这与生成图 5.29(c)、( f ) 和 (i) 中 Wiener 滤波器结果的步骤相同。通过比较约束最小二乘和Wiener滤波器的结果，我们发现前者在高噪声和中噪声情况下获得了更好的结果(尤其是在降噪方面)，而两种滤波器在低噪声情况下产生的结果基本相同。这并不奇怪，因为公式 (5-89) 中的参数 γ 是一个真正的标量，而公式 (5-85) 中的 K 值是两个未知频域函数之比的标量近似值，这两个函数的大小为 M × N 。因此，基于手动选择 γ 的结果理所当然地会是对未退化图像的更精确估计。与示例 5.11 一样，本例的结果比实际应用中通常得到的结果要好。我们这里关注的是噪声模糊对恢复的影响。如前所述，你可能会遇到恢复方案与原始图像不太接近的情况，正如我们在这两个示例中所展示的那样。

-------------------------------图 5.30：约束最小二乘滤波的结果。将 (a)、(b) 和 (c) 分别与图 5.29( c )、( f ) 和 ( i ) 中的维纳滤波结果进行比较。-----------------------------------------------------

如上例所述，可以交互地调整参数 γ，直到获得可接受的结果。但是，如果我们追求数学最优性，则必须调整该参数以满足公式 (5-88) 中的约束条件。通过迭代计算 γ 的过程如下。

基于公式 (5-89) ，我们定义一个“留数(residual)”向量 r 为

(5-91) $\mathbf{r} =\mathbf{g} - \mathbf{H}\hat{f}$

我们看到， $\hat{F}(u,v)$ ( 并意味着 $\hat{f}$ ) 是 γ 的一个函数。则可以推导出 r 也是这个参数的一个函数。可以证明(Hunt [1973]、Gonzalez 和 Woods [1992])

(5-92) $\phi( \gamma ) = \mathbf{r}^{T} \mathbf{r} = \Vert \mathbf{r} \Vert^{2}$

是 γ 的单调递增函数。我们想做的事就是调整 γ 使

(5-93) ${\Vert}\mathbf{r}\Vert^{2} = {\Vert}\mathbf{\eta}\Vert^{2} {\pm}\alpha$

其中， α 是一个精确因子。根据公式 (5-91), 若 ${\Vert}\mathbf{r}\Vert^{2} = {\Vert}\mathbf{\eta}\Vert^{2}$ ，则(5-88) 中的这个约束将得到严格满足。

因为 𝜙( γ ) 是单调的，因此不难求得 γ 的预期值。一种方法为

(1) 指定 γ 的一个初始值。

(2) 计算 $\Vert \mathbf{r} \Vert^{2}$ 的值。

(3) 若满足公式 (5-93) 则中止计算；否则，若 $\Vert\mathbf{r}\Vert^{2} < ( \Vert\mathbf{\eta}\Vert^{2} - \alpha )$ 则递增 γ ，或若 $\Vert\mathbf{r}\Vert^{2} > ( \Vert\mathbf{\eta}\Vert^{2} - \alpha )$ 则递降 γ ，然后回到第 (2) 步。使用 γ 的新值在公式 (5-89) 中计算最优估计 $\hat{F}(u,v)$ 。

可以使用其他程序(例如 Newton-Raphson算法)来提高收敛速度。

为了使用这种算法，我们需要量 $\Vert\mathbf{r}\Vert^{2}$ 和 $\Vert\mathbf{\eta}\Vert^{2}$ 。为了计算 $\Vert\mathbf{r}\Vert^{2}$ ，我们从公式 (5-91) 注意到，

(5-94) $R(u ,v) = G(u , v) - H(u , v)F(u , v)$

通过计算 R(u ,v) 的逆 Fourier 变换，我们得到 r (x, y)。然后，根据 Euclid 范数的定义，可以推断出

(5-95) $\displaystyle \Vert \mathbf{r}\Vert^{2} = \mathbf{r}^{T}\mathbf{r} = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}r^{2} f(x,y)$

$\Vert \mathbf{\eta}\Vert^{2}$ 的计算导出了一个有趣的结论。首先，考虑整个图像的噪声方差，我们使用表达式

(5-96) $\displaystyle \sigma_{\eta}^{2} = \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\left [\eta(x,y)-\overline{\eta} \right ]^{2}$

从样本中估计出来，其中，

(5-97) $\displaystyle \overline{\eta} = \frac{1}{MN}\sum_{x=0}^{M-1}\sum_{y=0}^{N-1}\eta(x,y)$

是采样均值。就公式 (5-95) 的形式而言，我们注意到，公式 (5-96) 的双重求和与 $\Vert\mathbf{\eta}\Vert^{2}$ 的平方成正比。这就导出了表达式

(5-98) $\displaystyle \Vert\mathbf{\eta}\Vert^{2}=MN[\sigma_{\eta}^{2}+\overline{\eta}^{2}]$

这是一个非常有用的结果。其告诉我们，我们仅需具备噪声均值和方差的知识，就可以估算未知量 $\Vert\mathbf{\eta}\Vert^{2}$ 。并不难估算这些量的值(见 5.2 节)，假设噪声和图像强度值不相关。这是本章讨论的所有方法的假设。

例子 5.13：包含最优约束最小二乘滤波器的迭代估计。

图 5.31(a) 显示了使用刚才描述的算法估算用于恢复图 5.25(b) 的最佳滤波器所获得的结果。γ 的初始值为 $10^{-5}$ ，用于调整 γ 的校正因子为 $10^{-6}$ ，α 的值为 0.25。指定的噪声参数与生成图 5.25(a) 时使用的相同：噪声方差为 $10^{-5}$ ，均值为零。恢复的结果与图 5.28(c) 相当，后者是通过 Wiener 滤波获得的，其中 K 是手动指定的，以获得最佳视觉效果。图 5.31(b) 显示了如果使用错误的噪声参数估计值可能发生的情况。在这种情况下，指定的噪声方差为 $10^{-2}$ ，均值为 0。在这种情况下的结果会更加模糊。

------------------------------图 5.31：(a)使用正确的噪声参数，迭代确定的约束最小平方法复原图 5.25(b)。(b)使用错误的噪声参数获得的结果。--------------------

5.10 几何均值滤波(Geometric mean filter)

可以对第 5.8 节讨论的 Wiener 滤波器进行一些推广。推广后的形式是所谓的几何均值滤波器：

(5-99) $\displaystyle \hat{F}(u,v)=\left [ \frac{H^{*}(u,v)}{|H(u,v)|^{2}} \right ]^{\alpha}\left [ \frac{H^{*}(u,v)}{|H(u,v)|^{2}+\beta\left [ \frac{S_{\eta}(u,v)}{S_{f}(u,v)} \right ]} \right ]^{1-\alpha}G(u,v)$

其中，α 和 β 是非负实常量。几何均值滤波器转移函数由方括号中的两个分别取 α 次方和 1 - α 次方的表达式构成。

当 α = 1 时，几何均值简化为逆滤波器。当 α = 0 时，此滤波器成为所谓的参数化Wiener 滤波器 ，当 β = 1 时，其简化为“标准”Wiener 滤波器。若 α = 1/2，则滤波器变为两个同次幂的量之积，这正是几何平均值的定义，因此该滤波器得其名。当 β = 1 时，随着 α 增大到 1/2 以上，滤波器的性能将更接近逆滤波器。类似地，当 α 减小到 1/2 以下时，滤波器的行为更接近 Wiener 滤波器。当 α = 1/2 且 β = 1 时，该滤波器通常被称为频谱均衡滤波器(spectrum equalization filter)。公式 (5-99) 在实现恢复滤波器时非常实用，因为它将一族滤波器组合成一个单一表达式。

5.11 基于投影的图像重建(Image reconstruction from projections)

本章前几节讨论了图像退化修复技术。本节将探讨如何从一系列投影图像中重建图像，重点关注 X 射线计算机分片扫描(CT)。X 射线计算机分片扫描是最早也是目前应用最广泛的 CT 类型，也是数字图像处理在医学领域的主要应用之一。

5.11.1 简介(Introduction)

重建问题原理很简单，可以用直接、直观的方式进行定性解释，无需使用公式(我们将在本节后面讨论数学部分)。首先，考虑图 5.32(a)，图中只有一个物体位于均匀的背景上。为了使以下解释具有物理意义，我们假设该图像是人体三维区域的横截面(cross-section)。同时假设图像中的背景代表柔软均匀的组织，而圆形物体是一个肿瘤，也是均匀的，但具有更高的 X 射线吸收特性。

接下来，假设我们从左到右(穿过图像平面)穿过一束细而平坦的 X 射线，如图 5.32(b) 所示，并假设物体对光束能量的吸收比背景吸收更多(通常情况下确实如此)。在区域的另一侧使用一排 X 射线吸收探测器，即可得到如图所示的信号(吸收曲线)，其振幅(强度)与吸收成正比。(注：本文不涉及 X 射线源和探测器的物理原理，而是侧重于 CT 图像处理方面。有关 X 射线成像物理原理的精彩介绍，请参阅 Prince 和 Links [2006]的著作。) 我们可以将信号中的任意一点视为光束中与该点空间位置对应的单条射线吸收值的总和(这种总和通常称为射线和(raysum))。此时，我们所掌握的关于物体的所有信息就是这个一维吸收信号。

我们无法仅凭一次投影来判断我们处理的是单个物体还是光束路径上的多个物体，但我们首先基于此信息创建图像来进行重建。该方法是将一维信号沿光束入射方向的反方向投影，如图 5.32(c) 所示。将一维信号反投影到二维区域的过程有时称为将投影涂抹(smearing)到该区域。就数字图像而言，这意味着在图像上复制相同的一维信号，方向垂直于光束方向。例如，图 5.32(c) 是通过在重建图像的所有列中复制一维信号而创建的。显然，上述方法称为反投影法(backprojection)。

接下来，假设我们将源探测器对的位置旋转 90°，如图 5.32(d)所示。重复上一段所述的步骤，即可得到垂直方向的反投影图像，如图5.32(e)所示。我们继续进行重建，将此结果与之前的反投影图像相加，得到图5.32( f )。现在，我们开始怀疑目标物体位于图中所示的正方形区域内，该区域的振幅是各个反投影振幅的两倍，这是因为信号被叠加了。我们应该能够通过以刚才描述的方式获取更多视图来了解更多关于目标物体形状的信息，如图 5.33 所示。随着投影数量的增加，不相交反投影的振幅强度相对于多个反投影相交区域的强度会降低。最终结果是，较亮的区域将主导结果，而几乎没有或完全没有交叉点的反投影将在图像缩放显示时逐渐融入背景。

------------------------图 5.32：(a) 包含单个物体的平坦区域。(b) 平行光束、探测器条带和所探测到的一维吸收信号轮廓。(c) 吸收轮廓反投影的结果。(d) 光束和探测器旋转 90°。(e) 反投影。( f ) (c) 和 (e) 的强度之和。反投影相交处的强度是各个反投影强度的两倍。---------------------------------------------------------------

图 5.33( f ) 由 32 个反投影构成，阐释了这一概念。然而，需要注意的是，虽然该重建图像能够较好地近似原始物体的形状，但图像仍受到“光晕(halo)”效应的模糊影响，其形成过程在图 5.33 中逐步展现。例如，图 5.33(e) 中的光晕呈现为一个“星形”，其强度低于物体，但高于背景。随着视图数量的增加，光晕的形状逐渐变为圆形，如图 5.33(f) 所示。CT 重建中的模糊问题是一个重要问题，其解决方案将在本节后续部分讨论。最后，我们从图 5.32 和图 5.33 的讨论中得出结论：相差 180°的反投影互为镜像，因此我们只需考虑圆周上半圈的角度，即可生成重建所需的所有反投影。

----------------------图 5.33：(a) 同图 5.32 (a)。(b)-(e)分别使用间隔 45°的 1 , 2 , 3 和 4 个反投影进行重建。( f )使用间隔 5.625°的 32 个反投影进行重建(注意模糊)。------------------------------------------

例子 5.14：包含两个物体的平面区域的反投影。

图 5.34 展示了使用反投影对包含两个吸收特性不同的物体(较大物体的吸收率更高)的区域进行重建的结果。图 5.34(b) 显示了使用一次反投影的结果。图中从下到上有三个主要特征：一条细长的水平灰色带，对应于较小物体未遮挡的部分；其上方一条较亮(吸收率更高)的带，对应于两个物体共享的区域；以及上方一条带，对应于椭圆形物体的其余部分。图 5.34(c) 和 (d) 分别显示了使用间隔 90° 的两个投影和间隔 45° 的四个投影进行的重建结果。这些图的解释与图 5.33(c) 至 (e) 的讨论类似。图 5.34( e ) 和 ( f ) 分别显示了使用 32 次和 64 次反投影进行的更精确的重建结果。后两个结果在视觉上非常接近，并且都存在前面提到的模糊问题。

--------------------------图 5.34：(a)两个具有不同吸收特性的物体。(b)-(d)分别使用间隔 45° 的 1 次、2 次和 4 次反投影进行重建。(e)使用间隔 5.625° 的 32 次反投影进行重建。( f )使用间隔 2.8125° 的 64 次反投影进行重建。-------------

5.11.2 X射线计算机分片扫描(CT)原理(Principles of X-ray computed tomography (CT))

与上一章讨论的 Fourier 变换一样，CT 所需的基本数学概念早在数字计算机普及之前就已存在多年。CT 的理论基础可以追溯到 Vienna 数学家 Johann Radon，他在1917年推导出了一种将二维物体沿平行光线投影的方法，这是他在线积分研究的一部分(该方法现在被为 Radon 变换，我们稍后会讨论)。45年后，Tufts 大学的物理学家 Allan M. Cormack 部分地“重新发现”了这些概念，并将其应用于 CT 。Cormack 在1963 年和 1964 年发表了他的初步研究成果，并展示了如何利用不同角度方向拍摄的 X 射线图像重建人体横截面图像。他给出了重建所需的数学公式，并构建了一个 CT 原型机来验证其想法的可行性。电气工程师 Godfrey N. Hounsfield 和他在伦敦 EMI 公司的同事们独立研究出类似的解决方案，并制造出了第一台医用 CT 机。Cormack 和Hounsfield 因对分片扫描技术在医学应用方面的贡献，共同获得了1979年 Nobel 生理学或医学奖。

X 射线计算机分片扫描的目标是通过从多个不同方向对物体进行 X 射线照射，获得物体内部结构的三维图像。想象一下传统的胸部 X 光检查，将受检者放置在 X 射线感光板上，并用锥形 X 射线束照射受检者。X 射线感光板上会生成一幅图像，图像中某一点的强度与 X 射线穿过受检者后照射到该点的能量成正比。这幅图像相当于我们在上一节讨论的投影的二维图像。我们可以对整个图像进行反投影，从而创建一个三维体数据。重复这个过程，从多个角度进行扫描，并将反投影结果相加，即可得到胸腔结构的三维图像。计算机分片扫描(tomography)尝试通过对人体进行分片扫描来获取相同的信息(或相同信息的局部信息)。然后，通过堆叠这些分片，即可获得三维图像。由于获得高分辨率分片所需的探测器数量远少于生成相同分辨率的完整二维投影所需的探测器数量，因此 CT 成像更加经济。计算负担和 X 射线剂量也相应降低，使得一维投影 CT 成为一种更实用的方法。(译注：tomography ，指“一个预先确定的平面的放射成像，截面 X 射线成像”，1935 年，源自希腊语 tomos “分片(slice)，截面(section)”(参见 tome)+ -graphy。因此，tomography 更确切说法是“分片扫描，截面扫描”。 )

第一代 (G1) CT 扫描仪采用“笔形”X 射线束和单个探测器，如图 5.35(a) 所示。对于给定的旋转角度，源/探测器对(pair)沿所示的线性方向逐步平移。通过测量探测器在每一个平移增量处的输出，生成一个投影(如图 5.32 所示)。完成一次完整的线性平移后，旋转源/探测器组件，并重复上述步骤以生成不同角度的另一个投影。对于 [0°, 180°] 范围内的所有所需角度，重复上述步骤以生成一组完整的投影图像，最终从中获得一个截面图像(三维物体的切片)(译注：可能是横截面，也可能是纵截面)，如前一节所述。通过让受试者在每次完整扫描后逐步移动，使其越过源/探测器平面(受试者头部上的十字标记s标识垂直于源/探测器平面的运动方向)，即可生成一组截面图像(切片)。将这些图像通过计算堆叠，即可生成人体某一部分的三维体积。G1 扫描仪已不再用于医学成像，但由于它们产生平行光束(如图 5.32 所示)，因此其几何结构主要用于介绍 CT 成像的基本原理，并作为推导从投影重建图像所需方程的起点。

-------------------------图 5.35：四代 CT 扫描仪。虚线箭头表示线性运动的增量。虚线箭头弧表示旋转的增量。受试者头部上的十字标记表示垂直于纸面的线性运动。(a)和(b)中的双箭头表示源/探测器单元先平移，然后再返回到其原始位置。----------------

第二代 (G2) CT 扫描仪 [图 5.35(b)] 的工作原理与第一代 (G1) 扫描仪相同，但其使用的射线束呈扇形。这使得可以使用多个探测器，从而减少了源/探测器对的平移次数。

第三代 (G3) 扫描仪相比前两代 CT几何结构有了显著改进。如图 5.35(c) 所示，G3 扫描仪采用足够长的探测器阵列(约 1000 个独立探测器)，足以覆盖更宽射线束的整个视野。因此，每一次角度增量都会产生一个完整的投影，无需像 G1 和 G2 扫描仪那样平移源/探测器对。

第四代 (G4) 扫描仪更进一步。通过使用环形探测器阵列(约 5000 个独立探测器)，只需旋转射线源即可。G3 和 G4 扫描仪的主要优势在于速度；主要缺点是成本较高且 X 射线散射较大。后者意味着，为了达到与 G1 和 G2 扫描仪相当的信噪比特征，需要更高的 X 射线剂量。

新型扫描方式正逐渐采用。例如，第五代 (G5) CT 扫描仪，也称为电子束计算机分片扫描(EBCT)扫描仪，通过使用电磁控制的电子束，消除了所有机械运动。这些电子束撞击环绕患者的钨阳极，产生 X 射线，然后将其整形为扇形束，穿过患者并激发环形探测器，这与 G4 扫描仪类似。

传统的 CT 成像方式是让患者在扫描期间保持静止，直至生成一张图像。扫描随后暂停，同时使用电动床沿垂直于成像平面的方向移动患者的位置。然后获取下一张图像，并重复此过程，直至覆盖身体的特定区域。虽然获取一张图像可能只需不到一秒，但有些检查(例如腹部和胸部扫描)需要患者在图像采集过程中屏住呼吸。完成这些检查，例如采集30 张图像，可能需要几分钟。螺旋(helical) CT(有时也称为第六代(G6)CT )的使用正在不断增加。在这种方法中，G3 或 G4 扫描仪配置了所谓的滑环，从而无需在源/探测器和处理单元之间连接电缆和信号线。源/探测器对随后连续旋转360°，同时患者沿垂直于扫描方向的轴线以恒定速度移动。由此产生连续的螺旋形数据，然后对这些数据进行处理以获得各个分片的图像。

第七代(G7)扫描仪(也称多层螺旋 CT 扫描仪)正在兴起，它采用“厚”扇形束，并结合平行探测器组，同时采集容积 CT 数据。也就是说，每次 X 射线脉冲生成的是三维截面“层”，而非单个截面图像。除了显著提高图像细节外，这种方法还具有更经济地利用 X 射线管的优势，从而降低成本并有可能减少辐射剂量。

在接下来的讨论中，我们将开发公式化图像投影和重建算法所需的数学工具。我们的重点是支撑所有已讨论的 CT 方法的图像处理基础知识。有关 CT 系统的机械特性和源/探测器特性的信息，请参阅本章末尾列出的参考文献。

5.11.3 投影和Radon变换(Projections and the Radon transform)

接下来，我们将详细阐述 X 射线计算机分片扫描背景下图像重建所需的数学原理。同样的原理也适用于其他 CT 成像方式，例如 SPECT(单光子发射分片扫描)、PET(正电子发射分布扫描)、MRI(磁共振成像)以及某些超声成像方式。

笛卡尔坐标系中的一条直线可以按其斜斜截距来描述( y = ax + b )，如图 5.36 所示，也可以按其常规表示描述：

(5-100) $x \cos(\theta) + y \sin( \theta ) = \rho$

-------------------------------图 5-36：一条直线的常规表示--------------------

如图 5.37 所示，平行光束的投影可以用一组这样的线来建模。在投影剖面中坐标 $( \rho_{j} , \theta_{k} )$ 处的任意点由沿直线 $x\cos({\theta}_{k}) + y\sin{(\theta_{k})}= \rho_{j}$ 的射线和给出。目前只处理连续量，射线和是一个线积分，即

(5-101) $\displaystyle g( \rho_{j} , \theta_{k} ) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta \bigg (x \cos({\theta}_{k} )+ y\sin(\theta_{k} )-\rho_{j} \bigg ) dxdy$

其中，我们使用了 4.5 节中所描述的脉冲 δ 的属性。换言之，公式 (5-101)的右边是零(除非 δ 参数是零)，这表明积分仅沿着直线 $x\cos({\theta}_{k}) + y\sin({\theta}_{k}) = \rho_{j}$ 进行计算。若我们考虑和的所有值，则前述公式可推广为

(5-102) $\displaystyle g( \rho , \theta) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta \bigg (x \cos({\theta})+ y\sin(\theta)-\rho \bigg ) dxdy$

这个公式( 其给出了 f ( x , y ) 沿 xy 平面中任意直线的投影(线积分))是上述提到的 Radon 变换。有时候，在公式 (5-102) 中使用符号 ℜ(f ( x , y ))或 ℜ( f )替代 g( ρ , θ ) 以表明它是 f ( x , y ) 的 Radon 变换，但公式 (5-102) 中使用的这种类型的符号更为常见(传统)。正如接下来的讨论将要阐明的那样，Radon 变换是投影重建的基石，而计算机分片扫描是其在图像处理邻域的主要应用。

-------------------------------图 5.37：平行光束的几何形状。-------------------

在离散情况下(注：在第四章中，我们非常谨慎地用 ( t, z ) 表示连续图像坐标，用 (x , y) 表示离散坐标。当时，这种区分非常重要，因为我们正在构建从连续量过渡到采样量的基本概念。然而，在本文的讨论中，我们会在连续坐标和离散坐标之间频繁切换，因此坚持之前的约定可能会造成不必要的混淆。鉴于此，同时也为了与该领域已发表的文献(例如，参见 Prince 和 Links [2006])保持一致，我们将根据上下文来判断坐标 (x , y)是连续的还是离散的。当坐标是连续的时，你会看到积分；否则，你会看到求和), 公式 (5-102) 所描述的 Radon 变换就成为

(5-103) $\displaystyle g( \rho , \theta) = \sum_{x=0}^{M-1}\sum_{y=0}^{N-1}f(x,y)\delta \bigg (x \cos({\theta})+ y\sin(\theta)-\rho \bigg )$

其中，x 和 y 在当前表示离散变量，而 M 和 N 是变换在其上进行的矩形区域的维数。若我们固定 θ 并允许 ρ 变动，我们将看到，(5-103) 只是对 f (x , y) 沿着由这两个参数值所指定的值定义的直线的像素强度值求和。遍历覆盖 M × N 区域所需的所有 ρ 值(θ 固定不变)，即可得到一个投影。改变 θ 并重复此过程，即可得到另一个投影，依此类推。图 5.32-5.34 中的投影正是以此方式生成的。

例子 5.15 利用 Radon 变换获取圆形区域的投影。

在继续之前，我们先演示如何使用 Radon 变换来获得图 5.38(a) 中圆形物体的投影的解析表达式：

$\displaystyle f(x,y)=\left \{ \begin{array}{l} A \hspace{0.2cm}(x^{2}+y^{2} \le r^{2}) \\ \\ 0 \hspace{0.2cm}({\rm{otherwise}}) \end{array} \right .$

其中，A 是常数，r 是物体的半径。我们假设圆心位于 xy 平面的原点。由于物体具有圆对称性，因此其投影在所有角度下都相同，所以我们只需计算 θ = 0° 时的投影即可。公式 (5-102) 变为：

$\displaystyle g(\rho,\theta ) = \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(x-\rho)dxdy = \int_{-\infty}^{\infty}f(\rho,y)dy$

其中第二个表达式由公式 (4-13) 得出。正如前面所述，这是一个线积分(在本例中沿曲线 L( ρ , 0) 进行积分 )。此外，请注意，当 |ρ| > r 时，g( ρ , θ ) = 0 。而当 |ρ| ≤ r 时，此积分沿从 $y = -(r^{2}-{\rho}^{2} )^{1/2}$ 到 $y = (r^{2}-{\rho}^{2} )^{1/2}$ 进行计算。因此，

$\displaystyle g(\rho,\theta)=\int_{-\sqrt{r^{2-\rho^{2}}}}^{\sqrt{r^{2-\rho^{2}}}}f(\rho,y)dy=\int_{-\sqrt{r^{2-\rho^{2}}}}^{\sqrt{r^{2-\rho^{2}}}}Ady$

计算这个积分得到

$\displaystyle g(\rho,\theta)=g(\rho)=\left \{ \begin{array}{l} 2A \sqrt{r^{2}+\rho^{2}}\hspace{0.2cm}(|\rho|\le r) \\ \\ 0 \hspace{0.2cm}({\rm{otherwise}}) \end{array} \right .$

这里我们利用了当 r > r 时 g( ρ , θ ) = 0 这一事实。图 5.38(b) 显示了这一结果的图像。注意，g( ρ , θ ) = g( ρ )；即，g 与 θ 无关，因为该物体关于原点对称。

--------------------图 5.38：(a) 一个圆盘，以及 (b) 其 Radon 变换的图像，该变换是通过解析方法推导出来的。这里我们能够绘制出变换的图像，因为它只依赖于一个变量。当 g 依赖于 ρ 和 θ 两个变量时，Radon 变换就变成了一幅图像，其坐标轴分别为 ρ 和 θ ，图像中每一个像素的强度与该像素位置处 g 的值成正比。------------------

当将 Radon 变换 g( ρ , θ ) 以 ρ 和 θ 为直角坐标显示为图像时，所得结果称为正弦图，其概念类似于显示 Fourier 频谱。与 Fourier 变换类似，正弦图包含重建 f ( x, y ) 所需的数据。然而，与Fourier变换不同的是，g( ρ , θ ) 始终是一个实函数。与Fourier 频谱的显示类似，对于简单的区域，正弦图很容易解释，但随着投影区域变得越来越复杂，正弦图也越来越难以“解读”。例如，图 5.39(b) 是左侧矩形的正弦图。垂直轴和水平轴分别对应于 θ 和 ρ 。因此，最下面一行是矩形在水平方向(即 θ = 0°)的投影，中间一行是矩形在垂直方向( θ = 90°) 的投影。最下面一行非零部分的长度小于中间一行非零部分的长度，这表明物体在水平方向上比在垂直方向上更窄。正弦图关于图像中心在两个方向上都对称，这表明我们处理的是一个对称且平行于 x 轴和 y 轴的物体。最后，正弦图是平滑的，表明物体的强度是均匀的。除了这些一般性的观察之外，我们无法对这个正弦图做出更多解释。

图 5.39(c) 是 Shepp-Logan 模型的图像(Shepp 和 Logan [1974])，这是一种广泛使用的合成图像，旨在模拟大脑主要区域(包括小肿瘤)的吸收情况。如图 5.39(d) 所示，该图像的正弦图更难解释。我们仍然可以推断出一些对称性特征，但除此之外，我们能说的就很少了。对正弦图进行视觉分析的实际应用价值有限，但它们在算法开发等任务中可能会有所帮助。

--------------------------图 5.39：两幅图像及其正弦图( Radon 变换)。正弦图的每一行代表沿垂直轴上相应角度的投影。(请注意，正弦图的水平轴表示 ρ 值。) 图像 (c) 称为 Shepp-Logan 幻影。其原始形式的对比度相当低。此处显示的是增强后的图像，以便于观察。---------------------------------------

5.11.4 反向投影(Backprojections)

为了获得 Radon 变换反投影图像的正式表达式，对于一个固定的旋转值 $\theta_{k}$ ，我们从完整投影 $g(\rho, \theta_{k})$ 的一个点 $g(\rho_{j} ,\theta_{k})$ 开始( 见图 5.37 )。通过反投影将这一个点构成图像的一部分，实际上就是将直线 $L(\rho_{j} ,\theta_{k} )$ 复制到图像上，其中，在那条直线上的每一个点的值(强度)都是 $g(\rho_{j} ,\theta_{k})$ 。对投影信号中所有 $\rho_{j}$ 值重复此过程( 但保持 θ 的值固定为 $\theta_{k}$ ) ，即可得到以下表达式：

$f(\theta_{k} )( x , y ) = g( \theta , \theta_{k}) = g( x\cos(\theta_{k}) + y \sin(\theta_k) , \theta_k)$

通过固定角度 $\theta_{k}$ 然后进行反向投影而获得的图像像如图 5.32(b) 所示。这个公式对 $\theta_{k}$ 的所有取值均成立，因此，我们通常可以将在 θ 角处进行单次反向投影而获得的图像写成

(5-104) $f_{\theta}( x , y ) = g( x\cos(\theta) + y\sin(\theta) , \theta )$

我们通过对所有反向投影图像进行积分来形成最终图像：

(5-105) $\displaystyle f( x , y ) = \int_{0}^{\pi}f_{\theta}( x ,y )d\theta$

在离散情况下，积分变成所有反投影图像的总和：

(5-106) $\displaystyle f ( x , y ) = \sum_{\theta=0}^{\pi}f_{\theta}( x ,y )$

其中，x , y , 和 θ 是离散量，如前所述，0°和 180°方向的投影互为镜像，因此求和运算只需进行到 180°之前的最后一个角度增量即可。例如，如果使用 0.5°的增量，则求和范围为0°到179.5°，以0.5°为步长。以这种方式形成的逆投影图像有时被称为层析图像。需要明确的是，层析图像只是生成投影的原始图像的近似值，这一点将在以下示例中进行说明。

例子 5.16 从正弦图中获取反投影图像。

图 5.32 至图 5.34 中的反投影图像是使用公式 (5-106) 生成的，而投影数据则来自公式 (5-103)。类似地，图 5.40(a) 和 (b) 也是使用这些公式生成的，它们分别显示了与图 5.39(b) 和 (d) 中的正弦图对应的反投影图像。与之前的图像一样，我们注意到图像存在明显的模糊，因此很明显，直接使用公式 (5-103) 和 (5-106) 无法获得令人满意的结果。早期的实验性 CT 系统就是基于这些公式构建的。然而，正如我们稍后讨论的那样，通过改进反投影方法，可以显著提高图像重建质量。

-----------------------图 5.40：图 5.39 中正弦图的反投影图像。----------------

5.11.5 Fourier分片定理(The Fourier-slice theorem)

在本节中，我们将推导出一个基本方程，该方程建立了投影的一维 Fourier 变换与获取该投影的区域的二维 Fourier 变换之间的关系。这种关系是解决我们迄今为止遇到的模糊问题的重建方法的基础。

一个关于 ρ 投影的一维 Fourier 变换是

(5-107) $\displaystyle G ( \omega , \theta ) =\int_{-\infty}^{\infty}g(\rho,\theta)e^{-j2{\pi}{\omega}\rho}d\rho$

其中，ω 是角频率变量，可以认为该表达式基于固定的 θ 值。代入 (5-102) 中的 g( ρ , θ ) ，我们得到

(5-108)

$\begin{array}{rlc} G (\omega ,\theta ) &\displaystyle=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(x\cos(\theta) + y\sin(\theta)-\rho) e^{-j2\pi\omega\rho} dxdyd\rho \\ \\ &\displaystyle= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\left [\int_{-\infty}^{\infty}\delta(x\cos(\theta) + y\sin(\theta)-\rho)e^{-j2\pi\omega\rho} \right ] dxdy \\ \\ &\displaystyle= \int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)e^{-j2{\pi}{\omega}(x\cos(\theta) + y\sin(\theta))}dxdy \end{array}$

其中，最后一步是根据第 4 章中讨论的脉冲函数的筛选属性得出的。通过令 u = ωcos(θ) ，v = ωsin(θ) ，我们可以将公式 (5-108) 写成

(5-109) $\displaystyle G ( \omega , \theta ) = \left [\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)e^{ -j2{\pi}(ux + vy ))} dxdy\right ]_{u = \omega\cos(\theta);v = \omega\sin(\theta)}$

我们可以将此表达式视为 f ( x , y ) [参见公式 (4-59)] 在指定 u 和 v 值处的二维 Fourier 变换。即，

(5-110) $\begin{array}{rlc}\displaystyle G ( \omega , \theta ) =\bigg [F(u,v)\bigg ]_{u = \omega\cos(\theta);v = \omega\sin(\theta)} = F ( \omega\cos(\theta) , \omega\sin(\theta))\end{array}$

照例，F(u,v) 表示函数 f ( x , y ) 的二维 Fourier 变换。

公式 (5-110) 中的结果称为 Fourier 分片定理(或投影分片定理)。它表明，投影的Fourier 变换是该投影所来自区域的二维 Fourier 变换的一个服片。这种术语的由来可以用图 5.41 来解释。如图所示，任意投影的一维 Fourier 变换是通过提取 F(u,v) 沿一条直线上的值得到的，这条直线的方向与生成投影时使用的角度相同。

在原则上，我们可以通过对 F(u,v) 进行 Fourier 逆变换来获得 f ( x , y )。然而，这种方法计算成本很高，因为它需要计算二维变换的逆变换。下一节讨论的方法效率要高得多。

5.11.6 使用平行束滤波反投影进行图像重建(Reconstruction using parallel-beam filtered backprojections)

如图 5.33、5.34 和 5.40 所示，直接进行反投影会产生模糊不清的结果。幸运的是，这个问题有一个简单的解决方法，即在计算反投影之前对投影进行滤波。根据公式 (4-60)，F(u,v) 的二维逆 Fourier 变换为

(5-111) $\displaystyle f ( x , y ) =\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}F(u,v)e^{j2\pi(ux + vy )}dudv$

若同 (5-109) 和 (5-110) 一样，令 u = ωcos(θ) ，v = ωsin(θ) ，且我们将公式 (5-111) 用极坐标进行描述：

(5-112) $\displaystyle f ( x , y ) =\int_{0}^{2\pi}\int_{0}^{\infty}F(\omega\cos{\theta},\omega\sin{\theta})e^{j2\pi\omega(x\cos(\theta) + y\sin(\theta) )}{\omega}{d\omega}d\theta$

则，使用 Fourier 分片定理，有

(5-113) $\displaystyle f ( x , y ) =\int_{0}^{2\pi}\int_{0}^{\infty}G(\omega,{\theta})e^{j2\pi\omega(x\cos(\theta) + y\sin(\theta) )}{\omega}{d\omega}d\theta$

通过将这个积分分裂成两个表达式，一部分 θ 位于 0 到 180° ，另一部分其位于 180° 到 360°，利用 G ( ω , θ + 180°) = G ( -ω , θ ) 的事实( 见问题 5.46 )，我们可以将公式 (5-113) 表示为

(5-114) $\displaystyle f ( x , y ) =\int_{0}^{\pi}\int_{-\infty}^{\infty}|\omega|G(\omega,{\theta})e^{j2\pi\omega(x\cos(\theta) + y\sin(\theta) )}{d\omega}d\theta$

项 x cos(θ) + y sin(θ) 是一个关于 ω 的常量，我们将其视为公式 (5-100) 中的 ρ 。因此，我们可以将公式 (5-114) 写成

(5-115) $\displaystyle f ( x , y ) =\int_{0}^{\pi}\left [\int_{-\infty}^{\infty}|\omega|G(\omega,{\theta})e^{j2\pi\omega\rho}{d\omega} \right ]_{\rho=x\cos(\theta) + y\sin(\theta)}d\theta$

内部表达式的形式为一维 Fourier 逆变换[参见公式(4-21)]，其加入项 |ω| 根据第 4.7节的讨论，我们可知其为一维滤波器转移函数。注意，|ω| 是一个斜坡函数[参见图5.42(a)]。由于其幅值在两个方向上都延伸到无穷大，因此该函数不可积，Fourier 逆变换也因此未定义。理论上，可以通过使用所谓的广义 δ 函数等方法来解决这个问题。实际上，常用的方法是对斜坡函数进行加窗处理，使其在定义的频率区间之外为零。也就是说，加窗操作对斜坡滤波器转移函数进行了频带限制。

-----------------------------------图 5.41：Fourier 分片定理示意图。一维投影的Fourier 变换是该投影所来自区域的二维 Fourier 变换的一个分片。请注意在两个图中角度 θ 的对应关系。-----------------------------------------------------------

对函数进行带宽限制的最简单方法是在频域中使用矩形窗函数。然而，正如我们在图 4.4 中看到的，矩形窗函数具有不理想的振铃特性。图 5.42(b) 和 (c) 展示了这一点。前者显示了斜坡转移函数经过矩形窗函数进行带宽限制后的曲线，后者显示了通过计算其逆Fourier变换得到的空间域表示。正如预期的那样，所得的加窗滤波器在空间域中表现出明显的振铃现象。我们从第 4 章中得知，频域滤波等效于空间域卷积，因此使用具有振铃特性的函数进行空间滤波也会产生受振铃污染的结果。使用平滑函数进行加窗处理有助于改善这种情况。用于一维 FFT 实现的常用 M 点离散窗函数为

(5-116) $\displaystyle H(\omega)=\left \{ \begin{array}{l} c+(c-1)\cos{\frac{2\pi\omega}{M}} \hspace{0.2cm} (0 \le \omega \le M-1) \\ \\ 0 \hspace{0.2cm}(\rm{otherwise}) \end{array} \right .$

当 c = 0.54 时，此函数称为 Hamming 窗(以 Richard Hamming 命名)；当 c = 0.5 时，则称为 Hann 窗(以 Julius von Hann 命名 )。Hamming 和 Hann 窗的主要区别在于后者两端的函数值为零。在图像处理应用中，两者的差异通常在视觉上难以察觉。

图 5.42(d) 是 Hamming 窗的曲线图，图 5.42(e) 显示了该窗口与图 5.42(b) 所示的带限斜坡滤波器转移函数的乘积。图 5.42( f ) 显示了该乘积在空间域的表示，通常通过计算逆快速 Fourier 变换 ( IFFT ) 获得。通过比较此图和图 5.42(c) 可以明显看出，加窗斜坡的振铃现象有所减少(图 5.42(c) 和 ( f ) 中峰谷比分别为 2.5 和 3.4)。另一方面，由于图 5.42( f ) 中主瓣的宽度略宽于图 5.42(c) 中的主瓣宽度，因此我们预计使用汉明窗进行反投影会减少振铃，但图像会略微模糊。正如下面的示例 5.17 所示，实际情况确实如此。

-------------------图 5.42：(a) 频率域斜坡滤波器转移函数。(b) 经矩形滤波器进行带限后的函数。(c) 空间域表示。(d) Hamming 窗函数。(e) 窗化斜坡滤波器，由 (b) 和 (d) 相乘得到。( f ) 乘积的空间域表示。(注意振铃现象有所减弱。)--------------

回顾一下公式 (5-107)，G ( ω , θ ) 是 g( ρ , θ ) 的一维 Fourier 变换，其中 g( ρ , θ ) 是在固定角度 u 下获得的单个投影。公式 (5-115) 表明，完整的反投影图像 f ( x, y ) 可以通过以下方式获得：

(1) 计算每一个投影的一维 Fourier 变换。

(2) 将每个一维 Fourier 变换乘以滤波器转移函数 |ω|，如上所述，该函数已乘以合适的(例如，Hamming)窗函数。

(3) 对每一个滤波后的变换结果进行一维 Fourier 逆变换。

(4) 对步骤 (3) 中得到的所有一维逆变换结果进行积分(求和)。

由于使用了滤波函数，这种图像重建方法贴切地称为滤波反投影法(filtered backprojection)。实际上，数据是离散的，因此所有频域计算都使用一维快速 Fourier 变换(FFT)算法进行，滤波则采用与第 4 章中解释的二维函数基本相同的步骤来实现。或者，我们也可以像后面解释的那样，在空间域中使用卷积来实现滤波。

前面的讨论探讨了滤波反投影的窗函数方面。与任何采样数据系统一样，我们也需要关注采样率。我们从第 4 章了解到，采样率的选择对图像处理结果有着深远的影响。在本次讨论中，有两个采样方面的考虑。第一个是使用的射线数量，它决定了每个投影中的样本数量。第二个是旋转角度增量的数量，它决定了重建图像的数量(这些图像的总和构成最终图像)。欠采样会导致混叠，正如我们在第 4 章中所看到的，混叠会在图像中表现为伪影，例如条纹。我们将在后续讨论中更详细地探讨 CT 扫描的采样问题。

例子 5.17 使用滤波反投影法进行图像重建。

本例的重点是展示使用滤波反投影进行图像重建，首先使用矩形窗函数进行滤波，然后使用 Hamming 窗函数进行滤波。这些滤波反投影的结果将与图 5.40 中的“原始”反投影结果进行比较。为了仅关注滤波带来的差异，本例中的结果是使用 0.5° 的旋转增量生成的，与生成图 5.40 时使用的增量相同。两种情况下，射线之间的间隔均为一个像素。两个示例中的图像尺寸均为 600 × 600 像素，因此对角线长度约为 $\sqrt{2} \times 600 \approx 849$ 。因此，当旋转角度为 45° 和 135° 时，使用了 849 条射线来覆盖整个区域。

图 5.43(a) 显示了使用矩形窗限带的斜坡函数重建的矩形图像。该结果最显著的特点是没有任何肉眼可见的模糊。然而，正如预期的那样，图像中存在振铃效应，表现为一些细微的线条，尤其是在矩形的角落周围。这些线条在图 5.43(c) 的放大图中更加明显。如图 5.43(b) 和 (d) 所示，使用汉明窗对斜坡函数进行处理显著改善了振铃问题，但代价是图像略微模糊。与图 5.40(a) 相比，这些改进(即使是使用矩形窗处理的斜坡函数)也是显而易见的。由于体模图像的过渡不像矩形那样尖锐和突出，因此即使使用矩形窗处理的斜坡函数，振铃效应在这种情况下也几乎察觉不到，如图 5.44(a) 所示。如图 5.44(b) 所示，使用汉明窗处理后图像更加平滑。这两个结果都比图 5.40(b) 有了显著的改进，再次说明了滤波反投影方法固有的显著优势。

在大多数 CT 应用中(尤其是在医学领域)，伪影(例如振铃伪影)是一个严重的问题，因此人们投入了大量精力来最大限度地减少这些伪影。调整滤波算法以及如前所述使用大量探测器都是有助于减少这些影响的设计考虑因素。

----------图 5.43：使用 (a) 斜坡滤波器和 (b) Hamming 窗加权斜坡滤波器对矩形进行滤波反投影。第二行显示了第一行图像的放大细节。请与图 5.40(a) 进行比较。------

--------------图 5.44：使用以下方法对头部模型进行滤波反投影：(a) 斜坡滤波器，以及 (b) Hamming 窗加权的斜坡滤波器。与图 5.40(b) 进行比较。------------------

前面的讨论基于通过快速 Fourier 变换(FFT)实现获得滤波后的反投影。然而，我们从第 4 章的卷积定理可知，使用空间卷积也可以获得等效的结果。具体来说，请注意，公式 (5-115) 中括号内的项是两个频域函数乘积的逆 Fourier 变换，根据卷积定理，我们知道它等于这两个函数的空间表示(逆 Fourier 变换)的卷积。换言之，令 s(ρ) 表示 ω 的逆 Fourier 变换，(注：如果使用窗函数(例如 Hamming 窗)，则对加窗后的斜坡信号进行Fourier 逆变换。) 我们可以将式 (5-115) 写成

(5-117) $\begin{array}{rl} f ( x , y )&=\displaystyle \int_{0}^{\pi}\left [\int_{-\infty}^{\infty}|\omega|G(\omega,{\theta})e^{j2\pi\omega\rho}{d\omega} \right ]_{\rho=x\cos(\theta) + y\sin(\theta)}d\theta \\ \\ &=\displaystyle \int_{0}^{\pi}\bigg [s(\rho){\bigstar}g(\rho,\theta)e^{j2\pi\omega\rho} \bigg ]_{\rho=x\cos(\theta) + y\sin(\theta)}d\theta \\ \\ &=\displaystyle \int_{0}^{\pi}\bigg [\int_{-\infty}^{\infty}g(\rho,\theta)s(x\cos(\theta) + y\sin(\theta)-\rho)d\rho \bigg ]d\theta \end{array}$

其中，与第 4 章一样，“★”表示卷积。第二行是根据上一段解释的原因从第一行推导出来的。第三行(包括 -ρ )是根据公式 (4-24) 中卷积的定义得出的。

公式 (5-117) 的最后两行表达的是同一个意思：角度为 θ 的单个反投影可以通过对相应的投影 g( ρ , θ ) 和斜坡滤波器转移函数的傅里叶逆变换 s(ρ) 进行卷积得到。与之前一样，完整的反投影图像是通过对所有单个反投影图像进行积分(求和)得到的。除了计算中的舍入误差之外，使用卷积方法得到的结果与使用快速 Fourier 变换 (FFT) 方法得到的结果是相同的。在实际的 CT 实现中，卷积通常在计算上更有效率，因此大多数现代 CT 系统都采用这种方法。Fourier 变换在理论公式推导和算法开发中确实发挥着核心作用(例如，MATLAB 中的 CT 图像处理就是基于 FFT 的)。此外，我们注意到在重建过程中无需存储所有反投影图像。

相反，只需用最新的反投影图像更新一个累加和即可。程序结束时，该累加和将等于所有反投影图像的总和。

最后，我们指出，由于斜坡滤波器(即使经过加窗处理)会在频域中使直流分量为零，因此每一个反投影图像的平均值都为零(参见图 4.29)。这意味着每一个反投影图像中的像素值既有正值也有负值。当所有反投影图像叠加形成最终图像时，一些负值位置可能会变为正值，平均值也可能不再为零，但通常情况下，最终图像中仍然会存在负像素值。

解决这个问题有几种方法。当不知道平均值应该是什么时，最简单的方法是接受负值是这种方法固有的事实，并使用公式 (2-31) 和 (2-32) 中描述的程序对结果进行缩放。本文采用的就是这种方法。如果已知“典型”平均值应该是什么，则可以将该值添加到频域中的滤波器转移函数中，从而抵消斜坡并防止直流项归零(参见图 4.30(c))。在空间域中使用卷积时，截断空间滤波器核(斜坡的 Fourier 逆变换)的长度本身就可以防止其平均值为零，从而完全避免了归零问题。

5.11.7 使用扇形束滤波反投影进行图像重建(Reconstruction using fan-beam filtered backprojections)

到目前为止，我们讨论的重点是平行光束。由于其简单直观，这种成像几何结构传统上被用于介绍计算机分片扫描技术。然而，更现代的 CT 系统使用扇形光束几何结构(参见图 5.35 )，这也是接下来讨论的主题。

图 5.45 显示了一种基本的扇形束成像几何结构，其中探测器排列在圆弧上，并且假设光源的角增量相等。令 p( α ,β ) 表示扇形束投影，其中 α 是特定探测器相对于中心射线的角位置，β 是光源相对于 y 轴的角位移，如图所示。我们还在图 5.45 中注意到，扇形束中的一条射线可以用标准形式的直线 L( ρ , θ ) 来表示，这与我们之前讨论的平行束成像几何结构中表示射线的方法相同。这使我们能够利用平行束的结果作为起点，推导出扇形束几何结构的相应方程。接下来，我们将基于卷积推导扇形束滤波反投影算法，以证明这一点。(注：Fourier 分片定理是针对平行光束几何结构推导出来的，不能直接应用于扇形光束。然而，公式 (5-118) 和 (5-119) 为将扇形光束几何结构转换为平行光束几何结构提供了基础，从而使我们能够使用上一节中开发的滤波平行反投影方法，而该方法正是基于分片定理。我们将在本节末尾更详细地讨论这一点。)

我们首先注意到，在图 5.45 中，直线 L( ρ , θ ) 的参数与扇形光束的参数之间存在以下关系：

(5-118) $\theta = \beta + \alpha$

和

(5-119) $\rho = D\sin(\alpha)$

其中 D 是光源中心到 xy 平面原点的距离。

-------------------------图 5.45：基本的扇形射线几何结构。穿过光源中心和原点（此处假定原点是光源的旋转中心）的直线称为中心射线。-------------------------

平行光束成像几何的卷积反投影公式由式 (5-117) 给出。不失一般性，假设我们关注以 xy 平面原点为中心、半径为 T 的圆形区域内的物体。那么，当 r > T 时，g( ρ , θ ) = 0，式 (5-117) 变为

(5-120) $\displaystyle f ( x,y ) = \frac{1}{2}\int_{0}^{2\pi}\int_{-T}^{T}g(\rho,\theta )s(x\cos(\theta)+ y\sin(\theta)-\rho)d{\rho}d\theta$

这里我们利用了前面提到的一个事实，即相隔 180° 的投影互为镜像。这样，方程 (5-120) 中外层积分的积分限就覆盖了一个完整的圆周，这符合扇形束排列的要求，因为探测器是呈圆形排列的。

我们感兴趣的是关于 α 和 β 的积分。为此，我们改用极坐标 ( r , φ )。也就是说，我们令 $x=r\cos{\varphi}$ 和 $y=r\sin{\varphi}$ ，由此可得

(5-121) $x\cos(\theta)+y\sin{\theta}=r\cos(\varphi)\cos(\theta)+r\sin(\varphi)\sin(\theta)=r\cos(\theta-\varphi)$

使用这个结果，我们可以将公式 (5-120) 写成

(5-122) $\displaystyle f ( x , y ) = \frac{1}{2} \int_{0}^{2\pi}\int_{-T}^{T}g(\rho ,\theta )s(r\cos(\theta-\varphi) -\rho)d{\rho}d\theta$

这个表达式不过是用极坐标表示的平行光束重建公式。然而，积分仍然是关于 ρ 和 θ 的积分。要对 α 和 β 进行积分，需要使用公式 (5-118) 和 (5-119) 进行坐标变换：

(5-123)

$\begin{array}{rl}\displaystyle f (r, \varphi) = \frac{1}{2}\int_{-\alpha}^{2\pi-\alpha}\!\! \int_{\sin^{-1 }(-T/D)}^{\sin^{-1}(T/D)}g(D\sin(\alpha) ,\alpha+ \beta)s(r\cos(\beta+\alpha-\varphi)-D\sin(\alpha))\cos{\varphi}d{\alpha}d{\beta} \end{array}$

其中，我们使用了 $d{\rho}d{\theta} = D\cos{\varphi}d{\alpha}d\beta$ [ 见公式 (5-112) 的解释 ]。

这个公式可以进一步简化。首先，注意到变量 β 的积分限从 -α 到 2π - α 涵盖了 360°的整个范围。由于所有关于 β 的函数都是周期为 2π 的周期函数，因此外层积分的上下限可以分别替换为 0 和 2π 。项 $\sin^{-1}(T/D)$ 存在一个最大值 $\alpha_{m}$ ，对应于 |β| > T 的情况，在此之后 g = 0 (参见图 5.46 )，因此我们可以将内层积分的上下限分别替换为 $-\alpha_{m}$ 和 $\alpha_{m}$ 。最后，考虑图 5.45 中的直线 L( ρ , θ ) 。沿这条直线的扇形束射线积分一定等于沿同一条直线的平行束射线积分。这是因为射线积分是沿直线所有值的总和，因此对于给定的射线，无论使用何种坐标系表示，结果都必须相同。对于 ( α ， β ) 和 ( ρ , θ ) 的对应值，任何射线积分都满足此条件。因此，令 p(α ， β) 表示扇形束投影，则有 p( α ， β ) = g( ρ , θ )，并且根据公式 (5-118) 和 (5-119)， $p(\alpha,\beta) = g(D\sin(\alpha), \alpha + \beta )$ 。将这些观察结果代入公式 (5-123) 即可得到表达式

(5-124) $\begin{array}{rlc} \displaystyle f(r,\varphi) = \frac{1}{2}\int_{0}^{2\pi}{\!\!\!}\int_{- \alpha_m}^{\alpha_m}p(\alpha ,\beta)s \bigg [ r\cos(\beta+\alpha-\varphi)-D\sin(\alpha) \bigg ]D\cos(\alpha)d{\alpha}d{\beta}\end{array}$

这是基于滤波反投影的扇形束重建基本公式。

方程 (5-124) 可以进一步变换为更常见的卷积形式。参考图 5.47，可以证明(参见习题 5.47)

(5-125) $\displaystyle r\cos(\beta + \alpha - \varphi ) = R\sin( \alpha^{'} - \alpha )$

中，R 是从光源到扇形光束中任意一点的距离，而 $\alpha^{'}$ 是这条光束与中心光束之间的夹角。注意， R 和 $\alpha^{'}$ 由 r ，φ 和 β 所确定，将公式 (5-125) 代入公式 (5-124) 得到

(5-126) $\begin{array}{rlc} \displaystyle f(r,\varphi) = \frac{1}{2}\int_{0}^{2\pi}{\!\!\!}\int_{- \alpha_m}^{\alpha_m}p(\alpha ,\beta) \bigg [ R\sin(\alpha^{'}-\alpha)\bigg ]D\cos(\alpha)d{\alpha}d{\beta}\end{array}$

-------------------图 5.46：覆盖兴趣区所需的 α 的最大值。---------------------

-------------------图 5.47：扇形光束射线上任意点的极坐标表示。-----------------

可以证明(参见问题 5.48)：

(5-127) $\displaystyle s( R\sin(\alpha)) = \left [\frac{\alpha}{R\sin(\alpha)}\right ]^{2}s(\alpha )$

利用这个公式，我们可以将公式 (5-126) 写成

(5-128) $\displaystyle f ( r , \varphi) = \frac{1}{2}\int_{0}^{2\pi}\frac{1}{R^{2}} \left [\int_{-{\alpha}_{m}}^{{\alpha}_{m}}q(\alpha,\beta)h(\alpha^{'} -\alpha)d{\alpha}\right ]d\beta$

其中，

(5-129) $\displaystyle h(\alpha) = \frac{1}{2}\left [\frac{\alpha}{\sin(\alpha )}\right ]^{2} s(\alpha )$

以及

(5-130) $q(\alpha, \beta) = p(\alpha,\beta )D\cos(\alpha)$

我们将式 (5-128) 中的内层积分视为卷积表达式，从而表明公式 (5-124) 中的图像重建公式可以实现为函数 q( α ， β ) 和 h(α) 的卷积。与平行投影的重建公式不同，基于扇形束投影的重建包含一个项 $1/R^{2}$ ，这是一个与距离光源的距离成反比的加权因子。公式 (5-128) 的计算细节超出了本文的讨论范围(有关此主题的详细论述，请参见 Kak 和 Slaney [2001])。

与直接应用公式 (5-128) 不同，一种常用的方法(尤其是在软件仿真中)是：(1) 使用公式 (5-118) 和 (5-119) 将扇形束几何结构转换为平行束几何结构；(2) 使用前面开发的平行束重建方法。本节最后将给出一个示例来说明如何实现这一点。如前所述，在角度 β 处获取的扇形束投影 p 对应于在相应角度 θ 处获取的平行束投影 g ，因此，

(5-131) $p( \alpha ,\beta ) = g( \rho , \theta ) = g( D\sin( \alpha ) , \alpha + \beta )$

其中最后一行由公式 (5-118) 和(5-119) 得出。

令 Δβ 表示连续扇形光束投影之间的角度增量，令 Δα 表示光线之间的角度增量，它决定了每次投影中的样本数。我们施加限制

(5-132) ${\Delta}\beta = {\Delta}\alpha = \gamma$

则，对于某两个整数 m 和 n ，有 β = mγ 且 α = nγ ，从而我们可以将公式 (5-131) 写成

(5-133) $p( n\gamma , m\gamma ) = g( D\sin(n\gamma) , (n + m)\gamma)$

该方程表明，第 m 个径向投影中的第 n 条射线等于第 (m + n) 个平行投影中的第 n 条射线。公式 (5-133) 右侧的 D sin(nγ ) 项表明，从扇形束投影转换而来的平行投影并非均匀采样，如果采样间隔 Δα 和 Δβ过大，则可能导致模糊、振铃和混叠伪影，如下例所示。

例子 5.18 ：使用滤波扇形反投影进行图像重建。

图 5.48(a) 显示了以下步骤的结果：(1) 生成矩形图像的扇形投影，其中 Δα = Δβ = 1°；(2) 使用公式 (5-133) 将每条扇形射线转换为相应的平行射线；(3) 使用前面针对平行射线开发的滤波反投影方法。图 5.48(b) 至 (d) 显示了使用 0.5°，0.25° 和 0.125° 的 Δα 和 Δβ 增量时的结果。所有情况下均使用了 Hamming 窗。我们使用不同的角度增量是为了说明欠采样的影响。

图 5.48(a) 的结果清楚地表明，1°的增量太粗糙，因为模糊和振铃现象非常明显。图 5.48(b) 的结果很有意思，因为它与图 5.43(b) 的结果相差甚远，而图 5.43(b) 是我们使用相同的 0.5° 角度增量生成的。事实上，如图 5.48(c) 所示，即使使用 0.25° 的角度增量，重建结果仍然不如图 5.43(b) 的效果好。如图 5.48(d) 所示，我们需要使用大约 0.125° 的角度增量，才能使两个结果变得相似。这种角度增量会产生 180 × (1/0.125) = 1440 个样本的投影，这几乎是示例 5.17 中平行投影使用的 849 条射线数量的两倍。因此，当使用 Δα = 0.125° 时，结果在外观上接近也就不足为奇了。

使用头部模型也获得了类似的结果，只不过在这种情况下，混叠现象以正弦波干扰的形式更加明显。我们在图 5.49(c) 中可以看到，即使 Δα = Δβ = 0.25°，仍然存在明显的失真，尤其是在椭圆的边缘区域。与矩形的情况类似，最终使用 0.125°的增量才获得了与头部模型的反投影图像相当的结果。这些结果解释了现代 CT 系统扇形束几何结构中需要使用数千个探测器来减少混叠伪影的主要原因之一。

--------------------------------图 5.48：从滤波后的扇形反投影重建矩形图像。(a) α 和 β 以 1°为增量。(b) 0.5°为增量。(c) 0.25°为增量。(d) 0.125°为增量。将 (d) 与图 5.43(b) 进行比较。--------------------------------------------------

--------------------------------图 5.49：从滤波扇形反投影重建头部模型图像。(a) α 和 β 步长为 1°。(b) 步长为 0.5°。(c) 步长为 0.25°。(d) 步长为 0.125°。将图 (d) 与图 5.44(b) 进行比较。----------------------------------------------