Reconstruction and visualization

This is a brief tutorial that demonstrates how to do a basic reconstruction and visualization of a data set in 2D and 3D, and explains some of the parameters necessary for the reconstruction to work correctly. The reconstruction only requires the usage of CPU resources. The simulated data used in the rendered version of this tutorial can be obtained using, e.g., wget with

import h5py
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
from ipywidgets import interact
from mumott.data_handling import DataContainer
from mumott.output_handling import ProjectionViewer
from mumott.output_handling.saving import dict_to_h5
from mumott.methods.basis_sets import SphericalHarmonics
from mumott.methods.projectors import SAXSProjector
from mumott.methods.residual_calculators import GradientResidualCalculator
from mumott.optimization.loss_functions import SquaredLoss
from mumott.optimization.optimizers import LBFGS
from mumott.optimization.regularizers import Laplacian
INFO:Setting the number of threads to 8
INFO:Setting numba log level to WARNING.

Loading data

First, we need to locate our input data. The path, name, and input type are defined as strings. The path can be either relative or absolute. The file name needs to include the file ending. The only allowed input format is 'h5' (which requires a specially formated hdf5 file).

Initializing the DataContainer

The input files must contain most of the necessary metadata, so creating a DataContainer object simply requires passing the file name as an argument.

data_container = DataContainer('saxstt_dataset_M.h5')
INFO:Rotation matrix generated from inner and outer angles, along with inner and outer rotation axis vectors. Rotation and tilt angles assumed to be in radians.
mumott/data_handling/ DeprecationWarning: Entry name rotations is deprecated. Use inner_angle instead.
mumott/data_handling/ DeprecationWarning: Entry name tilts is deprecated. Use outer_angle instead.
mumott/data_handling/ DeprecationWarning: Entry name offset_j is deprecated. Use j_offset instead.
mumott/data_handling/ DeprecationWarning: Entry name offset_k is deprecated. Use k_offset instead.
INFO:No sample geometry information was found. Default mumott geometry assumed.
INFO:No detector geometry information was found. Default mumott geometry assumed.

We can now take a quick look at the container object to see that the parameters make sense.



Number of projections 417
Corrected for transmission False

We can use the ProjectionViewer to check out projections. For interactive usage, it would be necessary to use the setting %matplotlib widget, but for the purpose of correctly rendering the widget in HMTL, this is left out.

p = ProjectionViewer(data_container, orientation_symmetry='transversal')

@interact(i=(0, len(data_container), 1))
def change_projection(i=0):
mumott/output_handling/ MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  self._colormap = get_cmap(colormap)
mumott/output_handling/ MatplotlibDeprecationWarning: The get_cmap function was deprecated in Matplotlib 3.7 and will be removed two minor releases later. Use ``matplotlib.colormaps[name]`` or ``matplotlib.colormaps.get_cmap(obj)`` instead.
  self._phase_colormap_lut = cm.get_cmap(self._phase_colormap, 256)
colorspacious/ RuntimeWarning: invalid value encountered in divide
  t = (C
ipywidgets/widgets/ DeprecationWarning: `ipykernel.pylab.backend_inline` is deprecated, directly use `matplotlib_inline.backend_inline`
  from ipykernel.pylab.backend_inline import flush_figures

We now need to choose a basis set (we choose SphericalHarmonics with ell_max=6), a projector (we use SAXSProjector, which is slow but runs on a CPU), a ResidualCalculator (we use GradientResidualCalculator), a loss function (we use a standard SquaredLoss), and an optimizer (we use LBFGS).

We also add a regularizer (the Laplacian) and append it to the loss function with a weight of 1e-2.

basis_set = SphericalHarmonics(ell_max=6)


Maximum "ell" 16
"ell" indices 28[0 2 ... 6 6]
"emm" indices 28[ 0 -2 ... 5 6]
Coefficient projection matrix (1, 1, 28)[[[ 1. 0. ... -0. 2.42]]]
Hash 175340c5
projector = SAXSProjector(data_container.geometry)


is_dirty 1False
hash 17e3d402
ResidualCalculator = GradientResidualCalculator(


BasisSet 1SphericalHarmonics
Projector 21SAXSProjector
Is dirty 1False
probed_coordinates (417, 8, 11, 3)947a95 (hash)
Hash 175fd7ee
loss_function = SquaredLoss(ResidualCalculator)


ResidualCalculator 1GradientResidualCalculator
use_weights 1False
preconditioner_hash 1None
residual_norm_multiplier 11
Function of residual r 1$L(r) = r^2$
Hash 1758feb0
regularizer = Laplacian()


Function of coefficients 1$R(\vec{x}) = \lambda \frac{\Vert \nabla^2 \vec{x} \Vert^2}{2}$
Hash 181df556
optimizer = LBFGS(loss_function, maxiter=10)


LossFunction 1SquaredLoss
Hash 176f47cd
options[maxiter] 110

To suppress excessive output from running the optimization, we will get the 'mumott.optimization.optimizer' logger and set the level to logging.WARNING, which will suppress INFO outputs.

results = optimizer.optimize()
100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 10/10 [00:42<00:00,  4.20s/it]

We use the basis_set to get useful output from our optimization coefficients in the form of a dictionary. For example, we plot the power spectrum for ell = 2.

output = basis_set.get_output(ResidualCalculator.coefficients)
plt.imshow(output['spherical_functions']['power_spectra'][:, 25, :, 1], cmap='cet_gouldian')
dict_keys(['means', 'variances', 'eigenvectors', 'r2_tensors', 'tensor_to_matrix_indices', 'eigenvalues', 'main_orientations', 'main_orientation_symmetries', 'normalized_standard_deviations', 'power_spectra', 'power_spectra_ell'])
<matplotlib.image.AxesImage at 0x7f4f4f9b3ca0>

We now use the function dict_to_h5, which will recursively save entries in dictionaries containing numpy arrays or array-like elements to h5.

dict_to_h5(dict_to_output=output, filename='output.h5', overwrite=True)
INFO:File output.h5 saved successfully!

Visualizing the output

We want to take a look at the orientations of our sample, so let’s load normalized_standard_deviations, means, and main_orientations for scale, color and orientation vector respectively.

with h5py.File(f'output.h5') as file:
    orientation = np.copy(file['spherical_functions/main_orientations'])
    scale = np.copy(file['spherical_functions/means'])
    color = np.copy(file['spherical_functions/normalized_standard_deviations'])

To do a 3D render, we use mayavi, a VTK frontend. This is a convenient package for viewing 3D data in Python, but not mandatory, as there are many other alternatives. For notebook rendering, we will use the png backend to make the render persist in the notebook; for interactive plots, ipy is preferable.

import colorcet
from matplotlib import colormaps as cm
(50, 50, 50, 3)

We want to threshold the scale and color, so we will plot a slice of each.

f, ax = plt.subplots(2, figsize=(5, 9))
image_0 = ax[0].imshow(scale[:, 25, :], cmap='cet_gouldian')
plt.colorbar(image_0, ax=ax[0])
image_1 = ax[1].imshow(color[:, 25, :], cmap='cet_fire', norm='log')
plt.colorbar(image_1, ax=ax[1])
<matplotlib.image.AxesImage at 0x7f4f4f93fd60>
<matplotlib.image.AxesImage at 0x7f4f7ac8d730>

Based on these images, we will choose 0.4 and 1.0 as our cutoffs for the scale, and 0.2 and 1 as our cutoffs for the color. We check our reconstruction out using matplotlib.pyplot.quiver, which can also be a useful alternative to 3D rendering.

scale_cutoff = (0.4, 1.0)
color_cutoff = (0.2, 1.)
test_scale = scale[5:45, 25, 5:45].clip(0, scale_cutoff[1])
test_scale[test_scale < scale_cutoff[0]] = np.nan
test_scale -= scale_cutoff[0]
test_scale /= (scale_cutoff[1] - scale_cutoff[0])
test_color = color[5:45, 25, 5:45].clip(*color_cutoff)
f, ax = plt.subplots(1, figsize=(9, 8))
x, y = np.meshgrid(np.arange(5, 45),
                   np.arange(5, 45),
test_orientation_x = orientation[5:45, 25, 5:45, 0] * test_scale
test_orientation_y = orientation[5:45, 25, 5:45, 2] * test_scale * np.sign(test_orientation_x)
test_orientation_x = abs(test_orientation_x)
q = ax.quiver(x, y, test_orientation_x, test_orientation_y,
              test_color, pivot='mid', angles='xy',
              headwidth=1, headlength=0, scale=30, width=0.005,
              cmap='cet_rainbow4', clim=color_cutoff)
plt.colorbar(q, ax=ax)
<matplotlib.quiver.Quiver at 0x7f4f7ada3eb0>

This concludes the tutorial.