ASTRA Toolbox学习笔记-Python

GitHub网址:https://github.com/astra-toolbox/astra-toolbox
官方网址:http://www.astra-toolbox.com/
文章网址:Fast and flexible X-ray tomography using the ASTRA toolbox



一、简介

ASTRA Toolbox 是一个用于处理 X 射线相干散射和 X 射线传输成像的开源算法工具箱。本文档概述了该工具箱的基本概念。

几何

ASTRA Toolbox 支持以下几种类型的投影几何:

  • 平行束二维投影
  • 锥束二维投影
  • 平行束三维投影
  • 锥束三维投影
  • 向量二维投影
  • 向量三维投影

每种几何类型都可以通过 ASTRA 提供的 create_*_proj_geom 函数创建。

数据和算法

ASTRA Toolbox 为 X 射线相干散射和 X 射线传输成像提供了许多不同的算法,包括以下几种:

  • Filtered Backprojection (FBP)
  • Simultaneous Algebraic ReconstructionTechnique (SART)
  • Simultaneous Iterative Reconstruction Technique (SIRT)
  • Conjugate Gradient (CG)
  • Total Variation (TV)

ASTRA Toolbox 提供了创建和管理算法对象的 API,用于指定算法的输入和输出数据以及算法的参数。有关如何使用 API 的详细信息,请参阅 ASTRA Toolbox 文档。

二、使用

1.创建几何(Geometry)

首先需要制定实验所需的几何形状,包括两部分:
(1)Volume geometries(体积几何)
指定了样品或幻影在二维或三维空间中的位置。它采用2D矩形或3D盒子的形式,通常以原点为中心。

创建一个3D的体积几何:

vol_geom = astra_create_vol_geom(y,x,z); 

astra_create_vol_geom(y, x, z)用于创建一个 3D 体积几何。它接收三个参数,分别表示体积在三个维度(y, x, z)上的大小。这个函数返回一个用于描述体积几何的对象,这个对象将在之后的投影和重建操作中使用。

例如,如果你想创建一个具有尺寸 128x128x128 的立方体,你可以使用以下代码:

vol_geom = astra_create_vol_geom(128, 128, 128)

这将创建一个表示 128x128x128 大小的 3D 体积几何的对象。在之后的投影和重建操作中,这个对象将被用于描述样品或模拟样品 的位置和尺寸。

(2)a projection geometry(投影几何)
指定了射线源和探测器的位置和轨迹。可用的投影几何类型有2D平行光束、2D扇形光束、3D平行光束和3D锥形光束。

proj_geom = astra_create_proj_geom('cone',  det_spacing_x, det_spacing_y, det_row_count, det_col_count, angles, source_origin, origin_det);

创建一个圆锥投影几何(cone)。函数的参数如下:

‘cone’:代表几何类型,这里是圆锥投影。
det_spacing_x:探测器像素在x方向上的间距。
det_spacing_y:探测器像素在y方向上的间距。
det_row_count:探测器面板的行数。
det_col_count:探测器面板的列数。
angles:投影角度的数组,单位通常是弧度。
source_origin:源点到物体原点的距离。
origin_det:物体原点到探测器的距离。

函数返回一个代表投影几何的对象,它包含了上述参数的信息。这个投影几何对象可用于后续的投影、反投影等操作。

2.创建投影数据对象(data)

使用astra.data3d模块来创建3D数据对象。主要使用的函数为以下两个:astra.data3d.create(data_type, geometry, data=None)astra.data3d.delete(data_object_ids),用于创建和删除数据对象。以下是该模块中的主要函数及其参数和用法:


1. astra.data3d.create(data_type, geometry, data=None):用于创建一个ASTRA数据对象。这个函数可以用于创建三维的投影数据或者体积数据。

  • data_type 是一个字符串,表示数据类型。对于三维投影数据,使用 ‘-proj3d’。对于三维体积数据,使用 ‘-vol’。
  • geometry 是一个包含投影几何或体积几何的字典。对于投影数据,使用由 astra.create_proj_geom() 创建的投影几何。对于体积数据,使用由 astra.create_vol_geom() 创建的体积几何。
  • data 是可选的 NumPy 数组,包含初始数据。如果未指定,将使用零填充数据。

返回值:返回一个数据对象ID,这个ID可以用于后续的ASTRA函数调用。

sino_id = astra.data3d.create('-sino', proj_geom, sino_data)
vol_id = astra.data3d.create('-vol', vol_geom, vol_data)

2. astra.data3d.get(dataset_id): 用于获取三维数据集对象的信息。具体来说,可以获取三维数据集对象的属性,例如数据类型、大小、数据范围等。
参数:

  • dataset_id: 数据对象的ID。

返回值:一个包含数据对象内容的NumPy数组。

import numpy as np
import astra

# 创建一个体积几何
vol_geom = astra.create_vol_geom(128, 128, 128)

# 使用现有的 NumPy 数组创建一个三维体积数据对象
data = np.random.rand(128, 128, 128)
vol_data_id = astra.data3d.create('-vol', vol_geom, data)

# 获取数据对象中的数据
retrieved_data = astra.data3d.get(vol_data_id)

# 检查原始数据和检索到的数据是否相同
print(np.allclose(data, retrieved_data))  # 如果数据相同,则输出 True

3. astra.data3d.store(data_id, data): 将数据对象的内容替换为新数据。此函数在ASTRA Toolbox 1.9.9及更高版本中已被弃用,推荐使用astra.data3d.set。

  • data_id:要存储数据的 ASTRA 数据对象的 ID。
  • data:一个 NumPy 数组,包含要存储到指定数据对象的数据。对于投影数据对象,它应该具有形状(det_row_count, num_angles, det_col_count);对于体积数据对象,它应该具有形状 (GridRowCount, GridColCount, GridSliceCount)。
  • 此函数没有返回值。在调用此函数之后,指定的数据对象将包含提供的数据。

4. astra.data3d.delete(data_object_ids):用于删除在 ASTRA Toolbox 中创建的一个或多个 3D 数据对象。此函数可以释放这些数据对象所占用的内存。当不再需要这些数据对象时,通常使用此函数进行清理。
参数:

  • data_object_ids:一个整数或整数列表,表示要删除的一个或多个 3D 数据对象的 ID。
import astra

# 创建一个投影几何对象
proj_geom = astra.create_proj_geom(...)

# 创建一个体积几何对象
vol_geom = astra.create_vol_geom(...)

# 创建投影数据对象
proj_id = astra.data3d.create('-proj3d', proj_geom)

# 创建体积数据对象
vol_id = astra.data3d.create('-vol', vol_geom)

# 删除投影数据对象和体积数据对象
astra.data3d.delete([proj_id, vol_id])

5. astra.data3d.info: 获取数据对象的信息,如数据类型、尺寸等。
参数:

  • id: 数据对象的ID。

返回值:一个字典,包含数据对象的信息。

sino_info = astra.data3d.info(sino_id)
vol_info = astra.data3d.info(vol_id)

这些函数为处理3D数据提供了方便的接口。在使用它们时,请确保在操作完成后删除创建的对象,以避免内存泄漏。

3.创建算法对象(Algorithm)

astra.algorithm.create函数用于在ASTRA Toolbox中创建算法对象。该函数需要一个包含算法配置信息的字典作为输入参数。通常,通过astra.astra_dict类创建这个字典。
以下是使用astra.algorithm.create函数创建SIRT(Simultaneous Iterative Reconstruction Technique)算法对象的示例:

import astra

# 创建astra_dict对象以存储算法配置信息
cfg = astra.astra_dict("SIRT")
cfg["ProjectionDataId"] = sino_id  # 投影数据对象ID
cfg["ReconstructionDataId"] = rec_id  # 重建数据对象ID(实际上是体数据对象ID)
cfg["ProjectorId"] = proj_id  # 投影仪对象ID
cfg["option"] = {"MinConstraint": 0}  # 设置其他选项,例如最小约束

# 使用astra.algorithm.create函数创建算法对象
alg_id = astra.algorithm.create(cfg)

  1. “type”(必需):算法的类型。例如,“SIRT”(Simultaneous Iterative Reconstruction
    Technique)、“FBP”(Filtered
    Backprojection)等。通过astra.astra_dict函数传递算法名称时,该键会自动设置。
  2. “ProjectionDataId”(必需):投影数据对象ID。这是一个整数,表示在ASTRA Toolbox中创建的投影数据对象。
  3. “ReconstructionDataId”(必需):重建数据对象ID。这是一个整数,表示在ASTRA
    Toolbox中创建的体数据对象,用于存储重建结果。
  4. “ProjectorId”(可选):投影仪对象ID。这是一个整数,表示在ASTRA
    Toolbox中创建的投影仪对象。对于某些算法,例如FBP,需要提供投影仪对象。投影仪对象可以通过astra.create_projector函数创建。
  5. “option”(可选):一个包含算法特定选项的字典。例如,对于SIRT算法,您可以设置最小约束和最大约束:
cfg["option"] = {"MinConstraint": 0, "MaxConstraint": 1}

对于FBP算法,您可以设置滤波器类型和滤波器参数:

cfg["option"] = {"FilterType": "Ram-Lak", "FilterParameter": 0.5}
  1. “DetectorSuperSampling”(可选):一个整数,表示检测器的超采样因子。用于提高投影和反投影操作的精度。仅适用于GPU算法。
  2. “VolumeSuperSampling”(可选):一个整数,表示体数据的超采样因子。用于提高投影和反投影操作的精度。仅适用于GPU算法。

三、使用示例

1.通过模拟的几何图形产生投影

模拟正方体模型:
在这里插入图片描述
投影数据:
在这里插入图片描述

Python程序进行投影:

import numpy as np
import astra
import matplotlib.pyplot as plt
import os
from skimage.io import imread

# 创建一个简单的3D幻灯片(立方体)
vol_geom = astra.create_vol_geom(128, 128, 128)
phantom = np.zeros(vol_geom['GridRowCount'] * vol_geom['GridColCount'] * vol_geom['GridSliceCount'], dtype='float32')
phantom = phantom.reshape(vol_geom['GridRowCount'], vol_geom['GridColCount'], vol_geom['GridSliceCount'])
phantom[30:100, 30:100, 30:100] = 1

# 创建一个圆锥光投影几何
angles = np.linspace(0, 2 * np.pi, 360, endpoint=False)
det_col_count = 200  # 列数
det_row_count = 200  # 行数---->最外层
source_origin = 200
origin_det = 200
proj_geom = astra.create_proj_geom('cone', 1, 1, det_row_count, det_col_count, angles, source_origin, origin_det)

# 创建一个投影数据对象
proj_id = astra.data3d.create('-proj3d', proj_geom)

# 计算投影数据
vol_id = astra.data3d.create('-vol', vol_geom, phantom)


alg_cfg = astra.astra_dict('FP3D_CUDA')
alg_cfg['ProjectionDataId'] = proj_id
alg_cfg['VolumeDataId'] = vol_id
alg_id = astra.algorithm.create(alg_cfg)
astra.algorithm.run(alg_id)

# 获取投影数据
projections = astra.data3d.get(proj_id)
# 可视化投影数据
plt.imshow(projections[1], cmap='gray')
plt.colorbar()
plt.show()
# *********************************************保存为jpg图片数据
# output_dir = 'projections'
# if not os.path.exists(output_dir):
#     os.makedirs(output_dir)

# # 归一化投影数据
# normalized_projections = projections / np.max(projections)
#
# # 将每个投影保存为JPG文件
# for i, projection in enumerate(normalized_projections):
#     output_file = os.path.join(output_dir, f'projection_{i:03d}.jpg')
#     plt.imsave(output_file, projection, cmap='gray')

#  ***************************************保存为二进制文件
output_file = 'projections.npy'
np.save(output_file, projections)


# 清理
astra.algorithm.delete(alg_id)
astra.data3d.delete(proj_id)
astra.data3d.delete(vol_id)

2.通过通过投影数据重建图像

重建切片图:
在这里插入图片描述
重建3D图:
在这里插入图片描述

Python程序进行重建:

import numpy as np
import astra
import os
import matplotlib.pyplot as plt
from skimage.io import imread

# 读取保存的投影数据
output_dir = 'projections'
num_projections = 360
det_row_count = 200
det_col_count = 200

# 从jpg获取投影数据
# projections = np.zeros((det_row_count, num_projections, det_col_count), dtype=np.float32)
#
# for i in range(det_row_count):
#     input_file = os.path.join(output_dir, f'projection_{i:03d}.jpg')
#     projection = imread(input_file, as_gray=True)
#     projections[i] = projection
#
# # 将投影数据恢复到原始范围
# projections *= np.max(projections)

# 获取二进制投影数据
input_file = 'projections.npy'
projections = np.load(input_file)
print(projections.shape)
# 投影几何
angles = np.linspace(0, 2 * np.pi, num_projections, endpoint=False)
source_origin = 200
origin_det = 200
proj_geom = astra.create_proj_geom('cone', 1, 1, det_row_count, det_col_count, angles, source_origin, origin_det)

# 体积几何
vol_geom = astra.create_vol_geom(128, 128, 128)


# 创建投影数据和重建数据对象
proj_id = astra.data3d.create('-proj3d', proj_geom, projections)
recon_id = astra.data3d.create('-vol', vol_geom)

# 配置并执行SIRT算法
# cfg = astra.astra_dict('SIRT3D_CUDA')
cfg = astra.astra_dict('FDK_CUDA')

cfg['ProjectionDataId'] = proj_id
cfg['ReconstructionDataId'] = recon_id
alg_id = astra.algorithm.create(cfg)
astra.algorithm.run(alg_id, 1)  # 迭代次数

# 获取并显示重建结果
reconstruction = astra.data3d.get(recon_id)
plt.imshow(reconstruction[64], cmap='gray')  # 显示中间切片
plt.colorbar()
plt.show()

# 清理
astra.algorithm.delete(alg_id)
astra.data3d.delete(proj_id)
astra.data3d.delete(recon_id)

四、更新中… …


*注意

1.astra.create_sino3d_gpu()和astra.create_backprojection3d_gpu()

astra.create_sino3d_gpu()astra.create_backprojection3d_gpu()是ASTRA工具箱Python接口中的两个函数,用于计算三维投影和反投影。这两个函数都利用了GPU(CUDA)加速计算。

  1. astra.create_sino3d_gpu():该函数用于计算三维体积数据的投影数据(正投影)。它接受一个三维数据(体积)、一个投影几何和一个体积几何。函数返回一个包含投影数据的NumPy数组。
proj_id, projections = astra.create_sino3d_gpu(data=phantom, proj_geom=proj_geom, vol_geom=vol_geom)
  1. astra.create_backprojection3d_gpu():该函数用于计算投影数据的反投影(即重建过程的一部分)。它接受投影数据、投影几何和体积几何。函数返回一个包含重建体积数据的NumPy数组。
rec_id, reconstruction = astra.create_backprojection3d_gpu(data=projections, proj_geom=proj_geom, vol_geom=vol_geom)

2.与astra.algorithm.create(cfg)的主要区别在于它们的用途和执行方式。

astra.algorithm.create(cfg):

  • 这个函数用于创建基于配置字典(cfg)的ASTRA算法对象。配置字典指定了要执行的操作(如正投影、反投影或重建算法)以及相关参数。
  • 创建算法对象后,需要使用astra.algorithm.run(alg_id)来执行算法。执行完成后,可以通过astra.data2d.get()astra.data3d.get()获取结果。
  • 这种方法允许更灵活地配置和执行各种操作,如多个迭代的重建算法。

astra.create_sino3d_gpu()和astra.create_backprojection3d_gpu():

  • 这两个函数专门用于快速计算三维投影和反投影,它们直接返回结果而无需创建和执行算法对象。
  • 它们利用了GPU加速,并为特定任务提供了简化的接口。
  • 这些函数更适用于直接进行正投影和反投影操作,而不涉及复杂的重建算法。

总之,astra.algorithm.create(cfg)astra.algorithm.run(alg_id)提供了一种更通用且灵活的方式来执行各种操作(包括正投影、反投影和重建算法),而astra.create_sino3d_gpu()astra.create_backprojection3d_gpu()则提供了一种快速且简便的方法来执行三维投影和反投影操作。

本图文内容来源于网友网络收集整理提供,作为学习参考使用,版权属于原作者。
THE END
分享
二维码
< <上一篇
下一篇>>