GPU-Accelerated reconstruction of holography stacks

This notebooks demonstrates the whole reconstruction process, from loading raw stacks of holograms and reference waves, to the final complex image. This notebook requires a CUDA-capable GPU. The sample imaged was a Fe3O4 nanoparticle, the reference waves were acquired in vacuum.

Note

Check the README for installation instructions including GPU support.

Packages

[1]:

from IPython import get_ipython
from pathlib import Path
import os
testdata_base_path = os.environ.get("TESTDATA_BASE_PATH")
if testdata_base_path is not None:
    testdata_base_path = Path(testdata_base_path)

[2]:

# Matplotlib widget magic
if os.environ.get("GENERATING_NOTEBOOKS"):
    get_ipython().run_line_magic("matplotlib", "inline")
else:
    get_ipython().run_line_magic("matplotlib", "widget")

[3]:

import cupy as cp

[4]:

import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import tqdm
from cupyx.scipy.ndimage import shift as shiftcp

[5]:

from libertem_holo.base.io import Results, InputData
from libertem_holo.base.reconstr import reconstruct_frame, reconstruct_bf, phase_offset_correction
from libertem_holo.base.utils import HoloParams
from libertem_holo.base.align import ImageCorrelator
from libertem_holo.base.filters import clipped
from libertem_holo.base.unwrap import phase_unwrap

Prepare input data

Before continuing, download the files called stack_ref_minus.dm4 and stack_obj_minus.dm4 from https://zenodo.org/uploads/21108532 and put it into the same folder as this notebook. For example using wget:

wget -O stack_ref_minus.dm4 'https://zenodo.org/records/21108532/files/stack_ref_minus.dm4?download=1&preview=1&token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6ImE2ZDgyOGQzLWZlMTgtNDRhMy1hOGUxLTBlOWU2NzVhZDFhMyIsImRhdGEiOnt9LCJyYW5kb20iOiIyZGU0NjIwMGQ2ZGRhNzlkNzEzODY4OWJjMGRhNGZkMSJ9.rQjUanP6uggj6JLxqUlWyjYvIC8eZP3yd1i1vPumju8roGieF7NV19hA0nA-JolhxKVyhxq-2zikUHK-3i3l7w'

wget -O stack_obj_minus.dm4 'https://zenodo.org/records/21108532/files/stack_obj_minus.dm4?download=1&preview=1&token=eyJhbGciOiJIUzUxMiJ9.eyJpZCI6ImE2ZDgyOGQzLWZlMTgtNDRhMy1hOGUxLTBlOWU2NzVhZDFhMyIsImRhdGEiOnt9LCJyYW5kb20iOiIyZGU0NjIwMGQ2ZGRhNzlkNzEzODY4OWJjMGRhNGZkMSJ9.rQjUanP6uggj6JLxqUlWyjYvIC8eZP3yd1i1vPumju8roGieF7NV19hA0nA-JolhxKVyhxq-2zikUHK-3i3l7w'

[6]:

# only for our internal testing:
if testdata_base_path is not None:
    ref_path = testdata_base_path / "dm/holo/reconstr/stack_ref_minus.dm4"
    obj_path = testdata_base_path / "dm/holo/reconstr/stack_obj_minus.dm4"
else:
    ref_path = Path(".") / "stack_ref_minus.dm4"
    obj_path = Path(".") / "stack_obj_minus.dm4"

Load hologram stacks

Data is loaded from dm3 or dm4 files directly into a ‘data’ object which contains the stack as well as metadata information such as the pixel size in meters, the total exposure time in seconds and the GMS tags.

[7]:

# Reference stack
data_ref = InputData.load_from_dm(ref_path)

# Object stack
data_obj = InputData.load_from_dm(obj_path)

[8]:

data_obj.pixelsize, data_obj.exposure_time, data_obj.shape

[8]:

(1.4522670209407806e-10, np.float64(60.0), (30, 3838, 3710))

[9]:

# Show one frame of the stack

frame = data_obj.data[0]
extent=(
    -data_obj.pixelsize*frame.shape[1]/2*1e9, data_obj.pixelsize*frame.shape[1]/2*1e9,
    -data_obj.pixelsize*frame.shape[0]/2*1e9, data_obj.pixelsize*frame.shape[0]/2*1e9,
)
plt.figure()
im = plt.imshow(
    frame, vmin=np.min(clipped(frame)), vmax=np.max(clipped(frame)),
    extent=extent, cmap='gray',
)
plt.colorbar()
plt.xlabel('nm')
plt.ylabel('nm')
plt.title("First Hologram")

[9]:

Text(0.5, 1.0, 'First Hologram')

../_images/usage_reconstruction_13_1.png

Hologram reconstruction parameters

The hologram reconstruction parameters are:

- sideband position

- sideband size

- output shape

- line filter for Fresnel fringes (length and width)

- Butterworth filter order

The parameters are stored in the ‘HoloParams’ object. A convenient function ‘from_hologram’ calculates all parameters automatically with default values, but every parameter is editable.

[10]:

# One frame is passed as argument for the parameters determination, for example the first frame of the reference stack.

holoparams = HoloParams.from_hologram(data_ref.data[0], xp=cp)

/cachedata/alex/source/modules/LiberTEM-holo/.tox/notebooks_gen/lib/python3.12/site-packages/numba_cuda/numba/cuda/cudadrv/driver.py:381: UserWarning: cuda-bindings was built for CUDA major version 13, but the NVIDIA driver only supports up to CUDA 12. Some cuda-bindings features may not work correctly. Consider updating your NVIDIA driver, or using a cuda-bindings version built for CUDA 12. (Set CUDA_PYTHON_DISABLE_MAJOR_VERSION_WARNING=1 to suppress this warning.)
  dev = Device(devnum)

[11]:

holoparams.sb_size, holoparams.sb_position, holoparams.out_shape

[11]:

(121.94465957966344, (3639.0, 3569.0), (959, 927))

If necessary, we can force a sideband size, or output shape in the holoparams if the output shape does not respect the aspect ratio of the original hologram, a warning is raised to warn the user that the pixel size will no longer be valid for both x and y dimensions, continue at your own risk

[12]:

holoparams_override = HoloParams.from_hologram(data_ref.data[0], xp=cp, sb_size=200, out_shape=(512, 512))

[13]:

holoparams_override.sb_size, holoparams_override.sb_position, holoparams_override.out_shape

[13]:

(200, (3639.0, 3569.0), (512, 512))

We can visualise the aperture created to filter the FFT:

[14]:

fig, ax = plt.subplots(ncols=2)
aperture = cp.fft.fftshift(holoparams.aperture).get()
fft = cp.fft.fftshift(cp.fft.fft2(cp.asarray(data_obj.data[0]))).get()
fft_rolled = np.roll(fft, np.array(holoparams.sb_position_int).astype(np.int64), axis=(0, 1))
fft_rolled = fft_rolled[holoparams.slice_fft]

ax[0].imshow(np.log1p(np.abs(fft))[holoparams.slice_fft], cmap='gray')
ax[0].set_title('FFT of first hologram')
ax[0].axis("off")
ax[1].imshow(aperture * np.log1p(np.abs(fft_rolled)), cmap='gray')
ax[1].set_title('Aperture over selected sideband')
ax[1].axis("off")

plt.show()

../_images/usage_reconstruction_22_0.png

Reconstruction of one frame

Before reconstructing the whole stack, it can be useful to check the reconstruction parameters by reconstructing one single frame and plotting the amplitude and phase of the complex wave.

[15]:

# First the object wave is reconstructed
wave_obj = reconstruct_frame(
    data_obj.data[0], sb_pos=holoparams.sb_position, aperture=holoparams.aperture,
    slice_fft=holoparams.slice_fft, xp=cp
)

# Then the reference wave...
wave_ref = reconstruct_frame(
    data_ref.data[0], sb_pos=holoparams.sb_position, aperture=holoparams.aperture,
    slice_fft=holoparams.slice_fft, xp=cp
)

# ... and the reference wave is removed from the object wave to reduce distorsions
wave = wave_obj / wave_ref

If you are using a GPU and the cupy library instead of numpy, the final wave must be transferred to CPU before plotting

[16]:

wave = wave.get()

The amplitude and phase images can then be extracted and plotted

[17]:

amp = np.abs(wave)
phase = np.angle(wave)

fig, ax = plt.subplots(ncols=2)
ax[0].imshow(amp, vmin=np.min(clipped(amp)), vmax=np.max(clipped(amp)))
ax[0].axis("off")
ax[0].set_title("Amplitude")
ax[1].imshow(phase)
ax[1].axis("off")
ax[1].set_title("Wrapped phase")

[17]:

Text(0.5, 1.0, 'Wrapped phase')

../_images/usage_reconstruction_28_1.png

The bright field (centerband filtering) of the hologram can also be reconstructed. The same ‘holoparams’ can be used as well.

[18]:

bf = reconstruct_bf(data_obj.data[0], aperture=holoparams.aperture_bf, slice_fft=holoparams.slice_fft, xp=cp)

[19]:

# if using cupy
bf = bf.get()

[20]:

fig, ax = plt.subplots(ncols=2, figsize=(10, 5))
ax[0].imshow(bf)
ax[0].set_title("Brightfield reconstructed from centerband", fontsize=10)
ax[0].axis('off')

ax[1].imshow(np.fft.fftshift(holoparams.aperture_bf.get()))
ax[1].set_title("Aperture used on centerband with line filters\nto remove Fresnel fringes", fontsize=10)
ax[1].axis('off')

[20]:

(np.float64(-0.5), np.float64(926.5), np.float64(958.5), np.float64(-0.5))

../_images/usage_reconstruction_32_1.png

Reconstruction stack

Now that we are satisfied with our reconstruction parameters, we can reconstruct a stack.

[21]:

from libertem_holo.base.convenience import reconstruct_stack

[22]:

res = reconstruct_stack(stack=data_obj, stack_ref=data_ref, holoparams=None)

100%|████████████████████████████████████████████████████████████████████████████| 30/30 [00:00<00:00, 42.83it/s]
100%|███████████████████████████████████████████████████████████████████████████| 30/30 [00:00<00:00, 109.35it/s]

[23]:

phase = phase_unwrap(np.angle(res.complex_wave))

[24]:

res.unwrapped_phase = phase

res contains complex_wave, unwrapped_phase and brightfield. It also contains metadata which include the effective pixel size after hologram reconstruction.

[25]:

res.metadata.keys()

[25]:

dict_keys(['stack_shape', 'exposure_time', 'effective_pixelsize', 'acquisition_timestamp', 'drifts_x', 'drifts_y'])

Additionnal metadata can be added as desired.

[26]:

res.metadata['temperature'] = 300
res.metadata['sample'] = 'Magnetite nanocube'

[27]:

res.metadata

[27]:

{'stack_shape': [30, 3838, 3710],
 'exposure_time': 60.0,
 'effective_pixelsize': 1.1636327403696693e-09,
 'acquisition_timestamp': '2025-12-02T14:44:32.458073+00:00',
 'drifts_x': [np.float64(0.0),
  np.float64(0.0),
  np.float64(0.19999999999998863),
  np.float64(0.0),
  np.float64(0.0),
  np.float64(0.19999999999998863),
  np.float64(-0.4000000000000057),
  np.float64(-0.4000000000000057),
  np.float64(-0.4000000000000057),
  np.float64(-1.4000000000000057),
  np.float64(-1.5999999999999943),
  np.float64(-1.0),
  np.float64(-1.4000000000000057),
  np.float64(-1.4000000000000057),
  np.float64(-1.4000000000000057),
  np.float64(-1.0),
  np.float64(-1.5999999999999943),
  np.float64(-1.5999999999999943),
  np.float64(-1.5999999999999943),
  np.float64(-4.599999999999994),
  np.float64(-4.599999999999994),
  np.float64(-4.400000000000006),
  np.float64(-4.199999999999989),
  np.float64(-3.4000000000000057),
  np.float64(-3.5999999999999943),
  np.float64(-3.0),
  np.float64(-3.1999999999999886),
  np.float64(-4.199999999999989),
  np.float64(-3.8000000000000114),
  np.float64(-3.5999999999999943)],
 'drifts_y': [np.float64(0.0),
  np.float64(0.0),
  np.float64(-0.19999999999998863),
  np.float64(-0.19999999999998863),
  np.float64(-0.19999999999998863),
  np.float64(-0.19999999999998863),
  np.float64(0.0),
  np.float64(0.0),
  np.float64(-0.19999999999998863),
  np.float64(0.5999999999999943),
  np.float64(0.5999999999999943),
  np.float64(0.0),
  np.float64(0.0),
  np.float64(0.19999999999998863),
  np.float64(0.5999999999999943),
  np.float64(0.0),
  np.float64(0.4000000000000057),
  np.float64(0.0),
  np.float64(0.0),
  np.float64(2.4000000000000057),
  np.float64(2.1999999999999886),
  np.float64(1.5999999999999943),
  np.float64(1.5999999999999943),
  np.float64(0.5999999999999943),
  np.float64(0.8000000000000114),
  np.float64(0.4000000000000057),
  np.float64(0.5999999999999943),
  np.float64(2.4000000000000057),
  np.float64(2.1999999999999886),
  np.float64(2.4000000000000057)],
 'temperature': 300,
 'sample': 'Magnetite nanocube'}

Phase ramp

Sometimes if the sideband is not perfectly centered or there are some distorsions from the microscope, a phase ramp can be observed across the phase image.

A function to remove this ramp can be used:

[28]:

from libertem_holo.base.utils import remove_phase_ramp

A region-of-interest where no phase gradient is expected can be added as input.

[29]:

def make_rect_fig(ax, rois):
    for roi in rois:
        ax.add_patch(plt.Rectangle(
            (roi[1].start, roi[0].start), roi[1].stop-roi[1].start, roi[0].stop-roi[0].start,
            edgecolor='red', facecolor='none', lw=2,
        ))

[30]:

fig, ax = plt.subplots()
roi = np.s_[50:150, 250:400]
ax.imshow(res.unwrapped_phase)
ax.set_title('Selected ROI for phaseramp removal')
make_rect_fig(ax, rois=[roi])

../_images/usage_reconstruction_50_0.png

[31]:

phase_t, _, _ = remove_phase_ramp(res.unwrapped_phase, roi=roi, method='gradient')
phase_t -= np.mean(phase_t[roi])

[32]:

plt.figure()
plt.plot(res.unwrapped_phase[200, ...], label='Original phase')
plt.plot(phase_t[200, ...], label='After ramp removal')
plt.title('Profile of phase before and after ramp removal')
plt.xlabel('Distance (px)')
plt.ylabel('Phase (rad)')
plt.legend()

[32]:

<matplotlib.legend.Legend at 0x7f5e961941d0>

../_images/usage_reconstruction_52_1.png

The res object can be modified to accept the new phase easily.

[33]:

res.unwrapped_phase = phase_t

Finally, the dataset can be saved as a .npz file.

[34]:

if False:  # enable here if you want to save
    path = ""
    filename = ""
    res.save(path + filename + ".npz")

Phase unwrapping

Phase unwrapping can introduce a lot of artefacts due to low signal-to-noise ratio, hard edges, large phase gradient… The basic phase_unwrap function is directly coming from the skimage package and can be used on a complex wave or a wrapped phase. Additionnal algorithms are also implement in the unwrap module based on quality guided unwrapping and least-square fitting unwrapping. However, the latter still introduces artefacts and should not be used for quantitative analysis of off-axis electron holography.

Quality based unwrapping

[35]:

from libertem_holo.base.unwrap import quality_unwrap, derivative_variance

[36]:

phase_wrapped = np.angle(res.complex_wave)
quality_map = derivative_variance(phase_wrapped)
phase_unwrapped = quality_unwrap(phase_wrapped, quality=quality_map)

# remove phase ramp as above:
phase_unwrapped, _, _ = remove_phase_ramp(phase_unwrapped, roi=roi, method='gradient')
phase_unwrapped -= np.mean(phase_unwrapped[roi])

[37]:

plt.figure()
plt.imshow(phase_unwrapped)
plt.title('Unwrapped phase using quality-guided unwrapping')

[37]:

Text(0.5, 1.0, 'Unwrapped phase using quality-guided unwrapping')

../_images/usage_reconstruction_62_1.png

Laplacian-based unwrapping

The laplacian unwrapping method is a lot faster than the quaity-based unwrap and the unwrapping function included with skimage, but currently only gives a qualitative result. It can be used for a fast visualization, or as a starting point for a quantitative unwrapping.

[38]:

from libertem_holo.base.unwrap import unwrap_phase_laplacian

[39]:

phase_unwrapped_lap = unwrap_phase_laplacian(wrapped_phase=cp.asarray(phase_wrapped), xp=cp).get()

# remove phase ramp as above:
phase_unwrapped_lap, _, _ = remove_phase_ramp(phase_unwrapped_lap, roi=roi, method='gradient')
phase_unwrapped_lap -= np.mean(phase_unwrapped_lap[roi])

[40]:

plt.figure()
plt.imshow(phase_unwrapped_lap)
plt.title('Unrwapped phase using laplacian-based unwrapping')

[40]:

Text(0.5, 1.0, 'Unrwapped phase using laplacian-based unwrapping')

../_images/usage_reconstruction_66_1.png

[41]:

fig, ax = plt.subplots(ncols=2)
ax[0].imshow(np.cos(phase_unwrapped_lap * 8))
ax[0].set_title('Laplacian-based unwrapping')
ax[0].axis("off")
ax[1].imshow(np.cos(phase_unwrapped * 8))
ax[1].set_title('Quality-guided unwrapping')
ax[1].axis("off")
fig.suptitle('Comparison of both unwrapping techniques—\nnotice the artifacts on the edge of the image')

[41]:

Text(0.5, 0.98, 'Comparison of both unwrapping techniques—\nnotice the artifacts on the edge of the image')

../_images/usage_reconstruction_67_1.png

Load and visualise a npz reconstructed dataset

[42]:

if testdata_base_path is not None:
    result_path = testdata_base_path / "dm/holo/align/minus.npz"
else:
    result_path = Path(".") / "minus.npz"

[43]:

data = Results.load(result_path)

[44]:

wave = data.complex_wave
amp = np.abs(wave)
phase = data.unwrapped_phase
bf = data.brightfield
pixel_size = data.metadata['effective_pixelsize'] # in meters !

[45]:

extent = (
    0, phase.shape[1] * pixel_size*1e9,
    0, phase.shape[0] * pixel_size*1e9
)

[46]:

fig, ax = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(13, 5))
idx = 50
roi=np.s_[idx:-idx, idx:-idx]

im = ax[0].imshow(amp[roi], extent=extent, vmin=np.min(clipped(amp)), vmax=np.max(clipped(amp)))
ax[0].set_title('Amplitude')
plt.colorbar(im, fraction=0.046, label='a.u.')
im = ax[1].imshow(phase[roi], extent=extent)
ax[1].set_title('Phase')
plt.colorbar(im, fraction=0.046, label='rad.')
im = ax[2].imshow(bf[roi], extent=extent)
ax[2].set_title('Brightfield')
plt.colorbar(im, fraction=0.046, label='a.u.')

[46]:

<matplotlib.colorbar.Colorbar at 0x7f5e28112180>

../_images/usage_reconstruction_73_1.png

[ ]: