diff --git a/changelog.d/20250903_155423_markiewicz_axiscodes.md b/changelog.d/20250903_155423_markiewicz_axiscodes.md new file mode 100644 index 00000000..7dcab4ee --- /dev/null +++ b/changelog.d/20250903_155423_markiewicz_axiscodes.md @@ -0,0 +1,48 @@ + + +### Added + +- Added support for extracting image orientation from NIfTI headers, + added to the BIDS schema in 1.10.1. + + + + + + + diff --git a/src/files/nifti.test.ts b/src/files/nifti.test.ts index a4d00aff..e3fed4bf 100644 --- a/src/files/nifti.test.ts +++ b/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']) + }) +}) diff --git a/src/files/nifti.ts b/src/files/nifti.ts index 774ee6fa..5639bb2d 100644 --- a/src/files/nifti.ts +++ b/src/files/nifti.ts @@ -65,8 +65,91 @@ export async function loadHeader(file: BIDSFile): Promise { }, 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]) +} diff --git a/src/schema/applyRules.ts b/src/schema/applyRules.ts index 1ca37b8b..28337ae4 100644 --- a/src/schema/applyRules.ts +++ b/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, diff --git a/src/tests/local/nifti_rules.test.ts b/src/tests/local/nifti_rules.test.ts new file mode 100644 index 00000000..376d20d6 --- /dev/null +++ b/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) + }) +}) diff --git a/tests/data/AIR.nii.gz b/tests/data/AIR.nii.gz new file mode 100644 index 00000000..c52148e5 Binary files /dev/null and b/tests/data/AIR.nii.gz differ diff --git a/tests/data/RAS.nii.gz b/tests/data/RAS.nii.gz new file mode 100644 index 00000000..8e8f40e9 Binary files /dev/null and b/tests/data/RAS.nii.gz differ diff --git a/tests/data/SPL.nii.gz b/tests/data/SPL.nii.gz new file mode 100644 index 00000000..24097180 Binary files /dev/null and b/tests/data/SPL.nii.gz differ