滤波反投影重建算法(FBP)实现及应用(matlab)
2017-09-28 22:24
507 查看
滤波反投影重建算法实现及应用(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. 滤波器(滤波函数)和内插函数的选取
由于直接使用反投影算法会存在两个对实验结果影响很不好的因素: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 && 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.
相关文章推荐
- 滤波反投影图像重建算法
- 【图】二分图最大匹配算法的应用及Matlab实现
- 图像算法之十二:非局部均值滤波及其Matlab实现
- 数字图像处理,经典滤波算法去噪对比实验(Matlab实现)
- otsu自适应阈值分割的算法描述和opencv实现,及其在肤色检测中的应用
- Matlab中图像处理实例:灰度变换,空域滤波,频域滤波,傅里叶变换的实现
- C#数据结构和算法学习系列六----堆栈、堆栈的实现和应用
- 基于肤色和眼睛定位的人脸检测算法——MATLAB实现
- PHP实现的简单排列组合算法应用示例
- 粒子滤波(Particle filter)算法简介及MATLAB实现
- 算法应用:用递归结构实现穷举服务器文件
- 神经网络之感知器算法简单介绍和MATLAB简单实现
- decode函数应用,对角线算法实现表中数据的行列转换
- matlab黄金分割算法实现与解析
- 【算法导论】八皇后问题的算法实现(C、MATLAB、Python版)
- [置顶] 【算法 机器学习】MATLAB、R、python三种编程语言实现简单线性回归算法比较
- matlab实现图割算法中的最大流最小割Max-flow/min-cut问题(一)
- 【原创】机器学习之PageRank算法应用与C#实现(1)算法介绍
- MATLAB中实现图像的空间域滤波和频率域滤波
- 卡尔曼滤波简介说明及其算法MATLAB实现代码