analysis

lacuna.analysis

Analysis modules for lesion decoding.

This module provides the base infrastructure and specific analysis implementations for processing lesion data.

Functions:

Name	Description
list_analyses	List all available analysis classes discovered via auto-discovery.
get_analysis	Get an analysis class by name.

Classes:

Name	Description
BaseAnalysis	Abstract base class for all analyses.
FunctionalNetworkMapping	Functional lesion network mapping analysis.
ParcelAggregation	Aggregate voxel-wise results to parcels.
RegionalDamage	Regional damage analysis.
StructuralNetworkMapping	Structural lesion network mapping analysis.

Examples:

>>> from lacuna.analysis import list_analyses, get_analysis
>>> for name, cls in list_analyses():
...     print(f"{name}: {cls.batch_strategy}")
FunctionalNetworkMapping: vectorized
ParcelAggregation: parallel
RegionalDamage: parallel
StructuralNetworkMapping: parallel

>>> FNM = get_analysis("FunctionalNetworkMapping")
>>> analysis = FNM(connectome_name="GSP1000")

AcceleratedFunctionalNetworkMapping

Bases: BaseAnalysis

Accelerated functional LNM via matrix multiplication (m @ C).

Parameters:

Name	Type	Description	Default
matrix_path	str or Path	Path to the group-FC matrix TSV produced by lacuna parcellate --modality functional. Must be a square matrix with region labels as its row/column index.	required
parcel_names	list[str]	Registered parcellation name(s) matching the atlas used to build C. Exactly one atlas must be supplied — its labels must align with the TSV row/column labels.	None
lesion_weighting	('fractional', 'binary', 'voxel_count')	Scheme for the m row vector: fractional: each touched region gets 1 / n_regions_touched binary: 1 if the region is touched, else 0 voxel_count: fraction of the parcel's voxels covered by the lesion	"fractional"
verbose	bool	If True, print progress messages.	False
keep_intermediate	bool	If True, include the m weight vector in results for inspection.	False

Source code in src/lacuna/analysis/accelerated_functional_network_mapping.py

class AcceleratedFunctionalNetworkMapping(BaseAnalysis):
    """Accelerated functional LNM via matrix multiplication (``m @ C``).

    Parameters
    ----------
    matrix_path : str or Path
        Path to the group-FC matrix TSV produced by
        ``lacuna parcellate --modality functional``.
        Must be a square matrix with region labels as its row/column index.
    parcel_names : list[str]
        Registered parcellation name(s) matching the atlas used to build ``C``.
        Exactly one atlas must be supplied — its labels must align with the
        TSV row/column labels.
    lesion_weighting : {"fractional", "binary", "voxel_count"}, default="fractional"
        Scheme for the ``m`` row vector:

        - ``fractional``: each touched region gets ``1 / n_regions_touched``
        - ``binary``: 1 if the region is touched, else 0
        - ``voxel_count``: fraction of the parcel's voxels covered by the lesion
    verbose : bool, default=False
        If True, print progress messages.
    keep_intermediate : bool, default=False
        If True, include the ``m`` weight vector in results for inspection.
    """

    batch_strategy = "parallel"
    TARGET_SPACE = None
    TARGET_RESOLUTION = None

    def __init__(
        self,
        matrix_path: str | Path,
        parcel_names: list[str] | None = None,
        lesion_weighting: Literal["fractional", "binary", "voxel_count"] = "fractional",
        verbose: bool = False,
        keep_intermediate: bool = False,
    ):
        super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

        if lesion_weighting not in _LESION_WEIGHTINGS:
            raise ValueError(
                f"lesion_weighting must be one of {_LESION_WEIGHTINGS}, got "
                f"'{lesion_weighting}'"
            )

        self.matrix_path = Path(matrix_path)
        self.lesion_weighting = lesion_weighting

        if not parcel_names:
            raise ValueError(
                "parcel_names is required for AcceleratedFunctionalNetworkMapping "
                "(pass --parcel-atlases or --custom-parcellation)."
            )
        if isinstance(parcel_names, str):
            parcel_names = [parcel_names]
        if len(parcel_names) != 1:
            raise ValueError(
                "AcceleratedFunctionalNetworkMapping currently supports exactly one "
                f"parcellation (got {len(parcel_names)}). Run it per atlas."
            )
        registered = {p.name for p in list_parcellations()}
        if parcel_names[0] not in registered:
            raise KeyError(
                f"Parcellation '{parcel_names[0]}' is not registered. "
                f"Use list_parcellations() or register via --custom-parcellation."
            )
        self.parcel_names = parcel_names

        self.logger = ConsoleLogger(verbose=verbose, width=70)

        self._c_matrix: np.ndarray | None = None
        self._c_labels: list[str] | None = None
        self._c_metadata: dict | None = None

    def _get_parameters(self) -> dict:
        params = super()._get_parameters()
        params.update(
            {
                "matrix_path": str(self.matrix_path),
                "parcel_names": list(self.parcel_names),
                "lesion_weighting": self.lesion_weighting,
            }
        )
        return params

    def _load_c_matrix(self) -> None:
        if self._c_matrix is not None:
            return
        if not self.matrix_path.exists():
            raise FileNotFoundError(f"Matrix TSV not found: {self.matrix_path}")
        df = pd.read_csv(self.matrix_path, sep="\t", index_col=0)
        if df.shape[0] != df.shape[1]:
            raise ValueError(
                f"Expected a square matrix in {self.matrix_path}, got shape {df.shape}."
            )
        if list(df.index) != list(df.columns):
            raise ValueError(
                "Matrix TSV row/column labels disagree; the file does not look like "
                "a ConnectivityMatrix produced by 'lacuna parcellate'."
            )
        self._c_matrix = df.values.astype(np.float64)
        self._c_labels = [str(x) for x in df.index]

        sidecar = self.matrix_path.with_suffix(".json")
        if sidecar.exists():
            import json

            with open(sidecar) as f:
                self._c_metadata = json.load(f)
        else:
            self._c_metadata = {}

    def _validate_inputs(self, mask_data: SubjectData) -> None:
        self._load_c_matrix()
        if mask_data.space is None:
            raise ValueError(
                "AcceleratedFunctionalNetworkMapping requires input with a known "
                "coordinate space."
            )

    def _build_weight_vector(
        self,
        mask_array: np.ndarray,
        atlas_values: np.ndarray,
        region_ids: list[int],
    ) -> tuple[np.ndarray, dict[int, int]]:
        """Return (m, voxel_counts_per_region) for the chosen weighting scheme."""
        flat_mask = mask_array.astype(bool).ravel()
        flat_atlas = atlas_values.ravel()

        voxels_per_region: dict[int, int] = {}
        hit_per_region: dict[int, int] = {}

        for rid in region_ids:
            region_mask = flat_atlas == rid
            voxels_per_region[rid] = int(region_mask.sum())
            hit_per_region[rid] = int(np.logical_and(region_mask, flat_mask).sum())

        touched = [rid for rid, hits in hit_per_region.items() if hits > 0]
        n_touched = len(touched)

        m = np.zeros(len(region_ids), dtype=np.float64)
        if n_touched == 0:
            return m, voxels_per_region

        if self.lesion_weighting == "binary":
            for col, rid in enumerate(region_ids):
                m[col] = 1.0 if hit_per_region[rid] > 0 else 0.0
        elif self.lesion_weighting == "fractional":
            w = 1.0 / n_touched
            for col, rid in enumerate(region_ids):
                if hit_per_region[rid] > 0:
                    m[col] = w
        else:  # voxel_count
            for col, rid in enumerate(region_ids):
                nvox = voxels_per_region[rid]
                if nvox > 0 and hit_per_region[rid] > 0:
                    m[col] = hit_per_region[rid] / nvox
        return m, voxels_per_region

    def _align_atlas_and_labels(
        self, mask_data: SubjectData
    ) -> tuple[np.ndarray, list[int], list[str]]:
        """Load atlas, resample to mask grid, and align with C's label order."""
        atlas_name = self.parcel_names[0]
        parc = load_parcellation(atlas_name)
        atlas_img = parc.image

        ref = nib.Nifti1Image(
            np.zeros(mask_data.mask_img.shape, dtype=np.int16),
            mask_data.mask_img.affine,
        )
        resampled = resample_to_img(
            atlas_img,
            ref,
            interpolation="nearest",
            force_resample=True,
            copy_header=True,
        )
        atlas_values = np.round(resampled.get_fdata()).astype(np.int64)

        label_to_rid: dict[str, int] = {}
        for rid, label in parc.labels.items():
            label_to_rid[str(label)] = int(rid)

        assert self._c_labels is not None
        region_ids: list[int] = []
        ordered_labels: list[str] = []
        missing: list[str] = []
        for lab in self._c_labels:
            if lab in label_to_rid:
                region_ids.append(label_to_rid[lab])
                ordered_labels.append(lab)
            else:
                missing.append(lab)
        if missing:
            raise ValueError(
                f"Atlas '{atlas_name}' is missing {len(missing)} label(s) present in "
                f"the C matrix (first few: {missing[:3]}). The atlas and matrix must "
                "match."
            )
        return atlas_values, region_ids, ordered_labels

    def _run_analysis(self, mask_data: SubjectData) -> dict[str, DataContainer]:
        self._load_c_matrix()
        assert self._c_matrix is not None and self._c_labels is not None

        atlas_values, region_ids, ordered_labels = self._align_atlas_and_labels(mask_data)
        mask_array = mask_data.mask_img.get_fdata()

        m, voxels_per_region = self._build_weight_vector(mask_array, atlas_values, region_ids)

        afnmap = m @ self._c_matrix  # (N,)
        atlas_name = self.parcel_names[0]

        results: dict[str, DataContainer] = {}

        n_touched = int((m > 0).sum())
        meta_base = {
            "atlas": atlas_name,
            "matrix_path": str(self.matrix_path),
            "lesion_weighting": self.lesion_weighting,
            "n_regions": len(ordered_labels),
            "n_regions_touched": n_touched,
        }

        afnmap_data = {lab: float(afnmap[i]) for i, lab in enumerate(ordered_labels)}
        afnmap_parc = ParcelData(
            name="afnmap",
            data=afnmap_data,
            region_labels=ordered_labels,
            parcel_names=[atlas_name],
            aggregation_method=f"afnm_{self.lesion_weighting}",
            metadata=dict(meta_base, description="Accelerated functional network map (m @ C)"),
        )
        results[
            build_result_key(
                atlas=atlas_name,
                source="AcceleratedFunctionalNetworkMapping",
                desc="afnmap",
            )
        ] = afnmap_parc

        if self.keep_intermediate:
            m_data = {lab: float(m[i]) for i, lab in enumerate(ordered_labels)}
            results[
                build_result_key(
                    atlas=atlas_name,
                    source="AcceleratedFunctionalNetworkMapping",
                    desc="afnmweights",
                )
            ] = ParcelData(
                name="afnmweights",
                data=m_data,
                region_labels=ordered_labels,
                parcel_names=[atlas_name],
                aggregation_method=f"afnm_{self.lesion_weighting}",
                metadata=dict(
                    meta_base,
                    description="Lesion-to-parcel weights m used for m @ C",
                    voxels_per_region={str(k): v for k, v in voxels_per_region.items()},
                ),
            )

        return results

BaseAnalysis

Bases: ABC

Abstract base class for all analysis modules.

This class defines the contract that all analysis implementations must follow, enabling plug-and-play extensibility. Subclasses must implement two abstract methods: - _validate_inputs: Check that input data meets analysis requirements - _run_analysis: Perform the actual analysis computation

The public run() method orchestrates the workflow and cannot be overridden.

Attributes:

Name	Type	Description
TARGET_SPACE	str or None	Coordinate space where computations are performed (e.g., "MNI152NLin6Asym"). Can be defined as a class attribute (static) or instance attribute (dynamic). For connectome-based analyses, this is typically set in init based on the registered connectome's space. Set to None or "atlas" for analyses that adapt to input data spaces (e.g., ParcelAggregation).
TARGET_RESOLUTION	float or None	Resolution in mm where computations are performed (e.g., 1 or 2). Can be defined as a class attribute or instance attribute. Set to None for analyses that adapt to input data resolution.
batch_strategy	str	Preferred batch processing strategy ("parallel", or "vectorized"). Default is "parallel". Subclasses should override this if they benefit from a different strategy (e.g., vectorized for network mapping analyses).

Examples:

>>> class LesionVolume(BaseAnalysis):
...     # Declare computation space - masks will be transformed to this space
...     TARGET_SPACE = "MNI152NLin6Asym"
...     TARGET_RESOLUTION = 2
...     batch_strategy = "parallel"  # Process subjects in parallel
...
...     def __init__(self):
...         super().__init__()
...
...     def _validate_inputs(self, mask_data):
...         # Validation happens AFTER automatic transformation to TARGET_SPACE
...         space = mask_data.get_coordinate_space()
...         if space != self.TARGET_SPACE:
...             raise ValueError(f"Expected {self.TARGET_SPACE}, got {space}")
...
...     def _run_analysis(self, mask_data):
...         # Lesion is guaranteed to be in TARGET_SPACE @ TARGET_RESOLUTION
...         volume = mask_data.get_volume_mm3()
...         return {"volume": volume}
...
>>> analysis = LesionVolume()
>>> result = analysis.run(mask_data)
>>> print(result.results["LesionVolume"])
{"volume": 523.5}

Source code in src/lacuna/analysis/base.py

class BaseAnalysis(ABC):
    """
    Abstract base class for all analysis modules.

    This class defines the contract that all analysis implementations must follow,
    enabling plug-and-play extensibility. Subclasses must implement two abstract
    methods:
    - `_validate_inputs`: Check that input data meets analysis requirements
    - `_run_analysis`: Perform the actual analysis computation

    The public `run()` method orchestrates the workflow and cannot be overridden.

    Attributes
    ----------
    TARGET_SPACE : str or None
        Coordinate space where computations are performed (e.g., "MNI152NLin6Asym").
        Can be defined as a class attribute (static) or instance attribute (dynamic).
        For connectome-based analyses, this is typically set in __init__ based on
        the registered connectome's space. Set to None or "atlas" for analyses
        that adapt to input data spaces (e.g., ParcelAggregation).
    TARGET_RESOLUTION : float or None
        Resolution in mm where computations are performed (e.g., 1 or 2).
        Can be defined as a class attribute or instance attribute.
        Set to None for analyses that adapt to input data resolution.
    batch_strategy : str
        Preferred batch processing strategy ("parallel", or "vectorized").
        Default is "parallel". Subclasses should override this if they benefit from
        a different strategy (e.g., vectorized for network mapping analyses).

    Examples
    --------
    >>> class LesionVolume(BaseAnalysis):
    ...     # Declare computation space - masks will be transformed to this space
    ...     TARGET_SPACE = "MNI152NLin6Asym"
    ...     TARGET_RESOLUTION = 2
    ...     batch_strategy = "parallel"  # Process subjects in parallel
    ...
    ...     def __init__(self):
    ...         super().__init__()
    ...
    ...     def _validate_inputs(self, mask_data):
    ...         # Validation happens AFTER automatic transformation to TARGET_SPACE
    ...         space = mask_data.get_coordinate_space()
    ...         if space != self.TARGET_SPACE:
    ...             raise ValueError(f"Expected {self.TARGET_SPACE}, got {space}")
    ...
    ...     def _run_analysis(self, mask_data):
    ...         # Lesion is guaranteed to be in TARGET_SPACE @ TARGET_RESOLUTION
    ...         volume = mask_data.get_volume_mm3()
    ...         return {"volume": volume}
    ...
    >>> analysis = LesionVolume()
    >>> result = analysis.run(mask_data)
    >>> print(result.results["LesionVolume"])
    {"volume": 523.5}
    """

    #: Preferred batch processing strategy (default: parallel)
    batch_strategy: str = "parallel"

    def __init__(self, verbose: bool = False, keep_intermediate: bool = False) -> None:
        """
        Initialize the analysis module.

        Parameters
        ----------
        verbose : bool, default=False
            If True, print progress messages. If False, run silently.
        keep_intermediate : bool, default=False
            If True, include intermediate results (e.g., warped mask images)
            in the output. Useful for debugging and quality control.

        Notes
        -----
        Subclasses should override this to accept analysis-specific parameters
        and store them as instance attributes for provenance tracking.
        Always call super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)
        when overriding.
        """
        self.verbose = verbose
        self.keep_intermediate = keep_intermediate

    def __repr__(self) -> str:
        """
        Return detailed string representation of the analysis object.

        Returns
        -------
        str
            String in format "ClassName(param1=value1, param2=value2, ...)"

        Examples
        --------
        >>> analysis = FunctionalNetworkMapping(connectome_name='GSP1000')
        >>> repr(analysis)
        "FunctionalNetworkMapping(connectome_name='GSP1000', ...)"
        """
        params = self._get_parameters()
        class_name = self.__class__.__name__

        if not params:
            return f"{class_name}()"

        # Format parameters - truncate long values
        param_strs = []
        for key, value in params.items():
            if isinstance(value, str) and len(value) > 50:
                value_str = f"'{value[:47]}...'"
            elif isinstance(value, str):
                value_str = f"'{value}'"
            else:
                value_str = str(value)
            param_strs.append(f"{key}={value_str}")

        params_formatted = ", ".join(param_strs)
        return f"{class_name}({params_formatted})"

    def __str__(self) -> str:
        """
        Return user-friendly string representation of the analysis.

        Returns
        -------
        str
            Human-readable description of the analysis configuration.

        Examples
        --------
        >>> analysis = FunctionalNetworkMapping(method='pearson')
        >>> print(analysis)
        FunctionalNetworkMapping Analysis
        Configuration:
          - method: pearson
          - connectome_path: /path/to/connectome.h5
          - t_threshold: 9.0
        """
        class_name = self.__class__.__name__
        params = self._get_parameters()

        if not params:
            return f"{class_name} Analysis (no parameters)"

        lines = [f"{class_name} Analysis", "Configuration:"]
        for key, value in params.items():
            # Truncate long strings for readability
            if isinstance(value, str) and len(value) > 60:
                value_str = f"{value[:57]}..."
            else:
                value_str = str(value)
            lines.append(f"  - {key}: {value_str}")

        return "\n".join(lines)

    @final
    def run(self, mask_data: SubjectData) -> SubjectData:
        """
        Execute the analysis on a SubjectData object.

        This is the ONLY public method users should call. It orchestrates
        the complete analysis workflow:
        1. Transform to TARGET_SPACE if defined
        2. Validate inputs via _validate_inputs()
        3. Run analysis via _run_analysis()
        4. Namespace results under the analysis class name
        5. Create new SubjectData with updated results
        6. Record provenance
        7. Return new SubjectData instance

        The input SubjectData is never modified (immutability principle).

        Parameters
        ----------
        mask_data : SubjectData
            Input data containing lesion mask, metadata, and any prior results.

        Returns
        -------
        SubjectData
            A NEW SubjectData instance with analysis results added to the
            .results dictionary under a namespace key (the analysis class name).

        Raises
        ------
        ValueError
            If input validation fails (via _validate_inputs).
        RuntimeError
            If analysis computation fails (via _run_analysis).

        Notes
        -----
        This method is marked @final to prevent subclasses from overriding it.
        All customization must happen in the protected abstract methods.

        When keep_intermediate=True, the mask used for analysis is included
        in results as 'analysis_mask'. This VoxelMap includes metadata
        indicating whether transformation occurred (was_transformed=True/False).

        Examples
        --------
        >>> analysis = FunctionalNetworkMapping(connectome='GSP1000')
        >>> result = analysis.run(mask_data)
        >>> print(result.results['FunctionalNetworkMapping']['rmap'])
        """
        # Track original input space info for analyses that need to transform back
        original_space = mask_data.space
        original_resolution = mask_data.resolution

        # Step 1: Transform to target space if TARGET_SPACE is defined
        transformed_data = self._ensure_target_space(mask_data)

        # Store original input info in metadata for _run_analysis to access
        # This allows analyses to transform results back to input space if requested
        if transformed_data is not mask_data:
            updated_metadata = transformed_data.metadata.copy()
            updated_metadata["_original_input_space"] = original_space
            updated_metadata["_original_input_resolution"] = original_resolution
            transformed_data = SubjectData(
                mask_img=transformed_data.mask_img,
                space=transformed_data.space,
                resolution=transformed_data.resolution,
                metadata=updated_metadata,
                provenance=transformed_data.provenance,
                results=transformed_data.results,
            )

        # Step 2: Validate inputs
        self._validate_inputs(transformed_data)

        # Step 3: Run analysis computation
        analysis_results = self._run_analysis(transformed_data)

        # Step 4: Namespace results under class name
        results_dict = analysis_results

        # Add analysis mask to intermediates if keep_intermediate=True
        # Note: This is a fallback that stores the space-transformed mask.
        # Analyses that resample further (e.g., FNM to connectome grid) should
        # store their own 'analysis_mask' in _run_analysis with the actual
        # mask used for computation. That will override this default.
        if self.keep_intermediate and "analysis_mask" not in results_dict:
            from lacuna.core.data_types import VoxelMap

            was_transformed = (
                original_space != transformed_data.space
                or original_resolution != transformed_data.resolution
            )
            analysis_mask = VoxelMap(
                name="analysis_mask",
                data=transformed_data.mask_img,
                space=transformed_data.space,
                resolution=transformed_data.resolution,
                metadata={
                    "description": (
                        "Mask transformed to analysis target space"
                        if was_transformed
                        else "Mask used for analysis (no transformation needed)"
                    ),
                    "was_transformed": was_transformed,
                    "original_space": original_space,
                    "original_resolution": original_resolution,
                    "analysis_space": transformed_data.space,
                    "analysis_resolution": transformed_data.resolution,
                },
            )
            results_dict["analysis_mask"] = analysis_mask

        namespace_key = self.__class__.__name__
        updated_results = transformed_data.results.copy()
        updated_results[namespace_key] = results_dict

        # Step 4: Create new SubjectData with updated results
        # Create a new instance with updated results (manual approach for namespace overwriting)
        result_mask_data = SubjectData(
            mask_img=transformed_data.mask_img,
            space=transformed_data.space,
            resolution=transformed_data.resolution,
            metadata=transformed_data.metadata,
            provenance=transformed_data.provenance,
            results=updated_results,
        )

        # Step 5: Record provenance
        provenance_record = create_provenance_record(
            function=f"{self.__class__.__module__}.{self.__class__.__name__}",
            parameters=self._get_parameters(),
            version=self._get_version(),
        )
        result_mask_data = result_mask_data.add_provenance(provenance_record)

        return result_mask_data

    @abstractmethod
    def _validate_inputs(self, mask_data: SubjectData) -> None:
        """
        Validate that mask_data meets the requirements for this analysis.

        Parameters
        ----------
        mask_data : SubjectData
            Input data to validate.

        Raises
        ------
        ValueError
            If validation fails. Error message should clearly explain what
            requirement was not met and how to fix it.

        Notes
        -----
        Common validations include:
        - Checking coordinate space (e.g., must be in MNI152)
        - Verifying lesion mask is binary
        - Ensuring required metadata fields are present
        - Checking for prerequisite results from other analyses

        This method is called automatically by run() before _run_analysis().

        Examples
        --------
        >>> def _validate_inputs(self, mask_data: SubjectData) -> None:
        ...     if mask_data.get_coordinate_space() != "MNI152Nlin6Asym":
        ...         raise ValueError(
        ...             "ExampleAnalysis requires data in MNI152Nlin6Asym space. "
        ...         )
        ...
        ...     data = mask_data.mask_img.get_fdata()
        ...     if not np.all(np.isin(data, [0, 1])):
        ...         raise ValueError("Mask data must be binary (0s and 1s).")
        """
        pass

    @abstractmethod
    def _run_analysis(self, mask_data: SubjectData) -> list["DataContainer"]:
        """
        Perform the core analysis computation.

        Parameters
        ----------
        mask_data : SubjectData
            Validated input data.

        Returns
        -------
        list[DataContainer]
            Analysis results as a list of DataContainer objects. Each result
            represents a distinct output (voxel map, ROI data, matrix, etc.).

        Raises
        ------
        RuntimeError
            If analysis computation fails.

        Notes
        -----
        This method contains the scientific logic of your analysis.
        It is called automatically by run() after validation succeeds.

        The returned list will be automatically namespaced under
        self.__class__.__name__ in the output SubjectData.results attribute.

        Do NOT modify the input mask_data object. Extract what you need,
        perform computations, and return results as a list of DataContainer objects.

        Examples
        --------
        >>> from lacuna.core.data_types import VoxelMap, ScalarMetric
        >>> def _run_analysis(self, mask_data: SubjectData) -> list[DataContainer]:
        ...     mask_array = mask_data.mask_img.get_fdata()
        ...
        ...     # Create voxel map result
        ...     correlation_img = self._compute_correlation_map(mask_array)
        ...     voxel_result = VoxelMap(
        ...         name="rmap",
        ...         data=correlation_img,
        ...         output_space=self.computation_space,
        ...         lesion_space=mask_data.coordinate_space
        ...     )
        ...
        ...     # Create summary statistics result
        ...     summary_result = ScalarMetric(
        ...         name="summarystatistics",
        ...         data={"sum": float(np.sum(mask_array))}
        ...     )
        ...
        ...     return [voxel_result, summary_result]
        """
        pass

    def _get_parameters(self) -> dict[str, Any]:
        """
        Get analysis parameters for provenance tracking.

        Returns
        -------
        Dict[str, Any]
            Dictionary of parameter names and values.

        Notes
        -----
        Override this method if your analysis has parameters that should
        be recorded in provenance. The base implementation returns verbose.
        Subclasses should call super()._get_parameters() and merge with their
        own parameters.

        Examples
        --------
        >>> def _get_parameters(self) -> Dict[str, Any]:
        ...     params = super()._get_parameters()  # Get verbose
        ...     params.update({
        ...         'threshold': self.threshold,
        ...         'method': self.method,
        ...         'connectome': self.connectome
        ...     })
        ...     return params
        """
        return {"verbose": self.verbose}

    def _ensure_target_space(self, mask_data: SubjectData) -> SubjectData:
        """
        Automatically transform lesion data to TARGET_SPACE if defined.

        This method is called automatically by run() before validation and analysis.
        If TARGET_SPACE and TARGET_RESOLUTION are defined (as class or instance
        attributes), the lesion data will be transformed to that space.

        Special cases:
        - If TARGET_SPACE is None or "atlas", no transformation is performed
          (analysis adapts to input space)
        - If TARGET_RESOLUTION is None, current resolution is preserved

        Note: Instance attributes take precedence over class attributes, allowing
        analyses like FunctionalNetworkMapping and StructuralNetworkMapping to
        dynamically set TARGET_SPACE based on connectome metadata.

        Parameters
        ----------
        mask_data : SubjectData
            Input lesion data

        Returns
        -------
        SubjectData
            Transformed lesion data (or original if no transformation needed)
        """
        # Check if this analysis defines a target space
        # Use getattr(self, ...) to pick up instance attributes (e.g., from connectome)
        # as well as class attributes
        target_space = getattr(self, "TARGET_SPACE", None)
        target_resolution = getattr(self, "TARGET_RESOLUTION", None)

        # Skip transformation if no target space defined or if set to "atlas" (adaptive)
        if target_space is None or target_space == "atlas":
            return mask_data

        # Get current space
        current_space = mask_data.space
        current_resolution = mask_data.resolution

        if current_space is None:
            raise ValueError(
                f"{self.__class__.__name__} requires lesion data with 'space' metadata. "
                f"Expected space: {target_space}"
            )

        # Validate that resolution is present when space is specified
        # This prevents silently ignoring resolution mismatches
        if current_resolution is None:
            raise ValueError(
                f"{self.__class__.__name__} requires lesion data with 'resolution' metadata. "
                f"Resolution is required when space is specified. "
                f"Got space='{current_space}' but resolution=None"
            )

        # Import here to avoid circular imports
        from lacuna.core.spaces import REFERENCE_AFFINES, CoordinateSpace, spaces_are_equivalent
        from lacuna.spatial.transform import transform_mask_data

        # Check if transformation needed (use space equivalence, not raw string match)
        needs_space_transform = not spaces_are_equivalent(current_space, target_space)
        needs_resolution_change = (
            target_resolution is not None
            and current_resolution is not None
            and current_resolution != target_resolution
        )

        if not needs_space_transform and not needs_resolution_change:
            # Already in target space
            return mask_data

        # Determine target resolution (use current if not specified)
        final_resolution = (
            target_resolution if target_resolution is not None else current_resolution
        )

        # Create target space object
        target_space_obj = CoordinateSpace(
            identifier=target_space,
            resolution=final_resolution,
            reference_affine=REFERENCE_AFFINES.get(
                (target_space, final_resolution), mask_data.affine
            ),
        )

        # Transform (logging handled by transform_mask_data)
        return transform_mask_data(
            mask_data, target_space_obj, image_name="mask", verbose=self.verbose
        )

    def _get_version(self) -> str:
        """
        Get analysis version for provenance tracking.

        Returns
        -------
        str
            Version string (e.g., "0.1.0").

        Notes
        -----
        Returns the lacuna package version for consistent provenance tracking.
        Override this method if you need custom version information.

        Examples
        --------
        >>> def _get_version(self) -> str:
        ...     from .. import __version__
        ...     return __version__
        """
        from .. import __version__

        return __version__

    def _validate_and_transform_space(
        self,
        mask_data: SubjectData,
        required_space: str,
        required_resolution: float | None = None,
    ) -> SubjectData:
        """Validate coordinate space and auto-transform if needed.

        This helper method provides a standard pattern for analysis modules
        to validate spatial requirements and automatically transform data
        to the required space if needed.

        Parameters
        ----------
        mask_data : SubjectData
            Input lesion data
        required_space : str
            Required coordinate space identifier (e.g., 'MNI152NLin2009cAsym')
        required_resolution : float, optional
            Required resolution in mm. If None, any resolution accepted.

        Returns
        -------
        SubjectData
            Original data (if already in required space) or transformed data

        Raises
        ------
        ValueError
            If space cannot be determined or transformation not available

        Examples
        --------
        >>> def _validate_inputs(self, mask_data: SubjectData) -> None:
        ...     # Ensure data is in MNI152NLin2009cAsym space at 2mm
        ...     mask_data = self._validate_and_transform_space(
        ...         mask_data,
        ...         required_space='MNI152NLin2009cAsym',
        ...         required_resolution=2
        ...     )
        ...     return mask_data
        """
        # Get current space from metadata
        current_space = mask_data.space
        current_resolution = mask_data.resolution

        if current_space is None:
            raise ValueError(
                "Cannot determine coordinate space from lesion data. "
                "Ensure metadata contains 'space' key."
            )

        # Check if transformation needed
        needs_space_transform = current_space != required_space
        needs_resolution_change = (
            required_resolution is not None and current_resolution != required_resolution
        )

        if not needs_space_transform and not needs_resolution_change:
            # Already in required space - no transformation needed
            return mask_data

        # Import transformation utilities
        from lacuna.core.spaces import REFERENCE_AFFINES, CoordinateSpace
        from lacuna.spatial.transform import transform_mask_data

        # Create target space
        target_resolution = (
            required_resolution if required_resolution is not None else current_resolution
        )
        target_space = CoordinateSpace(
            identifier=required_space,
            resolution=target_resolution,
            reference_affine=REFERENCE_AFFINES.get(
                (required_space, target_resolution), mask_data.affine
            ),
        )

        # Transform data (logging handled by transform_mask_data)
        return transform_mask_data(mask_data, target_space, image_name="mask", verbose=self.verbose)

__init__(verbose=False, keep_intermediate=False)

Initialize the analysis module.

Parameters:

Name	Type	Description	Default
verbose	bool	If True, print progress messages. If False, run silently.	False
keep_intermediate	bool	If True, include intermediate results (e.g., warped mask images) in the output. Useful for debugging and quality control.	False

Notes

Subclasses should override this to accept analysis-specific parameters and store them as instance attributes for provenance tracking. Always call super().init(verbose=verbose, keep_intermediate=keep_intermediate) when overriding.

Source code in src/lacuna/analysis/base.py

def __init__(self, verbose: bool = False, keep_intermediate: bool = False) -> None:
    """
    Initialize the analysis module.

    Parameters
    ----------
    verbose : bool, default=False
        If True, print progress messages. If False, run silently.
    keep_intermediate : bool, default=False
        If True, include intermediate results (e.g., warped mask images)
        in the output. Useful for debugging and quality control.

    Notes
    -----
    Subclasses should override this to accept analysis-specific parameters
    and store them as instance attributes for provenance tracking.
    Always call super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)
    when overriding.
    """
    self.verbose = verbose
    self.keep_intermediate = keep_intermediate

__repr__()

Return detailed string representation of the analysis object.

Returns:

Type	Description
str	String in format "ClassName(param1=value1, param2=value2, ...)"

Examples:

>>> analysis = FunctionalNetworkMapping(connectome_name='GSP1000')
>>> repr(analysis)
"FunctionalNetworkMapping(connectome_name='GSP1000', ...)"

Source code in src/lacuna/analysis/base.py

def __repr__(self) -> str:
    """
    Return detailed string representation of the analysis object.

    Returns
    -------
    str
        String in format "ClassName(param1=value1, param2=value2, ...)"

    Examples
    --------
    >>> analysis = FunctionalNetworkMapping(connectome_name='GSP1000')
    >>> repr(analysis)
    "FunctionalNetworkMapping(connectome_name='GSP1000', ...)"
    """
    params = self._get_parameters()
    class_name = self.__class__.__name__

    if not params:
        return f"{class_name}()"

    # Format parameters - truncate long values
    param_strs = []
    for key, value in params.items():
        if isinstance(value, str) and len(value) > 50:
            value_str = f"'{value[:47]}...'"
        elif isinstance(value, str):
            value_str = f"'{value}'"
        else:
            value_str = str(value)
        param_strs.append(f"{key}={value_str}")

    params_formatted = ", ".join(param_strs)
    return f"{class_name}({params_formatted})"

__str__()

Return user-friendly string representation of the analysis.

Returns:

Type	Description
str	Human-readable description of the analysis configuration.

Examples:

>>> analysis = FunctionalNetworkMapping(method='pearson')
>>> print(analysis)
FunctionalNetworkMapping Analysis
Configuration:
  - method: pearson
  - connectome_path: /path/to/connectome.h5
  - t_threshold: 9.0

Source code in src/lacuna/analysis/base.py

def __str__(self) -> str:
    """
    Return user-friendly string representation of the analysis.

    Returns
    -------
    str
        Human-readable description of the analysis configuration.

    Examples
    --------
    >>> analysis = FunctionalNetworkMapping(method='pearson')
    >>> print(analysis)
    FunctionalNetworkMapping Analysis
    Configuration:
      - method: pearson
      - connectome_path: /path/to/connectome.h5
      - t_threshold: 9.0
    """
    class_name = self.__class__.__name__
    params = self._get_parameters()

    if not params:
        return f"{class_name} Analysis (no parameters)"

    lines = [f"{class_name} Analysis", "Configuration:"]
    for key, value in params.items():
        # Truncate long strings for readability
        if isinstance(value, str) and len(value) > 60:
            value_str = f"{value[:57]}..."
        else:
            value_str = str(value)
        lines.append(f"  - {key}: {value_str}")

    return "\n".join(lines)

run(mask_data)

Execute the analysis on a SubjectData object.

This is the ONLY public method users should call. It orchestrates the complete analysis workflow: 1. Transform to TARGET_SPACE if defined 2. Validate inputs via _validate_inputs() 3. Run analysis via _run_analysis() 4. Namespace results under the analysis class name 5. Create new SubjectData with updated results 6. Record provenance 7. Return new SubjectData instance

The input SubjectData is never modified (immutability principle).

Parameters:

Name	Type	Description	Default
mask_data	SubjectData	Input data containing lesion mask, metadata, and any prior results.	required

Returns:

Type	Description
SubjectData	A NEW SubjectData instance with analysis results added to the .results dictionary under a namespace key (the analysis class name).

Raises:

Type	Description
ValueError	If input validation fails (via _validate_inputs).
RuntimeError	If analysis computation fails (via _run_analysis).

Notes

This method is marked @final to prevent subclasses from overriding it. All customization must happen in the protected abstract methods.

When keep_intermediate=True, the mask used for analysis is included in results as 'analysis_mask'. This VoxelMap includes metadata indicating whether transformation occurred (was_transformed=True/False).

Examples:

>>> analysis = FunctionalNetworkMapping(connectome='GSP1000')
>>> result = analysis.run(mask_data)
>>> print(result.results['FunctionalNetworkMapping']['rmap'])

Source code in src/lacuna/analysis/base.py

@final
def run(self, mask_data: SubjectData) -> SubjectData:
    """
    Execute the analysis on a SubjectData object.

    This is the ONLY public method users should call. It orchestrates
    the complete analysis workflow:
    1. Transform to TARGET_SPACE if defined
    2. Validate inputs via _validate_inputs()
    3. Run analysis via _run_analysis()
    4. Namespace results under the analysis class name
    5. Create new SubjectData with updated results
    6. Record provenance
    7. Return new SubjectData instance

    The input SubjectData is never modified (immutability principle).

    Parameters
    ----------
    mask_data : SubjectData
        Input data containing lesion mask, metadata, and any prior results.

    Returns
    -------
    SubjectData
        A NEW SubjectData instance with analysis results added to the
        .results dictionary under a namespace key (the analysis class name).

    Raises
    ------
    ValueError
        If input validation fails (via _validate_inputs).
    RuntimeError
        If analysis computation fails (via _run_analysis).

    Notes
    -----
    This method is marked @final to prevent subclasses from overriding it.
    All customization must happen in the protected abstract methods.

    When keep_intermediate=True, the mask used for analysis is included
    in results as 'analysis_mask'. This VoxelMap includes metadata
    indicating whether transformation occurred (was_transformed=True/False).

    Examples
    --------
    >>> analysis = FunctionalNetworkMapping(connectome='GSP1000')
    >>> result = analysis.run(mask_data)
    >>> print(result.results['FunctionalNetworkMapping']['rmap'])
    """
    # Track original input space info for analyses that need to transform back
    original_space = mask_data.space
    original_resolution = mask_data.resolution

    # Step 1: Transform to target space if TARGET_SPACE is defined
    transformed_data = self._ensure_target_space(mask_data)

    # Store original input info in metadata for _run_analysis to access
    # This allows analyses to transform results back to input space if requested
    if transformed_data is not mask_data:
        updated_metadata = transformed_data.metadata.copy()
        updated_metadata["_original_input_space"] = original_space
        updated_metadata["_original_input_resolution"] = original_resolution
        transformed_data = SubjectData(
            mask_img=transformed_data.mask_img,
            space=transformed_data.space,
            resolution=transformed_data.resolution,
            metadata=updated_metadata,
            provenance=transformed_data.provenance,
            results=transformed_data.results,
        )

    # Step 2: Validate inputs
    self._validate_inputs(transformed_data)

    # Step 3: Run analysis computation
    analysis_results = self._run_analysis(transformed_data)

    # Step 4: Namespace results under class name
    results_dict = analysis_results

    # Add analysis mask to intermediates if keep_intermediate=True
    # Note: This is a fallback that stores the space-transformed mask.
    # Analyses that resample further (e.g., FNM to connectome grid) should
    # store their own 'analysis_mask' in _run_analysis with the actual
    # mask used for computation. That will override this default.
    if self.keep_intermediate and "analysis_mask" not in results_dict:
        from lacuna.core.data_types import VoxelMap

        was_transformed = (
            original_space != transformed_data.space
            or original_resolution != transformed_data.resolution
        )
        analysis_mask = VoxelMap(
            name="analysis_mask",
            data=transformed_data.mask_img,
            space=transformed_data.space,
            resolution=transformed_data.resolution,
            metadata={
                "description": (
                    "Mask transformed to analysis target space"
                    if was_transformed
                    else "Mask used for analysis (no transformation needed)"
                ),
                "was_transformed": was_transformed,
                "original_space": original_space,
                "original_resolution": original_resolution,
                "analysis_space": transformed_data.space,
                "analysis_resolution": transformed_data.resolution,
            },
        )
        results_dict["analysis_mask"] = analysis_mask

    namespace_key = self.__class__.__name__
    updated_results = transformed_data.results.copy()
    updated_results[namespace_key] = results_dict

    # Step 4: Create new SubjectData with updated results
    # Create a new instance with updated results (manual approach for namespace overwriting)
    result_mask_data = SubjectData(
        mask_img=transformed_data.mask_img,
        space=transformed_data.space,
        resolution=transformed_data.resolution,
        metadata=transformed_data.metadata,
        provenance=transformed_data.provenance,
        results=updated_results,
    )

    # Step 5: Record provenance
    provenance_record = create_provenance_record(
        function=f"{self.__class__.__module__}.{self.__class__.__name__}",
        parameters=self._get_parameters(),
        version=self._get_version(),
    )
    result_mask_data = result_mask_data.add_provenance(provenance_record)

    return result_mask_data

FunctionalNetworkMapping

Bases: BaseAnalysis

Functional connectivity-based lesion network mapping.

This analysis maps functional connectivity disruption patterns by correlating a lesion's timeseries with whole-brain connectome data. Requires a pre-computed functional connectome in HDF5 format.

Computation space is determined dynamically from the registered connectome's metadata (space and resolution). Input masks are automatically transformed to match the connectome space before analysis.

Memory-efficient batch processing: Connectomes can be stored as single HDF5 files or directories with multiple batched files. All batches are processed sequentially to minimize memory usage.

Parameters:

Name	Type	Description	Default
connectome_name	str	Name of registered functional connectome (e.g., "GSP1000"). Use list_functional_connectomes() to see available connectomes. The connectome must be pre-registered via register_functional_connectome(). Each HDF5 file must contain: - 'timeseries': (n_subjects, n_timepoints, n_voxels) array - 'mask_indices': (3, n_voxels) or (n_voxels, 3) brain mask coordinates - 'mask_affine': (4, 4) affine transformation matrix - 'mask_shape': Tuple stored in attributes	required
method	(boes, pini)	Timeseries extraction method: - "boes": Mean timeseries across all lesion voxels - "pini": PCA-based selection of representative voxels	"boes"
pini_percentile	int	For PINI method: percentile threshold for PC1 loadings (0-100). Higher values select fewer, more representative voxels.	20
n_jobs	int	Number of parallel jobs for batch processing (not yet implemented).	1

Attributes:

Name	Type	Description
batch_strategy	str	Set to "vectorized" for optimized batch processing.

Methods:

Name	Description
run	Inherited from BaseAnalysis. Computes functional network mapping.

Examples:

>>> from lacuna import SubjectData
>>> from lacuna.analysis import FunctionalNetworkMapping
>>> from lacuna.assets.connectomes import (
...     list_functional_connectomes,
...     register_functional_connectome,
... )
>>>
>>> # Register a connectome (do this once)
>>> register_functional_connectome(
...     name="GSP1000",
...     space="MNI152NLin6Asym",
...     resolution=2.0,
...     data_path="/data/gsp1000_connectome.h5",
...     n_subjects=1000,
...     description="GSP1000 voxel-wise connectome"
... )
>>>
>>> # List available connectomes
>>> list_functional_connectomes()
>>>
>>> # Use registered connectome
>>> lesion = SubjectData.from_nifti("lesion_mni.nii.gz")
>>> analysis = FunctionalNetworkMapping(
...     connectome_name="GSP1000",
...     method="boes"
... )
>>> result = analysis.run(lesion)
>>> correlation_map = result.results["FunctionalNetworkMapping"]["rmap"]
>>> z_map = result.results["FunctionalNetworkMapping"]["zmap"]

Notes

The connectome HDF5 file(s) must follow this structure: - Datasets: timeseries, mask_indices, mask_affine - Attributes: mask_shape - All voxels should be in MNI152 space

When using multiple batch files, all files must have the same mask structure (mask_indices, mask_affine, mask_shape).

References

• Boes et al. (2015): https://doi.org/10.1093/brain/awv228
• Pini et al. (2020): https://doi.org/10.1093/braincomms/fcab259

Source code in src/lacuna/analysis/functional_network_mapping.py

class FunctionalNetworkMapping(BaseAnalysis):
    """Functional connectivity-based lesion network mapping.

    This analysis maps functional connectivity disruption patterns by
    correlating a lesion's timeseries with whole-brain connectome data.
    Requires a pre-computed functional connectome in HDF5 format.

    Computation space is determined dynamically from the registered connectome's
    metadata (space and resolution). Input masks are automatically transformed to
    match the connectome space before analysis.

    Memory-efficient batch processing: Connectomes can be stored as single
    HDF5 files or directories with multiple batched files. All batches are
    processed sequentially to minimize memory usage.

    Parameters
    ----------
    connectome_name : str
        Name of registered functional connectome (e.g., "GSP1000").
        Use list_functional_connectomes() to see available connectomes.
        The connectome must be pre-registered via register_functional_connectome().
        Each HDF5 file must contain:
        - 'timeseries': (n_subjects, n_timepoints, n_voxels) array
        - 'mask_indices': (3, n_voxels) or (n_voxels, 3) brain mask coordinates
        - 'mask_affine': (4, 4) affine transformation matrix
        - 'mask_shape': Tuple stored in attributes
    method : {"boes", "pini"}, default="boes"
        Timeseries extraction method:
        - "boes": Mean timeseries across all lesion voxels
        - "pini": PCA-based selection of representative voxels
    pini_percentile : int, default=20
        For PINI method: percentile threshold for PC1 loadings (0-100).
        Higher values select fewer, more representative voxels.
    n_jobs : int, default=1
        Number of parallel jobs for batch processing (not yet implemented).

    Attributes
    ----------
    batch_strategy : str
        Set to "vectorized" for optimized batch processing.

    Methods
    -------
    run(mask_data: SubjectData) -> SubjectData
        Inherited from BaseAnalysis. Computes functional network mapping.

    Examples
    --------
    >>> from lacuna import SubjectData
    >>> from lacuna.analysis import FunctionalNetworkMapping
    >>> from lacuna.assets.connectomes import (
    ...     list_functional_connectomes,
    ...     register_functional_connectome,
    ... )
    >>>
    >>> # Register a connectome (do this once)
    >>> register_functional_connectome(
    ...     name="GSP1000",
    ...     space="MNI152NLin6Asym",
    ...     resolution=2.0,
    ...     data_path="/data/gsp1000_connectome.h5",
    ...     n_subjects=1000,
    ...     description="GSP1000 voxel-wise connectome"
    ... )
    >>>
    >>> # List available connectomes
    >>> list_functional_connectomes()
    >>>
    >>> # Use registered connectome
    >>> lesion = SubjectData.from_nifti("lesion_mni.nii.gz")
    >>> analysis = FunctionalNetworkMapping(
    ...     connectome_name="GSP1000",
    ...     method="boes"
    ... )
    >>> result = analysis.run(lesion)
    >>> correlation_map = result.results["FunctionalNetworkMapping"]["rmap"]
    >>> z_map = result.results["FunctionalNetworkMapping"]["zmap"]

    Notes
    -----
    The connectome HDF5 file(s) must follow this structure:
    - Datasets: timeseries, mask_indices, mask_affine
    - Attributes: mask_shape
    - All voxels should be in MNI152 space

    When using multiple batch files, all files must have the same mask
    structure (mask_indices, mask_affine, mask_shape).

    References
    ----------
    - Boes et al. (2015): https://doi.org/10.1093/brain/awv228
    - Pini et al. (2020): https://doi.org/10.1093/braincomms/fcab259
    """

    # Class attribute for batch processing strategy
    batch_strategy = "vectorized"

    def __init__(
        self,
        connectome_name: str,
        method: str = "boes",
        pini_percentile: int = 20,
        n_jobs: int = 1,
        verbose: bool = False,
        compute_p_map: bool = True,
        fdr_alpha: float | None = 0.05,
        t_threshold: float | None = None,
        return_in_input_space: bool = True,
        output_resolution: int | None = None,
        keep_intermediate: bool = False,
    ):
        """Initialize functional network mapping analysis.

        Parameters
        ----------
        connectome_name : str
            Name of registered functional connectome (e.g., "GSP1000").
            Use list_functional_connectomes() to see available options.
        method : {"boes", "pini"}, default="boes"
            Timeseries extraction method.
        pini_percentile : int, default=20
            Percentile threshold for PINI method (0-100).
        n_jobs : int, default=1
            Number of parallel jobs for post-processing (result aggregation
            and spatial resampling). Set to -1 to use all available CPUs.
        verbose : bool, default=False
            If True, print progress messages. If False, run silently.
        compute_p_map : bool, default=True
            If True, compute p-value map (two-tailed) from t-statistics.
        fdr_alpha : float, optional, default=0.05
            If provided, compute FDR-corrected p-value map using Benjamini-Hochberg
            procedure at the specified alpha level. Set to None to disable FDR correction.
            Requires compute_p_map=True.
        t_threshold : float, optional
            If provided, create binary mask of voxels with |t| > threshold.
        return_in_input_space : bool, default=True
            If True, transform VoxelMap outputs back to the original input mask space.
            If False, outputs remain in the connectome space (e.g., MNI152NLin6Asym).
            Requires input SubjectData to have valid space metadata.
        output_resolution : int, optional
            Final output resolution in mm (1 or 2). Controls the resolution of VoxelMap outputs.
            If None (default), matches the input mask resolution when return_in_input_space=True,
            or uses analysis resolution when return_in_input_space=False.
            Set explicitly to ensure consistent output resolution across analyses.
        keep_intermediate : bool, default=False
            If True, include intermediate results (e.g., warped mask images)
            in the output. Useful for debugging and quality control.

        Raises
        ------
        ValueError
            If method is not 'boes' or 'pini'.
        KeyError
            If connectome_name not found in registry.
        """
        super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

        # Validate method parameter
        if method not in ("boes", "pini"):
            msg = f"method must be 'boes' or 'pini', got '{method}'"
            raise ValueError(msg)

        # Load connectome from registry
        try:
            connectome = load_functional_connectome(connectome_name)
        except KeyError as e:
            available = [c.name for c in list_functional_connectomes()]
            raise KeyError(
                f"Connectome '{connectome_name}' not found in registry. "
                f"Available connectomes: {', '.join(available)}. "
                f"Use register_functional_connectome() to add new connectomes."
            ) from e

        # Store connectome information
        self.connectome_name = connectome_name
        self.connectome_path = connectome.data_path
        self.output_space = connectome.metadata.space
        self.output_resolution = connectome.metadata.resolution
        self._is_batch_dir = connectome.is_batched

        # Set TARGET_SPACE dynamically from connectome (used by BaseAnalysis._ensure_target_space)
        self.TARGET_SPACE = connectome.metadata.space
        self.TARGET_RESOLUTION = connectome.metadata.resolution

        # Analysis parameters
        self.method = method
        self.pini_percentile = pini_percentile
        if not (1 <= pini_percentile <= 100):
            raise ValueError(f"pini_percentile must be between 1 and 100, got {pini_percentile}")
        self.n_jobs = n_jobs
        self.compute_p_map = compute_p_map
        self.fdr_alpha = fdr_alpha
        if fdr_alpha is not None and not (0 < fdr_alpha <= 1):
            raise ValueError(
                f"fdr_alpha must be between 0 (exclusive) and 1 (inclusive), got {fdr_alpha}"
            )
        self.t_threshold = t_threshold
        self.return_in_input_space = return_in_input_space
        self.final_output_resolution = output_resolution  # User-specified, None means auto

        # Initialize logger
        self.logger = ConsoleLogger(verbose=verbose, width=70)

        # Internal state
        self._batch_files = None
        self._mask_info = None

    def _format_subject_id(self, mask_data: SubjectData) -> str:
        """Format a human-readable identifier for a subject.

        Combines subject_id, session_id, and label into a compact string.

        Parameters
        ----------
        mask_data : SubjectData
            Subject data with metadata.

        Returns
        -------
        str
            Formatted identifier like 'sub-001/ses-01/lesion'
        """
        parts = []
        metadata = mask_data.metadata

        subject_id = metadata.get("subject_id", "unknown")
        parts.append(subject_id)

        session_id = metadata.get("session_id")
        if session_id:
            parts.append(session_id)

        label = metadata.get("label")
        if label:
            parts.append(label)

        return "/".join(parts)

    def _get_connectome_files(self) -> list[Path]:
        """Get list of HDF5 connectome files to process.

        Returns
        -------
        list[Path]
            List of HDF5 file paths, sorted alphabetically.

        Raises
        ------
        ValidationError
            If no HDF5 files found.
        """
        if self._batch_files is not None:
            return self._batch_files

        try:
            self._batch_files = list_connectome_batch_files(self.connectome_path)
        except FileNotFoundError as e:
            raise ValidationError(str(e)) from e

        return self._batch_files

    def _validate_inputs(self, mask_data: SubjectData) -> None:
        """Validate inputs for functional network mapping.

        This method validates that the mask data is ready for FNM analysis.
        By the time this is called, BaseAnalysis.run() has already transformed
        the mask to TARGET_SPACE (the connectome space) via _ensure_target_space().

        Parameters
        ----------
        mask_data : SubjectData
            Mask data to validate (already transformed to connectome space).

        Raises
        ------
        ValidationError
            If connectome file(s) don't exist or mask space doesn't match
            the expected connectome space.

        Notes
        -----
        Binary mask validation is handled by SubjectData.__init__, so we don't
        need to duplicate that check here.

        Space transformation is handled by BaseAnalysis._ensure_target_space(),
        so by the time we get here, mask_data.space should equal self.TARGET_SPACE.
        """
        # Check connectome path exists
        if not self.connectome_path.exists():
            msg = f"Connectome path not found: {self.connectome_path}"
            raise ValidationError(msg)

        # Validate that we have connectome files
        _ = self._get_connectome_files()  # Raises ValidationError if no files found

        # Validate coordinate space matches connectome space
        # (should already be transformed by _ensure_target_space)
        if mask_data.space != self.TARGET_SPACE:
            msg = (
                f"Mask space '{mask_data.space}' does not match connectome space "
                f"'{self.TARGET_SPACE}'. This is unexpected - space transformation "
                f"should have been handled by BaseAnalysis._ensure_target_space()."
            )
            raise ValidationError(msg)

    def _run_analysis(self, mask_data: SubjectData) -> dict[str, "AnalysisResult"]:
        """Execute functional network mapping analysis.

        Processes connectome batches sequentially to minimize memory usage.
        Accumulates z-transformed correlation maps across all batches and
        performs final aggregation.

        Parameters
        ----------
        mask_data : SubjectData
            Validated lesion data in MNI152 space.

        Returns
        -------
        dict[str, AnalysisResult]
            Dictionary containing:
            - 'correlation_map': VoxelMapResult for correlation (r values)
            - 'z_map': VoxelMapResult for Fisher z-transformed
            - 't_map': VoxelMapResult for t-statistics
            - 't_threshold_map': VoxelMapResult (if t_threshold provided)
            - 'summary_statistics': MiscResult for summary statistics

        Notes
        -----
        This method implements a memory-efficient FNM pipeline:
        1. Load mask info once (shared across batches)
        2. For each connectome batch:
           a. Load timeseries data
           b. Extract lesion timeseries
           c. Compute correlation maps
           d. Fisher z-transform and accumulate
           e. Free memory
        3. Aggregate statistics across all batches
        4. Convert to 3D volumes and create NIfTI images
        """
        # Load mask information once
        if self._mask_info is None:
            self.logger.info("Loading mask information from connectome...")
            self._load_mask_info()
            mask_shape = self._mask_info["mask_shape"]
            n_voxels = len(self._mask_info["mask_indices"][0])
            self.logger.success(
                "Mask loaded", details={"shape": str(mask_shape), "n_voxels": n_voxels}
            )

        # Empty masks: produce zero-valued output maps
        if mask_data.is_empty_mask:
            self.logger.warning("Empty mask detected — producing zero-valued network maps")
            return self._build_empty_mask_results()

        # Get mask voxel indices and resampled mask (computed once, reused for all batches)
        self.logger.info("Computing mask-connectome overlap...")
        mask_voxel_indices, resampled_mask_img = self._get_mask_voxel_indices(mask_data)

        if len(mask_voxel_indices) == 0:
            self.logger.warning(
                "No mask voxels overlap with connectome brain mask after resampling "
                "— producing zero-valued network maps"
            )
            return self._build_empty_mask_results()

        self.logger.success(f"Found {len(mask_voxel_indices):,} overlapping mask voxels")

        # Store resampled mask as analysis_mask if keep_intermediate=True
        # This is stored as an instance variable to be added to results at the end
        self._resampled_mask_img = resampled_mask_img

        # Initialize accumulators for batch processing
        all_z_maps = []
        total_subjects = 0

        # Process each connectome batch sequentially
        connectome_files = self._get_connectome_files()
        n_batches = len(connectome_files)

        if n_batches == 1:
            self.logger.info("Processing connectome...")
        else:
            self.logger.info(f"Processing {n_batches} connectome batches...")

        for batch_idx, batch_file in enumerate(connectome_files, 1):
            self.logger.progress(f"Loading {batch_file.name}", current=batch_idx, total=n_batches)

            # Load this batch's timeseries
            with h5py.File(batch_file, "r") as hf:
                batch_timeseries = hf["timeseries"][:]
                batch_n_subjects = batch_timeseries.shape[0]

            self.logger.debug(
                f"Extracting mask timeseries ({batch_n_subjects} subjects)", indent_level=1
            )

            # Extract mask timeseries for this batch
            if self.method == "boes":
                mask_ts = self._extract_lesion_timeseries_boes_batch(
                    batch_timeseries, mask_voxel_indices
                )
            else:  # pini
                mask_ts = self._extract_lesion_timeseries_pini_batch(
                    batch_timeseries, mask_voxel_indices
                )

            self.logger.debug("Computing correlation maps", indent_level=1)

            # Compute correlation maps for this batch
            batch_r_maps = self._compute_correlation_maps_batch(mask_ts, batch_timeseries)

            self.logger.debug("Applying Fisher z-transform", indent_level=1)

            # Fisher z-transform
            batch_z_maps = np.arctanh(batch_r_maps)
            batch_z_maps = np.nan_to_num(batch_z_maps, nan=0, posinf=10, neginf=-10)

            # Accumulate
            all_z_maps.append(batch_z_maps)
            total_subjects += batch_timeseries.shape[0]

            # Explicitly free memory
            del batch_timeseries, mask_ts, batch_r_maps, batch_z_maps

        self.logger.info(f"Aggregating results across {total_subjects} subjects...")

        # Concatenate all z-maps
        all_z_maps_array = np.vstack(all_z_maps)

        # Aggregate across all subjects
        mean_z_map = np.mean(all_z_maps_array, axis=0)
        mean_r_map = np.tanh(mean_z_map)

        # Compute t-statistics (always computed)
        self.logger.info("Computing t-statistics...")
        std_z_map = np.std(all_z_maps_array, axis=0, ddof=1)
        with np.errstate(divide="ignore", invalid="ignore"):
            std_error_map_flat = std_z_map / np.sqrt(total_subjects)
            t_map_flat = np.zeros_like(mean_z_map)
            np.divide(
                mean_z_map,
                std_error_map_flat,
                out=t_map_flat,
                where=(std_error_map_flat != 0),
            )

        # Compute p-values from t-statistics (two-tailed)
        p_map_flat = None
        p_fdr_map_flat = None
        if self.compute_p_map:
            self.logger.info("Computing p-values...")
            df = total_subjects - 1  # degrees of freedom
            # Two-tailed p-value: 2 * (1 - CDF(|t|))
            p_map_flat = 2 * stats.t.sf(np.abs(t_map_flat), df)
            p_map_flat = np.nan_to_num(p_map_flat, nan=1.0)

            # Compute FDR-corrected p-values if requested
            if self.fdr_alpha is not None:
                self.logger.info(f"Computing FDR-corrected p-values (alpha={self.fdr_alpha})...")
                # Benjamini-Hochberg FDR correction
                n_voxels = len(p_map_flat)
                sorted_indices = np.argsort(p_map_flat)
                sorted_pvals = p_map_flat[sorted_indices]

                # Compute FDR-adjusted p-values
                # p_adj[i] = min(p[i] * n / (rank[i]), 1.0)
                ranks = np.arange(1, n_voxels + 1)
                adjusted = sorted_pvals * n_voxels / ranks

                # Enforce monotonicity (cumulative minimum from the right)
                adjusted = np.minimum.accumulate(adjusted[::-1])[::-1]
                adjusted = np.clip(adjusted, 0, 1)

                # Map back to original order
                p_fdr_map_flat = np.empty_like(p_map_flat)
                p_fdr_map_flat[sorted_indices] = adjusted

                # Count significant voxels
                n_significant_fdr = np.sum(p_fdr_map_flat < self.fdr_alpha)
                pct_significant_fdr = (n_significant_fdr / n_voxels) * 100
                self.logger.info(
                    f"Found {n_significant_fdr:,} voxels ({pct_significant_fdr:.2f}%) "
                    f"significant at FDR q < {self.fdr_alpha}",
                    indent_level=1,
                )

        # Free memory
        del all_z_maps, all_z_maps_array

        self.logger.info("Creating 3D output volumes...")

        # Convert flat arrays to 3D brain volumes
        mask_shape = self._mask_info["mask_shape"]
        mask_indices = self._mask_info["mask_indices"]
        mask_affine = self._mask_info["mask_affine"]

        # Create 3D volumes
        # mask_indices is tuple of (x_coords, y_coords, z_coords)
        correlation_map_3d = np.zeros(mask_shape, dtype=np.float32)
        correlation_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_r_map

        z_map_3d = np.zeros(mask_shape, dtype=np.float32)
        z_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_z_map

        # Create NIfTI images
        correlation_map_nifti = nib.Nifti1Image(correlation_map_3d, mask_affine)
        z_map_nifti = nib.Nifti1Image(z_map_3d, mask_affine)

        # Create t-map (always computed)
        t_map_3d = np.zeros(mask_shape, dtype=np.float32)
        t_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = t_map_flat
        t_map_nifti = nib.Nifti1Image(t_map_3d, mask_affine)

        # Create thresholded binary map if threshold provided
        t_threshold_map_nifti = None
        if self.t_threshold is not None:
            self.logger.info(f"Creating thresholded map (|t| > {self.t_threshold})")
            t_threshold_mask = np.abs(t_map_flat) > self.t_threshold
            n_significant = np.sum(t_threshold_mask)
            pct_significant = (n_significant / len(t_map_flat)) * 100

            threshold_map_3d = np.zeros(mask_shape, dtype=np.uint8)
            threshold_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = (
                t_threshold_mask.astype(np.uint8)
            )
            t_threshold_map_nifti = nib.Nifti1Image(threshold_map_3d, mask_affine)

            self.logger.info(
                f"Found {n_significant:,} voxels ({pct_significant:.2f}%) above threshold",
                indent_level=1,
            )

        # Create p-value maps if computed
        p_map_nifti = None
        p_fdr_map_nifti = None
        if self.compute_p_map and p_map_flat is not None:
            p_map_3d = np.zeros(mask_shape, dtype=np.float32)
            p_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = p_map_flat
            # Set background to 1.0 (not significant) for clarity
            p_map_nifti = nib.Nifti1Image(p_map_3d, mask_affine)

            # Create FDR-corrected p-value map if computed
            if p_fdr_map_flat is not None:
                p_fdr_map_3d = np.zeros(mask_shape, dtype=np.float32)
                p_fdr_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = p_fdr_map_flat
                p_fdr_map_nifti = nib.Nifti1Image(p_fdr_map_3d, mask_affine)

        # Compute summary statistics
        mean_correlation = float(np.mean(mean_r_map))
        std_correlation = float(np.std(mean_r_map))
        max_correlation = float(np.max(mean_r_map))
        min_correlation = float(np.min(mean_r_map))

        # Success summary
        summary_details = {
            "mean_correlation": mean_correlation,
            "std_correlation": std_correlation,
            "correlation_range": f"[{min_correlation:.4f}, {max_correlation:.4f}]",
            "n_subjects": total_subjects,
        }

        if t_map_nifti is not None:
            t_min = float(np.min(t_map_flat))
            t_max = float(np.max(t_map_flat))
            summary_details["t_range"] = f"[{t_min:.2f}, {t_max:.2f}]"

        self.logger.success("Analysis complete", details=summary_details)

        # Create result objects as dict with descriptive keys
        results = {}

        # Correlation map (r values)
        correlation_result = VoxelMap(
            name="rmap",
            data=correlation_map_nifti,
            space=self.output_space,
            resolution=self.output_resolution,
            metadata={
                "method": self.method,
                "n_subjects": total_subjects,
                "n_batches": len(connectome_files),
                "statistic": "pearson_correlation_coefficient",
            },
        )
        results["rmap"] = correlation_result

        # Z-map (Fisher z-transformed correlations)
        z_result = VoxelMap(
            name="zmap",
            data=z_map_nifti,
            space=self.output_space,
            resolution=self.output_resolution,
            metadata={
                "method": self.method,
                "n_subjects": total_subjects,
                "n_batches": len(connectome_files),
                "statistic": "fisher_z",
            },
        )
        results["zmap"] = z_result

        # Summary statistics
        summary_dict = {
            "mean": mean_correlation,
            "std": std_correlation,
            "max": max_correlation,
            "min": min_correlation,
            "n_subjects": total_subjects,
            "n_batches": len(connectome_files),
        }

        # Add t-map results if computed
        if t_map_nifti is not None:
            t_result = VoxelMap(
                name="tmap",
                data=t_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "t_statistic",
                },
            )
            results["tmap"] = t_result
            summary_dict["t_min"] = float(np.min(t_map_flat))
            summary_dict["t_max"] = float(np.max(t_map_flat))

        if t_threshold_map_nifti is not None:
            threshold_result = VoxelMap(
                name="tthresholdmap",
                data=t_threshold_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "threshold": self.t_threshold,
                    "statistic": "thresholded_t",
                },
            )
            results["tthresholdmap"] = threshold_result
            summary_dict["n_significant_voxels"] = int(n_significant)
            summary_dict["pct_significant_voxels"] = float(pct_significant)

        # Add p-value map if computed
        if p_map_nifti is not None:
            p_result = VoxelMap(
                name="pmap",
                data=p_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "p_value_two_tailed",
                    "degrees_of_freedom": total_subjects - 1,
                },
            )
            results["pmap"] = p_result
            summary_dict["p_min"] = float(np.min(p_map_flat))
            summary_dict["p_max"] = float(np.max(p_map_flat))

        # Add FDR-corrected p-value map if computed
        if p_fdr_map_nifti is not None:
            p_fdr_result = VoxelMap(
                name="pfdrmap",
                data=p_fdr_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "p_value_fdr_corrected",
                    "fdr_alpha": self.fdr_alpha,
                    "correction_method": "benjamini_hochberg",
                },
            )
            results["pfdrmap"] = p_fdr_result
            summary_dict["n_significant_fdr"] = int(n_significant_fdr)
            summary_dict["pct_significant_fdr"] = float(pct_significant_fdr)
            summary_dict["fdr_alpha"] = self.fdr_alpha

            # Create binary significance mask at FDR threshold
            fdr_sig_mask = (p_fdr_map_flat < self.fdr_alpha).astype(np.uint8)
            fdr_sig_map_3d = np.zeros(mask_shape, dtype=np.uint8)
            fdr_sig_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = fdr_sig_mask
            fdr_sig_map_nifti = nib.Nifti1Image(fdr_sig_map_3d, mask_affine)

            fdr_threshold_result = VoxelMap(
                name="pfdrthresholdmap",
                data=fdr_sig_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "fdr_significant_binary",
                    "fdr_alpha": self.fdr_alpha,
                    "n_significant": n_significant_fdr,
                },
            )
            results["pfdrthresholdmap"] = fdr_threshold_result

        # Add summary statistics as ScalarMetric
        summary_result = ScalarMetric(
            name="summarystatistics",
            data=summary_dict,
            metadata={
                "method": self.method,
                "n_subjects": total_subjects,
            },
        )
        results["summarystatistics"] = summary_result

        # Add analysis_mask if keep_intermediate=True
        # This is the mask that was actually used for computation (resampled to connectome grid)
        if self.keep_intermediate and hasattr(self, "_resampled_mask_img"):
            analysis_mask = VoxelMap(
                name="analysis_mask",
                data=self._resampled_mask_img,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "description": "Mask resampled to connectome space for analysis",
                    "original_space": mask_data.metadata.get(
                        "_original_input_space", mask_data.space
                    ),
                    "original_resolution": mask_data.metadata.get(
                        "_original_input_resolution", mask_data.resolution
                    ),
                    "analysis_space": self.output_space,
                    "analysis_resolution": self.output_resolution,
                },
            )
            results["analysis_mask"] = analysis_mask

        # Transform VoxelMap results back to input space if requested
        if self.return_in_input_space:
            results = self._transform_results_to_input_space(results, mask_data)

        self.logger.success(f"Analysis complete ({len(results)} results)")
        return results

    def _build_empty_mask_results(self) -> dict[str, "AnalysisResult"]:
        """Build zero-valued results for an empty mask.

        Returns the same structure as a normal analysis run (rmap, zmap,
        summarystatistics) but with all-zero NIfTI maps.
        """
        # Use connectome mask info for output shape/affine
        if self._mask_info is None:
            self._load_mask_info()

        mask_shape = self._mask_info["mask_shape"]
        mask_affine = self._mask_info["mask_affine"]

        zero_vol = np.zeros(mask_shape, dtype=np.float32)

        results: dict[str, AnalysisResult] = {}

        results["rmap"] = VoxelMap(
            name="rmap",
            data=nib.Nifti1Image(zero_vol.copy(), mask_affine),
            space=self.output_space,
            resolution=self.output_resolution,
            metadata={"method": self.method, "empty_mask": True},
        )

        results["zmap"] = VoxelMap(
            name="zmap",
            data=nib.Nifti1Image(zero_vol.copy(), mask_affine),
            space=self.output_space,
            resolution=self.output_resolution,
            metadata={"method": self.method, "empty_mask": True},
        )

        results["summarystatistics"] = ScalarMetric(
            name="summarystatistics",
            data={
                "mean": 0.0,
                "std": 0.0,
                "max": 0.0,
                "min": 0.0,
                "n_subjects": 0,
                "n_batches": 0,
                "empty_mask": True,
            },
            metadata={"method": self.method},
        )

        return results

    def _load_mask_info(self) -> tuple:
        """Load mask information from first connectome file.

        Mask info is shared across all batch files and only needs to be
        loaded once. Sets self._mask_info.

        Returns
        -------
        tuple
            (mask_indices, mask_affine, mask_shape) tuple

        Raises
        ------
        ValidationError
            If HDF5 file structure is invalid.
        """
        connectome_files = self._get_connectome_files()
        self._mask_info = read_mask_info(connectome_files[0])
        return (
            self._mask_info["mask_indices"],
            self._mask_info["mask_affine"],
            self._mask_info["mask_shape"],
        )

    def _extract_lesion_timeseries_boes_batch(
        self, batch_timeseries: np.ndarray, lesion_voxel_indices: np.ndarray
    ) -> np.ndarray:
        """Extract mean timeseries across all lesion voxels (BOES method).

        Memory-efficient version that works on a single batch.

        Parameters
        ----------
        batch_timeseries : np.ndarray
            Shape (n_subjects, n_timepoints, n_voxels). Connectome batch data.
        lesion_voxel_indices : np.ndarray
            1D array of voxel indices within the connectome mask.

        Returns
        -------
        np.ndarray
            Shape (n_subjects, n_timepoints). Mean timeseries for each subject.
        """
        # Extract timeseries for lesion voxels
        # Shape: (n_subjects, n_timepoints, n_lesion_voxels)
        lesion_ts = batch_timeseries[:, :, lesion_voxel_indices]

        # Compute mean across voxels
        # Shape: (n_subjects, n_timepoints)
        lesion_mean_ts = np.mean(lesion_ts, axis=2)

        return lesion_mean_ts

    def _extract_lesion_timeseries_pini_batch(
        self, batch_timeseries: np.ndarray, lesion_voxel_indices: np.ndarray
    ) -> np.ndarray:
        """Extract representative timeseries using PCA (PINI method).

        Memory-efficient version that works on a single batch.

        Uses PCA to identify most representative voxels based on their
        correlation with the mean timeseries, then extracts mean from
        these selected voxels.

        Parameters
        ----------
        batch_timeseries : np.ndarray
            Shape (n_subjects, n_timepoints, n_voxels). Connectome batch data.
        lesion_voxel_indices : np.ndarray
            1D array of voxel indices within the connectome mask.

        Returns
        -------
        np.ndarray
            Shape (n_subjects, n_timepoints). Representative timeseries.
        """
        # Extract timeseries for lesion voxels
        lesion_ts = batch_timeseries[:, :, lesion_voxel_indices]

        # Compute initial mean timeseries
        mean_ts_across_voxels = np.mean(lesion_ts, axis=2)

        # Center timeseries
        mean_ts_centered = mean_ts_across_voxels - mean_ts_across_voxels.mean(axis=1, keepdims=True)
        voxel_ts_centered = lesion_ts - lesion_ts.mean(axis=1, keepdims=True)

        # Compute correlation between mean and each voxel
        # Using einsum for efficiency
        covariance = np.einsum(
            "it,itv->iv",
            mean_ts_centered,
            voxel_ts_centered,
            dtype=np.float64,
            optimize="optimal",
        )

        std_mean = np.sqrt(np.sum(mean_ts_centered**2, axis=1))
        std_voxels = np.sqrt(np.sum(voxel_ts_centered**2, axis=1))

        with np.errstate(divide="ignore", invalid="ignore"):
            pca_input_matrix = covariance / (std_mean[:, np.newaxis] * std_voxels)

        pca_input_matrix = np.nan_to_num(pca_input_matrix)

        # Apply PCA to find principal component
        pca = PCA(n_components=1)
        pca.fit(pca_input_matrix)
        pc1_loadings = pca.components_[0, :]

        # Select voxels above percentile threshold
        threshold = np.percentile(np.abs(pc1_loadings), self.pini_percentile)
        suprathreshold_indices = np.where(np.abs(pc1_loadings) >= threshold)[0]

        # Extract refined timeseries from selected voxels
        if len(suprathreshold_indices) == 0:
            # Fallback to all voxels if none selected
            refined_lesion_ts = lesion_ts
        else:
            refined_lesion_ts = lesion_ts[:, :, suprathreshold_indices]

        # Return mean timeseries from selected voxels
        return np.mean(refined_lesion_ts, axis=2)

    def _get_mask_voxel_indices(self, mask_data: SubjectData) -> tuple[np.ndarray, nib.Nifti1Image]:
        """Get indices of mask voxels within connectome mask (vectorized O(N) version).

        This uses a lookup array for O(N) complexity instead of O(N×M) nested loops,
        providing massive speedup for large masks.

        **Performance**:
        - Speedup: 15-2000x vs. legacy implementation (increases with mask size)

        **Benchmark Results** (MNI152 @ 2mm, ~335K brain voxels):
        - 100 voxels: 113ms → 7ms (15.7x speedup)
        - 1,000 voxels: 1,078ms → 4.7ms (228x speedup)
        - 10,000 voxels: 9,965ms → 4.9ms (2,025x speedup)

        **Implementation**:
        Uses a 3D lookup array (shape=mask_shape, dtype=int32) that maps
        spatial coordinates directly to flat indices in the connectome mask.
        This eliminates the need for searching through mask_indices for each
        mask voxel.

        Automatically resamples mask to connectome space if dimensions don't match.

        Parameters
        ----------
        mask_data : SubjectData
            Mask data in target space.

        Returns
        -------
        tuple[np.ndarray, nib.Nifti1Image]
            Tuple containing:
            - 1D array of indices into the connectome's voxel dimension
            - The mask image resampled to connectome space (for analysis_mask storage)

        Notes
        -----
        For batch processing scenarios, consider caching the lookup array to avoid
        rebuilding it for every subject
        """
        # Get connectome mask info
        mask_shape = self._mask_info["mask_shape"]
        mask_indices = self._mask_info["mask_indices"]
        mask_affine = self._mask_info["mask_affine"]

        # Get mask image
        mask_img = mask_data.mask_img
        input_shape = mask_img.shape

        # Check if resampling is needed
        if input_shape != mask_shape:
            # Resample mask to connectome space
            from nilearn.image import resample_to_img

            # Create template image in connectome space
            template_img = nib.Nifti1Image(np.zeros(mask_shape), mask_affine)

            # Resample mask to match connectome
            mask_img_resampled = resample_to_img(
                mask_img,
                template_img,
                interpolation="nearest",
                force_resample=True,
                copy_header=True,
            )
            resampled_mask = mask_img_resampled.get_fdata().astype(bool)
        else:
            mask_img_resampled = mask_img
            resampled_mask = mask_img.get_fdata().astype(bool)

        # Get mask coordinates
        mask_coords = np.where(resampled_mask)

        # Build lookup array: 3D array mapping coordinates to flat indices
        lookup = np.full(mask_shape, -1, dtype=np.int32)
        lookup[mask_indices] = np.arange(len(mask_indices[0]), dtype=np.int32)

        # Direct O(N) indexing to get flat indices
        flat_indices = lookup[mask_coords]

        # Filter out voxels not in connectome mask (value = -1)
        valid_indices = flat_indices[flat_indices >= 0]

        return valid_indices.astype(int), mask_img_resampled

    def _compute_correlation_maps_batch(
        self, lesion_timeseries: np.ndarray, batch_timeseries: np.ndarray
    ) -> np.ndarray:
        """Compute correlation maps between lesion and whole-brain timeseries.

        Memory-efficient version that works on a single batch.

        Parameters
        ----------
        lesion_timeseries : np.ndarray
            Shape (n_subjects, n_timepoints). Lesion timeseries.
        batch_timeseries : np.ndarray
            Shape (n_subjects, n_timepoints, n_voxels). Connectome batch data.

        Returns
        -------
        np.ndarray
            Shape (n_subjects, n_voxels). Correlation values for each voxel.
        """
        # Center timeseries
        brain_ts_centered = batch_timeseries - batch_timeseries.mean(axis=1, keepdims=True)
        lesion_ts_centered = lesion_timeseries - lesion_timeseries.mean(axis=1, keepdims=True)

        # Compute covariance using einsum
        # (n_subjects, n_timepoints) @ (n_subjects, n_timepoints, n_voxels)
        cov = np.einsum(
            "it,itv->iv",
            lesion_ts_centered,
            brain_ts_centered,
            dtype=np.float64,
            optimize="optimal",
        )

        # Compute standard deviations
        lesion_std = np.sqrt(np.sum(lesion_ts_centered**2, axis=1))
        brain_std = np.sqrt(np.sum(brain_ts_centered**2, axis=1))

        # Compute correlation (cov / (std_lesion * std_brain))
        with np.errstate(divide="ignore", invalid="ignore"):
            r_maps = cov / (lesion_std[:, np.newaxis] * brain_std)

        # Clean up NaN and inf values
        r_maps = np.nan_to_num(r_maps, nan=0, posinf=1, neginf=-1)

        return r_maps.astype(np.float32)

    def run_batch(self, mask_data_list: list[SubjectData]) -> list[SubjectData]:
        """Process multiple lesions together using vectorized operations.

        This is 10-50x faster than sequential processing because it:
        1. Processes all lesions through each connectome batch together
        2. Uses vectorized einsum: "lit,itv->liv" for batch correlations
        3. Minimizes loop overhead and leverages optimized BLAS operations

        This method is automatically called when using VectorizedStrategy
        for batch processing.

        Parameters
        ----------
        mask_data_list : list[SubjectData]
            Batch of lesions to process together

        Returns
        -------
        list[SubjectData]
            Processed lesions with results added

        Examples
        --------
        >>> from lacuna.batch import VectorizedStrategy
        >>> from lacuna.analysis import FunctionalNetworkMapping
        >>>
        >>> analysis = FunctionalNetworkMapping(...)
        >>> strategy = VectorizedStrategy()
        >>> results = strategy.execute(mask_data_list, analysis)
        """
        self.logger.info(f"Vectorized batch processing: {len(mask_data_list)} masks")

        # Validate all lesions first
        for mask_data in mask_data_list:
            self._validate_inputs(mask_data)

        # Load mask info once (shared across all batches)
        mask_indices, mask_affine, mask_shape = self._load_mask_info()
        connectome_files = self._get_connectome_files()

        self.logger.success(
            "Loaded connectome metadata",
            details={
                "connectome_batches": len(connectome_files),
                "mask_shape": mask_shape,
                "n_masks": len(mask_data_list),
            },
        )

        # Track empty masks — they get zero-valued results without batch processing
        empty_mask_indices: dict[int, dict] = {}
        for i, mask_data in enumerate(mask_data_list):
            if mask_data.is_empty_mask:
                subject_id = self._format_subject_id(mask_data)
                self.logger.warning(
                    f"Empty mask for {subject_id} — will produce zero-valued network maps"
                )
                empty_mask_indices[i] = self._build_empty_mask_results()

        mask_batch = []

        for i, mask_data in enumerate(mask_data_list):
            if i in empty_mask_indices:
                continue
            subject_id = self._format_subject_id(mask_data)

            voxel_indices, _ = self._get_mask_voxel_indices(mask_data)

            if len(voxel_indices) == 0:
                self.logger.warning(
                    f"No overlap for {subject_id}: mask has no voxels within connectome "
                    f"brain mask after resampling to {self.TARGET_SPACE} "
                    f"— will produce zero-valued network maps"
                )
                empty_mask_indices[i] = self._build_empty_mask_results()
                continue

            mask_batch.append(
                {
                    "mask_data": mask_data,
                    "voxel_indices": voxel_indices,
                    "index": i,
                }
            )

        # Process through all connectome batches (VECTORIZED)
        if not mask_batch:
            # All masks were empty or had no overlap — skip batch processing
            self.logger.info("All masks empty or non-overlapping — skipping connectome processing")
            results = [
                mask_data_list[idx].add_result(self.__class__.__name__, empty_result)
                for idx, empty_result in sorted(empty_mask_indices.items())
            ]
            return results

        self.logger.info("Processing connectome batches...")

        # Get number of voxels from first connectome batch
        with h5py.File(connectome_files[0], "r") as hf:
            n_voxels = hf["timeseries"].shape[2]

        # Initialize streaming aggregators for each mask (MEMORY OPTIMIZED)
        # Instead of storing all correlation maps, we accumulate statistics
        aggregators = []
        for _ in range(len(mask_batch)):
            aggregators.append(
                {
                    "sum_z": np.zeros(n_voxels, dtype=np.float64),  # Need higher precision for sums
                    "sum_z2": np.zeros(n_voxels, dtype=np.float64),
                    "n": 0,
                }
            )

        total_subjects = 0
        batch_times = []  # Track timing for each batch

        for batch_idx, connectome_path in enumerate(connectome_files):
            import time

            batch_start_time = time.time()

            with h5py.File(connectome_path, "r") as hf:
                timeseries_data = hf["timeseries"][:]  # (n_subj, n_time, n_vox)
                n_subjects = timeseries_data.shape[0]
                total_subjects += n_subjects

            # Vectorized processing for ALL masks at once
            batch_r_maps = self._compute_batch_correlations_vectorized(
                mask_batch, timeseries_data
            )  # (n_masks, n_subjects, n_voxels)

            # Convert to Fisher z-scores and update running statistics
            # This is the KEY optimization: we don't store full maps!
            with np.errstate(divide="ignore", invalid="ignore"):
                batch_z_maps = np.arctanh(batch_r_maps)
                batch_z_maps = np.nan_to_num(batch_z_maps, nan=0.0, posinf=3.0, neginf=-3.0)

            # Update aggregators with streaming statistics
            for i in range(len(mask_batch)):
                # Sum across subjects in this batch
                aggregators[i]["sum_z"] += np.sum(batch_z_maps[i], axis=0)
                aggregators[i]["sum_z2"] += np.sum(batch_z_maps[i] ** 2, axis=0)
                aggregators[i]["n"] += n_subjects

            # Memory cleanup - immediately free large arrays
            del timeseries_data, batch_r_maps, batch_z_maps

            # Display timing information
            batch_elapsed = time.time() - batch_start_time
            batch_times.append(batch_elapsed)

            # Show progress with time estimates
            if len(batch_times) > 2:
                avg_time = sum(batch_times) / len(batch_times)
                remaining_batches = len(connectome_files) - (batch_idx + 1)
                est_remaining = avg_time * remaining_batches
                self.logger.progress(
                    f"Batch completed in {batch_elapsed:.2f}s (est. {est_remaining:.1f}s remaining)",
                    current=batch_idx + 1,
                    total=len(connectome_files),
                )
            else:
                self.logger.progress(
                    f"Batch completed in {batch_elapsed:.2f}s",
                    current=batch_idx + 1,
                    total=len(connectome_files),
                )

        total_batch_time = sum(batch_times)
        avg_batch_time = total_batch_time / len(batch_times) if batch_times else 0
        self.logger.success(
            "Batch processing complete",
            details={
                "n_batches": len(connectome_files),
                "total_time": f"{total_batch_time:.2f}s",
                "avg_time_per_batch": f"{avg_batch_time:.2f}s",
            },
        )

        # Compute final statistics from aggregated values
        self.logger.info("Aggregating results...")

        # Define per-subject aggregation work
        def _aggregate_one(i, mask_info):
            """Aggregate results for a single mask (thread-safe)."""
            n = aggregators[i]["n"]
            mean_z = aggregators[i]["sum_z"] / n
            mean_r = np.tanh(mean_z).astype(np.float32)

            # Sample standard deviation with Bessel's correction
            var_z_population = (aggregators[i]["sum_z2"] / n) - (mean_z**2)
            var_z_sample = (n / (n - 1)) * var_z_population
            std_z = np.sqrt(np.maximum(var_z_sample, 0))

            mask_copy = mask_info["mask_data"].copy()

            result = self._aggregate_results_from_statistics(
                mask_copy,
                mean_r,
                mean_z,
                std_z,
                mask_indices,
                mask_affine,
                mask_shape,
                total_subjects,
            )
            return mask_info["index"], result

        effective_n_jobs = self.n_jobs
        use_parallel = effective_n_jobs != 1 and len(mask_batch) > 1

        if use_parallel:
            from joblib import Parallel, delayed, parallel_backend

            self.logger.info(
                f"Aggregating {len(mask_batch)} subjects in parallel "
                f"(n_jobs={effective_n_jobs})"
            )
            # inner_max_num_threads=1 prevents fork-after-BLAS-init deadlock on
            # many-core nodes: caps OMP/MKL/OpenBLAS thread pools to 1 in each
            # worker so no inherited mutex can be left locked by a thread that
            # didn't survive fork().
            with parallel_backend("loky", inner_max_num_threads=1):
                pairs = Parallel(n_jobs=effective_n_jobs)(
                    delayed(_aggregate_one)(i, mask_info) for i, mask_info in enumerate(mask_batch)
                )
            processed_results = dict(pairs)
        else:
            processed_results = {}
            for i, mask_info in enumerate(mask_batch):
                subject_id = self._format_subject_id(mask_info["mask_data"])
                self.logger.info(f"Aggregating results for: {subject_id}", indent_level=1)
                idx, result = _aggregate_one(i, mask_info)
                processed_results[idx] = result

        # Merge empty-mask results back in at their original indices
        # Wrap raw result dicts into SubjectData to match _aggregate_one output
        for idx, empty_result in empty_mask_indices.items():
            processed_results[idx] = mask_data_list[idx].add_result(
                self.__class__.__name__, empty_result
            )

        # Build final results list in original input order
        results = [processed_results[k] for k in sorted(processed_results)]

        self.logger.success(
            "Batch processing complete",
            details={
                "n_masks_processed": len(processed_results),
            },
        )

        return results

    def _compute_batch_correlations_vectorized(
        self,
        mask_batch: list[dict],
        timeseries_data: np.ndarray,
    ) -> np.ndarray:
        """Compute correlations for ALL masks at once (vectorized).

        Uses einsum "lit,itv->liv" to compute correlations for all masks simultaneously,
        reducing overhead and enabling optimized BLAS operations.

        Parameters
        ----------
        mask_batch : list[dict]
            List of mask dictionaries with 'voxel_indices' and 'mask_data'
        timeseries_data : np.ndarray
            Shape (n_subjects, n_timepoints, n_voxels). Connectome batch data.

        Returns
        -------
        np.ndarray
            Shape (n_masks, n_subjects, n_voxels). Correlation maps for all masks.
        """
        # Extract and process timeseries for all masks
        mask_mean_ts_list = []
        for mask_info in mask_batch:
            voxel_indices = mask_info["voxel_indices"]

            # Extract mask timeseries: (n_subjects, n_timepoints, n_mask_voxels)
            mask_ts = timeseries_data[:, :, voxel_indices]

            if self.method == "boes":
                # Simple mean across voxels
                mask_mean_ts = np.mean(mask_ts, axis=2)

            elif self.method == "pini":
                # PINI: PCA-based selection
                mask_mean_ts = self._compute_pini_timeseries_batch(mask_ts)

            mask_mean_ts_list.append(mask_mean_ts)

        # Stack into (n_masks, n_subjects, n_timepoints)
        mask_mean_ts_batch = np.stack(mask_mean_ts_list, axis=0)

        # Center data
        brain_ts_centered = timeseries_data - timeseries_data.mean(axis=1, keepdims=True)
        mask_ts_centered = mask_mean_ts_batch - mask_mean_ts_batch.mean(axis=2, keepdims=True)

        # VECTORIZED CORRELATION: Process all masks at once!
        # einsum: "lit,itv->liv"
        #   l = masks, i = subjects, t = timepoints, v = voxels
        # Use float32 for memory efficiency (sufficient precision for correlations)
        cov = np.einsum(
            "lit,itv->liv",
            mask_ts_centered.astype(np.float32),
            brain_ts_centered.astype(np.float32),
            dtype=np.float32,
            optimize="optimal",
        )

        # Compute standard deviations
        mask_std = np.sqrt(np.sum(mask_ts_centered**2, axis=2))  # (n_masks, n_subjects)
        brain_std = np.sqrt(np.sum(brain_ts_centered**2, axis=1))  # (n_subjects, n_voxels)

        # Compute correlations: cov / (mask_std * brain_std)
        with np.errstate(divide="ignore", invalid="ignore"):
            all_r_maps = cov / (mask_std[:, :, np.newaxis] * brain_std[np.newaxis, :, :])

        # Clean up NaN and inf values
        all_r_maps = np.nan_to_num(all_r_maps, nan=0, posinf=1, neginf=-1)

        return all_r_maps.astype(np.float32)

    def _compute_pini_timeseries_batch(self, lesion_ts: np.ndarray) -> np.ndarray:
        """Compute PINI timeseries for a batch of subjects.

        Parameters
        ----------
        lesion_ts : np.ndarray
            Shape (n_subjects, n_timepoints, n_lesion_voxels)

        Returns
        -------
        np.ndarray
            Shape (n_subjects, n_timepoints). PINI-refined timeseries.
        """
        n_subjects, n_timepoints, n_voxels = lesion_ts.shape

        # Compute mean timeseries first
        mean_ts_across_voxels = np.mean(lesion_ts, axis=2)  # (n_subjects, n_timepoints)

        # Center timeseries
        mean_ts_centered = mean_ts_across_voxels - mean_ts_across_voxels.mean(axis=1, keepdims=True)
        voxel_ts_centered = lesion_ts - lesion_ts.mean(axis=1, keepdims=True)

        # Compute covariance between mean and each voxel
        covariance = np.einsum(
            "it,itv->iv",
            mean_ts_centered,
            voxel_ts_centered,
            dtype=np.float64,
            optimize="optimal",
        )

        # Compute standard deviations
        std_mean = np.sqrt(np.sum(mean_ts_centered**2, axis=1))
        std_voxels = np.sqrt(np.sum(voxel_ts_centered**2, axis=1))

        # Compute correlation matrix for PCA
        with np.errstate(divide="ignore", invalid="ignore"):
            pca_input_matrix = covariance / (std_mean[:, np.newaxis] * std_voxels)

        pca_input_matrix = np.nan_to_num(pca_input_matrix)

        # Apply PCA to find principal component
        pca = PCA(n_components=1)
        pca.fit(pca_input_matrix)
        pc1_loadings = pca.components_[0, :]

        # Select voxels based on percentile threshold
        threshold = np.percentile(np.abs(pc1_loadings), self.pini_percentile)
        suprathreshold_indices = np.where(np.abs(pc1_loadings) >= threshold)[0]

        # Use refined voxel set or fall back to all voxels
        if len(suprathreshold_indices) == 0:
            refined_lesion_ts = lesion_ts
        else:
            refined_lesion_ts = lesion_ts[:, :, suprathreshold_indices]

        # Return mean of refined voxels
        return np.mean(refined_lesion_ts, axis=2)

    def _aggregate_results_from_statistics(
        self,
        mask_data: SubjectData,
        mean_r_map: np.ndarray,
        mean_z_map: np.ndarray,
        std_z_map: np.ndarray,
        mask_indices: tuple,
        mask_affine: np.ndarray,
        mask_shape: tuple,
        total_subjects: int,
    ) -> SubjectData:
        """Aggregate results from pre-computed statistics (memory-optimized).

        This method is used by vectorized batch processing with streaming
        aggregation. Instead of storing all individual correlation maps,
        it accepts pre-computed mean and standard deviation, reducing memory usage.

        Parameters
        ----------
        mask_data : SubjectData
            Original lesion data to add results to
        mean_r_map : np.ndarray
            Shape (n_voxels,). Mean correlation map (already Fisher z-averaged).
        mean_z_map : np.ndarray
            Shape (n_voxels,). Mean Fisher z-transformed map.
        std_z_map : np.ndarray
            Shape (n_voxels,). Standard deviation of z-maps (for t-statistics).
        mask_indices : tuple
            Brain mask voxel indices
        mask_affine : np.ndarray
            Affine transformation matrix
        mask_shape : tuple
            3D mask shape
        total_subjects : int
            Total number of subjects processed

        Returns
        -------
        SubjectData
            Lesion data with analysis results added
        """
        # Compute t-statistics (always computed)
        with np.errstate(divide="ignore", invalid="ignore"):
            std_error_map = std_z_map / np.sqrt(total_subjects)
            t_map_flat = np.zeros_like(mean_z_map)
            np.divide(
                mean_z_map,
                std_error_map,
                out=t_map_flat,
                where=(std_error_map != 0),
            )

        # Compute p-values from t-statistics (two-tailed)
        p_map_flat = None
        p_fdr_map_flat = None
        n_significant_fdr = 0
        pct_significant_fdr = 0.0

        if self.compute_p_map:
            df = total_subjects - 1  # degrees of freedom
            # Two-tailed p-value: 2 * (1 - CDF(|t|))
            p_map_flat = 2 * stats.t.sf(np.abs(t_map_flat), df)
            p_map_flat = np.nan_to_num(p_map_flat, nan=1.0)

            # Compute FDR-corrected p-values if requested
            if self.fdr_alpha is not None:
                # Benjamini-Hochberg FDR correction
                n_voxels = len(p_map_flat)
                sorted_indices = np.argsort(p_map_flat)
                sorted_pvals = p_map_flat[sorted_indices]

                # Compute FDR-adjusted p-values
                # p_adj[i] = min(p[i] * n / (rank[i]), 1.0)
                ranks = np.arange(1, n_voxels + 1)
                adjusted = sorted_pvals * n_voxels / ranks

                # Enforce monotonicity (cumulative minimum from the right)
                adjusted = np.minimum.accumulate(adjusted[::-1])[::-1]
                adjusted = np.clip(adjusted, 0, 1)

                # Map back to original order
                p_fdr_map_flat = np.empty_like(p_map_flat)
                p_fdr_map_flat[sorted_indices] = adjusted

                # Count significant voxels
                n_significant_fdr = int(np.sum(p_fdr_map_flat < self.fdr_alpha))
                pct_significant_fdr = (n_significant_fdr / n_voxels) * 100

        # Create 3D volumes
        correlation_map_3d = np.zeros(mask_shape, dtype=np.float32)
        correlation_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_r_map
        correlation_map_nifti = nib.Nifti1Image(correlation_map_3d, mask_affine)

        z_map_3d = np.zeros(mask_shape, dtype=np.float32)
        z_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_z_map.astype(np.float32)
        z_map_nifti = nib.Nifti1Image(z_map_3d, mask_affine)

        # Build results dictionary with snake_case keys (matching _run_analysis)
        # Wrap NIfTI images in VoxelMap for consistent unwrap behavior
        results = {
            "rmap": VoxelMap(
                name="rmap",
                data=correlation_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "pearson_correlation_coefficient",
                },
            ),
            "zmap": VoxelMap(
                name="zmap",
                data=z_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "fisher_z",
                },
            ),
            "summarystatistics": ScalarMetric(
                name="summarystatistics",
                data={
                    "mean": float(np.mean(mean_r_map)),
                    "std": float(np.std(mean_r_map)),
                    "max": float(np.max(mean_r_map)),
                    "min": float(np.min(mean_r_map)),
                    "n_subjects": total_subjects,
                    "n_batches": len(self._get_connectome_files()),
                },
                metadata={"method": self.method},
            ),
        }

        # Add t-map results if computed
        if t_map_flat is not None:
            t_map_3d = np.zeros(mask_shape, dtype=np.float32)
            t_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = t_map_flat
            t_map_nifti = nib.Nifti1Image(t_map_3d, mask_affine)

            results["tmap"] = VoxelMap(
                name="tmap",
                data=t_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "t_statistic",
                },
            )
            # Update summary statistics with t-map info
            results["summarystatistics"].data["t_min"] = float(np.min(t_map_flat))
            results["summarystatistics"].data["t_max"] = float(np.max(t_map_flat))

            # Create thresholded t-map if threshold provided
            if self.t_threshold is not None:
                t_threshold_mask = np.abs(t_map_flat) > self.t_threshold
                threshold_map_3d = np.zeros(mask_shape, dtype=np.uint8)
                threshold_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = (
                    t_threshold_mask.astype(np.uint8)
                )
                t_threshold_map_nifti = nib.Nifti1Image(threshold_map_3d, mask_affine)
                results["tthresholdmap"] = VoxelMap(
                    name="tthresholdmap",
                    data=t_threshold_map_nifti,
                    space=self.output_space,
                    resolution=self.output_resolution,
                    metadata={
                        "method": self.method,
                        "threshold": self.t_threshold,
                        "statistic": "thresholded_t",
                    },
                )
                results["summarystatistics"].data["t_threshold"] = self.t_threshold
                results["summarystatistics"].data["n_significant_voxels"] = int(
                    np.sum(t_threshold_mask)
                )

        # Add p-value map if computed
        if p_map_flat is not None:
            p_map_3d = np.zeros(mask_shape, dtype=np.float32)
            p_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = p_map_flat
            p_map_nifti = nib.Nifti1Image(p_map_3d, mask_affine)

            results["pmap"] = VoxelMap(
                name="pmap",
                data=p_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "p_value_two_tailed",
                    "degrees_of_freedom": total_subjects - 1,
                },
            )
            results["summarystatistics"].data["p_min"] = float(np.min(p_map_flat))
            results["summarystatistics"].data["p_max"] = float(np.max(p_map_flat))

        # Add FDR-corrected p-value map if computed
        if p_fdr_map_flat is not None:
            p_fdr_map_3d = np.zeros(mask_shape, dtype=np.float32)
            p_fdr_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = p_fdr_map_flat
            p_fdr_map_nifti = nib.Nifti1Image(p_fdr_map_3d, mask_affine)

            results["pfdrmap"] = VoxelMap(
                name="pfdrmap",
                data=p_fdr_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "p_value_fdr_corrected",
                    "fdr_alpha": self.fdr_alpha,
                    "correction_method": "benjamini_hochberg",
                },
            )
            results["summarystatistics"].data["n_significant_fdr"] = n_significant_fdr
            results["summarystatistics"].data["pct_significant_fdr"] = pct_significant_fdr
            results["summarystatistics"].data["fdr_alpha"] = self.fdr_alpha

            # Create binary significance mask at FDR threshold
            fdr_sig_mask = (p_fdr_map_flat < self.fdr_alpha).astype(np.uint8)
            fdr_sig_map_3d = np.zeros(mask_shape, dtype=np.uint8)
            fdr_sig_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = fdr_sig_mask
            fdr_sig_map_nifti = nib.Nifti1Image(fdr_sig_map_3d, mask_affine)

            results["pfdrthresholdmap"] = VoxelMap(
                name="pfdrthresholdmap",
                data=fdr_sig_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "fdr_significant_binary",
                    "fdr_alpha": self.fdr_alpha,
                    "n_significant": n_significant_fdr,
                },
            )

        # Transform VoxelMap results back to input space if requested
        if self.return_in_input_space:
            results = self._transform_results_to_input_space(results, mask_data)

        # Add results to mask data (returns new instance with results)
        batch_results = {
            "rmap": results["rmap"],
            "zmap": results["zmap"],
            "summarystatistics": results["summarystatistics"],
        }
        # Add optional results if present
        if "tmap" in results:
            batch_results["tmap"] = results["tmap"]
        if "tthresholdmap" in results:
            batch_results["tthresholdmap"] = results["tthresholdmap"]
        if "pmap" in results:
            batch_results["pmap"] = results["pmap"]
        if "pfdrmap" in results:
            batch_results["pfdrmap"] = results["pfdrmap"]
        if "pfdrthresholdmap" in results:
            batch_results["pfdrthresholdmap"] = results["pfdrthresholdmap"]

        mask_data_with_results = mask_data.add_result(self.__class__.__name__, batch_results)

        return mask_data_with_results

    def _aggregate_results(
        self,
        mask_data: SubjectData,
        all_r_maps: np.ndarray,
        mask_indices: tuple,
        mask_affine: np.ndarray,
        mask_shape: tuple,
        total_subjects: int,
    ) -> SubjectData:
        """Aggregate correlation maps across all subjects into final results.

        This method is reused by both single and batch processing.

        Parameters
        ----------
        mask_data : SubjectData
            Original lesion data to add results to
        all_r_maps : np.ndarray
            Shape (n_subjects, n_voxels). All correlation maps.
        mask_indices : tuple
            Brain mask voxel indices
        mask_affine : np.ndarray
            Affine transformation matrix
        mask_shape : tuple
            3D mask shape
        total_subjects : int
            Total number of subjects processed

        Returns
        -------
        SubjectData
            Lesion data with analysis results added
        """
        # Fisher z-transform
        all_z_maps = np.arctanh(all_r_maps)
        all_z_maps = np.nan_to_num(all_z_maps, nan=0, posinf=10, neginf=-10)

        # Compute statistics
        mean_z_map = np.mean(all_z_maps, axis=0)
        mean_r_map = np.tanh(mean_z_map)

        # Compute t-statistics (always computed)
        std_z_map = np.std(all_z_maps, axis=0, ddof=1)

        with np.errstate(divide="ignore", invalid="ignore"):
            std_error_map = std_z_map / np.sqrt(total_subjects)
            t_map_flat = np.zeros_like(mean_z_map)
            np.divide(
                mean_z_map,
                std_error_map,
                out=t_map_flat,
                where=(std_error_map != 0),
            )

        # Create 3D volumes
        correlation_map_3d = np.zeros(mask_shape, dtype=np.float32)
        correlation_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_r_map
        correlation_map_nifti = nib.Nifti1Image(correlation_map_3d, mask_affine)

        z_map_3d = np.zeros(mask_shape, dtype=np.float32)
        z_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = mean_z_map
        z_map_nifti = nib.Nifti1Image(z_map_3d, mask_affine)

        # Build results dictionary
        # Wrap NIfTI images in VoxelMap for consistent unwrap behavior
        results = {
            "rmap": VoxelMap(
                name="rmap",
                data=correlation_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "pearson_correlation_coefficient",
                },
            ),
            "network_map": correlation_map_nifti,  # Alias for backward compat (raw nifti)
            "zmap": VoxelMap(
                name="zmap",
                data=z_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "fisher_z",
                },
            ),
            "mean_correlation": float(np.mean(mean_r_map)),
            "summarystatistics": ScalarMetric(
                name="summarystatistics",
                data={
                    "mean": float(np.mean(mean_r_map)),
                    "std": float(np.std(mean_r_map)),
                    "max": float(np.max(mean_r_map)),
                    "min": float(np.min(mean_r_map)),
                    "n_subjects": total_subjects,
                    "n_batches": len(self._get_connectome_files()),
                },
                metadata={"method": self.method},
            ),
        }

        # Add t-map results if computed
        if t_map_flat is not None:
            t_map_3d = np.zeros(mask_shape, dtype=np.float32)
            t_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = t_map_flat
            t_map_nifti = nib.Nifti1Image(t_map_3d, mask_affine)

            results["tmap"] = VoxelMap(
                name="tmap",
                data=t_map_nifti,
                space=self.output_space,
                resolution=self.output_resolution,
                metadata={
                    "method": self.method,
                    "n_subjects": total_subjects,
                    "statistic": "t_statistic",
                },
            )
            results["summarystatistics"].data["t_min"] = float(np.min(t_map_flat))
            results["summarystatistics"].data["t_max"] = float(np.max(t_map_flat))

            # Create thresholded t-map if threshold provided
            if self.t_threshold is not None:
                t_threshold_mask = np.abs(t_map_flat) > self.t_threshold
                threshold_map_3d = np.zeros(mask_shape, dtype=np.uint8)
                threshold_map_3d[mask_indices[0], mask_indices[1], mask_indices[2]] = (
                    t_threshold_mask.astype(np.uint8)
                )
                t_threshold_map_nifti = nib.Nifti1Image(threshold_map_3d, mask_affine)
                results["tthresholdmap"] = VoxelMap(
                    name="tthresholdmap",
                    data=t_threshold_map_nifti,
                    space=self.output_space,
                    resolution=self.output_resolution,
                    metadata={
                        "method": self.method,
                        "threshold": self.t_threshold,
                        "statistic": "thresholded_t",
                    },
                )
                results["summarystatistics"].data["t_threshold"] = self.t_threshold
                results["summarystatistics"].data["n_significant_voxels"] = int(
                    np.sum(t_threshold_mask)
                )

        # Add results to lesion data (returns new instance with results)
        # Note: Using individual keys to match _run_analysis() structure
        batch_results = {
            "rmap": results["rmap"],
            "zmap": results["zmap"],
            "summarystatistics": results["summarystatistics"],
        }
        # Add optional results if present
        if "tmap" in results:
            batch_results["tmap"] = results["tmap"]
        if "tthresholdmap" in results:
            batch_results["tthresholdmap"] = results["tthresholdmap"]

        mask_data_with_results = mask_data.add_result(self.__class__.__name__, batch_results)

        return mask_data_with_results

    def _get_parameters(self) -> dict:
        """Get analysis parameters for provenance and display.

        Returns
        -------
        dict
            Dictionary of parameter names and values.
        """
        return {
            "connectome_name": self.connectome_name,
            "method": self.method,
            "pini_percentile": self.pini_percentile,
            "n_jobs": self.n_jobs,
            "compute_p_map": self.compute_p_map,
            "fdr_alpha": self.fdr_alpha,
            "t_threshold": self.t_threshold,
            "return_in_input_space": self.return_in_input_space,
            "output_resolution": self.final_output_resolution,
            "keep_intermediate": self.keep_intermediate,
            "analysis_space": self.output_space,
            "analysis_resolution": self.output_resolution,
            "verbose": self.verbose,
        }

    def _transform_results_to_input_space(self, results: dict, mask_data: SubjectData) -> dict:
        """Transform VoxelMap results back to original input space.

        Parameters
        ----------
        results : dict
            Dictionary of result objects
        mask_data : SubjectData
            Input mask data with space/resolution metadata. If the mask was
            transformed by BaseAnalysis, it will have _original_input_space
            and _original_input_resolution in metadata.

        Returns
        -------
        dict
            Results with transformed VoxelMap objects

        Raises
        ------
        ValueError
            If mask_data lacks space or resolution metadata
        """
        from lacuna.core.spaces import REFERENCE_AFFINES, CoordinateSpace
        from lacuna.spatial.transform import transform_image

        # Get original input space from metadata (if available) or current space
        target_space_id = mask_data.metadata.get("_original_input_space", mask_data.space)

        # Determine target resolution:
        # 1. If user specified output_resolution, use that
        # 2. Otherwise, match the original input resolution
        if self.final_output_resolution is not None:
            target_resolution = self.final_output_resolution
        else:
            target_resolution = mask_data.metadata.get(
                "_original_input_resolution", mask_data.resolution
            )

        # If we're already in the target space/resolution, no transformation needed
        if target_space_id == self.output_space and target_resolution == self.output_resolution:
            return results

        # Get reference affine for target space
        target_key = (target_space_id, target_resolution)
        if target_key not in REFERENCE_AFFINES:
            raise ValueError(
                f"No reference affine available for {target_space_id}@{target_resolution}mm. "
                f"Available spaces: {list(REFERENCE_AFFINES.keys())}"
            )

        target_space = CoordinateSpace(
            identifier=target_space_id,
            resolution=target_resolution,
            reference_affine=REFERENCE_AFFINES[target_key],
        )

        self.logger.info(
            f"Transforming VoxelMap outputs from {self.output_space}@{self.output_resolution}mm "
            f"to {target_space.identifier}@{target_space.resolution}mm"
        )

        transformed_results = {}
        for key, result in results.items():
            # Only transform VoxelMap results
            if isinstance(result, VoxelMap):
                # Auto-detect interpolation method based on data type
                # Use nearest for binary maps (thresholdmaps), linear for continuous
                data = result.data.get_fdata()
                unique_vals = np.unique(data[~np.isnan(data)])
                is_binary = len(unique_vals) <= 2 and set(unique_vals).issubset({0, 1})
                interpolation = "nearest" if is_binary else "linear"

                # Transform the image
                transformed_img = transform_image(
                    img=result.data,
                    source_space=self.output_space,
                    target_space=target_space,
                    source_resolution=int(self.output_resolution),
                    interpolation=interpolation,
                    verbose=self.verbose,
                )

                # Create new VoxelMap with updated space
                transformed_result = VoxelMap(
                    name=result.name,
                    data=transformed_img,
                    space=target_space.identifier,
                    resolution=target_space.resolution,
                    metadata={
                        **result.metadata,
                        "transformed_from": f"{self.output_space}@{self.output_resolution}mm",
                        "transformed_to": f"{target_space.identifier}@{target_space.resolution}mm",
                    },
                )
                transformed_results[key] = transformed_result
            else:
                # Keep non-VoxelMap results as-is
                transformed_results[key] = result

        return transformed_results

__init__(connectome_name, method='boes', pini_percentile=20, n_jobs=1, verbose=False, compute_p_map=True, fdr_alpha=0.05, t_threshold=None, return_in_input_space=True, output_resolution=None, keep_intermediate=False)

Initialize functional network mapping analysis.

Parameters:

Name	Type	Description	Default
connectome_name	str	Name of registered functional connectome (e.g., "GSP1000"). Use list_functional_connectomes() to see available options.	required
method	(boes, pini)	Timeseries extraction method.	"boes"
pini_percentile	int	Percentile threshold for PINI method (0-100).	20
n_jobs	int	Number of parallel jobs for post-processing (result aggregation and spatial resampling). Set to -1 to use all available CPUs.	1
verbose	bool	If True, print progress messages. If False, run silently.	False
compute_p_map	bool	If True, compute p-value map (two-tailed) from t-statistics.	True
fdr_alpha	float	If provided, compute FDR-corrected p-value map using Benjamini-Hochberg procedure at the specified alpha level. Set to None to disable FDR correction. Requires compute_p_map=True.	0.05
t_threshold	float	If provided, create binary mask of voxels with |t| > threshold.	None
return_in_input_space	bool	If True, transform VoxelMap outputs back to the original input mask space. If False, outputs remain in the connectome space (e.g., MNI152NLin6Asym). Requires input SubjectData to have valid space metadata.	True
output_resolution	int	Final output resolution in mm (1 or 2). Controls the resolution of VoxelMap outputs. If None (default), matches the input mask resolution when return_in_input_space=True, or uses analysis resolution when return_in_input_space=False. Set explicitly to ensure consistent output resolution across analyses.	None
keep_intermediate	bool	If True, include intermediate results (e.g., warped mask images) in the output. Useful for debugging and quality control.	False

Raises:

Type	Description
ValueError	If method is not 'boes' or 'pini'.
KeyError	If connectome_name not found in registry.

Source code in src/lacuna/analysis/functional_network_mapping.py

def __init__(
    self,
    connectome_name: str,
    method: str = "boes",
    pini_percentile: int = 20,
    n_jobs: int = 1,
    verbose: bool = False,
    compute_p_map: bool = True,
    fdr_alpha: float | None = 0.05,
    t_threshold: float | None = None,
    return_in_input_space: bool = True,
    output_resolution: int | None = None,
    keep_intermediate: bool = False,
):
    """Initialize functional network mapping analysis.

    Parameters
    ----------
    connectome_name : str
        Name of registered functional connectome (e.g., "GSP1000").
        Use list_functional_connectomes() to see available options.
    method : {"boes", "pini"}, default="boes"
        Timeseries extraction method.
    pini_percentile : int, default=20
        Percentile threshold for PINI method (0-100).
    n_jobs : int, default=1
        Number of parallel jobs for post-processing (result aggregation
        and spatial resampling). Set to -1 to use all available CPUs.
    verbose : bool, default=False
        If True, print progress messages. If False, run silently.
    compute_p_map : bool, default=True
        If True, compute p-value map (two-tailed) from t-statistics.
    fdr_alpha : float, optional, default=0.05
        If provided, compute FDR-corrected p-value map using Benjamini-Hochberg
        procedure at the specified alpha level. Set to None to disable FDR correction.
        Requires compute_p_map=True.
    t_threshold : float, optional
        If provided, create binary mask of voxels with |t| > threshold.
    return_in_input_space : bool, default=True
        If True, transform VoxelMap outputs back to the original input mask space.
        If False, outputs remain in the connectome space (e.g., MNI152NLin6Asym).
        Requires input SubjectData to have valid space metadata.
    output_resolution : int, optional
        Final output resolution in mm (1 or 2). Controls the resolution of VoxelMap outputs.
        If None (default), matches the input mask resolution when return_in_input_space=True,
        or uses analysis resolution when return_in_input_space=False.
        Set explicitly to ensure consistent output resolution across analyses.
    keep_intermediate : bool, default=False
        If True, include intermediate results (e.g., warped mask images)
        in the output. Useful for debugging and quality control.

    Raises
    ------
    ValueError
        If method is not 'boes' or 'pini'.
    KeyError
        If connectome_name not found in registry.
    """
    super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

    # Validate method parameter
    if method not in ("boes", "pini"):
        msg = f"method must be 'boes' or 'pini', got '{method}'"
        raise ValueError(msg)

    # Load connectome from registry
    try:
        connectome = load_functional_connectome(connectome_name)
    except KeyError as e:
        available = [c.name for c in list_functional_connectomes()]
        raise KeyError(
            f"Connectome '{connectome_name}' not found in registry. "
            f"Available connectomes: {', '.join(available)}. "
            f"Use register_functional_connectome() to add new connectomes."
        ) from e

    # Store connectome information
    self.connectome_name = connectome_name
    self.connectome_path = connectome.data_path
    self.output_space = connectome.metadata.space
    self.output_resolution = connectome.metadata.resolution
    self._is_batch_dir = connectome.is_batched

    # Set TARGET_SPACE dynamically from connectome (used by BaseAnalysis._ensure_target_space)
    self.TARGET_SPACE = connectome.metadata.space
    self.TARGET_RESOLUTION = connectome.metadata.resolution

    # Analysis parameters
    self.method = method
    self.pini_percentile = pini_percentile
    if not (1 <= pini_percentile <= 100):
        raise ValueError(f"pini_percentile must be between 1 and 100, got {pini_percentile}")
    self.n_jobs = n_jobs
    self.compute_p_map = compute_p_map
    self.fdr_alpha = fdr_alpha
    if fdr_alpha is not None and not (0 < fdr_alpha <= 1):
        raise ValueError(
            f"fdr_alpha must be between 0 (exclusive) and 1 (inclusive), got {fdr_alpha}"
        )
    self.t_threshold = t_threshold
    self.return_in_input_space = return_in_input_space
    self.final_output_resolution = output_resolution  # User-specified, None means auto

    # Initialize logger
    self.logger = ConsoleLogger(verbose=verbose, width=70)

    # Internal state
    self._batch_files = None
    self._mask_info = None

run_batch(mask_data_list)

Process multiple lesions together using vectorized operations.

This is 10-50x faster than sequential processing because it: 1. Processes all lesions through each connectome batch together 2. Uses vectorized einsum: "lit,itv->liv" for batch correlations 3. Minimizes loop overhead and leverages optimized BLAS operations

This method is automatically called when using VectorizedStrategy for batch processing.

Parameters:

Name	Type	Description	Default
mask_data_list	list[SubjectData]	Batch of lesions to process together	required

Returns:

Type	Description
list[SubjectData]	Processed lesions with results added

Examples:

>>> from lacuna.batch import VectorizedStrategy
>>> from lacuna.analysis import FunctionalNetworkMapping
>>>
>>> analysis = FunctionalNetworkMapping(...)
>>> strategy = VectorizedStrategy()
>>> results = strategy.execute(mask_data_list, analysis)

Source code in src/lacuna/analysis/functional_network_mapping.py

def run_batch(self, mask_data_list: list[SubjectData]) -> list[SubjectData]:
    """Process multiple lesions together using vectorized operations.

    This is 10-50x faster than sequential processing because it:
    1. Processes all lesions through each connectome batch together
    2. Uses vectorized einsum: "lit,itv->liv" for batch correlations
    3. Minimizes loop overhead and leverages optimized BLAS operations

    This method is automatically called when using VectorizedStrategy
    for batch processing.

    Parameters
    ----------
    mask_data_list : list[SubjectData]
        Batch of lesions to process together

    Returns
    -------
    list[SubjectData]
        Processed lesions with results added

    Examples
    --------
    >>> from lacuna.batch import VectorizedStrategy
    >>> from lacuna.analysis import FunctionalNetworkMapping
    >>>
    >>> analysis = FunctionalNetworkMapping(...)
    >>> strategy = VectorizedStrategy()
    >>> results = strategy.execute(mask_data_list, analysis)
    """
    self.logger.info(f"Vectorized batch processing: {len(mask_data_list)} masks")

    # Validate all lesions first
    for mask_data in mask_data_list:
        self._validate_inputs(mask_data)

    # Load mask info once (shared across all batches)
    mask_indices, mask_affine, mask_shape = self._load_mask_info()
    connectome_files = self._get_connectome_files()

    self.logger.success(
        "Loaded connectome metadata",
        details={
            "connectome_batches": len(connectome_files),
            "mask_shape": mask_shape,
            "n_masks": len(mask_data_list),
        },
    )

    # Track empty masks — they get zero-valued results without batch processing
    empty_mask_indices: dict[int, dict] = {}
    for i, mask_data in enumerate(mask_data_list):
        if mask_data.is_empty_mask:
            subject_id = self._format_subject_id(mask_data)
            self.logger.warning(
                f"Empty mask for {subject_id} — will produce zero-valued network maps"
            )
            empty_mask_indices[i] = self._build_empty_mask_results()

    mask_batch = []

    for i, mask_data in enumerate(mask_data_list):
        if i in empty_mask_indices:
            continue
        subject_id = self._format_subject_id(mask_data)

        voxel_indices, _ = self._get_mask_voxel_indices(mask_data)

        if len(voxel_indices) == 0:
            self.logger.warning(
                f"No overlap for {subject_id}: mask has no voxels within connectome "
                f"brain mask after resampling to {self.TARGET_SPACE} "
                f"— will produce zero-valued network maps"
            )
            empty_mask_indices[i] = self._build_empty_mask_results()
            continue

        mask_batch.append(
            {
                "mask_data": mask_data,
                "voxel_indices": voxel_indices,
                "index": i,
            }
        )

    # Process through all connectome batches (VECTORIZED)
    if not mask_batch:
        # All masks were empty or had no overlap — skip batch processing
        self.logger.info("All masks empty or non-overlapping — skipping connectome processing")
        results = [
            mask_data_list[idx].add_result(self.__class__.__name__, empty_result)
            for idx, empty_result in sorted(empty_mask_indices.items())
        ]
        return results

    self.logger.info("Processing connectome batches...")

    # Get number of voxels from first connectome batch
    with h5py.File(connectome_files[0], "r") as hf:
        n_voxels = hf["timeseries"].shape[2]

    # Initialize streaming aggregators for each mask (MEMORY OPTIMIZED)
    # Instead of storing all correlation maps, we accumulate statistics
    aggregators = []
    for _ in range(len(mask_batch)):
        aggregators.append(
            {
                "sum_z": np.zeros(n_voxels, dtype=np.float64),  # Need higher precision for sums
                "sum_z2": np.zeros(n_voxels, dtype=np.float64),
                "n": 0,
            }
        )

    total_subjects = 0
    batch_times = []  # Track timing for each batch

    for batch_idx, connectome_path in enumerate(connectome_files):
        import time

        batch_start_time = time.time()

        with h5py.File(connectome_path, "r") as hf:
            timeseries_data = hf["timeseries"][:]  # (n_subj, n_time, n_vox)
            n_subjects = timeseries_data.shape[0]
            total_subjects += n_subjects

        # Vectorized processing for ALL masks at once
        batch_r_maps = self._compute_batch_correlations_vectorized(
            mask_batch, timeseries_data
        )  # (n_masks, n_subjects, n_voxels)

        # Convert to Fisher z-scores and update running statistics
        # This is the KEY optimization: we don't store full maps!
        with np.errstate(divide="ignore", invalid="ignore"):
            batch_z_maps = np.arctanh(batch_r_maps)
            batch_z_maps = np.nan_to_num(batch_z_maps, nan=0.0, posinf=3.0, neginf=-3.0)

        # Update aggregators with streaming statistics
        for i in range(len(mask_batch)):
            # Sum across subjects in this batch
            aggregators[i]["sum_z"] += np.sum(batch_z_maps[i], axis=0)
            aggregators[i]["sum_z2"] += np.sum(batch_z_maps[i] ** 2, axis=0)
            aggregators[i]["n"] += n_subjects

        # Memory cleanup - immediately free large arrays
        del timeseries_data, batch_r_maps, batch_z_maps

        # Display timing information
        batch_elapsed = time.time() - batch_start_time
        batch_times.append(batch_elapsed)

        # Show progress with time estimates
        if len(batch_times) > 2:
            avg_time = sum(batch_times) / len(batch_times)
            remaining_batches = len(connectome_files) - (batch_idx + 1)
            est_remaining = avg_time * remaining_batches
            self.logger.progress(
                f"Batch completed in {batch_elapsed:.2f}s (est. {est_remaining:.1f}s remaining)",
                current=batch_idx + 1,
                total=len(connectome_files),
            )
        else:
            self.logger.progress(
                f"Batch completed in {batch_elapsed:.2f}s",
                current=batch_idx + 1,
                total=len(connectome_files),
            )

    total_batch_time = sum(batch_times)
    avg_batch_time = total_batch_time / len(batch_times) if batch_times else 0
    self.logger.success(
        "Batch processing complete",
        details={
            "n_batches": len(connectome_files),
            "total_time": f"{total_batch_time:.2f}s",
            "avg_time_per_batch": f"{avg_batch_time:.2f}s",
        },
    )

    # Compute final statistics from aggregated values
    self.logger.info("Aggregating results...")

    # Define per-subject aggregation work
    def _aggregate_one(i, mask_info):
        """Aggregate results for a single mask (thread-safe)."""
        n = aggregators[i]["n"]
        mean_z = aggregators[i]["sum_z"] / n
        mean_r = np.tanh(mean_z).astype(np.float32)

        # Sample standard deviation with Bessel's correction
        var_z_population = (aggregators[i]["sum_z2"] / n) - (mean_z**2)
        var_z_sample = (n / (n - 1)) * var_z_population
        std_z = np.sqrt(np.maximum(var_z_sample, 0))

        mask_copy = mask_info["mask_data"].copy()

        result = self._aggregate_results_from_statistics(
            mask_copy,
            mean_r,
            mean_z,
            std_z,
            mask_indices,
            mask_affine,
            mask_shape,
            total_subjects,
        )
        return mask_info["index"], result

    effective_n_jobs = self.n_jobs
    use_parallel = effective_n_jobs != 1 and len(mask_batch) > 1

    if use_parallel:
        from joblib import Parallel, delayed, parallel_backend

        self.logger.info(
            f"Aggregating {len(mask_batch)} subjects in parallel "
            f"(n_jobs={effective_n_jobs})"
        )
        # inner_max_num_threads=1 prevents fork-after-BLAS-init deadlock on
        # many-core nodes: caps OMP/MKL/OpenBLAS thread pools to 1 in each
        # worker so no inherited mutex can be left locked by a thread that
        # didn't survive fork().
        with parallel_backend("loky", inner_max_num_threads=1):
            pairs = Parallel(n_jobs=effective_n_jobs)(
                delayed(_aggregate_one)(i, mask_info) for i, mask_info in enumerate(mask_batch)
            )
        processed_results = dict(pairs)
    else:
        processed_results = {}
        for i, mask_info in enumerate(mask_batch):
            subject_id = self._format_subject_id(mask_info["mask_data"])
            self.logger.info(f"Aggregating results for: {subject_id}", indent_level=1)
            idx, result = _aggregate_one(i, mask_info)
            processed_results[idx] = result

    # Merge empty-mask results back in at their original indices
    # Wrap raw result dicts into SubjectData to match _aggregate_one output
    for idx, empty_result in empty_mask_indices.items():
        processed_results[idx] = mask_data_list[idx].add_result(
            self.__class__.__name__, empty_result
        )

    # Build final results list in original input order
    results = [processed_results[k] for k in sorted(processed_results)]

    self.logger.success(
        "Batch processing complete",
        details={
            "n_masks_processed": len(processed_results),
        },
    )

    return results

ParcelAggregation

Bases: BaseAnalysis

Atlas aggregation analysis.

Computations performed in input data space (atlases transformed to match input).

Aggregate voxel-level maps to ROI-level statistics using atlases.

This is a composable analysis that can: 1. Compute regional damage from lesion masks (percent overlap, volume) 2. Aggregate connectivity maps from network analyses (mean, sum, etc.) 3. Extract any voxel-level map to atlas ROI statistics

The analysis discovers all atlases in the specified directory and computes the specified aggregation method for each region in each atlas.

Computation Space: Atlases are automatically transformed to match the input data's coordinate space and resolution (parsed from metadata or BIDS-style filenames: tpl-{SPACE}res-{RES}...). If an atlas is already in the input space, no transformation is performed. After transformation, nilearn resamples the atlas to precisely match the input resolution for exact alignment.

Attributes:

Name	Type	Description
TARGET_SPACE	None	Space is determined from the input data. Atlases are transformed to match the input data's coordinate space automatically.
TARGET_RESOLUTION	None	Resolution is determined from the input data. Atlases are transformed to match the input resolution, then nilearn resamples for precise alignment.
batch_strategy	str	Batch processing strategy. Set to "sequential" to avoid race conditions with threading backends when accessing shared atlas resources.

Parameters:

Name	Type	Description	Default
source	str or list[str] or dict[str, str | list[str]]	Source of data to aggregate. Accepts multiple formats: String format: - "maskimg": Use the lesion mask directly - "{AnalysisName}.{result_key}": Use result from previous analysis Example: "FunctionalNetworkMapping.correlation_map" List format: - List of strings in the above formats for multi-source aggregation Example: ["SubjectData.mask_img", "FunctionalNetworkMapping.correlation_map"] Dict format (recommended for multi-source): - Mapping of analysis namespace to result key(s) Example: {"FunctionalNetworkMapping": "rmap"} Example: {"FunctionalNetworkMapping": ["rmap", "zmap"]} Example: {"SubjectData": "maskimg", "FunctionalNetworkMapping": ["rmap", "zmap"]}	"maskimg"
aggregation	str	Aggregation method to use. Options: - "mean": Mean value across ROI voxels - "sum": Sum of values across ROI voxels - "percent": Percentage of ROI voxels that are non-zero (for binary masks) - "volume": Volume (in mm³) of non-zero voxels in ROI - "median": Median value across ROI voxels - "std": Standard deviation across ROI voxels	"mean"
parcel_names	list of str or None	Names of atlases from the registry to process (e.g., "schaefer2018parcels100networks7"). If None, all registered atlases are processed. Use register_parcellation() or register_parcellationes_from_directory() to add custom atlases. If None, all parcellations found in atlas_dir will be processed. Example: ["schaefer2018parcels100networks7", "tian2020parcels16"]	None

Raises:

Type	Description
ValueError	If atlas_dir doesn't exist, is empty, or aggregation method is invalid.
FileNotFoundError	If specified atlas directory doesn't exist.

Notes

• Both 3D and 4D atlases support automatic resampling to match source data spatial resolution via nilearn
• 3D atlases: integer labels, use NiftiLabelsMasker with nearest-neighbor interpolation to preserve labels
• 4D atlases: automatically detect binary (0/1) vs probabilistic (0.0-1.0)
• Binary: use nearest-neighbor interpolation to preserve binary masks
• Probabilistic: use continuous interpolation for probability values
• For 3D atlases: regions defined by integer labels (automatically rounded)
• For 4D atlases: each volume is a binary or probability map for one region
• 4D probabilistic maps are thresholded at threshold parameter if provided
• Results stored in SubjectData.results["ParcelAggregation"] as dict mapping parcellation_name_region_name -> aggregated_value

Examples:

>>> # Use all bundled/registered atlases
>>> analysis = ParcelAggregation(
...     source="maskimg",
...     aggregation="percent"
... )
>>>
>>> # Use specific registered atlases
>>> analysis = ParcelAggregation(
...     source="maskimg",
...     aggregation="percent",
...     parcel_names=["schaefer2018parcels100networks7", "tian2020parcels16"]
... )
>>>
>>> # Register custom atlases first, then use them
>>> from lacuna.assets.parcellations import register_parcellations_from_directory
>>> register_parcellationes_from_directory("/data/my_atlases")
>>> analysis = ParcelAggregation(
...     source="maskimg",
...     aggregation="percent"
... )
>>>
>>> # Average functional connectivity per ROI
>>> analysis = ParcelAggregation(
...     source="FunctionalNetworkMapping.network_map",
...     aggregation="mean"
... )

See Also

RegionalDamage : Convenience wrapper for lesion overlap analysis BaseAnalysis : Parent class defining analysis interface

Source code in src/lacuna/analysis/parcel_aggregation.py

class ParcelAggregation(BaseAnalysis):
    """Atlas aggregation analysis.

    Computations performed in input data space (atlases transformed to match input).

    Aggregate voxel-level maps to ROI-level statistics using atlases.

    This is a composable analysis that can:
    1. Compute regional damage from lesion masks (percent overlap, volume)
    2. Aggregate connectivity maps from network analyses (mean, sum, etc.)
    3. Extract any voxel-level map to atlas ROI statistics

    The analysis discovers all atlases in the specified directory and computes
    the specified aggregation method for each region in each atlas.

    **Computation Space:**
    Atlases are automatically transformed to match the input data's coordinate space
    and resolution (parsed from metadata or BIDS-style filenames: tpl-{SPACE}_res-{RES}_...).
    If an atlas is already in the input space, no transformation is performed.
    After transformation, nilearn resamples the atlas to precisely match the input
    resolution for exact alignment.

    Attributes
    ----------
    TARGET_SPACE : None
        Space is determined from the input data. Atlases are transformed to match
        the input data's coordinate space automatically.
    TARGET_RESOLUTION : None
        Resolution is determined from the input data. Atlases are transformed to
        match the input resolution, then nilearn resamples for precise alignment.
    batch_strategy : str
        Batch processing strategy. Set to "sequential" to avoid race conditions
        with threading backends when accessing shared atlas resources.

    Parameters
    ----------
    source : str or list[str] or dict[str, str | list[str]], default="maskimg"
        Source of data to aggregate. Accepts multiple formats:

        **String format:**
        - "maskimg": Use the lesion mask directly
        - "{AnalysisName}.{result_key}": Use result from previous analysis
          Example: "FunctionalNetworkMapping.correlation_map"

        **List format:**
        - List of strings in the above formats for multi-source aggregation
          Example: ["SubjectData.mask_img", "FunctionalNetworkMapping.correlation_map"]

        **Dict format (recommended for multi-source):**
        - Mapping of analysis namespace to result key(s)
          Example: {"FunctionalNetworkMapping": "rmap"}
          Example: {"FunctionalNetworkMapping": ["rmap", "zmap"]}
          Example: {"SubjectData": "maskimg", "FunctionalNetworkMapping": ["rmap", "zmap"]}

    aggregation : str, default="mean"
        Aggregation method to use. Options:
        - "mean": Mean value across ROI voxels
        - "sum": Sum of values across ROI voxels
        - "percent": Percentage of ROI voxels that are non-zero (for binary masks)
        - "volume": Volume (in mm³) of non-zero voxels in ROI
        - "median": Median value across ROI voxels
        - "std": Standard deviation across ROI voxels
    parcel_names : list of str or None, default=None
        Names of atlases from the registry to process (e.g., "schaefer2018parcels100networks7").
        If None, all registered atlases are processed.
        Use register_parcellation() or register_parcellationes_from_directory() to add custom atlases.
        If None, all parcellations found in atlas_dir will be processed.
        Example: ["schaefer2018parcels100networks7", "tian2020parcels16"]

    Raises
    ------
    ValueError
        If atlas_dir doesn't exist, is empty, or aggregation method is invalid.
    FileNotFoundError
        If specified atlas directory doesn't exist.

    Notes
    -----
    - Both 3D and 4D atlases support automatic resampling to match source data
      spatial resolution via nilearn
    - 3D atlases: integer labels, use NiftiLabelsMasker with nearest-neighbor
      interpolation to preserve labels
    - 4D atlases: automatically detect binary (0/1) vs probabilistic (0.0-1.0)
      * Binary: use nearest-neighbor interpolation to preserve binary masks
      * Probabilistic: use continuous interpolation for probability values
    - For 3D atlases: regions defined by integer labels (automatically rounded)
    - For 4D atlases: each volume is a binary or probability map for one region
    - 4D probabilistic maps are thresholded at `threshold` parameter if provided
    - Results stored in SubjectData.results["ParcelAggregation"] as dict
      mapping parcellation_name_region_name -> aggregated_value

    Examples
    --------
    >>> # Use all bundled/registered atlases
    >>> analysis = ParcelAggregation(
    ...     source="maskimg",
    ...     aggregation="percent"
    ... )
    >>>
    >>> # Use specific registered atlases
    >>> analysis = ParcelAggregation(
    ...     source="maskimg",
    ...     aggregation="percent",
    ...     parcel_names=["schaefer2018parcels100networks7", "tian2020parcels16"]
    ... )
    >>>
    >>> # Register custom atlases first, then use them
    >>> from lacuna.assets.parcellations import register_parcellations_from_directory
    >>> register_parcellationes_from_directory("/data/my_atlases")
    >>> analysis = ParcelAggregation(
    ...     source="maskimg",
    ...     aggregation="percent"
    ... )
    >>>
    >>> # Average functional connectivity per ROI
    >>> analysis = ParcelAggregation(
    ...     source="FunctionalNetworkMapping.network_map",
    ...     aggregation="mean"
    ... )

    See Also
    --------
    RegionalDamage : Convenience wrapper for lesion overlap analysis
    BaseAnalysis : Parent class defining analysis interface
    """

    #: Space is determined from the input data
    TARGET_SPACE = None
    #: Resolution is determined from the input data
    TARGET_RESOLUTION = None
    #: Preferred batch processing strategy (sequential to avoid threading race conditions)
    batch_strategy: str = "sequential"

    VALID_AGGREGATIONS = ["mean", "sum", "percent", "volume", "median", "std"]
    VALID_SOURCES = ["maskimg"]

    def __init__(
        self,
        source: str | list[str] | dict[str, str | list[str]] = "maskimg",
        aggregation: str = "mean",
        parcel_names: list[str] | None = None,
        verbose: bool = False,
        keep_intermediate: bool = False,
    ):
        """Initialize ParcelAggregation analysis.

        Parameters
        ----------
        source : str or list[str] or dict, default="maskimg"
            Source of data to aggregate.
        aggregation : str, default="mean"
            Aggregation method to use.
        parcel_names : list of str or None, default=None
            Names of atlases from the registry to process.
        verbose : bool, default=False
            If True, print progress messages.
        keep_intermediate : bool, default=False
            If True, include intermediate results (e.g., warped mask images)
            in the output. Useful for debugging and quality control.
        """
        super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

        # Initialize logger for warnings and info messages
        from lacuna.utils.logging import ConsoleLogger

        self.logger = ConsoleLogger(verbose=verbose)

        # Normalize and validate source parameter
        self.sources = self._normalize_sources(source)
        self.source = source  # Keep original for compatibility
        self.aggregation = aggregation
        self.parcel_names = parcel_names

        # Validate aggregation method
        if aggregation not in self.VALID_AGGREGATIONS:
            from lacuna.utils.suggestions import format_suggestions, suggest_similar

            suggestions = suggest_similar(aggregation, list(self.VALID_AGGREGATIONS))
            hint = format_suggestions(suggestions)
            msg = (
                f"Invalid aggregation method: '{aggregation}'\n"
                f"Valid options: {', '.join(self.VALID_AGGREGATIONS)}"
            )
            if hint:
                msg = f"{msg}\n{hint}"
            raise ValueError(msg)

        # Threshold validation removed - accepts any float value (T061)
        # This allows for flexible thresholding (e.g., negative z-scores, arbitrary cutoffs)

        # Validate parcel_names if provided
        if parcel_names is not None:
            if not isinstance(parcel_names, list):
                raise TypeError(
                    f"parcel_names must be a list of strings or None, got {type(parcel_names).__name__}"
                )
            if not all(isinstance(name, str) for name in parcel_names):
                raise TypeError("All items in parcel_names must be strings")
            if not parcel_names:
                raise ValueError(
                    "parcel_names cannot be an empty list (use None to process all atlases)"
                )

        # Will be populated in _validate_inputs (thread-safe)
        self.atlases = []
        self._atlases_lock = threading.Lock()
        # Cache for loaded+transformed atlas images, keyed by
        # (atlas_path, input_space, input_resolution) to avoid redundant
        # disk I/O and spatial transformations across subjects
        self._atlas_cache: dict[tuple, nib.Nifti1Image] = {}

    def __getstate__(self):
        """Exclude non-picklable lock from serialization for multiprocessing."""
        state = self.__dict__.copy()
        # Remove the lock - it can't be pickled
        state.pop("_atlases_lock", None)
        # Don't serialize the atlas cache - rebuild in new process
        state.pop("_atlas_cache", None)
        return state

    def __setstate__(self, state):
        """Recreate lock after unpickling for multiprocessing."""
        self.__dict__.update(state)
        # Recreate the lock in the new process
        self._atlases_lock = threading.Lock()
        # Recreate empty atlas cache in the new process
        self._atlas_cache = {}

    def _normalize_sources(self, source: str | list[str] | dict[str, str | list[str]]) -> list[str]:
        """
        Normalize source parameter to a list of sources.

        Parameters
        ----------
        source : str or list[str] or dict[str, str | list[str]]
            Source specification in one of these formats:
            - str: Single source like "maskimg" or "FunctionalNetworkMapping.correlation_map"
            - list[str]: Multiple sources as strings
            - dict: Mapping of namespace to key(s), e.g.,
              {"FunctionalNetworkMapping": "rmap"} or
              {"FunctionalNetworkMapping": ["rmap", "zmap"]}

        Returns
        -------
        list[str]
            Normalized list of source strings in "Namespace.key" format.

        Raises
        ------
        TypeError
            If source is not str, list[str], or dict.
        ValueError
            If source list/dict is empty.

        Examples
        --------
        >>> agg._normalize_sources("maskimg")
        ['mask_img']
        >>> agg._normalize_sources({"FunctionalNetworkMapping": "rmap"})
        ['FunctionalNetworkMapping.rmap']
        >>> agg._normalize_sources({"FunctionalNetworkMapping": ["rmap", "zmap"]})
        ['FunctionalNetworkMapping.rmap', 'FunctionalNetworkMapping.zmap']
        """
        if isinstance(source, str):
            return [source]
        elif isinstance(source, dict):
            if not source:
                raise ValueError("source dict cannot be empty")
            sources = []
            for namespace, keys in source.items():
                if not isinstance(namespace, str):
                    raise TypeError(f"Source namespace must be str, got {type(namespace).__name__}")
                if isinstance(keys, str):
                    # Single key: {"FunctionalNetworkMapping": "rmap"}
                    sources.append(f"{namespace}.{keys}")
                elif isinstance(keys, list):
                    # Multiple keys: {"FunctionalNetworkMapping": ["rmap", "zmap"]}
                    if not keys:
                        raise ValueError(f"Source keys for '{namespace}' cannot be empty")
                    for key in keys:
                        if not isinstance(key, str):
                            raise TypeError(f"Source key must be str, got {type(key).__name__}")
                        sources.append(f"{namespace}.{key}")
                else:
                    raise TypeError(
                        f"Source value must be str or list[str], got {type(keys).__name__}"
                    )
            return sources
        elif isinstance(source, list):
            if not source:
                raise ValueError("source cannot be empty list")
            if not all(isinstance(s, str) for s in source):
                raise TypeError("All items in source list must be strings")
            return source
        else:
            raise TypeError(f"source must be str, list[str], or dict, got {type(source).__name__}")

    def run(
        self, data: "SubjectData | nib.Nifti1Image | list[nib.Nifti1Image]"
    ) -> "SubjectData | ParcelData | list[ParcelData]":
        """
        Execute atlas aggregation analysis on various input types.

        Supports flexible input types with matching return types:
        - SubjectData -> SubjectData (with results attached)
        - nibabel.Nifti1Image -> ParcelData
        - list[nibabel.Nifti1Image] -> list[ParcelData]

        Parameters
        ----------
        data : SubjectData or nibabel.Nifti1Image or list[nibabel.Nifti1Image]
            Input data to aggregate:
            - SubjectData: Standard workflow, returns SubjectData with results
            - nibabel.Nifti1Image: Single image, returns ParcelData
            - list[nibabel.Nifti1Image]: Batch processing, returns list of results

        Returns
        -------
        SubjectData or ParcelData or list[ParcelData]
            Results matching input type:
            - SubjectData input: New SubjectData instance with results in .results dict
            - nibabel input: Single ParcelData
            - list input: List of ParcelData objects (one per input image)

        Raises
        ------
        ValueError
            If input validation fails or source data not found.
        TypeError
            If input type is not supported.

        Notes
        -----
        This method overrides BaseAnalysis.run() to support flexible input types.
        The base class run() is designed for SubjectData only.

        Examples
        --------
        >>> # SubjectData input
        >>> mask_data = SubjectData(mask_img, space='MNI152NLin6Asym', resolution=2)
        >>> analysis = ParcelAggregation(aggregation='percent')
        >>> result = analysis.run(mask_data)
        >>> isinstance(data, SubjectData)
        True

        >>> # Nibabel image input
        >>> import nibabel as nib
        >>> img = nib.load('mask.nii.gz')
        >>> result = analysis.run(img)
        >>> isinstance(result, ParcelData)
        True

        >>> # List of images
        >>> images = [nib.load(f'mask_{i}.nii.gz') for i in range(5)]
        >>> results = analysis.run(images)
        >>> len(results) == 5
        True
        """
        from lacuna.core.data_types import VoxelMap

        # Detect input type and delegate to appropriate handler
        if isinstance(data, SubjectData):
            # Standard SubjectData workflow - use base class run()
            return super().run(data)

        elif isinstance(data, VoxelMap):
            # VoxelMap - run directly without SubjectData wrapper
            return self._run_voxelmap(data)

        elif isinstance(data, nib.Nifti1Image):
            # Single nibabel image - return ParcelData
            return self._run_single_image(data)

        elif isinstance(data, list):
            # List of images or VoxelMaps - return list of results
            if not data:
                raise ValueError("Empty list provided - at least one image required")

            # Check if all are VoxelMaps or all are Images
            if all(isinstance(item, VoxelMap) for item in data):
                # Process VoxelMaps directly
                return [self._run_voxelmap(vm) for vm in data]

            elif all(isinstance(img, nib.Nifti1Image) for img in data):
                return self._run_batch_images(data)

            else:
                raise TypeError(
                    "When providing a list, all items must be of the same type: "
                    "either all VoxelMap or all nibabel.Nifti1Image objects"
                )

        else:
            raise TypeError(
                f"Unsupported input type: {type(data).__name__}\n"
                "Supported types: SubjectData, VoxelMap, nibabel.Nifti1Image, "
                "list[VoxelMap], list[nibabel.Nifti1Image]"
            )

    def _run_single_image(self, img: nib.Nifti1Image) -> "ParcelData":
        """
        Run aggregation on a single nibabel image.

        This method auto-detects space and resolution from the image header,
        then runs aggregation directly without requiring a SubjectData wrapper.
        This allows processing of continuous-valued images (not just binary masks).

        Parameters
        ----------
        img : nibabel.Nifti1Image
            Input image to aggregate

        Returns
        -------
        ParcelData
            Aggregation result combining all atlas aggregations
        """
        # Load atlases using same logic as _run_voxelmap
        if not hasattr(self, "atlases") or not self.atlases:
            self.atlases = self._load_parcellations_from_registry()

        # Auto-detect space and resolution from image header
        input_space = SubjectData._detect_space_from_image(img)
        input_resolution = SubjectData._detect_resolution_from_image(img)

        # Use detected space or fall back with warning
        if input_space is None:
            input_space = "MNI152NLin6Asym"
            self.logger.warning(
                "Could not auto-detect coordinate space from image header. "
                "Assuming MNI152NLin6Asym. For explicit control, use SubjectData wrapper: "
                "SubjectData(img, space='...', resolution=...)"
            )
        if input_resolution is None:
            input_resolution = float(round(abs(img.affine[0, 0])))
            self.logger.info(
                f"Could not detect resolution from image header, using voxel size: {input_resolution}mm"
            )

        # Calculate voxel volume from source data
        voxel_volume_mm3 = np.abs(np.linalg.det(img.affine[:3, :3]))

        # Collect all ROI results across atlases
        all_roi_data = {}

        # Process each atlas
        for atlas_info in self.atlases:
            parcellation_name = atlas_info["name"]
            atlas_space = atlas_info.get("space")
            atlas_resolution = atlas_info.get("resolution")

            # Load atlas image
            atlas_img = nib.load(atlas_info["atlas_path"])

            # Transform atlas to match input data space if needed
            atlas_img = self._ensure_atlas_matches_input_space(
                atlas_img=atlas_img,
                atlas_space=atlas_space,
                atlas_resolution=atlas_resolution,
                input_space=input_space,
                input_resolution=input_resolution,
                input_affine=img.affine,
                parcellation_name=parcellation_name,
            )

            labels = atlas_info["labels"]
            atlas_data = atlas_img.get_fdata()

            if atlas_data.ndim == 3:
                # 3D integer-labeled atlas
                atlas_results = self._aggregate_3d_atlas(img, atlas_img, labels, voxel_volume_mm3)
            elif atlas_data.ndim == 4:
                # 4D probabilistic atlas
                atlas_results = self._aggregate_4d_atlas(img, atlas_img, labels, voxel_volume_mm3)
            else:
                continue

            # Merge results from this atlas
            all_roi_data.update(atlas_results)

        # Return single ParcelData with all ROI results
        from lacuna.core.data_types import ParcelData

        return ParcelData(
            name=f"{self.aggregation}_aggregation",
            data=all_roi_data,
            parcel_names=(
                self.parcel_names if self.parcel_names else [a["name"] for a in self.atlases]
            ),
            aggregation_method=self.aggregation,
            metadata={
                "source": "Nifti1Image",
                "n_regions": len(all_roi_data),
                "space": input_space,
                "resolution": input_resolution,
            },
        )

    def _run_batch_images(self, images: list[nib.Nifti1Image]) -> list["ParcelData"]:
        """
        Run aggregation on a batch of nibabel images.

        Parameters
        ----------
        images : list[nibabel.Nifti1Image]
            List of images to aggregate

        Returns
        -------
        list[ParcelData]
            List of aggregation results (one per input image)
        """
        results = []
        for img in images:
            result = self._run_single_image(img)
            results.append(result)

        return results

    def _run_voxelmap(self, voxel_map: "VoxelMap") -> "ParcelData":
        """
        Run aggregation on a VoxelMap directly.

        This bypasses SubjectData validation since VoxelMaps can contain
        continuous values (e.g., correlation maps, z-scores).

        Parameters
        ----------
        voxel_map : VoxelMap
            VoxelMap containing the data to aggregate

        Returns
        -------
        ParcelData
            Aggregation result combining all atlas aggregations
        """
        # Load atlases using same logic as _load_parcellations_from_registry
        if not hasattr(self, "atlases") or not self.atlases:
            self.atlases = self._load_parcellations_from_registry()

        # Get space and resolution from VoxelMap
        input_space = voxel_map.space
        input_resolution = voxel_map.resolution
        source_img = voxel_map.data

        # Calculate voxel volume from source data
        voxel_volume_mm3 = np.abs(np.linalg.det(source_img.affine[:3, :3]))

        # Collect all ROI results across atlases
        all_roi_data = {}

        # Process each atlas
        for atlas_info in self.atlases:
            parcellation_name = atlas_info["name"]
            atlas_space = atlas_info.get("space")
            atlas_resolution = atlas_info.get("resolution")

            # Load atlas image
            atlas_img = nib.load(atlas_info["atlas_path"])

            # Transform atlas to match input data space if needed
            atlas_img = self._ensure_atlas_matches_input_space(
                atlas_img=atlas_img,
                atlas_space=atlas_space,
                atlas_resolution=atlas_resolution,
                input_space=input_space,
                input_resolution=input_resolution,
                input_affine=source_img.affine,
                parcellation_name=parcellation_name,
            )

            labels = atlas_info["labels"]
            atlas_data = atlas_img.get_fdata()

            if atlas_data.ndim == 3:
                # 3D integer-labeled atlas
                atlas_results = self._aggregate_3d_atlas(
                    source_img, atlas_img, labels, voxel_volume_mm3
                )
            elif atlas_data.ndim == 4:
                # 4D probabilistic atlas
                atlas_results = self._aggregate_4d_atlas(
                    source_img, atlas_img, labels, voxel_volume_mm3
                )
            else:
                continue

            # Merge results from this atlas
            all_roi_data.update(atlas_results)

        # Return single ParcelData with all ROI results
        from lacuna.core.data_types import ParcelData

        return ParcelData(
            name=f"{self.aggregation}_aggregation",
            data=all_roi_data,
            parcel_names=(
                self.parcel_names if self.parcel_names else [a["name"] for a in self.atlases]
            ),
            aggregation_method=self.aggregation,
            metadata={
                "source": "VoxelMap",
                "source_name": voxel_map.name,
                "n_regions": len(all_roi_data),
                "space": input_space,
                "resolution": input_resolution,
            },
        )

    def _validate_inputs(self, mask_data: SubjectData) -> None:
        """
        Validate lesion data and load atlases from registry.

        Parameters
        ----------
        mask_data : SubjectData
            Lesion data to validate

        Raises
        ------
        ValueError
            If lesion data is invalid or source data not found
        """
        # Build list of available sources
        available = ["SubjectData.mask_img"]
        if mask_data.results:
            for analysis_name, analysis_results in mask_data.results.items():
                for key in analysis_results.keys():
                    available.append(f"{analysis_name}.{key}")

        # Validate each source exists
        missing_sources = []
        for src in self.sources:
            source_img = self._get_source_image_for_source(mask_data, src)
            if source_img is None:
                missing_sources.append(src)

        if missing_sources:
            from lacuna.utils.suggestions import format_suggestions, suggest_similar

            suggestions = []
            for missing in missing_sources:
                similar = suggest_similar(missing, available)
                if similar:
                    suggestions.extend(similar)

            error_msg = (
                f"Source data not found: {missing_sources}\n"
                "Check that the source exists in SubjectData.\n"
                f"Available sources: {', '.join(available)}"
            )
            if suggestions:
                error_msg += f"\n\nDid you mean: {format_suggestions(suggestions)}?"

            raise ValueError(error_msg)

        # Load atlases from registry (thread-safe)
        # Use lock to prevent race condition where multiple threads
        # simultaneously check 'if not self.atlases' and all try to load
        with self._atlases_lock:
            if not self.atlases:
                self.atlases = self._load_parcellations_from_registry()

        if not self.atlases:
            if self.parcel_names is not None:
                raise ValueError(
                    f"No matching parcellations found for specified names: {self.parcel_names}\n"
                    "Available parcellations in registry: check list_parcellations()\n"
                    "Use register_parcellation() or register_parcellationes_from_directory() to add atlases"
                )
            else:
                raise ValueError(
                    "No valid parcellations found in registry\n"
                    "Use register_parcellation() or register_parcellationes_from_directory() to add atlases"
                )

        # Warn if some requested atlases weren't found
        if self.parcel_names is not None:
            found_names = {atlas["name"] for atlas in self.atlases}
            missing_names = set(self.parcel_names) - found_names
            if missing_names:
                self.logger.warning(
                    f"Some requested parcellations were not found: {sorted(missing_names)}. "
                    f"Found: {sorted(found_names)}"
                )

    def _load_parcellations_from_registry(self) -> list[dict]:
        """
        Load atlases from the registry (bundled or user-registered).

        Returns
        -------
        list[dict]
            List of atlas dictionaries with keys: name, image, labels, space, resolution
        """
        from lacuna.assets.parcellations.loader import BUNDLED_PARCELLATIONS_DIR

        # Get atlases from registry (filter by names if provided)
        if self.parcel_names is not None:
            # Load specific atlases by name
            atlases_data = []
            for name in self.parcel_names:
                try:
                    atlas = load_parcellation(name)

                    # Resolve paths (absolute or relative to bundled dir)
                    atlas_filename_path = Path(atlas.metadata.parcellation_filename)
                    if atlas_filename_path.is_absolute():
                        atlas_path = atlas_filename_path
                    else:
                        atlas_path = (
                            BUNDLED_PARCELLATIONS_DIR / atlas.metadata.parcellation_filename
                        )

                    labels_filename_path = Path(atlas.metadata.labels_filename)
                    if labels_filename_path.is_absolute():
                        labels_path = labels_filename_path
                    else:
                        labels_path = BUNDLED_PARCELLATIONS_DIR / atlas.metadata.labels_filename

                    atlases_data.append(
                        {
                            "name": name,
                            "atlas_path": atlas_path,
                            "labels_path": labels_path,
                            "labels": atlas.labels,
                            "space": atlas.metadata.space,
                            "resolution": atlas.metadata.resolution,
                            "is_4d": getattr(atlas.metadata, "is_4d", False),
                        }
                    )
                except KeyError:
                    # Atlas not in registry - will be caught by validation
                    pass
        else:
            # Load all registered atlases
            atlas_metadatas = list_parcellations()
            atlases_data = []
            for metadata in atlas_metadatas:
                atlas = load_parcellation(metadata.name)

                # Resolve paths (absolute or relative to bundled dir)
                atlas_filename_path = Path(atlas.metadata.parcellation_filename)
                if atlas_filename_path.is_absolute():
                    atlas_path = atlas_filename_path
                else:
                    atlas_path = BUNDLED_PARCELLATIONS_DIR / atlas.metadata.parcellation_filename

                labels_filename_path = Path(atlas.metadata.labels_filename)
                if labels_filename_path.is_absolute():
                    labels_path = labels_filename_path
                else:
                    labels_path = BUNDLED_PARCELLATIONS_DIR / atlas.metadata.labels_filename

                atlases_data.append(
                    {
                        "name": metadata.name,
                        "atlas_path": atlas_path,
                        "labels_path": labels_path,
                        "labels": atlas.labels,
                        "space": metadata.space,
                        "resolution": metadata.resolution,
                        "is_4d": getattr(metadata, "is_4d", False),
                    }
                )

        return atlases_data

    def _ensure_atlas_matches_input_space(
        self,
        atlas_img: nib.Nifti1Image,
        atlas_space: str,
        atlas_resolution: int,
        input_space: str,
        input_resolution: int,
        input_affine: np.ndarray,
        parcellation_name: str | None = None,
    ) -> nib.Nifti1Image:
        """
        Transform atlas to match input data space if spaces don't match.

        This allows ParcelAggregation to work with any voxel-level image,
        not just lesion data, by transforming the atlas to the input space.

        Parameters
        ----------
        atlas_img : nib.Nifti1Image
            Atlas image to potentially transform
        atlas_space : str
            Atlas coordinate space (e.g., 'MNI152NLin6Asym')
        atlas_resolution : int
            Atlas resolution in mm (e.g., 1 or 2)
        input_space : str
            Input data coordinate space
        input_resolution : int
            Input data resolution in mm
        input_affine : np.ndarray
            Input data affine matrix

        Returns
        -------
        nib.Nifti1Image
            Atlas in input space (transformed if needed, original if already matching)
        """
        # If atlas doesn't specify space, assume it matches
        if atlas_space is None:
            return atlas_img

        # Validate declared space against image header (affine and shape)
        from lacuna.core.spaces import (
            REFERENCE_SHAPES,
            detect_space_from_header,
            spaces_are_equivalent,
        )

        detected = detect_space_from_header(atlas_img)
        if detected is None:
            # Affine check failed (e.g. flipped data strides) — fall back to
            # shape + voxel-size matching, which is orientation-independent.
            img_shape = atlas_img.shape[:3]
            voxel_size = round(float(atlas_img.header.get_zooms()[0]), 1)
            shape_to_space = {
                shape: space
                for (space, res), shape in REFERENCE_SHAPES.items()
                if res == voxel_size
            }
            if img_shape in shape_to_space:
                detected = (shape_to_space[img_shape], voxel_size)

        if detected is not None:
            detected_space, _ = detected
            if not spaces_are_equivalent(detected_space, atlas_space):
                raise ValueError(
                    f"Parcellation '{parcellation_name}': declared space is "
                    f"'{atlas_space}' but image header matches "
                    f"'{detected_space}'. Check that the correct coordinate "
                    f"space was specified for this atlas."
                )

        if spaces_are_equivalent(atlas_space, input_space):
            # Same space or equivalent alias - no coordinate transformation needed
            # (nilearn will handle resolution resampling during aggregation)
            return atlas_img

        # Need to transform atlas to input space
        from lacuna.core.spaces import CoordinateSpace
        from lacuna.spatial.transform import transform_image

        # Create target space matching input data
        target_space = CoordinateSpace(
            identifier=input_space,
            resolution=input_resolution,
            reference_affine=input_affine,
        )

        # Transform atlas using nearest neighbor to preserve labels
        # Logging is handled by transform_image
        return transform_image(
            img=atlas_img,
            source_space=atlas_space,
            target_space=target_space,
            source_resolution=atlas_resolution,
            interpolation="nearest",  # Preserve integer labels
            image_name=f"atlas '{parcellation_name}'" if parcellation_name else "atlas",
            verbose=self.verbose,
        )

    def _run_analysis(self, mask_data: SubjectData) -> dict[str, "DataContainer"]:
        """
        Compute ROI-level aggregation for all atlases and sources.

        Parameters
        ----------
        mask_data : SubjectData
            Validated lesion data

        Returns
        -------
        dict[str, DataContainer]
            Dictionary mapping BIDS-style keys to ParcelData objects.
            Keys follow the pattern: parc-{atlas}_source-{SourceClass}_desc-{key}
        """
        # Log analysis start with atlas names
        n_atlases = len(self.atlases) if hasattr(self, "atlases") and self.atlases else 0
        n_sources = len(self.sources)
        atlas_names = [a["name"] for a in self.atlases] if self.atlases else []
        self.logger.info(
            f"Aggregating {n_sources} source(s) across {n_atlases} atlas(es): {', '.join(atlas_names)}"
        )

        # Get input data space/resolution once
        input_space = mask_data.space
        input_resolution = mask_data.resolution

        # Collect results with BIDS-style keys
        all_results: dict[str, DataContainer] = {}

        # Process each source
        for source in self.sources:
            # Parse source string to extract source class and key
            if "." in source:
                # Cross-analysis source: "AnalysisName.result_key"
                source_class, source_key = source.split(".", 1)
            else:
                # Direct source: "maskimg" -> from SubjectData
                source_class = "SubjectData"
                source_key = source

            # Get source image for this source
            source_img = self._get_source_image_for_source(mask_data, source)

            # Calculate voxel volume from source data
            voxel_volume_mm3 = np.abs(np.linalg.det(source_img.affine[:3, :3]))

            # Process each atlas
            for atlas_info in self.atlases:
                parcellation_name = atlas_info["name"]
                atlas_space = atlas_info.get("space")
                atlas_resolution = atlas_info.get("resolution")

                # Load and transform atlas (cached across subjects)
                cache_key = (
                    atlas_info["atlas_path"],
                    input_space,
                    input_resolution,
                )
                if cache_key in self._atlas_cache:
                    atlas_img = self._atlas_cache[cache_key]
                else:
                    atlas_img = nib.load(atlas_info["atlas_path"])
                    atlas_img = self._ensure_atlas_matches_input_space(
                        atlas_img=atlas_img,
                        atlas_space=atlas_space,
                        atlas_resolution=atlas_resolution,
                        input_space=input_space,
                        input_resolution=input_resolution,
                        input_affine=source_img.affine,
                        parcellation_name=parcellation_name,
                    )
                    self._atlas_cache[cache_key] = atlas_img

                # Store warped atlas as intermediate if requested
                if self.keep_intermediate:
                    from lacuna.core.data_types import VoxelMap

                    # Build unique key for this atlas + source combination
                    intermediate_key = f"warped_atlas_{parcellation_name}_{source_key}"
                    warped_atlas = VoxelMap(
                        name=f"warped_{parcellation_name}",
                        data=atlas_img,
                        space=input_space,
                        resolution=input_resolution,
                        metadata={
                            "original_space": atlas_space,
                            "original_resolution": atlas_resolution,
                            "parcellation_name": parcellation_name,
                            "source": source,
                            "description": (
                                f"Atlas '{parcellation_name}' transformed from "
                                f"{atlas_space}@{atlas_resolution}mm to "
                                f"{input_space}@{input_resolution}mm"
                            ),
                        },
                    )
                    all_results[intermediate_key] = warped_atlas

                labels = atlas_info["labels"]
                atlas_data = atlas_img.get_fdata()

                # Warn if nilearn will resample atlas to match source resolution
                atlas_shape = atlas_data.shape[:3]  # Handle 4D atlases
                source_shape = source_img.get_fdata().shape
                if source_shape != atlas_shape:
                    self.logger.info(
                        f"Resampling parcellation '{parcellation_name}' to match source data "
                        f"(source: {source_shape}, parcellation: {atlas_shape})"
                    )

                if atlas_data.ndim == 3:
                    # 3D integer-labeled atlas - use nilearn NiftiLabelsMasker
                    atlas_results = self._aggregate_3d_atlas(
                        source_img, atlas_img, labels, voxel_volume_mm3
                    )
                elif atlas_data.ndim == 4:
                    # 4D probabilistic atlas - use nilearn resampling
                    atlas_results = self._aggregate_4d_atlas(
                        source_img, atlas_img, labels, voxel_volume_mm3
                    )
                else:
                    self.logger.warning(
                        f"Skipping parcellation '{parcellation_name}': "
                        f"unexpected dimensions {atlas_data.ndim}D"
                    )
                    continue

                # Create ParcelData for this atlas + source combination
                roi_result = ParcelData(
                    name=parcellation_name,
                    data=atlas_results,
                    parcel_names=[parcellation_name],
                    aggregation_method=self.aggregation,
                    metadata={
                        "source": source,
                        "source_class": source_class,
                        "source_key": source_key,
                        "n_regions": len(atlas_results),
                    },
                )

                # Build BIDS-style result key
                result_key = build_result_key(
                    atlas=parcellation_name,
                    source=source_class,
                    desc=source_key,
                )

                all_results[result_key] = roi_result

        self.logger.success(f"Aggregation complete ({len(all_results)} results)")
        return all_results

    def _aggregate_3d_atlas(
        self,
        source_img: nib.Nifti1Image,
        atlas_img: nib.Nifti1Image,
        labels: dict[int, str],
        voxel_volume_mm3: float,
    ) -> dict[str, float]:
        """
        Aggregate source data for 3D integer-labeled atlas using nilearn.

        Uses nilearn's NiftiLabelsMasker for robust extraction with automatic
        resampling, masking, and efficient computation.

        Note: Suppresses nilearn's verbose label removal warnings when verbose is False.

        Parameters
        ----------
        source_img : nib.Nifti1Image
            Source image to aggregate
        atlas_img : nib.Nifti1Image
            3D atlas with integer labels
        labels : dict[int, str]
            Mapping from region ID to region name
        voxel_volume_mm3 : float
            Volume of one voxel in mm³ (for volume aggregation)

        Returns
        -------
        dict[str, float]
            Mapping from region name to aggregated value
        """
        import warnings

        # Suppress nilearn's verbose label removal warnings unless in verbose mode
        # These warnings come from sklearn's set_output and are too verbose for standard use
        if not self.verbose:
            warnings.filterwarnings(
                "ignore",
                message=".*following labels were removed.*",
                category=UserWarning,
                module="sklearn",
            )

        # Map our aggregation methods to nilearn strategies
        strategy_map = {
            "mean": "mean",
            "sum": "sum",
            "median": "median",
            "std": "standard_deviation",
            "percent": "mean",  # Will multiply by 100
            "volume": "sum",  # Will multiply by voxel_volume_mm3
        }

        if self.aggregation not in strategy_map:
            raise ValueError(f"Unknown aggregation method: {self.aggregation}")

        strategy = strategy_map[self.aggregation]

        # Create label names list (NiftiLabelsMasker expects ordered list)
        # Background (0) should not be included
        atlas_data = atlas_img.get_fdata()

        # Round atlas values to ensure integer labels
        # This handles edge cases where resampling or data type conversion
        # might introduce small floating point values
        atlas_data_rounded = np.round(atlas_data).astype(int)

        region_ids = np.unique(atlas_data_rounded)
        region_ids = region_ids[
            region_ids > 0
        ]  # Exclude background        # Build ordered list of label names
        label_names = [labels.get(int(rid), f"Region{int(rid)}") for rid in sorted(region_ids)]

        # Initialize NiftiLabelsMasker with appropriate settings
        masker = NiftiLabelsMasker(
            labels_img=atlas_img,
            labels=label_names,
            background_label=0,
            strategy=strategy,
            resampling_target="data",  # Resample atlas to match source data
            standardize=False,  # Don't normalize for static maps
            detrend=False,  # No detrending for static maps
            memory=None,  # No caching for now
            verbose=0,
            keep_masked_labels=False,  # Remove empty region signals (future nilearn default)
        )

        # Extract values - nilearn expects 4D input (add time dimension if needed)
        if source_img.ndim == 3:
            # Add a dummy 4th dimension for time
            source_data_4d = source_img.get_fdata()[..., np.newaxis]
            source_img_4d = nib.Nifti1Image(source_data_4d, source_img.affine)
        else:
            source_img_4d = source_img

        # Transform: returns (n_timepoints, n_regions) array
        region_values = masker.fit_transform(source_img_4d)

        # Squeeze to get (n_regions,) for single timepoint
        if region_values.shape[0] == 1:
            region_values = region_values.squeeze(axis=0)

        # Apply post-processing based on aggregation type
        if self.aggregation == "percent":
            # Convert mean (0-1) to percentage (0-100)
            region_values = region_values * 100
        elif self.aggregation == "volume":
            # Convert count to volume (mm³)
            region_values = region_values * voxel_volume_mm3

        # Build results dict
        # Note: region_values length might not match label_names if regions are lost during resampling
        # We zip without strict=True and handle the mismatch
        results = {}
        for i, value in enumerate(region_values):
            if i < len(label_names):
                label_name = label_names[i]
            else:
                # Fallback if we get more regions than expected
                label_name = f"Region{i}"
            results[label_name] = float(value)

        return results

    def _aggregate_4d_atlas(
        self,
        source_img: nib.Nifti1Image,
        atlas_img: nib.Nifti1Image,
        labels: dict[int, str],
        voxel_volume_mm3: float,
    ) -> dict[str, float]:
        """
        Aggregate source data for 4D atlas with automatic resampling.

        Uses nilearn's resample_to_img for automatic spatial alignment of atlas
        to source data. Detects binary vs probabilistic atlases and uses appropriate
        interpolation method ('nearest' for binary, 'continuous' for probabilistic).

        Parameters
        ----------
        source_img : nib.Nifti1Image
            Source image to aggregate
        atlas_img : nib.Nifti1Image
            4D atlas (x, y, z, n_regions) with binary or probability maps
        labels : dict[int, str]
            Mapping from region ID to region name
        voxel_volume_mm3 : float
            Volume of one voxel in mm³ (for volume aggregation)

        Returns
        -------
        dict[str, float]
            Mapping from region name to aggregated value
        """
        # Detect if atlas is binary (only 0s and 1s) or probabilistic
        atlas_data_orig = atlas_img.get_fdata()
        unique_values = np.unique(atlas_data_orig)
        is_binary = np.all(np.isin(unique_values, [0.0, 1.0]))

        # Use appropriate interpolation based on atlas type
        # 'nearest' for binary to preserve 0/1 values
        # 'continuous' for probabilistic to interpolate between probability values
        interpolation = "nearest" if is_binary else "continuous"

        # Resample atlas to match source data spatial resolution
        atlas_resampled = resample_to_img(
            atlas_img,
            source_img,
            interpolation=interpolation,
            copy=True,
            force_resample=True,
            copy_header=True,
        )

        source_data = source_img.get_fdata()
        atlas_data = atlas_resampled.get_fdata()

        results = {}
        n_regions = atlas_data.shape[3]

        # Get sorted label IDs to map volume indices to label IDs
        # Volume index i corresponds to the i-th label ID in sorted order
        sorted_label_ids = sorted(labels.keys())

        # Validate that we have the right number of labels
        if len(sorted_label_ids) != n_regions:
            self.logger.warning(
                f"Number of volumes ({n_regions}) does not match number of labels "
                f"({len(sorted_label_ids)}). Using available labels."
            )

        for region_idx in range(n_regions):
            # Get probability map for this region
            prob_map = atlas_data[:, :, :, region_idx]

            # Create binary mask from non-zero probability values
            region_mask = prob_map > 0

            # Get values in this region
            region_values = source_data[region_mask]

            # Compute aggregation
            value = self._compute_aggregation(region_values, region_mask, voxel_volume_mm3)

            # Map volume index to label ID using sorted label IDs
            # Volume 0 → sorted_label_ids[0] (could be 0, 1, or any starting ID)
            # Volume 1 → sorted_label_ids[1], etc.
            if region_idx < len(sorted_label_ids):
                region_id = sorted_label_ids[region_idx]
                region_name = labels[region_id]
            else:
                # Fallback if more volumes than labels
                region_name = f"Region{region_idx}"

            results[region_name] = value

        return results

    def _compute_aggregation(
        self,
        region_values: np.ndarray,
        region_mask: np.ndarray,
        voxel_volume_mm3: float,
    ) -> float:
        """
        Compute specified aggregation method on region values.

        Parameters
        ----------
        region_values : np.ndarray
            Values within the region
        region_mask : np.ndarray
            Boolean mask for the region
        voxel_volume_mm3 : float
            Volume of one voxel in mm³

        Returns
        -------
        float
            Aggregated value
        """
        if len(region_values) == 0:
            return 0.0

        if self.aggregation == "mean":
            return float(np.mean(region_values))

        elif self.aggregation == "sum":
            return float(np.sum(region_values))

        elif self.aggregation == "median":
            return float(np.median(region_values))

        elif self.aggregation == "std":
            return float(np.std(region_values))

        elif self.aggregation == "percent":
            # Percentage of ROI voxels that are non-zero
            n_total = np.sum(region_mask)
            n_nonzero = np.sum(region_values > 0)
            return (n_nonzero / n_total * 100) if n_total > 0 else 0.0

        elif self.aggregation == "volume":
            # Volume of non-zero voxels in mm³
            n_nonzero = np.sum(region_values > 0)
            return n_nonzero * voxel_volume_mm3

        else:
            raise ValueError(f"Unknown aggregation method: {self.aggregation}")

    def _get_source_image_for_source(
        self, mask_data: SubjectData, source: str
    ) -> nib.Nifti1Image | None:
        """
        Get source image from SubjectData for a specific source string.

        Parameters
        ----------
        mask_data : SubjectData
            Lesion data containing source
        source : str
            Source specification (e.g., "SubjectData.mask_img", "FunctionalNetworkMapping.correlation_map")

        Returns
        -------
        nib.Nifti1Image or None
            Source image, or None if not found
        """
        # Handle "SubjectData.mask_img" or just "maskimg"
        if source == "maskimg" or source == "SubjectData.mask_img":
            return mask_data.mask_img

        # Result from previous analysis: "AnalysisName.result_key"
        if "." in source:
            analysis_name, result_key = source.split(".", 1)

            # Handle SubjectData prefix
            if analysis_name == "SubjectData":
                if result_key == "maskimg":
                    return mask_data.mask_img
                return None

            if analysis_name in mask_data.results:
                analysis_results = mask_data.results[analysis_name]

                if result_key in analysis_results:
                    result = analysis_results[result_key]

                    # If it's a NIfTI image, return it
                    if isinstance(result, nib.Nifti1Image):
                        return result

                    # If it's a VoxelMap, return the underlying image
                    from lacuna.core.data_types import VoxelMap

                    if isinstance(result, VoxelMap):
                        return result.data

                    # If it's a path, load it
                    if isinstance(result, (str, Path)):
                        result_path = Path(result)
                        if result_path.exists():
                            return nib.load(result_path)

        return None

    def _get_source_image(self, mask_data: SubjectData) -> nib.Nifti1Image | None:
        """
        Get source image from SubjectData based on first source in sources list.

        This is a compatibility method for single-source usage.

        Parameters
        ----------
        mask_data : SubjectData
            Lesion data containing source

        Returns
        -------
        nib.Nifti1Image or None
            Source image, or None if not found
        """
        if self.sources:
            return self._get_source_image_for_source(mask_data, self.sources[0])
        return None

    def _get_parameters(self) -> dict:
        """Get analysis parameters for provenance and display.

        Returns
        -------
        dict
            Dictionary of parameter names and values.
        """
        return {
            "source": self.source,
            "aggregation": self.aggregation,
            "parcel_names": self.parcel_names,
            "num_atlases": len(self.atlases) if hasattr(self, "atlases") else None,
            "verbose": self.verbose,
        }

__getstate__()

Exclude non-picklable lock from serialization for multiprocessing.

Source code in src/lacuna/analysis/parcel_aggregation.py

def __getstate__(self):
    """Exclude non-picklable lock from serialization for multiprocessing."""
    state = self.__dict__.copy()
    # Remove the lock - it can't be pickled
    state.pop("_atlases_lock", None)
    # Don't serialize the atlas cache - rebuild in new process
    state.pop("_atlas_cache", None)
    return state

__init__(source='maskimg', aggregation='mean', parcel_names=None, verbose=False, keep_intermediate=False)

Initialize ParcelAggregation analysis.

Parameters:

Name	Type	Description	Default
source	str or list[str] or dict	Source of data to aggregate.	"maskimg"
aggregation	str	Aggregation method to use.	"mean"
parcel_names	list of str or None	Names of atlases from the registry to process.	None
verbose	bool	If True, print progress messages.	False
keep_intermediate	bool	If True, include intermediate results (e.g., warped mask images) in the output. Useful for debugging and quality control.	False

Source code in src/lacuna/analysis/parcel_aggregation.py

def __init__(
    self,
    source: str | list[str] | dict[str, str | list[str]] = "maskimg",
    aggregation: str = "mean",
    parcel_names: list[str] | None = None,
    verbose: bool = False,
    keep_intermediate: bool = False,
):
    """Initialize ParcelAggregation analysis.

    Parameters
    ----------
    source : str or list[str] or dict, default="maskimg"
        Source of data to aggregate.
    aggregation : str, default="mean"
        Aggregation method to use.
    parcel_names : list of str or None, default=None
        Names of atlases from the registry to process.
    verbose : bool, default=False
        If True, print progress messages.
    keep_intermediate : bool, default=False
        If True, include intermediate results (e.g., warped mask images)
        in the output. Useful for debugging and quality control.
    """
    super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

    # Initialize logger for warnings and info messages
    from lacuna.utils.logging import ConsoleLogger

    self.logger = ConsoleLogger(verbose=verbose)

    # Normalize and validate source parameter
    self.sources = self._normalize_sources(source)
    self.source = source  # Keep original for compatibility
    self.aggregation = aggregation
    self.parcel_names = parcel_names

    # Validate aggregation method
    if aggregation not in self.VALID_AGGREGATIONS:
        from lacuna.utils.suggestions import format_suggestions, suggest_similar

        suggestions = suggest_similar(aggregation, list(self.VALID_AGGREGATIONS))
        hint = format_suggestions(suggestions)
        msg = (
            f"Invalid aggregation method: '{aggregation}'\n"
            f"Valid options: {', '.join(self.VALID_AGGREGATIONS)}"
        )
        if hint:
            msg = f"{msg}\n{hint}"
        raise ValueError(msg)

    # Threshold validation removed - accepts any float value (T061)
    # This allows for flexible thresholding (e.g., negative z-scores, arbitrary cutoffs)

    # Validate parcel_names if provided
    if parcel_names is not None:
        if not isinstance(parcel_names, list):
            raise TypeError(
                f"parcel_names must be a list of strings or None, got {type(parcel_names).__name__}"
            )
        if not all(isinstance(name, str) for name in parcel_names):
            raise TypeError("All items in parcel_names must be strings")
        if not parcel_names:
            raise ValueError(
                "parcel_names cannot be an empty list (use None to process all atlases)"
            )

    # Will be populated in _validate_inputs (thread-safe)
    self.atlases = []
    self._atlases_lock = threading.Lock()
    # Cache for loaded+transformed atlas images, keyed by
    # (atlas_path, input_space, input_resolution) to avoid redundant
    # disk I/O and spatial transformations across subjects
    self._atlas_cache: dict[tuple, nib.Nifti1Image] = {}

__setstate__(state)

Recreate lock after unpickling for multiprocessing.

Source code in src/lacuna/analysis/parcel_aggregation.py

def __setstate__(self, state):
    """Recreate lock after unpickling for multiprocessing."""
    self.__dict__.update(state)
    # Recreate the lock in the new process
    self._atlases_lock = threading.Lock()
    # Recreate empty atlas cache in the new process
    self._atlas_cache = {}

run(data)

Execute atlas aggregation analysis on various input types.

Supports flexible input types with matching return types: - SubjectData -> SubjectData (with results attached) - nibabel.Nifti1Image -> ParcelData - list[nibabel.Nifti1Image] -> list[ParcelData]

Parameters:

Name	Type	Description	Default
data	SubjectData or Nifti1Image or list[Nifti1Image]	Input data to aggregate: - SubjectData: Standard workflow, returns SubjectData with results - nibabel.Nifti1Image: Single image, returns ParcelData - list[nibabel.Nifti1Image]: Batch processing, returns list of results	required

Returns:

Type	Description
SubjectData or ParcelData or list[ParcelData]	Results matching input type: - SubjectData input: New SubjectData instance with results in .results dict - nibabel input: Single ParcelData - list input: List of ParcelData objects (one per input image)

Raises:

Type	Description
ValueError	If input validation fails or source data not found.
TypeError	If input type is not supported.

Notes

This method overrides BaseAnalysis.run() to support flexible input types. The base class run() is designed for SubjectData only.

Examples:

>>> # SubjectData input
>>> mask_data = SubjectData(mask_img, space='MNI152NLin6Asym', resolution=2)
>>> analysis = ParcelAggregation(aggregation='percent')
>>> result = analysis.run(mask_data)
>>> isinstance(data, SubjectData)
True

>>> # Nibabel image input
>>> import nibabel as nib
>>> img = nib.load('mask.nii.gz')
>>> result = analysis.run(img)
>>> isinstance(result, ParcelData)
True

>>> # List of images
>>> images = [nib.load(f'mask_{i}.nii.gz') for i in range(5)]
>>> results = analysis.run(images)
>>> len(results) == 5
True

Source code in src/lacuna/analysis/parcel_aggregation.py

def run(
    self, data: "SubjectData | nib.Nifti1Image | list[nib.Nifti1Image]"
) -> "SubjectData | ParcelData | list[ParcelData]":
    """
    Execute atlas aggregation analysis on various input types.

    Supports flexible input types with matching return types:
    - SubjectData -> SubjectData (with results attached)
    - nibabel.Nifti1Image -> ParcelData
    - list[nibabel.Nifti1Image] -> list[ParcelData]

    Parameters
    ----------
    data : SubjectData or nibabel.Nifti1Image or list[nibabel.Nifti1Image]
        Input data to aggregate:
        - SubjectData: Standard workflow, returns SubjectData with results
        - nibabel.Nifti1Image: Single image, returns ParcelData
        - list[nibabel.Nifti1Image]: Batch processing, returns list of results

    Returns
    -------
    SubjectData or ParcelData or list[ParcelData]
        Results matching input type:
        - SubjectData input: New SubjectData instance with results in .results dict
        - nibabel input: Single ParcelData
        - list input: List of ParcelData objects (one per input image)

    Raises
    ------
    ValueError
        If input validation fails or source data not found.
    TypeError
        If input type is not supported.

    Notes
    -----
    This method overrides BaseAnalysis.run() to support flexible input types.
    The base class run() is designed for SubjectData only.

    Examples
    --------
    >>> # SubjectData input
    >>> mask_data = SubjectData(mask_img, space='MNI152NLin6Asym', resolution=2)
    >>> analysis = ParcelAggregation(aggregation='percent')
    >>> result = analysis.run(mask_data)
    >>> isinstance(data, SubjectData)
    True

    >>> # Nibabel image input
    >>> import nibabel as nib
    >>> img = nib.load('mask.nii.gz')
    >>> result = analysis.run(img)
    >>> isinstance(result, ParcelData)
    True

    >>> # List of images
    >>> images = [nib.load(f'mask_{i}.nii.gz') for i in range(5)]
    >>> results = analysis.run(images)
    >>> len(results) == 5
    True
    """
    from lacuna.core.data_types import VoxelMap

    # Detect input type and delegate to appropriate handler
    if isinstance(data, SubjectData):
        # Standard SubjectData workflow - use base class run()
        return super().run(data)

    elif isinstance(data, VoxelMap):
        # VoxelMap - run directly without SubjectData wrapper
        return self._run_voxelmap(data)

    elif isinstance(data, nib.Nifti1Image):
        # Single nibabel image - return ParcelData
        return self._run_single_image(data)

    elif isinstance(data, list):
        # List of images or VoxelMaps - return list of results
        if not data:
            raise ValueError("Empty list provided - at least one image required")

        # Check if all are VoxelMaps or all are Images
        if all(isinstance(item, VoxelMap) for item in data):
            # Process VoxelMaps directly
            return [self._run_voxelmap(vm) for vm in data]

        elif all(isinstance(img, nib.Nifti1Image) for img in data):
            return self._run_batch_images(data)

        else:
            raise TypeError(
                "When providing a list, all items must be of the same type: "
                "either all VoxelMap or all nibabel.Nifti1Image objects"
            )

    else:
        raise TypeError(
            f"Unsupported input type: {type(data).__name__}\n"
            "Supported types: SubjectData, VoxelMap, nibabel.Nifti1Image, "
            "list[VoxelMap], list[nibabel.Nifti1Image]"
        )

RegionalDamage

Bases: ParcelAggregation

Compute lesion overlap with atlas regions.

This is a convenience wrapper around ParcelAggregation that: - Sets source="maskimg" (analyze the lesion mask) - Sets aggregation="percent" (compute overlap percentages)

This provides a simpler interface for the common use case of computing how much of each brain region is damaged by a lesion.

Attributes:

Name	Type	Description
batch_strategy	str	Batch processing strategy. Set to "sequential" to avoid race conditions with threading backends when accessing shared atlas resources.

Parameters:

Name	Type	Description	Default
parcel_names	list of str or None	Names of atlases from the registry to process (e.g., "schaefer2018parcels100networks7"). If None, all registered atlases are processed. Use list_parcellations() to see available atlases.	None

Raises:

Type	Description
ValueError	If parcel_names contains non-existent atlas names.

Notes

• Results show percentage of each region overlapping with mask
• For more control (e.g., computing volume instead of percent), use ParcelAggregation directly

Examples:

>>> # Use all registered atlases
>>> from lacuna import SubjectData
>>> from lacuna.analysis import RegionalDamage
>>>
>>> mask = SubjectData.from_nifti("mask.nii.gz")
>>> analysis = RegionalDamage()  # Uses all registered atlases
>>> result = analysis.run(mask)
>>>
>>> # Results are in RegionalDamage namespace
>>> overlap_pcts = result.results["RegionalDamage"]
>>> for region, pct in overlap_pcts.items():
...     if pct > 10:  # Show regions with >10% damage
...         print(f"{region}: {pct:.1f}%")
>>>
>>> # Process only specific atlases
>>> analysis = RegionalDamage(
...     parcel_names=["schaefer2018parcels100networks7"]
... )
>>> result = analysis.run(mask)

See Also

ParcelAggregation : More flexible aggregation with custom sources/methods

Source code in src/lacuna/analysis/regional_damage.py

class RegionalDamage(ParcelAggregation):
    """
    Compute lesion overlap with atlas regions.

    This is a convenience wrapper around ParcelAggregation that:
    - Sets source="maskimg" (analyze the lesion mask)
    - Sets aggregation="percent" (compute overlap percentages)

    This provides a simpler interface for the common use case of computing
    how much of each brain region is damaged by a lesion.

    Attributes
    ----------
    batch_strategy : str
        Batch processing strategy. Set to "sequential" to avoid race conditions
        with threading backends when accessing shared atlas resources.

    Parameters
    ----------
    parcel_names : list of str or None, default=None
        Names of atlases from the registry to process (e.g., "schaefer2018parcels100networks7").
        If None, all registered atlases are processed.
        Use list_parcellations() to see available atlases.

    Raises
    ------
    ValueError
        If parcel_names contains non-existent atlas names.

    Notes
    -----
    - Results show percentage of each region overlapping with mask
    - For more control (e.g., computing volume instead of percent),
      use ParcelAggregation directly

    Examples
    --------
    >>> # Use all registered atlases
    >>> from lacuna import SubjectData
    >>> from lacuna.analysis import RegionalDamage
    >>>
    >>> mask = SubjectData.from_nifti("mask.nii.gz")
    >>> analysis = RegionalDamage()  # Uses all registered atlases
    >>> result = analysis.run(mask)
    >>>
    >>> # Results are in RegionalDamage namespace
    >>> overlap_pcts = result.results["RegionalDamage"]
    >>> for region, pct in overlap_pcts.items():
    ...     if pct > 10:  # Show regions with >10% damage
    ...         print(f"{region}: {pct:.1f}%")
    >>>
    >>> # Process only specific atlases
    >>> analysis = RegionalDamage(
    ...     parcel_names=["schaefer2018parcels100networks7"]
    ... )
    >>> result = analysis.run(mask)

    See Also
    --------
    ParcelAggregation : More flexible aggregation with custom sources/methods
    """

    #: Preferred batch processing strategy (sequential to avoid threading race conditions)
    batch_strategy: str = "sequential"

    def __init__(
        self,
        parcel_names: list[str] | None = None,
        verbose: bool = False,
        keep_intermediate: bool = False,
    ):
        """
        Initialize RegionalDamage analysis.

        This is equivalent to:
        ParcelAggregation(source="maskimg",
                        aggregation="percent",
                        parcel_names=parcel_names,
                        verbose=verbose,
                        keep_intermediate=keep_intermediate)

        Parameters
        ----------
        parcel_names : list[str] | None, optional
            List of specific parcellation names to use. If None, uses all available.
        verbose : bool, default=False
            If True, print progress messages. If False, run silently.
        keep_intermediate : bool, default=False
            If True, include intermediate results (e.g., warped mask images)
            in the output. Useful for debugging and quality control.
        """
        super().__init__(
            source="maskimg",
            aggregation="percent",
            parcel_names=parcel_names,
            verbose=verbose,
            keep_intermediate=keep_intermediate,
        )

    def _validate_inputs(self, mask_data) -> None:
        """
        Validate inputs for regional damage analysis.

        Extends parent validation to ensure mask is binary.

        Parameters
        ----------
        mask_data : SubjectData
            Mask data to validate

        Raises
        ------
        ValueError
            If mask is not binary (contains values other than 0 and 1)
        """
        # Run parent validation first
        super()._validate_inputs(mask_data)

        # Check that mask is binary
        import numpy as np

        mask_data_arr = mask_data.mask_img.get_fdata()
        unique_vals = np.unique(mask_data_arr)

        # Binary mask should only have 0 and 1 (or just 0, or just 1)
        if not np.all(np.isin(unique_vals, [0, 1])):
            raise ValueError(
                f"RegionalDamage requires binary mask (0 and 1 only).\n"
                f"Found values: {unique_vals}\n"
                f"Use thresholding or binarization to convert continuous maps."
            )

    def _run_analysis(self, mask_data: SubjectData) -> dict[str, DataContainer]:
        """
        Compute percentage and binary regional damage for all atlases.

        Extends parent to add binary damage results (1 if region has any
        overlap with the mask, 0 otherwise) alongside percentage results.

        Parameters
        ----------
        mask_data : SubjectData
            Validated lesion data.

        Returns
        -------
        dict[str, DataContainer]
            Dictionary with both percentage and binary ParcelData objects.
        """
        from lacuna.core.data_types import ParcelData

        # Get percentage results from parent (keyed as InputMask)
        parent_results = super()._run_analysis(mask_data)

        # Re-key percentage results and derive binary results
        results = {}
        for key, result in parent_results.items():
            if not isinstance(result, ParcelData):
                # Keep non-ParcelData results (e.g., intermediates) as-is
                results[key] = result
                continue

            # Re-key percentage result: atlas-{atlas}_source-RegionalDamage_desc-damagepct
            pct_key = build_result_key(
                atlas=result.name,
                source="RegionalDamage",
                desc="damagepct",
            )
            results[pct_key] = result

            # Derive binary result: >0 → 1, else → 0
            binary_data = {
                region: 1.0 if value > 0 else 0.0 for region, value in result.data.items()
            }

            binary_parcel = ParcelData(
                name=result.name,
                data=binary_data,
                parcel_names=result.parcel_names,
                aggregation_method="binary",
                metadata={
                    **result.metadata,
                    "derived_from": "percent",
                    "description": "Binary damage indicator (1 = any overlap, 0 = none)",
                },
            )

            binary_key = build_result_key(
                atlas=result.name,
                source="RegionalDamage",
                desc="damagebin",
            )
            results[binary_key] = binary_parcel

        return results

    def _get_parameters(self) -> dict:
        """Get analysis parameters for provenance and display.

        Returns
        -------
        dict
            Dictionary of parameter names and values.
        """
        params = super()._get_parameters()
        # RegionalDamage is a specific configuration of ParcelAggregation
        # Override the source and aggregation to reflect the simplified API
        params.update(
            {
                "analysis_type": "RegionalDamage",
                "threshold": params.get("threshold"),
            }
        )
        return params

__init__(parcel_names=None, verbose=False, keep_intermediate=False)

Initialize RegionalDamage analysis.

This is equivalent to: ParcelAggregation(source="maskimg", aggregation="percent", parcel_names=parcel_names, verbose=verbose, keep_intermediate=keep_intermediate)

Parameters:

Name	Type	Description	Default
parcel_names	list[str] | None	List of specific parcellation names to use. If None, uses all available.	None
verbose	bool	If True, print progress messages. If False, run silently.	False
keep_intermediate	bool	If True, include intermediate results (e.g., warped mask images) in the output. Useful for debugging and quality control.	False

Source code in src/lacuna/analysis/regional_damage.py

def __init__(
    self,
    parcel_names: list[str] | None = None,
    verbose: bool = False,
    keep_intermediate: bool = False,
):
    """
    Initialize RegionalDamage analysis.

    This is equivalent to:
    ParcelAggregation(source="maskimg",
                    aggregation="percent",
                    parcel_names=parcel_names,
                    verbose=verbose,
                    keep_intermediate=keep_intermediate)

    Parameters
    ----------
    parcel_names : list[str] | None, optional
        List of specific parcellation names to use. If None, uses all available.
    verbose : bool, default=False
        If True, print progress messages. If False, run silently.
    keep_intermediate : bool, default=False
        If True, include intermediate results (e.g., warped mask images)
        in the output. Useful for debugging and quality control.
    """
    super().__init__(
        source="maskimg",
        aggregation="percent",
        parcel_names=parcel_names,
        verbose=verbose,
        keep_intermediate=keep_intermediate,
    )

StructuralNetworkMapping

Bases: BaseAnalysis

Structural lesion network mapping using tractography-based disconnection.

This analysis quantifies white matter disconnection caused by a lesion by: 1. Filtering a whole-brain tractogram to streamlines passing through the lesion 2. Computing track density images (TDI) for both lesion and whole-brain tracts 3. Calculating disconnection probability as the ratio of lesion TDI to whole-brain TDI 4. Optionally computing parcellated connectivity matrices (if atlas provided)

Outputs: - Voxel-wise disconnection map (always): 3D NIfTI showing disconnection probability per voxel - Connectivity matrices (optional): Parcellated edge-wise disconnection when atlas is provided

Computation Space: All computations are performed in the tractogram's native space (typically MNI152NLin2009cAsym @ 1mm for high-resolution structural connectivity). Lesions are transformed to this high-resolution template space for accurate white matter fiber filtering.

The analysis requires MRtrix3 to be installed and available in the system PATH.

Attributes:

Name	Type	Description
batch_strategy	str	Batch processing strategy. Set to "parallel" as each subject is independent and tractogram filtering is compute-intensive.

Parameters:

Name	Type	Description	Default
connectome_name	str	Name of registered structural connectome (e.g., "dTOR985"). Use list_structural_connectomes() to see available connectomes. The connectome must be pre-registered via register_structural_connectome().	required
parcellation_name	str	Name of registered atlas for parcellated connectivity matrices. Use list_parcellations() to see available atlases.	None
compute_disconnectivity_matrix	bool	If True and parcellation_name provided, compute disconnectivity matrices including per-ROI disconnection percentages.	False
output_resolution	(1, 2)	Output resolution in mm (must match connectome resolution).	1
cache_tdi	bool	If True, cache computed TDI maps for reuse.	True
n_jobs	int	Number of threads for MRtrix3 processing.	1
keep_intermediate	bool	If True, include intermediate results in output (lesion tractogram, lesion TDI, warped atlas if transformed). Useful for debugging and quality control. Intermediate files are exported to the output directory alongside other results; temporary working directories are always cleaned up.	False
check_dependencies	bool	If True, checks for MRtrix3 availability.	True
verbose	bool	If True, print progress messages. If False, run silently.	True
show_mrtrix_output	bool	If True, display MRtrix3 command outputs. If False, suppress verbose MRtrix3 messages for cleaner output.	False

Raises:

Type	Description
MRtrixError	If MRtrix3 is not installed or not available in PATH
FileNotFoundError	If tractogram_path or whole_brain_tdi files don't exist

Examples:

Register and use structural connectome:

>>> from lacuna import SubjectData
>>> from lacuna.analysis import StructuralNetworkMapping
>>> from lacuna.assets.connectomes import (
...     list_structural_connectomes,
...     register_structural_connectome,
... )
>>>
>>> # Register connectome (do this once)
>>> register_structural_connectome(
...     name="dTOR985",
...     space="MNI152NLin2009cAsym",
...     tractogram_path="/data/dtor/dtor985_tractogram.tck",
...     description="HCP dTOR tractogram - TDI computed on-demand"
... )
>>>
>>> # List available connectomes
>>> list_structural_connectomes()
>>>
>>> # Load lesion data
>>> lesion = SubjectData.from_nifti("lesion.nii.gz")
>>>
>>> # Interactive analysis
>>> analysis = StructuralNetworkMapping(
...     connectome_name="dTOR985",
...     n_jobs=8,
... )
>>> result = analysis.run(lesion)
>>> disconn_map = result.results["StructuralNetworkMapping"]["disconnection_pct"]
>>> disconn_map.orthoview()

Batch processing:

>>> from lacuna import batch_process
>>>
>>> analysis = StructuralNetworkMapping(
...     connectome_name="dTOR985",
...     n_jobs=8,
... )
>>> results = batch_process(lesions, analysis, n_jobs=2)
>>>
>>> # Save results
>>> for result in results:
...     disconn_map = result.results["StructuralNetworkMapping"]["disconnection_pct"]
...     nib.save(disconn_map, f"output/{subject_id}_disconn.nii.gz")

Notes

• Requires MRtrix3: https://www.mrtrix.org/download/
• Processing time scales with lesion size and tractogram density
• For large tractograms, processing can take several minutes per subject

See Also

FunctionalNetworkMapping : Functional connectivity-based lesion network mapping RegionalDamage : Atlas-based regional overlap quantification

Source code in src/lacuna/analysis/structural_network_mapping.py

class StructuralNetworkMapping(BaseAnalysis):
    """
    Structural lesion network mapping using tractography-based disconnection.

    This analysis quantifies white matter disconnection caused by a lesion by:
    1. Filtering a whole-brain tractogram to streamlines passing through the lesion
    2. Computing track density images (TDI) for both lesion and whole-brain tracts
    3. Calculating disconnection probability as the ratio of lesion TDI to whole-brain TDI
    4. Optionally computing parcellated connectivity matrices (if atlas provided)

    **Outputs:**
    - **Voxel-wise disconnection map** (always): 3D NIfTI showing disconnection probability per voxel
    - **Connectivity matrices** (optional): Parcellated edge-wise disconnection when atlas is provided

    **Computation Space:**
    All computations are performed in the tractogram's native space (typically
    MNI152NLin2009cAsym @ 1mm for high-resolution structural connectivity).
    Lesions are transformed to this high-resolution template space for accurate
    white matter fiber filtering.

    The analysis requires MRtrix3 to be installed and available in the system PATH.

    Attributes
    ----------
    batch_strategy : str
        Batch processing strategy. Set to "parallel" as each subject is independent
        and tractogram filtering is compute-intensive.

    Parameters
    ----------
    connectome_name : str
        Name of registered structural connectome (e.g., "dTOR985").
        Use list_structural_connectomes() to see available connectomes.
        The connectome must be pre-registered via register_structural_connectome().
    parcellation_name : str, optional
        Name of registered atlas for parcellated connectivity matrices.
        Use list_parcellations() to see available atlases.
    compute_disconnectivity_matrix : bool, default=False
        If True and parcellation_name provided, compute disconnectivity matrices
        including per-ROI disconnection percentages.
    output_resolution : {1, 2}, default=2
        Output resolution in mm (must match connectome resolution).
    cache_tdi : bool, default=True
        If True, cache computed TDI maps for reuse.
    n_jobs : int, default=1
        Number of threads for MRtrix3 processing.
    keep_intermediate : bool, default=False
        If True, include intermediate results in output (lesion tractogram, lesion TDI,
        warped atlas if transformed). Useful for debugging and quality control.
        Intermediate files are exported to the output directory alongside other results;
        temporary working directories are always cleaned up.
    check_dependencies : bool, default=True
        If True, checks for MRtrix3 availability.
    verbose : bool, default=True
        If True, print progress messages. If False, run silently.
    show_mrtrix_output : bool, default=False
        If True, display MRtrix3 command outputs. If False, suppress verbose
        MRtrix3 messages for cleaner output.

    Raises
    ------
    MRtrixError
        If MRtrix3 is not installed or not available in PATH
    FileNotFoundError
        If tractogram_path or whole_brain_tdi files don't exist

    Examples
    --------
    **Register and use structural connectome:**

    >>> from lacuna import SubjectData
    >>> from lacuna.analysis import StructuralNetworkMapping
    >>> from lacuna.assets.connectomes import (
    ...     list_structural_connectomes,
    ...     register_structural_connectome,
    ... )
    >>>
    >>> # Register connectome (do this once)
    >>> register_structural_connectome(
    ...     name="dTOR985",
    ...     space="MNI152NLin2009cAsym",
    ...     tractogram_path="/data/dtor/dtor985_tractogram.tck",
    ...     description="HCP dTOR tractogram - TDI computed on-demand"
    ... )
    >>>
    >>> # List available connectomes
    >>> list_structural_connectomes()
    >>>
    >>> # Load lesion data
    >>> lesion = SubjectData.from_nifti("lesion.nii.gz")
    >>>
    >>> # Interactive analysis
    >>> analysis = StructuralNetworkMapping(
    ...     connectome_name="dTOR985",
    ...     n_jobs=8,
    ... )
    >>> result = analysis.run(lesion)
    >>> disconn_map = result.results["StructuralNetworkMapping"]["disconnection_pct"]
    >>> disconn_map.orthoview()

    **Batch processing:**

    >>> from lacuna import batch_process
    >>>
    >>> analysis = StructuralNetworkMapping(
    ...     connectome_name="dTOR985",
    ...     n_jobs=8,
    ... )
    >>> results = batch_process(lesions, analysis, n_jobs=2)
    >>>
    >>> # Save results
    >>> for result in results:
    ...     disconn_map = result.results["StructuralNetworkMapping"]["disconnection_pct"]
    ...     nib.save(disconn_map, f"output/{subject_id}_disconn.nii.gz")

    Notes
    -----
    - Requires MRtrix3: https://www.mrtrix.org/download/
    - Processing time scales with lesion size and tractogram density
    - For large tractograms, processing can take several minutes per subject

    See Also
    --------
    FunctionalNetworkMapping : Functional connectivity-based lesion network mapping
    RegionalDamage : Atlas-based regional overlap quantification
    """

    #: Preferred batch processing strategy - sequential because MRtrix3's tckedit
    #: uses internal parallelization (-nthreads) and running multiple instances
    #: in parallel causes resource contention and memory-mapping conflicts
    batch_strategy: str = "sequential"

    def __init__(
        self,
        connectome_name: str,
        parcellation_name: str | list[str] | None = None,
        compute_disconnectivity_matrix: bool = False,
        compute_roi_disconnection: bool = False,
        include_streamline_counts: bool = False,
        output_resolution: Literal[1, 2] = 2,
        cache_tdi: bool = True,
        n_jobs: int = 1,
        keep_intermediate: bool = False,
        check_dependencies: bool = True,
        verbose: bool = False,
        show_mrtrix_output: bool = False,
        return_in_input_space: bool = True,
    ):
        """Initialize StructuralNetworkMapping analysis.

        Parameters
        ----------
        connectome_name : str
            Name of registered structural connectome (e.g., "HCP842_dTOR").
            Use list_structural_connectomes() to see available options.
        parcellation_name : str, optional
            Name of registered atlas for parcellated connectivity matrices.
        compute_disconnectivity_matrix : bool, default=False
            If True and parcellation_name provided, compute disconnectivity matrices
            (full, mask, and disconnectivity percentage connectivity matrices).
        compute_roi_disconnection : bool, default=False
            If True and parcellation_name provided, compute per-ROI disconnection
            percentages and intact (post-disconnection) connectivity matrix.
        include_streamline_counts : bool, default=False
            If True, include raw streamline count outputs (disconnection_tdi VoxelMap).
            By default, only percentage-based outputs are produced.
        output_resolution : {1, 2}, default=2
            Output resolution in mm (must match connectome resolution).
        cache_tdi : bool, default=True
            If True, cache computed TDI maps.
        n_jobs : int, default=1
            Number of threads for MRtrix3.
        keep_intermediate : bool, default=False
            If True, include intermediate results in output (lesion tractogram,
            lesion TDI, warped atlas). Useful for debugging and QC.
        check_dependencies : bool, default=True
            If True, checks for MRtrix3 availability.
        verbose : bool, default=True
            If True, print progress messages. If False, run silently.
        show_mrtrix_output : bool, default=False
            If True, display MRtrix3 command outputs. If False, suppress verbose
            MRtrix3 messages for cleaner output.
        return_in_input_space : bool, default=True
            If True, transform VoxelMap outputs back to the original input mask space.
            If False, outputs remain in the connectome space.
            Requires input SubjectData to have valid space/resolution metadata.

        Raises
        ------
        MRtrixError
            If MRtrix3 is not available and check_dependencies=True.
        KeyError
            If connectome_name not found in registry.
        ValueError
            If output_resolution is not 1 or 2.
        """
        super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

        # Validate output_resolution
        if output_resolution not in (1, 2):
            raise ValueError(f"output_resolution must be 1 or 2, got: {output_resolution}")

        # Load connectome from registry
        try:
            connectome = load_structural_connectome(connectome_name)
        except KeyError as e:
            available = [c.name for c in list_structural_connectomes()]
            raise KeyError(
                f"Connectome '{connectome_name}' not found in registry. "
                f"Available connectomes: {', '.join(available)}. "
                f"Use register_structural_connectome() to add new connectomes."
            ) from e

        # Store connectome information
        self.connectome_name = connectome_name
        self.tractogram_path = connectome.tractogram_path
        self.tractogram_space = connectome.metadata.space
        self.template = connectome.template_path  # May be None

        # Store analysis parameters — normalize parcellation_name(s) to list
        if parcellation_name is None:
            self.parcellation_names: list[str] = []
        elif isinstance(parcellation_name, str):
            self.parcellation_names = [parcellation_name]
        else:
            self.parcellation_names = list(parcellation_name)
        # Keep singular attribute for backward compatibility
        self.parcellation_name = self.parcellation_names[0] if self.parcellation_names else None

        self.compute_disconnectivity_matrix = compute_disconnectivity_matrix
        self.compute_roi_disconnection = compute_roi_disconnection
        self.include_streamline_counts = include_streamline_counts

        # Validate: compute flags require parcellation_name
        if (
            compute_disconnectivity_matrix or compute_roi_disconnection
        ) and not self.parcellation_names:
            flags = []
            if compute_disconnectivity_matrix:
                flags.append("compute_disconnectivity_matrix")
            if compute_roi_disconnection:
                flags.append("compute_roi_disconnection")
            raise ValueError(
                f"{' and '.join(flags)} requires parcellation_name to be specified "
                f"(via --parcel-atlases or --custom-parcellation). "
                f"Use list_parcellations() to see available options."
            )

        # Validate parcellation names early (full loading deferred to _validate_inputs)
        if self.parcellation_names:
            available_names = [a.name for a in list_parcellations()]
            for name in self.parcellation_names:
                if name not in available_names:
                    raise ValueError(
                        f"Atlas '{name}' not found in registry. "
                        f"Available parcellations: {', '.join(available_names[:10])}... "
                        f"Use list_parcellations() to see all options."
                    )

        self.output_resolution = output_resolution
        self.cache_tdi = cache_tdi
        self.n_jobs = n_jobs
        if n_jobs != -1 and n_jobs < 1:
            raise ValueError(f"n_jobs must be -1 (all CPUs) or >= 1, got {n_jobs}")
        self.keep_intermediate = keep_intermediate

        self.show_mrtrix_output = show_mrtrix_output
        self.return_in_input_space = return_in_input_space

        # Target space matches tractogram space, resolution from output_resolution
        # (used by BaseAnalysis._ensure_target_space)
        self.TARGET_SPACE = self.tractogram_space
        self.TARGET_RESOLUTION = self.output_resolution

        # Initialize logger
        self.logger = ConsoleLogger(verbose=verbose, width=70)

        # Internal state
        self.whole_brain_tdi = None  # Will be set during validation
        # Per-atlas state: list of dicts with keys: name, image, labels, resolved_path, was_transformed, etc.
        self._atlases: list[dict] = []
        self._full_connectivity_matrices: dict[str, np.ndarray] = {}  # Cached per atlas name
        # Keep legacy attributes for backward compat (set to first atlas after validation)
        self._atlas_image = None
        self._atlas_labels = None
        self._parcellation_resolved = None
        self._full_connectivity_matrix = None
        self._cached_tdi_path = None

        # Check MRtrix3 availability if requested
        if check_dependencies:
            try:
                check_mrtrix_available()
            except MRtrixError as e:
                raise MRtrixError(
                    f"MRtrix3 is required for StructuralNetworkMapping but is not available.\n{e}"
                ) from e

    def _get_tdi_cache_path(self) -> Path:
        """Get deterministic cache path for whole-brain TDI.

        Returns
        -------
        Path
            Path to cached TDI file, unique to tractogram and resolution
        """
        # Create deterministic hash from tractogram path and resolution
        tractogram_str = str(self.tractogram_path.resolve())
        hash_input = f"{tractogram_str}_{self.output_resolution}mm"
        file_hash = hashlib.md5(hash_input.encode()).hexdigest()[:12]

        # Use unified cache directory
        cache_dir = get_tdi_cache_dir()

        cache_filename = f"tdi_{file_hash}_{self.output_resolution}mm.nii.gz"
        return cache_dir / cache_filename

    def _compute_tdi_to_path(self, output_path: Path) -> None:
        """Compute whole-brain TDI and save to specified path.

        Parameters
        ----------
        output_path : Path
            Where to save the computed TDI
        """
        compute_tdi_map(
            tractogram_path=self.tractogram_path,
            template=self.template,
            output_path=output_path,
            n_jobs=self.n_jobs,
            verbose=self.show_mrtrix_output,
        )

    def _ensure_tdi_cached(self, cache_path: Path) -> None:
        """Ensure TDI is computed and cached, with file locking for parallel safety.

        Uses file locking to prevent race conditions when multiple parallel
        workers try to compute the same TDI simultaneously. The first worker
        to acquire the lock computes the TDI; others wait and use the cached file.

        Parameters
        ----------
        cache_path : Path
            Cache file path from _get_tdi_cache_path()
        """
        import fcntl

        # Check if already cached (fast path, no lock needed)
        if cache_path.exists():
            self.logger.info(f"Using cached TDI: {cache_path}")
            return

        # Use a lock file to coordinate parallel workers
        lock_path = cache_path.with_suffix(".lock")
        lock_path.parent.mkdir(parents=True, exist_ok=True)

        with open(lock_path, "w") as lock_file:
            self.logger.debug(f"Acquiring TDI cache lock: {lock_path}")
            fcntl.flock(lock_file.fileno(), fcntl.LOCK_EX)
            try:
                # Re-check after acquiring lock (another worker may have computed it)
                if cache_path.exists():
                    self.logger.info(
                        f"Using cached TDI (computed by another process): {cache_path}"
                    )
                    return

                # We have the lock and TDI doesn't exist - compute it
                self.logger.info(
                    f"Computing whole-brain TDI at {self.output_resolution}mm resolution..."
                )
                self._compute_tdi_to_path(cache_path)
                self.logger.info(f"Cached TDI to: {cache_path}")
            finally:
                fcntl.flock(lock_file.fileno(), fcntl.LOCK_UN)

        # Clean up lock file (best effort, ignore errors)
        try:
            lock_path.unlink(missing_ok=True)
        except OSError:
            pass

    def _compute_and_cache_tdi(self, cache_path: Path) -> None:
        """Compute whole-brain TDI and save to cache.

        Parameters
        ----------
        cache_path : Path
            Cache file path from _get_tdi_cache_path()
        """
        self._compute_tdi_to_path(cache_path)
        self.logger.info(f"Cached TDI to: {cache_path}")

    def run(self, mask_data: SubjectData) -> SubjectData:
        """Run structural network mapping analysis.

        Automatically transforms lesion to tractogram space if needed.

        Parameters
        ----------
        mask_data : SubjectData
            Lesion data to analyze (can be in any MNI152 space)

        Returns
        -------
        SubjectData
            Analysis results

        Raises
        ------
        ValueError
            If lesion is in native space (cannot transform)
        """
        # Check for native space (cannot transform)
        space = mask_data.get_coordinate_space()

        if space.lower() == "native":
            raise ValueError(
                "Native space lesions are not supported for structural network mapping. "
                f"Lesions must be in a standard space. Tractogram is in {self.tractogram_space}."
            )

        # Transform to tractogram space (handled by base class)
        # The base class will handle space equivalence and transformations
        return super().run(mask_data)

    def _validate_inputs(self, mask_data: SubjectData) -> None:
        """
        Validate that lesion data meets requirements for structural network mapping.

        This method validates that the mask data is ready for SNM analysis and
        performs essential setup (template loading, TDI computation, atlas loading).
        By the time this is called, BaseAnalysis.run() has already transformed
        the mask to TARGET_SPACE (the tractogram space) via _ensure_target_space().

        Parameters
        ----------
        mask_data : SubjectData
            Lesion data to validate (already transformed to tractogram space).

        Raises
        ------
        ValidationError
            If mask space doesn't match the expected tractogram space.
        FileNotFoundError
            If required input files don't exist (tractogram, template).

        Notes
        -----
        Binary mask validation is handled by SubjectData.__init__, so we don't
        need to duplicate that check here.

        Space transformation is handled by BaseAnalysis._ensure_target_space(),
        so by the time we get here, mask_data.space should equal self.TARGET_SPACE.
        """
        # Validate that required files exist
        if not self.tractogram_path.exists():
            raise FileNotFoundError(f"Tractogram file not found: {self.tractogram_path}")

        # Validate coordinate space matches tractogram space
        # (should already be transformed by _ensure_target_space)
        from lacuna.core.spaces import spaces_are_equivalent

        if not spaces_are_equivalent(mask_data.space, self.TARGET_SPACE):
            raise ValueError(
                f"Mask space '{mask_data.space}' does not match tractogram space "
                f"'{self.TARGET_SPACE}'. This is unexpected - space transformation "
                f"should have been handled by BaseAnalysis._ensure_target_space()."
            )

        # Load template from asset management if not provided (MUST BE DONE BEFORE TDI COMPUTATION)
        if self.template is None:
            # Use output_resolution for template (not lesion resolution)
            # This ensures TDI and template match
            template_name = f"{self.TARGET_SPACE}_res-{self.output_resolution}"

            try:
                template_path = load_template(template_name)
                self.template = template_path
            except (KeyError, FileNotFoundError) as e:
                raise FileNotFoundError(
                    f"Could not load template for {self.TARGET_SPACE} at {self.output_resolution}mm resolution. "
                    f"Template '{template_name}' not found in registry."
                ) from e
        else:
            self.template = Path(self.template)

        # Compute or load TDI with caching (template must be set first!)
        if self.cache_tdi:
            tdi_cache_path = self._get_tdi_cache_path()
            # Use file locking to prevent race conditions in parallel processing
            # Multiple workers may try to compute the same TDI simultaneously
            self._ensure_tdi_cached(tdi_cache_path)
            self.whole_brain_tdi = tdi_cache_path
        else:
            # Compute TDI without caching (temporary file)
            temp_tdi = get_temp_dir(prefix="tdi_") / "whole_brain_tdi.nii.gz"
            self.logger.info(
                f"Computing whole-brain TDI at {self.output_resolution}mm resolution..."
            )
            self._compute_tdi_to_path(temp_tdi)
            self.whole_brain_tdi = temp_tdi

        # Verify template exists
        if not self.template.exists():
            raise FileNotFoundError(f"Template not found: {self.template}")

        # Load atlases from registry
        self._atlases = []
        for parc_name in self.parcellation_names:
            try:
                atlas = load_parcellation(parc_name)
                atlas_info = {
                    "name": parc_name,
                    "image": atlas.image,
                    "labels": atlas.labels,
                    "resolved_path": None,
                    "was_transformed": False,
                    "original_space": None,
                    "original_resolution": None,
                }

                # Check if atlas space matches tractogram space
                atlas_space = atlas.metadata.space
                atlas_resolution = atlas.metadata.resolution

                if atlas_space != self.tractogram_space:
                    self.logger.info(
                        f"Atlas space ({atlas_space}) differs from tractogram space "
                        f"({self.tractogram_space}). Transforming atlas '{parc_name}'..."
                    )

                    # Transform atlas to tractogram space
                    from lacuna.core.spaces import CoordinateSpace
                    from lacuna.spatial.transform import transform_image
                    from lacuna.utils.cache import get_cache_dir

                    # Define target space matching tractogram
                    template_img = (
                        nib.load(self.template)
                        if isinstance(self.template, (str, Path))
                        else self.template
                    )
                    target_space = CoordinateSpace(
                        identifier=self.tractogram_space,
                        resolution=self.output_resolution,
                        reference_affine=template_img.affine,
                    )

                    transformed_atlas_img = transform_image(
                        img=atlas_info["image"],
                        source_space=atlas_space,
                        target_space=target_space,
                        source_resolution=atlas_resolution,
                        interpolation="nearest",
                        image_name=f"atlas '{parc_name}'",
                        verbose=self.verbose,
                    )

                    # Save transformed atlas to cache
                    atlas_cache_dir = get_cache_dir() / "atlases"
                    atlas_cache_dir.mkdir(exist_ok=True, parents=True)

                    atlas_hash = hashlib.md5(
                        f"{parc_name}_{self.tractogram_space}_{self.output_resolution}".encode()
                    ).hexdigest()[:12]
                    transformed_atlas_path = atlas_cache_dir / f"atlas_{atlas_hash}.nii.gz"

                    nib.save(transformed_atlas_img, transformed_atlas_path)
                    atlas_info["resolved_path"] = transformed_atlas_path
                    atlas_info["image"] = transformed_atlas_img
                    atlas_info["was_transformed"] = True
                    atlas_info["original_space"] = atlas_space
                    atlas_info["original_resolution"] = atlas_resolution

                    self.logger.info(f"Atlas transformed and cached to: {transformed_atlas_path}")
                else:
                    # No transformation needed - use original atlas file
                    from lacuna.assets.parcellations.loader import BUNDLED_PARCELLATIONS_DIR

                    atlas_filename_path = Path(atlas.metadata.parcellation_filename)
                    if atlas_filename_path.is_absolute():
                        atlas_info["resolved_path"] = atlas_filename_path
                    else:
                        atlas_info["resolved_path"] = (
                            BUNDLED_PARCELLATIONS_DIR / atlas.metadata.parcellation_filename
                        )

                    if not atlas_info["resolved_path"].exists():
                        raise FileNotFoundError(
                            f"Atlas file not found: {atlas_info['resolved_path']}"
                        )

                self._atlases.append(atlas_info)

            except KeyError as e:
                available = [a.name for a in list_parcellations()]
                raise ValueError(
                    f"Atlas '{parc_name}' not found in registry. "
                    f"Available parcellations: {', '.join(available[:5])}... "
                    f"Use list_parcellations() to see all options."
                ) from e

        # Set legacy single-atlas attributes from first atlas (backward compat)
        if self._atlases:
            first = self._atlases[0]
            self._atlas_image = first["image"]
            self._atlas_labels = first["labels"]
            self._parcellation_resolved = first["resolved_path"]
            self._atlas_was_transformed = first["was_transformed"]
            if first["was_transformed"]:
                self._original_atlas_space = first["original_space"]
                self._original_atlas_resolution = first["original_resolution"]

    def _run_analysis(self, mask_data: SubjectData) -> dict[str, "AnalysisResult"]:
        """
        Execute structural network mapping analysis.

        Parameters
        ----------
        mask_data : SubjectData
            Input lesion data

        Returns
        -------
        dict[str, AnalysisResult]
            Dictionary mapping result names to results:
            - 'disconnection_pct': VoxelMap for disconnection percentage map
            - 'disconnection_tdi': VoxelMap for raw streamline count (lesion TDI)
            - 'summary_statistics': ScalarMetric for summary statistics
            - 'mask_tractogram': Tractogram (if keep_intermediate=True)
            - Connectivity results (if atlas provided): see _compute_connectivity_matrices

        Notes
        -----
        Processing steps:
        1. Filter tractogram to streamlines passing through mask (tckedit)
        2. Compute TDI from filtered tractogram (tckmap)
        3. Compute disconnection as ratio of mask TDI to whole-brain TDI (mrcalc)
        """
        # Get subject ID for informative output
        subject_id = mask_data.metadata.get("subject_id", "unknown")

        # Handle empty masks: produce zero-valued output maps
        if mask_data.is_empty_mask:
            self.logger.warning("Empty mask detected — producing zero-valued disconnection maps")
            return self._build_empty_mask_results(mask_data)

        # Log analysis start
        self.logger.info("Filtering tractogram by mask...")

        # Create temporary directory for intermediate files
        temp_dir_path = get_temp_dir(prefix=f"slnm_{subject_id}_")

        if self.keep_intermediate:
            self.logger.info(f"Intermediate files will be saved to: {temp_dir_path}")

        try:
            # Step 1: Filter tractogram by mask
            mask_tck_path = temp_dir_path / "mask_streamlines.tck"
            filter_tractogram_by_mask(
                tractogram_path=self.tractogram_path,
                mask=mask_data.mask_img,
                output_path=mask_tck_path,
                n_jobs=self.n_jobs,
                force=True,
                verbose=self.show_mrtrix_output,
            )
            self.logger.success("Tractogram filtered")

            # Step 2: Compute TDI from mask-filtered tractogram
            self.logger.info("Computing track-density image (TDI)...")
            # Use anatomical template to define output grid
            mask_tdi_path = temp_dir_path / "mask_tdi.nii.gz"
            compute_tdi_map(
                tractogram_path=mask_tck_path,
                template=self.template,  # Use anatomical template
                output_path=mask_tdi_path,
                n_jobs=self.n_jobs,
                force=True,
                verbose=self.show_mrtrix_output,
            )

            # Step 3: Compute disconnection map
            self.logger.info("Computing disconnection map...")
            disconn_map_path = temp_dir_path / "disconnection_map.nii.gz"
            compute_disconnection_map(
                mask_tdi=mask_tdi_path,
                whole_brain_tdi=self.whole_brain_tdi,
                output_path=disconn_map_path,
                force=True,
                verbose=self.show_mrtrix_output,
            )

            # Load results
            # Use memory-mapped loading for computing statistics efficiently
            disconn_map = nib.load(disconn_map_path, mmap=True)

            # Compute summary statistics (this will load data temporarily but release it)
            disconn_array = np.nan_to_num(disconn_map.get_fdata(), nan=0.0, posinf=0.0, neginf=0.0)
            positive_mask = disconn_array > 0
            mean_disconnection = (
                float(np.mean(disconn_array[positive_mask])) if np.any(positive_mask) else 0.0
            )

            # Free memory immediately after computing statistics
            del disconn_array

            # Count streamlines in mask tractogram (from TDI sum)
            mask_tdi = nib.load(mask_tdi_path, mmap=True)
            mask_streamline_count = int(np.sum(mask_tdi.get_fdata()))

            # Load disconnection map into memory
            # This ensures results are independent of temp directory lifecycle
            disconn_data = nib.load(disconn_map_path).get_fdata()
            # Replace non-finite values (NaN/Inf from division by zero) with 0
            disconn_data = np.nan_to_num(disconn_data, nan=0.0, posinf=0.0, neginf=0.0)
            final_disconn_map = nib.Nifti1Image(
                disconn_data, disconn_map.affine, disconn_map.header
            )

            # Build results dict
            results = {}

            # VoxelMapResult for disconnection map
            disconnection_result = VoxelMap(
                name="disconnection_pct",
                data=final_disconn_map,
                space=self.tractogram_space,
                resolution=float(self.output_resolution),
                metadata={
                    "tractogram": str(self.tractogram_path),
                    "whole_brain_tdi": str(self.whole_brain_tdi),
                    "template": str(self.template),
                    "n_jobs": self.n_jobs,
                    "keep_intermediate": self.keep_intermediate,
                },
            )
            results["disconnection_pct"] = disconnection_result

            # VoxelMap for raw streamline count (lesion TDI) - optional
            if self.include_streamline_counts:
                mask_tdi_data = nib.load(mask_tdi_path).get_fdata()
                final_mask_tdi = nib.Nifti1Image(
                    mask_tdi_data, disconn_map.affine, disconn_map.header
                )
                disconnection_tdi_result = VoxelMap(
                    name="disconnection_tdi",
                    data=final_mask_tdi,
                    space=self.tractogram_space,
                    resolution=float(self.output_resolution),
                    metadata={
                        "description": "Track density image for mask-filtered tractogram (raw streamline counts)",
                        "tractogram": str(self.tractogram_path),
                    },
                )
                results["disconnection_tdi"] = disconnection_tdi_result

            # MiscResult for summary statistics
            summary_result = ScalarMetric(
                name="summarystatistics",
                data={
                    "mean_disconnection": mean_disconnection,
                    "mask_streamline_count": mask_streamline_count,
                },
                metadata={
                    "tractogram": str(self.tractogram_path),
                },
            )
            results["summarystatistics"] = summary_result

            # Add intermediate results if keep_intermediate=True
            # Copy files to a persistent staging dir so they survive temp cleanup
            if self.keep_intermediate:
                import shutil

                staging_dir = get_temp_dir(prefix=f"slnm_intermediates_{subject_id}_")
                staged_tck = staging_dir / mask_tck_path.name
                shutil.copy2(mask_tck_path, staged_tck)

                mask_tractogram_result = Tractogram(
                    name="mask_tractogram",
                    tractogram_path=staged_tck,
                    metadata={
                        "description": "Tractogram filtered by mask",
                    },
                )
                results["mask_tractogram"] = mask_tractogram_result

                # Add warped atlases if transformed
                for atlas_info in self._atlases:
                    if atlas_info["was_transformed"]:
                        warped_atlas_result = VoxelMap(
                            name=f"warped_atlas_{atlas_info['name']}",
                            data=atlas_info["image"],
                            space=self.tractogram_space,
                            resolution=self.output_resolution,
                            metadata={
                                "description": (
                                    f"Atlas '{atlas_info['name']}' transformed from "
                                    f"{atlas_info['original_space']}@{atlas_info['original_resolution']}mm "
                                    f"to {self.tractogram_space}@{self.output_resolution}mm"
                                ),
                                "original_space": atlas_info["original_space"],
                                "original_resolution": atlas_info["original_resolution"],
                                "parcellation_name": atlas_info["name"],
                                "cached_path": str(atlas_info["resolved_path"]),
                            },
                        )
                        results[f"warped_atlas_{atlas_info['name']}"] = warped_atlas_result

            # Optional: Compute parcellated connectivity matrices per atlas
            if self._atlases and (
                self.compute_disconnectivity_matrix or self.compute_roi_disconnection
            ):
                for atlas_info in self._atlases:
                    self.logger.info(
                        f"Computing connectivity matrices for '{atlas_info['name']}'..."
                    )
                    connectivity_results = self._compute_connectivity_matrices(
                        mask_data=mask_data,
                        mask_tck_path=mask_tck_path,
                        temp_dir_path=temp_dir_path,
                        subject_id=subject_id,
                        atlas_info=atlas_info,
                    )
                    results.update(connectivity_results)

                # Save intact tractogram as intermediate (atlas-independent, only one copy)
                if self.keep_intermediate and self.compute_roi_disconnection:
                    intact_tck_path = temp_dir_path / f"{subject_id}_intact.tck"
                    if intact_tck_path.exists():
                        staged_intact = staging_dir / intact_tck_path.name
                        shutil.copy2(intact_tck_path, staged_intact)
                        intact_tractogram_result = Tractogram(
                            name="intact_tractogram",
                            tractogram_path=staged_intact,
                            metadata={
                                "description": "Intact (post-disconnection) tractogram",
                            },
                        )
                        results["intact_tractogram"] = intact_tractogram_result

            # Transform VoxelMap results back to input space if requested
            if self.return_in_input_space:
                results = self._transform_results_to_input_space(results, mask_data)

            self.logger.success(f"Analysis complete ({len(results)} results)")
            return results

        finally:
            # Always clean up temporary working directory
            import shutil

            shutil.rmtree(temp_dir_path, ignore_errors=True)

    def _build_empty_mask_results(self, mask_data: SubjectData) -> dict[str, "AnalysisResult"]:
        """Build zero-valued results for an empty mask.

        Parameters
        ----------
        mask_data : SubjectData
            Empty-mask subject data (already in tractogram space).

        Returns
        -------
        dict[str, AnalysisResult]
            Zero-valued VoxelMaps and summary statistics matching the
            structure returned by the normal analysis path.
        """
        # Load template to get output shape and affine
        template_img = nib.load(self.template)
        template_shape = template_img.shape[:3]
        template_affine = template_img.affine

        zero_vol = np.zeros(template_shape, dtype=np.float32)

        results: dict[str, AnalysisResult] = {}

        results["disconnection_pct"] = VoxelMap(
            name="disconnection_pct",
            data=nib.Nifti1Image(zero_vol.copy(), template_affine),
            space=self.tractogram_space,
            resolution=float(self.output_resolution),
            metadata={
                "tractogram": str(self.tractogram_path),
                "whole_brain_tdi": str(self.whole_brain_tdi),
                "template": str(self.template),
                "empty_mask": True,
            },
        )

        if self.include_streamline_counts:
            results["disconnection_tdi"] = VoxelMap(
                name="disconnection_tdi",
                data=nib.Nifti1Image(zero_vol.copy(), template_affine),
                space=self.tractogram_space,
                resolution=float(self.output_resolution),
                metadata={
                    "description": "Track density image for mask-filtered tractogram (raw streamline counts)",
                    "tractogram": str(self.tractogram_path),
                    "empty_mask": True,
                },
            )

        results["summarystatistics"] = ScalarMetric(
            name="summarystatistics",
            data={
                "mean_disconnection": 0.0,
                "mask_streamline_count": 0,
                "empty_mask": True,
            },
            metadata={"tractogram": str(self.tractogram_path)},
        )

        # Compute connectivity matrices with zeros if atlas provided
        if self._atlases and (
            self.compute_disconnectivity_matrix or self.compute_roi_disconnection
        ):
            for atlas_info in self._atlases:
                parc_name = atlas_info["name"]
                labels = atlas_info["labels"]
                n_labels = len(labels)

                # Get atlas labels as ordered list
                atlas_labels = [f"region_{i}" for i in range(n_labels)]
                if labels is not None:
                    sorted_regions = sorted(labels.items())
                    atlas_labels = [name for region_id, name in sorted_regions]

                if self.compute_disconnectivity_matrix:
                    zero_matrix = np.zeros((n_labels, n_labels), dtype=np.float32)

                    # mask_connectivity_matrix
                    mask_conn_key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="mask_connectivity_matrix",
                    )
                    results[mask_conn_key] = ConnectivityMatrix(
                        name="mask_connectivity_matrix",
                        matrix=zero_matrix.copy(),
                        region_labels=atlas_labels,
                        matrix_type="structural",
                        metadata={
                            "atlas": parc_name,
                            "empty_mask": True,
                        },
                    )

                    # disconnectionpct
                    disconn_key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="disconnectionpct",
                    )
                    results[disconn_key] = ConnectivityMatrix(
                        name="disconnectionpct",
                        matrix=zero_matrix.copy(),
                        region_labels=atlas_labels,
                        matrix_type="structural",
                        metadata={
                            "atlas": parc_name,
                            "description": "Percentage of streamlines disconnected by mask",
                            "empty_mask": True,
                        },
                    )

                    # full_connectivity_matrix
                    full_conn_key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="full_connectivity_matrix",
                    )
                    results[full_conn_key] = ConnectivityMatrix(
                        name="full_connectivity_matrix",
                        matrix=zero_matrix.copy(),
                        region_labels=atlas_labels,
                        matrix_type="structural",
                        metadata={
                            "atlas": parc_name,
                            "description": "Full brain connectivity matrix (reference)",
                            "empty_mask": True,
                        },
                    )

                    # intact_connectivity_matrix
                    intact_conn_key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="intact_connectivity_matrix",
                    )
                    results[intact_conn_key] = ConnectivityMatrix(
                        name="intact_connectivity_matrix",
                        matrix=zero_matrix.copy(),
                        region_labels=atlas_labels,
                        matrix_type="structural",
                        metadata={
                            "atlas": parc_name,
                            "description": "Intact connectivity excluding disconnected streamlines",
                            "empty_mask": True,
                        },
                    )

                    # matrix_statistics
                    matrix_stats_key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="matrix_statistics",
                    )
                    results[matrix_stats_key] = ScalarMetric(
                        name="matrix_statistics",
                        data={
                            "n_parcels": n_labels,
                            "n_edges_total": 0,
                            "n_edges_affected": 0,
                            "percent_edges_affected": 0.0,
                            "mean_disconnection_percent": 0.0,
                            "max_disconnection_percent": 0.0,
                            "mean_degree_reduction": 0.0,
                            "max_degree_reduction": 0.0,
                            "most_affected_parcel": atlas_labels[0] if atlas_labels else "N/A",
                        },
                        metadata={
                            "atlas": parc_name,
                            "empty_mask": True,
                        },
                    )

                if self.compute_roi_disconnection:
                    zero_values = dict.fromkeys(atlas_labels, 0.0)
                    key = build_result_key(
                        atlas=parc_name,
                        source="StructuralNetworkMapping",
                        desc="roidisconnectionpct",
                    )
                    results[key] = ParcelData(
                        name="roidisconnectionpct",
                        data=zero_values,
                        region_labels=atlas_labels,
                        aggregation_method="percent",
                        metadata={
                            "atlas": parc_name,
                            "description": "Percentage of streamlines disconnected per ROI",
                            "unit": "percent",
                            "empty_mask": True,
                        },
                    )

        # Transform to input space if requested
        if self.return_in_input_space:
            results = self._transform_results_to_input_space(results, mask_data)

        return results

    def _compute_connectivity_matrices(
        self,
        mask_data: SubjectData,
        mask_tck_path: Path,
        temp_dir_path: Path,
        subject_id: str,
        atlas_info: dict | None = None,
    ) -> dict[str, "AnalysisResult"]:
        """Compute parcellated connectivity matrices for a single atlas.

        Parameters
        ----------
        mask_data : SubjectData
            Subject data with mask image
        mask_tck_path : Path
            Path to mask-filtered tractogram
        temp_dir_path : Path
            Temporary directory for intermediate files
        subject_id : str
            SubjectData identifier for file naming
        atlas_info : dict, optional
            Atlas info dict with keys: name, image, labels, resolved_path.
            If None, falls back to legacy self._parcellation_resolved.

        Returns
        -------
        dict[str, AnalysisResult]
            Dictionary containing per-atlas results keyed by BIDS-style names.
        """
        # Resolve atlas info (support legacy single-atlas callers)
        if atlas_info is None:
            atlas_info = {
                "name": self.parcellation_name,
                "labels": self._atlas_labels,
                "resolved_path": self._parcellation_resolved,
            }

        parc_name = atlas_info["name"]
        parc_labels = atlas_info["labels"]
        parc_resolved = atlas_info["resolved_path"]

        # Step 1: Compute full-brain connectivity matrix (cached per atlas)
        if parc_name not in self._full_connectivity_matrices:
            self.logger.info(
                f"Computing full-brain connectivity matrix for '{parc_name}' (will be cached)",
                indent_level=1,
            )
            self._full_connectivity_matrices[parc_name] = self._compute_connectivity_matrix(
                tractogram_path=self.tractogram_path,
                matrix_name=f"full_connectivity_{parc_name}",
                parcellation_path=parc_resolved,
            )
        else:
            self.logger.info(
                f"Using cached full-brain connectivity matrix for '{parc_name}'", indent_level=1
            )

        full_matrix = self._full_connectivity_matrices[parc_name]

        # Step 2: Compute mask connectivity matrix
        self.logger.info("Computing mask connectivity matrix", indent_level=1)
        mask_matrix = self._compute_connectivity_matrix(
            tractogram_path=mask_tck_path,
            matrix_name=f"{subject_id}_mask_connectivity_{parc_name}",
            parcellation_path=parc_resolved,
        )

        # Step 3: Compute disconnectivity percentage
        self.logger.info("Computing disconnectivity percentage", indent_level=1)
        with np.errstate(divide="ignore", invalid="ignore"):
            disconn_pct = (mask_matrix / full_matrix) * 100

        # Handle division by zero
        disconn_pct = np.nan_to_num(disconn_pct, nan=0.0, posinf=0.0, neginf=0.0)

        # Step 4: Optional - compute intact (post-disconnection) connectivity
        intact_matrix = None
        if self.compute_roi_disconnection:
            # The intact tractogram (tckedit -exclude) is atlas-independent,
            # so reuse if already computed for this subject
            intact_tck_path = temp_dir_path / f"{subject_id}_intact.tck"
            if not intact_tck_path.exists():
                self.logger.info("Computing intact (post-disconnection) tractogram", indent_level=1)
                exclude_mask_path = temp_dir_path / f"{subject_id}_exclude_mask.nii.gz"
                nib.save(mask_data.mask_img, exclude_mask_path)

                command = [
                    "tckedit",
                    str(self.tractogram_path),
                    str(intact_tck_path),
                    "-exclude",
                    str(exclude_mask_path),
                    "-force",
                ]
                command.extend(_get_nthreads_args(self.n_jobs))
                run_mrtrix_command(command, verbose=self.show_mrtrix_output)

            self.logger.info("Computing intact connectivity matrix", indent_level=1)
            intact_matrix = self._compute_connectivity_matrix(
                tractogram_path=intact_tck_path,
                matrix_name=f"{subject_id}_intact_connectivity_{parc_name}",
                parcellation_path=parc_resolved,
            )

        # Step 5: Compute summary statistics
        matrix_stats = self._compute_matrix_statistics(
            full_matrix=full_matrix,
            mask_matrix=mask_matrix,
            disconn_pct=disconn_pct,
            intact_matrix=intact_matrix,
        )

        # Build results dict
        results = {}

        # Get atlas labels for ConnectivityMatrixResult
        # Convert dict[int, str] to list[str] ordered by region ID
        atlas_labels = [f"region_{i}" for i in range(mask_matrix.shape[0])]
        if parc_labels is not None:
            sorted_regions = sorted(parc_labels.items())
            atlas_labels = [name for region_id, name in sorted_regions]

        # ConnectivityMatrixResult for mask connectivity
        mask_connectivity_result = ConnectivityMatrix(
            name="mask_connectivity_matrix",
            matrix=mask_matrix,
            region_labels=atlas_labels,
            matrix_type="structural",
            metadata={
                "atlas": parc_name,
                "tractogram": str(self.tractogram_path),
            },
        )
        mask_conn_key = build_result_key(
            atlas=parc_name,
            source="StructuralNetworkMapping",
            desc="mask_connectivity_matrix",
        )
        results[mask_conn_key] = mask_connectivity_result

        # ConnectivityMatrixResult for disconnection percentage
        disconn_result = ConnectivityMatrix(
            name="disconnectionpct",
            matrix=disconn_pct,
            region_labels=atlas_labels,
            matrix_type="structural",
            metadata={
                "atlas": parc_name,
                "description": "Percentage of streamlines disconnected by mask",
            },
        )
        disconn_pct_key = build_result_key(
            atlas=parc_name,
            source="StructuralNetworkMapping",
            desc="disconnectionpct",
        )
        results[disconn_pct_key] = disconn_result

        # ConnectivityMatrixResult for full connectivity (reference)
        full_connectivity_result = ConnectivityMatrix(
            name="full_connectivity_matrix",
            matrix=full_matrix,
            region_labels=atlas_labels,
            matrix_type="structural",
            metadata={
                "atlas": parc_name,
                "description": "Full brain connectivity matrix (reference)",
            },
        )
        full_conn_key = build_result_key(
            atlas=parc_name,
            source="StructuralNetworkMapping",
            desc="full_connectivity_matrix",
        )
        results[full_conn_key] = full_connectivity_result

        # Optional: intact (post-disconnection) connectivity matrix
        if intact_matrix is not None:
            intact_result = ConnectivityMatrix(
                name="intact_connectivity_matrix",
                matrix=intact_matrix,
                region_labels=atlas_labels,
                matrix_type="structural",
                metadata={
                    "atlas": parc_name,
                    "description": "Intact connectivity excluding disconnected streamlines",
                },
            )
            intact_conn_key = build_result_key(
                atlas=parc_name,
                source="StructuralNetworkMapping",
                desc="intact_connectivity_matrix",
            )
            results[intact_conn_key] = intact_result

            # Compute per-ROI disconnection percentage
            full_roi_streamlines = np.sum(full_matrix, axis=1)
            mask_roi_streamlines = np.sum(mask_matrix, axis=1)

            with np.errstate(divide="ignore", invalid="ignore"):
                roi_disconnection_pct = (mask_roi_streamlines / full_roi_streamlines) * 100

            roi_disconnection_pct = np.nan_to_num(
                roi_disconnection_pct, nan=0.0, posinf=0.0, neginf=0.0
            )

            roi_disconnection_data = {
                label: float(roi_disconnection_pct[i]) for i, label in enumerate(atlas_labels)
            }
            roi_disconnection_result = ParcelData(
                name="roidisconnectionpct",
                data=roi_disconnection_data,
                region_labels=atlas_labels,
                aggregation_method="percent",
                metadata={
                    "atlas": parc_name,
                    "description": "Percentage of streamlines disconnected per ROI",
                    "unit": "percent",
                },
            )
            roi_disconnection_key = build_result_key(
                atlas=parc_name,
                source="StructuralNetworkMapping",
                desc="roidisconnectionpct",
            )
            results[roi_disconnection_key] = roi_disconnection_result

        # MiscResult for matrix statistics
        stats_result = ScalarMetric(
            name="matrix_statistics",
            data=matrix_stats,
            metadata={
                "atlas": parc_name,
            },
        )
        matrix_stats_key = build_result_key(
            atlas=parc_name,
            source="StructuralNetworkMapping",
            desc="matrix_statistics",
        )
        results[matrix_stats_key] = stats_result

        return results

    def _compute_connectivity_matrix(
        self,
        tractogram_path: Path,
        matrix_name: str,
        parcellation_path: Path | None = None,
    ) -> np.ndarray:
        """Compute connectivity matrix from tractogram using tck2connectome.

        Parameters
        ----------
        tractogram_path : Path
            Path to tractogram file
        matrix_name : str
            Name for temporary CSV file
        parcellation_path : Path, optional
            Path to atlas NIfTI file. Defaults to self._parcellation_resolved.

        Returns
        -------
        np.ndarray
            Connectivity matrix (n_parcels x n_parcels)
        """
        from lacuna.utils.cache import make_temp_file

        if parcellation_path is None:
            parcellation_path = self._parcellation_resolved

        with make_temp_file(suffix=".csv", delete=True, mode="w") as tmp_csv:
            output_csv = Path(tmp_csv.name)

            command = [
                "tck2connectome",
                str(tractogram_path),
                str(parcellation_path),
                str(output_csv),
                "-symmetric",
                "-zero_diagonal",
                "-force",
            ]
            command.extend(_get_nthreads_args(self.n_jobs))

            run_mrtrix_command(command, verbose=self.show_mrtrix_output)

            # Load matrix
            matrix = np.loadtxt(output_csv, delimiter=",")

        return matrix

    def _compute_matrix_statistics(
        self,
        full_matrix: np.ndarray,
        mask_matrix: np.ndarray,
        disconn_pct: np.ndarray,
        intact_matrix: np.ndarray | None,
    ) -> dict:
        """Compute summary statistics for connectivity matrices.

        Parameters
        ----------
        full_matrix : np.ndarray
            Full connectivity matrix
        mask_matrix : np.ndarray
            Connectivity matrix from streamlines passing through mask
        disconn_pct : np.ndarray
            Disconnectivity percentage matrix
        intact_matrix : np.ndarray | None
            Intact connectivity matrix (excluding disconnected streamlines)

        Returns
        -------
        dict
            Summary statistics
        """
        # Edge-wise statistics
        n_edges = int(np.sum(full_matrix > 0))
        n_affected_edges = int(np.sum(mask_matrix > 0))
        mean_disconnection_pct = (
            float(np.mean(disconn_pct[full_matrix > 0])) if n_edges > 0 else 0.0
        )

        # Node-wise statistics (degree)
        full_degree = np.sum(full_matrix > 0, axis=1)
        mask_degree = np.sum(mask_matrix > 0, axis=1)
        degree_reduction = full_degree - mask_degree

        stats = {
            "n_parcels": full_matrix.shape[0],
            "n_edges_total": n_edges,
            "n_edges_affected": n_affected_edges,
            "percent_edges_affected": (
                float(n_affected_edges / n_edges * 100) if n_edges > 0 else 0.0
            ),
            "mean_disconnection_percent": mean_disconnection_pct,
            "max_disconnection_percent": float(np.max(disconn_pct)),
            "mean_degree_reduction": float(np.mean(degree_reduction)),
            "max_degree_reduction": int(np.max(degree_reduction)),
            "most_affected_parcel": int(np.argmax(degree_reduction)),
        }

        # Add intact matrix statistics if computed
        if intact_matrix is not None:
            intact_degree = np.sum(intact_matrix > 0, axis=1)
            stats["intact_mean_degree"] = float(np.mean(intact_degree))

            # Quality control: mask + intact should approximately equal full
            combined = mask_matrix + intact_matrix
            preservation = np.sum(combined > 0) / n_edges if n_edges > 0 else 0
            stats["connectivity_preservation_ratio"] = float(preservation)

        return stats

    def _get_version(self) -> str:
        """Get analysis version for provenance tracking."""
        from .. import __version__

        return __version__

    def _get_parameters(self) -> dict:
        """Get analysis parameters for provenance and display.

        Returns
        -------
        dict
            Dictionary of parameter names and values.
        """
        return {
            "connectome_name": self.connectome_name,
            "parcellation_name": str(self.parcellation_name) if self.parcellation_name else None,
            "compute_disconnectivity_matrix": self.compute_disconnectivity_matrix,
            "output_resolution": self.output_resolution,
            "n_jobs": self.n_jobs,
            "keep_intermediate": self.keep_intermediate,
            "show_mrtrix_output": self.show_mrtrix_output,
            "return_in_input_space": self.return_in_input_space,
            "verbose": self.verbose,
        }

    def _transform_results_to_input_space(self, results: dict, mask_data: SubjectData) -> dict:
        """Transform VoxelMap results back to original input mask space.

        Parameters
        ----------
        results : dict
            Dictionary of result objects
        mask_data : SubjectData
            Input mask data with space/resolution metadata. If the mask was
            transformed by BaseAnalysis, it will have _original_input_space
            and _original_input_resolution in metadata.

        Returns
        -------
        dict
            Results with transformed VoxelMap objects

        Raises
        ------
        ValueError
            If mask_data lacks space or resolution metadata

        Notes
        -----
        Transforms to original input space but keeps output_resolution for
        consistent output resolution. The output_resolution parameter controls
        both TDI generation and final output resolution.
        """
        from lacuna.core.spaces import REFERENCE_AFFINES, CoordinateSpace
        from lacuna.spatial.transform import transform_image

        # Get original input space from metadata (set by BaseAnalysis before transformation)
        # Fall back to current mask_data.space if not available (already in input space)
        original_space = mask_data.metadata.get("_original_input_space", mask_data.space)

        # Use output_resolution for final output (not input resolution)
        # This ensures consistent output resolution as specified by user
        target_resolution = self.output_resolution

        # Get reference affine for target space
        target_key = (original_space, target_resolution)
        if target_key not in REFERENCE_AFFINES:
            raise ValueError(
                f"No reference affine available for {original_space}@{target_resolution}mm. "
                f"Available spaces: {list(REFERENCE_AFFINES.keys())}"
            )

        target_space = CoordinateSpace(
            identifier=original_space,
            resolution=target_resolution,
            reference_affine=REFERENCE_AFFINES[target_key],
        )

        # Check if transformation is actually needed
        if original_space == self.TARGET_SPACE and target_resolution == self.output_resolution:
            return results

        self.logger.info(
            f"Transforming VoxelMap outputs from {self.TARGET_SPACE}@{self.output_resolution}mm "
            f"to {target_space.identifier}@{target_space.resolution}mm"
        )

        transformed_results = {}
        for key, result in results.items():
            # Only transform VoxelMap results
            from lacuna.core.data_types import VoxelMap

            if isinstance(result, VoxelMap):
                # Auto-detect interpolation method based on data type
                # Use nearest for binary maps (thresholdmaps), linear for continuous
                data = result.data.get_fdata()
                unique_vals = np.unique(data[~np.isnan(data)])
                is_binary = len(unique_vals) <= 2 and set(unique_vals).issubset({0, 1})
                interpolation = "nearest" if is_binary else "linear"

                # Transform the image
                transformed_img = transform_image(
                    img=result.data,
                    source_space=self.TARGET_SPACE,
                    target_space=target_space,
                    source_resolution=int(self.output_resolution),
                    interpolation=interpolation,
                    verbose=self.verbose,
                )

                # Create new VoxelMap with updated space
                transformed_result = VoxelMap(
                    name=result.name,
                    data=transformed_img,
                    space=target_space.identifier,
                    resolution=target_space.resolution,
                    metadata={
                        **result.metadata,
                        "transformed_from": f"{self.TARGET_SPACE}@{self.output_resolution}mm",
                        "transformed_to": f"{target_space.identifier}@{target_space.resolution}mm",
                    },
                )
                transformed_results[key] = transformed_result
            else:
                # Keep non-VoxelMap results as-is
                transformed_results[key] = result

        return transformed_results

__init__(connectome_name, parcellation_name=None, compute_disconnectivity_matrix=False, compute_roi_disconnection=False, include_streamline_counts=False, output_resolution=2, cache_tdi=True, n_jobs=1, keep_intermediate=False, check_dependencies=True, verbose=False, show_mrtrix_output=False, return_in_input_space=True)

Initialize StructuralNetworkMapping analysis.

Parameters:

Name	Type	Description	Default
connectome_name	str	Name of registered structural connectome (e.g., "HCP842_dTOR"). Use list_structural_connectomes() to see available options.	required
parcellation_name	str	Name of registered atlas for parcellated connectivity matrices.	None
compute_disconnectivity_matrix	bool	If True and parcellation_name provided, compute disconnectivity matrices (full, mask, and disconnectivity percentage connectivity matrices).	False
compute_roi_disconnection	bool	If True and parcellation_name provided, compute per-ROI disconnection percentages and intact (post-disconnection) connectivity matrix.	False
include_streamline_counts	bool	If True, include raw streamline count outputs (disconnection_tdi VoxelMap). By default, only percentage-based outputs are produced.	False
output_resolution	(1, 2)	Output resolution in mm (must match connectome resolution).	1
cache_tdi	bool	If True, cache computed TDI maps.	True
n_jobs	int	Number of threads for MRtrix3.	1
keep_intermediate	bool	If True, include intermediate results in output (lesion tractogram, lesion TDI, warped atlas). Useful for debugging and QC.	False
check_dependencies	bool	If True, checks for MRtrix3 availability.	True
verbose	bool	If True, print progress messages. If False, run silently.	True
show_mrtrix_output	bool	If True, display MRtrix3 command outputs. If False, suppress verbose MRtrix3 messages for cleaner output.	False
return_in_input_space	bool	If True, transform VoxelMap outputs back to the original input mask space. If False, outputs remain in the connectome space. Requires input SubjectData to have valid space/resolution metadata.	True

Raises:

Type	Description
MRtrixError	If MRtrix3 is not available and check_dependencies=True.
KeyError	If connectome_name not found in registry.
ValueError	If output_resolution is not 1 or 2.

Source code in src/lacuna/analysis/structural_network_mapping.py

def __init__(
    self,
    connectome_name: str,
    parcellation_name: str | list[str] | None = None,
    compute_disconnectivity_matrix: bool = False,
    compute_roi_disconnection: bool = False,
    include_streamline_counts: bool = False,
    output_resolution: Literal[1, 2] = 2,
    cache_tdi: bool = True,
    n_jobs: int = 1,
    keep_intermediate: bool = False,
    check_dependencies: bool = True,
    verbose: bool = False,
    show_mrtrix_output: bool = False,
    return_in_input_space: bool = True,
):
    """Initialize StructuralNetworkMapping analysis.

    Parameters
    ----------
    connectome_name : str
        Name of registered structural connectome (e.g., "HCP842_dTOR").
        Use list_structural_connectomes() to see available options.
    parcellation_name : str, optional
        Name of registered atlas for parcellated connectivity matrices.
    compute_disconnectivity_matrix : bool, default=False
        If True and parcellation_name provided, compute disconnectivity matrices
        (full, mask, and disconnectivity percentage connectivity matrices).
    compute_roi_disconnection : bool, default=False
        If True and parcellation_name provided, compute per-ROI disconnection
        percentages and intact (post-disconnection) connectivity matrix.
    include_streamline_counts : bool, default=False
        If True, include raw streamline count outputs (disconnection_tdi VoxelMap).
        By default, only percentage-based outputs are produced.
    output_resolution : {1, 2}, default=2
        Output resolution in mm (must match connectome resolution).
    cache_tdi : bool, default=True
        If True, cache computed TDI maps.
    n_jobs : int, default=1
        Number of threads for MRtrix3.
    keep_intermediate : bool, default=False
        If True, include intermediate results in output (lesion tractogram,
        lesion TDI, warped atlas). Useful for debugging and QC.
    check_dependencies : bool, default=True
        If True, checks for MRtrix3 availability.
    verbose : bool, default=True
        If True, print progress messages. If False, run silently.
    show_mrtrix_output : bool, default=False
        If True, display MRtrix3 command outputs. If False, suppress verbose
        MRtrix3 messages for cleaner output.
    return_in_input_space : bool, default=True
        If True, transform VoxelMap outputs back to the original input mask space.
        If False, outputs remain in the connectome space.
        Requires input SubjectData to have valid space/resolution metadata.

    Raises
    ------
    MRtrixError
        If MRtrix3 is not available and check_dependencies=True.
    KeyError
        If connectome_name not found in registry.
    ValueError
        If output_resolution is not 1 or 2.
    """
    super().__init__(verbose=verbose, keep_intermediate=keep_intermediate)

    # Validate output_resolution
    if output_resolution not in (1, 2):
        raise ValueError(f"output_resolution must be 1 or 2, got: {output_resolution}")

    # Load connectome from registry
    try:
        connectome = load_structural_connectome(connectome_name)
    except KeyError as e:
        available = [c.name for c in list_structural_connectomes()]
        raise KeyError(
            f"Connectome '{connectome_name}' not found in registry. "
            f"Available connectomes: {', '.join(available)}. "
            f"Use register_structural_connectome() to add new connectomes."
        ) from e

    # Store connectome information
    self.connectome_name = connectome_name
    self.tractogram_path = connectome.tractogram_path
    self.tractogram_space = connectome.metadata.space
    self.template = connectome.template_path  # May be None

    # Store analysis parameters — normalize parcellation_name(s) to list
    if parcellation_name is None:
        self.parcellation_names: list[str] = []
    elif isinstance(parcellation_name, str):
        self.parcellation_names = [parcellation_name]
    else:
        self.parcellation_names = list(parcellation_name)
    # Keep singular attribute for backward compatibility
    self.parcellation_name = self.parcellation_names[0] if self.parcellation_names else None

    self.compute_disconnectivity_matrix = compute_disconnectivity_matrix
    self.compute_roi_disconnection = compute_roi_disconnection
    self.include_streamline_counts = include_streamline_counts

    # Validate: compute flags require parcellation_name
    if (
        compute_disconnectivity_matrix or compute_roi_disconnection
    ) and not self.parcellation_names:
        flags = []
        if compute_disconnectivity_matrix:
            flags.append("compute_disconnectivity_matrix")
        if compute_roi_disconnection:
            flags.append("compute_roi_disconnection")
        raise ValueError(
            f"{' and '.join(flags)} requires parcellation_name to be specified "
            f"(via --parcel-atlases or --custom-parcellation). "
            f"Use list_parcellations() to see available options."
        )

    # Validate parcellation names early (full loading deferred to _validate_inputs)
    if self.parcellation_names:
        available_names = [a.name for a in list_parcellations()]
        for name in self.parcellation_names:
            if name not in available_names:
                raise ValueError(
                    f"Atlas '{name}' not found in registry. "
                    f"Available parcellations: {', '.join(available_names[:10])}... "
                    f"Use list_parcellations() to see all options."
                )

    self.output_resolution = output_resolution
    self.cache_tdi = cache_tdi
    self.n_jobs = n_jobs
    if n_jobs != -1 and n_jobs < 1:
        raise ValueError(f"n_jobs must be -1 (all CPUs) or >= 1, got {n_jobs}")
    self.keep_intermediate = keep_intermediate

    self.show_mrtrix_output = show_mrtrix_output
    self.return_in_input_space = return_in_input_space

    # Target space matches tractogram space, resolution from output_resolution
    # (used by BaseAnalysis._ensure_target_space)
    self.TARGET_SPACE = self.tractogram_space
    self.TARGET_RESOLUTION = self.output_resolution

    # Initialize logger
    self.logger = ConsoleLogger(verbose=verbose, width=70)

    # Internal state
    self.whole_brain_tdi = None  # Will be set during validation
    # Per-atlas state: list of dicts with keys: name, image, labels, resolved_path, was_transformed, etc.
    self._atlases: list[dict] = []
    self._full_connectivity_matrices: dict[str, np.ndarray] = {}  # Cached per atlas name
    # Keep legacy attributes for backward compat (set to first atlas after validation)
    self._atlas_image = None
    self._atlas_labels = None
    self._parcellation_resolved = None
    self._full_connectivity_matrix = None
    self._cached_tdi_path = None

    # Check MRtrix3 availability if requested
    if check_dependencies:
        try:
            check_mrtrix_available()
        except MRtrixError as e:
            raise MRtrixError(
                f"MRtrix3 is required for StructuralNetworkMapping but is not available.\n{e}"
            ) from e

run(mask_data)

Run structural network mapping analysis.

Automatically transforms lesion to tractogram space if needed.

Parameters:

Name	Type	Description	Default
mask_data	SubjectData	Lesion data to analyze (can be in any MNI152 space)	required

Returns:

Type	Description
SubjectData	Analysis results

Raises:

Type	Description
ValueError	If lesion is in native space (cannot transform)

Source code in src/lacuna/analysis/structural_network_mapping.py

def run(self, mask_data: SubjectData) -> SubjectData:
    """Run structural network mapping analysis.

    Automatically transforms lesion to tractogram space if needed.

    Parameters
    ----------
    mask_data : SubjectData
        Lesion data to analyze (can be in any MNI152 space)

    Returns
    -------
    SubjectData
        Analysis results

    Raises
    ------
    ValueError
        If lesion is in native space (cannot transform)
    """
    # Check for native space (cannot transform)
    space = mask_data.get_coordinate_space()

    if space.lower() == "native":
        raise ValueError(
            "Native space lesions are not supported for structural network mapping. "
            f"Lesions must be in a standard space. Tractogram is in {self.tractogram_space}."
        )

    # Transform to tractogram space (handled by base class)
    # The base class will handle space equivalence and transformations
    return super().run(mask_data)

get_analysis(name)

Get an analysis class by name.

Parameters:

Name	Type	Description	Default
name	str	Name of the analysis class (e.g., "FunctionalNetworkMapping").	required

Returns:

Type	Description
type[BaseAnalysis]	The analysis class.

Raises:

Type	Description
KeyError	If the analysis name is not found.

Examples:

>>> from lacuna.analysis import get_analysis
>>> FNM = get_analysis("FunctionalNetworkMapping")
>>> analysis = FNM(connectome_name="GSP1000")

Source code in src/lacuna/analysis/registry.py

def get_analysis(name: str) -> type[BaseAnalysis]:
    """
    Get an analysis class by name.

    Parameters
    ----------
    name : str
        Name of the analysis class (e.g., "FunctionalNetworkMapping").

    Returns
    -------
    type[BaseAnalysis]
        The analysis class.

    Raises
    ------
    KeyError
        If the analysis name is not found.

    Examples
    --------
    >>> from lacuna.analysis import get_analysis
    >>> FNM = get_analysis("FunctionalNetworkMapping")
    >>> analysis = FNM(connectome_name="GSP1000")
    """
    return AnalysisRegistry.get(name)

list_analyses()

List all available analysis classes.

This function returns all discovered analyses that: - Are subclasses of BaseAnalysis - Are concrete (not abstract) - Are not private (name doesn't start with '_')

Returns:

Type	Description
list of tuple	List of (name, class) tuples for all available analyses, sorted alphabetically by name.

Examples:

>>> from lacuna.analysis import list_analyses
>>> for name, cls in list_analyses():
...     print(f"{name}: {cls.batch_strategy}")
FunctionalNetworkMapping: vectorized
ParcelAggregation: parallel
RegionalDamage: parallel
StructuralNetworkMapping: parallel

Source code in src/lacuna/analysis/registry.py

def list_analyses() -> list[tuple[str, type[BaseAnalysis]]]:
    """
    List all available analysis classes.

    This function returns all discovered analyses that:
    - Are subclasses of BaseAnalysis
    - Are concrete (not abstract)
    - Are not private (name doesn't start with '_')

    Returns
    -------
    list of tuple
        List of (name, class) tuples for all available analyses,
        sorted alphabetically by name.

    Examples
    --------
    >>> from lacuna.analysis import list_analyses
    >>> for name, cls in list_analyses():
    ...     print(f"{name}: {cls.batch_strategy}")
    FunctionalNetworkMapping: vectorized
    ParcelAggregation: parallel
    RegionalDamage: parallel
    StructuralNetworkMapping: parallel
    """
    return AnalysisRegistry.list_analyses()