| 滤波反投影重建算法（FBP）实现及应用（matlab） | 您所在的位置：网站首页 › ct图像处理技术有哪些 › 滤波反投影重建算法（FBP）实现及应用（matlab） |
滤波反投影重建算法（FBP）实现及应用（matlab）
|
滤波反投影重建算法实现及应用（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 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 |
【本文地址】