滤波反投影重建算法实现及应用（matlab）
 
1. 滤波反投影重建算法原理
 
 
 
滤波反投影重建算法常用在CT成像重建中，背后的数学原理是傅立叶变换：对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换。（傅立叶中心切片定理）
 
 
CT重建算法大致分为解析重建算法和迭代重建算法，随着CT技术的发展，重建算法也变得多种多样，各有各的有特点。本文使用目前应用最广泛的重建算法——滤波反投影算法（FBP）作为模型的基础算法。FBP算法是在傅立叶变换理论基础之上的一种空域处理技术。它的特点是在反投影前将每一个采集投影角度下的投影进行卷积处理，从而改善点扩散函数引起的形状伪影，重建的图像质量较好。
 
 
上图应可以清晰的描述傅立叶中心切片定理的过程：对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换
 
傅立叶切片定理的意义在于，通过投影上执行傅立叶变换，可以从每个投影中得到二维傅立叶变换。从而投影图像重建的问题，可以按以下方法进行求解：采集不同时间下足够多的投影（一般为180次采集），求解各个投影的一维傅立叶变换，将上述切片汇集成图像的二维傅立叶变换，再利用傅立叶反变换求得重建图像。
 
 
 
投影相关知识请参考fbp的matlab实现
 
 
2. 滤波反投影重建算法过程（以平行束为例）
 
投影重建的过程是，先把投影由线阵探测器上获得的投影数据进行一次一维傅立叶变换，再与滤波器函数进行卷积运算，得到各个方向卷积滤波后的投影数据；然后把它们沿各个方向进行反投影，即按其原路径平均分配到每一矩阵单元上，进行重叠后得到每一矩阵单元的CT值；再经过适当处理后得到被扫描物体的断层图像
 算法步骤如下：
 1. 将原始投影进行一次一维傅立叶变换
 2. 设计合适的滤波器，在φ_i的角度下将得到原始投影p(x_r,φ_i)进行卷积滤波，得到滤波后的投影。
 3. 将滤波后的投影进行反投影，得到满足x_r=r cos⁡（(θ - φ_i)）方向上的原图像的密度。
 4. 将所有反投影进行叠加，得到重建后的投影。
 
3. 滤波器(滤波函数)和内插函数的选取
 
由于直接使用反投影算法会存在两个对实验结果影响很不好的因素：
 
不准确的数据重建图像就会产生各种伪影。
投影的数据是天然离散的，处理不当的话会产生很大的误差。
 
常见的滤波器有R-S滤波函数和S-L滤波函数。R-L滤波函数滤波计算简单，避免了大量的正弦、余弦计算，得到的采样序列分是分段现行的，并没有明显的降低图像质量，所以重建图像轮廓清楚，空间分辨率高。
 
常见的插值方法有最近邻插值和双线插值，最近邻插值即将离散点中间的缺失值用离它最近的整数处的投影值来替代。
 
4. FBP的matlab实现
 
使用R-L滤波器和最近邻插值方法
 
clc,clear;
%% 各个参数信息
%重建后图片像素个数
M=512;%建议先和接收传感器的个数一样

%旋转的角度 180次旋转
theta=xlsread('angles_180_3.xlsx','Sheet1','a2:a181');

%投影为512行，180列的数据，列的数据对应每一次旋转，512个接收传感器
R=xlsread('Accessory_A.xls','附件3');

%由投影进行反变换
size(R,1);%512

% 设置快速傅里叶变换的宽度
width = 2^nextpow2(size(R,1));  

%% 对投影做快速傅里叶变换并滤波
%傅立叶变换
proj_fft = fft(R, width);

% filter 滤波
% R-L是一种基础的滤波算法
filter = 2*[0:(width/2-1), width/2:-1:1]'/width;

% 滤波后结果 proj_filtered
proj_filtered = zeros(width,180);
for i = 1:180
    proj_filtered(:,i) = proj_fft(:,i).*filter;
end
figure
subplot(1,2,1),imshow(proj_fft),title('傅立叶变换')
subplot(1,2,2),imshow(proj_filtered),title('傅立叶变换+滤波')

%% 逆快速傅里叶变换并反投影
% 逆快速傅里叶变换 proj_ifft
proj_ifft = real(ifft(proj_filtered)); 
figure,imshow(proj_ifft),title('逆傅立叶变换')

%反投影到x轴，y轴
fbp = zeros(M); % 假设初值为0
for i = 1:180
    rad = theta(i);%弧度， %这个rad 是投影角，不是投影线与x轴夹角，他们之间相差 pi/2
    for x = 1:M
        for y = 1:M
            %{
            %最近邻插值法
            t = round((x-M/2)*cos(rad)-(y-M/2)*sin(rad));%将每个元素舍X入到最接近的整数。
            if t<size(R,1)/2 && t>-size(R,1)/2
                fbp(x,y)=fbp(x,y)+proj_ifft(round(t+size(R,1)/2),i);
            end
            %}
            t_temp = (x-M/2) * cos(rad) - (y-M/2) * sin(rad)+M/2  ;
             %最近邻插值法
            t = round(t_temp) ;
            if t>0 && t<=512
                fbp(x,y)=fbp(x,y)+proj_ifft(t,i);
            end
        end
    end
end
fbp = (fbp*pi)/180;%512x512 原图像每个像素位置的密度

xlswrite('problem2_origin.xlsx',fbp,'Sheet1');%将得到的重建后的图像数据写入
% 显示结果
figure,imshow(fbp),title('反投影变换后的图像')


 
其中每一个文件都有作用的说明，如需源文件请留言哈。（如有错误请指正）
 fbp算法实现案例
 
github代码下载 别忘了给star哟~
 
参考文献：
 
 
 
【1】 范慧赟.CT 图像滤波反投影重建算法的研究[D].硕士学位论文,西北工业大学,2007.
 【2】 余晓锷,龚剑,马建华等.CT 原理与技术[M].北京:科学出版社,2013,95-97.
 【3】 毛小渊. 二维CT图像重建算法研究[D].南昌航空大学,2016.
 【4】 洪虹. CT中金属伪影的校正研究[D].南方医科大学,2013.
 【5】 范慧赟. CT图像滤波反投影重建算法的研究[D].西北工业大学,2007.
 

 
摘要：X-CT(即X射线计算机断层成像)图像重建的卷积反投影原理对医学影像专业学生较为抽象难懂，本文在简化条件下采用图解法较为系统地描述该原理，达到直观易懂的目的。
 
1　X-CT图像重建的卷积反投射原理
 
X-CT成像原理在不少刊物上有详细介绍[1-3]。所谓CT图像重建是按照采集后的数据，求解图像矩阵中像素，然后重新构造图像的过程；而图像矩阵的求解由计算机完成，在X-CT图像重建的解析法中，当前最常用的是采用卷积运算的滤波反投影法，也称卷积反射影法。
 
卷积反投影重建图像时，先把由检测器上获得的原始数据与一个滤波函数进行了卷积运算，得到各方向卷积的投影函数；然后再把它们从各方向进行反投影，即按其原路径平均分配到每一矩阵元上，进行叠加后得到每一矩阵元的CT值；再经过适当处理后就可以得到被扫描物体的断层图像，卷积反投影可消除单纯的反投影产生的边缘失锐效应，补偿投影中的高频成分和降低投影中心密度，并保证重建图像边缘清晰和内部分布均匀。
 
卷积反投影法的数学表达式可由雷当投影定理推导出[2]；下面直接给出表达式。
 
如图1所示，取物体某断层面为x-y坐标系，X射线管A发出的单能窄束X射线束穿过衰减系数为μ(x,y)的物体断层，沿S路径X线束的投影可写为：
 

 
式中I0是X线管输出强度，Im是检测器测得的X线束强度。
 

 
图1　被测物体断面
 
由上式可知，投影值P(x,y)就是衰减系数μ(x,y)在断层平面的二维积分，若测出投影值P(x,y)，可求出该断层平面上衰减系数μ(x,y)的二维分布。
 
在X线束的平移-旋转扫描中，X线的投影值P(x,y)总是与X线束路径S有关的，在图1中，引入一个新坐标系R-θ来描述X线束路径S在x-y坐标系中的位置；取R为X线束路径S到坐标中心的距离，与y轴夹角为θ，则用直线方程：x
 
cosθ+y
 
sinθ=R描述路径S。若取δ(t)-函数为筛选因子，则μ(x,y)δ(t)为指定某一路径的衰减系数，这时，某一角度θ相对应的投影可写为：
 

 
式中的Pθ(R,θ)仅是R的函数，记为P(R)；若取滤波函数为Φ(t)，为了消除边缘失锐效应，需要用Φ(t)对投影Pθ(R，θ)进行有效地滤波，用滤波函数与反投影信号相加，使各投影信号在其两旁出现正和负的成分，形成所谓滤波反投影信号，由于各次信号迭加时，这些正和负的微小脉冲几乎互相抵消，相当于把失锐部分消去，使图像与实物相似；这种滤波等效于将Pθ(R，θ)与Φ(t)进行卷积运算。
 

 
式中P(R)*Φ(R)表示两者的卷积。
 
对应角度θ的反投影可写为：
 

 
将角度θ从0变化到π相应的所有反射影值相加，可得到断层图像重建的衰减系数：
 

 
2　卷积反投影的图解法
 
图像重建的卷积反投影原理及解析式的推导，对于医学院校医学影像专业的学生较为抽象难懂；在一定简化条件下，下面通过图解法系统地描述卷积反投影原理，达到直观易懂的目的，获得较好的课堂教学效果。
 
如图2(a)所示，设X线束穿过一半径为r衰减系数为μ的均匀圆形物体，X线束分别沿y轴方向(θ=0°)和X轴方向(θ=90°)各取得五次投影，如图2(b)，沿Y轴投影：
 

 

 
(C)D4=D8=6.0、D5=D3=9.0、D6=10.0
 
图2　沿Y轴的五个投影值(a～c)
 
上式得到五个投影值分别为(令r=1):1.20,1.83,2.0,1.83,1.20；为简化取整，令μ=5，则这五个投影值为：6.0，9.0，10.0，9.0，6.0；如图2(c)所示，沿X轴投影时，同理得到同样的五个投影值，如果将上述投影值单纯反投影，进行叠加后所得图像矩阵，如图3所示。图像均匀度按下式计算[4]：
 
W=[1-(μmax-μmin)/(μmax+μmin)]×100%
 
式中μmax、μmin分别为测量区域内矩阵元最大值和最小值，对于图3，W=98%，图像均匀度较低；从图3可看出沿投影方向有伪迹，图像中心密度增高。
 

 
图3　单纯反投影所得图像矩阵
 
设滤波函数Φ(a,n)为：
 

 
式中a为图像每条平行射线间距，n为图像矩阵单列或单行的长度；Φ(a,n)函数图像如图4所示。
 
由于X-CT数据采集为离散的，且滤波函数也取离散形式，故卷积公式可写成离散形式，并由积分式化成求和式，有：
 

 
式中P(n)、Φ(n)分别是投影值和滤波离散函数，N是图像矩阵单列或单行的长度，这里N=12；用Dn表示P(n)的各投影值，用Fn表示Φ(n)的各函数值，当n=0,1,2......11时，对于P(n):
 
D4=6.0,D5=9.0,D6=10.0,D7=9.0,
 
D8=6.0，其余Dn=0。如图2(c)所示；而对于Φ(n):F0=1.27,F1=-0.42,F2=-0.08,F3=0.04，其余Fn=0。见图4。
 

 
图4　滤波函数图像
 
设CDn为卷积值，P(n)与Φ(n)的卷积运算由图5(a)-(g)给出图解过程，卷积时可看成二个函数先相向后背向运动：在它产未相交前，卷积值CD0为零；第一点相交时，可把两个函数的相对应部分相乘，得出第一个卷积值CD1；然后滤波函数再向右移动一位，把二个交叉值分别相乘再相加，求得第二个卷积值CD2；再向右移动，重复上述过程可求出第三个卷积值CD3……，直到两个函数完全分开，卷积值又为零，卷积运算后得到的各个原始数值CDn如图6所示。
 

 
图5　P(n)与Φ(n)的卷积运算的图解法(a～g)
 

 
图6　卷积后得到的原始数值CDn
 
将卷积后的原始数值反投影，进行叠加后的图像短阵如图7所示，其图像均匀度W=91%，可见单纯珠反投影图像相比，其图像均匀度有所提高；同时由图7可见，卷积反投影后图像边缘的高频分量增高，图像中心宽度有所降低，若增加反投影方向数量，所得图像的均匀度更高，例如：X线束沿θ=0°，θ=45°和θ=90°投影时，卷积反投影后图像均匀度W=92%。
 

 
图7　卷积反投影后得到的图像矩阵
 
在课堂教学当中，将上述图解过程制作成幻灯片投影，能使学生对X-CT图像重建的卷积反投影原理有较直观透彻的理解；特别是对解释卷积反投影法消除边缘失锐，降低中心密度的作用，有较好的教学效果。

 
濾波反投影重建算法實現及應用(matlab)
 
1. 濾波反投影重建算法原理
 
濾波反投影重建算法常用在CT成像重建中，背后的數學原理是傅立葉變換：對投影的一維傅立葉變換等效於對原圖像進行二維的傅立葉變換。(傅立葉中心切片定理)
 
CT重建算法大致分為解析重建算法和迭代重建算法，隨着CT技術的發展，重建算法也變得多種多樣，各有各的有特點。本文使用目前應用最廣泛的重建算法——濾波反投影算法(FBP)作為模型的基礎算法。FBP算法是在傅立葉變換理論基礎之上的一種空域處理技術。它的特點是在反投影前將每一個采集投影角度下的投影進行卷積處理，從而改善點擴散函數引起的形狀偽影，重建的圖像質量較好。
 

 
上圖應可以清晰的描述傅立葉中心切片定理的過程：對投影的一維傅立葉變換等效於對原圖像進行二維的傅立葉變換
 
傅立葉切片定理的意義在於，通過投影上執行傅立葉變換，可以從每個投影中得到二維傅立葉變換。從而投影圖像重建的問題，可以按以下方法進行求解：采集不同時間下足夠多的投影(一般為180次采集)，求解各個投影的一維傅立葉變換，將上述切片匯集成圖像的二維傅立葉變換，再利用傅立葉反變換求得重建圖像。
 
2. 濾波反投影重建算法過程(以平行束為例)
 
投影重建的過程是，先把投影由線陣探測器上獲得的投影數據進行一次一維傅立葉變換，再與濾波器函數進行卷積運算，得到各個方向卷積濾波后的投影數據；然后把它們沿各個方向進行反投影，即按其原路徑平均分配到每一矩陣單元上，進行重疊后得到每一矩陣單元的CT值；再經過適當處理后得到被掃描物體的斷層圖像
 
算法步驟如下：
 
1. 將原始投影進行一次一維傅立葉變換
 
2. 設計合適的濾波器，在φ_i的角度下將得到原始投影p(x_r,φ_i)進行卷積濾波，得到濾波后的投影。
 
3. 將濾波后的投影進行反投影，得到滿足x_r=r cos⁡((θ - φ_i))方向上的原圖像的密度。
 
4. 將所有反投影進行疊加，得到重建后的投影。
 
3. 濾波器(濾波函數)和內插函數的選取
 
由於直接使用反投影算法會存在兩個對實驗結果影響很不好的因素：
 
1. 不准確的數據重建圖像就會產生各種偽影。
 
2. 投影的數據是天然離散的，處理不當的話會產生很大的誤差。
 
常見的濾波器有R-S濾波函數和S-L濾波函數。R-L濾波函數濾波計算簡單，避免了大量的正弦、余弦計算，得到的采樣序列分是分段現行的，並沒有明顯的降低圖像質量，所以重建圖像輪廓清楚，空間分辨率高。
 
常見的插值方法有最近鄰插值和雙線插值，最近鄰插值即將離散點中間的缺失值用離它最近的整數處的投影值來替代。
 
4. FBP的matlab實現
 
使用R-L濾波器和最近鄰插值方法
 
clc,clear;
 
%% 各個參數信息
 
%重建后圖片像素個數
 
M=512;%建議先和接收傳感器的個數一樣
 
%旋轉的角度 180次旋轉
 
theta=xlsread('angles_180_3.xlsx','Sheet1','a2:a181');
 
%投影為512行，180列的數據，列的數據對應每一次旋轉，512個接收傳感器
 
R=xlsread('Accessory_A.xls','附件3');
 
%由投影進行反變換
 
size(R,1);%512
 
% 設置快速傅里葉變換的寬度
 
width = 2^nextpow2(size(R,1));
 
%% 對投影做快速傅里葉變換並濾波
 
%傅立葉變換
 
proj_fft = fft(R, width);
 
% filter 濾波
 
% R-L是一種基礎的濾波算法
 
filter = 2*[0:(width/2-1), width/2:-1:1]'/width;
 
% 濾波后結果 proj_filtered
 
proj_filtered = zeros(width,180);
 
for i = 1:180
 
proj_filtered(:,i) = proj_fft(:,i).*filter;
 
end
 
figure
 
subplot(1,2,1),imshow(proj_fft),title('傅立葉變換')
 
subplot(1,2,2),imshow(proj_filtered),title('傅立葉變換+濾波')
 
%% 逆快速傅里葉變換並反投影
 
% 逆快速傅里葉變換 proj_ifft
 
proj_ifft = real(ifft(proj_filtered));
 
figure,imshow(proj_ifft),title('逆傅立葉變換')
 
%反投影到x軸，y軸
 
fbp = zeros(M); % 假設初值為0
 
for i = 1:180
 
rad = theta(i);%弧度， %這個rad 是投影角，不是投影線與x軸夾角，他們之間相差 pi/2
 
for x = 1:M
 
for y = 1:M
 
%{
 
%最近鄰插值法
 
t = round((x-M/2)*cos(rad)-(y-M/2)*sin(rad));%將每個元素舍X入到最接近的整數。
 
if t-size(R,1)/2
 
fbp(x,y)=fbp(x,y)+proj_ifft(round(t+size(R,1)/2),i);
 
end
 
%}
 
t_temp = (x-M/2) * cos(rad) - (y-M/2) * sin(rad)+M/2 ;
 
%最近鄰插值法
 
t = round(t_temp) ;
 
if t>0 && t<=512
 
fbp(x,y)=fbp(x,y)+proj_ifft(t,i);
 
end
 
end
 
end
 
end
 
fbp = (fbp*pi)/180;%512x512 原圖像每個像素位置的密度
 
xlswrite('problem2_origin.xlsx',fbp,'Sheet1');%將得到的重建后的圖像數據寫入
 
% 顯示結果
 
figure,imshow(fbp),title('反投影變換后的圖像')
 
其中每一個文件都有作用的說明，如需源文件請留言哈。(如有錯誤請指正)
 
fbp算法實現案例
 
github代碼下載 別忘了給star喲~
 
參考文獻：
 
【1】 范慧贇.CT 圖像濾波反投影重建算法的研究[D].碩士學位論文,西北工業大學,2007.
 
【2】 余曉鍔,龔劍,馬建華等.CT 原理與技術[M].北京:科學出版社,2013,95-97.
 
【3】 毛小淵. 二維CT圖像重建算法研究[D].南昌航空大學,2016.
 
【4】 洪虹. CT中金屬偽影的校正研究[D].南方醫科大學,2013.
 
【5】 范慧贇. CT圖像濾波反投影重建算法的研究[D].西北工業大學,2007.

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

      CT图像的形成和重构，在数学上的描述分别为拉东变换和拉东反变换。简单点说，拉动变换是用于描述通过X射线扫描物体形成CT图像这个过程的一种数学表示，而拉东反变换描述的则是对CT投影数据进行重构，还原成物体图像的一种数学方法。本文章主要讲的是CT图像重构的方法，也就是拉东反变换的具体实现方法，当然可能也有小伙伴想了解一下什么是拉东变换，这里贴上传送门：拉东变换。
 
其中拉东反变换在数学上的表达式如下：
 
 
由公式可看出拉东逆变换在数学上实际上是一个二维傅里叶逆变换，但由于二维傅里叶逆变换计算量大，耗时长，因此计算机对于拉东逆变换的具体实现方法主要为反投影法，反投影法又分为直接反投影法和滤波反投影法。由此可见CT重构方法主要为以下三种：
 
（1）傅里叶逆变换法
 
（2）直接反投影法
 
（3）滤波反投影法
 
下面将具体对这三种图像重构方法进行解释。
 
 

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

  
  
 
 
 
 

  
  
 
 
根据该理论，傅里叶逆变换法可以简单分成以下步骤：
 
① 假设每旋转1°就扫描一次，当对物体扫描了180°之后，我们就能得到180个投影信号（就是180根投影线）→在临床上，若使用平行扫描CT，我们拿到手的数据就是这个（在数学上，就是对图像进行拉东变换）
 
② 对180个投影信号进行一维傅里叶变换
 
③ 对②得到的180个一维频域信号，根据相对应的扫描角度，在空间中旋转排列，拼成一个二维频域空间（如图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。
 
      
 
 
 

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

  
  
 
 
这个方法的好处是，把两次二维傅里叶变换变成了两次一维傅里叶变换，计算速度大大提升。于是滤波反投影法的核心问题就变成了如何选择一个合适的滤波器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)，程序部分已修改。
 
 
 
%% CT反投影法
%%==========================Part 1==========================
%% 1、直接反投影法
% 绘制一张简单的图，进行直接反投影法
clc,clear,close all
% ================================================
% 简单图
% I=zeros(512,512);
% for i=1:512
%     for j=1:512
%         if ((i-256)*(i-256)+(j-256)*(j-256))<1024
%             I(i,j)=1;
%         end
%     end
% end
% figure,imshow(I,[]);
% ===============================================

%============================================
% 复杂图
I = phantom(512); 
figure,imshow(I,[]);
%============================================

% 绘制代表角度的反投影图像
R=radon(I,0:179);  %对物体CT扫描180°的数据，每1°扫描一次
figure,imshow(R,[]);title('CT图像');
% 绘制1、45、90、135、180的投影信号和反投影图像
theta=[1 45 90 135];
for i=1:length(theta)
    r=R(:,theta(i));
    II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;%回抹（反投影）之前不滤波
    II1{i}=II;
    figure,
    subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信号']);
    subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']);
end
II2=II1{1}+II1{2}+II1{3}+II1{4};
figure,imshow(II2,[]);

% 对反投影图像进行叠加
r=R(:,1);
II_r=iradon([r r],[1 1])/2;
k=[30 15 10 1];
for j=1:4
    A=II_r;
    for i=2:k(j):180
        r=R(:,i);
        II=iradon([r r],[i i],'linear','none')/2;%回抹（反投影）之前不滤波
        A=A+II;   
    end
%     figure,imshow(A,[min(min(I)) max(max(I))]);
    figure,imshow((A),[]);
end



%% =========================Part 2=============================
%% 2、滤波反投影法
clc,clear,close all
% 复杂图
I = phantom(512); 
figure,imshow(I,[]);

% 对投影做傅里叶变换
R=radon(I,0:179);
width=length(R);
% 设计R-L滤波器
filter=2*[0:round(width/2-1), width/2:-1:1]'/width;
% 展示滤波前的CT投影信号
figure,imshow(R,[]),title('滤波前的CT投影信号');

% 每一列做傅里叶变换
r_fft=fft(R,729); 
% 每一列做傅里叶变换后滤波
r_fft_filter=r_fft.*filter; 
% 滤波后反变换(real取实部)
r_fft_filter_v=real(ifft(r_fft_filter));
% 展示滤波后的CT投影信号
% figure,imshow(r_fft_filter_v,[]),title('滤波后的CT投影信号');
figure,imshow(r_fft_filter_v),title('滤波后的CT投影信号');

% 绘制1、45、90、135、180的投影信号和反投影图像
theta=[1 45 90 135];
for i=1:length(theta)
    r=r_fft_filter_v(:,theta(i));
    II=iradon([r r],[theta(i) theta(i)],'linear','none')/2;
    %由于iradon默认有滤波，因此关掉默认的滤波选项（none），直接用上面已经手动滤波的数据进行反投影
    II1{i}=II;
    figure,
    subplot(1,2,1),imshow(r,[]);title([num2str(theta(i)) '°投影信号']);
    subplot(1,2,2),imshow(II,[]);title(['往' num2str(theta(i)) '°方向反投影']);
end
II2=II1{1}+II1{2}+II1{3}+II1{4};
figure,imshow(II2,[]);

% 对反投影图像进行叠加
r=r_fft_filter_v(:,1);
II_r=iradon([r r],[1 1])/2;
k=[30 15 10 1];
for j=1:4
    A=II_r;
    for i=2:k(j):180
        r=r_fft_filter_v(:,i);
        II=iradon([r r],[i i],'linear','none')/2;
        %由于iradon默认有滤波，因此关掉默认的滤波选项（none），直接用上面已经手动滤波的数据进行反投影
        A=A+II;   
    end
%     figure,imshow(A,[min(min(I)) max(max(I))]);
    figure,imshow((A),[]);
end

% 上面是反投影叠加的分步计算，下面直接输入全部投影数据到iradon函数中进行验证
B=iradon(r_fft_filter_v,0:179,'linear','none');
figure,
imshow(B,[]);title('验证：滤波后图像');


 
 
 
 

以上就是对CT图像重构的整理结果，更多是自己在学习过程中的理解，若有表述不当的地方欢迎指出。
 
 
 
------------------------------------------------------- 2021.05.27  笔记补充 -----------------------------------------------------------
 
 
重构图像（514×514）与原图像（512×512）尺寸不一致的原因分析：
 
首先需要感谢楼下 study_art 同学提出了这个问题，以及 xjch1900 同学的提醒，一开始确实没考虑过这个问题，后来查看源代码找到了原因，在此更新一下笔记。
 
 
①radon函数的matlab文件解释：
 
% Grandfathered syntax
%   R = RADON(I,THETA,N) returns a Radon transform with the
%   projection computed at N points. R has N rows. If you do not
%   specify N, the number of points the projection is computed at
%   is:
%
%        2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3
%
%   This number is sufficient to compute the projection at unit
%   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为：
 
size(I)=[512 512];
N=2*ceil(norm(size(I)-floor((size(I)-1)/2)-1))+3

N =
   
   729
 
 
这比自己的理论值725大了一些，但实际上这是为了涵盖所有图像情况造成的，这并不会对最后的结果带来太大影响。
 
在正文最后的matlab程序中，我并没有专门定义这个N，因此radon函数自动把N通过公式计算出来了，大小为729，同时还定义了180个扫描角度，所以得到的所有投影数据的矩阵大小为729×180.
 
 
②iradon函数的matlab文件解释：
 
%   I = IRADON(R,THETA,INTERPOLATION,FILTER,FREQUENCY_SCALING,OUTPUT_SIZE)   
%   OUTPUT_SIZE is a scalar that specifies the number of rows and columns in
%   the reconstructed image.  If OUTPUT_SIZE is not specified, the size is
%   determined from the length of the projections:
%
%       OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))
%
%   If you specify OUTPUT_SIZE, IRADON reconstructs a smaller or larger
%   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公式后就直接算出：
 
size(R,1)=729;
OUTPUT_SIZE = 2*floor(size(R,1)/(2*sqrt(2)))

OUTPUT_SIZE =
   
           514
 
所以最后得到的重构图像大小就变成了514×514.
 
 
③解决方法
 
 使用radon函数时，自己定义N。如在正文的例子中，直接在radon函数中输入N的理论计算值725，这样在重构iradon时就不需要设置output_size了，输出结果直接为512×512
使用radon函数时，不定义N，让函数自己计算图像的N，然后在重构iradon时设置output_size的大小为原始图像大小512，输出结果同样为512×512
 备注：若radon和iradon函数都不专门设置时，会出现重构图像与原图像尺寸大小不一致的情况，但这并不会改变重构图像本身的缩放比例，因此若没有专门的尺寸大小必须前后一致的严格要求，可以不用专门设置。