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.
[1]:
# %pylab inline
[2]:
from skimage.data import shepp_logan_phantom
from skimage.transform import radon
import numpy
Create projection frames#
[3]:
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)
[4]:
# imshow(phantom)
[5]:
radios = numpy.swapaxes(sinograms, 2, 0)
radios = numpy.swapaxes(radios, 1, 2)
[6]:
# imshow(radios[0])
Create dark frames#
[7]:
max_shape = max(phantom.shape)
[8]:
dark = numpy.zeros((20, max_shape))
[9]:
# imshow(dark)
Create flat frames#
[10]:
flat = numpy.ones((20, max_shape))
Add noise to radiographs#
[11]:
tmp_radios = []
for radio in radios:
tmp_radios.append(dark + radio * (flat - dark))
radios = numpy.asarray(tmp_radios)
[12]:
# imshow(radios[0])
Create alignment frames#
To keep things simple we reuse one of the previously generated radiographs.
[13]:
alignment = radios[200]
alignment_angle = proj_rotation_angles[200]
Build an NXtomo that matches the sequence#
[ ]:
[14]:
from nxtomo import NXtomo
[15]:
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.
[16]:
# 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#
[17]:
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#
[18]:
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.
[19]:
my_nxtomo.instrument.detector.field_of_view = "Full"
Pixel size#
[20]:
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#
[21]:
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.
[22]:
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.
[23]:
from h5glance import H5Glance
H5Glance(nx_tomo_file_path)
[23]:
A good practice is to verify the frames, image keys, and rotation angles to ensure the values are consistent.
[24]:
# ! 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:
[25]:
# ! 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.
[26]:
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.
[27]:
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.
[28]:
my_nxtomo.sample.name = "my sample"
[29]:
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)
[30]:
my_nxtomo.save(
file_path=nx_tomo_file_path,
data_path="entry",
overwrite=True,
)