Merge pull request #112 from effigies/feat/axiscodes

feat: Extract image orientation from NIfTI header

Author: effigies

changelog.d/20250903_155423_markiewicz_axiscodes.md

@@ -0,0 +1,48 @@
+<!--
+A new scriv changelog fragment.
+
+Uncomment the section that is right (remove the HTML comment wrapper).
+For top level release notes, leave all the headers commented out.
+-->
+
+### Added
+
+- Added support for extracting image orientation from NIfTI headers,
+  added to the BIDS schema in 1.10.1.
+
+<!--
+### Changed
+
+- A bullet item for the Changed category.
+
+-->
+<!--
+### Fixed
+
+- A bullet item for the Fixed category.
+
+-->
+<!--
+### Deprecated
+
+- A bullet item for the Deprecated category.
+
+-->
+<!--
+### Removed
+
+- A bullet item for the Removed category.
+
+-->
+<!--
+### Security
+
+- A bullet item for the Security category.
+
+-->
+<!--
+### Infrastructure
+
+- A bullet item for the Infrastructure category.
+
+-->

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([])
@@ -73,3 +73,26 @@ Deno.test('Test loading nifti header', async (t) => {
     })
   })
 })
+
+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

@@ -65,8 +65,91 @@ export async function loadHeader(file: BIDSFile): Promise<NiftiHeader> {
       },
       qform_code: header.qform_code,
       sform_code: header.sform_code,
+      axis_codes: axisCodes(header.affine),
     } as NiftiHeader
   } catch (err) {
     throw { code: '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])
+}

src/schema/applyRules.ts

@@ -146,7 +146,7 @@ function mapEvalCheck(statements: string[], context: BIDSContext): boolean {
  * Classic rules interpreted like selectors. Examples in specification:
  * schema/rules/checks/*
  */
-function evalRuleChecks(
+export function evalRuleChecks(
   rule: GenericRule,
   context: BIDSContext,
   schema: GenericSchema,

src/tests/local/nifti_rules.test.ts

@@ -0,0 +1,111 @@
+import { assertEquals, assertObjectMatch } from '@std/assert'
+import type { BIDSFile, FileTree } from '../../types/filetree.ts'
+import type { GenericSchema, GenericRule } from '../../types/schema.ts'
+import { loadHeader } from '../../files/nifti.ts'
+import { BIDSFileDeno } from '../../files/deno.ts'
+import { BIDSContextDataset, BIDSContext } from '../../schema/context.ts'
+import { loadSchema } from '../../setup/loadSchema.ts'
+import { evalRuleChecks } from '../../schema/applyRules.ts'
+// import { applyRules } from '../../schema/applyRules.ts'
+import type { Context, NiftiHeader } from '@bids/schema/context'
+import type { Schema } from '@bids/schema/metaschema'
+
+import { expressionFunctions } from '../../schema/expressionLanguage.ts'
+
+function prepContext(header: NiftiHeader, dir: string, pedir: string): BIDSContext {
+  const fullContext = {
+    dataset: new BIDSContextDataset({}),
+    nifti_header: header,
+    entities: {direction: dir},
+    sidecar: {PhaseEncodingDirection: pedir},
+  } as unknown as BIDSContext
+  Object.assign(fullContext, expressionFunctions)
+  return fullContext
+}
+
+Deno.test('Test NIFTI-specific rules', async (t) => {
+  const RAS = await loadHeader(new BIDSFileDeno('', 'tests/data/RAS.nii.gz'))
+  const SPL = await loadHeader(new BIDSFileDeno('', 'tests/data/SPL.nii.gz'))
+  const AIR = await loadHeader(new BIDSFileDeno('', 'tests/data/AIR.nii.gz'))
+
+  const schema = await loadSchema() as Schema
+  const NiftiPEDir = schema.rules?.checks?.nifti?.NiftiPEDir as GenericRule
+
+  await t.step('Test reading NIfTI axis codes' , async () => {
+    assertEquals(RAS.axis_codes, ['R', 'A', 'S'])
+    assertEquals(SPL.axis_codes, ['S', 'P', 'L'])
+    assertEquals(AIR.axis_codes, ['A', 'I', 'R'])
+  })
+
+  await t.step('Test rules.checks.nifti.NiftiPEDir' , async () => {
+    const schemaPath = 'rules.checks.nifti.NiftiPEDir'
+
+    let context = prepContext(RAS, 'PA', 'j')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(RAS, 'AP', 'j-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(RAS, 'LR', 'i')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(RAS, 'RL', 'i-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(RAS, 'IS', 'k')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(RAS, 'SI', 'k-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+
+    // Common flips
+    context = prepContext(RAS, 'AP', 'j')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'PA', 'j-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'RL', 'i')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'LR', 'i-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+
+    // Wrong axes
+    context = prepContext(RAS, 'PA', 'i')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'AP', 'k')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'LR', 'j')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+    context = prepContext(RAS, 'RL', 'k-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({code: 'NIFTI_PE_DIRECTION_CONSISTENCY' }).length, 1)
+
+    // A couple checks on SPL and AIR
+    context = prepContext(SPL, 'IS', 'i')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(SPL, 'AP', 'j')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(SPL, 'RL', 'k')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+
+    context = prepContext(AIR, 'IS', 'j-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(AIR, 'AP', 'i-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+    context = prepContext(AIR, 'RL', 'k-')
+    evalRuleChecks(NiftiPEDir, context, {} as GenericSchema, schemaPath)
+    assertEquals(context.dataset.issues.get({}).length, 0)
+  })
+})