Creating an NXtomo file

Creating an NXtomo file#

HTTomo can automatically find tomography data in a NeXus NXtomo file. This tutorial shows how to pack either NumPy arrays or stacks of TIFF images into that format.

Download the complete create_nxtomo.py script before following the examples. The script is self-contained and may also be imported as a Python module.

Install the packages used by the script in the environment where the packing will run:

python -m pip install numpy h5py tifffile

tifffile is only needed for the TIFF example.

Input layout#

Projection, flat-field, and dark-field data must have the axis order (frames, detector_y, detector_x). For projections, frames is the rotation-angle axis. The angle array must be one-dimensional, measured in degrees, and contain exactly one value per projection.

A single flat or dark image may be supplied as a two-dimensional array. The script adds its frame axis automatically. Flats and darks are optional, but operations that perform flat/dark correction need meaningful calibration images.

From NumPy arrays#

Import write_nxtomo() from the downloaded script. Arrays can come from numpy.load, another Python library, or calculations performed in the same program:

import numpy as np

from create_nxtomo import write_nxtomo

projections = np.load("projections.npy")  # (n_angles, detector_y, detector_x)
angles = np.load("angles.npy")             # (n_angles,), in degrees
flats = np.load("flats.npy")               # (n_flats, detector_y, detector_x)
darks = np.load("darks.npy")               # (n_darks, detector_y, detector_x)

write_nxtomo(
    "scan.nxs",
    projections,
    angles,
    flats=flats,
    darks=darks,
    sample_name="my sample",
    compression="gzip",
)

Use overwrite=True only when an existing output file should be replaced. Set compression=None for faster writing and a larger file, or compression="lzf" for lightweight compression. The default is compression="gzip".

From TIFF stacks#

Keep each image type in its own directory, for example:

scan/
├── projections/
│   ├── projection_0000.tif
│   ├── projection_0001.tif
│   └── ...
├── flats/
│   ├── flat_0000.tif
│   └── ...
├── darks/
│   ├── dark_0000.tif
│   └── ...
└── angles.npy

Then run the script. Quote each glob so that the script, rather than the shell, receives and sorts the complete file list:

python create_nxtomo.py scan.nxs \
    --projections "scan/projections/*.tif" \
    --angles scan/angles.npy \
    --flats "scan/flats/*.tif" \
    --darks "scan/darks/*.tif" \
    --sample-name "my sample"

The angles may instead be stored as a one-column .txt file or a comma-separated .csv file. Omit --flats or --darks when that image type is unavailable. Use --overwrite to replace an existing output file.

The TIFF loader sorts filenames naturally, so projection_2.tif precedes projection_10.tif. The filenames still need to encode the acquisition order correctly. Every matched TIFF must be a single two-dimensional grayscale image, and all images must have the same height and width.

What the script writes#

The output stores darks, then flats, then projections in one three-dimensional detector dataset. The image_key gives every frame its NXtomo meaning:

  • 0 – projection;

  • 1 – flat field; and

  • 2 – dark field.

Calibration frames are assigned the first projection angle because NXtomo requires one rotation angle per frame. HTTomo uses image_key to select only the projection frames and their corresponding angles.

The important part of the resulting file is:

/entry                              NXentry
├── definition                     "NXtomo"
├── instrument                     NXinstrument
│   └── detector                   NXdetector
│       ├── data                   (frames, detector_y, detector_x)
│       └── image_key              (frames,)
├── sample                         NXsample
│   └── rotation_angle             (frames,), units="deg"
└── data                           NXdata
    ├── data                       link to detector/data
    ├── image_key                  link to detector/image_key
    └── rotation_angle             link to sample/rotation_angle

The links in /entry/data do not duplicate the arrays. They expose the standard NXtomo paths that HTTomo’s automatic discovery uses.

Loading the result in HTTomo#

Use the standard tomography loader and set all discoverable paths to auto:

- method: standard_tomo
  module_path: httomo.data.hdf.loaders
  parameters:
    data_path: auto
    image_key_path: auto
    rotation_angles: auto

Pass scan.nxs as the input data file when running HTTomo. If the file has no flats or darks, HTTomo supplies dummy calibration arrays. Alternatively, set flats: ignore or darks: ignore explicitly when the corresponding correction should not use stored calibration images; see Darks and flats.

Reusable NeXus writer code#

The core NumPy writer is included below for reference. The downloadable script also contains TIFF loading, validation, filename sorting, and its command-line interface.

Reusable NeXus writer code
  1def write_nxtomo(
  2    output_file: str | Path,
  3    projections: np.ndarray,
  4    angles: np.ndarray,
  5    *,
  6    flats: np.ndarray | None = None,
  7    darks: np.ndarray | None = None,
  8    sample_name: str = "sample",
  9    compression: str | None = "gzip",
 10    overwrite: bool = False,
 11) -> Path:
 12    """Write NumPy arrays to an HTTomo-compatible NXtomo file.
 13
 14    Parameters
 15    ----------
 16    output_file
 17        Destination ``.nxs`` or ``.h5`` file.
 18    projections
 19        Projection stack with shape ``(angles, detector_y, detector_x)``.
 20    angles
 21        One rotation angle in degrees per projection.
 22    flats, darks
 23        Optional calibration stacks. A single 2D image is also accepted.
 24    sample_name
 25        Descriptive sample name stored in the NXsample group.
 26    compression
 27        HDF5 compression: ``"gzip"``, ``"lzf"``, or ``None``.
 28    overwrite
 29        Replace ``output_file`` if it already exists.
 30    """
 31    projections = _as_frame_stack(projections, "projections")
 32    angles = np.asarray(angles)
 33    if angles.ndim != 1 or len(angles) != len(projections):
 34        raise ValueError("angles must be 1D with one value per projection")
 35    if (
 36        not (
 37            np.issubdtype(angles.dtype, np.integer)
 38            or np.issubdtype(angles.dtype, np.floating)
 39        )
 40        or not np.isfinite(angles).all()
 41    ):
 42        raise ValueError("angles must contain finite numeric values")
 43
 44    optional_stacks = []
 45    for name, stack in (("darks", darks), ("flats", flats)):
 46        if stack is None:
 47            optional_stacks.append(None)
 48            continue
 49        stack = _as_frame_stack(stack, name)
 50        if stack.shape[1:] != projections.shape[1:]:
 51            raise ValueError(
 52                f"{name} frame shape {stack.shape[1:]} does not match "
 53                f"projection frame shape {projections.shape[1:]}"
 54            )
 55        optional_stacks.append(stack)
 56    darks, flats = optional_stacks
 57
 58    compression = None if compression == "none" else compression
 59    if compression not in (None, "gzip", "lzf"):
 60        raise ValueError("compression must be 'gzip', 'lzf', or None")
 61
 62    stacks = [stack for stack in (darks, flats, projections) if stack is not None]
 63    output_dtype = np.result_type(*(stack.dtype for stack in stacks))
 64    dark_count = 0 if darks is None else len(darks)
 65    flat_count = 0 if flats is None else len(flats)
 66    projection_count = len(projections)
 67    frame_count = dark_count + flat_count + projection_count
 68    detector_y, detector_x = projections.shape[1:]
 69
 70    image_key = np.concatenate(
 71        (
 72            np.full(dark_count, 2, dtype=np.int8),
 73            np.full(flat_count, 1, dtype=np.int8),
 74            np.zeros(projection_count, dtype=np.int8),
 75        )
 76    )
 77    # NXtomo requires an angle for every frame. Calibration frames use the
 78    # first projection angle; HTTomo selects projection angles via image_key.
 79    frame_angles = np.concatenate(
 80        (np.full(dark_count + flat_count, angles[0]), angles)
 81    ).astype(np.float32, copy=False)
 82
 83    output_file = Path(output_file)
 84    output_file.parent.mkdir(parents=True, exist_ok=True)
 85    mode = "w" if overwrite else "w-"
 86    now = datetime.now(timezone.utc).isoformat()
 87
 88    with h5py.File(output_file, mode) as nexus_file:
 89        nexus_file.attrs.update(
 90            default="entry", file_name=output_file.name, file_time=now
 91        )
 92
 93        entry = _nx_group(nexus_file, "entry", "NXentry")
 94        entry.attrs["default"] = "data"
 95        entry.create_dataset("definition", data="NXtomo")
 96        entry.create_dataset("title", data=sample_name)
 97        entry.create_dataset("start_time", data=now)
 98
 99        instrument = _nx_group(entry, "instrument", "NXinstrument")
100        detector = _nx_group(instrument, "detector", "NXdetector")
101        data = detector.create_dataset(
102            "data",
103            shape=(frame_count, detector_y, detector_x),
104            dtype=output_dtype,
105            chunks=(1, min(detector_y, 256), min(detector_x, 256)),
106            compression=compression,
107        )
108        data.attrs.update(
109            interpretation="image",
110            axes="frame,detector_y,detector_x",
111        )
112        key = detector.create_dataset("image_key", data=image_key)
113        key.attrs["meaning"] = "0=projection; 1=flat; 2=dark"
114
115        first = 0
116        for stack in (darks, flats, projections):
117            if stack is not None:
118                data[first : first + len(stack)] = stack
119                first += len(stack)
120
121        sample = _nx_group(entry, "sample", "NXsample")
122        sample.create_dataset("name", data=sample_name)
123        rotation_angle = sample.create_dataset("rotation_angle", data=frame_angles)
124        rotation_angle.attrs["units"] = "deg"
125
126        # NXdata contains hard links, not copies. These paths are also the
127        # locations used by HTTomo's automatic NXtomo discovery.
128        nx_data = _nx_group(entry, "data", "NXdata")
129        nx_data.attrs["signal"] = "data"
130        nx_data.attrs["axes"] = np.asarray(
131            ["rotation_angle", ".", "."], dtype=h5py.string_dtype()
132        )
133        nx_data["data"] = data
134        nx_data["image_key"] = key
135        nx_data["rotation_angle"] = rotation_angle
136
137        entry.create_dataset("end_time", data=datetime.now(timezone.utc).isoformat())
138
139    return output_file
140
141