滤波反投影重建算法实现及应用（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:180proj_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:180rad = theta(i);%弧度， %这个rad 是投影角，不是投影线与x轴夹角，他们之间相差 pi/2for x = 1:Mfor 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)/2fbp(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<=512fbp(x,y)=fbp(x,y)+proj_ifft(t,i);endendend
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.

滤波反投影重建算法（FBP）实现及应用（matlab）相关推荐

matlab fbp fan arc,滤波反投影重建算法（FBP）实现及应用（matlab）
滤波反投影重建算法实现及应用(matlab) 1. 滤波反投影重建算法原理滤波反投影重建算法常用在CT成像重建中,背后的数学原理是傅立叶变换:对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换 ...
【youcans 的 OpenCV 例程 200 篇】112. 滤波反投影重建图像
欢迎关注『youcans 的 OpenCV 例程 200 篇』系列,持续更新中欢迎关注『youcans 的 OpenCV学习课』系列,持续更新中 [youcans 的 OpenCV 例程 2 ...
【OpenCV 例程 300 篇】112. 滤波反投影重建图像
专栏地址:『youcans 的 OpenCV 例程 300篇 - 总目录』 [第 7 章:图像复原与重建] 110. 投影和雷登变换 111. 雷登变换反投影重建图像 112. 滤波反投影重建图像 [ ...
直接反投影 matlab,濾波反投影重建算法（FBP）實現及應用（matlab）
濾波反投影重建算法實現及應用(matlab) 1. 濾波反投影重建算法原理濾波反投影重建算法常用在CT成像重建中,背后的數學原理是傅立葉變換:對投影的一維傅立葉變換等效於對原圖像進行二維的傅立葉變換 ...
FBP 滤波反投影重建
一.为什么光有反投影不足以得到质量较好的图像,下面这个图像能很好的说明问题,图像的边缘模糊二-中心切片定理二维图像的中心切片定理指出:二维函数 f(x, y) 的投影 p(s) 之傅里叶变换 P( ...
运用滤波反投影的方法对图像进行重建matlab仿真
目录 1.算法描述 2.仿真效果预览 3.MATLAB部分代码预览 4.完整MATLAB程序 1.算法描述直接由正弦图得到反投影图像,会存在严重的模糊,这是早期 CT 系统所存在的问题.傅立叶中心切 ...
CT 系统参数标定及反投影重建成像-2017数模国赛论文A298编程分析
CT 系统参数标定及反投影重建成像-2017数模国赛论文A298编程分析之前的同学已经讲解清楚了这篇论文建模的主要思路,我主要讲解代码对建模思路的实现. 本文提到的论文下载地址:http://dxs ...
基于matlab的医学成像技术滤波反投影仿真,包括直接反投影,S-L滤波,R-L滤波,Lewitt滤波
目录 1.算法描述 2.仿真效果预览 3.MATLAB部分代码预览 4.完整MATLAB程序 1.算法描述医学成像技术滤波反投影含R-L滤波, R-S滤波,Lewitt滤波重建后图像清晰. ...
第5章 Python 数字图像处理(DIP) - 图像复原与重建17 - 由投影重建图像、雷登变换、投影、反投影、反投影重建
标题由投影重建图像投影和雷登变换 Johann Radon 反投影滤波反投影重建由投影重建图像本由投影重建图像,主要是雷登变换与雷登把变换的应用,所以也没有太多的研究,只为了保持完整性,而添 ...

滤波反投影重建算法（FBP）实现及应用（matlab）