Support ODX format for fixels/sh

src/modelarrayio/cli/export_results.py

@@ -14,6 +14,7 @@
 from modelarrayio.cli.h5_to_cifti import h5_to_cifti
 from modelarrayio.cli.h5_to_mif import h5_to_mif
 from modelarrayio.cli.h5_to_nifti import h5_to_nifti
+from modelarrayio.cli.h5_to_odx import h5_to_odx
 from modelarrayio.cli.parser_utils import _is_file, add_log_level_arg
 from modelarrayio.utils.misc import detect_modality_from_path
 
@@ -87,7 +88,7 @@ def export_results(
                 f'The template file appears to be MIF/fixel ({template_path!r}). '
                 'For MIF/fixel results, supply --index-file and --directions-file.'
             )
-        modality = 'cifti'
+        modality = 'odx' if detected == 'odx' else 'cifti'
     else:
         raise ValueError(
             'Cannot determine modality. Provide --mask (NIfTI), '
@@ -128,6 +129,21 @@ def export_results(
         )
         return 0
 
+    if modality == 'odx':
+        if cohort_file is None and example_file is None:
+            raise ValueError('One of --cohort-file or --example-file is required for ODX results.')
+        if example_file is None:
+            logger.warning('No example ODX provided; using first source_file from cohort.')
+            example_file = pd.read_csv(cohort_file)['source_file'].iloc[0]
+        h5_to_odx(
+            example_odx=example_file,
+            in_file=in_file,
+            analysis_name=analysis_name,
+            compress=compress,
+            output_dir=output_path,
+        )
+        return 0
+
     # cifti
     if cohort_file is None and example_file is None:
         raise ValueError('One of --cohort-file or --example-file is required for CIFTI results.')

src/modelarrayio/cli/h5_to_odx.py

@@ -0,0 +1,70 @@
+"""Convert ModelArray HDF5 results back to an ODX file.
+
+The ODX analogue of :mod:`modelarrayio.cli.h5_to_mif`. Reads the
+``results/<analysis_name>/results_matrix`` written by ModelArray and paints each
+result metric (e.g. ``<term>.estimate``, ``<term>.p.value.fdr``) onto a template
+ODX's group-fixel geometry as a per-fixel (DPF) array, producing a single ODX
+that can be visualized (e.g. in trxviz). For each ``p.value`` metric a
+``1m.p.value`` (1 - p) companion is also written, matching ``h5_to_mif``.
+"""
+
+from __future__ import annotations
+
+import logging
+from pathlib import Path
+
+import h5py
+import numpy as np
+
+from modelarrayio.cli import utils as cli_utils
+from modelarrayio.utils.odx import write_odx_results
+
+logger = logging.getLogger(__name__)
+
+
+def h5_to_odx(example_odx, in_file, analysis_name, compress, output_dir):
+    """Write ModelArray results from an HDF5 file onto a template ODX.
+
+    Parameters
+    ----------
+    example_odx : path-like
+        Template ODX whose group-fixel geometry the results are painted onto
+        (any of the per-subject ODX used to build the HDF5, or the group ODX).
+    in_file : path-like
+        HDF5 file containing ``results/<analysis_name>/results_matrix``.
+    analysis_name : str
+        Name of the ModelArray results group inside the HDF5 file.
+    compress : bool
+        Accepted for signature parity with the other ``h5_to_*`` writers; the
+        ``.odx`` archive is already compressed.
+    output_dir : path-like
+        Directory where ``<analysis_name>.odx`` is written.
+
+    Returns
+    -------
+    status : :obj:`int`
+        0 if successful.
+    """
+    del compress  # ODX archives are self-compressed; kept for API parity.
+    output_path = Path(output_dir)
+    output_path.mkdir(parents=True, exist_ok=True)
+
+    with h5py.File(in_file, 'r') as h5_data:
+        results_matrix = h5_data[f'results/{analysis_name}/results_matrix']
+        results_names = cli_utils.read_result_names(
+            h5_data, analysis_name, results_matrix, logger=logger
+        )
+        results: dict[str, np.ndarray] = {}
+        for result_col, result_name in enumerate(results_names):
+            safe = cli_utils.sanitize_result_name(result_name)
+            column = np.asarray(results_matrix[result_col, :], dtype=np.float32)
+            results[f'{analysis_name}_{safe}'] = column
+            if 'p.value' in safe:
+                results[f'{analysis_name}_{safe.replace("p.value", "1m.p.value")}'] = (
+                    1.0 - column
+                ).astype(np.float32)
+
+    out_odx = output_path / f'{analysis_name}.odx'
+    write_odx_results(example_odx, results, out_odx)
+    logger.info('Wrote %d result DPF arrays to %s', len(results), out_odx)
+    return 0

src/modelarrayio/cli/odx_to_h5.py

@@ -0,0 +1,117 @@
+"""Convert per-subject ODX fixel data to an HDF5 (or TileDB) ModelArray file.
+
+The ODX analogue of :mod:`modelarrayio.cli.mif_to_h5`. Each cohort row points
+to one per-subject ODX (single-column DPF) produced by
+``odx combine --per-subject-odx DIR``; the shared group-fixel geometry is read
+directly from the first ODX, so no separate index/directions files are needed.
+"""
+
+from __future__ import annotations
+
+import logging
+from pathlib import Path
+
+import h5py
+
+from modelarrayio.cli import utils as cli_utils
+from modelarrayio.utils.odx import gather_fixels_from_odx, load_cohort_odx
+
+logger = logging.getLogger(__name__)
+
+
+def odx_to_h5(
+    cohort_long,
+    backend='hdf5',
+    output=Path('fixelarray.h5'),
+    storage_dtype='float32',
+    compression='gzip',
+    compression_level=4,
+    shuffle=True,
+    chunk_voxels=0,
+    target_chunk_mb=2.0,
+    workers=1,
+    s3_workers=1,
+    split_outputs=False,
+):
+    """Load per-subject ODX fixel data and write a ModelArray store.
+
+    Parameters mirror :func:`modelarrayio.cli.mif_to_h5.mif_to_h5` minus the
+    ``index_file``/``directions_file`` arguments — ODX carries its own geometry.
+
+    Returns
+    -------
+    status : :obj:`int`
+        0 if successful, non-zero if an expected output was not written.
+    """
+    if cohort_long.empty:
+        raise ValueError('Cohort file does not contain any ODX scalar entries.')
+
+    # Every template-space ODX shares the group-fixel geometry; read it once.
+    first_source = str(cohort_long.iloc[0]['source_file'])
+    fixel_table, voxel_table = gather_fixels_from_odx(first_source)
+
+    logger.info('Extracting ODX data...')
+    scalars, sources_lists = load_cohort_odx(cohort_long, s3_workers)
+    if not sources_lists:
+        raise ValueError('Unable to derive scalar sources from cohort file.')
+
+    scalar_names = list(sources_lists.keys())
+
+    if backend == 'hdf5':
+        if split_outputs:
+            outputs: list[Path] = []
+            for scalar_name in scalar_names:
+                scalar_output = cli_utils.prepare_output_parent(
+                    cli_utils.prefixed_output_path(output, scalar_name)
+                )
+                with h5py.File(scalar_output, 'w') as h5_file:
+                    cli_utils.write_table_dataset(h5_file, 'fixels', fixel_table)
+                    cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table)
+                    cli_utils.write_hdf5_scalar_matrices(
+                        h5_file,
+                        {scalar_name: scalars[scalar_name]},
+                        {scalar_name: sources_lists[scalar_name]},
+                        storage_dtype=storage_dtype,
+                        compression=compression,
+                        compression_level=compression_level,
+                        shuffle=shuffle,
+                        chunk_voxels=chunk_voxels,
+                        target_chunk_mb=target_chunk_mb,
+                    )
+                outputs.append(scalar_output)
+            return int(not all(path.exists() for path in outputs))
+
+        output = cli_utils.prepare_output_parent(output)
+        with h5py.File(output, 'w') as h5_file:
+            cli_utils.write_table_dataset(h5_file, 'fixels', fixel_table)
+            cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table)
+            cli_utils.write_hdf5_scalar_matrices(
+                h5_file,
+                scalars,
+                sources_lists,
+                storage_dtype=storage_dtype,
+                compression=compression,
+                compression_level=compression_level,
+                shuffle=shuffle,
+                chunk_voxels=chunk_voxels,
+                target_chunk_mb=target_chunk_mb,
+            )
+        return int(not output.exists())
+
+    # tiledb backend
+    for scalar_name in scalar_names:
+        scalar_output = (
+            cli_utils.prefixed_output_path(output, scalar_name) if split_outputs else output
+        )
+        cli_utils.write_tiledb_scalar_matrices(
+            scalar_output,
+            {scalar_name: scalars[scalar_name]},
+            {scalar_name: sources_lists[scalar_name]},
+            storage_dtype=storage_dtype,
+            compression=compression,
+            compression_level=compression_level,
+            shuffle=shuffle,
+            chunk_voxels=chunk_voxels,
+            target_chunk_mb=target_chunk_mb,
+        )
+    return 0

src/modelarrayio/cli/to_modelarray.py

@@ -11,6 +11,7 @@
 from modelarrayio.cli.cifti_to_h5 import cifti_to_h5
 from modelarrayio.cli.mif_to_h5 import mif_to_h5
 from modelarrayio.cli.nifti_to_h5 import nifti_to_h5
+from modelarrayio.cli.odx_to_h5 import odx_to_h5
 from modelarrayio.cli.parser_utils import _is_file, add_log_level_arg
 from modelarrayio.utils.misc import load_and_normalize_cohort
 
@@ -37,7 +38,7 @@ def to_modelarray(
 ):
     """Load neuroimaging data and write to an HDF5 or TileDB modelarray file.
 
-    The modality (NIfTI, CIFTI, or MIF/fixel) is autodetected from the source
+    The modality (NIfTI, CIFTI, MIF/fixel, or ODX) is autodetected from the source
     file extensions listed in the cohort file.
 
     Parameters
@@ -91,6 +92,10 @@ def to_modelarray(
             )
         return mif_to_h5(index_file=index_file, directions_file=directions_file, **common_kwargs)
 
+    if modality == 'odx':
+        # ODX carries its own group-fixel geometry, so no index/directions files.
+        return odx_to_h5(**common_kwargs)
+
     # cifti
     return cifti_to_h5(**common_kwargs)
 
@@ -106,7 +111,7 @@ def _parse_to_modelarray():
     parser = argparse.ArgumentParser(
         description=(
             'Convert neuroimaging data to a modelarray HDF5 file. '
-            'The modality (NIfTI, CIFTI, or MIF/fixel) is autodetected from '
+            'The modality (NIfTI, CIFTI, MIF/fixel, or ODX) is autodetected from '
             'the source file extensions in the cohort file.'
         ),
         formatter_class=argparse.ArgumentDefaultsHelpFormatter,

src/modelarrayio/utils/misc.py

@@ -14,7 +14,7 @@
 
 
 def detect_modality_from_path(path: str) -> str:
-    """Return ``'cifti'``, ``'mif'``, or ``'nifti'`` based on file extension.
+    """Return ``'cifti'``, ``'mif'``, ``'odx'``, or ``'nifti'`` based on file extension.
 
     Parameters
     ----------
@@ -36,11 +36,13 @@ def detect_modality_from_path(path: str) -> str:
         return 'cifti'
     if path.endswith(('.mif.gz', '.mif')):
         return 'mif'
+    if path.rstrip('/').endswith('.odx'):
+        return 'odx'
     if path.endswith(('.nii.gz', '.nii')):
         return 'nifti'
     raise ValueError(
         f'Cannot detect modality from file extension: {path!r}. '
-        'Expected .mif, .nii, .nii.gz, or a CIFTI compound extension '
+        'Expected .mif, .odx, .nii, .nii.gz, or a CIFTI compound extension '
         '(e.g. .dscalar.nii, .pscalar.nii).'
     )