feat: Calculate axis codes from affines

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,111 @@ export async function loadHeader(file: BIDSFile): Promise<NiftiHeader> {
     throw { key: 'NIFTI_HEADER_UNREADABLE' }
   }
 }
+
+function add(a: number[], b: number[]): number[] {
+  return a.map((x, i) => x + b[i])
+}
+
+function sub(a: number[], b: number[]): number[] {
+  return a.map((x, i) => x - b[i])
+}
+
+function scale(vec: number[], scalar: number): number[] {
+  return vec.map((x) => x * scalar)
+}
+
+function dot(a: number[], b: number[]): number {
+  return a.map((x, i) => x * b[i]).reduce((acc, x) => acc + x, 0)
+}
+
+function extractRotation(affine: number[][]): number[][] {
+  // 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)
+  //
+  // To explain the conventions of the s{xyz}* parameters:
+  //
+  // The upper left 3x3 of the affine is a matrix we will call RZS which can be decomposed
+  //
+  //    RZS = R * Z * S
+  //
+  // where R is a 3x3 rotation matrix, Z is a diagonal matrix of scalings
+  //
+  //    Z = diag([sx, xy, sz])
+  //
+  // and S is a shear matrix with the form
+  //
+  //    S = [[1, sxy, sxz],
+  //         [0,   1, syz],
+  //         [0,   0,   1]]
+  //
+  // Note that this function does not return scales, shears or translations, and
+  // does not guarantee a right-handed rotation matrix, as that is not necessary for our use.
+
+  // 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]))
+
+  const sx = Math.sqrt(dot(cosX, cosX))
+  const normX = cosX.map((x) => x / sx) // Unit vector
+
+  // Orthogonalize cosY with respect to normX
+  const sx_sxy = dot(normX, cosY)
+  const orthY = sub(cosY, scale(normX, sx_sxy))
+  const sy = Math.sqrt(dot(orthY, orthY))
+  const normY = orthY.map((y) => y / sy)
+
+  // Orthogonalize cosZ with respect to normX and normY
+  const sx_sxz = dot(normX, cosZ)
+  const sy_syz = dot(normY, cosZ)
+  const orthZ = sub(cosZ, add(scale(normX, sx_sxz), scale(normY, sy_syz)))
+  const sz = Math.sqrt(dot(orthZ, orthZ))
+  const normZ = orthZ.map((z) => z / sz)
+
+  // Transposed normalized cosines
+  return [normX, normY, normZ]
+}
+
+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[] {
+  // Note that rotation is transposed
+  const rotations = extractRotation(affine)
+  const maxIndices = rotations.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][rotations[i][idx] > 0 ? 0 : 1])
+}