由投影重建图像：滤波反投影、FDK、TFDK三维重建算法理论基础

1. 基础理论从：
[1] RafaelC.Gonzalez, RichardE.Woods, Gonzalez,等. 阮秋琦等译.数字图像处理(第三版)[M]. 电子工业出版社, 2011.P232

[2] RafaelC.Gonzalez, RichardE.Woods, StevenL.Eddins. 阮秋琦译.数字图像处理:MATLAB版:本科教学版[M]. 电子工业出版社, 2014. (第二版)P102

参考文献[1]从P228-P245

了解：

1.    由投影重建图像

1.     计算机断层（CT）的原理

2.     投影和Radon（雷登）变换（radontransform）

3.     正弦图（sinogram）和Shepp-Logan幻影

4.     傅里叶切片定理

5.     使用平行射线束滤波反投影的重建

6.     Ram-Lak滤波器与Hamming、Hann等窗函数

7.     使用扇形射线束滤波反投影的重建

1. 由投影重建图像

利用Matlab来复现参考文献[1]中的处理流程：

首先生成一幅原始图像：

%代码如下：
I=zeros(512,512);
for i=1:512for j=1:512if ((i-256)*(i-256)+(j-256)*(j-256))<1024I(i,j)=1;endend
end
imshow(I,[]);

得到1°和90°两个方向上的投影（一维），并回抹得到反投影图像（二维）

         

这两个方向的反投影相叠加：

%代码如下：（这里参考了MATLAB自带的函数说明输入：doc iradon打开该参考说明）
%就是怎么得到特定角度下的一个投影，并形成反投影图像：
R=radon(I,0:179);
r1=R(:,2);
II1=iradon([r1 r1],[1 1])/2; %得到1°角度下的投影图像
imshow(II1);r90=R(:,91);
II90=iradon([r90 r90],[90 90])/2; %得到90°角度下的投影图像
imshow(II90);II1_90=II1+II90;%二者叠加
imshow(II1_90);

同理，得到45°和135°方向的反投影图像，并将这四幅反投影图像相叠加：

%代码如下：
r45=R(:,46);
II45=iradon([r45 r45],[45 45])/2;
imshow(II45);r135=R(:,136);
II135=iradon([r135 r135],[135 135])/2;
imshow(II135);II_1_90_45_135=II1+II45+II90+II135;%四个反投影图像相叠加！
imshow(II_1_90_45_135);

当反投影的角度采样增多时，这里从0°到179°每隔1°共采180个反投影，相叠加后形成反投影图像：（不使用滤波器时的反投影重建图像如下：）

%这里用了iradon函数，线性插值，不使用滤波器。
ii2=iradon(R,0:179,'linear','none');
imshow(ii1,[]);

ii1=iradon(R,0:179,'linear','Ram-Lak');%默认使用线性插值，Ram-Lak滤波器！

该结果如下：

将重建后的图像与原图像进行对比：（左边是原图、右边是重建得到的图像，注意模糊与振铃现象）

  

我们再使用MATLAB自带的一头部幻影图像Shepp-Logan:

原图：

未使用滤波器的重建图像：

使用Ram-Lak滤波器  ：

使用Shepp-Logan滤波器：

%生成以上图像的代码如下：
P=phantom(512);
theta=0:179;
[R,xp]=radon(P,theta);
I1P=iradon(R,0:179,'linear','none');
I2P=iradon(R,0:179,'linear','Ram-Lak');
I3P=iradon(R,0:179,'linear','Hamming');

以下内容参考文献[1]将滤波反投影的基础知识过一遍：

笛卡尔坐标系的一条直线由它的斜截式描述： ；或由其法线方程来描述：

平行射线束的投影可以由一组这样的直线建模。如图5.37所示，投影信号中的任意一点由沿着直线 的射线和给出。连续变量的情况下，线求个变为线积分，由下式给出：

该式是沿xy平面内任意一条直线的f(x,y)的投影（线积分）的公式，就是雷登(Radon)变换。

符号R{f(x,y)}或R{f}有时用于代替(3)式中的 来表示f的雷登变换。雷登变换是由投影重建图像的基石，计算机断层(CT)是其在图像处理领域的主要应用。在离散情况下，（3）式变为

正弦图与雷登变换：

    

正弦图包含了重建图像f(x,y)所需的信息。

正弦图的视觉分析仅限于实际应用、但有时对于算法开发是有帮助的。

 

CT的关键目的是从投影得到物体的三维表示。其方法是反投影每一个投影，然后对反投影求和以产生一幅图像（切片），再堆积所有的结果以产生三维物体的再现。（这里是指用二维切片堆积成三维体，与FDK/TFDK直接就是三维重建是不同的！）

5.11.4傅里叶切片定理

 

傅里叶切片定理即投影的一维傅里叶变换和被投影区域图像的二维傅里叶变换间的关系。

投影 的一维傅里叶变换为：

式（11）就是著名的傅里叶切片定理（或投影切片定理）。它说明了一个投影(一维)的傅里叶是得到这个投影的二维区域f(x,y)的二维傅里叶变换对应角度下的一个切片。正如图5.41所示，任意一个投影的一维傅里叶变换可以沿着一个角度提取一条直线的F(u,v)的值来得到，而该角度就是投影时所用的角度。

下面推导滤波反投影公式，将用到傅里叶切片定理。

F(u,v)的反傅里叶变换为：

当c=0.54时，该函数称为汉明窗（RichardHamming）；

当c=0.5时，称为韩窗（Juliusvon Hann）

加了窗函数的滤波器在空域的振铃现象减弱。

我们可以预期，由于使用汉明窗的反投影有较小的振铃，但稍微模糊一点。见下图：

原图  ；使用Ram-Lak滤波器  ；  使用Hamming窗加窗后的滤波器

        

在CT的多数应用中（特别是医学上），像振铃这样的人为缺陷有严重的厉害关系，使其最小化是有意义的工作。调整滤波算法（可以做文章的地方，一些硕士论文就这样自创新的滤波器，发现效果有所改进，好，成文！比如文献

[1] 张銮. 基于平板探测器的锥束CT重建技术研究[D]. 中北大学, 2010.）、和硬件制造方面的改进（如，提供探测器的检测细腻度，即采样粒度）

因为斜坡滤波器（甚至在被加窗时）在频率域的直流项为零，故每一幅反投影图像的均值将为零。这将意味着，每一幅反投影图像都将有正像素和负像素值，当所有的反投影图像相加形成最终的重建图像时，一些负像素值位置可能变成正像素，而平均值可能不为零，但是，典型地，最终的图像将还是有负像素值。

当有关平均值的知识未知时，就使用标定的方法将图像的像素值都归一化到一个区间[0,255]。当典型的平均值的知识是可用的时候，可将该值加到频率域的滤波器上，从而抵消斜波并防止直流项为零。当在空域中使用卷积时，截断空间滤波器的长度（斜坡的反傅里叶变换）的真正效果都将防止其有零均值，这样就完全避免的迫零问题。

滤波反投影算法：

1.     计算每一个投影的一维傅里叶变换；

2.     用滤波函数 乘以每一个傅里叶变换；

3.     得到每一个滤波后的变换的一维反傅里叶变换；

4.     对步骤3得到的所以一维反变换积分（求和）