From e94150cad8187d8a13e01ab94aa404fc5b91297c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sat, 1 Nov 2025 10:34:30 -0400 Subject: [PATCH 01/13] ENH: Add transform analysis utils Add transform analysis utils. Transfer contents from the `NiFreeze` projects so that hey can be reused across projects requiring transform analysis: https://github.com/nipreps/nifreeze/blob/d27ba7552bbd9095c3c13b46443d87a4b5504c4c/src/nifreeze/analysis/motion.py https://github.com/nipreps/nifreeze/blob/d27ba7552bbd9095c3c13b46443d87a4b5504c4c/src/nifreeze/data/utils.py Add a fixture to be able to reuse a random number generator across tests. Co-authored-by: Oscar Esteban --- nitransforms/analysis/__init__.py | 0 nitransforms/analysis/utils.py | 188 ++++++++++++++++++++++++++++ nitransforms/tests/conftest.py | 11 ++ nitransforms/tests/test_analysis.py | 156 +++++++++++++++++++++++ 4 files changed, 355 insertions(+) create mode 100644 nitransforms/analysis/__init__.py create mode 100644 nitransforms/analysis/utils.py create mode 100644 nitransforms/tests/conftest.py create mode 100644 nitransforms/tests/test_analysis.py diff --git a/nitransforms/analysis/__init__.py b/nitransforms/analysis/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py new file mode 100644 index 00000000..8b627cd5 --- /dev/null +++ b/nitransforms/analysis/utils.py @@ -0,0 +1,188 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Utilities to aid in performing and evaluating image registration. + +This module provides functions to compute displacements of image coordinates +under a transformation, useful for assessing the accuracy of image registration +processes. + +""" + +from __future__ import annotations + +from itertools import product +from typing import Tuple + +import nibabel as nb +import numpy as np +from scipy.stats import zscore + +from nitransforms.base import TransformBase + + +RADIUS = 50.0 +"""Typical radius (in mm) of a sphere mimicking the size of a typical human brain.""" + + +def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = RADIUS) -> np.ndarray: + """Compute framewise displacement (FD) from motion parameters. + + Each row in the motion parameters represents one frame, and columns + represent each coordinate axis ``x``, `y``, and ``z``. Translation + parameters are followed by rotation parameters column-wise. + + Parameters + ---------- + motion_parameters : :obj:`~numpy.ndarray` + Motion parameters. + radius : :obj:`float`, optional + Radius (in mm) of a sphere mimicking the size of a typical human brain. + + Returns + ------- + :obj:`~numpy.ndarray` + The framewise displacement (FD) as the sum of absolute differences + between consecutive frames. + """ + + translations = motion_parameters[:, :3] + rotations_deg = motion_parameters[:, 3:] + rotations_rad = np.deg2rad(rotations_deg) + + # Compute differences between consecutive frames + d_translations = np.vstack([np.zeros((1, 3)), np.diff(translations, axis=0)]) + d_rotations = np.vstack([np.zeros((1, 3)), np.diff(rotations_rad, axis=0)]) + + # Convert rotations from radians to displacement on a sphere + rotation_displacement = d_rotations * radius + + # Compute FD as sum of absolute differences + return np.sum(np.abs(d_translations) + np.abs(rotation_displacement), axis=1) + + +def compute_fd_from_transform( + img: nb.spatialimages.SpatialImage, + test_xfm: TransformBase, + radius: float = RADIUS, +) -> float: + """ + Compute the framewise displacement (FD) for a given transformation. + + Parameters + ---------- + img : :obj:`~nibabel.spatialimages.SpatialImage` + The reference image. Used to extract the center coordinates. + test_xfm : :obj:`~nitransforms.base.TransformBase` + The transformation to test. Applied to coordinates around the image center. + radius : :obj:`float`, optional + The radius (in mm) of the spherical neighborhood around the center of the image. + + Returns + ------- + :obj:`float` + The average framewise displacement (FD) for the test transformation. + + """ + affine = img.affine + # Compute the center of the image in voxel space + center_ijk = 0.5 * (np.array(img.shape[:3]) - 1) + # Convert to world coordinates + center_xyz = nb.affines.apply_affine(affine, center_ijk) + # Generate coordinates of points at radius distance from center + fd_coords = np.array(list(product(*((radius, -radius),) * 3))) + center_xyz + # Compute the average displacement from the test transformation + return np.mean(np.linalg.norm(test_xfm.map(fd_coords) - fd_coords, axis=-1)) + + +def displacements_within_mask( + mask_img: nb.spatialimages.SpatialImage, + test_xfm: TransformBase, + reference_xfm: TransformBase | None = None, +) -> np.ndarray: + """ + Compute the distance between voxel coordinates mapped through two transforms. + + Parameters + ---------- + mask_img : :obj:`~nibabel.spatialimages.SpatialImage` + A mask image that defines the region of interest. Voxel coordinates + within the mask are transformed. + test_xfm : :obj:`~nitransforms.base.TransformBase` + The transformation to test. This transformation is applied to the + voxel coordinates. + reference_xfm : :obj:`~nitransforms.base.TransformBase`, optional + A reference transformation to compare with. If ``None``, the identity + transformation is assumed (no transformation). + + Returns + ------- + :obj:`~numpy.ndarray` + An array of displacements (in mm) for each voxel within the mask. + + """ + # Mask data as boolean (True for voxels inside the mask) + maskdata = np.asanyarray(mask_img.dataobj) > 0 + # Convert voxel coordinates to world coordinates using affine transform + xyz = nb.affines.apply_affine( + mask_img.affine, + np.argwhere(maskdata), + ) + # Apply the test transformation + targets = test_xfm.map(xyz) + + # Compute the difference (displacement) between the test and reference transformations + diffs = targets - xyz if reference_xfm is None else targets - reference_xfm.map(xyz) + return np.linalg.norm(diffs, axis=-1) + + +def extract_motion_parameters(affine: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: + """Extract translation (mm) and rotation (degrees) parameters from an affine matrix. + + Parameters + ---------- + affine : :obj:`~numpy.ndarray` + The affine transformation matrix. + + Returns + ------- + :obj:`tuple` + Extracted translation and rotation parameters. + """ + + translation = affine[:3, 3] + rotation_rad = np.arctan2( + [affine[2, 1], affine[0, 2], affine[1, 0]], [affine[2, 2], affine[0, 0], affine[1, 1]] + ) + rotation_deg = np.rad2deg(rotation_rad) + return *translation, *rotation_deg + + +def identify_spikes(fd: np.ndarray, threshold: float = 2.0) -> Tuple[np.ndarray, np.ndarray]: + """Identify motion spikes in framewise displacement data. + + Identifies high-motion frames as timepoint exceeding a given threshold value + based on z-score normalized framewise displacement (FD) values. + + Parameters + ---------- + fd : :obj:`~numpy.ndarray` + Framewise displacement data. + threshold : :obj:`float`, optional + Threshold value to determine motion spikes. + + Returns + ------- + indices : :obj:`~numpy.ndarray` + Indices of identified motion spikes. + mask : :obj:`~numpy.ndarray` + Mask of identified motion spikes. + """ + + # Normalize (z-score) + fd_norm = zscore(fd) + + mask = fd_norm > threshold + indices = np.where(mask)[0] + + return indices, mask diff --git a/nitransforms/tests/conftest.py b/nitransforms/tests/conftest.py new file mode 100644 index 00000000..380536c1 --- /dev/null +++ b/nitransforms/tests/conftest.py @@ -0,0 +1,11 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +import numpy as np +import pytest + + +@pytest.fixture(autouse=True) +def random_number_generator(request): + """Automatically set a fixed-seed random number generator for all tests.""" + request.node.rng = np.random.default_rng(1234) diff --git a/nitransforms/tests/test_analysis.py b/nitransforms/tests/test_analysis.py new file mode 100644 index 00000000..685ba57b --- /dev/null +++ b/nitransforms/tests/test_analysis.py @@ -0,0 +1,156 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: + +import numpy as np +import nibabel as nb +import pytest + +import nitransforms as nt + +from nitransforms.analysis.utils import ( + compute_fd_from_motion, + compute_fd_from_transform, + displacements_within_mask, + extract_motion_parameters, + identify_spikes, +) + + +@pytest.fixture +def identity_affine(): + return np.eye(4) + + +@pytest.fixture +def simple_mask_img(identity_affine): + # 3x3x3 mask with center voxel as 1, rest 0 + data = np.zeros((3, 3, 3), dtype=np.uint8) + data[1, 1, 1] = 1 + return nb.Nifti1Image(data, identity_affine) + + +@pytest.fixture +def translation_transform(): + # Simple translation of (1, 2, 3) mm + return nt.linear.Affine(map=np.array([ + [1, 0, 0, 1], + [0, 1, 0, 2], + [0, 0, 1, 3], + [0, 0, 0, 1], + ])) + + +@pytest.fixture +def rotation_transform(): + # 90 degree rotation around z axis + angle = np.pi / 2 + rot = np.array([ + [np.cos(angle), -np.sin(angle), 0, 0], + [np.sin(angle), np.cos(angle), 0, 0], + [0, 0, 1, 0], + [0, 0, 0, 1], + ]) + return nt.linear.Affine(map=rot) + + +@pytest.mark.parametrize( + "test_xfm, reference_xfm, expected", + [ + (nt.linear.Affine(np.eye(4)), None, np.zeros(1)), + (nt.linear.Affine(np.array([ + [1, 0, 0, 1], + [0, 1, 0, 2], + [0, 0, 1, 3], + [0, 0, 0, 1], + ])), None, [np.linalg.norm([1, 2, 3])]), + (nt.linear.Affine(np.array([ + [1, 0, 0, 1], + [0, 1, 0, 2], + [0, 0, 1, 3], + [0, 0, 0, 1], + ])), nt.linear.Affine(np.eye(4)), [np.linalg.norm([1, 2, 3])]), + ], +) +def test_displacements_within_mask(simple_mask_img, test_xfm, reference_xfm, expected): + disp = displacements_within_mask(simple_mask_img, test_xfm, reference_xfm) + np.testing.assert_allclose(disp, expected) + + +@pytest.mark.parametrize( + "test_xfm, expected", + [ + (nt.linear.Affine(np.eye(4)), 0), + (nt.linear.Affine(np.array([ + [1, 0, 0, 1], + [0, 1, 0, 2], + [0, 0, 1, 3], + [0, 0, 0, 1], + ])), np.linalg.norm([1, 2, 3])), + ], +) +def test_compute_fd_from_transform(simple_mask_img, test_xfm, expected): + fd = compute_fd_from_transform(simple_mask_img, test_xfm) + assert np.isclose(fd, expected) + + +@pytest.mark.parametrize( + "motion_params, radius, expected", + [ + (np.zeros((5, 6)), 50, np.zeros(5)), # 5 frames, 3 trans, 3 rot + ( + np.array([ + [0,0,0,0,0,0], + [2,0,0,0,0,0], # 2mm translation in x at frame 1 + [2,0,0,90,0,0], + ]), # 90deg rotation in x at frame 2 + 50, + [0, 2, abs(np.deg2rad(90)) * 50] + ), # First frame: 0, Second: translation 2mm, Third: rotation (pi/2)*50 + ], +) +def test_compute_fd_from_motion(motion_params, radius, expected): + fd = compute_fd_from_motion(motion_params, radius=radius) + np.testing.assert_allclose(fd, expected, atol=1e-4) + + +@pytest.mark.parametrize( + "affine, expected_trans, expected_rot", + [ + (np.eye(4) + np.array([[0,0,0,10],[0,0,0,15],[0,0,0,20],[0,0,0,0]]), # translation only + [10, 15, 20], [0, 0, 0]), + (np.array([ + [1, 0, 0, 0], + [0, np.cos(np.deg2rad(30)), -np.sin(np.deg2rad(30)), 0], + [0, np.sin(np.deg2rad(30)), np.cos(np.deg2rad(30)), 0], + [0, 0, 0, 1], # rotation only + ]), [0, 0, 0], [30, 0, 0]), # Only one rot will be close to 30 + ], +) +def test_extract_motion_parameters(affine, expected_trans, expected_rot): + params = extract_motion_parameters(affine) + assert np.allclose(params[:3], expected_trans) + # For rotation case, at least one value close to 30 + if np.any(np.abs(expected_rot)): + assert np.any(np.isclose(np.abs(params[3:]), 30, atol=1e-4)) + else: + assert np.allclose(params[3:], expected_rot) + + +def test_identify_spikes(request): + rng = request.node.rng + + n_samples = 450 + + fd = rng.normal(0, 5, n_samples) + threshold = 2.0 + + expected_indices = np.asarray( + [5, 57, 85, 100, 127, 180, 191, 202, 335, 393, 409] + ) + expected_mask = np.zeros(n_samples, dtype=bool) + expected_mask[expected_indices] = True + + obtained_indices, obtained_mask = identify_spikes(fd, threshold=threshold) + + assert np.array_equal(obtained_indices, expected_indices) + assert np.array_equal(obtained_mask, expected_mask) From fe6816265cb3d9616daac26554314cba5718bde2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sat, 3 Jan 2026 14:37:52 -0500 Subject: [PATCH 02/13] Apply suggestions from code review Co-authored-by: Oscar Esteban --- nitransforms/analysis/utils.py | 17 +++++++++++------ nitransforms/tests/test_analysis.py | 4 ++-- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 8b627cd5..fbe5999c 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -21,11 +21,16 @@ from nitransforms.base import TransformBase -RADIUS = 50.0 -"""Typical radius (in mm) of a sphere mimicking the size of a typical human brain.""" +DEFAULT_FD_RADIUS = 50.0 +""" +Default radius (in mm) of a sphere where framewise displacements are calculated. +The choice was proposed by +`Power et al. (2015) `__, and it +represents approximately the mean distance from the cerebral cortex to the center of the head. +""" -def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = RADIUS) -> np.ndarray: +def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = DEFAULT_FD_RADIUS) -> np.ndarray: """Compute framewise displacement (FD) from motion parameters. Each row in the motion parameters represents one frame, and columns @@ -158,7 +163,7 @@ def extract_motion_parameters(affine: np.ndarray) -> Tuple[np.ndarray, np.ndarra return *translation, *rotation_deg -def identify_spikes(fd: np.ndarray, threshold: float = 2.0) -> Tuple[np.ndarray, np.ndarray]: +def identify_spikes(fd: np.ndarray, z_threshold: float = 2.0) -> Tuple[np.ndarray, np.ndarray]: """Identify motion spikes in framewise displacement data. Identifies high-motion frames as timepoint exceeding a given threshold value @@ -168,7 +173,7 @@ def identify_spikes(fd: np.ndarray, threshold: float = 2.0) -> Tuple[np.ndarray, ---------- fd : :obj:`~numpy.ndarray` Framewise displacement data. - threshold : :obj:`float`, optional + z_threshold : :obj:`float`, optional Threshold value to determine motion spikes. Returns @@ -182,7 +187,7 @@ def identify_spikes(fd: np.ndarray, threshold: float = 2.0) -> Tuple[np.ndarray, # Normalize (z-score) fd_norm = zscore(fd) - mask = fd_norm > threshold + mask = fd_norm > z_threshold indices = np.where(mask)[0] return indices, mask diff --git a/nitransforms/tests/test_analysis.py b/nitransforms/tests/test_analysis.py index 685ba57b..3307b8b4 100644 --- a/nitransforms/tests/test_analysis.py +++ b/nitransforms/tests/test_analysis.py @@ -142,7 +142,7 @@ def test_identify_spikes(request): n_samples = 450 fd = rng.normal(0, 5, n_samples) - threshold = 2.0 + z_threshold = 2.0 expected_indices = np.asarray( [5, 57, 85, 100, 127, 180, 191, 202, 335, 393, 409] @@ -150,7 +150,7 @@ def test_identify_spikes(request): expected_mask = np.zeros(n_samples, dtype=bool) expected_mask[expected_indices] = True - obtained_indices, obtained_mask = identify_spikes(fd, threshold=threshold) + obtained_indices, obtained_mask = identify_spikes(fd, z_threshold=z_threshold) assert np.array_equal(obtained_indices, expected_indices) assert np.array_equal(obtained_mask, expected_mask) From 7016e0b42300eff2eb58a89fcc2b77924e0cd2ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sat, 3 Jan 2026 14:50:57 -0500 Subject: [PATCH 03/13] DOC: Cite Power et al. in FD computation function docstring Cite Power et al. in FD computation function docstring. --- nitransforms/analysis/utils.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index fbe5999c..5bfa5769 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -33,6 +33,10 @@ def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = DEFAULT_FD_RADIUS) -> np.ndarray: """Compute framewise displacement (FD) from motion parameters. + The framewise displacement is the sum of the magnitudes of the translational + and rotational motion, computed from the frame-to-frame differences along + the three spatial axes (`Power et al. (2015) `__). + Each row in the motion parameters represents one frame, and columns represent each coordinate axis ``x``, `y``, and ``z``. Translation parameters are followed by rotation parameters column-wise. From 9c48c5faba32c1bab97fb0924da8fbe7c5668667 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sat, 3 Jan 2026 14:54:09 -0500 Subject: [PATCH 04/13] Update nitransforms/analysis/utils.py --- nitransforms/analysis/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 5bfa5769..54babcde 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -73,7 +73,7 @@ def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = DEFAUL def compute_fd_from_transform( img: nb.spatialimages.SpatialImage, test_xfm: TransformBase, - radius: float = RADIUS, + radius: float = DEFAULT_FD_RADIUS, ) -> float: """ Compute the framewise displacement (FD) for a given transformation. From 43cafbc187b75bf57b859a632056ac0977ecd4ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sat, 3 Jan 2026 15:01:17 -0500 Subject: [PATCH 05/13] STY: Fix style errors Fix style errors. --- nitransforms/analysis/utils.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 54babcde..2449cc86 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -30,12 +30,16 @@ """ -def compute_fd_from_motion(motion_parameters: np.ndarray, radius: float = DEFAULT_FD_RADIUS) -> np.ndarray: +def compute_fd_from_motion( + motion_parameters: np.ndarray, + radius: float = DEFAULT_FD_RADIUS, +) -> np.ndarray: """Compute framewise displacement (FD) from motion parameters. The framewise displacement is the sum of the magnitudes of the translational and rotational motion, computed from the frame-to-frame differences along - the three spatial axes (`Power et al. (2015) `__). + the three spatial axes ( + `Power et al. (2015) `__). Each row in the motion parameters represents one frame, and columns represent each coordinate axis ``x``, `y``, and ``z``. Translation From ee16315e14af22326ed8d2a9ac4477c79a6a6d47 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Sun, 4 Jan 2026 16:09:47 +0100 Subject: [PATCH 06/13] fix: generalize sampling on unit sphere computation --- nitransforms/analysis/utils.py | 175 +++++++++++++++++++++++++++++++-- 1 file changed, 167 insertions(+), 8 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 2449cc86..7d3c0ffa 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -7,10 +7,18 @@ under a transformation, useful for assessing the accuracy of image registration processes. +References +---------- +.. [Power2012] Power, JD. et al. (2012). "Spurious but systematic correlations in functional + connectivity MRI networks arise from subject motion." NeuroImage, 59(3):2142-2154. + doi:`10.1016/j.neuroimage.2011.10.018 `__. + """ from __future__ import annotations + +import math from itertools import product from typing import Tuple @@ -24,9 +32,8 @@ DEFAULT_FD_RADIUS = 50.0 """ Default radius (in mm) of a sphere where framewise displacements are calculated. -The choice was proposed by -`Power et al. (2015) `__, and it -represents approximately the mean distance from the cerebral cortex to the center of the head. +The choice was proposed by [Power2012]_, and it represents approximately the mean +distance from the cerebral cortex to the center of the head. """ @@ -38,8 +45,7 @@ def compute_fd_from_motion( The framewise displacement is the sum of the magnitudes of the translational and rotational motion, computed from the frame-to-frame differences along - the three spatial axes ( - `Power et al. (2015) `__). + the three spatial axes [Power2012]_. Each row in the motion parameters represents one frame, and columns represent each coordinate axis ``x``, `y``, and ``z``. Translation @@ -78,10 +84,15 @@ def compute_fd_from_transform( img: nb.spatialimages.SpatialImage, test_xfm: TransformBase, radius: float = DEFAULT_FD_RADIUS, + n_vertices: int = 8, ) -> float: """ Compute the framewise displacement (FD) for a given transformation. + This implementation varies with respect to the original formulation by [Power2012]_ + in that the FD is computed as the average across a number of vertices sample of the + sphere. + Parameters ---------- img : :obj:`~nibabel.spatialimages.SpatialImage` @@ -90,6 +101,8 @@ def compute_fd_from_transform( The transformation to test. Applied to coordinates around the image center. radius : :obj:`float`, optional The radius (in mm) of the spherical neighborhood around the center of the image. + n_vertices : :obj:`int`, optional + The number of vertices to sample on the sphere. Returns ------- @@ -103,7 +116,7 @@ def compute_fd_from_transform( # Convert to world coordinates center_xyz = nb.affines.apply_affine(affine, center_ijk) # Generate coordinates of points at radius distance from center - fd_coords = np.array(list(product(*((radius, -radius),) * 3))) + center_xyz + fd_coords = sample_unit_sphere(n_points=n_vertices) * radius + center_xyz # Compute the average displacement from the test transformation return np.mean(np.linalg.norm(test_xfm.map(fd_coords) - fd_coords, axis=-1)) @@ -165,13 +178,159 @@ def extract_motion_parameters(affine: np.ndarray) -> Tuple[np.ndarray, np.ndarra translation = affine[:3, 3] rotation_rad = np.arctan2( - [affine[2, 1], affine[0, 2], affine[1, 0]], [affine[2, 2], affine[0, 0], affine[1, 1]] + [affine[2, 1], affine[0, 2], affine[1, 0]], + [affine[2, 2], affine[0, 0], affine[1, 1]], ) rotation_deg = np.rad2deg(rotation_rad) return *translation, *rotation_deg -def identify_spikes(fd: np.ndarray, z_threshold: float = 2.0) -> Tuple[np.ndarray, np.ndarray]: +def sample_unit_sphere(n_points: int = 8) -> np.ndarray: + """Returns :math:`N` evenly distributed points on a unit-radius sphere. + + This function returns a **deterministic**, **quasi-uniform** point set on the + surface of the unit sphere :math:`S^2 \\subset \\mathbb{R}^3`. + + Notes + ----- + - There is no unique notion of "evenly distributed" for arbitrary `N` on a sphere. + This function uses: + * **Platonic solids** for certain small `N` (high symmetry; e.g. `N=6` gives + the ±axis points). + * A **Fibonacci / golden-angle spiral** otherwise (fast, simple, good coverage). + + Parameters + ---------- + n_points : int + Number of points on the sphere. + + Returns + ------- + numpy.ndarray + An array of shape ``(n_points, 3)`` whose rows have unit norm. + + Examples + -------- + Basic shape + unit norm: + + >>> import numpy as np + >>> X = sample_unit_sphere(10) + >>> X.shape + (10, 3) + >>> bool(np.allclose(np.linalg.norm(X, axis=1), 1.0)) + True + + Edge case: + + >>> sample_unit_sphere(0) + Traceback (most recent call last): + ... + ValueError: n_points must be a positive integer + + For N=6, return the ±axis points (octahedron vertices): + + >>> X = sample_unit_sphere(6) + >>> # Each row has exactly one coordinate with magnitude 1, others 0 + >>> bool(np.all((np.abs(X) == 1.0).sum(axis=1) == 1)) + True + >>> bool(np.all((np.abs(X) == 1.0).sum(axis=0) == 2)) # each axis appears twice (±) + True + + For N=4, the tetrahedron has constant pairwise dot product -1/3 off-diagonal: + + >>> X = sample_unit_sphere(4) + >>> D = X @ X.T + >>> off = D[~np.eye(4, dtype=bool)] + >>> bool(np.allclose(off, -1/3)) + True + + For a quasi-uniform set, the second moment matrix is close to I/3, and the + minimum angular separation is non-trivial: + + >>> X = sample_unit_sphere(200) + >>> M = (X.T @ X) / len(X) + >>> bool(np.allclose(M, np.eye(3) / 3, atol=1e-3)) + True + >>> def min_angle_rad(Y): + ... dots = np.clip(Y @ Y.T, -1.0, 1.0) + ... np.fill_diagonal(dots, 1.0) + ... ang = np.arccos(dots) + ... np.fill_diagonal(ang, np.inf) + ... return float(ang.min()) + >>> min_angle_rad(X) > 0.18 + True + + """ + if isinstance(n_points, (bool, np.bool_)) or not isinstance( + n_points, (int, np.integer) + ): + raise TypeError("n_points must be an integer") + if n_points < 1: + raise ValueError("n_points must be a positive integer") + + def _normalize(X): + X = np.asarray(X, dtype=float) + X /= np.linalg.norm(X, axis=1, keepdims=True) + return X + + # --- Highly symmetric small-N cases (Platonic solids / degeneracies) --- + if n_points == 1: + return np.array([[0.0, 0.0, 1.0]]) + if n_points == 2: + return np.array([[0.0, 0.0, 1.0], [0.0, 0.0, -1.0]]) + if n_points == 4: # tetrahedron (4 cube corners with even number of minus signs) + X = np.array( + [[1.0, 1.0, 1.0], [1.0, -1.0, -1.0], [-1.0, 1.0, -1.0], [-1.0, -1.0, 1.0]] + ) + return _normalize(X) + if n_points == 6: # octahedron: ±axes + return np.vstack([np.eye(3), -np.eye(3)]).astype(float) + if n_points in (8, 20): # hexahedron (cube) & base for dodecahedron + X = np.array(list(product([-1.0, 1.0], repeat=3)), dtype=float) + if n_points == 8: + return _normalize(X) + + # Dodecahedron, add 12 vertices + extra = [] + phi = (1.0 + math.sqrt(5.0)) / 2.0 + invphi = 1.0 / phi + for a in (-1.0, 1.0): + for b in (-1.0, 1.0): + extra.extend( + [ + [0.0, a * invphi, b * phi], + [a * invphi, b * phi, 0.0], + [a * phi, 0.0, b * invphi], + ] + ) + X = np.vstack([X, np.asarray(extra, dtype=float)]) + return _normalize(X) + + if n_points == 12: # icosahedron + phi = (1.0 + math.sqrt(5.0)) / 2.0 + X = [] + for a in (-1.0, 1.0): + for b in (-1.0, 1.0): + X.append([0.0, a, b * phi]) + X.append([a, b * phi, 0.0]) + X.append([a * phi, 0.0, b]) + return _normalize(X) + + # --- General N: Fibonacci / golden-angle spiral (deterministic quasi-uniform) --- + i = np.arange(n_points, dtype=float) + golden_angle = math.pi * (3.0 - math.sqrt(5.0)) + + z = 1.0 - 2.0 * (i + 0.5) / n_points + r = np.sqrt(np.maximum(0.0, 1.0 - z * z)) + theta = golden_angle * i + + X = np.column_stack((r * np.cos(theta), r * np.sin(theta), z)) + return _normalize(X) + + +def identify_spikes( + fd: np.ndarray, z_threshold: float = 2.0 +) -> Tuple[np.ndarray, np.ndarray]: """Identify motion spikes in framewise displacement data. Identifies high-motion frames as timepoint exceeding a given threshold value From e8b57834392142d2696bb8de67e755bd4c29fb18 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Sun, 4 Jan 2026 16:21:19 +0100 Subject: [PATCH 07/13] fix: defensive increase of tolerance in FD test --- nitransforms/tests/test_analysis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/tests/test_analysis.py b/nitransforms/tests/test_analysis.py index 3307b8b4..85a6f77c 100644 --- a/nitransforms/tests/test_analysis.py +++ b/nitransforms/tests/test_analysis.py @@ -90,7 +90,7 @@ def test_displacements_within_mask(simple_mask_img, test_xfm, reference_xfm, exp ) def test_compute_fd_from_transform(simple_mask_img, test_xfm, expected): fd = compute_fd_from_transform(simple_mask_img, test_xfm) - assert np.isclose(fd, expected) + assert np.isclose(fd, expected, atol=1e-4, rtol=1e-6) @pytest.mark.parametrize( From f9744b9482b9fbec443c761fbc90c8a4eeb547e0 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Sun, 4 Jan 2026 16:21:38 +0100 Subject: [PATCH 08/13] doc: add missing modules to apidoc --- docs/_api/analysis.rst | 6 ++++++ docs/_api/resampling.rst | 6 ++++++ docs/api.rst | 6 ++++-- 3 files changed, 16 insertions(+), 2 deletions(-) create mode 100644 docs/_api/analysis.rst create mode 100644 docs/_api/resampling.rst diff --git a/docs/_api/analysis.rst b/docs/_api/analysis.rst new file mode 100644 index 00000000..5d3b146d --- /dev/null +++ b/docs/_api/analysis.rst @@ -0,0 +1,6 @@ +======== +Analysis +======== + +.. automodule:: nitransforms.analysis.utils + :members: diff --git a/docs/_api/resampling.rst b/docs/_api/resampling.rst new file mode 100644 index 00000000..7488e3de --- /dev/null +++ b/docs/_api/resampling.rst @@ -0,0 +1,6 @@ +========== +Resampling +========== + +.. automodule:: nitransforms.resampling + :members: diff --git a/docs/api.rst b/docs/api.rst index a57d6836..4e384142 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -5,11 +5,13 @@ Information on specific functions, classes, and methods for developers. .. toctree:: :maxdepth: 1 + _api/analysis _api/base + _api/interp _api/io _api/linear _api/manip _api/nonlinear - _api/surface - _api/interp _api/patched + _api/resampling + _api/surface From c437e382e3d83061bf90a3f07df6a983461b52a9 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Sun, 4 Jan 2026 17:06:54 +0100 Subject: [PATCH 09/13] tst: exercise all branches --- nitransforms/analysis/utils.py | 37 +++++++++++++++++++++------------- 1 file changed, 23 insertions(+), 14 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 7d3c0ffa..e8bdb7d6 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -214,18 +214,15 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: Basic shape + unit norm: >>> import numpy as np - >>> X = sample_unit_sphere(10) - >>> X.shape - (10, 3) - >>> bool(np.allclose(np.linalg.norm(X, axis=1), 1.0)) - True - - Edge case: - - >>> sample_unit_sphere(0) - Traceback (most recent call last): - ... - ValueError: n_points must be a positive integer + >>> for n_points in (1, 2, 8, 10, 12, 20): + ... X = sample_unit_sphere(n_points) + ... X.shape, bool(np.allclose(np.linalg.norm(X, axis=1), 1.0)) + ((1, 3), True) + ((2, 3), True) + ((8, 3), True) + ((10, 3), True) + ((12, 3), True) + ((20, 3), True) For N=6, return the ±axis points (octahedron vertices): @@ -260,13 +257,25 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: >>> min_angle_rad(X) > 0.18 True + Improper inputs: + + >>> sample_unit_sphere(True) + Traceback (most recent call last): + ... + TypeError: n_points must be a positive integer + + >>> sample_unit_sphere(0) + Traceback (most recent call last): + ... + ValueError: n_points must be 1 or greater + """ if isinstance(n_points, (bool, np.bool_)) or not isinstance( n_points, (int, np.integer) ): - raise TypeError("n_points must be an integer") + raise TypeError("n_points must be a positive integer") if n_points < 1: - raise ValueError("n_points must be a positive integer") + raise ValueError("n_points must be 1 or greater") def _normalize(X): X = np.asarray(X, dtype=float) From bacdc83791c0237c90da7f3670e4815093d328a3 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Sun, 4 Jan 2026 17:19:49 +0100 Subject: [PATCH 10/13] fix: address bug in ``compute_fd_from_transform`` --- nitransforms/analysis/utils.py | 27 +++++++++++++++++---------- 1 file changed, 17 insertions(+), 10 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index e8bdb7d6..229b50e2 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -27,6 +27,7 @@ from scipy.stats import zscore from nitransforms.base import TransformBase +from nitransforms.linear import Affine DEFAULT_FD_RADIUS = 50.0 @@ -82,7 +83,8 @@ def compute_fd_from_motion( def compute_fd_from_transform( img: nb.spatialimages.SpatialImage, - test_xfm: TransformBase, + xfm: TransformBase, + xfm_prev: TransformBase | None = None, radius: float = DEFAULT_FD_RADIUS, n_vertices: int = 8, ) -> float: @@ -97,8 +99,11 @@ def compute_fd_from_transform( ---------- img : :obj:`~nibabel.spatialimages.SpatialImage` The reference image. Used to extract the center coordinates. - test_xfm : :obj:`~nitransforms.base.TransformBase` + xfm : :obj:`~nitransforms.base.TransformBase` The transformation to test. Applied to coordinates around the image center. + xfm_prev : :obj:`~nitransforms.base.TransformBase`, optional + A previous transformation to compare with. If ``None``, the identity + transformation is assumed (no transformation). radius : :obj:`float`, optional The radius (in mm) of the spherical neighborhood around the center of the image. n_vertices : :obj:`int`, optional @@ -110,6 +115,8 @@ def compute_fd_from_transform( The average framewise displacement (FD) for the test transformation. """ + xfm_prev = Affine() if xfm_prev is None else xfm_prev + affine = img.affine # Compute the center of the image in voxel space center_ijk = 0.5 * (np.array(img.shape[:3]) - 1) @@ -118,13 +125,13 @@ def compute_fd_from_transform( # Generate coordinates of points at radius distance from center fd_coords = sample_unit_sphere(n_points=n_vertices) * radius + center_xyz # Compute the average displacement from the test transformation - return np.mean(np.linalg.norm(test_xfm.map(fd_coords) - fd_coords, axis=-1)) + return np.mean(np.linalg.norm(xfm.map(fd_coords) - xfm_prev.map(fd_coords), axis=-1)) def displacements_within_mask( mask_img: nb.spatialimages.SpatialImage, - test_xfm: TransformBase, - reference_xfm: TransformBase | None = None, + xfm: TransformBase, + xfm_prev: TransformBase | None = None, ) -> np.ndarray: """ Compute the distance between voxel coordinates mapped through two transforms. @@ -134,11 +141,11 @@ def displacements_within_mask( mask_img : :obj:`~nibabel.spatialimages.SpatialImage` A mask image that defines the region of interest. Voxel coordinates within the mask are transformed. - test_xfm : :obj:`~nitransforms.base.TransformBase` + xfm : :obj:`~nitransforms.base.TransformBase` The transformation to test. This transformation is applied to the voxel coordinates. - reference_xfm : :obj:`~nitransforms.base.TransformBase`, optional - A reference transformation to compare with. If ``None``, the identity + xfm_prev : :obj:`~nitransforms.base.TransformBase`, optional + A previous (reference) transformation to compare with. If ``None``, the identity transformation is assumed (no transformation). Returns @@ -155,10 +162,10 @@ def displacements_within_mask( np.argwhere(maskdata), ) # Apply the test transformation - targets = test_xfm.map(xyz) + targets = xfm.map(xyz) # Compute the difference (displacement) between the test and reference transformations - diffs = targets - xyz if reference_xfm is None else targets - reference_xfm.map(xyz) + diffs = targets - xyz if xfm_prev is None else targets - xfm_prev.map(xyz) return np.linalg.norm(diffs, axis=-1) From 6799c9b6b450dfe28c1228dc3857825d2061bac1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sun, 4 Jan 2026 15:50:22 -0500 Subject: [PATCH 11/13] Apply suggestions from code review --- nitransforms/analysis/utils.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index 229b50e2..e2e12a68 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -92,7 +92,7 @@ def compute_fd_from_transform( Compute the framewise displacement (FD) for a given transformation. This implementation varies with respect to the original formulation by [Power2012]_ - in that the FD is computed as the average across a number of vertices sample of the + in that the FD is computed as the average across a number of vertices sampled over the sphere. Parameters @@ -200,9 +200,9 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: Notes ----- - - There is no unique notion of "evenly distributed" for arbitrary `N` on a sphere. + - There is no unique notion of "evenly distributed" for arbitrary ``N`` on a sphere. This function uses: - * **Platonic solids** for certain small `N` (high symmetry; e.g. `N=6` gives + * **Platonic solids** for certain small ``N`` (high symmetry; e.g. ``N=6`` gives the ±axis points). * A **Fibonacci / golden-angle spiral** otherwise (fast, simple, good coverage). @@ -231,7 +231,7 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: ((12, 3), True) ((20, 3), True) - For N=6, return the ±axis points (octahedron vertices): + For ``N=6``, return the ±axis points (octahedron vertices): >>> X = sample_unit_sphere(6) >>> # Each row has exactly one coordinate with magnitude 1, others 0 @@ -240,7 +240,7 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: >>> bool(np.all((np.abs(X) == 1.0).sum(axis=0) == 2)) # each axis appears twice (±) True - For N=4, the tetrahedron has constant pairwise dot product -1/3 off-diagonal: + For ``N=4``, the tetrahedron has constant pairwise dot product -1/3 off-diagonal: >>> X = sample_unit_sphere(4) >>> D = X @ X.T @@ -248,7 +248,7 @@ def sample_unit_sphere(n_points: int = 8) -> np.ndarray: >>> bool(np.allclose(off, -1/3)) True - For a quasi-uniform set, the second moment matrix is close to I/3, and the + For a quasi-uniform set, the second moment matrix is close to ``I/3``, and the minimum angular separation is non-trivial: >>> X = sample_unit_sphere(200) From bbf7bf9254122e3794dec2ea7f4d06439fbbb278 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=20Haitz=20Legarreta=20Gorro=C3=B1o?= Date: Sun, 4 Jan 2026 15:50:40 -0500 Subject: [PATCH 12/13] Apply suggestions from code review --- nitransforms/analysis/utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index e2e12a68..bca72536 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -93,7 +93,8 @@ def compute_fd_from_transform( This implementation varies with respect to the original formulation by [Power2012]_ in that the FD is computed as the average across a number of vertices sampled over the - sphere. + sphere. See :func:`~nitransforms.analysis.utils.sample_unit_sphere` for details + about the vertex sampling method. Parameters ---------- From dab55cf116718b569487a40e45ce0af2a4a30204 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 8 Jan 2026 08:24:24 +0100 Subject: [PATCH 13/13] Apply suggestions from code review Co-authored-by: Chris Markiewicz --- nitransforms/analysis/utils.py | 21 +++++++++------------ 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/nitransforms/analysis/utils.py b/nitransforms/analysis/utils.py index bca72536..a89b3bd3 100644 --- a/nitransforms/analysis/utils.py +++ b/nitransforms/analysis/utils.py @@ -67,18 +67,15 @@ def compute_fd_from_motion( """ translations = motion_parameters[:, :3] - rotations_deg = motion_parameters[:, 3:] - rotations_rad = np.deg2rad(rotations_deg) + rotations = np.deg2rad(motion_parameters[:, 3:]) + + displacements = np.hstack(( + np.diff(translations, axis=0, prepend=0), + np.diff(rotations * radius, axis=0, prepend=0) + )) - # Compute differences between consecutive frames - d_translations = np.vstack([np.zeros((1, 3)), np.diff(translations, axis=0)]) - d_rotations = np.vstack([np.zeros((1, 3)), np.diff(rotations_rad, axis=0)]) - - # Convert rotations from radians to displacement on a sphere - rotation_displacement = d_rotations * radius - - # Compute FD as sum of absolute differences - return np.sum(np.abs(d_translations) + np.abs(rotation_displacement), axis=1) + # FD is the L1 norm (sum of absolute values) + return np.linalg.norm(displacements, ord=1, axis=1) def compute_fd_from_transform( @@ -126,7 +123,7 @@ def compute_fd_from_transform( # Generate coordinates of points at radius distance from center fd_coords = sample_unit_sphere(n_points=n_vertices) * radius + center_xyz # Compute the average displacement from the test transformation - return np.mean(np.linalg.norm(xfm.map(fd_coords) - xfm_prev.map(fd_coords), axis=-1)) + return np.mean(np.linalg.norm(xfm.map(fd_coords) - xfm_prev.map(fd_coords), ord=1, axis=-1)) def displacements_within_mask(