Hyperspectral CT

Hyperspectral CT#

The hsnt module processes hyperspectral neutron transmission data: the transmission of every detector pixel in each of many wavelength bins. A sample made of a few materials has a low-rank attenuation X = W H, with a nonnegative spectrum per component in H and a map per component in W. The module estimates W and H from the data (dehydration), multiplies them back into denoised data (rehydration), estimates the number of components, and reads and writes the hsnt HDF5 layout. The scripts in experiments/hsnt are worked examples, and Command line runs the same steps on files.

Dehydration and rehydration#

Dehydration fits W and H by maximum likelihood: the Poisson likelihood of the counts, written in terms of the transmission T = exp(-X). The natural input is therefore the transmission, the counts divided by the open beam, and a zero count is simply a zero transmission:

transmission = counts / open_beam                       # normalize by the detector calibration
maps, spectra, dataset_type = hsnt.dehydrate(transmission, dataset_type='transmission')
denoised = hsnt.rehydrate([maps, spectra, dataset_type])

The dehydrated form [maps, spectra, dataset_type] is small: the maps have one entry per pixel and component, and the spectra one per component and bin. rehydrate returns the same quantity as the input, transmission or attenuation. When the number of components is not given, it is estimated by likelihood-ratio tests. The components span the materials’ spectra but need not be the materials themselves. Data too large for the device are fitted by chunks of pixels. Given a subspace_basis, dehydrate fits only the maps, so a large scan can be dehydrated piece by piece against spectra fitted on part of it.

mbirtorch.hsnt.dehydrate(data, dataset_type='attenuation', num_materials=None, *, subspace_basis=None, spectra='mle', dose=None, penalty='auto', max_steps=1000, rel_tol=1e-08, max_rank=6, device=None, compile_mode='auto', mode='auto', chunk_pixels=None, max_passes=5, verbose=1, **kwargs)[source]#

Dehydrate a hyperspectral dataset: factor its attenuation as X = W H, W, H >= 0, by maximum likelihood.

The spectral axis is last; the leading axes are kept. The fit maximizes the Poisson likelihood of the counts, written in terms of the transmission T = exp(-X), so zero counts are handled directly. When num_materials is not given, the rank is estimated by likelihood-ratio tests (see estimate_rank). The components are a nonnegative basis of the data, not necessarily the pure materials. Data too large for the device are fitted by chunks of pixels. Given a subspace_basis, only the maps W are fitted.

Parameters:

  • data (numpy.ndarray or torch.Tensor) – Hyperspectral data, spectral axis of length \(N_k\) last. NaN or infinite transmissions count as zero counts; negative transmissions are clipped at zero.

  • dataset_type (str, optional) – ‘attenuation’ or ‘transmission’, with attenuation = -log(transmission). Defaults to ‘attenuation’.

  • num_materials (int, optional) – Rank \(N_m\) of the factorization. Defaults to None, which estimates it.

  • subspace_basis (numpy.ndarray or torch.Tensor, optional) – Spectra of shape \((N_m, N_k)\) to hold fixed; only the maps are then fitted. Defaults to None.

  • spectra (str, optional) – ‘mle’, the maximum-likelihood spectra; ‘unconstrained’, a re-estimate that removes the bias of the nonnegativity constraint at low dose; ‘support’, which selects the components present in each pixel and refits (needs dose). Defaults to ‘mle’.

  • dose (float, optional) – Open-beam counts per pixel and bin, for spectra=’support’. Defaults to None.

  • penalty (str or float, optional) – Charge per component for spectra=’support’, as a multiple of log(N_k), or ‘auto’. Defaults to ‘auto’.

  • max_steps (int, optional) – Solver iteration cap. Defaults to 1000.

  • rel_tol (float, optional) – Stop when the relative loss change stays below this. Defaults to 1e-8.

  • max_rank (int, optional) – Largest rank the estimate considers. Defaults to 6.

  • device (str, optional) – Torch device. Defaults to None: CUDA if available, else CPU.

  • compile_mode (str, optional) – ‘auto’, ‘on’ or ‘off’ for torch.compile. Defaults to ‘auto’.

  • mode (str, optional) – ‘full’ fits on the device at once, ‘stream’ by chunks of pixels, ‘auto’ picks from the available memory. Defaults to ‘auto’.

  • chunk_pixels (int, optional) – Pixels per chunk when streamed. Defaults to None, from the available memory.

  • max_passes (int, optional) – Passes over the chunks when streamed. Defaults to 5.

  • verbose (int, optional) – 0 prints nothing; 1 a summary; 2 also the rank search. Defaults to 1.

Returns:

list – [subspace_data, subspace_basis, dataset_type]: W reshaped to the leading axes plus \(N_m\), H of shape \((N_m, N_k)\), both float32, and the input’s dataset_type, so that rehydrate returns the same quantity as the input.

mbirtorch.hsnt.rehydrate(dehydrated_data, hyperspectral_idx=None)[source]#

Rehydrate dehydrated hyperspectral data: multiply the maps by the spectra, for all or selected bins.

The method is described in M. S. N. Chowdhury et al., “Fast Hyperspectral Neutron Tomography,” IEEE Trans. Computational Imaging, vol. 11, pp. 663-677, 2025, doi:10.1109/TCI.2025.3567854.

Parameters:

  • dehydrated_data – [subspace_data, subspace_basis, dataset_type]: the maps, with the component axis of length \(N_s\) last; the spectra, shape \((N_s, N_k)\); and ‘attenuation’ or ‘transmission’.

  • hyperspectral_idx – Indices of the spectral bins to rehydrate. Defaults to None, all \(N_k\) bins.

Returns:

numpy.ndarray – Hyperspectral data with the shape of subspace_data, except that the last axis holds the selected bins; attenuation or transmission as dataset_type says.

mbirtorch.hsnt.hyper_denoise(data, dataset_type='attenuation', num_materials=None, **kwargs)[source]#

Denoise a hyperspectral dataset: dehydrate, then rehydrate.

Parameters:

  • data (numpy.ndarray or torch.Tensor) – Hyperspectral data with the spectral axis last.

  • dataset_type (str, optional) – ‘attenuation’ or ‘transmission’. Defaults to ‘attenuation’.

  • num_materials (int, optional) – Rank of the factorization. Defaults to None, which estimates it.

  • **kwargs – The keyword arguments of dehydrate().

Returns:

numpy.ndarray – The denoised data, float32, with the shape and dataset_type of the input.

mbirtorch.hsnt.estimate_rank(data, dataset_type='attenuation', max_rank=6, device=None, pool='auto', verbose=0)[source]#

Estimate the number of components of a hyperspectral dataset by likelihood-ratio tests.

Ranks 1 to max_rank are fitted in turn, and a component is kept while its gain in log-likelihood stands well above the gain a component that fits only noise would give. When the leading axes are (views, rows, cols) or (rows, cols), the test is also run on blocks of pooled pixels, which is more sensitive to faint materials at low dose, and the larger of the two ranks is returned. This is the rank dehydrate uses when num_materials is not given.

Parameters:

  • data (numpy.ndarray or torch.Tensor) – Hyperspectral data with any leading axes and the spectral axis last.

  • dataset_type (str, optional) – ‘attenuation’ or ‘transmission’. Defaults to ‘attenuation’.

  • max_rank (int, optional) – Largest rank considered. Defaults to 6.

  • device (str, optional) – Torch device. Defaults to None: CUDA if available, else CPU.

  • pool (str or int, optional) – Pooling block size in pixels; ‘auto’ chooses it from the dose, 0 disables pooling. Defaults to ‘auto’.

  • verbose (int, optional) – 1 prints each search’s gains and decision. Defaults to 0.

Returns:

(rank, note, detail) – the rank, a one-line account of how it was chosen, and a dict with the searches’ gains, effective doses and thresholds.

Import/Export#

mbirtorch.hsnt.import_hsnt_data_hdf5(filename)[source]#

Import a hyperspectral dataset and its metadata from an HDF5 file.

Parameters:

filename – Path to the HDF5 file.

Returns:

[data, metadata] – data is an ndarray with the spectral axis last (hyperspectral form), a list [subspace_data, subspace_basis, dataset_type] (dehydrated form), or None; metadata is a dict with the keys of create_hsnt_metadata.

mbirtorch.hsnt.create_hsnt_metadata(**kwargs)[source]#

Create the metadata dict of a hyperspectral neutron dataset. Keys not given are None.

Parameters:

  • dataset_name – Character string with the name of the dataset.

  • dataset_type – ‘attenuation’ or ‘transmission’.

  • dataset_modality – ‘hyperspectral neutron’.

  • wavelengths – Array of wavelength values in Angstroms.

  • alu_unit – Character string defining geometry unit (e.g., ‘mm’ or ‘cm’).

  • alu_value – Float that represents the value of 1 ALU in the defined unit.

  • delta_det_channel – Detector channel spacing in ALU.

  • delta_det_row – Detector row spacing in ALU.

  • dataset_geometry – ‘parallel’ or ‘cone’.

  • angles – Array of view angles in degrees.

  • det_channel_offset – Assumed offset between center of rotation and center of detector in ALU.

  • source_detector_dist – Distance from source to detector in ALU.

  • source_iso_dist – Distance from source to iso in ALU.

Returns:

dict – The metadata.

Example

>>> metadata = create_hsnt_metadata(dataset_name="sample1", dataset_type="transmission",
...                                 wavelengths=np.linspace(1.0, 5.0, 50), angles=np.linspace(0, 180, 10))
mbirtorch.hsnt.export_hsnt_data_hdf5(filename, data, metadata=None)[source]#

Export a hyperspectral dataset and its metadata to an HDF5 file.

Parameters:

  • filename – Path to the HDF5 file.

  • data – ndarray with the spectral axis last (hyperspectral form), or [subspace_data, subspace_basis, dataset_type] (dehydrated form).

  • metadata – The dict from create_hsnt_metadata. Defaults to None.

Synthetic data#

A simulated scan for trying the method: a phantom of three materials, three synthetic spectra, and the hyperspectral sinogram of the phantom through a tomography model, with Poisson counts. demo_7_hyperspectral.py dehydrates, reconstructs, and rehydrates one.

mbirtorch.hsnt.gen_material_phantom(recon_shape, density=0.08)[source]#

A phantom of three materials: three cylinders along the slice axis, spaced evenly around the center of the volume, one of each material.

Parameters:

  • recon_shape (tuple of int) – (num_rows, num_cols, num_slices).

  • density (float, optional) – The value inside each cylinder, which scales that material’s spectrum to the attenuation per voxel. Defaults to 0.08.

Returns:

numpy.ndarray – One map per material, float32, shape (3,) + recon_shape, density inside that material’s cylinder and 0 elsewhere.

mbirtorch.hsnt.synthetic_material_basis(num_wavelengths=200)[source]#

Three synthetic attenuation spectra, one per material, for simulated scans.

Each spectrum falls slowly with wavelength and drops sharply at one wavelength of its own, the shape of a Bragg edge. The values are of order one and are not those of any real material; a material map gives the attenuation per unit length by scaling them.

Parameters:

num_wavelengths (int, optional) – The number of wavelength bins. Defaults to 200.

Returns:

numpy.ndarray – The spectra, float32, shape (3, num_wavelengths), all values positive.

mbirtorch.hsnt.generate_hyper_sinogram(ct_model, material_maps, material_basis, dosage_rate=300, noisy=True, seed=None)[source]#

Simulate a hyperspectral scan of a volume of several materials.

Each material map is projected with the model, and the projections are combined with the spectra into the attenuation at every wavelength. The open beam is dosage_rate counts in every pixel and bin, and the counts through the sample are Poisson with that mean times the transmission. Divide the counts by the open beam to get the transmission that dehydrate() takes with dataset_type='transmission'.

Parameters:

  • ct_model (TomographyModel) – The model of the scan.

  • material_maps (numpy.ndarray) – One volume per material, shape (num_materials,) + recon_shape, the attenuation per unit length at a spectrum value of 1.

  • material_basis (numpy.ndarray) – The spectra, shape (num_materials, num_wavelengths).

  • dosage_rate (float, optional) – Open beam counts per pixel and bin. Defaults to 300.

  • noisy (bool, optional) – If True, the counts are Poisson; if False, their mean. Defaults to True.

  • seed (int, optional) – Seed of the Poisson draw. Defaults to None.

Returns:

tuple – (counts, open_beam, attenuation), each of shape (num_views, num_det_rows, num_det_channels, num_wavelengths), float32: the counts through the sample, the open beam, and the noiseless attenuation.

mbirtorch.hsnt.generate_hyper_data(material_basis, num_angles=1, detector_rows=64, detector_columns=64, dosage_rate=300, material_density=None, verbose=1, *, noisy=True)[source]#

Simulate hyperspectral neutron counts for \(N_m=3\) materials (Ni, Cu, Al) and \(N_k\) wavelength bins.

The open beam is noiseless, dosage_rate counts in every pixel and bin; the sample counts are Poisson. Divide the counts by the open beam to get the transmission that dehydrate takes with dataset_type=’transmission’; a zero count is a zero transmission, which the likelihood handles.

Parameters:

  • material_basis – ndarray of shape \((N_m, N_k)\), where rows are material linear attenuation coefficient spectra.

  • num_angles – Number of view angles \((N_v)\). Defaults to 1.

  • detector_rows – Number of rows in the detector \((N_r)\). Defaults to 64.

  • detector_columns – Number of columns in the detector \((N_c)\). Defaults to 64.

  • dosage_rate – Neutron dosage rate during hyperspectral data collection. Defaults to 300.

  • material_density – Material density (vol. fraction) for Ni, Cu, and Al. Defaults to {“Ni”: 0.2, “Cu”: 0.2, “Al”: 1.0}.

  • verbose – Verbosity level. If 0, prints nothing; if 1, prints details; if >1, also generates plots. Defaults to 1.

  • noisy – Whether to generate noisy data (keyword only). Defaults to True.

Returns:

A list in the form [counts, open_beam, angles, gt_hyper_projection].

  • counts: Poisson counts through the sample, shape \((N_v, N_r, N_c, N_k)\), in the dtype of material_basis.

  • open_beam: Open-beam counts of the same shape, dosage_rate everywhere.

  • angles: ndarray of view angles in radians.

  • gt_hyper_projection: Ground truth noiseless attenuation of the same shape.

Command line#

mbirtorch-hsnt runs the same steps from the shell, on an HDF5 file in the hsnt layout or on a directory of TIFF images, one per wavelength bin. Its subcommands are inspect, convert (TIFFs to HDF5, normalized by an open beam), dehydrate, rehydrate and denoise:

mbirtorch-hsnt convert sample_tifs/ --open-beam open_beam/ -o sample.h5
mbirtorch-hsnt denoise sample.h5 -o results/

mbirtorch-hsnt <subcommand> -h lists the options.