"""
Tomographic projector backends.
@author: Nicola VIGANÒ, Computational Imaging group, CWI, The Netherlands,
and ESRF - The European Synchrotron, Grenoble, France
"""
from abc import ABC, abstractmethod
from collections.abc import Mapping, Sequence
import numpy as np
import skimage
import skimage.transform as skt
from numpy.typing import ArrayLike, NDArray
from . import filters
from .models import ProjectionGeometry, VolumeGeometry
try:
import astra
has_astra = True
has_cuda = astra.astra.use_cuda()
if not has_cuda:
print(
"WARNING: CUDA is not available. Only 2D operations on CPU are available (scikit-image will be used as default)."
)
except ImportError:
has_astra = False
has_cuda = False
print("WARNING: ASTRA is not available. Only 2D operations on CPU are available (scikit-image will be used).")
if has_cuda:
try:
import astra.experimental
has_astra_direct = True
except ImportError:
has_astra_direct = False
print("WARNING: the experimental ASTRA direct interface is not available. The traditional interface will be used.")
[docs]
def compute_attenuation(vol: NDArray, angle_rad: float, invert: bool = False) -> NDArray:
"""
Compute the attenuation volume for the given local attenuation, and angle.
Parameters
----------
vol : NDArray
The local attenuation volume.
angle_rad : float
The angle along which to compute the attenuation.
invert : bool, optional
Whether to invert propagation direction. The default is False.
Returns
-------
NDArray
The attenuation volume.
"""
def pad_vol(vol: NDArray, edges: Sequence[int]):
paddings = np.zeros((len(vol.shape), 1), dtype=int)
paddings[-2], paddings[-1] = edges[0], edges[1]
return np.pad(vol, paddings, mode="constant")
def compute_cumsum(vol: NDArray, rot_angle_deg: float):
vol = skt.rotate(vol, -rot_angle_deg, order=1, clip=False)
vol += np.roll(vol, 1, axis=-2)
vol = np.cumsum(vol / 2, axis=-2)
return skt.rotate(vol, rot_angle_deg, order=1, clip=False)
size_lims = np.array(vol.shape[-2:])
min_size = np.ceil(np.sqrt(np.sum(size_lims**2)))
edges = np.ceil((min_size - size_lims) / 2).astype(int)
if invert:
angle_rad += np.pi
rot_angle_deg = np.rad2deg(angle_rad)
cum_arr = pad_vol(vol, edges)
if cum_arr.ndim > 2:
prev_shape = np.array(cum_arr.shape, ndmin=1)
num_slices = np.prod(prev_shape[:-2])
cum_arr = cum_arr.reshape([num_slices, *prev_shape[-2:]])
for ii in range(num_slices):
cum_arr[ii] = compute_cumsum(cum_arr[ii], rot_angle_deg)
cum_arr = cum_arr.reshape(prev_shape)
else:
cum_arr = compute_cumsum(cum_arr, rot_angle_deg)
cum_arr = cum_arr[..., edges[0] : -edges[0], edges[1] : -edges[1]]
cum_arr = np.exp(-cum_arr)
return cum_arr
[docs]
class ProjectorBackend(ABC):
"""Base abstract projector backend class. All backends should inherit from this class."""
vol_geom: VolumeGeometry
vol_shape_zxy: NDArray[np.integer]
angles_w_rad: NDArray[np.floating]
prj_shape_vwu: NDArray[np.integer]
prj_shape_vu: NDArray[np.integer]
is_initialized: bool
is_ready: bool
def __init__(self):
"""Initialize base abstract projector backend class."""
self.is_initialized = False
self.is_ready = False
[docs]
def initialize_geometry(
self,
vol_geom: VolumeGeometry,
angles_rot_rad: ArrayLike | NDArray,
rot_axis_shift_pix: ArrayLike | NDArray | None = None,
prj_geom: ProjectionGeometry | None = None,
create_single_projs: bool = False,
):
"""Initialize the projector geometry.
Parameters
----------
vol_geom : VolumeGeometry
The volume geometry.
angles_rot_rad : ArrayLike | NDArray
The projection angles.
rot_axis_shift_pix : float | NDArray | None, optional
Relative position of the rotation center with respect to the volume center. The default is None.
prj_geom : ProjectionGeometry, optional
The fully specified projection geometry.
When active, the rotation axis shift is ignored. The default is None.
create_single_projs : bool, optional
Whether to create projectors for single projections. Used for corrections and SART. The default is False.
"""
self.vol_geom = vol_geom
self.vol_shape_zxy = self.vol_geom.shape_zxy
self.angles_w_rad = np.array(angles_rot_rad, ndmin=1, dtype=np.float32)
self.vol_eff_shape_zxy = np.floor(self.vol_shape_zxy * self.vol_geom.vox_size).astype(int)
# Basic sizes, unless overridden
self.prj_shape_vwu = np.array(
[*self.vol_eff_shape_zxy[:-2], len(self.angles_w_rad), self.vol_eff_shape_zxy[-1]], dtype=int
)
self.prj_shape_vu = np.array([*self.vol_eff_shape_zxy[:-2], 1, self.vol_eff_shape_zxy[-1]], dtype=int)
self.has_individual_projs = create_single_projs
self.is_initialized = True
[docs]
def get_vol_shape(self) -> NDArray:
"""Return the expected and produced volume shape (in ZXY coordinates).
Returns
-------
NDArray
The volume shape.
"""
return self.vol_shape_zxy
[docs]
def get_prj_shape(self) -> NDArray:
"""Return the expected and produced projection shape (in VWU coordinates).
Returns
-------
NDArray
The projection shape.
"""
return self.prj_shape_vwu
[docs]
def make_ready(self) -> None:
"""Initialize the projector.
It should make sure that all the resources have been allocated.
"""
self.is_ready = True
[docs]
def dispose(self) -> None:
"""De-initialize the projector.
It should make sure that all the resources have been de-allocated.
"""
self.is_ready = False
[docs]
def __del__(self):
"""De-initialize projector on deletion."""
if self.is_ready:
self.dispose()
[docs]
def __repr__(self) -> str:
"""Build a string representation of the projector backend.
Returns
-------
str
The representation of the projector.
"""
class_name = f"{self.__class__.__name__}: "
if self.is_initialized:
return (
class_name
+ "{\n"
+ f" Shape vol ZXY: {self.get_vol_shape()}\n"
+ f" Shape prj VWU: {self.get_prj_shape()}\n"
+ f" Angles (deg): {np.rad2deg(self.angles_w_rad)}\n"
+ "}"
)
else:
return class_name + "{ Not initialized! }"
[docs]
@abstractmethod
def fp(self, vol: NDArray, angle_ind: int | None = None) -> NDArray:
"""Forward-project volume.
Forward-projection interface. Derived backends need to implement this method.
Parameters
----------
vol : NDArray
The volume to forward-project.
angle_ind : int, optional
The angle index to foward project. The default is None.
"""
[docs]
@abstractmethod
def bp(self, prj: NDArray, angle_ind: int | None = None) -> NDArray:
"""Back-project data.
Back-projection interface. Derived backends need to implement this method.
Parameters
----------
prj : NDArray
The sinogram to back-project or a single line.
angle_ind : int, optional
The angle index to foward project. The default is None.
"""
[docs]
class ProjectorBackendSKimage(ProjectorBackend):
"""Projector backend based on scikit-image."""
def __init__(self) -> None:
"""Initialize scikit-image projector backend class."""
super().__init__()
self.is_ready = True
[docs]
def initialize_geometry(
self,
vol_geom: VolumeGeometry,
angles_rot_rad: ArrayLike | NDArray,
rot_axis_shift_pix: ArrayLike | NDArray | None = None,
prj_geom: ProjectionGeometry | None = None,
create_single_projs: bool = False,
):
"""Initialize projector backend based on scikit-image.
Parameters
----------
vol_geom : VolumeGeometry
The volume shape.
angles_rot_rad : ArrayLike
The projection angles.
rot_axis_shift_pix : float | NDArray | None, optional
Relative position of the rotation center with respect to the volume center. The default is None.
NOT SUPPORTED: if anything else than None is passed, it will throw an error!
prj_geom : ProjectionGeometry, optional
The fully specified projection geometry.
When active, the rotation axis shift is ignored. The default is None.
NOT SUPPORTED: if anything else than None is passed, it will throw an error!
create_single_projs : bool, optional
Whether to create projectors for single projections. Used for corrections and SART. The default is False.
Raises
------
ValueError
In case the volume dimensionality is larger than 2D, and if a rotation axis shift is passed.
"""
if vol_geom.is_3D():
raise ValueError("With the scikit-image backend, only 2D volumes are allowed!")
if rot_axis_shift_pix is not None:
raise ValueError("With the scikit-image backend, rotation axis shifts are not supported!")
if prj_geom is not None:
raise ValueError("With the scikit-image backend, `ProjectionGeometry` is not supported!")
super().initialize_geometry(vol_geom, angles_rot_rad, create_single_projs=create_single_projs)
self.angles_w_deg = np.rad2deg(self.angles_w_rad)
[docs]
@staticmethod
def _set_filter_name(filt):
if skimage.__version__ >= "0.18":
return dict(filter_name=filt)
else:
return dict(filter=filt)
[docs]
@staticmethod
def _set_bpj_size(output_size):
return dict(circle=False, output_size=output_size)
[docs]
def fp(self, vol: NDArray, angle_ind: int | None = None) -> NDArray:
"""Forward-projection of the volume to the sinogram or a single sinogram line.
Parameters
----------
vol : NDArray
The volume to forward-project.
angle_ind : int, optional
The angle index to foward project. The default is None.
Returns
-------
NDArray
The forward-projected sinogram or sinogram line.
"""
if angle_ind is None:
prj = np.empty(self.prj_shape_vwu, dtype=vol.dtype)
for ii_a, a in enumerate(self.angles_w_deg):
prj[ii_a, :] = np.squeeze(skt.radon(vol, [a]))
return prj
else:
return np.squeeze(skt.radon(vol, self.angles_w_deg[angle_ind : angle_ind + 1 :]))
[docs]
def bp(self, prj: NDArray, angle_ind: int | None = None) -> NDArray:
"""Back-projection of a single sinogram line to the volume.
Parameters
----------
prj : NDArray
The sinogram to back-project or a single line.
angle_ind : int, optional
The angle index to foward project. The default is None.
Returns
-------
NDArray
The back-projected volume.
"""
filter_name = self._set_filter_name(None)
bpj_size = self._set_bpj_size(self.vol_shape_zxy[-1])
if angle_ind is None:
vol = np.empty([self.prj_shape_vwu[-2], *self.vol_shape_zxy], dtype=prj.dtype)
for ii_a, a in enumerate(self.angles_w_deg):
vol[ii_a, ...] = skt.iradon(prj[ii_a, :, np.newaxis], [a], **bpj_size, **filter_name)
vol = vol.sum(axis=0)
else:
vol = skt.iradon(prj[:, np.newaxis], self.angles_w_deg[angle_ind : angle_ind + 1 :], **bpj_size, **filter_name)
return vol * 2 / np.pi
[docs]
class ProjectorBackendASTRA(ProjectorBackend):
"""Projector backend based on astra-toolbox."""
proj_id: list
astra_vol_geom: Mapping
proj_geom_ind: Sequence[Mapping]
proj_geom_all: Mapping
def __init__(self, super_sampling: int = 1):
"""Initialize the ASTRA projector backend.
Parameters
----------
super_sampling : int, optional
Super sampling factor for the pixels and voxels, by default 1.
"""
super().__init__()
self.super_sampling = super_sampling
[docs]
def initialize_geometry(
self,
vol_geom: VolumeGeometry,
angles_rot_rad: ArrayLike | NDArray,
rot_axis_shift_pix: ArrayLike | NDArray | None = None,
prj_geom: ProjectionGeometry | None = None,
create_single_projs: bool = False,
):
"""Initialize geometry of projector backend based on astra-toolbox.
Parameters
----------
vol_geom : VolumeGeometry
The volume shape.
angles_rot_rad : ArrayLike
The projection angles.
rot_axis_shift_pix : float | NDArray | None, optional
Relative position of the rotation center with respect to the volume center. The default is None.
prj_geom : ProjectionGeometry, optional
The fully specified projection geometry.
When active, the rotation axis shift is ignored. The default is None.
create_single_projs : bool, optional
Whether to create projectors for single projections. Used for corrections and SART. The default is False.
Raises
------
ValueError
In case the volume dimensionality is larger than 2D and CUDA is not available.
"""
if vol_geom.is_3D() and not has_cuda:
raise ValueError("CUDA is not available: only 2D volumes are allowed!")
if not (rot_axis_shift_pix is None or isinstance(rot_axis_shift_pix, (int, float, Sequence, np.ndarray))):
raise ValueError(
"Rotation axis shift should either be None or one of the following: int, a float or a sequence of floats"
+ f" ({type(rot_axis_shift_pix)} given instead)."
)
super().initialize_geometry(vol_geom, angles_rot_rad, create_single_projs=create_single_projs)
self.proj_id = []
self.dispose()
num_angles = self.angles_w_rad.size
if self.vol_geom.is_3D():
self.astra_vol_geom = astra.create_vol_geom(*vol_geom.shape_xyz[list([1, 0, 2])], *self.vol_geom.extent)
if prj_geom is None:
prj_geom = ProjectionGeometry.get_default_parallel(geom_type="3d", rot_axis_shift_pix=rot_axis_shift_pix)
if prj_geom.det_shape_vu is None:
prj_geom.det_shape_vu = np.array(self.prj_shape_vwu[list([-3, -1])], dtype=int)
else:
# Here the projections are supposed to be larger or smaller than the sample size
self.prj_shape_vwu = np.array([*prj_geom.det_shape_vu[:-1], num_angles, prj_geom.det_shape_vu[-1]])
self.prj_shape_vu = np.array([*prj_geom.det_shape_vu[:-1], 1, prj_geom.det_shape_vu[-1]])
rot_geom = prj_geom.rotate(self.angles_w_rad)
vectors = np.empty([num_angles, 12])
# source / beam direction
vectors[:, 0:3] = rot_geom.get_field_scaled("src_pos_xyz")
# center of detector
vectors[:, 3:6] = rot_geom.get_field_scaled("det_pos_xyz")
# vector from detector pixel (0, 0) to (0, 1)
vectors[:, 6:9] = rot_geom.get_field_scaled("det_u_xyz")
# vector from detector pixel (0, 0) to (1, 0)
vectors[:, 9:12] = rot_geom.get_field_scaled("det_v_xyz")
geom_type_str = prj_geom.geom_type
else:
self.astra_vol_geom = astra.create_vol_geom(*vol_geom.shape_xyz[list([1, 0])], *self.vol_geom.extent)
if prj_geom is None:
prj_geom = ProjectionGeometry.get_default_parallel(geom_type="2d", rot_axis_shift_pix=rot_axis_shift_pix)
if prj_geom.det_shape_vu is None:
prj_geom.det_shape_vu = np.array(self.prj_shape_vwu[list([-1])], dtype=int)
else:
# Here the projections are supposed to be larger or smaller than the sample size
self.prj_shape_vwu = np.array([num_angles, prj_geom.det_shape_vu[-1]])
self.prj_shape_vu = np.array([1, prj_geom.det_shape_vu[-1]])
rot_geom = prj_geom.rotate(self.angles_w_rad, patch_astra_2d=True)
vectors = np.empty([num_angles, 6])
# source / beam direction
vectors[:, 0:2] = rot_geom.get_field_scaled("src_pos_xyz")
# center of detector
vectors[:, 2:4] = rot_geom.get_field_scaled("det_pos_xyz")
# vector from detector pixel 0 to 1
vectors[:, 4:6] = rot_geom.get_field_scaled("det_u_xyz")
geom_type_str = prj_geom.geom_type[:-2]
vectors /= vol_geom.vox_size
if self.has_individual_projs:
self.proj_geom_ind = [
astra.create_proj_geom(geom_type_str + "_vec", *prj_geom.det_shape_vu, vectors[ii : ii + 1 :, :])
for ii in range(num_angles)
]
self.proj_geom_all = astra.create_proj_geom(geom_type_str + "_vec", *prj_geom.det_shape_vu, vectors)
# def get_vol_shape(self) -> NDArray:
# """Return the expected and produced volume shape (in ZYX coordinates).
# Returns
# -------
# NDArray
# The volume shape.
# """
# return astra.functions.geom_size(self.astra_vol_geom)
# def get_prj_shape(self) -> NDArray:
# """Return the expected and produced projection shape (in VWU coordinates).
# Returns
# -------
# NDArray
# The projection shape.
# """
# return astra.functions.geom_size(self.proj_geom_all)
[docs]
def make_ready(self) -> None:
"""Initialize the ASTRA projectors."""
if not self.is_ready:
if self.vol_geom.is_3D():
projector_type = "cuda3d"
self.algo_type = "3D_CUDA"
self.data_mod = astra.data3d
else:
if has_cuda:
projector_type = "cuda"
self.algo_type = "_CUDA"
else:
projector_type = "linear"
self.algo_type = ""
self.data_mod = astra.data2d
voxel_sampling = int(self.super_sampling * np.fmax(1, self.vol_geom.vox_size))
pixel_sampling = int(self.super_sampling / np.fmin(1, self.vol_geom.vox_size))
opts = {"VoxelSuperSampling": voxel_sampling, "DetectorSuperSampling": pixel_sampling}
if self.has_individual_projs:
self.proj_id = [
astra.create_projector(projector_type, pg, self.astra_vol_geom, opts) for pg in self.proj_geom_ind
]
self.proj_id.append(astra.create_projector(projector_type, self.proj_geom_all, self.astra_vol_geom, opts))
super().make_ready()
[docs]
def _check_data(self, x: NDArray, expected_shape: Sequence[int] | NDArray[np.integer]) -> NDArray:
if x.dtype != np.float32:
x = x.astype(np.float32)
if not x.flags["C_CONTIGUOUS"]:
x = np.ascontiguousarray(x)
try:
return x.reshape(np.array(expected_shape, ndmin=1))
except ValueError:
print(f"Could not reshape input data of shape={x.shape} into expected shape={expected_shape}")
raise
[docs]
def _check_prj_shape(self, prj: NDArray) -> None:
inds_vu = [ind for ind in range(-len(self.prj_shape_vu), 0) if self.prj_shape_vu[ind] > 1]
expected_prj_shape_vu = self.prj_shape_vu[inds_vu]
given_prj_shape_vu = np.array(prj.shape)[inds_vu]
if np.any(given_prj_shape_vu != expected_prj_shape_vu):
raise ValueError(
f"Expected prj shape (VU: {expected_prj_shape_vu} - likely inferred from `vol_geom`),"
f" and passed prj shape (VU: {given_prj_shape_vu}) from the projection data, do not match.\n"
"If you intend to use projection shapes different from the volume shape, please adjust"
" the `prj_geom` object's `det_shape_vu` accordingly."
)
[docs]
def dispose(self) -> None:
"""De-initialize the ASTRA projectors."""
for p in self.proj_id:
astra.projector.delete(p)
self.proj_id = []
super().dispose()
[docs]
def fp(self, vol: NDArray, angle_ind: int | None = None) -> NDArray:
"""Apply forward-projection of the volume to the sinogram or a single sinogram line.
Parameters
----------
vol : NDArray
The volume to forward-project.
angle_ind : int | None, optional
The angle index to forward project. The default is None.
Returns
-------
NDArray
The forward-projected sinogram or sinogram line.
"""
self.make_ready()
vol = self._check_data(vol, self.vol_shape_zxy)
if angle_ind is None:
prj = np.empty(self.prj_shape_vwu, dtype=np.float32)
prj_geom = self.proj_geom_all
proj_id = self.proj_id[-1]
else:
if not self.has_individual_projs:
raise ValueError("Individual projectors not available!")
prj = np.empty(self.prj_shape_vu, dtype=np.float32)
prj_geom = self.proj_geom_ind[angle_ind]
proj_id = self.proj_id[angle_ind]
vid = self.data_mod.link("-vol", self.astra_vol_geom, vol)
try:
sid = self.data_mod.link("-sino", prj_geom, prj)
try:
cfg = astra.creators.astra_dict("FP" + self.algo_type)
cfg["ProjectionDataId"] = sid
cfg["VolumeDataId"] = vid
cfg["ProjectorId"] = proj_id
fp_id = astra.algorithm.create(cfg)
try:
astra.algorithm.run(fp_id)
finally:
astra.algorithm.delete(fp_id)
finally:
self.data_mod.delete([sid])
finally:
self.data_mod.delete([vid])
return np.squeeze(prj)
[docs]
def bp(self, prj: NDArray, angle_ind: int | None = None) -> NDArray:
"""Apply back-projection of a single sinogram line to the volume.
Parameters
----------
prj : NDArray
The sinogram to back-project or a single line.
angle_ind : int | None, optional
The angle index to forward project. The default is None.
Returns
-------
NDArray
The back-projected volume.
"""
self.make_ready()
self._check_prj_shape(prj)
vol = np.empty(self.vol_shape_zxy, dtype=np.float32)
if angle_ind is None:
prj = self._check_data(prj, self.prj_shape_vwu)
prj_geom = self.proj_geom_all
proj_id = self.proj_id[-1]
else:
if not self.has_individual_projs:
raise ValueError("Individual projectors not available!")
prj = self._check_data(prj, self.prj_shape_vu)
prj_geom = self.proj_geom_ind[angle_ind]
proj_id = self.proj_id[angle_ind]
vid = self.data_mod.link("-vol", self.astra_vol_geom, vol)
try:
sid = self.data_mod.link("-sino", prj_geom, prj)
try:
cfg = astra.creators.astra_dict("BP" + self.algo_type)
cfg["ProjectionDataId"] = sid
cfg["ReconstructionDataId"] = vid
cfg["ProjectorId"] = proj_id
bp_id = astra.algorithm.create(cfg)
try:
astra.algorithm.run(bp_id)
finally:
astra.algorithm.delete(bp_id)
finally:
self.data_mod.delete([sid])
finally:
self.data_mod.delete([vid])
return vol
[docs]
class ProjectorBackendDirectASTRA(ProjectorBackendASTRA):
"""Experimental astra-toolbox functions projector."""
astra_vol_shape: Sequence
astra_prj_shape: Sequence
astra_angle_prj_shape: Sequence
angle_prj_shape: Sequence
[docs]
def initialize_geometry(
self,
vol_geom: VolumeGeometry,
angles_rot_rad: ArrayLike | NDArray,
rot_axis_shift_pix: ArrayLike | NDArray | None = None,
prj_geom: ProjectionGeometry | None = None,
create_single_projs: bool = False,
):
"""Initialize projector backend based on experimental astra-toolbox functions.
Parameters
----------
vol_geom : VolumeGeometry
The volume shape.
angles_rot_rad : ArrayLike
The projection angles.
rot_axis_shift_pix : float, optional
Relative position of the rotation center with respect to the volume center. The default is 0.0.
geom : ProjectionGeometry, optional
The fully specified projection geometry.
When active, the rotation axis shift is ignored. The default is None.
create_single_projs : bool, optional
Whether to create projectors for single projections. Used for corrections and SART. The default is False.
Raises
------
ValueError
In case CUDA is not available.
"""
if not has_cuda:
raise ValueError("CUDA is not available, but it is required for the direct functions!")
if not has_astra_direct:
raise ValueError("ASTRA direct is not available, but it is required!")
if not (rot_axis_shift_pix is None or isinstance(rot_axis_shift_pix, (int, float, Sequence, np.ndarray))):
raise ValueError(
"Rotation axis shift should either be None or one of the following: int, a float or a sequence of floats"
+ f" ({type(rot_axis_shift_pix)} given instead)."
)
ProjectorBackend.initialize_geometry(self, vol_geom, angles_rot_rad, create_single_projs=create_single_projs)
self.proj_id = []
self.dispose()
num_angles = self.angles_w_rad.size
vol_geom_tmp3d = self.vol_geom.get_3d()
self.astra_vol_geom = astra.create_vol_geom(*vol_geom_tmp3d.shape_xyz[list([1, 0, 2])], *vol_geom_tmp3d.extent)
if prj_geom is None:
prj_geom = ProjectionGeometry.get_default_parallel(geom_type="3d", rot_axis_shift_pix=rot_axis_shift_pix)
else:
if prj_geom.det_shape_vu is not None:
# Here the projections are supposed to be larger or smaller than the sample size
# We use the original version because it makes sure that 2D is respected.
self.prj_shape_vwu = np.array([*prj_geom.det_shape_vu[:-1], num_angles, prj_geom.det_shape_vu[-1]])
self.prj_shape_vu = np.array([*prj_geom.det_shape_vu[:-1], 1, prj_geom.det_shape_vu[-1]])
prj_geom = prj_geom.get_3d()
if prj_geom.det_shape_vu is None:
prj_shape_vu = np.delete(self.prj_shape_vwu, obj=-2)
prj_geom.det_shape_vu = np.ones(2, dtype=int)
prj_geom.det_shape_vu[-len(prj_shape_vu) :] = prj_shape_vu
rot_geom = prj_geom.rotate(self.angles_w_rad)
vectors = np.empty([num_angles, 12])
# source / beam direction
vectors[:, 0:3] = rot_geom.get_field_scaled("src_pos_xyz")
# center of detector
vectors[:, 3:6] = rot_geom.get_field_scaled("det_pos_xyz")
# vector from detector pixel (0, 0) to (0, 1)
vectors[:, 6:9] = rot_geom.get_field_scaled("det_u_xyz")
# vector from detector pixel (0, 0) to (1, 0)
vectors[:, 9:12] = rot_geom.get_field_scaled("det_v_xyz")
geom_type_str = prj_geom.geom_type
if self.has_individual_projs:
self.proj_geom_ind = [
astra.create_proj_geom(geom_type_str + "_vec", *prj_geom.det_shape_vu, vectors[ii : ii + 1 :, :])
for ii in range(num_angles)
]
self.proj_geom_all = astra.create_proj_geom(geom_type_str + "_vec", *prj_geom.det_shape_vu, vectors)
self.astra_vol_shape = tuple(vol_geom_tmp3d.shape_zxy)
self.astra_prj_shape = (prj_geom.det_shape_vu[-2], num_angles, prj_geom.det_shape_vu[-1])
self.astra_angle_prj_shape = (prj_geom.det_shape_vu[-2], 1, prj_geom.det_shape_vu[-1])
self.angle_prj_shape = (prj_geom.det_shape_vu[-2], prj_geom.det_shape_vu[-1])
[docs]
def make_ready(self):
"""Initialize the ASTRA projectors."""
if not self.is_ready:
voxel_sampling = int(self.super_sampling * np.fmax(1, self.vol_geom.vox_size))
pixel_sampling = int(self.super_sampling / np.fmin(1, self.vol_geom.vox_size))
opts = {"VoxelSuperSampling": voxel_sampling, "DetectorSuperSampling": pixel_sampling}
if self.has_individual_projs:
self.proj_id = [astra.create_projector("cuda3d", pg, self.astra_vol_geom, opts) for pg in self.proj_geom_ind]
self.proj_id.append(astra.create_projector("cuda3d", self.proj_geom_all, self.astra_vol_geom, opts))
ProjectorBackend.make_ready(self)
[docs]
def fp(self, vol: NDArray, angle_ind: int | None = None):
"""Apply forward-projection of the volume to the sinogram or a single sinogram line.
Parameters
----------
vol : NDArray
The volume to forward-project.
angle_ind : int | None, optional
The angle index to forward project. The default is None.
Returns
-------
NDArray
The forward-projected sinogram or sinogram line.
"""
self.make_ready()
if angle_ind is None:
prj = np.zeros(self.astra_prj_shape, dtype=np.float32)
proj_id = self.proj_id[-1]
out_shape = self.prj_shape_vwu
else:
if not self.has_individual_projs:
raise ValueError("Individual projectors not available!")
prj = np.zeros(self.astra_angle_prj_shape, dtype=np.float32)
proj_id = self.proj_id[angle_ind]
out_shape = self.angle_prj_shape
vol = self._check_data(vol, self.astra_vol_shape)
astra.experimental.direct_FP3D(proj_id, vol, prj)
return prj.reshape(out_shape)
[docs]
def bp(self, prj: NDArray, angle_ind: int | None = None):
"""Apply back-projection of a single sinogram line to the volume.
Parameters
----------
prj : NDArray
The sinogram to back-project or a single line.
angle_ind : int | None, optional
The angle index to forward project. The default is None.
Returns
-------
NDArray
The back-projected volume.
"""
self.make_ready()
self._check_prj_shape(prj)
if angle_ind is None:
prj = self._check_data(prj, self.astra_prj_shape)
proj_id = self.proj_id[-1]
else:
if not self.has_individual_projs:
raise ValueError("Individual projectors not available!")
prj = self._check_data(prj, self.astra_angle_prj_shape)
proj_id = self.proj_id[angle_ind]
vol = np.zeros(self.astra_vol_shape, dtype=np.float32)
astra.experimental.direct_BP3D(proj_id, vol, prj)
return vol.reshape(self.vol_shape_zxy)