Use MifImage class instead of mrconvert

.github/workflows/tox.yml

@@ -44,6 +44,12 @@ jobs:
         with:
           fetch-depth: 0
 
+      - name: Install mrtrix3
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y --no-install-recommends mrtrix3
+          mrconvert -version | head -1
+
       - name: Install the latest version of uv
         uses: astral-sh/setup-uv@v8.0.0
 

docs/examples/plot_fixel_workflow.py

@@ -6,7 +6,7 @@
 use the ``modelarrayio to-modelarray`` command to convert the MIF files to the HDF5 format
 (``.h5``) used by **ModelArray**,
 and ``modelarrayio export-results`` to export results back to MIF.
-This guide assumes **ModelArrayIO** and **MRtrix** are already installed.
+This guide assumes **ModelArrayIO** is already installed.
 """
 
 # %%
@@ -132,8 +132,8 @@
 #
 # .. warning::
 #
-#    **Existing files are not overwritten.**  ``modelarrayio export-results`` calls ``mrconvert`` without
-#    ``-force``, so any ``.mif`` file already present in ``--output-dir`` with the same name
+#    **Existing files are not overwritten.** With ``modelarrayio export-results``, any ``.mif`` file
+#    already present in ``--output-dir`` with the same name
 #    will be left unchanged.  If ``--output-dir`` itself already exists you will see a
 #    ``WARNING: Output directory exists`` message, but no files will be deleted.  To start
 #    fresh, manually remove the output directory before re-running ``modelarrayio export-results``.

docs/installation.rst

@@ -14,16 +14,6 @@ If you want to use the most up-to-date version, you can install from the ``main`
     pip install git+https://github.com/PennLINC/ModelArrayIO.git
 
 
-MRtrix (required for fixel ``.mif`` only)
------------------------------------------
-
-Fixel-wise ``.mif`` conversion tools use MRtrix ``mrconvert``.
-Install MRtrix from `MRtrix's webpage <https://www.mrtrix.org/download/>`_ if needed.
-Run ``mrview`` in the terminal to verify the installation.
-
-If your data are in NIfTI or CIFTI format only, you can skip this step.
-
-
 What Next?
 ----------
 

pyproject.toml

@@ -182,6 +182,7 @@ exclude_lines = [
 [tool.pytest.ini_options]
 markers = [
   "s3: tests that require network access to public S3 (deselect with '-m not s3')",
+  "downloaded_data: tests that download external fixture data before running",
 ]
 log_cli = true
 

src/modelarrayio/cli/h5_to_mif.py

@@ -6,10 +6,10 @@
 from pathlib import Path
 
 import h5py
-import nibabel as nb
+import numpy as np
 
 from modelarrayio.cli import utils as cli_utils
-from modelarrayio.utils.mif import mif_to_nifti2, nifti2_to_mif
+from modelarrayio.utils.mif import MifImage, mif_to_image
 
 logger = logging.getLogger(__name__)
 
@@ -24,10 +24,9 @@ def h5_to_mif(example_mif, in_file, analysis_name, compress, output_dir):
     named ``results/has_names``. This data can be of any type and does not need to contain
     more than a single row of data. Instead, its attributes are read to get column names
     for the data represented in ``results/results_matrix``.
-    The function takes the example mif file and converts it to Nifti2 to get a header.
-    Then each column in ``results/results_matrix`` is extracted to fill the data of a
-    new Nifti2 file that gets converted to mif and named according to the corresponding
-    item in ``results/has_names``.
+    The function takes the example mif file as a template header. Then each column in
+    ``results/results_matrix`` is extracted to fill the data of a new ``MifImage`` and
+    named according to the corresponding item in ``results/has_names``.
 
     Parameters
     ----------
@@ -47,9 +46,10 @@ def h5_to_mif(example_mif, in_file, analysis_name, compress, output_dir):
     -------
     None
     """
-    # Get a template nifti image.
-    nifti2_img, _ = mif_to_nifti2(example_mif)
+    # Use the example MIF as the template so layout and metadata stay native to MIF.
+    template_img, _ = mif_to_image(example_mif)
     output_path = Path(output_dir)
+    output_path.mkdir(parents=True, exist_ok=True)
     ext = '.mif.gz' if compress else '.mif'
     with h5py.File(in_file, 'r') as h5_data:
         results_matrix = h5_data[f'results/{analysis_name}/results_matrix']
@@ -60,22 +60,50 @@ def h5_to_mif(example_mif, in_file, analysis_name, compress, output_dir):
         for result_col, result_name in enumerate(results_names):
             valid_result_name = cli_utils.sanitize_result_name(result_name)
             out_mif = output_path / f'{analysis_name}_{valid_result_name}{ext}'
-            temp_nifti2 = nb.Nifti2Image(
-                results_matrix[result_col, :].reshape(-1, 1, 1),
-                nifti2_img.affine,
-                header=nifti2_img.header,
+            write_mif(
+                arr=np.asarray(results_matrix[result_col, :], dtype=np.float32),
+                template_img=template_img,
+                out_file=out_mif,
             )
-            nifti2_to_mif(temp_nifti2, out_mif)
 
             if 'p.value' not in valid_result_name:
                 continue
 
             valid_result_name_1mpvalue = valid_result_name.replace('p.value', '1m.p.value')
             out_mif_1mpvalue = output_path / f'{analysis_name}_{valid_result_name_1mpvalue}{ext}'
-            output_mifvalues_1mpvalue = 1 - results_matrix[result_col, :]
-            temp_nifti2_1mpvalue = nb.Nifti2Image(
-                output_mifvalues_1mpvalue.reshape(-1, 1, 1),
-                nifti2_img.affine,
-                header=nifti2_img.header,
+            output_mifvalues_1mpvalue = np.asarray(
+                1 - results_matrix[result_col, :],
+                dtype=np.float32,
             )
-            nifti2_to_mif(temp_nifti2_1mpvalue, out_mif_1mpvalue)
+            write_mif(
+                arr=output_mifvalues_1mpvalue,
+                template_img=template_img,
+                out_file=out_mif_1mpvalue,
+            )
+
+    return 0
+
+
+def write_mif(arr, template_img, out_file):
+    """Write array to MIF file.
+
+    Parameters
+    ----------
+    arr : :obj:`numpy.ndarray`
+        Array to write to file.
+    template_img : :obj:`MifImage`
+        Template MIF image.
+    out_file : :obj:`pathlib.Path`
+        Output file to write.
+        If it already exists, this function will raise a warning and not overwrite.
+    """
+    if out_file.exists():
+        logger.warning('Output file already exists. Not overwriting. %s', out_file)
+        return
+
+    result_data = arr.reshape(template_img.shape)
+    result_header = template_img.header.copy()
+    result_header.set_data_shape(result_data.shape)
+    result_header.set_data_dtype(result_data.dtype)
+    result_img = MifImage(result_data, template_img.affine, header=result_header)
+    result_img.to_filename(out_file)

src/modelarrayio/utils/_mif_format.py

@@ -0,0 +1,175 @@
+"""Low-level MIF format parsing helpers."""
+
+import sys
+
+import numpy as np
+
+
+def _readline(fileobj) -> bytes:
+    """Read one newline-terminated line from *fileobj* using only ``read(1)``.
+
+    This works with any object that implements ``read(n)``, including
+    nibabel's ``ImageOpener`` and gzip file objects that lack ``readline``.
+    """
+    buf = bytearray()
+    while True:
+        ch = fileobj.read(1)
+        if not ch:
+            break
+        buf.extend(ch if isinstance(ch, (bytes, bytearray)) else ch.encode('latin-1'))
+        if buf[-1:] == b'\n':
+            break
+    return bytes(buf)
+
+
+_MIF_DTYPE_MAP: dict[str, str] = {
+    'Int8': 'i1',
+    'UInt8': 'u1',
+    'Int16': 'i2',
+    'UInt16': 'u2',
+    'Int32': 'i4',
+    'UInt32': 'u4',
+    'Int64': 'i8',
+    'UInt64': 'u8',
+    'Float32': 'f4',
+    'Float64': 'f8',
+    'CFloat32': 'c8',
+    'CFloat64': 'c16',
+}
+
+_NUMPY_TO_MIF_BASE: dict[tuple[str, int], str] = {
+    ('i', 1): 'Int8',
+    ('u', 1): 'UInt8',
+    ('i', 2): 'Int16',
+    ('u', 2): 'UInt16',
+    ('i', 4): 'Int32',
+    ('u', 4): 'UInt32',
+    ('i', 8): 'Int64',
+    ('u', 8): 'UInt64',
+    ('f', 4): 'Float32',
+    ('f', 8): 'Float64',
+    ('c', 8): 'CFloat32',
+    ('c', 16): 'CFloat64',
+}
+
+
+def _mif_parse_dtype(dtype_str: str) -> np.dtype:
+    """Convert a MIF datatype string (e.g. ``'Float32LE'``) to a numpy dtype."""
+    dtype_str = dtype_str.strip()
+    if dtype_str.endswith('LE'):
+        endian, base = '<', dtype_str[:-2]
+    elif dtype_str.endswith('BE'):
+        endian, base = '>', dtype_str[:-2]
+    else:
+        endian = '<' if sys.byteorder == 'little' else '>'
+        base = dtype_str
+
+    if base not in _MIF_DTYPE_MAP:
+        raise ValueError(f'Unknown MIF datatype: {dtype_str!r}')
+
+    type_char = _MIF_DTYPE_MAP[base]
+    if type_char in ('i1', 'u1'):  # single-byte types have no endianness
+        return np.dtype(type_char)
+    return np.dtype(endian + type_char)
+
+
+def _mif_dtype_to_str(dtype: np.dtype) -> str:
+    """Convert a numpy dtype to a MIF datatype string."""
+    dtype = np.dtype(dtype)
+    base_name = _NUMPY_TO_MIF_BASE.get((dtype.kind, dtype.itemsize))
+    if base_name is None:
+        raise ValueError(f'Cannot represent numpy dtype {dtype!r} in MIF format')
+    if dtype.itemsize == 1:
+        return base_name
+
+    byte_order = dtype.byteorder
+    if byte_order == '=':
+        byte_order = '<' if sys.byteorder == 'little' else '>'
+    elif byte_order == '|':
+        return base_name
+    return base_name + ('LE' if byte_order == '<' else 'BE')
+
+
+def _mif_parse_layout(layout_str: str, ndim: int) -> list[int]:
+    """Parse a MIF layout string to a list of symbolic strides (1-indexed, signed).
+
+    For example ``'-0,-1,+2'`` becomes ``[-1, -2, 3]``.  The absolute value
+    encodes ordering (1 = fastest-varying axis) and the sign encodes direction.
+    """
+    strides = []
+    for token in layout_str.strip().split(','):
+        token = token.strip()
+        if token.startswith('+'):
+            sign, val = 1, int(token[1:])
+        elif token.startswith('-'):
+            sign, val = -1, int(token[1:])
+        else:
+            sign, val = 1, int(token)
+        strides.append(sign * (val + 1))  # convert 0-indexed to 1-indexed
+    if len(strides) != ndim:
+        raise ValueError(f'Layout has {len(strides)} axes but dim has {ndim}: {layout_str!r}')
+    return strides
+
+
+def _mif_layout_to_str(layout: list[int]) -> str:
+    """Convert symbolic strides list to a MIF layout string."""
+    tokens = []
+    for s in layout:
+        sign = '+' if s >= 0 else '-'
+        val = abs(s) - 1  # convert 1-indexed back to 0-indexed
+        tokens.append(f'{sign}{val}')
+    return ','.join(tokens)
+
+
+def _mif_apply_layout(raw_flat: np.ndarray, shape: tuple, layout: list[int]) -> np.ndarray:
+    """Reorder flat MIF disk data into a canonical (positive-stride) numpy array.
+
+    MIF stores data with the axis whose ``|layout[i]|`` equals 1 varying
+    fastest on disk, and axes with a negative ``layout[i]`` stored reversed.
+    This function returns the array in the logical (mrtrix-canonical) order:
+    output index ``(i, j, k, ...)`` corresponds to mrtrix image coordinate
+    ``(i, j, k, ...)``, matching what ``mrconvert <in> <out> -strides 1,2,3,...``
+    would produce and what MRtrix tools (``fixelcfestats``, ``mrstats``, ...)
+    see when they access the image via the ``Image`` API.
+    """
+    ndim = len(shape)
+    # Sort axes from fastest (|layout|=1) to slowest
+    order = sorted(range(ndim), key=lambda i: abs(layout[i]))
+    # Disk layout in C-order: [slowest, ..., fastest]
+    disk_axes = list(reversed(order))
+    disk_shape = tuple(shape[i] for i in disk_axes)
+
+    data = raw_flat.reshape(disk_shape)
+
+    # Transpose: output axis i came from disk position inv_perm[i]
+    inv_perm = [0] * ndim
+    for disk_pos, orig_axis in enumerate(disk_axes):
+        inv_perm[orig_axis] = disk_pos
+    data = data.transpose(inv_perm)
+
+    # Flip any axis with a negative symbolic stride so the result is in
+    # mrtrix-canonical positive-stride order.
+    slicer = tuple(slice(None, None, -1) if layout[i] < 0 else slice(None) for i in range(ndim))
+    data = data[slicer]
+
+    return np.ascontiguousarray(data)
+
+
+def _mif_apply_layout_for_write(data: np.ndarray, layout: list[int]) -> np.ndarray:
+    """Reorder a canonical numpy array into MIF disk layout for writing.
+
+    Inverse of :func:`_mif_apply_layout`: first flip each axis whose
+    symbolic stride is negative, then transpose to the disk axis order
+    (``[slowest, ..., fastest]`` in C-order).
+    """
+    ndim = len(data.shape)
+
+    # Flip axes that will be stored reversed on disk
+    slicer = tuple(slice(None, None, -1) if layout[i] < 0 else slice(None) for i in range(ndim))
+    data = data[slicer]
+
+    # Transpose to disk order: [slowest, ..., fastest] in C-order
+    order = sorted(range(ndim), key=lambda i: abs(layout[i]))
+    disk_axes = list(reversed(order))
+    data = data.transpose(disk_axes)
+    return np.ascontiguousarray(data)