Python Quickstart

To deomstrate the Python bindings of your framework, we’ll use a simple example for Computed Tomography. This is pretty direct translation of the C++ quickstart example.

2D example

First, we import NumPy, elsa and matplotlib.

import numpy as np
import matplotlib.pyplot as plt

import pyelsa as elsa

The first thing, we have to do is set up some phantom data. Working with phantom data is quite useful for many applications, as we have a ground truth and compute error norms.

You can plot any data from elsa as follows (assume you got some DataContainer called dc):

plt.imshow(np.asarray(dc), cmap="gray")
plt.show()

Next, we have to define a trajectory. The trajectory defines the position of the source and the detector. In our case, this is a fan-beam geometry setup.

numAngles = 720
arc = 360
DSD = size[0] * 100
DSO = size[0] * 2
angles = np.linspace(0, arc, numAngles, endpoint=False)
sinoDescriptor = elsa.CircleTrajectoryGenerator.trajectoryFromAngles(
    angles, phantom.getDataDescriptor(), DSD, DSO)

The trajectory has 720 positions (numAngles) over the whole circle (arc). The two distances DSO and DSD are the distance from the source to the center of the volume and the distance from the center of the volume to the detector. In this case we make the distance from the center to detector and source dependent on the size of the volume. The values provided here are okay defaults. But usually, if you work with real data, these are of course given for your by your setup.

It’s usually a good idea to move the source further away from the volume, and the detector closer to it. This way also the edges of the phantom are covered nicely, and no artifacts appear during reconstruction (but don’t trust me, try it!).

Now, we need to create the sinogram. The sinogram is the result of projecting the phantom.

projector = elsa.JosephsMethod(phantom.getDataDescriptor(), sinoDescriptor)
sinogram = projector.apply(phantom)

The projector is the approximation of the Radon transform. Our current projectors are implemented based on ray traversals.

With this setup, we finally have everything to setup the tomographic problem.

h = elsa.IndicatorNonNegativity(phantom.getDataDescriptor())
solver = elsa.APGD(projector, sinogram, h)
reconstruction = solver.solve(50)

Without getting to into the math behind it (checkout the paper, if you want to dive into it), this sets up an optimization problem. We want to find an solution to this problem In this case, we’ll find it it using the Conjugate gradient method, or CGLS for short. It’s an iterative algorithm to solve. As this is still a quite simple example, we don’t need to let it run for too many iterations.

Now, you have the reconstruction! In the best case, this should already look quite similar to the original phantom. We can have a look at the difference and the L2-Norm:

plt.imshow(np.array(reconstruction), cmap="gray")
plt.colorbar()
plt.show()