Create an NXtomo (from scratch)

The goal of this tutorial is to build an NXtomo file from scratch, without converting it directly from a BLISS scan.

Example description

Let’s create an NXtomo that matches the following sequence:

Frame index	Rotation angle (degrees)	Frame type	Control image key	Image key
0	0.0	Dark	2	2
1	0.0	Flat	1	1
2-201	0.0 - 89.9	Projection	0	0
202	90.0	Flat	1	1
203-402	90.0 - 180.0	Projection	0	0
403	180.0	Flat	1	1
404	90.0	Alignment	-1	0

create dummy dataset

To simplify the setup we create a single dark frame and a single flat frame that we reuse across the tutorial. Remember that these are raw frames. In real acquisitions you usually record several dark and flat frames, and they vary depending on when they are acquired.

# %pylab inline

from skimage.data import shepp_logan_phantom
from skimage.transform import radon
import numpy

Create projection frames

phantom = shepp_logan_phantom()
projections = {}
proj_rotation_angles = numpy.linspace(0.0, 180.0, max(phantom.shape), endpoint=False)
sinogram = radon(phantom, theta=proj_rotation_angles)

sinograms = numpy.asarray([sinogram] * 20)

# imshow(phantom)

radios = numpy.swapaxes(sinograms, 2, 0)
radios = numpy.swapaxes(radios, 1, 2)

# imshow(radios[0])

Create dark frames

max_shape = max(phantom.shape)

dark = numpy.zeros((20, max_shape))

# imshow(dark)

Create flat frames

flat = numpy.ones((20, max_shape))

Add noise to radiographs

tmp_radios = []
for radio in radios:
    tmp_radios.append(dark + radio * (flat - dark))
radios = numpy.asarray(tmp_radios)

# imshow(radios[0])

Create alignment frames

To keep things simple we reuse one of the previously generated radiographs.

alignment = radios[200]
alignment_angle = proj_rotation_angles[200]

Build an NXtomo that matches the sequence

from nxtomo import NXtomo

my_nxtomo = NXtomo()

Provide mandatory data for contrast tomography

Mandatory information for contrast tomography includes:

• detector frames: raw data

• image-key (control): the frame type for each image (projections, flats, darks, alignment)

• rotation angles: the rotation angle in degrees for every frame

Detector frames

To follow the sequence described earlier we create the series: dark, flat, first half of the projections, flat, second half of the projections, flat, and finally the alignment frame.

All frames must be provided as a three-dimensional NumPy array.

# reshape dark, flat and alignment that need to be 3d when numpy.concatenate is called
darks_stack = dark.reshape(1, dark.shape[0], dark.shape[1])
flats_stack = flat.reshape(1, flat.shape[0], flat.shape[1])
alignment_stack = alignment.reshape(1, alignment.shape[0], alignment.shape[1])

assert darks_stack.ndim == 3
assert flats_stack.ndim == 3
assert alignment_stack.ndim == 3
assert radios.ndim == 3
print("radios shape is", radios.shape)
# create the array
data = numpy.concatenate(
    [
        darks_stack,
        flats_stack,
        radios[:200],
        flats_stack,
        radios[200:],
        flats_stack,
        alignment_stack,
    ]
)
assert data.ndim == 3
print(data.shape)
# then register the data to the detector
my_nxtomo.instrument.detector.data = data

radios shape is (400, 20, 400)
(405, 20, 400)

Image-key control

from nxtomo.nxobject.nxdetector import ImageKey

image_key_control = numpy.concatenate(
    [
        [ImageKey.DARK_FIELD] * 1,
        [ImageKey.FLAT_FIELD] * 1,
        [ImageKey.PROJECTION] * 200,
        [ImageKey.FLAT_FIELD] * 1,
        [ImageKey.PROJECTION] * 200,
        [ImageKey.FLAT_FIELD] * 1,
        [ImageKey.ALIGNMENT] * 1,
    ]
)

# insure with have the same number of frames and image key
assert len(image_key_control) == len(data)
# print position of flats in the sequence
print("flats indexes are", numpy.where(image_key_control == ImageKey.FLAT_FIELD))
# then register the image keys to the detector
my_nxtomo.instrument.detector.image_key_control = image_key_control

flats indexes are (array([  1, 202, 403]),)

Rotation angles

import pint

ureg = pint.get_application_registry()

rotation_angle = numpy.concatenate(
    [
        [
            0.0,
        ],
        [
            0.0,
        ],
        proj_rotation_angles[:200],
        [
            90.0,
        ],
        proj_rotation_angles[200:],
        [
            180.0,
        ],
        [
            90.0,
        ],
    ]
)
assert len(rotation_angle) == len(data)
# register rotation angle to the sample
my_nxtomo.sample.rotation_angle = rotation_angle * ureg.degree

Field of view

The field of view can be either Half or Full.

my_nxtomo.instrument.detector.field_of_view = "Full"

Pixel size

my_nxtomo.instrument.detector.x_pixel_size = (
    my_nxtomo.instrument.detector.y_pixel_size
) = (10 * ureg.micrometer)

When a unit is provided it is stored as a dataset attribute. The software that reads the NXtomo file must interpret that unit.

Save the NXtomo to disk

import os

os.makedirs("output", exist_ok=True)
nx_tomo_file_path = os.path.join("output", "nxtomo.nx")
my_nxtomo.save(file_path=nx_tomo_file_path, data_path="entry", overwrite=True)

Check the saved data

Use the tomoscan validator to ensure we have enough information for nabu to process the dataset.

try:
    import tomoscan #NOQA
except ImportError:
    has_tomoscan = False
else:
    from tomoscan.esrf import NXtomoScan
    from tomoscan.validator import ReconstructionValidator

    has_tomoscan = True

if has_tomoscan:
    scan = NXtomoScan(nx_tomo_file_path, entry="entry")
    validator = ReconstructionValidator(
        scan, check_phase_retrieval=False, check_values=True
    )
    assert validator.is_valid()

You can also inspect the file layout to confirm it looks correct.

from h5glance import H5Glance

H5Glance(nx_tomo_file_path)

• output/nxtomo.nx
  • entry
    • data
      • data → /entry/instrument/detector/data
      • image_key → /entry/instrument/detector/image_key
      • image_key_control → /entry/instrument/detector/image_key_control
      • rotation_angle → /entry/sample/rotation_angle
    • instrument
      • detector
        • data [📋]: 405 × 20 × 400 entries, dtype: float64
        • field_of_view [📋]: scalar entries, dtype: UTF-8 string
        • image_key [📋]: 405 entries, dtype: int64
        • image_key_control [📋]: 405 entries, dtype: int64
        • x_pixel_size [📋]: scalar entries, dtype: int64
        • y_pixel_size [📋]: scalar entries, dtype: int64
      • source
        • name [📋]: scalar entries, dtype: UTF-8 string
        • probe [📋]: scalar entries, dtype: UTF-8 string
        • type [📋]: scalar entries, dtype: UTF-8 string
    • sample
      • rotation_angle [📋]: 405 entries, dtype: float64
    • definition [📋]: scalar entries, dtype: UTF-8 string

A good practice is to verify the frames, image keys, and rotation angles to ensure the values are consistent.

# ! silx view output/nxtomo.nx

Reconstruct using nabu

Now that we have a valid NXtomo we can reconstruct it with nabu.

Create a nabu configuration file for contrast tomography, named nabu-ct.conf, to reconstruct one slice of the volume.

In the configuration file make sure to disable take_logarithm, because the dataset already contains processed values.

If nabu is installed you can run it:

# ! nabu nabu-cf.conf

Provide mandatory data for phase tomography

To perform phase tomography you must also record:

• incoming beam energy (in keV)

• sample-to-detector distance (in metres)

We can reuse the existing my_nxtomo and add the missing elements.

my_nxtomo.energy = 12.5 * ureg.keV  # in keV by default
my_nxtomo.instrument.detector.distance = 0.2 * ureg.meter

Then you can reconstruct it with phase retrieval by modifying the nabu configuration file accordingly.

Provide more metadata

You can also store the sample translations along x, y, and z during the acquisition.

my_nxtomo.sample.x_translation = [0, 12] * ureg.millimeter

Add useful contextual information such as the sample name, source details, and acquisition start and end times.

my_nxtomo.sample.name = "my sample"

from datetime import datetime

my_nxtomo.instrument.source.name = "ESRF"  # default value
my_nxtomo.instrument.source.type = "Synchrotron X-ray Source"  # default value
my_nxtomo.start_time = datetime.now()
my_nxtomo.end_time = datetime(2022, 2, 27)

my_nxtomo.save(
    file_path=nx_tomo_file_path,
    data_path="entry",
    overwrite=True,
)