图像重建图像重建之ASTRA算法

Posted 2021-05-18 Matlab走起

滤波反投影重建算法实现及应用（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滤波函数滤波计算简单，避免了大量的正弦、余弦计算，得到的采样序列分是分段现行的，并没有明显的降低图像质量，所以重建图像轮廓清楚，空间分辨率高。

常见的插值方法有最近邻插值和双线插值，最近邻插值即将离散点中间的缺失值用离它最近的整数处的投影值来替代。

 
% load a phantom image
im = phantom(256);
% and flatten it to a vector
x = im(:);
 
%% Setting up the geometry
% projection geometry
proj_geom = astra_create_proj_geom('parallel', 1, 256, linspace2(0,pi,180));
% object dimensions
vol_geom  = astra_create_vol_geom(256,256);
 
%% Generate projection data
% Create the Spot operator for ASTRA using the GPU.
W = opTomo('cuda', proj_geom, vol_geom);
 
p = W*x;
 
% reshape the vector into a sinogram
sinogram = reshape(p, W.proj_size);
imshow(sinogram, []);
 
 
%% Reconstruction
% We use a least squares solver lsqr from Matlab to solve the 
% equation W*x = p.
% Max number of iterations is 100, convergence tolerance of 1e-6.
y = lsqr(W, p, 1e-6, 100);
 
% the output is a vector, so we reshape it into an image
reconstruction = reshape(y, W.vol_size);
 
subplot(1,3,1);
imshow(reconstruction, []);
title('Reconstruction');
 
subplot(1,3,2);
imshow(im, []);
title('Ground truth');
 
% The transpose of the operator corresponds to the backprojection.
backProjection = W'*p;
subplot(1,3,3);
imshow(reshape(backProjection, W.vol_size), []);
title('Backprojection');

完整代码添加QQ1575304183