feat: Extract image orientation from NIfTI header

src/files/nifti.test.ts

@@ -1,8 +1,8 @@
-import { assert, assertObjectMatch } from '@std/assert'
+import { assert, assertEquals, assertObjectMatch } from '@std/assert'
 import { FileIgnoreRules } from './ignore.ts'
 import { BIDSFileDeno } from './deno.ts'
 
-import { loadHeader } from './nifti.ts'
+import { loadHeader, axisCodes } from './nifti.ts'
 
 Deno.test('Test loading nifti header', async (t) => {
   const ignore = new FileIgnoreRules([])
@@ -54,3 +54,26 @@ Deno.test('Test loading nifti header', async (t) => {
     assertObjectMatch(error, { key: 'NIFTI_HEADER_UNREADABLE' })
   })
 })
+
+Deno.test('Test extracting axis codes', async (t) => {
+  await t.step('Identify RAS', async () => {
+    const affine = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), ['R', 'A', 'S'])
+  })
+  await t.step('Identify LPS (flips)', async () => {
+    const affine = [[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), ['L', 'P', 'S'])
+  })
+  await t.step('Identify SPL (flips + swap)', async () => {
+    const affine = [[0, 0, -1, 0], [0, -1, 0, 0], [1, 0, 0, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), ['S', 'P', 'L'])
+  })
+  await t.step('Identify SLP (flips + rotate)', async () => {
+    const affine = [[0, -1, 0, 0], [0, 0, -1, 0], [1, 0, 0, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), ['S', 'L', 'P'])
+  })
+  await t.step('Identify ASR (rotate)', async () => {
+    const affine = [[0, 0, 1, 0], [1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1]]
+    assertEquals(axisCodes(affine), ['A', 'S', 'R'])
+  })
+})

src/files/nifti.ts

@@ -58,3 +58,85 @@ export async function loadHeader(file: BIDSFile): Promise<NiftiHeader> {
     throw { key: 'NIFTI_HEADER_UNREADABLE' }
   }
 }
+
+/** Vector addition */
+function add(a: number[], b: number[]): number[] {
+  return a.map((x, i) => x + b[i])
+}
+
+/** Vector subtraction */
+function sub(a: number[], b: number[]): number[] {
+  return a.map((x, i) => x - b[i])
+}
+
+/** Scalar multiplication */
+function scale(vec: number[], scalar: number): number[] {
+  return vec.map((x) => x * scalar)
+}
+
+/** Dot product */
+function dot(a: number[], b: number[]): number {
+  return a.map((x, i) => x * b[i]).reduce((acc, x) => acc + x, 0)
+}
+
+function argMax(arr: number[]): number {
+  return arr.reduce((acc, x, i) => (x > arr[acc] ? i : acc), 0)
+}
+
+/**
+ * Identify the nearest principle axes of an image affine.
+ *
+ * Affines transform indices in a data array into mm right, anterior and superior of
+ * an origin in "world coordinates". If moving along an axis in the positive direction
+ * predominantly moves right, that axis is labeled "R".
+ *
+ * @example The identity matrix is in "RAS" orientation:
+ *
+ * # Usage
+ *
+ * ```ts
+ * const affine = [[1, 0, 0, 0],
+ *                 [0, 1, 0, 0],
+ *                 [0, 0, 1, 0],
+ *                 [0, 0, 0, 1]]
+ *
+ * axisCodes(affine)
+ * ```
+ *
+ * # Result
+ * ```ts
+ * ['R', 'A', 'S']
+ * ```
+ *
+ * @returns character codes describing the orientation of an image affine.
+ */
+export function axisCodes(affine: number[][]): string[] {
+  // 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.
+
+  // 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]))
+
+  // Orthogonalize cosY with respect to cosX
+  const orthY = sub(cosY, scale(cosX, dot(cosX, cosY)))
+
+  // Orthogonalize cosZ with respect to cosX and orthY
+  const orthZ = sub(
+    cosZ, add(scale(cosX, dot(cosX, cosZ)), scale(orthY, dot(orthY, cosZ)))
+  )
+
+  const basis = [cosX, orthY, orthZ]
+  const maxIndices = basis.map((row) => argMax(row.map(Math.abs)))
+
+  // Check that indices are 0, 1 and 2 in some order
+  if (maxIndices.toSorted().some((idx, i) => idx !== i)) {
+    throw { key: 'AMBIGUOUS_AFFINE' }
+  }
+
+  // Positive/negative codes for each world axis
+  const codes = ['RL', 'AP', 'SI']
+  return maxIndices.map((idx, i) => codes[idx][basis[i][idx] > 0 ? 0 : 1])
+}