Algorithm background

This page summarises what fitscube-rs does when it combines a set of single-plane FITS images into a cube and when it extracts a plane back out. It is a Rust port of the fitscube Python package and follows the same conventions, so cubes produced by either tool are interchangeable.

Combining images into a cube

The input is a list of FITS images, each holding a single plane: one frequency channel (the default) or one time step (--time-domain). fitscube-rs reads the 2D data from each image and stacks them along a new leading axis, producing a 3D (or 4D, with a degenerate Stokes axis preserved) cube whose plane order follows the order the files are given.

Every input is checked for a consistent shape and pixel grid; a plane that does not match the others is an error rather than a silent reshape. The spatial WCS (CRVAL1/2, CDELT1/2, projection, …) is taken from the first image and carried onto the cube unchanged.

Frequency vs time domain

The new axis can be either spectral or temporal:

  • Frequency (default): the per-plane coordinate is the spectral value of each image, read from its spectral WCS axis (CTYPE/CRVAL/CUNIT) or supplied explicitly. The cube gains a FREQ axis.

  • Time (--time-domain): the per-plane coordinate is a time stamp and the cube gains a time axis instead.

In both cases the per-plane coordinates may be provided out-of-band via a spec file (--spec-file) or an inline list (--specs), which overrides what is read from the headers; --ignore-spec skips reading the coordinate entirely and writes a bare pixel axis.

Even vs uneven spacing detection

After collecting the per-plane coordinates, fitscube-rs decides how to describe the new axis in the output WCS:

  • Evenly spaced: if successive coordinates differ by a constant step (within a small tolerance), the axis is written as a linear WCS — a single reference value (CRVALn) and increment (CDELTn). This is the compact, standard form most downstream tools expect.

  • Unevenly spaced: if the step varies, no single CDELT can describe the axis without losing information. fitscube-rs instead writes the full list of per-plane coordinates as an explicit table alongside the cube, so the exact frequency (or time) of every plane is recoverable. The user is warned that the axis is non-linear.

Per-channel beams (BEAMS table)

Radio images frequently carry a restoring beam (BMAJ, BMIN, BPA) that differs from plane to plane. fitscube-rs preserves this:

  • If all input planes share one beam, the single BMAJ/BMIN/BPA is written to the primary header.

  • If the beams differ across planes, fitscube-rs writes a CASA-style BEAMS binary-table extension — one row per plane with that plane’s beam (and its channel/Stokes index) — and sets CASAMBM=T in the primary header. This is the multi-beam convention understood by CASA and astropy, so per-channel beam information survives the round trip into the cube.

Bounding-box trimming

With --bounding-box, fitscube-rs trims away the blank border shared by every plane. It computes, across all input planes, the smallest rectangle that contains all the valid (non-blank) pixels, then crops every plane to that common box and updates the spatial reference pixel (CRPIX1/2) so the WCS stays correct. This shrinks cubes that were padded out to a large common canvas without discarding any real data. --invalidate-zeros first treats exact-zero pixels as blank, so zero-padded borders are trimmed too.

Floating-point precision

--floating {32,64} selects the pixel data type of the output cube (float32 / BITPIX −32, or float64 / BITPIX −64). These are the only IEEE float widths the FITS standard defines, so other values are rejected. The default follows the inputs; downcasting to float32 is offered for cubes where storage matters more than the last bits of precision.

Plane extraction

extract is the inverse operation. Given a cube and a plane selector (--channel-index for spectral cubes, --time-index for time cubes), it reads that single plane, rebuilds a 2D image header from the cube WCS (dropping the combined axis and restoring the per-plane coordinate and beam where available), and writes a standalone FITS image. --hdu-index selects which HDU of the cube to read from when it is not the primary.