create an NXtomo (from scratch)

The goal of this tutorial is to create an NXtomo from scratch. In the sense without converting directly from a bliss scan.

description of the example

Let say we want to create an NXtomo matching the following sequence

frame index	Rotation Angle (in degree)	Frame Type	control Image Key	Image Key
0	0	Dark	2	2
1	0	Flat	1	1
2-201	0 - 89.9	Projection	0	0
202	90	Flat	1	1
203 - 402	90 - 180	Projection	0	0
403	180	Flat	1	1
404	90	Alignment	-1	0

create dummy dataset

In order to siplify the setup we will take only create one dark and one flat frame that we will reuse. But keep in mind that those are raw data. So in ‘real’ life it is expected that you will have several frame for dark and several frames for flats and there will differ of course depending on when you acquire them.

# %pylab inline

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

create some projection

phantom = shepp_logan_phantom()
projections = {}
proj_rotation_angles = numpy.linspace(0., 180., 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 some dark

max_shape = max(phantom.shape)

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

# imshow(dark)

create some flat

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

add some noise to radios

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

# imshow(radios[0])

create some alignment

In order to keep it simple we will pick one of the radio created

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

create an NXtomo that fits the sequence we want

import nxtomo

from nxtomo import NXtomo

my_nxtomo = NXtomo()

provide mandatory data for Contrast Tomography

Mandatory information for contrast tomography are: * detector frames: raw data * image-key (control): for each frame the “type” of frame (projections, flats, darks and alignment). * rotation angles: for each frame the rotation angle in degree

detector frames

to fit the sequence describe previously we need to create the following sequence: dark, flat, first half of the projections, flat, second half of the projections, flat and alignment frame.

And we need to provide them as a numpy array (3d)

# 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 angle

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

Field of view

field of view can either be 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 = 1e-7  # pixel size must be provided in SI: meter

but for attribute with a unit you can specify the unit the value should be “converted to” using the ‘unit’ attribute like:

my_nxtomo.instrument.detector.x_pixel_size = my_nxtomo.instrument.detector.y_pixel_size = 0.1
my_nxtomo.instrument.detector.x_pixel_size.unit = my_nxtomo.instrument.detector.y_pixel_size.unit = "micrometer"

When the unit is provided it will be stored as a property of the dataset. It must be interpreted by the software reading the NXtomo.

save the nx 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 data saved

We can use some validator from tomoscan to insure we have enought data to be treated by nabu

try:
    import tomoscan
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 check the layout of the file to insure it seems valid as well

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: float64
        • y_pixel_size [📋]: scalar entries, dtype: float64
      • 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 pratice is also to check frames, image_key and rotation angles to insure values seems valid.

# ! silx view output/nxtomo.nx

reconstruct using nabu

now that we have a valid nxtomo we are able to reconstruct it using nabu.

We create a nabu configuration file for contrast tomography reconstruction name nabu-ct.conf in order to reconstruct one slice of the volume.

on the configuration you must disable take_logarithm due to the dataset.

if nabu is installed you can run it:

# ! nabu nabu-cf.conf

provide mandatory data for Phase Tomography

in order to compute Phase Tomography you must also register: * incoming beam energy (in keV) * sample / detector distance (in meter)

we can take back the existing my_nxtomo and add it the missing elements

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

And then you can reconstruct it with phase retrieval from modifing the nabu configuration file.

provide more metadata

you can also provide x, y and z translation of the sample during the acquisition.

my_nxtomo.sample.x_translation = [0, 12]

as a sample name, source information, start and end time

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,
)