CT成像技术—20231205

平行束傅里叶切片定理（谢强，3.5）

静坐标系$x,y$；旋转坐标系$s,t$，$s$轴与射线同方向，$t$轴平行于探测器，$t$轴与$x$轴夹角为$\theta$。
$t=(x,y)\cdot (cos\theta ,sin\theta )$
$s=(-x,y)\cdot (sin\theta ,cos\theta )$
注：我写博客主要是文字，需要大家和书籍上的图片对应。

$f(x,y)$表示物体截面（静坐标系描述），其在旋转坐标系中的描述为$f^{\prime }(t,s)$
$p(t,\theta )$表示$f(x,y)$沿s方向的投影
$f(x,y)$的二维傅里叶变换：$F(u,v)$
$F(u,v)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-i2\pi w(xu+yv)}dxdy$

$p(t,\theta )$对变量$t$的一维傅里叶变换：$P(w,\theta )$
$p(t,\theta )={\int }_{-\infty }^{\infty }{f}^{\prime }(t,s)ds$
$P(w,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{f}^{\prime }(t,s)ds{e}^{-i2\pi wt}dt$
将$P(w,\theta )$的描述从旋转坐标系转化至静坐标系
$P(w,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)e^{-i2\pi w(xcos\theta +ysin\theta )}dxdy$
注：最后一步还用到了两个坐标系中面积元素的关系$dsdt=Jdxdy$，$J$为两行两列的行列式，其元素为$s,t$对$x,y$的偏微分

令$u=wcos\theta ,v=wsin\theta$，便可将$F(u,v)$与$P(w,\theta )$联系起来。（中心切片定理）
$F(wcos\theta ,wsin\theta )=P(w,\theta )$
$u,v$为频域中的静坐标系；$w,\theta$为频域中的极坐标系

依靠中心切片定理直接重建需要：
1、将投影变换至频域
2、将频域投影重新插值，从极坐标系变换至直角坐标系
3、进行二维傅里叶反变换
第二步中产生的局部误差，经过第三步反变换后会变成全局误差，这就是频域插值的敏感性。当然现在肯定有学者在研究如何克服这种问题。

另外，目标重建（感兴趣区域重建，ROI重建）时的图像矩阵大小与ROI的尺寸成反比。对于很小的ROI，则傅里叶变换困难。
注：ROI感兴趣区域，FOV视野。之后有时间的话，我会研究如何写代码进行ROI重建。

平行束滤波反投影算法（谢强，3.6）

如前所述，依据中心切片定理直接重建很困难，本节将其进一步转化为滤波反投影算法。

空域截面可以由其频域表述重建：
$f(x,y)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }F(u,v)e^{i2\pi w(xu+yv)}dxdy$
将$F(u,v)$的表述由直角坐标系变换至极坐标系$w,\theta$：
$f(x,y)={\int }_0^{2\pi }d\theta {\int }_0^{\infty }F(wcos\theta ,wsin\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}wdw$
由此我们可以将$P(w,\theta )$替换$F(wcos\theta ,wsin\theta )$，这便是用到了中心切片定理：
$f(x,y)={\int }_0^{2\pi }d\theta {\int }_0^{\infty }P(w,\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}wdw$
$f(x,y)={\int }_0^{\pi }d\theta {\int }_0^{\infty }P(w,\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}wdw+{\int }_0^{\pi }d\theta {\int }_0^{\infty }P(w,\theta +\pi )e^{-i2\pi w(xwcos\theta +ysin\theta )}wdw$
另外，由于$P(w,\theta +\pi )=P(-w,\theta )$，将其带入上式可得：
$f(x,y)={\int }_0^{\pi }d\theta {\int }_{-\infty }^{\infty }P(w,\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}|w|dw$
注：我现在距离学习高数已有8年，真是忘记的太多了。

这便是基于傅里叶中心切片定理进一步推断出的先滤波再反投影算法。

将滤波反投影公式中的后一部分记为$g(x,y)$，可将公式重写为：
$f(x,y)={\int }_0^{\pi }g(x,y)d\theta$
$g(x,y)$表示将探测器数据经过空域滤波（或频域滤波再反变换回空域）后，将数据沿射线方向进行涂抹。
对每个角度下的数据重复操作，直到180度为止。最后得到的便是精确的重建截面。
注：这个滤波器称为斜坡滤波器

滤波在频域中就是与一个滤波函数相乘，在空域中就是与一个滤波函数卷积。
频域斜坡滤波函数$|w|$
空域斜坡滤波函数$\xi (t)={\int }_{-\infty }^{\infty }|w|e^{i2\pi wt}dw$

滤波函数无法实现，在空域中：
$\xi (0)={\int }_{-\infty }^{\infty }|w|dw$，当$w\to \infty$时，$\xi (0)\to \infty$。
在频域中：
当$w\to \infty$时，$|w|\to \infty$。

我们假定$P(w,\theta )$是有限带宽的，即在$[-\Gamma ,\Gamma ]$以外$P(w,\theta )=0$。
$g(t,\theta )={\int }_{-\Gamma }^{\Gamma }P(w,\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}|w|dw$
注：截面，截面投影的傅里叶变换都只是一部分区域内值不为0，超出这部分区域外值为0。因此，我们只需要知道这部分区域下的滤波函数即可。我们认为物体具有有限空间紧支集。

为保证无混叠采样，$\Gamma =1/(2\delta )$，$\delta$为探测器两个像素中心的间距。
$\Gamma$单位$cycles/mm$，$\delta$单位$mm$。
注：空间频率指在每单位长度上，出现几个同样的几何结构。空间频率单位：$cycles/mm$或$LP/mm$。

有限带宽的空域滤波函数$h(t)={\int }_{-\Gamma }^{\Gamma }|w|e^{i2\pi wt}dw$
$h(t)=\frac{1}{2{\delta }^2}\frac{sin2\pi \Gamma t}{2\pi \Gamma t}-\frac{1}{4{\delta }^2}(\frac{sin\pi \Gamma t}{\pi \Gamma t}{)}^2$
$h(t)$及其傅里叶变换都是一个实偶函数。
我们只关心位于像素中心的滤波器函数值，即当$t=n\delta$时（$n$为整数）的函数值：
$h(n\delta )=\frac{1}{4{\delta }^2}$，$n=0$
$h(n\delta )=0$，$n$为偶数
$h(n\delta )=-\frac{1}{(n\pi \delta {)}^2}$，$n$为奇数

截面在静坐标系中表示$f(x,y)={\int }_0^{\pi }d\theta {\int }_{-\infty }^{\infty }P(w,\theta )e^{i2\pi w(xwcos\theta +ysin\theta )}|w|dw$
在旋转坐标系中表示$f^{\prime }(t,s)={\int }_0^{\pi }d\theta {\int }_{-\infty }^{\infty }P(w,\theta )e^{i2\pi wt}|w|dw$
$f^{\prime }(t,s)=f(x,y)$

对$P(w,\theta )$滤波后的傅里叶反变换，相当于直接在空域中对$p(t,\theta )$进行卷积滤波。
$f^{\prime }(t,s)={\int }_0^{\pi }d\theta {\int }_{-t_m}^{t_m}p(t^{\prime },\theta )h(t-t^{\prime })d{t}^{\prime }$
注：在$[-t_m,t_m]$之外$p(t^{\prime },\theta )=0$
对$g(t,\theta )$离散表述：$g(n\delta ,\theta )=\delta {\Sigma }_{k=0}^{N-1}h(n\delta -k\delta )p(k\delta ,\theta ),n=0,...,N-1$
$N$表示一维探测器的像素个数，当像素个数较多时，还是在频域滤波的速度较快。

在$[-\Gamma ,\Gamma ]$内，有限带宽滤波函数$h(t)$和$|w|$有些差异，主要体现在直流部分（0频部分）。

若要在频域中滤波，会导致warp-around卷绕效应（周期间干涉）从而产生干涉伪像。为了避免干涉伪像，需要对空域每个投影数据补$N-1$个0，之后将数据变换至频域进行滤波再反变换回空域。
注：时域补零相当于频域插值，增大了频域分辨率。继续深究，需要看信号处理书籍了。

斜坡滤波器突出高频成分，斜坡滤波可当作一个反卷积运算（微分运算），去掉反投影产生的模糊。

$h(t)$相当于在$|w|$上加了一个矩形窗，我们也可以加一个Hanning窗，sinc窗等。
通过加窗实现在空间分辨率与噪声之间的折中，矩形窗就是空间分辨率最好但噪声最大，通过如sinc函数可以降低空间分辨率但也降低噪声。
另外，除了改变窗形状，我们也可以修改窗的截止频率，以实现空间分辨率与噪声之间的折中。
针对软组织，肺，骨的滤波器是不一样的，这是CT商业公司的机密。
注：用其他窗，将截至频率提前，这两种做法都是抑制高频，可去除噪声，但将边缘变模糊了。

前面一直说的是滤波，现在说反投影，即如何将滤波后的数据进行反投影。
反投影两种方式：射线驱动，截面驱动。
射线驱动：计算出一条过探测器像素中心且与探测器垂直的直线，由该直线与截面的相交情况，将该滤波后的探测器像素值分给截面像素。
截面驱动：计算出一条过截面像素中心且与探测器垂直的直线，确定出该直线与探测器的交点位置。该点的值可通过滤波后的探测器像素插值出来，其值便作为截面像素值。（这是我自己想到的，书上写的内容我刚开始没读懂，现在看跟我想的一样，读书得静下心来。）
前者需要二维插值，后者需要一维插值，计算量小。
对于第二种方法，为加速反投影，还可以先对滤波后的数据进行高阶插值，再采用第二种方法求解截面像素值，此时就仅采用低阶插值。

好的插值方法会产生一条平坦的MTF曲线，这样系统响应独立于输入频率，空间分辨率更好。
若一种插值方法的所产生的截面空间分辨率高，其噪声肯定也增大了。这就是空间分辨率与噪声的折中，不一定某种插值方法是最好的，只选择适合自己任务的。

滤波反投影算法的完整流程：
1、对投影数据补零
2、将投影数据进行傅里叶变换
3、对频域投影数据进行滤波
4、对频域数据进行傅里叶反变换
5、对滤波投影数据进行高阶插值
6、用高阶插值数据，进行射线驱动的反投影，确定截面像素值时用低阶插值
在所有角度下，重复上述操作。

ROI重建

$50cm$的FOV，图像矩阵$512\ast 512$，则图像像素尺寸约为$50cm/512\approx 1mm$，图像最高频率为$1/(2\ast 1)LP/mm=5LP/cm$。如果希望最高频率为$20LP/cm$，可采用：
方法1，其余不变，只将图像矩阵设为$2048\times 2048$。$50cm$的FOV，每个图像像素尺寸约为$0.25mm$，则图像最高频率$1/(2\ast 0.25)LP/mm=20LP/cm$。
方法2，其余不变，只将FOV设为$12.5cm$。图像矩阵$512\ast 512$，每个图像像素尺寸约为$0.25mm$，则图像最高频率$1/(2\ast 0.25)LP/mm=20LP/cm$。

以上例子说明想检查具有更高空间分辨率成分的物体横截面，要么增大重建图像矩阵尺寸或减小重建视野。

ROI重建与直接对原图像进行插值放大不同，前者的空间分辨率更大。

前面是今天写的平行束重建内容，后面是12月5日写的。
接下来我会再写关于扇形束与锥束重建的博客。

历史，方法

成像就是物质波与被检测物质相互作用结果的可视化。

1895年，德国，伦琴发现X射线；1901年，伦琴获得诺贝尔物理学奖。
1917年，奥地利，Radon给出图像投影，以及由投影重建图像的数学理论。可惜并未受到重视。
1963年，美国，Cormack给出当代图像重建的数学方法。
1972年，英国EMI员工，Hounsfield研制出第一台临床用的CT。
1979年，Cormack与Hounsfield获得诺贝尔生理和医学奖。

图像重建方法分为两类：以Radon变换为基础的解析方法，以解方程为思想的迭代方法。

1、解析方法
二维图像重建基础是傅里叶中心切片定理，基于该定理可以得到两种方法：直接傅里叶重建算法、滤波反投影算法。对于扇形束需经过重拍或加权修正，便可进行滤波反投影重建。
三维图像重建算法分为近似与精确算法。前者包括：FDK（1984，Feldkamp），1998年Turbell提出螺旋锥束重建的PI方法。后者包括：Grangeat算法（以1961年Kirillov给出的锥束几何复值函数的逆变换公式为基础，Tuy，Grangeat，Smith分别提出三种算法，其中Grangeat算法易于实现，从而最出名），Katsevich算法，BPF算法。
另外，Tuy给出三维图像精确重建的条件。

2、迭代方法
代数迭代：1970年，Gordon，ART算法；1984年，Andersen，SART算法。
统计迭代：EM，最小范数法，MAP。

Python包

python包列表 https://tomopedia.github.io/software/，列举：
ASTRA
TomoPy
TIGRE
NeuTomPy 中子断层重建 https://neutompy-toolbox.readthedocs.io/en/latest/
PYRO-NN 深度学习重建 https://github.com/csyben/PYRO-NN

另外，可以参考Github学习资源：
扇束，锥束CT https://github.com/leehoy/CTReconstruction
一些论文 https://github.com/LoraLinH/Awesome-CT-Reconstruction
锥束CT https://github.com/Apoorva0607/CBCT
…

三维数据（曾更生，医学图像重建，5.1，5.2，5.3）

平行的线积分数据

三维线积分投影成像的中心切片定理：（将二维线积分投影成像的中心切片定理的表述进行升维便是三维的。）
$f(x,y,z)$的在二维探测器上投影的二维傅里叶变换等于$f(x,y,z)$的三维傅里叶变换被过原点且与探测器平行的平面所截取的片段

图像重建算法取决于方向向量$\theta$的轨迹不同而不同。（$\theta$表征射线方向，是一个单位向量）

Orlov条件（Orlov’s condition）：如果圆心在原点的所有单位圆都与$\theta$的轨迹有交点，则获得的就是一组完整的投影数据。

平行的面积分数据

Radon变换：对$n$维物体$f$做平行的$n-1$维超平面积分。对二维物体做平行的线积分，对三维物体做平行的面积分。

三维面积分投影成像的中心切片定理：$f(x,y,z)$的在一维探测器上的面积分投影的一维傅里叶变换等于$f(x,y,z)$的三维傅里叶变换被过原点且与一维探测器平行的直线所截取的片段。
即$f(x,y,z)$的在一维探测器上的面积分投影的一维傅里叶变换等于$f(x,y,z)$的三维傅里叶变换在一维探测器上的面积分投影。

锥形束数据

锥形束数据目前并没有中心切片定理，但1983年，Tuy给出的引理近似于锥形束中心切片定理。

Tuy给出锥形束数据充分条件：每一个与物体相交的平面都必须包含至少一个锥形束的焦点位置。
（扇形束数据充分条件：每一条与物体相交的直线都必须包含至少一个扇形束的焦点位置。）

螺旋轨道，圆圈加直线轨道满足Tuy条件；圆轨道不满足Tuy条件。