Reconstruction using pipelines¶
This tutorial demonstrates the pipelines for simultaneous iterative reconstruction technique (SIRT), the modular iterative tomographic reconstruction algorithm (MITRA), as well as discrete directions (DD), 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
[1]:
import warnings
from mumott.data_handling import DataContainer
from mumott.pipelines import run_sirt, run_mitra, run_discrete_directions
from mumott.methods.utilities.grids_on_the_sphere import football
from mumott.methods.basis_sets import NearestNeighbor, GaussianKernels
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), and os.environ["OPENBLAS_NUM_THREADS"] = f"{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 – simultaneous iterative reconstruction 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.
[2]:
%%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.74it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 100/100 [00:12<00:00, 8.05it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 500/500 [01:02<00:00, 8.02it/s]
CPU times: user 1min 10s, sys: 9.16 s, total: 1min 20s
Wall time: 1min 19s
[3]:
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 result 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 and can be straightforwardly extended to tensor-tomography using DD, better reconstructions can be obtained by using regularization, available with MITRA. 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.
[4]:
%%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, 8.11it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 20/20 [00:02<00:00, 8.13it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 8.16it/s]
CPU times: user 9.12 s, sys: 882 ms, total: 10 s
Wall time: 9.98 s
[5]:
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.
[6]:
%%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:03<00:00, 7.90it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.88it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.94it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.91it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.78it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.86it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.85it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.84it/s]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [00:03<00:00, 7.81it/s]
CPU times: user 38.4 s, sys: 4.03 s, total: 42.4 s
Wall time: 42.3 s
[7]:
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).
[8]:
%%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 [01:00<00:00, 2.02s/it]
CPU times: user 1min 33s, sys: 10min 41s, total: 12min 15s
Wall time: 1min 5s
To analyze the results we here define a convenience function that allows us to extract mean, standard deviation, and orientation data from the reconstruction.
[9]:
def get_tensor_properties(basis_set, reconstruction):
if type(basis_set) is not NearestNeighbor:
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, ...])
eigenvectors = spherical_functions['eigenvectors']
else:
mean = reconstruction.mean(-1)[28]
std = np.std(reconstruction, axis=-1)[28]
eigenvectors = basis_set.get_output(reconstruction)['eigenvectors']
orientation = np.arctan2(-eigenvectors[28, ..., 1, 0],
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
[10]:
mean, std, orientation, orientation_alpha = get_tensor_properties(basis_set, reconstruction)
Next we plot the results.
[11]:
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()])
[12]:
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).
[13]:
%%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:23<00:00, 2.79s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [01:21<00:00, 2.73s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [01:19<00:00, 2.64s/it]
CPU times: user 5min 35s, sys: 32min 44s, total: 38min 20s
Wall time: 4min 13s
[14]:
reconstructions = [r['result']['x'] for r in result]
basis_sets = [r['basis_set'] for r in result]
[15]:
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).
[16]:
%%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:22<00:00, 2.75s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [01:27<00:00, 2.91s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 30/30 [01:23<00:00, 2.78s/it]
CPU times: user 5min 48s, sys: 34min 11s, total: 40min
Wall time: 4min 22s
[17]:
reconstructions = [r['result']['x'] for r in result]
basis_sets = [r['basis_set'] for r in result]
[18]:
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.
DD - discrete directions¶
Rather than regularizing by imposing constraints on the whole solution, we can place the measured reciprocal space map components into bins on the sphere, and reconstruct each bin independently. This is what the DD pipeline does. Each reconstruction is that solved by SIRT. The drawback of this method is that it is not compatible with constraint-based regularization. However, due to the inherent regularization properties of SIRT, smooth reconstructions can be obtained by using a small number
of iterations in the reconstruction.
In the following cells, we run the reconstruction and plot the result. The reconstruction shown uses the same grid as the GaussianKernels of the MITRA pipelines, for easy comparison.
[19]:
%%time
theta, phi = GaussianKernels(grid_scale=5).grid
directions = np.concatenate(((np.sin(theta) * np.cos(phi))[:, None], (np.sin(theta) * np.sin(phi))[:, None], np.cos(theta)[:, None]), axis=1)
result = []
maximum_iterations = [10, 100, 300]
for maxiter in maximum_iterations:
result.append(run_discrete_directions(data_container, directions, maxiter=maxiter, use_gpu=True))
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 72/72 [02:11<00:00, 1.82s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 72/72 [02:22<00:00, 1.98s/it]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 72/72 [05:17<00:00, 4.41s/it]
CPU times: user 9min 31s, sys: 20 s, total: 9min 51s
Wall time: 9min 51s
[20]:
reconstructions = [r['result']['x'] for r in result]
basis_set = NearestNeighbor(directions)
We calculate a number of scalar quantities from the tensor tomography solution for plotting.
[21]:
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_set, 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'$maxiter = ' + f'{maximum_iterations[i]}' + 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)
From this reconstruction, we can see that about 100 iterations is a good number, and the quality of the reconstruction as far as can be assessed from the mean, standard deviation, and orientations of this slice, is comparable to that of the regularized GaussianKernels reconstruction.
[ ]: