ENH: Add gradient unit-sphere normalization function

Add gradient unit-sphere normalization function.

Author: jhlegarreta

src/nifreeze/data/dmri.py

@@ -51,6 +51,9 @@
 DEFAULT_LOWB_THRESHOLD = 50
 """The lower bound for the b-value so that the orientation is considered a DW volume."""
 
+DEFAULT_GRADIENT_EPS = 1e-8
+"""Epsilon value for b-vector normalization."""
+
 DEFAULT_HIGHB_THRESHOLD = 8000
 """A b-value cap for DWI data."""
 
@@ -180,7 +183,18 @@ def format_gradients(
         formatted = formatted.astype(int)
 
     # Transpose if column-major
-    return formatted.T if formatted.shape[0] == 4 and formatted.shape[1] != 4 else formatted
+    formatted = formatted.T if formatted.shape[0] == 4 and formatted.shape[1] != 4 else formatted
+
+    # Normalize b-vectors in-place
+    bvecs = formatted[:, :3]
+    norms = np.linalg.norm(bvecs, axis=1)
+    mask = norms > DEFAULT_GRADIENT_EPS
+    if np.any(mask):
+        formatted[mask, :3] = bvecs[mask] / norms[mask, None]  # Norm b-vectors
+        formatted[mask, 3] *= norms[mask]  # Scale b-values by norm
+        formatted[~mask, :] = 0.0  # Zero-out small b-vectors
+
+    return formatted
 
 
 def validate_gradients(
@@ -652,3 +666,81 @@ def transform_fsl_bvec(
     ijk2ijk_xfm = np.linalg.inv(imaffine) @ xfm @ imaffine
 
     return ijk2ijk_xfm[:3, :3] @ b_ijk[:3]
+
+
+def normalize_gradients(value: np.ndarray, eps: float = 1e-8, copy: bool = True) -> np.ndarray:
+    """Normalize b-vectors in arrays of common shapes.
+
+    Parameters
+    ----------
+    value : :obj:`~numpy.ndarray`
+        Input array with shape one of:
+        - (N, 3) : rows are b-vector components (e.g., [gx gy gz])
+        - (N, 4) : first 3 columns are b-vector components (e.g., [gx gy gz b])
+        - (3, N) : columns are b-vector components (e.g., [gx gy gz].T)
+        - (4, N) : first 3 rows are b-vector components (e.g., [gx gy gz b].T)
+        - (3,) or (1,3) or (3,1) : single b-vector
+        Columns are checked first to disambiguate Nx3/Nx4 cases.
+    eps : float, optional
+        Threshold below which a vector is considered zero and left unchanged.
+    copy : bool, optional
+        If ``True``, returns a new array; modify in-place otherwise.
+
+    Returns
+    -------
+    out : :obj:`~numpy.ndarray`
+        Array with the same shape as ``value`` with each 3-component b-vector
+        normalized.
+    """
+    arr = np.asarray(value, dtype=float)
+
+    # 1D single vector
+    if arr.ndim == 1:
+        if arr.size != 3:
+            raise ValueError(GRADIENT_NORMALIZATION_LENGTH_ERROR_MSG)
+        norm = np.linalg.norm(arr)
+        if norm > eps:
+            if copy:
+                return arr / norm
+            else:
+                # Perform in-place normalization on the array view
+                arr[:] = arr / norm
+                return arr
+        else:
+            return arr.copy() if copy else arr
+
+    if arr.ndim != 2:
+        raise ValueError(GRADIENT_NORMALIZATION_SHAPE_ERROR_MSG)
+
+    rows, cols = arr.shape
+
+    # Prepare output (copy or in-place)
+    normalized_arr = arr.copy() if copy else arr
+
+    # Determine where the 3-component vectors live and create a (N, 3) view
+    # Check columns first to make Nx3/Nx4 deterministic
+    if cols == 4:
+        # Nx4: first 3 columns are b-vectors components, last are b-values
+        vecs = normalized_arr[:, :3]  # shape (N, 3)
+    elif cols == 3:
+        # Nx3: rows are vectors
+        vecs = normalized_arr  # shape (N, 3)
+    elif rows == 4:
+        # 4xN: first 3 rows are b-vector components, last row are b-values
+        # Create a (N, 3) view by transposing first 3 rows
+        vecs = normalized_arr[:3, :].T  # shape (N, 3)
+    elif rows == 3:
+        # 3xN: columns are vectors: normalize per-column
+        vecs = normalized_arr.T  # shape (N, 3)
+    else:
+        raise ValueError(
+            GRADIENT_NORMALIZATION_UNRECOGNIZED_SHAPE_ERROR_MSG.format(shape=arr.shape)
+        )
+
+    # Normalize in-place on vecs (which is a view into output)
+    norms = np.linalg.norm(vecs, axis=1)
+    mask = norms > eps
+    if np.any(mask):
+        vecs[mask] = vecs[mask] / norms[mask, None]
+
+    return normalized_arr

test/test_data_dmri.py

@@ -24,6 +24,7 @@
 
 import re
 from pathlib import Path
+from string import Formatter
 
 import attrs
 import nibabel as nb
@@ -49,6 +50,28 @@
 from nifreeze.utils.ndimage import load_api
 
 
+def _template_has_field(template: str, field_name: str | None = None) -> bool:
+    """Return True if `template` contains a format field.
+    If `field_name` is provided, return True only if that named field appears.
+
+    This uses Formatter.parse() so it recognizes real format fields and
+    ignores literal substrings that merely look like "{shape}".
+    """
+    formatter = Formatter()
+    for _literal_text, field, _format_spec, _conversion in formatter.parse(template):
+        if field is None:
+            # no field in this segment
+            continue
+        # field can be '' (positional {}), 'shape', or complex like 'shape[0]' or 'obj.attr'
+        if field_name is None:
+            return True
+        # Compare the base name before any attribute/indexing syntax
+        base = field.split(".", 1)[0].split("[", 1)[0]
+        if base == field_name:
+            return True
+    return False
+
+
 def _dwi_data_to_nifti(
     dwi_dataobj,
     affine,