Reconstruction using pipelines

This tutorial demonstrates the pipelines for spherical integral geometric tensor tomography (SIRT) and the modular iterative tomographic reconstruction algorithm (MITRA), and how they can be used to get started with a reconstruction quickly.

The dataset of trabecular bone used in this tutorial is available via Zenodo, and can be downloaded either via a browser or the command line using, e.g., wget or curl,

wget https://zenodo.org/records/10074598/files/trabecular_bone_9.h5

import warnings
from mumott.data_handling import DataContainer
from mumott.pipelines import run_sirt, run_mitra
from mumott.optimization.regularizers import HuberNorm, TotalVariation
import matplotlib.pyplot as plt
import colorcet
import logging
import numpy as np
data_container = DataContainer('trabecular_bone_9.h5')

INFO:Setting the number of threads to 8. If your physical cores are fewer than this number, you may want to use numba.set_num_threads(n) to set the number of threads to the number of physical cores n.
INFO:Setting numba log level to WARNING.
INFO:Inner axis found in dataset base directory. This will override the default.
INFO:Outer axis found in dataset base directory. This will override the default.
INFO:Rotation matrices were loaded from the input file.
INFO:Sample geometry loaded from file.
INFO:Detector geometry loaded from file.

SIRT – spherical integral geometric tensor tomography

One of the most straightforward ways to obtain a good reconstruction of transmission data is by using SIRT, which is essentially a preconditioned least-squares formulation of tomography. Because of this preconditioning, the system can be relatively easily solved using gradient descent, and only requires a single parameter, the number of iterations, to be adjusted.

For illustration, in the following cell we run a SIRT reconstruction three times but using a different number of iterations.

%%time
logger = logging.getLogger('mumott.optimization.optimizers.gradient_descent')
logger.setLevel(logging.WARNING)

result_1 = run_sirt(data_container,
                    use_gpu=True,
                    maxiter=10)['result']['x']
result_2 = run_sirt(data_container,
                    use_gpu=True,
                    maxiter=100)['result']['x']
result_3 = run_sirt(data_container,
                    use_gpu=True,
                    maxiter=500)['result']['x']

100%|███████████████████████████████████████████| 10/10 [00:01<00:00,  7.11it/s]
100%|█████████████████████████████████████████| 100/100 [00:15<00:00,  6.36it/s]
100%|█████████████████████████████████████████| 500/500 [01:19<00:00,  6.26it/s]
CPU times: user 1min 31s, sys: 8.14 s, total: 1min 40s
Wall time: 1min 40s

imshow_kwargs=dict(vmin=0, vmax=0.05, cmap='cet_gouldian')
fig, ax = plt.subplots(1, 3, figsize=(11.6, 3.2), dpi=140, sharey=True)
plt.subplots_adjust(wspace=0)
im = ax[0].imshow(result_1[28, ..., 0], **imshow_kwargs);
im = ax[1].imshow(result_2[28, ..., 0], **imshow_kwargs);
im = ax[2].imshow(result_3[28, ..., 0], **imshow_kwargs);
ax[0].set_title('10 iterations')
ax[1].set_title('100 iterations')
ax[2].set_title('500 iterations')
plt.colorbar(im, ax=ax);

The results show that 10 iterations are insufficient to reach convergence, where the resul after 500 iterations indicates overfitting. We can thus conclude that the optimal number of iterations is somewhere between 100 and 500.

MITRA - modular iterative tomographic reconstruction algorithm

While SIRT is easy to get started with, better reconstructions can be obtained by leveraging the options available with MITRA, including regularizers. By default, MITRA also uses a Nestorov gradient in its descent optimization, which greatly improves the rate of convergence. It also uses the same preconditioner as SIRT.

Scalar reconstruction

In the following cell we first run the “standard” MITRA pipeline (without regularizers) for three different numbers of iterations. The results illustrate the faster convergence enabled by the Nestorov optimizer used here. By default (use_absorbances=True) the reconstruction is run using the absorbance data, corresponding to the reconstructiong of a scalar field.

%%time
result_1 = run_mitra(data_container,
                     use_gpu=True,
                     maxiter=10)['result']['x']
result_2 = run_mitra(data_container,
                     use_gpu=True,
                     maxiter=20)['result']['x']
result_3 = run_mitra(data_container,
                     use_gpu=True,
                     maxiter=30)['result']['x']

100%|███████████████████████████████████████████| 10/10 [00:01<00:00,  6.39it/s]
100%|███████████████████████████████████████████| 20/20 [00:03<00:00,  6.12it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.27it/s]
CPU times: user 11.7 s, sys: 931 ms, total: 12.6 s
Wall time: 12.6 s

fig, ax = plt.subplots(1, 3, figsize=(11.6, 3.2), dpi=140, sharey=True)
plt.subplots_adjust(wspace=0)
im = ax[0].imshow(result_1[28, ..., 0], **imshow_kwargs);
im = ax[1].imshow(result_2[28, ..., 0], **imshow_kwargs);
im = ax[2].imshow(result_3[28, ..., 0], **imshow_kwargs);
ax[0].set_title('10 iterations')
ax[1].set_title('20 iterations')
ax[2].set_title('30 iterations')
plt.colorbar(im, ax=ax);

Scalar reconstruction with regularization

This is only the start. We can also attach regularizers to MITRA, further improving the reconstruction. We will illustrate this with the HuberNorm and TotalVariation regularizers.

The HuberNorm regularizer suppresses noise while converging more easily than a traditional \(L_1\) norm. The latter is achieved by switching from the \(L_1\) to the \(L_2\) norm for residuals below a threshold value that is set by the delta parameter.

The TotalVariation regularizer adds smoothness while preserving edges, and like the HuberNorm, is spliced with a more easily converging function below a threshold set by its respective delta parameter.

Here, we set the delta parameter to 1e-2, representing a typical “small” value in the reconstruction.

%%time
def get_regularizers(alpha, beta):
    return [dict(name='hn', regularizer=HuberNorm(delta=1e-2), regularization_weight=alpha),
            dict(name='tv', regularizer=TotalVariation(delta=1e-2), regularization_weight=beta)]
alpha_values = np.geomspace(2e-3, 2e-5, 3)
beta_values = np.geomspace(1e-4, 1e-6, 3)
result = []
values = []
for alpha in alpha_values:
    for beta in beta_values:
        result.append(run_mitra(data_container,
                                use_gpu=True,
                                maxiter=30,
                                Regularizers=get_regularizers(alpha, beta))['result']['x'])
        values.append((alpha, beta))

100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.88it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.54it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.72it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.80it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.65it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.78it/s]
100%|███████████████████████████████████████████| 30/30 [00:05<00:00,  5.97it/s]
100%|███████████████████████████████████████████| 30/30 [00:04<00:00,  6.24it/s]
100%|███████████████████████████████████████████| 30/30 [00:04<00:00,  6.18it/s]
CPU times: user 49 s, sys: 3.82 s, total: 52.8 s
Wall time: 52.8 s

fig, ax = plt.subplots(3, 3, figsize=(11., 9.4), dpi=140, sharey='row', sharex='col')
plt.subplots_adjust(wspace=0.00, hspace=0.)
for i in range(3):
    ax[i, 0].set_ylabel(r'$\alpha = {:0.0e}$'.format(values[i*3][0]))
    for j in range(3):
        im = ax[i, j].imshow(result[i*3 + j][28, ..., 0], **imshow_kwargs);
        ax[2, j].set_xlabel(r'$\beta = {:0.0e}$'.format(values[i*3 + j][1]))
plt.colorbar(im, ax=ax, shrink=0.5, pad=0.0);

The top row has large parts of the expected reconstruction blanked out due to over-regularization. Conversely the bottom row has barely any noise suppression at all. Similarly, the highly regularized cases to the left become blurry, while the under-regularized cases to the right become grainy.

We are thus looking for a choice in the middle of these extremes. Based on this visual inspection, we should select values of \(\alpha = 2e-4\) and \(\beta = 1e-5\) for our regularization parameters.

Application to tensor tomography

MITRA is not limited to scalar reconstruction but be applied to tensor tomography as well. By default it will use GaussianKernels to represent functions on the sphere, and apply tensor SIRT weights.

For demonstration, we run the MITRA pipeline using the “full” data from the data container (i.e., we do not limit us to the absorbances).

%%time
results = run_mitra(data_container,
                    use_gpu=True,
                    maxiter=30,
                    use_absorbances=False)

reconstruction = results['result']['x']
basis_set = results['basis_set']

100%|███████████████████████████████████████████| 30/30 [00:54<00:00,  1.83s/it]
CPU times: user 47.8 s, sys: 11.4 s, total: 59.2 s
Wall time: 59.2 s

To analyze the results we here define a convenience function that allows us to extract mean, standard deviation, and orientation data from the reconstruction.

def get_tensor_properties(basis_set, reconstruction):
    with warnings.catch_warnings():
        warnings.filterwarnings('ignore', r'invalid value encountered in divide')
        spherical_functions = basis_set.get_output(reconstruction)['spherical_harmonic_analysis']['spherical_functions']
    mean = spherical_functions['means'][28, ...]
    std = np.sqrt(spherical_functions['variances'][28, ...])

    orientation = np.arctan2(-spherical_functions['eigenvectors'][28, ..., 1, 0],
                             spherical_functions['eigenvectors'][28, ..., 2, 0])
    orientation[orientation < 0] += np.pi
    orientation[orientation > np.pi] -= np.pi
    orientation *= 180 / np.pi
    orientation_alpha = np.clip(2 * std / std.max(), 0, 1)
    return mean, std, orientation, orientation_alpha

mean, std, orientation, orientation_alpha = get_tensor_properties(basis_set, reconstruction)

Next we plot the results.

colorbar_kwargs = dict(orientation='horizontal', shrink=0.75, pad=0.1)
mean_kwargs = dict(vmin=0, vmax=800, cmap='cet_gouldian')
std_kwargs = dict(vmin=0, vmax=400, cmap='cet_fire')
orientation_kwargs = dict(vmin=0, vmax=180, cmap='cet_CET_C10')

def config_cbar(cbar):
    cbar.set_ticks([i for i in np.linspace(0, 180, 5)])
    cbar.set_ticklabels([r'$' + f'{int(v):d}' + r'^{\circ}$' for v in bar.get_ticks()])

fig, ax = plt.subplots(1, 3, figsize=(10.3, 4.6), dpi=140, sharey=True)

im0 = ax[0].imshow(mean, **mean_kwargs);
im1 = ax[1].imshow(std, **std_kwargs);
im2 = ax[2].imshow(orientation, alpha=orientation_alpha, **orientation_kwargs);

ax[0].set_title('Mean amplitude')
ax[1].set_title('Anisotropic ampl.')
ax[2].set_title('Orientation')

plt.subplots_adjust(wspace=0)

plt.colorbar(im0, ax=ax[0], **colorbar_kwargs)
plt.colorbar(im1, ax=ax[1], **colorbar_kwargs)
bar = plt.colorbar(im2, ax=ax[2], **colorbar_kwargs)
config_cbar(bar)

Tensor reconstruction with regularization

We can regularize our reconstruction in a very similar manner as for the scalar case. To keep the plots clear, we start by only running over beta and thus limit ourselves to the TotalVariation regularizer. We set delta=1 as there is no single optimal value correspoding to a “small” number; one could reasonably set it anywhere between 0.8 and 80. A larger number will make the reconstruction converge more easily, but a smaller number will be more effective for removing noise. As a general guideline, one thus would like to use the smallest number that for which convergence can be achieved (with reasonable effort).

%%time
def get_regularizers(alpha=0, beta=0):
    return [dict(name='hn', regularizer=HuberNorm(delta=1), regularization_weight=alpha),
            dict(name='tv', regularizer=TotalVariation(delta=1), regularization_weight=beta)]
beta_values = np.geomspace(0.5, 5, 3)
result = []
values = []
for beta in beta_values:
    result.append(run_mitra(data_container,
                            use_gpu=True,
                            maxiter=30,
                            use_absorbances=False,
                            Regularizers=get_regularizers(beta=beta)))
    values.append(beta)

100%|███████████████████████████████████████████| 30/30 [01:18<00:00,  2.62s/it]
100%|███████████████████████████████████████████| 30/30 [01:17<00:00,  2.57s/it]
100%|███████████████████████████████████████████| 30/30 [01:11<00:00,  2.40s/it]
CPU times: user 3min 7s, sys: 47.7 s, total: 3min 55s
Wall time: 3min 55s

reconstructions = [r['result']['x'] for r in result]
basis_sets = [r['basis_set'] for r in result]

colorbar_kwargs = dict(orientation='horizontal', shrink=0.75, pad=0.033)
fig, axes = plt.subplots(3, 3, figsize=(9.4, 11.6), dpi=140, sharey=True, sharex=True)
for i in range(3):
    ax = axes[i]
    mean, std, orientation, orientation_alpha = get_tensor_properties(basis_sets[i], reconstructions[i])
    plt.subplots_adjust(wspace=0, hspace=0)
    im0 = ax[0].imshow(mean, **mean_kwargs);
    im1 = ax[1].imshow(std, **std_kwargs);
    im2 = ax[2].imshow(orientation, alpha=orientation_alpha, **orientation_kwargs);
    if i == 0:
        ax[0].set_title('Mean amplitude')
        ax[1].set_title('Anisotropic ampl.')
        ax[2].set_title('Orientation')
    ax[0].set_ylabel(r'$\beta = ' + f'{beta_values[i]:.2f}' + r'$')
    if i == 2:
        plt.colorbar(im0, ax=axes[:, 0], **colorbar_kwargs);
        plt.colorbar(im1, ax=axes[:, 1], **colorbar_kwargs);
        bar = plt.colorbar(im2, ax=axes[:, 2], **colorbar_kwargs);
        config_cbar(bar)

As before with increasing beta, the reconstructions become blurrier. The comparison suggests that a suitable choice for beta should be close to the middle row.

Next we scan different values for alpha (HuberNorm) with beta=2 (TotalVariation).

%%time
alpha_values = np.geomspace(1, 40, 3)
result = []
values = []
for alpha in alpha_values:
    result.append(run_mitra(data_container,
                  use_gpu=True,
                  maxiter=30,
                  use_absorbances=False,
                  Regularizers=get_regularizers(alpha=alpha, beta=2.)))
    values.append(alpha)

100%|███████████████████████████████████████████| 30/30 [01:09<00:00,  2.33s/it]
100%|███████████████████████████████████████████| 30/30 [01:08<00:00,  2.27s/it]
100%|███████████████████████████████████████████| 30/30 [01:06<00:00,  2.21s/it]
CPU times: user 2min 46s, sys: 45.4 s, total: 3min 32s
Wall time: 3min 32s

reconstructions = [r['result']['x'] for r in result]
basis_sets = [r['basis_set'] for r in result]

fig, axes = plt.subplots(3, 3, figsize=(9.4, 11.6), dpi=140, sharey=True, sharex=True)
for i in range(3):
    ax = axes[i]
    mean, std, orientation, orientation_alpha = get_tensor_properties(basis_sets[i], reconstructions[i])
    plt.subplots_adjust(wspace=0, hspace=0)
    im0 = ax[0].imshow(mean, **mean_kwargs);
    im1 = ax[1].imshow(std, **std_kwargs);
    im2 = ax[2].imshow(orientation, alpha=orientation_alpha, **orientation_kwargs);
    if i == 0:
        ax[0].set_title('Mean amplitude')
        ax[1].set_title('Anisotropic ampl.')
        ax[2].set_title('Orientation')
    ax[0].set_ylabel(r'$\alpha = ' + f'{alpha_values[i]:.2f}' + r'$')
    if i == 2:
        plt.colorbar(im0, ax=axes[:, 0],  **colorbar_kwargs);
        plt.colorbar(im1, ax=axes[:, 1],  **colorbar_kwargs);
        bar = plt.colorbar(im2, ax=axes[:, 2],  **colorbar_kwargs);
        config_cbar(bar)

We can now see the effect of the \(L_1\) regularizer. It has a broad dampening effect that first only acts on the background, but eventually becomes strong enough to dampen the entire reconstruction.

The comparison suggests a value for alpha somewhere in the range from 6 (middle row) to 40 (bottom row). Generally, one must start with a broader sweep than we have done in this tutorial, where we already start with good guesses. In practice a regularization parameter typically sweep might span 10 to 15 orders of magnitude.