【转】CT解析重建**

转自：CT解析重建 - 知乎

1、傅里叶变换（Fourier Transform）

白光可以分解成彩色光，彩色光也可合成白光；同样的通过傅里叶变换可将时域下的信号转变成傅里叶域的信号，通过傅里叶逆变换可转换回来。此外，很多问题在傅里叶域讨论会有一片新的天地。 $傅里叶变换F ( w ) = \int _ { - \infty } ^ { \infty } f ( t ) e ^ { - j 2 \pi w t } d t\\ 傅里叶逆变换f ( t ) = \int _ { - \infty } ^ { \infty }F( w ) e ^ { j 2 \pi w t } d w\\ 2D傅里叶变换F ( \mu,v ) = \int _ { - \infty } ^ { \infty } \int _ { - \infty } ^ { \infty } f ( x,y) e ^ { - j 2 \pi (\mu x+vy) } dxdy\\ 2D傅里叶逆变换f(x,y) = \int _ { - \infty } ^ { \infty } \int _ { \infty } ^ { \infty } F ( \mu,v) e ^ { j 2 \pi (\mu x+vy) } d \mu dv\\ 离散2D傅里叶变换F ( \mu , v ) = \sum _ { x = 0 } ^ { m - 1 } \sum _ { y = 0 } ^ { n - 1 } f ( x , y ) e ^ { - j 2 \pi ( \mu x / m + v y / v ) }\\ 离散2D傅里叶逆变换f ( x , y ) = \frac { 1 } { M N } \sum _ { x = 0 } ^ { m - 1 } \sum _ { y = 0 } ^ { n - 1 } F ( \mu , v ) e ^ { j 2 \pi ( \mu x / m + v y / v ) }$

上面一行的图是傅里叶域表示，中心是低频部分越往外表示频率越高，大部分能量都聚集在低频部分，(b)表示将低频部分置0，相当于高通滤波，保留图像的边缘等像素变化大的部分；(c)表示只保留低频成分，低通滤波，图像模糊。

傅里叶变换举例：

2、中心切片定理

简单来说就是 $p(s,\theta)的1D \ \ \ Fourier \ \ transform \ \ P_{\theta}(w)=\\ F(\mu,v)\ \ when\ \ \mu=wccos\theta,v=wsin\theta$

证明过程：

$\left. { P _ { \theta } ( w ) = \int _ { - \infty } ^ { \infty } f _ { \theta } ( s ) e ^ { - j 2 \pi \omega s } d s } 傅里叶变换\\ { = \int _ { - \infty } ^ { \infty } [ \int _ { - \infty } ^ { \infty } \int _ { - \infty } ^ { \infty } f ( x , y ) \delta ( x \cos \theta + y \sin \theta - s ) d x d y ] }e ^ { - j 2 \pi \omega s } d s拉东变换\\ { = \int _ { - \infty } ^ { \infty } \int _ { - \infty } ^ { \infty } [ \int _ { - \infty } ^ { \infty } f ( x , y ) \delta ( x \cos \theta + y \sin \theta - s ) e ^ { - j 2 \pi w s } d s ] d x d y } \right.积分换序\\ =\left. { \int _ { - \infty } ^ { \infty } \int _ { - \infty } ^ { \infty } f ( x , y ) e ^ { - j 2 \pi w } ( x \cos \theta + y \sin \theta ) d x d y } \\ = { \int _ { - \infty } ^ { + \infty } \int _ { - \infty } ^ { \infty } f ( x , y ) e ^ { - j 2 \pi ( w x \cos \theta + w y \sin \theta ) } d x d y } \\ ={ F ( w \cos \theta , w \sin \theta ) } \\ =\ { F ( \mu , v ) } \right|_{\mu=wcos\theta,v=wsin\theta}.$

3、一些重建方法

方法1：FBP

（1）求投影数据 p(s, θ) 的以 s 为变量的一维傅里叶变换，得到P(ω,θ ) 。

（2）对 P(ω,θ ) 乘以斜坡滤波器的传递函数 |ω|，得到 Q(ω, θ)。

（3）求 Q(ω, θ) 的以ω 为变量的一维傅里叶反变换，得到 q(s, θ)。

（4）反投影得到重建图像f(x,y)。

方法2：根据傅里叶变换理论，在 ω 域中做乘法等价于在 s 域中做卷积

（1）q(s,θ ) = p(s,θ ) ∗ h(s)；h(s) 是卷积积分中的卷积核，是H(ω)=|ω|的一维傅里叶反变换

（2）反投影得到重建图像f(x,y)。

方法3： $H(w)=|w|=i2\pi w \times \frac{1}{i2 \pi} \ sgn(w)$

傅里叶变换的两个性质:

性质 1：在傅里叶域 (即ω 域) 中乘以 i2πω 相当于在空间域 (即 s 域) 中求导数。

性质 2：函数 -i sgn(ω) 的傅里叶反变换是 1/(πs)。与 1/(πs) 做卷积叫做希尔伯特变换。

$q(s,\theta) = \frac{dp(s,\theta)}{ds} \ast \frac{-1}{2\pi^2s}$ 然后再反投影。

方法4：改变斜坡滤波和反投影的次序，先反投影后滤波。

（1）对反投影得到的图像 b(x, y) 求二维傅里叶变换，得到B(ω x ,ω y ) 。

（2）对 B(ω x ,ω y ) 乘以斜坡滤波器的传递函数 |ω|= ωx2 + ω y2 ，得到F(ω x ,ω y ) 。

（3）对 F(ω x ,ω y ) 求二维傅里叶反变换，得到 f (x, y) 。

方法5：求导，希尔伯特变换，和反投影可换序

（1）对投影数据 p(s,θ ) 以变量s求导(实际上是求偏导)，得到dp(s,θ ) / ds 。

（2）对 dp(s,θ ) / ds 做 180° 的反投影。

（3）对反投影得到的图像逐行的做(一维的)希尔伯特变换。其方向是与探测器在 90º 角的位置相平行。

希尔伯特变换可以在空间域中做卷积来实现，也可以在傅里叶域中做乘法来实现。除此以外，希尔伯特变换还可以在空间域中做积分来实现，这个积分并非卷积，而是在有限区间上的积分。这个有限积分的希尔伯特变换在处理不完整的(即截断的)投影数据时有着重要的应用。

4、卷积核

参考：毛小渊. 二维CT图像重建算法研究[D].南昌航空大学,2016.

上面介绍的滤波器H(w)=|w|是一个频带无限地滤波器，无法实现，所以考虑其替代。在实践滤波过程中，可以把一个信号的绝大部分用有用频率予以保留，丢弃无关紧要的频率，在实际的卷积过程中，投影数据的傅立叶变换是有限带宽的。也就是说在频率间隔(B,B)以外的能量0。可得： $q(s,\theta)= \int_{-B}^{B}|w|P _ { \theta } ( w )e^{j2\pi ws_{r}}dw$ 根据奈奎斯特采样定理，为了保证无混叠的采样，采样间隔必须不大于最高截止频率 2 倍的倒数，也就是： $B=1/2d$ 。

（1）R-L滤波器： $H_{R-L} (w) = |w|W(w) = |w|rect(w/2B)\\ 当|w|<B时，rect(w/2B)=1，其他情况为0\\ 其中，w为空域频率，W( \ )表示窗函数\\ 相应的时域卷积函数为：h_{R-L}(s_{r}) = 2B^2sinc(2s_{r}B)-B^2sinc^2(s_{r}B)\\ 将s_{r}=nd,B=1/2d代入得到离散形式：\\ h_{R-L}(nd) = 1/4d^2,\ n=0\\ h_{R-L}(nd) = 0,\ n=偶数\\ h_{R-L}(nd) = -1/n^2\pi ^2d^2,\ n=奇数\\$

R-L 滤波器的频域波形如图所示，其中截止频率 d=1。它在频域中的图像类似于斜坡，故也称为斜坡滤波器。R-L 滤波器形式简单实用，用它重建图像，轮廓清楚。缺点是有 Gibbs 现象，表现为明显的振荡响应。

（2）S-L滤波器： $H_{S-L}(w) = |w|sinc(w/2B)rect(w/2B)=|\frac{2B}{\pi}sin\frac{\pi w}{2B}|\\ h_{S-L}(s_{r}) = \frac{1}{2}(\frac{4B}{\pi})^2\frac{1-4Bs_{r}sin(\frac{\pi}{2}4Bs_{r})}{1-(4Bs_{r})^2}\\ 离散形式：h_{S-L}(nd) =\frac{-2}{\pi ^2 d^2(4n^2-1)}\\$

（3）Cosine滤波器：

窗函数： $H=cos(\frac{\pi w}{2B})rect(\frac{w}{2B})$

（4）Hanning 滤波器与 Hamming 滤波器

广义Hanning窗口： $H=\alpha+(1-\alpha)cos(\frac{\pi w}{B})rect(\frac{w}{2B})\\$

α为参数，[0.5,1),当α=0.54时，为 Hamming 窗函数： $H_{ham}=0.54+0.46cos(\frac{\pi w}{B})rect(\frac{w}{2B})\\$

5、扇束

6、锥束