CT图像重构方法详解——傅里叶逆变换法、直接反投影法、滤波反投影法

Good

绪

      在做CT图像处理的时候遇到很多问题，对于滤波反变换有许多细节存在疑问，经过多天查找资料和利用MATLAB程序一步步实现后终于豁然开朗，于是想要总结成文，作为笔记方便今后查看。文中若有错误欢迎指出！

目 录

（1）傅里叶逆变换法
（2）直接反投影法
（3）滤波反投影法
（4）MATLAB程序

---

      CT图像的形成和重构，在数学上的描述分别为拉东变换和拉东反变换。简单点说，拉动变换是用于描述通过X射线扫描物体形成CT图像这个过程的一种数学表示，而拉东反变换描述的则是对CT投影数据进行重构，还原成物体图像的一种数学方法。本文章主要讲的是CT图像重构的方法，也就是拉东反变换的具体实现方法，当然可能也有小伙伴想了解一下什么是拉东变换，这里贴上传送门：拉东变换。
其中拉东反变换在数学上的表达式如下：

由公式可看出拉东逆变换在数学上实际上是一个二维傅里叶逆变换，但由于二维傅里叶逆变换计算量大，耗时长，因此计算机对于拉东逆变换的具体实现方法主要为反投影法，反投影法又分为直接反投影法和滤波反投影法。由此可见CT重构方法主要为以下三种：
（1）傅里叶逆变换法
（2）直接反投影法
（3）滤波反投影法
下面将具体对这三种图像重构方法进行解释。

---

（1）傅里叶逆变换法

该方法基于一个重要的定理：中心切片定理。
该定理简单地理解就是：通过角度为θ扫描得到的投影，该投影的一维傅里叶变换，与对整个图像二维傅里叶变换后，二维频域中对应θ角度的一个切片信号是相同的，下面两个图理解起来更直观。
     

     图1        

     图2    根据该理论，傅里叶逆变换法可以简单分成以下步骤：
① 假设每旋转1°就扫描一次，当对物体扫描了180°之后，我们就能得到180个投影信号（就是180根投影线）→在临床上，若使用平行扫描CT，我们拿到手的数据就是这个（在数学上，就是对图像进行拉东变换）
② 对180个投影信号进行一维傅里叶变换
③ 对②得到的180个一维频域信号，根据相对应的扫描角度，在空间中旋转排列，拼成一个二维频域空间（如图3）
④ 由于数据是离散的，直接按照角度进行排列难以铺满整个二维空间，因此还需要对空缺的地方进行插值（一般三次样条插值效果最好），但插值会带来一定误差。另外，由于中心的信号密集，周围的信号稀疏，显然会损失一部分高频数据，造成高频信号失真，这就是采用傅里叶逆变换法重构图像时会使得图像边缘模糊的原因。
⑤ 对④中拼接而成的二维图像进行二维傅里叶逆变换，就可重构原图
     

     图3  傅里叶逆变换法    该方法的缺点有：
a/高频信号有所失真
b/在插值时还涉及到极坐标和直角坐标的变换，计算量大
c/需要用到二维傅里叶逆变换，总体耗费时间长
由于计算机处理二维傅里叶逆变换的计算量太大，因此很少直接使用该方法实现拉东逆变换。

（2）直接反投影法

反投影法的原理是将所测得的投影值，按照其原投影路径，平均地分配到经过的每一个点上，把各个方向的投影值都这样反投影后，在把每个角度的反投影图像进行累加，从而推断出原图。为方便理解，下面分别用两张图，通过MATLAB编程来还原一下该过程（程序在最后）：
原图如下：（图1简单，图2复杂）

①首先看一下图像经过拉东变换（即CT扫描）后的数据（图中每一列代表对应每一个角度的投影信号，在这里只是拼起来成一张图了）

②对每个角度的信号进行反投影（选了几个代表性角度）：

② 把以上角度发反投影图像叠加起来：

当越来越多的反投影图像加起来之后，可以看到：

到此，直接反投影法结束。值得注意的是，由于反投影图是离散叠加的，显然在中心处信号集中，边缘处信号稀疏，因此在最后需要在空缺的地方进行插值，才能得到最终的原图像。上面演示的过程中没有进行插值，仅仅是把反投影图叠加，因此能够看出每一张图都不是很光滑。
从以上过程可以很直观地看出直接反投影法出现伪迹的原因：
①在“回抹“过程中（就是把投影信号平均到每个二维空间点的过程中），会把原图像本来是0的像素点也”抹“上一个平均值，最终使得重构的图像中存在误差；
②插值过程中会带来误差，且周围信号稀疏，高频信号有所失真，导致图像边缘模糊；
③在反投影图不断叠加过程中，能够看到这种叠加方式会带来明显的“星状伪迹“，这是造成图像边缘模糊的最重要因素。在投影数据少的时候更明显（现实中投影数据的多少取决于机器）。
关于第③点，在数学上有研究表明，反投影重构后的图像fb(x,y)，和真实原图f(x,y)之间存在以下卷积关系：

如果要重构成真实的图像，就需要把这个由于反投影法本身带来的1/r效应去掉，由于1/r会使得图片变得模糊，因此1/r也称之为模糊因子。在实际操作中，把1/r影响去除的方法这就是下面这种最常用的（计算机、商业普遍使用）CT图像重构方法：滤波反投影法。

（3）滤波反投影法

从上面的式（1）可以看到，1/r跟原图像是卷积关系，通过傅里叶变换转变到频域上后会变成乘积关系，这说明直接在频率上处理会更加简便。

在频域中，式（3）中的r称之为权重因子，其实就相当于一个滤波器，想要去除模糊因子的影响，重构成原始图像，只需要三步：
①把重构图像fb(x,y)进行二维傅里叶变换得到Fb(x,y)
②在频域中把Fb(x,y)与滤波器r相乘
③把r×Fb(x,y)进行二维傅里叶逆变换即能够得到原图f(x,y)
具体过程见下图，图中的ρ就是公式中的r。
      
     

     图4  直接反投影法   
这种先反投影，后滤波的方法需要用到两次二维傅里叶变换，计算量太大，因此也不是一个特别理想的矫正方法。
而滤波反投影法，顾名思义就是先滤波，后反投影的重构方法。具体步骤如下：
①对180个投影信号先分别进行一维傅里叶变换
②在频域中对所有投影信号进行滤波，也就是乘上权重因子r
③把所有滤波后的信号再进行一维傅里叶逆变换，还原到时域
④对每一个已经滤波的投影信号进行反投影，最后叠加（这一步跟前面的直接反投影法步骤完全相同，只是投影信号已经过了滤波）
具体流程见下图：
     

     图5  滤波反投影法    这个方法的好处是，把两次二维傅里叶变换变成了两次一维傅里叶变换，计算速度大大提升。于是滤波反投影法的核心问题就变成了如何选择一个合适的滤波器r。数学中的滤波器r是一个理想化的滤波函数，现实中不存在，因此只能设计一个近似的滤波函数来代替r的作用。
下面常用的滤波函数有：R-L滤波函数和S-L滤波函数。
①R-L滤波函数（Ramp-Lak滤波器）

从频域图中可以看到该滤波器的作用是把低频信号减少，从而突出高频信号，因此经过该函数滤波后，重建图像的轮廓会更清晰，并且函数简单，实现起来更方便。但由于该函数在频域上有一个加窗处理（就是把高频部分直接截断了），因此在时域中会出现有许多振荡的小尾巴（截窗振荡效应），这将会让重构出来的图像存在Gibb’s现象（重构图像存在振荡，不连续）。
②S-L滤波函数（Shepp-Logan滤波器）

S-L滤波器在频域上不是直接加窗截断，而是通过一些比较平滑的窗函数对函数进行约束。因此经过S-L滤波后重构的图像振荡更少，重构的图像质量比R-L滤波器更好一些，但是因为S-L在高频段并没有直接截断，偏离了理想滤波器的效果，因此在高频段（轮廓）上的重构效果没有R-L滤波器好。
可以在下图更直观地看到两个滤波器之间的区别：

下面通过MATLAB编程来还原一下滤波反投影法的过程（程序在最后，只用复杂那个图进行展示）：
①对每一列投影信号分别在频域进行R-L滤波(高频信号明显突出了)

②对每个滤波后的投影信号反投影（这里开始与直接反投影法一样）

③ 把所有反投影图像叠加起来 

为了验证手动重构的效果，下面将使用matlab中的iradon函数直接对CT投影数据进行重构，把手动重构的图像和matlab自带函数重构的图像进行对比。
（iradon是matlab中的拉东逆变换函数，但是其实际采用的重构方法并不是对数据进行二维傅里叶逆变换，而是通过滤波反投影法实现的，且默认插值方法为“线性插值linear“，默认滤波方法为”R-L滤波函数“，这在matlab的help中可以直接看到。）

仔细观察两个重构结果，会发现iradon的重构结果更光滑一些，那是因为手动方法只是简单粗暴地把全部离散的反变换图片叠加起来，而iradon则会在最后对图像进行线性插值，使得图像更连续更光滑，重构效果更好。
最后把iradon重构的图像放大来看看图像细节：

可以看到重构的图像其实并不连续，甚至能够看出一些振荡的波纹，这就是R-L滤波函数截窗振荡效应所带来的Gibb’s现象。

（4）MATLAB程序

2022年1月11日修改：
十分感谢评论区 tdxt 指出的错误：第31行和第86行的 i 应改为 theta(i)，程序部分已修改。
  

1. %% CT反投影法
2. %%==========================Part 1==========================
3. %% 1、直接反投影法
4. % 绘制一张简单的图，进行直接反投影法
5. clc,clear,close all
6. % ================================================
7. % 简单图
8. % I=zeros(512,512);
9. % for i=1:512
10. %     for j=1:512
11. %         if ((i-256)*(i-256)+(j-256)*(j-256))<1024
12. %             I(i,j)=1;
13. %         end
14. %     end
15. % end
16. % figure,imshow(I,[]);
17. % ===============================================
19. %============================================
20. % 复杂图
21. I = phantom(512);
22. figure,imshow(I,[]);
23. %============================================
25. % 绘制代表角度的反投影图像
26. R=radon(I,0:179);  %对物体CT扫描180°的数据，每1°扫描一次
27. figure,imshow(R,[]);title('CT图像');
28. % 绘制1、45、90、135、180的投影信号和反投影图像
29. theta=[1 45 90 135];
30. for i=1:length(theta)
31.     r=R(:,theta(i));
32.     II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;%回抹（反投影）之前不滤波
33.     II1{i}=II;
34.     figure,
35.     subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信号']);
36.     subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']);
37. end
38. II2=II1{1}+II1{2}+II1{3}+II1{4};
39. figure,imshow(II2,[]);
41. % 对反投影图像进行叠加
42. r=R(:,1);
43. II_r=iradon([r r],[1 1])/2;
44. k=[30 15 10 1];
45. for j=1:4
46.     A=II_r;
47.     for i=2:k(j):180
48.         r=R(:,i);
49.         II=iradon([r r],[i i],'linear','none')/2;%回抹（反投影）之前不滤波
50.         A=A+II;   
51.     end
52. %     figure,imshow(A,[min(min(I)) max(max(I))]);
53.     figure,imshow((A),[]);
54. end
58. %% =========================Part 2=============================
59. %% 2、滤波反投影法
60. clc,clear,close all
61. % 复杂图
62. I = phantom(512);
63. figure,imshow(I,[]);
65. % 对投影做傅里叶变换
66. R=radon(I,0:179);
67. width=length(R);
68. % 设计R-L滤波器
69. filter=2*[0:round(width/2-1), width/2:-1:1]'/width;
70. % 展示滤波前的CT投影信号
71. figure,imshow(R,[]),title('滤波前的CT投影信号');
73. % 每一列做傅里叶变换
74. r_fft=fft(R,729);
75. % 每一列做傅里叶变换后滤波
76. r_fft_filter=r_fft.*filter;
77. % 滤波后反变换(real取实部)
78. r_fft_filter_v=real(ifft(r_fft_filter));
79. % 展示滤波后的CT投影信号
80. % figure,imshow(r_fft_filter_v,[]),title('滤波后的CT投影信号');
81. figure,imshow(r_fft_filter_v),title('滤波后的CT投影信号');
83. % 绘制1、45、90、135、180的投影信号和反投影图像
84. theta=[1 45 90 135];
85. for i=1:length(theta)
86.     r=r_fft_filter_v(:,theta(i));
87.     II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;
88.     %由于iradon默认有滤波，因此关掉默认的滤波选项（none），直接用上面已经手动滤波的数据进行反投影
89.     II1{i}=II;
90.     figure,
91.     subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信号']);
92.     subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']);
93. end
94. II2=II1{1}+II1{2}+II1{3}+II1{4};
95. figure,imshow(II2,[]);
97. % 对反投影图像进行叠加
98. r=r_fft_filter_v(:,1);
99. II_r=iradon([r r],[1 1])/2;
100. k=[30 15 10 1];
101. for j=1:4
102.     A=II_r;
103.     for i=2:k(j):180
104.         r=r_fft_filter_v(:,i);
105.         II=iradon([r r],[i i],'linear','none')/2;
106.         %由于iradon默认有滤波，因此关掉默认的滤波选项（none），直接用上面已经手动滤波的数据进行反投影
107.         A=A+II;   
108.     end
109. %     figure,imshow(A,[min(min(I)) max(max(I))]);
110.     figure,imshow((A),[]);
111. end
113. % 上面是反投影叠加的分步计算，下面直接输入全部投影数据到iradon函数中进行验证
114. B=iradon(r_fft_filter_v,0:179,'linear','none');
115. figure,
116. imshow(B,[]);title('验证：滤波后图像');

复制代码

  

---

以上就是对CT图像重构的整理结果，更多是自己在学习过程中的理解，若有表述不当的地方欢迎指出。


------------------------------------------------------- 2021.05.27  笔记补充 -----------------------------------------------------------

重构图像（514×514）与原图像（512×512）尺寸不一致的原因分析：
首先需要感谢楼下 study_art 同学提出了这个问题，以及 xjch1900 同学的提醒，一开始确实没考虑过这个问题，后来查看源代码找到了原因，在此更新一下笔记。

①radon函数的matlab文件解释：

1. % Grandfathered syntax
2. %   R = RADON(I,THETA,N) returns a Radon transform with the
3. %   projection computed at N points. R has N rows. If you do not
4. %   specify N, the number of points the projection is computed at
5. %   is:
6. %
7. %        2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
8. %
9. %   This number is sufficient to compute the projection at unit
10. %   intervals, even along the diagonal.

复制代码

我们一开始先使用了radon函数来模拟了CT的扫描过程，得到的结果是每个角度的投影数据。假设我们扫描的是一张正方形的图像，大小为512×512，如果我们要把扫描的数据都存储下来，就必须设置一个足够长的数组来记录这些数据，显然，当扫描方向为45°角时得到的扫描数据是最长的（如下图所示），需要的数组长度至少为 512×√2≈725 。

回到matlab的radon函数，从上面的文件解释中可看到radon的语法中的N便是用来定义这个最长长度的。但当我们没有对这个N有专门的限制时，该函数就会自动计算出这个“最长长度”，为了涵盖图像形状的不同情况（如有的图像为长方形），它设置了一个公用的计算公式：
N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
【其中ceil函数表示“向上取整”；floor表示“向下取整”；norm表示“范数计算”，用在这里其实就是计算正方形斜边的长度】
如果套入这个公式，计算出来的N为：

1. size(I)=[512 512];
2. N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
4. N =
5.    
6.    729

复制代码

这比自己的理论值725大了一些，但实际上这是为了涵盖所有图像情况造成的，这并不会对最后的结果带来太大影响。
在正文最后的matlab程序中，我并没有专门定义这个N，因此radon函数自动把N通过公式计算出来了，大小为729，同时还定义了180个扫描角度，所以得到的所有投影数据的矩阵大小为729×180.

②iradon函数的matlab文件解释：

1. %   I = IRADON(R,THETA,INTERPOLATION,FILTER,FREQUENCY_SCALING,OUTPUT_SIZE)   
2. %   OUTPUT_SIZE is a scalar that specifies the number of rows and columns in
3. %   the reconstructed image.  If OUTPUT_SIZE is not specified, the size is
4. %   determined from the length of the projections:
5. %
6. %       OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))
7. %
8. %   If you specify OUTPUT_SIZE, IRADON reconstructs a smaller or larger
9. %   portion of the image, but does not change the scaling of the data.

复制代码

在使用iradon函数进行图像重构时，同样有一个output_size的参数可以设置（表示输出图像的大小），我们可以直接设置为原始图像的大小（即512），这样得到的图像就跟原图像大小一模一样了。但在正文最后部分，我同样没有专门设置图像大小，所以该函数就会自动套入一条涵盖了所有图像形状的计算公式：
OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))
由于我们在第一步用radon函数模拟CT扫描时没有专门定义N，得到的投影数据的矩阵R大小为729×180，因此直接套入这个output_size公式后就直接算出：

1. size(R,1)=729;
2. OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))
4. OUTPUT_SIZE =
5.    
6.            514

复制代码

所以最后得到的重构图像大小就变成了514×514.

③解决方法

•  使用radon函数时，自己定义N。如在正文的例子中，直接在radon函数中输入N的理论计算值725，这样在重构iradon时就不需要设置output_size了，输出结果直接为512×512
  
• 使用radon函数时，不定义N，让函数自己计算图像的N，然后在重构iradon时设置output_size的大小为原始图像大小512，输出结果同样为512×512
  

 备注：若radon和iradon函数都不专门设置时，会出现重构图像与原图像尺寸大小不一致的情况，但这并不会改变重构图像本身的缩放比例，因此若没有专门的尺寸大小必须前后一致的严格要求，可以不用专门设置。

来源：https://blog.csdn.net/qq_33924470/article/details/106874678
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！