Geometry Description

In X-ray CT, the trajectory is a description of the measurement process. It defines from where each projection is taken and how the measurements are acquired. As such the trajectory is an essential part of the reconstruction process, playing a cruicial role in the overall quality of the reconstruction. This guide explores the essentials of defining trajectories using elsa, providing a basic building block for tomographic reconstructions.

Circular Trajectory for phantom data

There are two distinct use cases. First, if you are working with phantom data (e.g. the Shepp-Logan phantom), you need to create a trajectory, to create the measurements. Second, if you are working with CT data, you need to create the trajectory for that specific data. The trajectory is created in the same fashion, however, some of the code around will be slightly different. We will explore both ways in this guide.

If you are working with phantom data, you can access the DataDescriptor the following way:

# Create phantom
size = np.array([512] * 3)
phantom = elsa.phantoms.modifiedSheppLogan(size)

# Access the descriptor
volume_descriptor = phantom.getDataDescriptor()

Then you need to define some parameters describing the trajectory. For now, we’ll assume a circular trajectory. A circular trajectory is defined by a list of angles, the distance from the source to the center of rotation, and the distance from the center or rotation to the detector. Optionally, you can define an offset of the detector, an offset of the center of rotation and the size and spacing of the detector.

Here, is an example for a very simple circular trajectory with a 360 degree angular resolution of 1° degree increments:

angles = np.linspace(0, 360, 360)
dist_source_origin = 5000
dist_origin_detector = 700

sino_descriptor = elsa.CircleTrajectoryGenerator.trajectoryFromAngles(
    angles,
    volume_descriptor,
    dist_source_origin,
    dist_origin_detector,
)

The size of the detector is calculated automatically as int(np.max(vol.shape) * np.sqrt(2)), and it’s spacing is all ones. This ensure an good coverage of the edges of the volume.

Trajectory for CT data

In the case of real measurements, you need to know the setup the data was acquired with. Usually, that information is stored in a file, or as metadata. In this quickstart guide, we will setup a trajectory that will work with the walnut data from the Finish Inverse Problem Society, which can be found here. Hence, all values are specific to this dataset, but it should be straight forward to adapt to other datasets.

First, we assume you have loaded the data as an array of shape (nangles, u_size, v_size). Then, the dataset consists of 721 projection images, which were acquired in a circular trajectory over the complete circle with an equal distance of 0.5° increments. Next, we need the distance between the source and the center of rotation, which was 210.66mm. The distance between the center of rotation and detector is not given directly, but with the distance from source to detector it can be calculated as 553.74mm - 210.66mm = 343.08mm. Each detector pixel is of size 0.050mm x 0.050mm.

Define the volume

Okay, with that out of the way, wef first need to create the volume descriptor. In theory, we can define it any size we want it to. If you have a reference reconstruction, you can use these parameters. Here, we will use the magnification of the system, to determine a suitable size:

# Given parameters about distance
dist_source_detector = 553.74
dist_source_origin = 210.66
dist_origin_detector = dist_source_detector - dist_source_origin

# Shape and spacing of the detector
det_shape = (2240, 2368)
det_spacing = np.asarray([0.050, 0.050])

# Calculate system magnification
magnification = dist_source_detector / dist_source_origin

# Volume is as large as the largest detector axis
vol_shape = [np.max(det_shape)] * 3

# Spacing is based on the magnification
vol_spacing = detectorSpacing[0] / magnification

# Descriptor for the size and spacing of the volume
volume_descriptor = elsa.VolumeDescriptor(
    vol_shape,
    vol_spacing,
)

Note

For this dataset, all sizes are given in millimeters. That is quite typical but not necessarily the case. elsa by itself is completely unitless. As long as all units are in the same unit, the results will be the same. I.e. using miliimeters vs centimeters is not a problem.

This is cruicial to cover a wide variety of different setups, where some might Micro CT setups, with only a few µm spacing. Assuming a too large default unit, might create floating point problems in such scenarios. Hence, we provide this flexibility to the user.

Create the circular trajectory

With that out of the way, we can create the circular trajectory.

angles = np.linspace(0, 360, 721)

# Create a circular trajectory for the FIPS walnut data
sino_descriptor = elsa.CircleTrajectoryGenerator.trajectoryFromAngles(
    angles,
    volume_descriptor,
    dist_source_origin,
    dist_origin_detector,
    [0, 0],     # No offset of detector
    [0, 0, 0],  # No offset of center of rotation
    det_shape,
    det_spacing
)

Finally, we can create the DataContainer with the measurements:

sinogram = elsa.DataContainer(measurements, sino_descriptor)

And that is it, with that you have a sinogram, and perform any iterative reconstruction algorithm from our framework. See the next chapter on how to do that.

Arbitrary Trajectory

The circular trajectory is one of the most commonly used trajectories. However, if not, we got you covered! In this overview, we create a trajectory with a non-standard rotation axis. We’ll be using scipy’s rotation to help us, so we need to import it.

import numpy as np
import pyelsa as elsa
from scipy.spatial.transform import Rotation as R

Next, for this quickstart, we’ll define some example size and distances for the volume and detector

vol_shape = np.array([64] * 3)
vol_spacing = np.array([1.] * 3)
det_shape = np.array([100] * 3)
det_spacing = np.array([1.] * 2)

dist_source_origin = 1000.
dist_origin_detector = 100.

Okay, next we need to define the rotation axis. In this case, I’ve opted for a slightly tilted rotation axis. This is common for e.g. laminography. With a tilt, it can be possible to move object closer to the source, thus increasing the magnification. And objects might not conceal each other. But for these reasons, we choose a rotation axis, that slightly tilts towards the source:

rot_dir = np.array([0, 1, -0.3])
rot_dir = rot_dir / np.linalg.norm(rot_dir)

Note

We can take any rotation matrix, here we only opted for this specific example.

In this case, we’ll use our constructor for geometry, that takes in a rotation matrix. Here, the scipy rotation is used to generate the matrices:

angles = np.linspace(0, 360, 360, endpoint=False)
rot_matrices = [R.from_rotvec(
    angle * rot_dir, degrees=True).as_matrix() for angle in angles]

We create angles in 1° increments, and then compute the associated rotation matrix. The rotation matrix, is now a rotation matrix around the given direction. As a last step, we need to create a list of Geometry. We’ll be using the following way:

geom = []
for rot_mat in rot_matrices:
    g = elsa.Geometry(
        dist_source_origin, dist_origin_detector,
        vol_shape, det_shape,
        rot_mat,
        vol_spacing,  # optional, (default [1, 1, 1])
        det_spacing,  # optional, (default [1, 1])
        [0, 0],       # optional, detector shift (default [0, 0])
        [0, 0, 0]     # optional, center of rotation offset (default [0, 0, 0])
    )
    geom.append(g)

Again, you can see the default parameters explicitly, such that you can easily use them if you need them!

sino_desc = elsa.PlanarDetectorDescriptor(det_shape, geom)

The last step is only to create a flat detector with the associated geometry. That’s it, you’ve created a non-standard trajectory.