Tomosipo: fast, flexible, and convenient 3D tomography for complex scanning geometries in Python

: Tomography is a powerful tool for reconstructing the interior of an object from a series of projection images. Typically, the source and detector traverse a standard path (e


Introduction
Tomographic imaging enables the examination of the internal structure of an object. The object is typically placed between a source and detector, and its structure is reconstructed using projection images from a range of different positions. Collectively, the position information of the source, object, and detector determine the acquisition geometry.
Most common tomographic techniques rely on a selection of standard acquisition geometries, such as circular cone beam or single-axis parallel beam [1]. In recent years, several scientific and industrial applications have emerged whose needs are not met by the standard selection of paths. Such scientific applications include diffraction contrast tomography (DCT) [2] and X-ray scattering tensor tomography (XSTT) [3]. These techniques measure X-ray effects other than absorption, which necessarily give rise to more complex acquisition geometries. Complex geometries also arise in industrial applications like automotive and aerospace testing [4,5], as objects may be too large to fit in conventional scanners. Instead, a robot arm moves the source and detector along an irregular path around the object.
Efficient reconstruction algorithms exist for many common acquisition geometries [1,6]. Such filtered backprojection (FBP)-type algorithms are typically fast to compute [7], but require the source and detector to follow a regular path. Algorithms that permit flexible acquisition geometries, such as SIRT [8] and total variation minimization (TV-MIN) [9], typically follow an iterative reconstruction scheme. As iterative algorithms tend to be more computationally demanding than FBP-type algorithms, they benefit more from an efficient implementation.

Standard tomography problem
Common tomography setups expose a sample to a beam of high energy particles, e.g., photons, electrons, or neutrons, which are collected on a detector. Contrast in the measured projection images is generated by differences in attenuation, refraction index or scattering of the object (e.g. phase and diffraction, respectively), or the emission of secondary signals (e.g. X-ray fluorescence, Compton, Auger). Many of these problems can be modeled as a collection of line integrals through space where the ith measurement y i ∈ R is obtained as a line integral through a point s i ∈ R 3 with direction η i ∈ R 3 . The canonical case is absorption contrast tomography, which we describe here.
In absorption contrast tomography, the reconstruction problem can be posed as a linear discrete inverse problem. Suppose measurements y ∈ R N θ ×N 2 p are acquired from N θ positions using a square detector that is divided into N 2 p pixels. Define the cubic reconstruction volume x ∈ R N 3 v on a voxel grid and let A denote the projection matrix such that A ij describes the absorption by object voxel j of the ray to measurement i. The goal is to determine the value of x that gave rise to the measurement Ax = y.
The computation of the linear operator A depends strongly on the geometry of the acquisition. This includes the direction of the rays, the position and orientation of the reconstruction volume, and the position and orientation of the detector. In tomosipo, the operator A is not stored as a matrix. Instead, the GPU implementation provided by the ASTRA Toolbox uses ray-tracing based techniques.

Framework concepts
Three concepts are essential to the tomosipo package. These are geometries, geometric transformations, and the projection operator A. Geometries represent the position of the source, sample, and detector at each time step. The sample's position and orientation is represented by a volume geometry, and the X-ray source and flat panel detector are represented by a projection geometry, which can model both point sources (cone beam geometry) and parallel box beams (parallel beam geometry). All geometries have two representations: a simple representation that defines a standard trajectory, and a flexible representation that permits arbitrary movement and orientation. Volume and projection geometries are discussed in Section 3.2. Geometries can be manipulated using geometric transforms, as well as split and joined using subsampling and concatenation. In this way, complicated acquisition geometries can be assembled from simple geometries. This is described in Sections 3.3 and 3.5.
Together, a volume and projection geometry define the projection operator A. In tomosipo, the computation using A is GPU-accelerated using the ASTRA Toolbox. Most tomosipo geometries have an ASTRA counterpart, except for the flexible volume geometry whose movement and orientation is compensated for by exploiting the flexibility of ASTRA's projection geometries. The creation of projection operators is discussed in the next section, and the integration with the ASTRA Toolbox and Python array libraries in Section 3.4. Finally, visualization of geometries is discussed in Section 3.6. The main concepts of tomosipo and their relation to the ASTRA Toolbox and the physical geometry are summarized in Fig. 1.

Tomographic projection
In this section, we describe the creation and use of the projection operator A. Tomosipo provides a convenient representation of the projection operator A, offering an API that is similar to the Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.
In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator. Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.  Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.

Tomographic projection
In this section, we describe the creation and use of the projection operator A. Tomosipo provides a convenient representation of the projection operator A, offering an API that is similar to the OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Equation (2) can be obtained as follows: The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Equation (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps. In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator. Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.
OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Eq. (2) can be obtained as follows: The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Eq. (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows:  Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.

Tomographic projection
In this section, we describe the creation and use of the projection operator A. Tomosipo provides a convenient representation of the projection operator A, offering an API that is similar to the OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Equation (2)  The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: 5 x = np.ones(A.domain_shape, dtype=np.float32) # 32-bit floating point array The volume, containing all ones, can be forward projected and the result backprojected with: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Equation (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps. In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator. Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.
OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Eq. (2)  The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Eq. (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are  Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.

Tomographic projection
In this section, we describe the creation and use of the projection operator A. Tomosipo provides a convenient representation of the projection operator A, offering an API that is similar to the OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Equation (2)  The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Equation (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps. In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator. Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.
OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Eq. (2)  The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Eq. (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Eq. (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps   154  155  156  157   158  159  160  161  162   163  164  165  166  167  168   169  170  171  172  173  174   175  176  177  178  179  180   181  182  183  184  185   186  187  188  189  190  191   192  193  194  195  196  197   198  199  200  201  202  203 204  Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.

Tomographic projection
In this section, we describe the creation and use of the projection operator A. Tomosipo provides a convenient representation of the projection operator A, offering an API that is similar to the OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Equation (2)  The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Equation (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps. In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator. Tomosipo can be roughly divided in actions and concepts. The concepts describe the acquisition geometry and the X-ray projection and can be directly mapped onto ASTRA primitives, except for the flexible (moving) volume, which has no ASTRA counterpart. The actions provide the means to transform, recombine, and visualize tomosipo's geometry primitives.
OpTomo Spot operator in the ASTRA Toolbox [23]. Given a volume geometry vg and projection geometry pg, the linear operator A from Eq. (2) can be obtained as follows: The operator A is a stand-alone object. It has domain_shape and range_shape properties that facilitate the creation of data of the correct number of voxels and pixels in each dimension in its mathematical domain and range, i.e., image space and sinogram space. An array representing the volume is created as follows: The computation is performed on the GPU, and is handled by the ASTRA Toolbox. The operator A can be used to solve the inverse problem posed in Eq. (2). In the code below, this is demonstrated by computing a simple Landweber iteration [26] with step size eta in n steps In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are In the next section, we describe how to define the volume and projection geometries that are required to create a projection operator.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are defined in the ASTRA Toolbox as well, except for the volume vector geometry that can represent an arbitrarily oriented moving reconstruction grid. The geometries are illustrated in Fig. 2 with accompanying code.  Fig. 2. Creation of typical tomographic geometries. From left to right: a volume geometry, single-axis parallel beam geometry, and a circular cone beam geometry. Below, arbitrarily oriented vector geometries are shown. The parameters are specified using keyword-only arguments [27]. The pos parameter, for instance, determines the position of a volume, other parameters have accompanying labels in the diagrams. The diagrams are created using tomosipo.

Acquisition geometry primitives
Tomosipo provides three standard geometries: the fixed volume geometry, the single-axis parallel beam geometry, and the circular cone beam geometry. In addition, these geometries have a flexible counterpart that permits arbitrary orientation and movement. The flexible geometries are known as vector geometries, following the terminology of [20]. All geometry primitives are defined in the ASTRA Toolbox as well, except for the volume vector geometry that can represent an arbitrarily oriented moving reconstruction grid. The geometries are illustrated in Figure 2 with accompanying code.
In contrast to the standard projection geometries, whose movement is parameterized by the rotation angle, vector geometries move arbitrarily in time. We therefore refer to the state of the acquisition geometry at a specific time as a time step. Furthermore, geometries have a num_steps property that describes in how many time steps their movement is discretized.

Standard geometries
Volume geometry. A volume geometry is created with the ts.volume function and describes the position and size of an axis-aligned voxel grid on which the object is reconstructed. A volume Fig. 2. Creation of typical tomographic geometries. From left to right: a volume geometry, single-axis parallel beam geometry, and a circular cone beam geometry. Below, arbitrarily oriented vector geometries are shown. The parameters are specified using keyword-only arguments [27]. The pos parameter, for instance, determines the position of a volume, other parameters have accompanying labels in the diagrams. The diagrams are created using tomosipo.
In contrast to the standard projection geometries, whose movement is parameterized by the rotation angle, vector geometries move arbitrarily in time. We therefore refer to the state of the acquisition geometry at a specific time as a time step. Furthermore, geometries have a num_steps property that describes in how many time steps their movement is discretized.

Standard geometries
Volume geometry. A volume geometry is created with the ts.volume function and describes the position and size of an axis-aligned voxel grid on which the object is reconstructed. A volume geometry can be created with size, pos, and shape parameters, which define its physical size, center position, and the number of voxels in each direction. By default, the volume is centered on the origin, and if the size is not specified, it is set to equal the shape, causing the voxel size to equal 1. Other parametrizations, such as in terms of the volume's extents, are described in the documentation.
Single-axis parallel beam. In the parallel beam geometry, X-rays run along parallel lines and are collected on a flat panel detector that rotates around a single axis on the origin. It can be created with the ts.parallel function which has parameters size, shape, and angles, which define the detector's physical size, the number of pixels in each dimension, and the rotation angles. If an integer argument is provided for angles, equispaced rotation angles in the interval [0, π) are used. Otherwise, a provided array is interpreted as containing the rotation angles in radians.
Circular cone beam. Like the parallel beam geometry, the flat panel detector of a cone beam geometry rotates around an axis located on the origin. A cone beam geometry is created with the ts.cone function that also takes the angles, shape, and size parameters. In contrast to the parallel beam geometry, the rays in a cone beam geometry are emitted from a point source, and the source-to-origin distance and source-to-detector distances can be specified using the src_orig_dist and src_det_dist parameters. Also, when angles is provided as an integer, a rotation is performed along a full arc [0, 2π) as opposed to [0, π).

Flexible vector geometries
Any geometry g can be converted to a vector geometry by calling g.to_vec(). Vector geometries can also be created directly as described below.
Volume vector geometry. In contrast to a volume geometry, which is static, a volume vector geometry may move over time and the reconstruction grid may be arbitrarily oriented. It can be created with the ts.volume_vec function which takes as parameter the shape of the voxel grid and 3 vectors describing the local frame of reference of the grid at each point in time. In practice, a vector volume geometry is easier to obtain by applying a geometric transformation to a standard volume geometry.
The ASTRA Toolbox, tomosipo's computational back end, does not support non-axis-aligned volume geometries. Internally, tomosipo aligns the volume to the origin and moves the projection geometry with it. The transformed geometries are handed to ASTRA, causing the projection operation to be performed in the frame of reference of the object.
Arbitrarily oriented parallel beam. In a parallel vector geometry, the detector can be arbitrarily oriented and positioned. In addition, the direction of the incoming rays can be adjusted to a direction that is not necessarily orthogonal to the detector plane. It can be created with the ts.parallel_vec function by specifying a fixed detector shape and varying ray directions, detector positions, and detector orientations at each time step. The orientation is determined by parameters det_u and det_v that specify the vector from detector pixel (0, 0) to (0, 1) and (0, 0) to (1, 0), respectively. An example is the dual-axis parallel beam geometry, which is common in electron tomography [28].
Arbitrarily oriented cone beam. In a cone vector geometry, the detector can be arbitrarily oriented and the source can be placed in an arbitrary location. For instance, this geometry can represent a helical cone beam acquisition, as we show in Section 4.1. It can be created with the ts.cone_vec function that has almost the same parameters as ts.parallel_vec: instead of a ray direction, however, a source position must be provided for each time step. In the next section, we describe in more detail how vector geometries can be obtained as transformations of simple geometries.

Geometric transforms
Tomosipo defines geometric transforms that can rotate, translate, scale, and reflect the previously introduced geometries. In addition, the package provides a perspective transform to switch between different frames of reference. The transforms are stand-alone objects instead of functions that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4×4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4×4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector . This way, application of a geometric transform -notably translation -to points and vectors can be performed by matrix multiplication. That is, in homogeneous coordinates, the transformed vector equals M − → v and the transformed point equals M − → p . In code, a vector v and point p in Euclidean coordinates are transformed using a transform T as follows Application of a transform T to a geometry vg is expressed in code as 1

transformed_vg = T * vg
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition T = T 1 • T 2 of two transforms T 1 , T 2 represented by matrices M 1 , M 2 is equal to the matrix product of the matrices, i.e., M = M 1 M 2 . In code, this is expressed as 1

T = T1 * T2
Composition of transforms is demonstrated in Section 4.1, where a helical cone beam geometry is created.
Rigid and scaling transforms. Tomosipo provides the ts.translate, ts.rotate, ts.scale, and ts.reflect functions to create a translation, rotation, scaling, or reflection transform, respectively. These are illustrated in Fig. 3. A transform may change over time, i.e., at each time step it can define a different geometric transformation. The functions that create the transforms are designed to facilitate defining transforms that vary over time.
A translation transform is parameterized by an axis and an array alpha. The displacement vector at time step i is defined by alpha[i] * axis.
A rotation transform is created using the axis angle representation. The axis, pos, and angles parameters describe the orientation and location of the rotation axis, as well as the angle of rotation. Each of these parameters may be provided as an array to define the rotation at multiple time steps. The angles are expressed in radians and the direction of rotation is right-handed by default.
A scaling transform describes a scaling operation centered on a position. The scaling is not necessarily isotropic: some directions can be scaled more than others. An alpha parameter can be used to modulate the scaling at each time step.
A reflection transform describes a reflection in a plane that is parameterized by a position pos and a normal vector axis. Both can be specified as an array, defining a reflection in a moving plane at several time steps.
Perspective. The ts.from_perspective function creates a perspective transform. This function takes a volume and returns the transform that moves the volume back to the origin and rotates it back into a single axis-aligned orientation. All projection geometries have a to_vol() method that describes the frame of reference of the detector at each time step. This makes it easy to create a transform that converts to the detector's frame of reference. In the case of a circular cone beam trajectory, for example, the source and detector rotate around the volume, from the volume's perspective. From the perspective of the detector, on the other hand, the volume rotates.
Application of a transform T to a geometry vg is expressed in code as that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4×4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector . This way, application of a geometric transform -notably translation -to points and vectors can be performed by matrix multiplication. That is, in homogeneous coordinates, the transformed vector equals M − → v and the transformed point equals M − → p . In code, a vector v and point p in Euclidean coordinates are transformed using a transform T as follows Application of a transform T to a geometry vg is expressed in code as 1

transformed_vg = T * vg
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition T = T 1 • T 2 of two transforms T 1 , T 2 represented by matrices M 1 , M 2 is equal to the matrix product of the matrices, i.e., M = M 1 M 2 . In code, this is expressed as 1

T = T1 * T2
Composition of transforms is demonstrated in Section 4.1, where a helical cone beam geometry is created.
Rigid and scaling transforms. Tomosipo provides the ts.translate, ts.rotate, ts.scale, and ts.reflect functions to create a translation, rotation, scaling, or reflection transform, respectively. These are illustrated in Fig. 3. A transform may change over time, i.e., at each time step it can define a different geometric transformation. The functions that create the transforms are designed to facilitate defining transforms that vary over time.
A translation transform is parameterized by an axis and an array alpha. The displacement vector at time step i is defined by alpha[i] * axis.
A rotation transform is created using the axis angle representation. The axis, pos, and angles parameters describe the orientation and location of the rotation axis, as well as the angle of rotation. Each of these parameters may be provided as an array to define the rotation at multiple time steps. The angles are expressed in radians and the direction of rotation is right-handed by default.
A scaling transform describes a scaling operation centered on a position. The scaling is not necessarily isotropic: some directions can be scaled more than others. An alpha parameter can be used to modulate the scaling at each time step.
A reflection transform describes a reflection in a plane that is parameterized by a position pos and a normal vector axis. Both can be specified as an array, defining a reflection in a moving plane at several time steps.
Perspective. The ts.from_perspective function creates a perspective transform. This function takes a volume and returns the transform that moves the volume back to the origin and rotates it back into a single axis-aligned orientation. All projection geometries have a to_vol() method that describes the frame of reference of the detector at each time step. This makes it easy to create a transform that converts to the detector's frame of reference. In the case of a circular cone beam trajectory, for example, the source and detector rotate around the volume, from the volume's perspective. From the perspective of the detector, on the other hand, the volume rotates.
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition T = T 1 • T 2 of two transforms T 1 , T 2 represented by matrices M 1 , M 2 is equal to the matrix product of the matrices, i.e., M = M 1 M 2 . In code, this is expressed as that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4×4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector . This way, application of a geometric transform -notably translation -to points and vectors can be performed by matrix multiplication. That is, in homogeneous coordinates, the transformed vector equals M − → v and the transformed point equals M − → p . In code, a vector v and point p in Euclidean coordinates are transformed using a transform T as follows Application of a transform T to a geometry vg is expressed in code as 1

transformed_vg = T * vg
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition T = T 1 • T 2 of two transforms T 1 , T 2 represented by matrices M 1 , M 2 is equal to the matrix product of the matrices, i.e., M = M 1 M 2 . In code, this is expressed as Composition of transforms is demonstrated in Section 4.1, where a helical cone beam geometry is created.
Rigid and scaling transforms. Tomosipo provides the ts.translate, ts.rotate, ts.scale, and ts.reflect functions to create a translation, rotation, scaling, or reflection transform, respectively. These are illustrated in Fig. 3. A transform may change over time, i.e., at each time step it can define a different geometric transformation. The functions that create the transforms are designed to facilitate defining transforms that vary over time.
A translation transform is parameterized by an axis and an array alpha. The displacement vector at time step i is defined by alpha[i] * axis.
A rotation transform is created using the axis angle representation. The axis, pos, and angles parameters describe the orientation and location of the rotation axis, as well as the angle of rotation. Each of these parameters may be provided as an array to define the rotation at multiple time steps. The angles are expressed in radians and the direction of rotation is right-handed by default.
A scaling transform describes a scaling operation centered on a position. The scaling is not necessarily isotropic: some directions can be scaled more than others. An alpha parameter can be used to modulate the scaling at each time step.
A reflection transform describes a reflection in a plane that is parameterized by a position pos and a normal vector axis. Both can be specified as an array, defining a reflection in a moving plane at several time steps.
Perspective. The ts.from_perspective function creates a perspective transform. This function takes a volume and returns the transform that moves the volume back to the origin and rotates it back into a single axis-aligned orientation. All projection geometries have a to_vol() method that describes the frame of reference of the detector at each time step. This makes it easy to create a transform that converts to the detector's frame of reference. In the case of a circular cone beam trajectory, for example, the source and detector rotate around the volume, from the volume's perspective. From the perspective of the detector, on the other hand, the volume rotates.
Composition of transforms is demonstrated in Section 4.1, where a helical cone beam geometry is created.
Rigid and scaling transforms. Tomosipo provides the ts.translate, ts.rotate, ts.scale, and ts.reflect functions to create a translation, rotation, scaling, or reflection transform, respectively. These are illustrated in Fig. 3. A transform may change over time, i.e., at each time step it can define a different geometric transformation. The functions that create the transforms are designed to facilitate defining transforms that vary over time.
A translation transform is parameterized by an axis and an array alpha. The displacement vector at time step i is defined by alpha[i] * axis.
A rotation transform is created using the axis angle representation. The axis, pos, and angles parameters describe the orientation and location of the rotation axis, as well as the angle of rotation. Each of these parameters may be provided as an array to define the rotation at multiple time steps. The angles are expressed in radians and the direction of rotation is right-handed by default.
A scaling transform describes a scaling operation centered on a position. The scaling is not necessarily isotropic: some directions can be scaled more than others. An alpha parameter can be used to modulate the scaling at each time step.
A reflection transform describes a reflection in a plane that is parameterized by a position pos and a normal vector axis. Both can be specified as an array, defining a reflection in a moving plane at several time steps.
Perspective. The ts.from_perspective function creates a perspective transform. This function takes a volume and returns the transform that moves the volume back to the origin and rotates it back into a single axis-aligned orientation. All projection geometries have a to_vol() method that describes the frame of reference of the detector at each time step. This makes it easy to create a transform that converts to the detector's frame of reference. In the case of a circular cone beam trajectory, for example, the source and detector rotate around the volume, from the volume's perspective. From the perspective of the detector, on the other hand, the volume rotates. This change in perspective is illustrated in the last two panes of Fig. 3. Both perspectives yield the same projection operator A.

Interoperability and GPU-acceleration
In this section, we discuss tomosipo's interoperability with NumPy arrays [30] and GPUaccelerated Python packages. In addition, we discuss the performance benefits of using GPU-accelerated arrays and also some trade-offs in favor of CPU arrays.

Geometric transforms
Tomosipo defines geometric transforms that can rotate, translate, scale, and reflect the previously introduced geometries. In addition, the package provides a perspective transform to switch between different frames of reference. The transforms are stand-alone objects instead of functions that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4 × 4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector . This way, application of a geometric transform -notably translation -to points and vectors can be performed by matrix multiplication. That is, in homogeneous coordinates, the transformed vector equals Mì v and the transformed point equals Mì p. In code, a vector v and point p in Euclidean coordinates are transformed using a transform T as follows Application of a transform T to a geometry vg is expressed in code as 1

transformed_vg = T * vg
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition = 1 • 2 of two transforms 1 , 2 Projection operations are calculated using the ASTRA Toolbox. We have extended the ASTRA Toolbox API to enable direct operation on NumPy arrays. Before any ASTRA operation, the input arrays are automatically linked to the ASTRA runtime, and unlinked afterwards. This represents a substantial ergonomic improvement over the existing API.
Apart from NumPy arrays, array types from other Python packages can also be linked. Out of the box, tomosipo interoperates with PyTorch and CuPy [24,31]. More integrations can be added though an API, which can be used in the future to add interoperability with a variety of array libraries through the currently developing Python array API standard [32].
Integration with GPU array libraries can enable substantial performance improvements. In the snippet below, a NumPy array and a PyTorch array are forward projected using a projection operator A. The NumPy array is located in RAM attached to the CPU, and the PyTorch array is located on the GPU.

Geometric transforms
Tomosipo defines geometric transforms that can rotate, translate, scale, and reflect the previously introduced geometries. In addition, the package provides a perspective transform to switch between different frames of reference. The transforms are stand-alone objects instead of functions that act on geometries directly. We first discuss the internal representation of the transforms and then we introduce the built-in functions to create transforms.
Representation. Internally, the transform at each time step is represented by a 4 × 4 matrix M describing the transformation in homogeneous coordinates [29]. An orientation vector . This way, application of a geometric transform -notably translation -to points and vectors can be performed by matrix multiplication. That is, in homogeneous coordinates, the transformed vector equals Mì v and the transformed point equals Mì p. In code, a vector v and point p in Euclidean coordinates are transformed using a transform T as follows Application of a transform T to a geometry vg is expressed in code as 1

transformed_vg = T * vg
In the internal representation, the composition of two transforms is also computed by matrix multiplication. The matrix representation of the composition = 1 • 2 of two transforms 1 , 2 This change in perspective is illustrated in the last two panes of Fig. 3. Both perspectives yield the same projection operator A.

Interoperability and GPU-acceleration
In this section, we discuss tomosipo's interoperability with NumPy arrays [30] and GPUaccelerated Python packages. In addition, we discuss the performance benefits of using GPU-accelerated arrays and also some trade-offs in favor of CPU arrays.
Projection operations are calculated using the ASTRA Toolbox. We have extended the ASTRA Toolbox API to enable direct operation on NumPy arrays. Before any ASTRA operation, the input arrays are automatically linked to the ASTRA runtime, and unlinked afterwards. This represents a substantial ergonomic improvement over the existing API.
Apart from NumPy arrays, array types from other Python packages can also be linked. Out of the box, tomosipo interoperates with PyTorch and CuPy [24,31]. More integrations can be added though an API, which can be used in the future to add interoperability with a variety of array libraries through the currently developing Python array API standard [32].
Integration with GPU array libraries can enable substantial performance improvements. In the snippet below, a NumPy array and a PyTorch array are forward projected using a projection operator A. The NumPy array is located in RAM attached to the CPU, and the PyTorch array is located on the GPU. On line 1, the data is first moved to the GPU, the forward projection is calculated, and the data is moved back to the CPU. On line 2, no data movement takes place: the forward projection is calculated on the GPU. In iterative algorithms, where the forward and backprojection are repeatedly executed substeps of the algorithm, the latency imposed by CPU-GPU communication can dominate the computation time, as we demonstrate in Section 6. Note that PyTorch arrays can also be created on the CPU. In that case, the computation of the forward projection goes through exactly the same steps as a NumPy array would.
There are cases where it is beneficial to keep data on CPU. When data is too big to fit in GPU memory, the ASTRA Toolbox automatically splits data residing on CPU and performs data is moved back to the CPU. On line 2, no data movement takes place: the forward projection is calculated on the GPU. In iterative algorithms, where the forward and backprojection are repeatedly executed substeps of the algorithm, the latency imposed by CPU-GPU communication can dominate the computation time, as we demonstrate in Section 6. Note that PyTorch arrays can also be created on the CPU. In that case, the computation of the forward projection goes through exactly the same steps as a NumPy array would.
There are cases where it is beneficial to keep data on CPU. When data is too big to fit in GPU memory, the ASTRA Toolbox automatically splits data residing on CPU and performs the computation on the GPU in a streaming fashion. In this case, the user does not have to split up the data manually. When multiple GPUs are present on the system, they can be used automatically. In the code below, the ASTRA Toolbox is instructed to use four GPUs on line 1. The computation of the forward projection on line 2 is distributed over the four GPUs.

Splitting and joining geometries
The ability to split and join geometries in tomosipo's API allows users to customize their design easily. Tomosipo allows subsampling a geometry to obtain a sub-geometry. In addition, it allows joining the time steps of sequences of geometries into a single geometry.
First, we demonstrate subsampling of geometries. Here, we limit ourselves to projection geometries. Subsampling of volume geometries and geometric transforms works similarly and is described in the documentation. A projection geometry pg can be subsampled to obtain a geometry describing a subset of the detector surface. In the code below, the detector surface is cropped, removing a border of 100 pixels from each side. On the next line, the detector surface is subsampled, selecting every other row and column of pixels. Subsampling induces a slight shift in the detector's center, which is taken into account and described in detail in the documentation. Subsampling the angular dimension is also possible. In this dimension, subsampling supports both slicing and Boolean masks [30]. In the code below, the angular direction is subsampled, obtaining a geometry that contains every other projection angle. In the line below, angles are selected when a condition array equals True. In addition to subsampling, tomosipo includes the ts.concatenate function that concatenates geometries and transforms. The concatenation of multiple projection geometries combines their time steps into a single geometry. In the code below, two projection geometries pg1 and pg2 are combined. In the next line, a rotation R is repeatedly composed with different translations T1, T2, T3. The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.

Splitting and joining geometries
The ability to split and join geometries in tomosipo's API allows users to customize their design easily. Tomosipo allows subsampling a geometry to obtain a sub-geometry. In addition, it allows joining the time steps of sequences of geometries into a single geometry.
First, we demonstrate subsampling of geometries. Here, we limit ourselves to projection geometries. Subsampling of volume geometries and geometric transforms works similarly and is described in the documentation. A projection geometry pg can be subsampled to obtain a geometry describing a subset of the detector surface. In the code below, the detector surface is cropped, removing a border of 100 pixels from each side. On the next line, the detector surface is subsampled, selecting every other row and column of pixels. Subsampling induces a slight shift in the detector's center, which is taken into account and described in detail in the documentation.
is calculated on the GPU. In iterative algorithms, where the forward and backprojection are repeatedly executed substeps of the algorithm, the latency imposed by CPU-GPU communication can dominate the computation time, as we demonstrate in Section 6. Note that PyTorch arrays can also be created on the CPU. In that case, the computation of the forward projection goes through exactly the same steps as a NumPy array would.
There are cases where it is beneficial to keep data on CPU. When data is too big to fit in GPU memory, the ASTRA Toolbox automatically splits data residing on CPU and performs the computation on the GPU in a streaming fashion. In this case, the user does not have to split up the data manually. When multiple GPUs are present on the system, they can be used automatically. In the code below, the ASTRA Toolbox is instructed to use four GPUs on line 1. The computation of the forward projection on line 2 is distributed over the four GPUs.

Splitting and joining geometries
The ability to split and join geometries in tomosipo's API allows users to customize their design easily. Tomosipo allows subsampling a geometry to obtain a sub-geometry. In addition, it allows joining the time steps of sequences of geometries into a single geometry.
First, we demonstrate subsampling of geometries. Here, we limit ourselves to projection geometries. Subsampling of volume geometries and geometric transforms works similarly and is described in the documentation. A projection geometry pg can be subsampled to obtain a geometry describing a subset of the detector surface. In the code below, the detector surface is cropped, removing a border of 100 pixels from each side. On the next line, the detector surface is subsampled, selecting every other row and column of pixels. Subsampling induces a slight shift in the detector's center, which is taken into account and described in detail in the documentation. Subsampling the angular dimension is also possible. In this dimension, subsampling supports both slicing and Boolean masks [30]. In the code below, the angular direction is subsampled, obtaining a geometry that contains every other projection angle. In the line below, angles are selected when a condition array equals True. In addition to subsampling, tomosipo includes the ts.concatenate function that concatenates geometries and transforms. The concatenation of multiple projection geometries combines their time steps into a single geometry. In the code below, two projection geometries pg1 and pg2 are combined. In the next line, a rotation R is repeatedly composed with different translations T1, T2, T3. The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.
Subsampling the angular dimension is also possible. In this dimension, subsampling supports both slicing and Boolean masks [30]. In the code below, the angular direction is subsampled, obtaining a geometry that contains every other projection angle. In the line below, angles are selected when a condition array equals True. is calculated on the GPU. In iterative algorithms, where the forward and backprojection are repeatedly executed substeps of the algorithm, the latency imposed by CPU-GPU communication can dominate the computation time, as we demonstrate in Section 6. Note that PyTorch arrays can also be created on the CPU. In that case, the computation of the forward projection goes through exactly the same steps as a NumPy array would.
There are cases where it is beneficial to keep data on CPU. When data is too big to fit in GPU memory, the ASTRA Toolbox automatically splits data residing on CPU and performs the computation on the GPU in a streaming fashion. In this case, the user does not have to split up the data manually. When multiple GPUs are present on the system, they can be used automatically. In the code below, the ASTRA Toolbox is instructed to use four GPUs on line 1. The computation of the forward projection on line 2 is distributed over the four GPUs.

Splitting and joining geometries
The ability to split and join geometries in tomosipo's API allows users to customize their design easily. Tomosipo allows subsampling a geometry to obtain a sub-geometry. In addition, it allows joining the time steps of sequences of geometries into a single geometry.
First, we demonstrate subsampling of geometries. Here, we limit ourselves to projection geometries. Subsampling of volume geometries and geometric transforms works similarly and is described in the documentation. A projection geometry pg can be subsampled to obtain a geometry describing a subset of the detector surface. In the code below, the detector surface is cropped, removing a border of 100 pixels from each side. On the next line, the detector surface is subsampled, selecting every other row and column of pixels. Subsampling induces a slight shift in the detector's center, which is taken into account and described in detail in the documentation. Subsampling the angular dimension is also possible. In this dimension, subsampling supports both slicing and Boolean masks [30]. In the code below, the angular direction is subsampled, obtaining a geometry that contains every other projection angle. In the line below, angles are selected when a condition array equals True. In addition to subsampling, tomosipo includes the ts.concatenate function that concatenates geometries and transforms. The concatenation of multiple projection geometries combines their time steps into a single geometry. In the code below, two projection geometries pg1 and pg2 are combined. In the next line, a rotation R is repeatedly composed with different translations T1, T2, T3. The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.
In Section 4.3, Boolean masking is demonstrated in the case study of X-ray diffraction tomography, where diffraction only occurs in a subset of projection angles.
In addition to subsampling, tomosipo includes the ts.concatenate function that concatenates geometries and transforms. The concatenation of multiple projection geometries combines their time steps into a single geometry. In the code below, two projection geometries pg1 and pg2 are combined. In the next line, a rotation R is repeatedly composed with different translations T1, T2, T3.
is calculated on the GPU. In iterative algorithms, where the forward and backprojection are repeatedly executed substeps of the algorithm, the latency imposed by CPU-GPU communication can dominate the computation time, as we demonstrate in Section 6. Note that PyTorch arrays can also be created on the CPU. In that case, the computation of the forward projection goes through exactly the same steps as a NumPy array would.
There are cases where it is beneficial to keep data on CPU. When data is too big to fit in GPU memory, the ASTRA Toolbox automatically splits data residing on CPU and performs the computation on the GPU in a streaming fashion. In this case, the user does not have to split up the data manually. When multiple GPUs are present on the system, they can be used automatically. In the code below, the ASTRA Toolbox is instructed to use four GPUs on line 1. The computation of the forward projection on line 2 is distributed over the four GPUs.

Splitting and joining geometries
The ability to split and join geometries in tomosipo's API allows users to customize their design easily. Tomosipo allows subsampling a geometry to obtain a sub-geometry. In addition, it allows joining the time steps of sequences of geometries into a single geometry.
First, we demonstrate subsampling of geometries. Here, we limit ourselves to projection geometries. Subsampling of volume geometries and geometric transforms works similarly and is described in the documentation. A projection geometry pg can be subsampled to obtain a geometry describing a subset of the detector surface. In the code below, the detector surface is cropped, removing a border of 100 pixels from each side. On the next line, the detector surface is subsampled, selecting every other row and column of pixels. Subsampling induces a slight shift in the detector's center, which is taken into account and described in detail in the documentation. Subsampling the angular dimension is also possible. In this dimension, subsampling supports both slicing and Boolean masks [30]. In the code below, the angular direction is subsampled, obtaining a geometry that contains every other projection angle. In the line below, angles are selected when a condition array equals True. In addition to subsampling, tomosipo includes the ts.concatenate function that concatenates geometries and transforms. The concatenation of multiple projection geometries combines their time steps into a single geometry. In the code below, two projection geometries pg1 and pg2 are combined. In the next line, a rotation R is repeatedly composed with different translations T1, T2, T3. The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.
The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.

Visualization
Tomosipo provides extensive support for visualizing geometries. Animations can be saved as a video or as a scalable vector graphic (SVG). In addition, geometries can be investigated in a 3D-accelerated environment on the desktop, allowing the user to zoom, pan, and rotate the view. Finally, an interactive SVG animation can be shown in an online Jupyter notebook, allowing for quick inspection of intermediate results. These options are illustrated in Fig. 4. All other illustrations in this article have been generated using tomosipo. They were saved in the SVG format and extended using Inkscape.

Research Article
Vol. 29 The concatenation of transforms is demonstrated in the case study of X-ray scattering tensor tomography in Section 4.4, where it is used to define a repeated rotation at several tilt angles.

Visualization
Tomosipo provides extensive support for visualizing geometries. Animations can be saved as a video or as a scalable vector graphic (SVG). In addition, geometries can be investigated in a 3D-accelerated environment on the desktop, allowing the user to zoom, pan, and rotate the view. Finally, an interactive SVG animation can be shown in an online Jupyter notebook, allowing for quick inspection of intermediate results. These options are illustrated in Figure 4. All other illustrations in this article have been generated using tomosipo. They were saved in the SVG format and extended using Inkscape.

Case studies
In this section, the concepts developed in the previous section are put into practice. We describe two simple examples and two complex acquisition schemes that are in use at synchrotron tomography beamlines. The first example demonstrates how geometric transforms can be composed to create a helical cone beam geometry. The second example models single-axis parallel beam tomography with a non-standard center of rotation in the frame of reference of the laboratory. In the first case study, we describe X-ray diffraction contrast tomography, which demonstrates the use of reflection and subsampling using a Boolean mask. In the second case study, we describe X-ray scattering tensor tomography, which demonstrates the use of concatenation. The case studies demonstrate that X-ray diffraction and scattering can be modeled using tomosipo's projection operators.

Basic example: Helical cone beam geometry
As a demonstration of the composition of two primitive transforms, we define the helical cone beam geometry [34]. Here, the source and detector follow a helical path around the object. Using the notation of [34], we describe the helical geometry as a composition of translation and rotation in Figure 5. First, a static volume and a static cone beam geometry are defined. Next, a rotation and translation are defined, which rotate around and translate along the z-axis. The helical transform is defined such that it applies the rotation and then a translation at time step . When it is applied to the cone beam geometry, the resulting trajectory of the source and detector is helical. We note that the helical trajectory could have been obtained as a translation of a non-static cone beam geometry, effectively hiding the rotation in the cone beam geometry.

Case studies
In this section, the concepts developed in the previous section are put into practice. We describe two simple examples and two complex acquisition schemes that are in use at synchrotron tomography beamlines. The first example demonstrates how geometric transforms can be composed to create a helical cone beam geometry. The second example models single-axis parallel beam tomography with a non-standard center of rotation in the frame of reference of the laboratory. In the first case study, we describe X-ray diffraction contrast tomography, which demonstrates the use of reflection and subsampling using a Boolean mask. In the second case study, we describe X-ray scattering tensor tomography, which demonstrates the use of concatenation. The case studies demonstrate that X-ray diffraction and scattering can be modeled using tomosipo's projection operators.

Basic example: Helical cone beam geometry
As a demonstration of the composition of two primitive transforms, we define the helical cone beam geometry [34]. Here, the source and detector follow a helical path around the object. Using the notation of [34], we describe the helical geometry as a composition of translation and rotation in Fig. 5. First, a static volume and a static cone beam geometry are defined. Next, a rotation R and translation T are defined, which rotate around and translate along the z-axis. The helical transform H is defined such that it applies the rotation R i and then a translation T i at time step i. When it is applied to the cone beam geometry, the resulting trajectory of the source and detector is helical. We note that the helical trajectory could have been obtained as a translation of a non-static cone beam geometry, effectively hiding the rotation in the cone beam geometry.

Basic example: Parallel beam in the lab frame
Acquisition using the single-axis parallel beam geometry is common at synchrotron beam lines. The detector is often positioned at a fixed location and the sample is mounted on a movable rotation stage. Typically, it is assumed that the center of rotation and the center of the detector coincide. In many cases in practice, however, it is difficult to achieve this with the described setup. Therefore, the offset between the center of rotation and the center of the detector has to be taken into account in order to achieve an accurate reconstruction. This is possible in tomosipo by positioning the detector, volume, and rotation axis independently from each other.
In Fig. 6, the acquisition geometry is defined in the frame of reference of the laboratory. First, a static detector is translated from the origin to its final position by a transform T. Next, a static volume geometry is created at the initial position of the sample. A rotation is defined with a specific position of the rotation axis. Finally, the rotation is applied to the static volume, obtaining a rotating volume whose center rotates around the rotation axis.

Basic example: Parallel beam in the lab frame
Acquisition using the single-axis parallel beam geometry is common at synchrotron beam lines. The detector is often positioned at a fixed location and the sample is mounted on a movable rotation stage. Typically, it is assumed that the center of rotation and the center of the detector coincide. In many cases in practice, however, it is difficult to achieve this with the described setup. Therefore, the offset between the center of rotation and the center of the detector has to be taken into account in order to achieve an accurate reconstruction. This is possible in tomosipo by positioning the detector, volume, and rotation axis independently from each other. In Figure 6, the acquisition geometry is defined in the frame of reference of the laboratory. First, a static detector is translated from the origin to its final position by a transform T. Next, a static volume geometry is created at the initial position of the sample. A rotation is defined with a specific position of the rotation axis. Finally, the rotation is applied to the static volume, obtaining a rotating volume whose center rotates around the rotation axis.
There is an advantage to this formulation. In existing tomography packages, the position of the volume is commonly chosen to coincide with the rotation axis. However, this causes the reconstructed images to be translated when a different center of rotation is provided. This can cause problems when the determination of the correct center of rotation is based on the reconstructed images. In contrast, the proposed formulation opens up the possibility of determining the correct center of rotation by maximizing the auto-correlation in the reconstruction at several values of the center of rotation.  6. A single-axis parallel beam acquisition with a custom center of rotation. An object, whose changing position during rotation is indicated in blue, is rotated around a non-standard axis of rotation (in red). The center of the detector is indicated in green.

Basic example: Parallel beam in the lab frame
Acquisition using the single-axis parallel beam geometry is common at synchrotron beam lines. The detector is often positioned at a fixed location and the sample is mounted on a movable rotation stage. Typically, it is assumed that the center of rotation and the center of the detector coincide. In many cases in practice, however, it is difficult to achieve this with the described setup. Therefore, the offset between the center of rotation and the center of the detector has to be taken into account in order to achieve an accurate reconstruction. This is possible in tomosipo by positioning the detector, volume, and rotation axis independently from each other.
In Figure 6, the acquisition geometry is defined in the frame of reference of the laboratory. First, a static detector is translated from the origin to its final position by a transform T. Next, a static volume geometry is created at the initial position of the sample. A rotation is defined with a specific position of the rotation axis. Finally, the rotation is applied to the static volume, obtaining a rotating volume whose center rotates around the rotation axis.
There is an advantage to this formulation. In existing tomography packages, the position of the volume is commonly chosen to coincide with the rotation axis. However, this causes the reconstructed images to be translated when a different center of rotation is provided. This can cause problems when the determination of the correct center of rotation is based on the reconstructed images. In contrast, the proposed formulation opens up the possibility of determining the correct center of rotation by maximizing the auto-correlation in the reconstruction at several values of the center of rotation.  6. A single-axis parallel beam acquisition with a custom center of rotation. An object, whose changing position during rotation is indicated in blue, is rotated around a non-standard axis of rotation (in red). The center of the detector is indicated in green.

Fig. 6.
A single-axis parallel beam acquisition with a custom center of rotation. An object, whose changing position during rotation is indicated in blue, is rotated around a non-standard axis of rotation (in red). The center of the detector is indicated in green.
There is an advantage to this formulation. In existing tomography packages, the position of the volume is commonly chosen to coincide with the rotation axis. However, this causes the reconstructed images to be translated when a different center of rotation is provided. This can cause problems when the determination of the correct center of rotation is based on the reconstructed images. In contrast, the proposed formulation opens up the possibility of determining the correct center of rotation by maximizing the auto-correlation in the reconstruction at several values of the center of rotation.

Complex case study: X-ray diffraction contrast tomography
X-ray diffraction contrast tomography (DCT) [2] is an imaging technique used to investigate the internal structure of poly-crystalline materials. The crystal lattice is divided into grains, homogeneous regions where the lattice has a similar orientation. The orientation, size, shape, and arrangement of individual grains strongly influence macroscopic properties of the material. Therefore, mapping the orientation of grains is important to characterize a poly-crystalline material [35]. Here, we take as an example the three-dimensional DCT acquisition geometry as described in [2] to demonstrate specific features of tomosipo.
The goal of DCT is to reconstruct a vector field representing the intra-granular orientation of the crystal lattice. This is achieved by discretizing the orientation space on a regular grid that can be represented by unit vectorsô 1 , . . . ,ô N o . For each orientationô k , a scalar field, i.e., a volume, is reconstructed that represents the diffraction "intensity" at that orientation. A variational reconstruction algorithm ensures that neighboring voxels have similar orientations.
A crystal lattice reflects an incoming X-ray beam at specific incidence directions, characterized by the so-called Bragg angles. When the diffraction geometry of the material under investigation is known beforehand, the intra-granular orientation of the crystal lattice can be recovered from those projection images at which Bragg diffraction is expected to occur.
As shown in Fig. 7, the acquisition uses a monochromatic parallel box beam and the diffracted signal of the sample is measured on a flat-panel detector. As the sample is rotated, the reflection of the incoming beam in a voxel with local orientationô forms a figure of eight on the detector. Bragg diffraction only occurs in the instances where the beam and local lattice are in Bragg condition. Occurrence of the Bragg condition is relatively rare, as shown in Fig. 7(b). Both the reflection and its intermittent nature can be modeled in tomosipo.

Complex case study: X-ray diffraction contrast tomography
X-ray diffraction contrast tomography (DCT) [2] is an imaging technique used to investigate the internal structure of poly-crystalline materials. The crystal lattice is divided into grains, homogeneous regions where the lattice has a similar orientation. The orientation, size, shape, and arrangement of individual grains strongly influence macroscopic properties of the material. Therefore, mapping the orientation of grains is important to characterize a poly-crystalline material [35]. Here, we take as an example the three-dimensional DCT acquisition geometry as described in [2] to demonstrate specific features of tomosipo. The goal of DCT is to reconstruct a vector field representing the intra-granular orientation of the crystal lattice. This is achieved by discretizing the orientation space on a regular grid that can be represented by unit vectorsô 1 , . . . ,ô . For each orientationô , a scalar field, i.e., a volume, is reconstructed that represents the diffraction "intensity" at that orientation. A variational reconstruction algorithm ensures that neighboring voxels have similar orientations.
A crystal lattice reflects an incoming X-ray beam at specific incidence directions, characterized by the so-called Bragg angles. When the diffraction geometry of the material under investigation is known beforehand, the intra-granular orientation of the crystal lattice can be recovered from those projection images at which Bragg diffraction is expected to occur.
As shown in Figure 7, the acquisition uses a monochromatic parallel box beam and the diffracted signal of the sample is measured on a flat-panel detector. As the sample is rotated, the reflection of the incoming beam in a voxel with local orientationô forms a figure of eight on the detector. Bragg diffraction only occurs in the instances where the beam and local lattice are in Bragg condition. Occurrence of the Bragg condition is relatively rare, as shown in Figure 7 (b). Both the reflection and its intermittent nature can be modeled in tomosipo.
Bragg diffraction (reflection). The diffraction, i.e., reflection, of an incident X-ray beam can be represented in tomosipo. As shown in Figure 7, a parallel bundle of rays remains parallel after it has been reflected. Therefore, the measurement of a diffracted parallel beam can be modeled using a standard parallel-beam geometry with altered ray direction.  7. (a) In X-ray diffraction contrast tomography, a crystal sample is illuminated by a monochromatic X-ray box beam. The crystal sample is divided into grains, which have a minimal spread in local orientation. As the sample is rotated, the incoming beam is diffracted when its incident angle with the local lattice plane equals the Bragg angle. The intersection points of the reflected beam with the detector form a figure of eight, on which Bragg diffraction occurs only twice (marked in red) as the sample is rotated. (b) For a random sample of 20 orientations, the occurrence of Bragg diffraction at a rotation angle is indicated in black. Here, diffraction occurs in just 3.3% of the orientation-rotation combinations.
Bragg diffraction (reflection). The diffraction, i.e., reflection, of an incident X-ray beam can be represented in tomosipo. As shown in Fig. 7, a parallel bundle of rays remains parallel after it has been reflected. Therefore, the measurement of a diffracted parallel beam can be modeled using a standard parallel-beam geometry with altered ray direction.
The code below models the reflection of the incoming beam by a rotating crystal lattice. The orientation of the lattice is represented by a plane normal vector. First, the plane normal of the crystal lattice is rotated. Then, a reflection M is created in the rotating plane normal. The position of the reflection is arbitrary, as it is used to transform the direction of the beam and not its location. Finally, a static parallel beam geometry is modified such that the ray direction corresponds to that of the beam reflected in the rotating plane normal. The position and orientation of the detector remain static. Bragg condition (Boolean masking ). The occurrence of Bragg diffraction can be considered as a Boolean mask, an example of which is shown in Fig. 7(b). It is computed in the code below. First, the plane normal is rotated. Then, the Bragg condition is determined at each rotation angle. The created Boolean mask is used to select a subset of each projection geometry. For each orientation, the code below creates an operator that computes the forward projection only at rotation angles where Bragg diffraction occurs. The created Boolean mask is used to select a subset of each projection geometry. For each orientation, the code below creates an operator that computes the forward projection only at rotation angles where Bragg diffraction occurs. Multi-orientation tomography (sums of masked operators). The forward projection computes the diffraction pattern of x ∈ R N o ×N 3 v , representing all discretized plane orientations at N 3 v locations, onto y ∈ R N θ ×N 2 p , representing the N 2 p pixel intensities at N θ rotation angles. The operation is a linear combination of the masked operators defined above. The backprojection operation can be defined similarly as we show in Code 1 (Ref. [36]), containing the full code listing.
The created Boolean mask is used to select a subset of each projection geometry. For each orientation, the code below creates an operator that computes the forward projection only at rotation angles where Bragg diffraction occurs. Implementing the variational reconstruction algorithm described in [2] is outside of the scope of this manuscript. We have shown how the DCT geometry can be succinctly expressed using tomosipo's rotation and reflection transformations. In addition, we have used subsampling with a Boolean mask to limit the forward projection to the few instances where Bragg diffraction occurs.

4.4.
Complex case study: X-ray scattering tensor tomography X-ray scattering tensor tomography (XSTT) is an imaging technique used to investigate materials with micro-and nano-scale structures over an orders of magnitude larger volumetric field of view, compared to conventional tomographic modalities [36,37]. Here, we take the XSTT acquisition geometry that is described in [3] as an example to demonstrate specific features of tomosipo.
The goal of XSTT is to reconstruct a vector field representing the directional scattering intensity of a sample. This is achieved by reconstructing ŝ ≥ 6 scalar fields that represent the squared scattering coefficient along unit vectorŝ 1 , . . . ,ŝ ŝ at each voxel. After reconstruction, the directional scattering intensities are fine-tuned using per-voxel PCA (principal component analysis) [38]. XSTT has various biological and industrial applications [3]. As an example, the recovered local directional scattering intensities can be used to predict macroscopic properties of fibrous materials. These properties depend on the local fiber arrangement. Fibers scatter X-rays the least along their primary fiber orientation. Therefore, the local fiber orientation can be recovered from the shortest principal axis (smallest scattering magnitude) of the fitted scattering ellipsoid. The possibility to investigate these local structures over large enough volumes is valuable for the research and development of new materials.
We describe the acquisition process to obtain one of the scalar fields xŝ, representing the squared scattering coefficient along a unit vectorŝ. First, we discuss the forward model at a single Because the calculation is performed in the lab frame, the vectorŝ is rotated rather than the beam directionb or sensitivity vectorĝ, Scaled linear combinations. After ν ∈ R Nŝ×Nĝ×N tilt N rot is calculated for all values ofŝ 1 , . . . ,ŝ Nŝ , allĝ i , and all tilts and rotations, then the full projection can be calculated. Here, the measurement alongĝ i is the sum of the contributions of the Nŝ scalar fields representing the scattering coefficients of the sample, as calculated below. The backprojection is defined in Code 1 (Ref. [36]). Implementing the variational reconstruction algorithm described in [2] is outside of the scope of this manuscript. We have shown how the DCT geometry can be succinctly expressed using tomosipo's rotation and reflection transformations. In addition, we have used subsampling with a Boolean mask to limit the forward projection to the few instances where Bragg diffraction occurs.
4.4. Complex case study: X-ray scattering tensor tomography X-ray scattering tensor tomography (XSTT) is an imaging technique used to investigate materials with micro-and nano-scale structures over an orders of magnitude larger volumetric field of view, compared to conventional tomographic modalities [36,37]. Here, we take the XSTT acquisition geometry that is described in [3] as an example to demonstrate specific features of tomosipo.
The goal of XSTT is to reconstruct a vector field representing the directional scattering intensity of a sample. This is achieved by reconstructing ŝ ≥ 6 scalar fields that represent the squared scattering coefficient along unit vectorŝ 1 , . . . ,ŝ ŝ at each voxel. After reconstruction, the directional scattering intensities are fine-tuned using per-voxel PCA (principal component analysis) [38]. XSTT has various biological and industrial applications [3]. As an example, the recovered local directional scattering intensities can be used to predict macroscopic properties of fibrous materials. These properties depend on the local fiber arrangement. Fibers scatter Implementing the variational reconstruction algorithm described in [2] is outside of the scope of this manuscript. We have shown how the DCT geometry can be succinctly expressed using tomosipo's rotation and reflection transformations. In addition, we have used subsampling with a Boolean mask to limit the forward projection to the few instances where Bragg diffraction occurs.

Complex case study: X-ray scattering tensor tomography
X-ray scattering tensor tomography (XSTT) is an imaging technique used to investigate materials with micro-and nano-scale structures over an orders of magnitude larger volumetric field of view, compared to conventional tomographic modalities [37,38]. Here, we take the XSTT acquisition geometry that is described in [3] as an example to demonstrate specific features of tomosipo.
The goal of XSTT is to reconstruct a vector field representing the directional scattering intensity of a sample. This is achieved by reconstructing Nŝ ≥ 6 scalar fields that represent the squared scattering coefficient along unit vectorŝ 1 , . . . ,ŝ Nŝ at each voxel. After reconstruction, the directional scattering intensities are fine-tuned using per-voxel PCA (principal component analysis) [39]. XSTT has various biological and industrial applications [3]. As an example, the recovered local directional scattering intensities can be used to predict macroscopic properties of fibrous materials. These properties depend on the local fiber arrangement. Fibers scatter X-rays the least along their primary fiber orientation. Therefore, the local fiber orientation can be recovered from the shortest principal axis (smallest scattering magnitude) of the fitted scattering ellipsoid. The possibility to investigate these local structures over large enough volumes is valuable for the research and development of new materials.
We describe the acquisition process to obtain one of the scalar fields xŝ, representing the squared scattering coefficient along a unit vectorŝ. First, we discuss the forward model at a single orientation of the sample, i.e., without any rotation or tilting. Let xŝ ∈ R N 3 v represent the sample's squared scattering coefficient along a vectorŝ. The sample is illuminated by a monochromatic parallel X-ray beam. Before they are measured on a detector, the X-rays travel through a panel that is etched with a periodic array of multi-circular gratings [3], generating a reference pattern. The panel is placed at a fixed propagation distance from the detector to maximize the visibility of the patterns. Different types of gratings require different acquisition geometries. The acquisition discussed in this case study is specifically geared to circular gratings.
With the use of circular gratings, the pixels of the detector are grouped into 9×9 pixel unit cells. In each unit cell, a 2D directional intensity is measured along multiple unit vectorsĝ i . The measured intensity y i along the vectorĝ i on the detector for a single beam directionb is computed by scaling the forward projection with the scalar νb ,ŝ,ĝ i [37], defined by The scaling is the same for each unit cell on the detector and varies as the sample (and thusŝ) is rotated.
Rotation and tilt (concatenation). When using circular gratings, the sample must be measured with multiple tilted rotation axes [3,40]. In the XSTT acquisition described in [3], the sample stage is rotated, while the stage is tilted in steps. At each step, the stage makes a full rotation, as illustrated in Fig. 8. Each step can be represented in tomosipo by composing a single tilt operation with a full rotation. In the code below, the full motion is computed by concatenating each step. Implementing the variational reconstruction algorithm described in [2] is outside of the scope of this manuscript. We have shown how the DCT geometry can be succinctly expressed using tomosipo's rotation and reflection transformations. In addition, we have used subsampling with a Boolean mask to limit the forward projection to the few instances where Bragg diffraction occurs.
4.4. Complex case study: X-ray scattering tensor tomography X-ray scattering tensor tomography (XSTT) is an imaging technique used to investigate materials with micro-and nano-scale structures over an orders of magnitude larger volumetric field of view, compared to conventional tomographic modalities [36,37]. Here, we take the XSTT acquisition geometry that is described in [3] as an example to demonstrate specific features of tomosipo.
The goal of XSTT is to reconstruct a vector field representing the directional scattering intensity of a sample. This is achieved by reconstructing ŝ ≥ 6 scalar fields that represent the squared scattering coefficient along unit vectorŝ 1 , . . . ,ŝ ŝ at each voxel. After reconstruction, the directional scattering intensities are fine-tuned using per-voxel PCA (principal component analysis) [38]. XSTT has various biological and industrial applications [3]. As an example, the recovered local directional scattering intensities can be used to predict macroscopic properties of fibrous materials. These properties depend on the local fiber arrangement. Fibers scatter X-rays the least along their primary fiber orientation. Therefore, the local fiber orientation can be recovered from the shortest principal axis (smallest scattering magnitude) of the fitted scattering ellipsoid. The possibility to investigate these local structures over large enough volumes is valuable for the research and development of new materials.
We describe the acquisition process to obtain one of the scalar fields xŝ, representing the squared scattering coefficient along a unit vectorŝ. First, we discuss the forward model at a single G r a t in g D e te c to r 360° rotation U n it c e ll 0 -45° tilt Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9 × 9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorsĝ . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ are reconstructed. Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9×9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorŝ g i . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ k are reconstructed.

Voxel Beam
allĝ i , and all tilts and rotations, then the full projection can be calculated. Here, the measurement alongĝ i is the sum of the contributions of the Nŝ scalar fields representing the scattering coefficients of the sample, as calculated below. The backprojection is defined in Code 1 (Ref. [36] return y Data size. The reconstruction problem considered in [3] fits in memory on modern GPUs. The measured data consists of 46 tilt angles, 50 rotation angles, and 100×144 unit cells. Measuring along 8ĝ i vectors, the total number of measured unit cells equals 46 × 50 × 100 × 144 × 8 ≈ 25 × 10 6 , which requires approximately 1 GB when stored in 32 bit precision. The reconstruction volume consists of 44×71×71 voxels, repeated for each of Nŝ = 7 scattering directions. In total, it requires roughly 6 MB to store in 32 bit precision. The size of the scaling matrix ν is negligible in comparison. Modern data center GPUs range in memory size from 16 GB to 80 GB. Therefore, it is possible to run an iterative SIRT reconstruction of the full problem on GPU. Benchmarks comparing the performance on GPU versus CPU are provided in Section 6.

occurs.
4.4. Complex case study: X-ray scattering tensor tomography X-ray scattering tensor tomography (XSTT) is an imaging technique used to investigate materials with micro-and nano-scale structures over an orders of magnitude larger volumetric field of view, compared to conventional tomographic modalities [36,37]. Here, we take the XSTT acquisition geometry that is described in [3] as an example to demonstrate specific features of tomosipo.
The goal of XSTT is to reconstruct a vector field representing the directional scattering intensity of a sample. This is achieved by reconstructing ŝ ≥ 6 scalar fields that represent the squared scattering coefficient along unit vectorŝ 1 , . . . ,ŝ ŝ at each voxel. After reconstruction, the directional scattering intensities are fine-tuned using per-voxel PCA (principal component analysis) [38]. XSTT has various biological and industrial applications [3]. As an example, the recovered local directional scattering intensities can be used to predict macroscopic properties of fibrous materials. These properties depend on the local fiber arrangement. Fibers scatter X-rays the least along their primary fiber orientation. Therefore, the local fiber orientation can be recovered from the shortest principal axis (smallest scattering magnitude) of the fitted scattering ellipsoid. The possibility to investigate these local structures over large enough volumes is valuable for the research and development of new materials.
We describe the acquisition process to obtain one of the scalar fields xŝ, representing the squared scattering coefficient along a unit vectorŝ. First, we discuss the forward model at a single G r a t in g D e te c to r Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9 × 9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorsĝ . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ are reconstructed. Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9×9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorŝ g i . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ k are reconstructed. Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9 × 9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorsĝ . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ are reconstructed. Fig. 8. In X-ray scattering tensor tomography, a sample is illuminated by a box beam. The sample is repeatedly rotated at several tilt angles. The scattered signal passes through an array of gratings before being measured on a detector. The detector pixels are grouped into 9×9 pixel unit cells. In each unit cell, a directional intensity is measured along vectorŝ g i . In each voxel, scattering coefficients for multiple scattering sensitivity vectorsŝ k are reconstructed.

Voxel Beam
At each tilt and rotation angle, the scaling νb ,ŝ,ĝ i from Eq. (3) is calculated as follows: Because the calculation is performed in the lab frame, the vectorŝ is rotated rather than the beam directionb or sensitivity vectorĝ, Scaled linear combinations. After ν ∈ R Nŝ×Nĝ×N tilt N rot is calculated for all values ofŝ 1 , . . . ,ŝ Nŝ , allĝ i , and all tilts and rotations, then the full projection can be calculated. Here, the measurement alongĝ i is the sum of the contributions of the Nŝ scalar fields representing the scattering coefficients of the sample, as calculated below. The backprojection is defined in Code 1 (Ref. [36]).
Because the calculation is performed in the lab frame, the vectorŝ is rotated rather than the beam directionb or sensitivity vectorĝ, Scaled linear combinations. After ν ∈ R Nŝ×Nĝ×N tilt N rot is calculated for all values ofŝ 1 , . . . ,ŝ Nŝ , allĝ i , and all tilts and rotations, then the full projection can be calculated. Here, the measurement alongĝ i is the sum of the contributions of the Nŝ scalar fields representing the scattering coefficients of the sample, as calculated below. The backprojection is defined in Code 1 (Ref. [36] return y Data size. The reconstruction problem considered in [3] fits in memory on modern GPUs. The measured data consists of 46 tilt angles, 50 rotation angles, and 100×144 unit cells. Measuring along 8ĝ i vectors, the total number of measured unit cells equals 46 × 50 × 100 × 144 × 8 ≈ 25 × 10 6 , which requires approximately 1 GB when stored in 32 bit precision. The reconstruction volume consists of 44×71×71 voxels, repeated for each of Nŝ = 7 scattering directions. In total, it requires roughly 6 MB to store in 32 bit precision. The size of the scaling matrix ν is negligible in comparison. Modern data center GPUs range in memory size from 16 GB to 80 GB. Therefore, it is possible to run an iterative SIRT reconstruction of the full problem on GPU. Benchmarks comparing the performance on GPU versus CPU are provided in Section 6.

Experimental data
In this section, we show reconstructions of experimental data acquired using the standard circular cone beam and single-axis parallel beam trajectories, as well as a reconstruction of an X-ray scattering tensor tomography dataset. The reconstructions have been computed using the algorithms implemented in the separate ts_algorithms package. Circular cone beam. A laboratory micro-CT dataset of a bell pepper was acquired at the FleX-ray laboratory [41] at the CWI, Amsterdam, The Netherlands. A polychromatic microfocus X-ray point source with tube voltage and power of 90 kV and 49.5 W was used. The data consisted of 3600 projection images of 1512×1912 pixels, acquired over a 360°rotation. A reconstruction was computed on a grid of 1512×1912×1912 voxels using FDK, a backprojection-type algorithm [42]. An axial slice of the reconstruction is shown in Fig. 9(a). Single axis parallel beam. A 3D micro-tomography dataset of a fuel cell from the publicly available TomoBank [43] was used. This dataset (#81) was acquired at the TOMCAT beamline at Data size. The reconstruction problem considered in [3] fits in memory on modern GPUs. The measured data consists of 46 tilt angles, 50 rotation angles, and 100×144 unit cells. Measuring along 8ĝ i vectors, the total number of measured unit cells equals 46 × 50 × 100 × 144 × 8 ≈ 25 × 10 6 , which requires approximately 1 GB when stored in 32 bit precision. The reconstruction volume consists of 44×71×71 voxels, repeated for each of Nŝ = 7 scattering directions. In total, it requires roughly 6 MB to store in 32 bit precision. The size of the scaling matrix ν is negligible in comparison. Modern data center GPUs range in memory size from 16 GB to 80 GB. Therefore, it is possible to run an iterative SIRT reconstruction of the full problem on GPU. Benchmarks comparing the performance on GPU versus CPU are provided in Section 6.

Experimental data
In this section, we show reconstructions of experimental data acquired using the standard circular cone beam and single-axis parallel beam trajectories, as well as a reconstruction of an X-ray scattering tensor tomography dataset. The reconstructions have been computed using the algorithms implemented in the separate ts_algorithms package. Circular cone beam. A laboratory micro-CT dataset of a bell pepper was acquired at the FleX-ray laboratory [41] at the CWI, Amsterdam, The Netherlands. A polychromatic microfocus X-ray point source with tube voltage and power of 90 kV and 49.5 W was used. The data consisted of 3600 projection images of 1512×1912 pixels, acquired over a 360°rotation. A reconstruction was computed on a grid of 1512×1912×1912 voxels using FDK, a backprojection-type algorithm [42]. An axial slice of the reconstruction is shown in Fig. 9(a).
Single axis parallel beam. A 3D micro-tomography dataset of a fuel cell from the publicly available TomoBank [43] was used. This dataset (#81) was acquired at the TOMCAT beamline at the Swiss Light Source (SLS) at the Paul Scherrer Institut (PSI), Villigen, Switzerland [44]. The first 3600 projection images of 1100×1440 pixels were used to compute a reconstruction on an axial slice of 1400×1400 pixels. The reconstructions were computed using FBP (Ram-Lak filter), SIRT (200 iterations), and TV-MIN (500 iterations with λ = 2 × 10 −7 ), as shown in Fig. 9(b -d).
X-ray scattering tensor tomography. The same validation sample was used as in a previous publication [3], which was also acquired at the TOMCAT beamline. It consisted of a 4×4×4 mm 3 plastic box containing three orthogonally oriented bundles of carbon fiber with a 12 µm diameter. The pixel and resulting unit cell size was 11×11 µm 2 and 99×99 µm 2 , generating the dataset size   Fig. 10. X-ray scattering tensor tomography reconstruction of a validation sample. The sample contains three orthogonally oriented carbon fiber bundles. A reconstruction of the orientation map is shown in two axial slices.

Benchmarks
In this section, we give a demonstration of the computational speed of tomosipo. First, we compare an implementation of SIRT in tomosipo to the built-in implementation in the ASTRA Toolbox.
Using the tomosipo implementation, we also investigate the impact of storing intermediate data on the CPU rather than on the GPU. This comparison is run on the examples from Section 4 with data sizes that fit on a single GPU. We exclude DCT, as its reconstruction algorithm is out of the scope of this manuscript. We also benchmark a non-iterative algorithm on a circular cone beam dataset that does not fit on the GPU. Here, we compare the speed of the built-in FDK implementation of the ASTRA Toolbox to the FDK implementation in ts_algorithms, tomosipo's accompanying reconstruction algorithms package. We describe the algorithms, data size, and benchmark methodology. The SIRT reconstructions were computed in 50 iterations. The implementation in tomosipo used PyTorch and the ASTRA implementation used the SIRT3D_CUDA algorithm. The FDK benchmark compared the FDK implementation provided by ts_algorithms to ASTRA's built-in accumulate_FDK implementation. The dataset of the parallel beam and helical cone beam cases consisted of 768 × 768 pixel projection images acquired over 512 angles and was reconstructed on a 512 3 voxel volume. The sizes of the XSTT and circular cone beam dataset are described in Sections 4.4 and 5, respectively. The benchmarks were conducted on a dual-socket system containing 8-core Intel Xeon Silver 4110 CPUs at 2.10 GHz (Intel, Santa Clara, CA, USA) with 192 GB of RAM and four Nvidia GeForce GTX 1080 Ti GPUs (Nvidia, Santa Clara, CA, USA). Each benchmark was run once without measurement to minimize startup and caching effects. The mean and standard deviation of three trials are reported.
SIRT on GPU-sized problems. The results of the SIRT benchmark are shown in Figure 11. The ASTRA Toolbox and the Tomosipo implementation with intermediate data on the GPU are close in performance. They are are 2 -9× faster than the tomosipo implementation with intermediate data located on CPU memory. This indicates that CPU-GPU communication latency is non-negligible and that reconstruction algorithms benefit from being completely computed on the GPU. We note that in all three implementations the forward and backprojection are computed on the GPU using the ASTRA Toolbox. The native ASTRA SIRT implementation does not have described at the end of Section 4.4. An illustration of the validation sample and its reconstruction using tomosipo is shown in Fig. 10. The reconstructions show the orientation of the fibers after post-processing using PCA and a similar thresholding strategy as in [3]. Thresholding causes noise in the background to be suppressed, as the X-ray scattering induced by plastic container is known to be negligible.

Benchmarks
In this section, we give a demonstration of the computational speed of tomosipo. First, we compare an implementation of SIRT in tomosipo to the built-in implementation in the ASTRA Toolbox. Using the tomosipo implementation, we also investigate the impact of storing intermediate data on the CPU rather than on the GPU. This comparison is run on the examples from Section 4 with data sizes that fit on a single GPU. We exclude DCT, as its reconstruction algorithm is out of the scope of this manuscript. We also benchmark a non-iterative algorithm on a circular cone beam dataset that does not fit on the GPU. Here, we compare the speed of the built-in FDK implementation of the ASTRA Toolbox to the FDK implementation in ts_algorithms, tomosipo's accompanying reconstruction algorithms package.
We describe the algorithms, data size, and benchmark methodology. The SIRT reconstructions were computed in 50 iterations. The implementation in tomosipo used PyTorch and the ASTRA implementation used the SIRT3D_CUDA algorithm. The FDK benchmark compared the FDK implementation provided by ts_algorithms to ASTRA's built-in accumulate_FDK implementation. The dataset of the parallel beam and helical cone beam cases consisted of 768×768 pixel projection images acquired over 512 angles and was reconstructed on a 512 3 voxel volume. The sizes of the XSTT and circular cone beam dataset are described in Sections 4.4 and 5, respectively. The benchmarks were conducted on a dual-socket system containing 8-core Intel Xeon Silver 4110 CPUs at 2.10 GHz (Intel, Santa Clara, CA, USA) with 192 GB of RAM and four Nvidia GeForce GTX 1080 Ti GPUs (Nvidia, Santa Clara, CA, USA). Each benchmark was run once without measurement to minimize startup and caching effects. The mean and standard deviation of three trials are reported.
SIRT on GPU-sized problems. The results of the SIRT benchmark are shown in Fig. 11. The ASTRA Toolbox and the Tomosipo implementation with intermediate data on the GPU are close in performance. They are are 2 -9× faster than the tomosipo implementation with intermediate data located on CPU memory. This indicates that CPU-GPU communication latency is non-negligible and that reconstruction algorithms benefit from being completely computed on the GPU. We note that in all three implementations the forward and backprojection are computed on the GPU using the ASTRA Toolbox. The native ASTRA SIRT implementation does not have an option to store intermediate data on CPU memory. Fig. 11. Comparison of reconstruction times using SIRT on a GPU-sized problem and using FDK on a lab-CT-sized problem. The SIRT implementations are compared on a parallel, helical and X-ray scattering tensor tomography (XSTT) acquisition geometry. The XSTT reconstruction cannot be implemented using the built-in ASTRA SIRT API. Because the FDK dataset is too large, intermediate data cannot be stored on the GPU, and the ASTRA implementation is compared to a tomosipo implementation that performs the filtering step on the CPU. FDK on a lab-CT-sized problem. The FDK dataset is too big to fit in GPU memory. In tomosipo's implementation, the filtering step is performed on the CPU and the computation of the backprojection on chunks of projection data is distributed over multiple GPUs. The FDK implementation in the ASTRA Toolbox, on the other hand, first distributes chunks of projection data over available GPUs and performs the filtering and backprojection in a single step on each GPU. Figure 11 shows the results of the FDK benchmark using one and four GPUs. Using one GPU, tomosipo's implementation of FDK is faster than ASTRA's. This can be attributed to fast filtering on the CPU, which is implemented using the Fast Fourier Transform provided by PyTorch and is approximately as fast as filtering on a single GPU. Using four GPUs, the run times of both implementations are reduced, but the ASTRA implementation comes out ahead. When four GPUs are available, the ASTRA implementation distributes the computation of the filter step over four GPUs, whereas the tomosipo implementation still computes the filtering step on the same amount of CPU cores.
The results show that a naive implementation of an iterative algorithm in tomosipo is not necessarily slower than a native implementation in the ASTRA Toolbox. In addition, the results illustrate the substantial negative impact that CPU-GPU communication has on reconstruction speed. Finally, the FDK results illustrate the benefits of interoperability with fast array libraries, but highlight the need for effective APIs to address multi-device streaming computation.

Discussion and conclusion
In short, tomosipo provides the expressive power to quickly and naturally define complex geometries, thereby unlocking the flexibility provided by the ASTRA Toolbox. We have demonstrated the ease of making common adjustments to an acquisition geometry, such as changing the center of rotation. In addition, the design and implementation of more complex geometries, such as the demonstrated X-ray diffraction and scattering setups, is made considerably easier by using tomosipo, especially compared to entering the formulas for all directional vectors manually. Reconstructions of real-world data from synchrotron and laboratory micro-CT sources are shown, computed using several common reconstruction algorithms. Finally, bechmarks demonstrate that the package enables the user to write fully GPU-accelerated reconstruction algorithms in Python whose speed is on par with native implementations. Because of tomosipo's interoperability with GPU-accelerated Python array libraries, intermediate results can remain on the GPU, avoiding the latency imposed by CPU-GPU communication.
The tomosipo package follows best practices. It has a comprehensive unit test suite, it is installable through the Anaconda package manager, it follows semantic versioning, it is developed in the open on GitHub, and it has extensive documentation.
Future developments are expected to go hand in hand with improvements in the ASTRA Toolbox. This includes support for curved detectors and more fine grained control of streams on the GPU, allowing for concurrency through pipelining. In addition, we intend to extend the interoperability of tomosipo's projection operator to more optimization packages. We note that the integration of tomosipo's projection operator in deep learning-based reconstruction methods using PyTorch is possible and is described in the documentation.
Compared to existing tomographic software packages, two features set tomosipo apart. First, the facilities to manipulate geometries significantly simplify defining complex acquisition geometries such as those in the described case studies. Although other packages including the ASTRA Toolbox and the Core Imaging Library (CIL) [14] can represent these acquisition geometries, they do not provide tools to define them. Specifically, the geometric transforms, subsampling, concatenation, and visualization features are not provided by the ASTRA Toolbox. Second, the extensible integration of tomosipo with GPU-accelerated Python array libraries provides two advantages. It enables the user to write custom reconstruction algorithms in Python that are comparable in computational efficiency to a native implementation. In addition, it enables integrating tomographic operators in deep learning-based reconstruction methods. This is technically possible using the ASTRA Toolbox, but the APIs that it exposes are designed to be wrapped by a user-friendly library, such as tomosipo.
We stress that tomosipo aims to be a building block in a larger system. Therefore, other software packages may be preferable for many purposes. Facilities for loading of various file formats, preprocessing of tomographic data, or post-processing of reconstructed images are present in TomoPy, Savu, and CIL [10,12,14]. Packages such as PyLops, CIL, and JUDI [14,16,17] provide building blocks and built-in optimization algorithms that enable rapid prototyping of variational reconstruction methods, among others. The reconstruction algorithms show-cased in this manuscript, on the other hand, are implemented in a separate package [25]. An advantage of the focused scope of tomosipo, is that it has only two required dependencies (NumPy and the ASTRA Toolbox), making it easy to install on various platforms, but contains several integrations with third-party packages, making it easy integrate into an existing system.

Research Article
Vol. 29, No. 24 / 22 Nov 2021 / Optics Express 40512 In summary, tomosipo provides scientists with an excellent tool to model and visualize complex tomographic acquisition geometries while maintaining and extending the fast reconstruction capabilities of the ASTRA Toolbox.