找回密码
 注册
关于网站域名变更的通知
查看: 642|回复: 3
打印 上一主题 下一主题

基于matlab ASTRA算法图像重建

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-5-13 13:59 | 只看该作者 |只看大图 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
( a, Y$ D: \# B2 D! K/ l# m7 n4 V
一、简介" U9 D0 d9 }) B
滤波反投影重建算法实现及应用(matlab)
& g, Z7 `5 _* u4 `$ v: c! Q( \) D2 _. N0 h. X2 F
  • 滤波反投影重建算法原理/ f/ w+ O. o9 E4 S, t2 |
      滤波反投影重建算法常用在CT成像重建中,背后的数学原理是傅立叶变换:对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换。(傅立叶中心切片定理): [% s- m- e% c
CT重建算法大致分为解析重建算法和迭代重建算法,随着CT技术的发展,重建算法也变得多种多样,各有各的有特点。本文使用目前应用最广泛的重建算法——滤波反投影算法(FBP)作为模型的基础算法。FBP算法是在傅立叶变换理论基础之上的一种空域处理技术。它的特点是在反投影前将每一个采集投影角度下的投影进行卷积处理,从而改善点扩散函数引起的形状伪影,重建的图像质量较好。
/ ~* Y5 O3 i5 w" j9 A8 i' d4 n5 X; `* A! S3 {! [/ z
  g7 y, x0 |3 S! f8 q

/ f* q, j" y0 s3 P" Z" r/ P7 ]- v' ?+ A: d
上图应可以清晰的描述傅立叶中心切片定理的过程:对投影的一维傅立叶变换等效于对原图像进行二维的傅立叶变换( X- Y# M) p% ~) N/ A& y9 a2 Y
9 C/ c% X  G  W! D' w' }
傅立叶切片定理的意义在于,通过投影上执行傅立叶变换,可以从每个投影中得到二维傅立叶变换。从而投影图像重建的问题,可以按以下方法进行求解:采集不同时间下足够多的投影(一般为180次采集),求解各个投影的一维傅立叶变换,将上述切片汇集成图像的二维傅立叶变换,再利用傅立叶反变换求得重建图像。
7 A4 v' \4 Y7 c7 {& O# t. |* e2. 滤波反投影重建算法过程(以平行束为例)
4 `6 ~  x7 t- e: N投影重建的过程是,先把投影由线阵探测器上获得的投影数据进行一次一维傅立叶变换,再与滤波器函数进行卷积运算,得到各个方向卷积滤波后的投影数据;然后把它们沿各个方向进行反投影,即按其原路径平均分配到每一矩阵单元上,进行重叠后得到每一矩阵单元的CT值;再经过适当处理后得到被扫描物体的断层图像- X& e7 u* y5 F" |- }" P: f
1 S, {* W, H  ~" Y; d3 m+ Y
% ^! L& \2 d6 R
算法步骤如下:
1 j# r% J! h1 o( F' d! T' K; l
1 ]0 o& j( V+ ]* q8 W0 H$ S- W
  • 将原始投影进行一次一维傅立叶变换
  • 设计合适的滤波器,在φ_i的角度下将得到原始投影p(x_r,φ_i)进行卷积滤波,得到滤波后的投影。
  • 将滤波后的投影进行反投影,得到满足x_r=r cos⁡((θ - φ_i))方向上的原图像的密度。
  • 将所有反投影进行叠加,得到重建后的投影。
  • 滤波器(滤波函数)和内插函数的选取
    9 l- Z; C8 O. n, Y. e4 m8 ]
      由于直接使用反投影算法会存在两个对实验结果影响很不好的因素:
, Y. q! U# g) [4 \- \5 {( ]% y+ r. B
不准确的数据重建图像就会产生各种伪影。
* \8 q/ X$ M$ p7 e投影的数据是天然离散的,处理不当的话会产生很大的误差。( g8 p3 z" I2 f" W/ ]/ g
常见的滤波器有R-S滤波函数和S-L滤波函数。R-L滤波函数滤波计算简单,避免了大量的正弦、余弦计算,得到的采样序列分是分段现行的,并没有明显的降低图像质量,所以重建图像轮廓清楚,空间分辨率高。  o$ @# L+ {% c

: g4 I, N$ b% s: X6 y3 a" G常见的插值方法有最近邻插值和双线插值,最近邻插值即将离散点中间的缺失值用离它最近的整数处的投影值来替代。
0 [4 t% E9 P- C' M; v
2 ~! q4 I. [1 _! ^' ?# L
2 s5 Y! @6 u' x  M二、源代码
* w& @  Q4 u* U6 Q, P2 P
; S1 M3 D6 o$ O
  • % 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');
    $ g. ^) A) E6 j5 o& q2 Z* X7 S# M
          % f9 u3 n& P% s$ _: R' Z' ^" v

& H1 _* j! m) h8 P* x' [三、运行结果, G9 G: D  q0 w. L2 H- \5 A8 H
% ^) h5 k. `$ v$ K' s9 ^3 `! V/ P

# B$ q8 }* ~6 d  q" D1 k" u7 E

该用户从未签到

2#
发表于 2021-5-13 14:47 | 只看该作者
基于matlab ASTRA算法图像重建

该用户从未签到

3#
发表于 2021-5-13 17:05 | 只看该作者
基于matlab ASTRA算法图像重建

该用户从未签到

4#
发表于 2021-5-13 17:06 | 只看该作者
基于matlab ASTRA算法图像重建
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-10-6 13:15 , Processed in 0.156250 second(s), 26 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表