fix: Set nifti_header.axis_codes to null for bad qforms, rather than error

changelog.d/20251124_080614_markiewicz_bad_qform.md

@@ -0,0 +1,5 @@
+### Fixed
+
+- NIfTI files with bad qform matrices, resulting from non-normalized quaternions,
+  would previously raise a NIFTI_HEADER_UNREADABLE error. Now only the axis codes
+  are disabled, preventing orientation checks, but not raising errors.

src/files/nifti.test.ts

@@ -96,6 +96,10 @@ Deno.test('Test extracting axis codes', async (t) => {
     const affine = [[0, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]]
     assertEquals(axisCodes(affine), ['A', 'S', 'R'])
   })
+  await t.step('Fail gracefully on NaNs', async () => {
+    const affine = [[Number.NaN, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), null)
+  })
 })
 
 testAsyncFileAccess('Test file access errors for loadHeader', loadHeader)

src/files/nifti.ts

@@ -34,9 +34,9 @@ async function extract(buffer: Uint8Array, nbytes: number): Promise<Uint8Array<A
 
 export async function loadHeader(file: BIDSFile): Promise<NiftiHeader> {
   const buf = await readBytes(file, 1024)
+  let header
   try {
     const data = isCompressed(buf.buffer) ? await extract(buf, 540) : buf.slice(0, 540)
-    let header
     if (isNIFTI1(data.buffer)) {
       header = new NIFTI1()
       // Truncate to 348 bytes to avoid attempting to parse extensions
@@ -48,29 +48,30 @@ export async function loadHeader(file: BIDSFile): Promise<NiftiHeader> {
     if (!header) {
       throw { code: 'NIFTI_HEADER_UNREADABLE' }
     }
-    const ndim = header.dims[0]
-    return {
-      dim: header.dims,
-      // Hack: round pixdim to 3 decimal places; schema should add rounding function
-      pixdim: header.pixDims.map((pixdim) => Math.round(pixdim * 1000) / 1000),
-      shape: header.dims.slice(1, ndim + 1),
-      voxel_sizes: header.pixDims.slice(1, ndim + 1),
-      dim_info: {
-        freq: header.dim_info & 0x03,
-        phase: (header.dim_info >> 2) & 0x03,
-        slice: (header.dim_info >> 4) & 0x03,
-      },
-      xyzt_units: {
-        xyz: ['unknown', 'meter', 'mm', 'um'][header.xyzt_units & 0x03],
-        t: ['unknown', 'sec', 'msec', 'usec'][(header.xyzt_units >> 3) & 0x03],
-      },
-      qform_code: header.qform_code,
-      sform_code: header.sform_code,
-      axis_codes: axisCodes(header.affine),
-    } as NiftiHeader
   } catch (err) {
     throw { code: 'NIFTI_HEADER_UNREADABLE' }
   }
+
+  const ndim = header.dims[0]
+  return {
+    dim: header.dims,
+    // Hack: round pixdim to 3 decimal places; schema should add rounding function
+    pixdim: header.pixDims.map((pixdim) => Math.round(pixdim * 1000) / 1000),
+    shape: header.dims.slice(1, ndim + 1),
+    voxel_sizes: header.pixDims.slice(1, ndim + 1),
+    dim_info: {
+      freq: header.dim_info & 0x03,
+      phase: (header.dim_info >> 2) & 0x03,
+      slice: (header.dim_info >> 4) & 0x03,
+    },
+    xyzt_units: {
+      xyz: ['unknown', 'meter', 'mm', 'um'][header.xyzt_units & 0x03],
+      t: ['unknown', 'sec', 'msec', 'usec'][(header.xyzt_units >> 3) & 0x03],
+    },
+    qform_code: header.qform_code,
+    sform_code: header.sform_code,
+    axis_codes: axisCodes(header.affine),
+  } as NiftiHeader
 }
 
 /** Vector addition */
@@ -124,13 +125,18 @@ function argMax(arr: number[]): number {
  *
  * @returns character codes describing the orientation of an image affine.
  */
-export function axisCodes(affine: number[][]): string[] {
+export function axisCodes(affine: number[][]): string[] | null {
   // This function is an extract of the Python function transforms3d.affines.decompose44
   // (https://github.com/matthew-brett/transforms3d/blob/6a43a98/transforms3d/affines.py#L10-L153)
   //
   // As an optimization, this only orthogonalizes the basis,
   // and does not normalize to unit vectors.
 
+  // Bad qforms result in NaNs in the rotation matrix
+  if (affine.some((row) => row.some((val) => !Number.isFinite(val)))) {
+    return null
+  }
+
   // Operate on columns, which are the cosines that project input coordinates onto output axes
   const [cosX, cosY, cosZ] = [0, 1, 2].map((j) => [0, 1, 2].map((i) => affine[i][j]))
 
@@ -148,7 +154,7 @@ export function axisCodes(affine: number[][]): string[] {
 
   // Check that indices are 0, 1 and 2 in some order
   if (maxIndices.toSorted().some((idx, i) => idx !== i)) {
-    throw { key: 'AMBIGUOUS_AFFINE' }
+    throw { code: 'AMBIGUOUS_AFFINE' }
   }
 
   // Positive/negative codes for each world axis