diff --git a/.github/workflows/testing.yml b/.github/workflows/testing.yml index c65b24d2..73c94a7c 100644 --- a/.github/workflows/testing.yml +++ b/.github/workflows/testing.yml @@ -5,7 +5,6 @@ # .. code-block:: bash # # $ nosetests --with-coverage --cover-package=brainprep --verbosity=2 -# --with-doctest --doctest-options='+ELLIPSIS,+NORMALIZE_WHITESPACE' ### name: "testing[nosetests]" @@ -48,7 +47,7 @@ jobs: python -m pip install --progress-bar off ".[ci]" - name: Run unit tests run: | - nosetests --with-coverage --cover-package=brainprep --verbosity=2 --with-doctest --doctest-options='+ELLIPSIS,+NORMALIZE_WHITESPACE' + nosetests --with-coverage --cover-package=brainprep --verbosity=2 - name: Coveralls if: matrix.python-version == 3.12 env: diff --git a/CHANGELOG.rst b/CHANGELOG.rst index d08e6154..8a537dfb 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -10,12 +10,12 @@ NEW --- - :bdg-success:`Enhancement` Add the dmriprep workflow. -- :bdg-success:`Enhancement` Add the mrophologist workflow. +- :bdg-success:`Enhancement` Add the morphologist workflow. Fixes ----- -- :bdg-danger:`Deprecation` Fix the containers that are using mri_synthstrip. +- :bdg-danger:`Deprecation` Fix the containers that are using `mri_synthstrip`. - :bdg-danger:`Deprecation` The run mapping file has been moved to avoid conflicts with FreeSurfer. @@ -24,10 +24,15 @@ Enhancements - :bdg-success:`Enhancement` Add signature hook. - :bdg-success:`Enhancement` Add live comand line monitoring support. +- :bdg-success:`Enhancement` Support multi-modality in Quasi-Raw workflow. +- :bdg-success:`Enhancement` Add `quick` mode in Quasi-Raw workflow. +- :bdg-success:`Enhancement` Support multi-modality in Deface workflow. Changes ------- +- :bdg-danger:`Deprecation` Optimize the Quasi-Raw workflow steps. + 2.0.0 ===== diff --git a/brainprep/decorators.py b/brainprep/decorators.py index 72e10f1e..f1227274 100644 --- a/brainprep/decorators.py +++ b/brainprep/decorators.py @@ -65,16 +65,6 @@ class Hook: ``before_call`` returns the inputs unchanged, and ``after_call`` returns the outputs unchanged. - Methods - ------- - before_call(func, inputs) - Hook executed before the wrapped function is called. - Must return a dictionary of (possibly modified) inputs. - - after_call(func, outputs) - Hook executed after the wrapped function returns. - Must return the (possibly modified) output value. - Notes ----- Subclasses may override one or both methods. If a method is not @@ -87,14 +77,22 @@ def before_call( func: Callable, inputs: dict[str, Any], ) -> dict[str, Any]: - """Transform and inspect inputs before the function call.""" + """ + Hook executed before the wrapped function is called. + Transform and/or inspect inputs. + Must return a dictionary of (possibly modified) inputs. + """ return inputs def after_call( self, outputs: Any, ) -> Any: - """Transform and inspect outputs after the function call.""" + """ + Hook executed after the wrapped function returns. + Transform and/or inspect outputs. + Must return the (possibly modified) output value. + """ return outputs diff --git a/brainprep/interfaces/__init__.py b/brainprep/interfaces/__init__.py index d06edf9f..2c9e4ef9 100644 --- a/brainprep/interfaces/__init__.py +++ b/brainprep/interfaces/__init__.py @@ -79,6 +79,7 @@ fmriprep_metrics, incremental_pca, mask_overlap, + maskdiff, mean_correlation, mriqc_metrics, network_entropy, @@ -91,7 +92,6 @@ from .utils import ( anonfile, copyfiles, - maskdiff, movedir, ungzfile, write_uuid_mapping, diff --git a/brainprep/interfaces/ants.py b/brainprep/interfaces/ants.py index 9088ff47..24b464c6 100644 --- a/brainprep/interfaces/ants.py +++ b/brainprep/interfaces/ants.py @@ -40,7 +40,8 @@ def biasfield( image_file: File, mask_file: File, output_dir: Directory, - entities: dict) -> tuple[list[str], tuple[File]]: + entities: dict, + quick: bool = False) -> tuple[list[str], tuple[File]]: """ Bias field correction of a BIDS-compliant anatomical image using ANTs's `N4BiasFieldCorrection`. @@ -55,6 +56,10 @@ def biasfield( Directory where the reoriented image will be saved. entities : dict A dictionary of parsed BIDS entities including modality. + quick : bool + Increased shrink factor from `1` to `4`, which downsamples the image + before estimating the bias field. + Default False. Returns ------- @@ -73,7 +78,7 @@ def biasfield( "N4BiasFieldCorrection", "-d", "3", "-i", str(image_file), - "-s", "1", + "-s", "4" if quick else "1", "-b", "[1x1x1,3]", "-c", "[50x50x50x50,0.001]", "-t", "[0.15,0.01,200]", diff --git a/brainprep/interfaces/freesurfer.py b/brainprep/interfaces/freesurfer.py index 81d33505..6b89412d 100644 --- a/brainprep/interfaces/freesurfer.py +++ b/brainprep/interfaces/freesurfer.py @@ -77,25 +77,29 @@ def brainmask( command : list[str] Skull-stripping command-line. outputs : tuple[File] - - mask_file : File - Skull-stripped brain image file. + -brain_file : File - Skull-stripped brain image file. + -mask_file : File - Binary brain mask image file. References ---------- .. footbibliography:: """ - basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}_brainmask".format( + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}".format( **entities) - mask_file = output_dir / f"{basename}.nii.gz" + brain_file = output_dir / f"{basename}_brain.nii.gz" + mask_file = output_dir / f"{basename}_brainmask.nii.gz" command = [ "mri_synthstrip", "-i", str(image_file), + "-o", str(brain_file), "-m", str(mask_file), + "-f", "0", "--no-csf", ] - return command, (mask_file, ) + return command, (brain_file, mask_file, ) @step( diff --git a/brainprep/interfaces/fsl.py b/brainprep/interfaces/fsl.py index 53290f54..a3324201 100644 --- a/brainprep/interfaces/fsl.py +++ b/brainprep/interfaces/fsl.py @@ -12,6 +12,7 @@ """ import os +from pathlib import Path from ..decorators import ( CoerceparamsHook, @@ -61,7 +62,7 @@ def reorient( outputs : tuple[File] - reorient_image_file : File - Reoriented input image file. """ - basename = "sub-{sub}_ses-{ses}_run-{run}_mod-T1w_reorient".format( + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}_reorient".format( **entities) reorient_image_file = output_dir / f"{basename}.nii.gz" @@ -109,7 +110,8 @@ def deface( outputs : tuple[File | list[File]] - deface_file : File - Defaced input T1w image file. - mask_file : File - Defacing binary mask. - - vol_files : list[File] - Defacing 3d rendering. + - transform_file : File - Affine transformation from original space to + MNI152 space. Raises ------ @@ -126,17 +128,29 @@ def deface( **entities) deface_file = output_dir / f"{basename}.nii.gz" mask_file = output_dir / f"{basename}mask.nii.gz" - - command = [ - "fsl_deface", - str(t1_file), - str(deface_file), - "-d", str(mask_file), - "-f", "0.5", - "-B", + transform_file = output_dir / f"{basename}affine.mat" + + resource_dir = Path(__file__).parent.parent / "resources" + bigfov_transfrom_file = resource_dir / "MNI_BigFov_to_MNI.mat" + + commands = [ + [ + "fsl_deface", + str(t1_file), + str(deface_file), + "-d", str(mask_file), + "-m13", str(transform_file), + "-f", "0.5", + "-B", + ], + [ + "convert_xfm", + "-omat", str(transform_file), + "-concat", str(bigfov_transfrom_file), str(transform_file), + ] ] - return command, (deface_file, mask_file, ) + return commands, (deface_file, mask_file, transform_file) @step( @@ -206,7 +220,8 @@ def scale( image_file: File, scale: int, output_dir: Directory, - entities: dict) -> tuple[list[str], tuple[File]]: + entities: dict, + interpolation: str = "spline") -> tuple[list[str], tuple[File]]: """ Apply an isotropic resampling transformation to a BIDS-compliant image file using FSL's `flirt`. @@ -221,6 +236,10 @@ def scale( Directory where the scaled image will be saved. entities : dict A dictionary of parsed BIDS entities including modality. + interpolation: str + The interpolation method: 'trilinear', 'nearestneighbour', 'sinc', or + 'spline'. + Default 'spline'. Returns ------- @@ -240,6 +259,7 @@ def scale( "-in", str(image_file), "-ref", str(image_file), "-applyisoxfm", str(scale), + "-interp", interpolation, "-out", str(scaled_anatomical_file), "-omat", str(transform_file), "-verbose", "1", @@ -263,7 +283,9 @@ def affine( anatomical_file: File, template_file: File, output_dir: Directory, - entities: dict) -> tuple[list[str], tuple[File]]: + entities: dict, + rigid: bool = False, + quick: bool = False) -> tuple[list[str], tuple[File]]: """ Affinely register a BIDS-compliant anatomical image to a template file using FSL's `flirt`. @@ -278,6 +300,15 @@ def affine( Directory where the affine transformation will be saved. entities : dict A dictionary of parsed BIDS entities including modality. + rigid : bool + Estimate a 6 DOF transformation that maintains the original size and + shape of the brain. By default a 9 DOF transformation allows for + additional scaling in the x, y, and z directions, adjusting the size + and shape of the brain during the alignment process. + Default False. + quick : bool + Restricted rotation search range to +/-30° on all three axes. + Default False. Returns ------- @@ -301,11 +332,17 @@ def affine( "-anglerep", "euler", "-bins", "256", "-interp", "trilinear", - "-dof", "9", + "-dof", "6" if rigid else "9", "-out", str(aligned_anatomical_file), "-omat", str(transform_file), "-verbose", "1" ] + if quick: + command += [ + "-searchrx", "-30", "30", + "-searchry", "-30", "30", + "-searchrz", "-30", "30", + ] return command, (aligned_anatomical_file, transform_file) @@ -346,7 +383,8 @@ def applyaffine( A dictionary of parsed BIDS entities including modality. interpolation: str The interpolation method: 'trilinear', 'nearestneighbour', 'sinc', or - 'spline'. Default 'spline'. + 'spline'. + Default 'spline'. Returns ------- @@ -364,7 +402,7 @@ def applyaffine( "-in", str(image_file), "-ref", str(template_file), "-init", str(transform_file), - "-interp", str(interpolation), + "-interp", interpolation, "-applyxfm", "-out", str(aligned_image_file), ] diff --git a/brainprep/interfaces/plotting.py b/brainprep/interfaces/plotting.py index afd6ff44..ea059c58 100644 --- a/brainprep/interfaces/plotting.py +++ b/brainprep/interfaces/plotting.py @@ -13,6 +13,7 @@ import itertools import warnings +import matplotlib.lines as mlines import matplotlib.pyplot as plt import nibabel import numpy as np @@ -145,7 +146,7 @@ def plot_defacing_mosaic( mosaic_file : File Path to the saved mosaic image. """ - basename = "sub-{sub}_ses-{ses}_run-{run}_mod-T1w_deface".format( + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}_deface".format( **entities) mosaic_file = output_dir / f"{basename}mosaic.png" @@ -193,6 +194,7 @@ def plot_histogram( col_name: str, output_dir: Directory, bar_coords: list[float] | None = None, + suffix: str | None = None, dryrun: bool = False) -> tuple[File]: """ Generates a histogram image with optional vertical bars. @@ -206,7 +208,11 @@ def plot_histogram( output_dir : Directory Directory where the image with the histogram will be saved. bar_coords: list[float] | None - Coordianates of vertical lines to be displayed in red. Default None. + Coordianates of vertical lines to be displayed in red. + Default None. + suffix : str | None + Suffix added to the generated PNG file. + Default None. dryrun : bool If True, skip actual computation and file writing. Default False. @@ -215,7 +221,7 @@ def plot_histogram( histogram_file : File Generated image with the histogram. """ - histogram_file = output_dir / f"histogram_{col_name}.png" + histogram_file = output_dir / f"histogram_{col_name}{suffix or ''}.png" if dryrun: return (histogram_file, ) @@ -242,7 +248,6 @@ def plot_histogram( ax.axvline(x=x_coord, color="red") ax.spines["right"].set_visible(False) ax.spines["top"].set_visible(False) - ax.legend() plt.savefig(histogram_file) @@ -362,6 +367,7 @@ def plot_brainparc( def plot_pca( pca_file: File, output_dir: Directory, + suffix: str | None = None, dryrun: bool = False) -> tuple[File]: """ Plot the two first PCA components. @@ -374,6 +380,9 @@ def plot_pca( and ``run``. output_dir : Directory Directory where the result image will be saved. + suffix : str | None + Suffix added to the generated PNG file. + Default None. dryrun : bool If True, skip actual computation and file writing. Default False. @@ -382,7 +391,7 @@ def plot_pca( pca_image_file : File Generated image with the two first PCA components. """ - pca_image_file = output_dir / f"pca.png" + pca_image_file = output_dir / f"pca{suffix or ''}.png" if dryrun: return (pca_image_file, ) @@ -391,15 +400,35 @@ def plot_pca( fig, ax = plt.subplots(figsize=(20, 10)) ax.scatter(df.pc1, df.pc2) - for idx in range(len(df)): - ax.annotate( - f"{df.participant_id[idx]}-{df.session[idx]}-{df.run[idx]}", - xy=(df.pc1[idx], df.pc2[idx]), + df.apply( + lambda row: ax.annotate( + f"{row.participant_id}-{row.session}-{row.run}", + xy=(row.pc1, row.pc2), xytext=(4, 4), - textcoords="offset pixels" - ) - plt.xlabel(f"PC1 (var={df.explained_variance_ratio_pc1[0]:.2f})") - plt.ylabel(f"PC2 (var={df.explained_variance_ratio_pc2[1]:.2f})") + textcoords="offset pixels", + fontsize=9, + ), + axis=1, + ) + annotation_desc = mlines.Line2D( + [], + [], + color="none", + label="Participant - Session - Run", + ) + ax.legend( + handles=[annotation_desc], + loc="upper right", + frameon=True, + facecolor="#f9f9f9", + edgecolor="gray", + ) + plt.xlabel( + fr"$\mathbf{{PC1}}$ (var={df.explained_variance_ratio_pc1[0]:.2f})" + ) + plt.ylabel( + fr"$\mathbf{{PC2}}$ (var={df.explained_variance_ratio_pc2[1]:.2f})" + ) plt.axis("equal") ax.spines["right"].set_visible(False) ax.spines["top"].set_visible(False) diff --git a/brainprep/interfaces/qualcheck.py b/brainprep/interfaces/qualcheck.py index fe84a9ad..0a750f84 100644 --- a/brainprep/interfaces/qualcheck.py +++ b/brainprep/interfaces/qualcheck.py @@ -160,6 +160,7 @@ def mask_overlap( maskdiff_files_regex: str, output_dir: Directory, overlap_threshold: float = 0.05, + suffix: str | None = None, dryrun: bool = False) -> tuple[File]: """ Compute overlap ratios between mask pairs from `maskdiff` summary files. @@ -180,9 +181,14 @@ def mask_overlap( Directory where a TSV file containing the mean correlation values is created. overlap_threshold : float - Quality control threshold applied on the overalp score. Default 0.05. + Quality control threshold applied on the overalp score. + Default 0.05. + suffix : str | None + Suffix added to the generated TSV file. + Default None. dryrun : bool - If True, skip actual computation and file writing. Default False. + If True, skip actual computation and file writing. + Default False. Returns ------- @@ -199,7 +205,7 @@ def mask_overlap( indicating whether the overlap score do not exceeds the threshold: ``qc = 1`` if ``overlap < overlap_threshold``, otherwise ``qc = 0``. """ - overlap_file = output_dir / "mask_overlap.tsv" + overlap_file = output_dir / f"mask_overlap{suffix or ''}.tsv" if dryrun: return (overlap_file, ) @@ -265,29 +271,44 @@ def mask_overlap( ] ) def mean_correlation( - image_files_regex: str, - atlas_file: File, + data_files_regex: str, + atlas_file: File | None, output_dir: Directory, - correlation_threshold: float = 0.5, + correlation_threshold: float | None = 0.5, + suffix: str | None = None, dryrun: bool = False) -> tuple[File]: """ Compute the mean Pearson correlation between a reference image and a list - of other images. + of input images. + + It can use pre-computed correlation data for each subject or compute the + correlation internally. + If the correlaton is computed internally, individual pre-computed + correlation data are saved in the ``subjects`` directory. + Additionally, it performs quality control based on a specified correlation + threshold. Parameters ---------- - image_files_regex : str + data_files_regex : str A REGEX to image files, each representing an image of the same shape - and geometry as `atlas_file`. - atlas_file : File - An file representing the reference image. + and geometry as `atlas_file` or pre-computed TSV correlation data. + atlas_file : File | None + A file representing the reference image. If None, expect pre-computed + correlation data. output_dir : Directory Directory where a TSV file containing the mean correlation values is created. - correlation_threshold : float - Quality control threshold on the correlation score. Default 0.5. + correlation_threshold : float | None + Quality control threshold on the correlation score. + If None do not add the ``qc`` column. + Default 0.5. + suffix : str | None + Suffix added to the generated TSV file. + Default None. dryrun : bool - If True, skip actual computation and file writing. Default False. + If True, skip actual computation and file writing. + Default False. Returns ------- @@ -302,6 +323,7 @@ def mean_correlation( ------ ValueError If the atlas and an image have incompatible shape or geometry. + If invalid pre-computed correlation data are provided. Notes ----- @@ -310,17 +332,18 @@ def mean_correlation( ``qc = 1`` if ``mean_correlation > correlation_threshold``, otherwise ``qc = 0``. """ - correlations_file = output_dir / "mean_correlations.tsv" + correlations_file = output_dir / f"mean_correlations{suffix or ''}.tsv" if dryrun: return (correlations_file, ) image_files = coerce_to_path( - glob.glob(str(image_files_regex)), + glob.glob(str(data_files_regex)), expected_type=list[File], ) - atlas_im = nibabel.load(atlas_file) - atlas_arr = atlas_im.get_fdata() + if atlas_file is not None: + atlas_im = nibabel.load(atlas_file) + atlas_arr = atlas_im.get_fdata() scores = pd.DataFrame( columns=( @@ -332,30 +355,55 @@ def mean_correlation( ) for path in image_files: entities = parse_bids_keys(path) - im = nibabel.load(path) - arr = atlas_im.get_fdata() - if atlas_arr.shape != arr.shape: - raise ValueError( - f"Atlas and image have incompatible shape: {path}" - ) - if not np.allclose(atlas_im.affine, im.affine): - raise ValueError( - f"Atlas and image have incompatible orientation: {path}" + if atlas_file is not None: + im = nibabel.load(path) + arr = atlas_im.get_fdata() + if atlas_arr.shape != arr.shape: + raise ValueError( + f"Atlas and image have incompatible shape: {path}" + ) + if not np.allclose(atlas_im.affine, im.affine): + raise ValueError( + f"Atlas and image have incompatible orientation: {path}" + ) + corr, _ = pearsonr( + atlas_arr.flatten(), + arr.flatten(), ) - corr, _ = pearsonr( - atlas_arr.flatten(), - arr.flatten(), - ) + else: + df_ = pd.read_csv(path, sep="\t") + if len(df_) != 1: + raise ValueError( + f"Invalid pre-computed correlation data: {path}" + ) + corr = df_.iloc[0]["mean_correlation"] scores.loc[len(scores)] = [ entities["sub"], entities["ses"], entities["run"], - corr, + float(corr), ] + if atlas_file is not None: + individual_score = scores.iloc[[-1]] + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}".format( + **entities + ) + individual_score_file = ( + path.parent / + "quality_check" / + f"{basename}_corr.tsv" + ) + individual_score_file.parent.mkdir(parents=True, exist_ok=True) + individual_score.to_csv( + individual_score_file, + sep="\t", + index=False, + ) - scores["qc"] = ( - scores["mean_correlation"] > correlation_threshold - ).astype(int) + if correlation_threshold is not None: + scores["qc"] = ( + scores["mean_correlation"] > correlation_threshold + ).astype(int) scores = scores.sort_values(by=["participant_id", "session", "run"]) scores.to_csv( correlations_file, @@ -366,6 +414,114 @@ def mean_correlation( return (correlations_file, ) +@step( + hooks=[ + CoerceparamsHook(), + OutputdirHook( + quality_check=True + ), + LogRuntimeHook( + bunched=False + ), + PythonWrapperHook(), + SignatureHook(), + ] +) +def maskdiff( + mask1_file: File, + mask2_file: File, + output_dir: Directory, + entities: dict, + inv_mask1: bool = False, + inv_mask2: bool = False, + dryrun: bool = False) -> tuple[File]: + """ + Compute summary statistics comparing two binary masks. + + This function loads two binary mask images, verifies that they share + the same spatial dimensions and affine transformation, computes their + voxel-wise intersection, and writes a summary table containing voxel + counts and physical volumes (in mm³) for each mask and their intersection. + + Parameters + ---------- + mask1_file : File + Path to the first binary mask image. + mask2_file : File + Path to the second binary mask image. + output_dir : Directory + Directory where the defacing mask will be saved. + entities : dict + A dictionary of parsed BIDS entities including modality. + inv_mask1 : bool + If True, the first mask is inverted before comparison. This is + useful when the mask represents an exclusion region rather than an + inclusion region. Default False. + inv_mask2 : bool + If True, the second mask is inverted before comparison. This is + useful when the mask represents an exclusion region rather than an + inclusion region. Default False. + dryrun : bool + If True, skip actual computation and file writing. Default False. + + Returns + ------- + summary_file : File + Path to the generated summary TSV file. + + Raises + ------ + ValueError + If both masks have not identical shapes and affines. + """ + basename = "sub-{sub}_ses-{ses}_run-{run}_mod-{mod}_maskdiff".format( + **entities) + summary_file = output_dir / f"{basename}.tsv" + + if not dryrun: + + mask1_im = nibabel.load(mask1_file) + mask2_im = nibabel.load(mask2_file) + mask1 = mask1_im.get_fdata().astype(bool) + mask2 = mask2_im.get_fdata().astype(bool) + + if inv_mask1: + mask1 = ~mask1 + if inv_mask2: + mask1 = ~mask2 + + if mask1.shape != mask2.shape: + raise ValueError( + f"Mask shapes differ: {mask1.shape} vs {mask2.shape}. " + "Resampling is required." + ) + if not np.allclose(mask1_im.affine, mask2_im.affine): + raise ValueError( + "Mask affines differ. Resampling is required before " + "intersection." + ) + + intersection = np.logical_and(mask1, mask2) + voxel_volume = np.abs(np.linalg.det(mask1_im.affine[:3, :3])) + + summary_df = pd.DataFrame({ + "mask": ["mask1", "mask2", "intersection"], + "voxels": [ + mask1.sum(), + mask2.sum(), + intersection.sum(), + ], + "volume_mm3": [ + mask1.sum() * voxel_volume, + mask2.sum() * voxel_volume, + intersection.sum() * voxel_volume, + ] + }) + summary_df.to_csv(summary_file, sep="\t", index=False) + + return (summary_file, ) + + @step( hooks=[ CoerceparamsHook(), @@ -384,6 +540,7 @@ def incremental_pca( image_files_regex: str, output_dir: Directory, batch_size: int = 10, + suffix: str | None = None, dryrun: bool = False) -> tuple[File]: """ Perform an Incremental PCA with 2 components on a collection of images @@ -410,6 +567,9 @@ def incremental_pca( batch_size : int Number of images to use in each batch. If None, a single batch is used. Default is 10. + suffix : str | None + Suffix added to the generated TSV file. + Default None. dryrun : bool If True, skip actual computation and file writing. Default False. @@ -425,7 +585,7 @@ def incremental_pca( If the dataset contains fewer than 2 images, which prevents PCA computation. """ - pca_file = output_dir / "pca.tsv" + pca_file = output_dir / f"pca{suffix or ''}.tsv" if dryrun: return (pca_file, ) diff --git a/brainprep/interfaces/utils.py b/brainprep/interfaces/utils.py index 61dbbe40..5598f047 100644 --- a/brainprep/interfaces/utils.py +++ b/brainprep/interfaces/utils.py @@ -15,8 +15,6 @@ import shutil import socket -import nibabel -import numpy as np import pandas as pd from ..decorators import ( @@ -37,112 +35,6 @@ ) -@step( - hooks=[ - CoerceparamsHook(), - OutputdirHook(), - LogRuntimeHook( - bunched=False - ), - PythonWrapperHook(), - SignatureHook(), - ] -) -def maskdiff( - mask1_file: File, - mask2_file: File, - output_dir: Directory, - entities: dict, - inv_mask1: bool = False, - inv_mask2: bool = False, - dryrun: bool = False) -> tuple[File]: - """ - Compute summary statistics comparing two binary masks. - - This function loads two binary mask images, verifies that they share - the same spatial dimensions and affine transformation, computes their - voxel-wise intersection, and writes a summary table containing voxel - counts and physical volumes (in mm³) for each mask and their intersection. - - Parameters - ---------- - mask1_file : File - Path to the first binary mask image. - mask2_file : File - Path to the second binary mask image. - output_dir : Directory - Directory where the defacing mask will be saved. - entities : dict - A dictionary of parsed BIDS entities including modality. - inv_mask1 : bool - If True, the first mask is inverted before comparison. This is - useful when the mask represents an exclusion region rather than an - inclusion region. Default False. - inv_mask2 : bool - If True, the second mask is inverted before comparison. This is - useful when the mask represents an exclusion region rather than an - inclusion region. Default False. - dryrun : bool - If True, skip actual computation and file writing. Default False. - - Returns - ------- - summary_file : File - Path to the generated summary TSV file. - - Raises - ------ - ValueError - If both masks have not identical shapes and affines. - """ - basename = "sub-{sub}_ses-{ses}_run-{run}_mod-T1w_defacemask".format( - **entities) - summary_file = output_dir / f"{basename}.tsv" - - if not dryrun: - - mask1_im = nibabel.load(mask1_file) - mask2_im = nibabel.load(mask2_file) - mask1 = mask1_im.get_fdata().astype(bool) - mask2 = mask2_im.get_fdata().astype(bool) - - if inv_mask1: - mask1 = ~mask1 - if inv_mask2: - mask1 = ~mask2 - - if mask1.shape != mask2.shape: - raise ValueError( - f"Mask shapes differ: {mask1.shape} vs {mask2.shape}. " - "Resampling is required." - ) - if not np.allclose(mask1_im.affine, mask2_im.affine): - raise ValueError( - "Mask affines differ. Resampling is required before " - "intersection." - ) - - intersection = np.logical_and(mask1, mask2) - voxel_volume = np.abs(np.linalg.det(mask1_im.affine[:3, :3])) - - summary_df = pd.DataFrame({ - "mask": ["mask1", "mask2", "intersection"], - "voxels": [ - mask1.sum(), - mask2.sum(), - intersection.sum(), - ], - "volume_mm3": [ - mask1.sum() * voxel_volume, - mask2.sum() * voxel_volume, - intersection.sum() * voxel_volume, - ] - }) - summary_df.to_csv(summary_file, sep="\t", index=False) - - return (summary_file, ) - - @step( hooks=[ CoerceparamsHook(), @@ -155,29 +47,38 @@ def maskdiff( ] ) def copyfiles( - source_image_files: list[File], - destination_image_files: list[File], + source_files: list[File], + destination_files: list[File], output_dir: Directory, + move_files: bool = False, dryrun: bool = False) -> None: """ - Copy input image files. + Copy or move input files to a specified destination. Parameters ---------- - source_image_files : list[File] - Path to the image to be copied. - destination_image_files : list[File] - Path to the locations where images will be copied. + source_files : list[File] + List of files to be copied or moved. + destination_files : list[File] + List of files representing the target locations for the copied or + moved files. output_dir : Directory - Directory where the images are copied. + The directory where the files will be copied or moved to. + move_files : bool + If True, move the input files instead of copying them. + Default False. dryrun : bool - If True, skip actual computation and file writing. Default False. + If True, skip actual computation and file writing. + Default False. """ if not dryrun: - for src_path, dest_path in zip(source_image_files, - destination_image_files, + for src_path, dest_path in zip(source_files, + destination_files, strict=True): - shutil.copy(src_path, dest_path) + if move_files: + shutil.move(src_path, dest_path) + else: + shutil.copy(src_path, dest_path) @step( diff --git a/brainprep/reporting/rst_reporting.py b/brainprep/reporting/rst_reporting.py index dbea8a79..d47754ff 100644 --- a/brainprep/reporting/rst_reporting.py +++ b/brainprep/reporting/rst_reporting.py @@ -62,6 +62,7 @@ class SingletonReport(type): >>> class Report(metaclass=SingletonReport): ... def __init__(self): ... self._registry = {} + ... self._commands = {} >>> r1 = Report() >>> r2 = Report() @@ -100,6 +101,7 @@ def __call__( if not is_reloadable: inst._count = 0 inst._registry.clear() + inst._commands.clear() if is_increment: inst._count += 1 inst._reloadable = is_reloadable @@ -214,8 +216,14 @@ def register( if identifier not in self._registry: self._registry[identifier] = Bunch() if name in self._registry[identifier]: + items_str = [ + f"- {name_}\n" + for name_ in self._registry[identifier] + ] raise ValueError( - "Duplicated name in registry." + f"Duplicated name in registry: {name}\n" + f">> {identifier}\n" + f"{''.join(items_str)}" ) if not (isinstance(data, Bunch) or (isinstance(data, str) and name in self._str_fields) diff --git a/brainprep/resources/MNI152_T1_1mm_brain.nii.gz b/brainprep/resources/MNI152_T1_1mm_brain.nii.gz index 3bd6af20..0bd66f5c 100644 Binary files a/brainprep/resources/MNI152_T1_1mm_brain.nii.gz and b/brainprep/resources/MNI152_T1_1mm_brain.nii.gz differ diff --git a/brainprep/resources/MNI152_T1_2mm_brain.nii.gz b/brainprep/resources/MNI152_T1_2mm_brain.nii.gz new file mode 100755 index 00000000..3fda6c09 Binary files /dev/null and b/brainprep/resources/MNI152_T1_2mm_brain.nii.gz differ diff --git a/brainprep/resources/MNI152_T2_1mm_brain.nii.gz b/brainprep/resources/MNI152_T2_1mm_brain.nii.gz new file mode 100644 index 00000000..8f61c183 Binary files /dev/null and b/brainprep/resources/MNI152_T2_1mm_brain.nii.gz differ diff --git a/brainprep/resources/MNI152_T2_2mm_brain.nii.gz b/brainprep/resources/MNI152_T2_2mm_brain.nii.gz new file mode 100755 index 00000000..84dcb0e3 Binary files /dev/null and b/brainprep/resources/MNI152_T2_2mm_brain.nii.gz differ diff --git a/brainprep/resources/MNI_BigFov_to_MNI.mat b/brainprep/resources/MNI_BigFov_to_MNI.mat new file mode 100644 index 00000000..fcb756da --- /dev/null +++ b/brainprep/resources/MNI_BigFov_to_MNI.mat @@ -0,0 +1,4 @@ +1 -0 -0 -50 +0 1 -0 -50 +0 0 1 -100 +0 0 0 1 diff --git a/brainprep/resources/MNI_to_MNI_BigFoV.mat b/brainprep/resources/MNI_to_MNI_BigFoV.mat new file mode 100644 index 00000000..063ab124 --- /dev/null +++ b/brainprep/resources/MNI_to_MNI_BigFoV.mat @@ -0,0 +1,4 @@ +1 0 0 50 +0 1 0 50 +0 0 1 100 +0 0 0 1 diff --git a/brainprep/tests/test_docstring.py b/brainprep/tests/test_docstring.py index 3a59fb4b..4c2e015d 100644 --- a/brainprep/tests/test_docstring.py +++ b/brainprep/tests/test_docstring.py @@ -12,19 +12,36 @@ import pkgutil import unittest +from brainprep.reporting import RSTReport -def load_tests(loader, tests, ignore): - for _, module_name, ispkg in pkgutil.walk_packages( - brainprep.__path__, - brainprep.__name__ + "."): - module = importlib.import_module(module_name) - tests.addTests( - doctest.DocTestSuite( + +def doctest_setup(test): + test.globs["report"] = RSTReport() + + +class TestDocString(unittest.TestCase): + + def test_doctests(self): + result = unittest.TestResult() + n_tests = 0 + for _, module_name, ispkg in pkgutil.walk_packages( + brainprep.__path__, + brainprep.__name__ + "."): + module = importlib.import_module(module_name) + suite = doctest.DocTestSuite( module, + setUp=doctest_setup, optionflags=doctest.NORMALIZE_WHITESPACE | doctest.ELLIPSIS ) - ) - return tests + suite(result) + n_tests += 1 + if not result.wasSuccessful(): + report = "" + n_errors = 0 + for test, err in result.failures + result.errors: + report += f"\nTest fail: {test}\n>> {err}" + n_errors += 1 + self.fail(f"Error in doctests: {n_errors}/{n_tests}\n{report}") if __name__ == "__main__": diff --git a/brainprep/tests/test_workflow.py b/brainprep/tests/test_workflow.py index 2d5765b5..af84b542 100644 --- a/brainprep/tests/test_workflow.py +++ b/brainprep/tests/test_workflow.py @@ -13,12 +13,15 @@ import runpy from pathlib import Path +from brainprep.reporting import RSTReport + class TestGalleryExamples(unittest.TestCase): def setUp(self, test_interfaces=True): self.test_interfaces = test_interfaces self.examples_dir = Path(__file__).parent.parent.parent / "examples" + self.report = RSTReport() @staticmethod def run_cmd(cmd): @@ -31,24 +34,25 @@ def run_cmd(cmd): ) return None except subprocess.CalledProcessError as e: - return f"Command failed: {cmd}" + return f"Command failed: {' '.join(cmd)}" def _test_interface_commands(self, env): if not self.test_interfaces: return outdir = Path(env["outdir"]) - commands = [] + commands, commands_files = [], [] for commands_file in outdir.rglob("commands_*.rst"): commands.extend( commands_file.read_text().splitlines() ) - commands = [cmd.split(" ") for cmd in commands] - print(f"Parsed: {outdir}") + commands_files.append(f"\n - {commands_file}") + commands = [[*cmd.split(" "), "--dryrun"] for cmd in commands] + print(f"Parsed: {''.join(commands_files)}") print(f"Interface commands: {len(commands)}") failures = [] - with ProcessPoolExecutor(max_workers=50) as pool: + with ProcessPoolExecutor(max_workers=20) as pool: for msg in pool.map(TestGalleryExamples.run_cmd, commands): if msg is not None: failures.append(msg) @@ -78,7 +82,7 @@ def test_quality_assurance(self): "plot_quality_assurance.py" ) env = runpy.run_path(str(script_path)) - self._test_interface_commands(env) + # self._test_interface_commands(env) def test_defacing(self): script_path = ( diff --git a/brainprep/utils/color.py b/brainprep/utils/color.py index 32907a01..deb15dbc 100644 --- a/brainprep/utils/color.py +++ b/brainprep/utils/color.py @@ -94,8 +94,9 @@ def print_color(color: str, text: str, end: str = "\n") -> None: Name of the color style to apply. text : str The message to display. - end : str, optional - String appended after the message. Default ``"\\n"``. + end : str + String appended after the message. + Default ``"\\n"``. Notes ----- diff --git a/brainprep/workflow/defacing.py b/brainprep/workflow/defacing.py index e786d4a8..4ea2376d 100644 --- a/brainprep/workflow/defacing.py +++ b/brainprep/workflow/defacing.py @@ -11,6 +11,7 @@ """ import shutil +from pathlib import Path import brainprep.interfaces as interfaces @@ -37,7 +38,7 @@ CoerceparamsHook(), BidsHook( process="defacing", - bids_file="t1_file", + bids_file="anatomical_file", add_subjects=True, container="neurospin/brainprep-deface" ), @@ -49,25 +50,29 @@ ] ) def brainprep_defacing( - t1_file: File, + anatomical_file: File, output_dir: Directory, keep_intermediate: bool = False, **kwargs: dict) -> Bunch: """ - Defacing pre-processing workflow for anatomical T1-weighted images. + Defacing pre-processing workflow for anatomical images. Applies FSL's `fsl_deface` tool :footcite:p:`almagro2018deface` with - default settings to remove facial features (face and ears) from the input - image. This includes: + default settings to remove facial features (face and ears) from an input + T1-weighted MRI image. Apply defacing mask to T2-weighted or FLAIR MRI + images. This includes: - 1) Reorient the T1w image to standard MNI152 template space. - 2) Deface the T1w image. - 3) Generate a mosaic image of the defaced T1w image. + 1) Reorient the anatomical image to standard MNI152 template space. + 2) Compute a brain mask using a skull-stripping tool. + 3) Deface the T1w image or apply defacing to T2w and FLAIR images using + coregistration. + 4) Compute brain mask and defacing mask intersection. + 5) Generate a mosaic image of the defaced anatomical image. Parameters ---------- - t1_file : File - Path to the input T1w anatomical image file. + anatomical_file : File + Path to the input image file: T1w, T2w or FLAIR. output_dir : Directory Directory where the defaced image and related outputs will be saved (i.e., the root of your dataset). @@ -83,22 +88,29 @@ def brainprep_defacing( Bunch A dictionary-like object containing: - - deface_t1_file : File - path to the defaced image. + - deface_anatomical_file : File - path to the defaced image. - mask_file : File - path to the defacing mask. - mosaic_file : File - path to defacing snapshots. - - summary_file : File - a TSV file containing voxel counts and + - maskdiff_file : File - a TSV file containing voxel counts and physical volumes (in mm³) for the brain/defacing masks and their intersection. + - correlations_file : File - a TSV file containing mean correlation + of aligned input image to the reference image. + - transform_file : File - path to the 12 dof (T1w) or 6 dof (T2w and + FLAIR coregistration) affine transformation. Raises ------ ValueError - If the T1w file do not follow BIDS convention. + If the input anatomical file do not follow BIDS convention. + If the input modality is not supported. + If a T1w image in the same session has not already been faced for + T2w or FLAIR processings. Notes ----- - This workflow assumes the input image is a valid T1-weighted anatomical - scan. + This workflow assumes a T1w image in the same session has already + been defaced for T2w or FLAIR processings. References ---------- @@ -108,13 +120,11 @@ def brainprep_defacing( Examples -------- >>> from brainprep.config import Config - >>> from brainprep.reporting import RSTReport >>> from brainprep.workflow import brainprep_defacing >>> >>> with Config(dryrun=True, verbose=False): - ... report = RSTReport() ... outputs = brainprep_defacing( - ... t1_file=( + ... anatomical_file=( ... "/tmp/dataset/rawdata/sub-01/ses-01/anat/" ... "sub-01_ses-01_run-01_T1w.nii.gz" ... ), @@ -122,49 +132,151 @@ def brainprep_defacing( ... ) >>> outputs Bunch( - deface_t1_file: PosixPath('...') + deface_anatomical_file: PosixPath('...') mask_file: PosixPath('...') mosaic_file: PosixPath('...') - summary_file: PosixPath('...') + maskdiff_file: PosixPath('...') + correlation_file: PosixPath('...') + transform_file: PosixPath('...') ) """ entities = kwargs.get("entities", {}) if len(entities) == 0: raise ValueError( - f"The T1w file '{t1_file}' is not BIDS-compliant." + f"Input file not BIDS-compliant: {anatomical_file}" ) + modality = entities["mod"] + if modality not in ("T1w", "T2w", "FLAIR"): + raise ValueError( + f"Modality not supported: {entities['mod']}" + ) + + resource_dir = Path(interfaces.__file__).parent.parent / "resources" + template_file = resource_dir / f"MNI152_T1_1mm_brain.nii.gz" workspace_dir = output_dir / f"workspace_{entities['run']}" workspace_dir.mkdir(parents=True, exist_ok=True) print_info(f"setting workspace directory: {workspace_dir}") - reoriented_t1_file = interfaces.reorient( - t1_file, - workspace_dir, + reoriented_anatomical_file = interfaces.reorient( + anatomical_file, + workspace_dir / "01-reorient", entities, ) - brainmask_file = interfaces.brainmask( - reoriented_t1_file, - workspace_dir, + _, brainmask_file = interfaces.brainmask( + reoriented_anatomical_file, + workspace_dir / "02-brainmask", entities, ) - deface_t1_file, mask_file = interfaces.deface( - reoriented_t1_file, + if modality == "T1w": + deface_anatomical_file, mask_file, transform_file = interfaces.deface( + reoriented_anatomical_file, + workspace_dir / "03-deface", + entities, + ) + aligned_anatomical_file = interfaces.applyaffine( + reoriented_anatomical_file, + template_file, + transform_file, + workspace_dir / "03-deface", + entities, + interpolation="spline", + ) + else: + t1_file = list( + output_dir.glob( + f"sub-{entities['sub']}_ses-{entities['ses']}_run-*_T1w.nii.gz" + ) + ) + mask_t1_file = list( + output_dir.glob( + f"sub-{entities['sub']}_ses-{entities['ses']}_run-*_mod-T1w_" + "defacemask.nii.gz" + ) + ) + if len(t1_file) != 1 or len(mask_t1_file) != 1: + raise ValueError( + f"No T1w defaced image found: {t1_file}, {mask_t1_file}" + ) + t1_file, mask_t1_file = t1_file[0], mask_t1_file[0] + print_info(f"using T1w: {t1_file}") + print_info(f"using defacing mask: {mask_t1_file}") + aligned_anatomical_file = reoriented_anatomical_file + template_file, transform_file = interfaces.affine( + t1_file, + reoriented_anatomical_file, + workspace_dir / "03-deface", + entities, + rigid=True, + quick=True, + ) + mask_file = interfaces.applyaffine( + mask_t1_file, + reoriented_anatomical_file, + transform_file, + workspace_dir / "03-deface", + entities, + interpolation="nearestneighbour", + ) + deface_anatomical_file = interfaces.applymask( + reoriented_anatomical_file, + mask_file, + workspace_dir / "03-deface", + entities, + ) + maskdiff_file = interfaces.maskdiff( + brainmask_file, + mask_file, output_dir, entities, + inv_mask2=True, + ) + correlation_file = interfaces.mean_correlation( + aligned_anatomical_file, + template_file, + output_dir, + correlation_threshold=None, + suffix=f"_{modality}", ) mosaic_file = interfaces.plot_defacing_mosaic( mask_file, - t1_file, + anatomical_file, output_dir, entities, ) - summary_file = interfaces.maskdiff( - brainmask_file, - mask_file, + + basename = "sub-{sub}_ses-{ses}_run-{run}".format(**entities) + out_deface_anatomical_file = output_dir / f"{basename}_{modality}.nii.gz" + out_mask_file = output_dir / f"{basename}_mod-{modality}_defacemask.nii.gz" + out_summary_file = output_dir / f"{basename}_mod-{modality}_maskinter.tsv" + out_transform_file = output_dir / f"{basename}_mod-{modality}_affine.txt" + out_correlation_file = ( + correlation_file.parent / f"{basename}_mod-{modality}_corr.tsv" + ) + interfaces.copyfiles( + [ + deface_anatomical_file, + mask_file, + transform_file, + correlation_file, + ], + [ + out_deface_anatomical_file, + out_mask_file, + out_transform_file, + out_correlation_file, + ], output_dir, - entities, - inv_mask2=True, + ) + interfaces.copyfiles( + [ + correlation_file, + ], + [ + out_correlation_file, + ], + output_dir, + move_files=True, ) if not keep_intermediate: @@ -172,10 +284,12 @@ def brainprep_defacing( shutil.rmtree(workspace_dir) return Bunch( - deface_t1_file=deface_t1_file, - mask_file=mask_file, + deface_anatomical_file=out_deface_anatomical_file, + mask_file=out_mask_file, mosaic_file=mosaic_file, - summary_file=summary_file, + maskdiff_file=maskdiff_file, + correlation_file=out_correlation_file, + transform_file=out_transform_file, ) @@ -194,26 +308,38 @@ def brainprep_defacing( ] ) def brainprep_group_defacing( + modality: str, output_dir: Directory, overlap_threshold: float = 0.05, + correlation_threshold: float = 0.5, keep_intermediate: bool = False) -> Bunch: """ - Group level defacing pre-processing. + Group-level defacing pre-processing. - Applies the following quality control procedure: + This function applies a quality control procedure to defaced images at + the group level. It includes the following steps: 1) Generate a TSV table containing the intersection between the brain and - defacing masks. - 2) Apply threshold-based quality checks on the selected quality metrics. - 3) Generate a histogram showing the distribution of these quality metrics. + defacing masks. The optimal scenario is when there is no intersection. + 2) Generate a TSV file containing the mean correlation of each image to + the reference image (MNI for T1w or T1w for T2w and FLAIR). The optimal + scenario is when the correlation is maximized. + 3) Apply threshold-based quality checks on the selected quality metrics. + 4) Generate a histogram showing the distribution of these quality metrics. Parameters ---------- + modality : str + Modality: T1w, T2w or FLAIR. output_dir : Directory - Directory where the quality assurance related outputs will be saved + Directory where the defacing related outputs will be saved (i.e., the root of your dataset). overlap_threshold : float - Quality control threshold on the overalp score. Default 0.05. + Quality control threshold on the overalp score. + Default 0.05. + correlation_threshold : float + Quality control threshold on the correlation score. + Default 0.5. keep_intermediate : bool If True, retains intermediate results (no effect on this workflow). Default False. @@ -223,57 +349,99 @@ def brainprep_group_defacing( Bunch A dictionary-like object containing: + - correlations_file : File - a TSV file containing mean correlation + of each input image to the reference image. + - correlation_histogram_file : File - a PNG file containing the + histogram of the computed mean correlations. - overalp_file : File - a TSV file containing brain/defacing masks - intersection quality check (QC) data. + intersections. - overalp_histogram_file : File - PNG file containing the histogram of the computed overlaps. + Raises + ------ + ValueError + If the input modality is not supported. + Notes ----- This workflow assumes the subject-level analyses have already been performed. - A ``qc`` column is added to the TSV QC output table. It contains a binary flag indicating whether the produced results should be kept: ``qc = 1`` if the result passes the thresholds, otherwise ``qc = 0``. - The associated PNG histograms help verify that the chosen thresholds are neither too restrictive nor too permissive. Examples -------- >>> from brainprep.config import Config - >>> from brainprep.reporting import RSTReport >>> from brainprep.workflow import brainprep_group_defacing >>> >>> with Config(dryrun=True, verbose=False): - ... report = RSTReport() ... outputs = brainprep_group_defacing( + ... modality="T1w", ... output_dir="/tmp/dataset/derivatives", ... ) >>> outputs Bunch( + correlations_file: PosixPath('...') + correlation_histogram_file: PosixPath('...') overlap_file: PosixPath('...') overalp_histogram_file: PosixPath('...') ) """ + if modality not in ("T1w", "T2w", "FLAIR"): + raise ValueError( + f"Modality not supported: {modality}" + ) + + correlations_file = interfaces.mean_correlation( + ( + output_dir / + "subjects" / + "sub-*" / + "ses-*" / + "quality_check" / + f"*mod-{modality}_corr.tsv" + ), + None, + output_dir, + correlation_threshold, + suffix=f"_{modality}", + ) + correlation_histogram_file = interfaces.plot_histogram( + correlations_file, + "mean_correlation", + output_dir, + bar_coords=[correlation_threshold], + suffix=f"_{modality}", + ) + overlap_file = interfaces.mask_overlap( ( output_dir / - "subjects" / "sub-*" / "ses-*" / - "*mod-T1w_defacemask.tsv" + "subjects" / + "sub-*" / + "ses-*" / + "quality_check" / + f"*mod-{modality}_maskdiff.tsv" ), output_dir, overlap_threshold, + suffix=f"_{modality}", ) overalp_histogram_file = interfaces.plot_histogram( overlap_file, "overlap", output_dir, bar_coords=[overlap_threshold], + suffix=f"_{modality}", ) return Bunch( + correlations_file=correlations_file, + correlation_histogram_file=correlation_histogram_file, overlap_file=overlap_file, overalp_histogram_file=overalp_histogram_file, ) diff --git a/brainprep/workflow/quasiraw.py b/brainprep/workflow/quasiraw.py index 1c9ec165..d741926a 100644 --- a/brainprep/workflow/quasiraw.py +++ b/brainprep/workflow/quasiraw.py @@ -52,35 +52,48 @@ def brainprep_quasiraw( anatomical_file: File, output_dir: Directory, + rigid: bool = False, + quick: bool = False, keep_intermediate: bool = False, **kwargs: dict) -> Bunch: """ Quasi-RAW pre-processing. Applies the Quasi-RAW pre-processing described in - :footcite:p:`dufumier2022openbhb`. This includes: + :footcite:p:`dufumier2022openbhb` to T1-weighted, T2-weighted and FLAIR + MRI images. This includes: 1) Reorient the anatomical image to standard MNI152 template space. 2) Compute a brain mask using a skull-stripping tool. - 3) Apply the brain mask to the anatomical image. + 3) Perform N4 bias field correction. 4) Resample the anatomical image to 1mm isotropic voxel size. - 5) Resample the brain mask image to 1mm isotropic voxel size. - 6) Perform N4 bias field correction. - 7) Linearly (9 dof) register the image to the MNI152 1mm template space. - 8) Apply the registration to the antomical image. - 9) Apply the registration to the brain mask image. - 10) Apply the brain mask to the registered anatomical image. + 5) Linearly register the image to the MNI152 1mm template space (6 or 9 + DOF). + 6) Apply the registration to the bias field corrected antomical image. + 7) Apply the registration to the brain mask image. Parameters ---------- - anatomical_file: File - Path to the input image file. - output_dir: Directory + anatomical_file : File + Path to the input image file: T1w, T2w or FLAIR. + output_dir : Directory Directory where the outputs will be saved (i.e., the root of your dataset). + rigid : bool + Estimate a 6 DOF transformation that maintains the original size and + shape of the brain. By default a 9 DOF transformation allows for + additional scaling in the x, y, and z directions, adjusting the size + of the brain during the alignment process. + Default False. + quick : bool + Speed up processing by applying optimizations that trade accuracy + for computational efficiency. This is particularly useful for + large-scale batch processing where speed is prioritized. + Default False. keep_intermediate : bool If True, retains intermediate results (i.e., the workspace); useful - for debugging. Default False. + for debugging. + Default False. **kwargs : dict entities: dict Dictionary of parsed BIDS entities. @@ -90,9 +103,9 @@ def brainprep_quasiraw( Bunch A dictionary-like object containing: - - aligned_anatomical_file : File - path to the aligned anatomical + - aligned_anatomical_file : File - path to the aligned 1 mm anatomical image - a Nifti file with the suffix "_T1w". - - aligned_mask_file : File - path to the aligned mask image - a + - aligned_mask_file : File - path to the aligned 1 mm mask image - a Nifti file with the suffix "_mod-T1w_brainmask". - transform_file : File - path to the 9 dof affine transformation - a text file with the suffix "_mod-T1w_affine". @@ -100,11 +113,20 @@ def brainprep_quasiraw( Raises ------ ValueError - If the input anatomical file is not BIDS-compliant. + If the input anatomical file is not BIDS-compliant or if the input + modality is not supported. Notes ----- - This workflow assumes the anatomical image is organized in BIDS. + This workflow assumes the anatomical image is organized in BIDS and applies + the following optimizations in `quick` mode: + + - **Use a coarser resolution**: Increase the shrink factor from `1` to `4` + to downsample the image before estimating the bias field, employ the + MNI152 2mm template as the reference image and scale data to a 2mm + space. + - **Use a Coarser Search Space**: Restricted rotation search range to + +/-30° on all three axes for the registration. References ---------- @@ -140,7 +162,14 @@ def brainprep_quasiraw( ) resource_dir = Path(interfaces.__file__).parent.parent / "resources" - template_file = resource_dir / "MNI152_T1_1mm_brain.nii.gz" + modality = entities["mod"] + if modality not in ("T1w", "T2w", "FLAIR"): + raise ValueError( + f"Modality not supported: {entities['mod']}" + ) + modality = "T2" if modality == "FLAIR" else modality[:-1] + template_file = resource_dir / f"MNI152_{modality}_1mm_brain.nii.gz" + lowres_template_file = resource_dir / f"MNI152_{modality}_2mm_brain.nii.gz" print_info(f"setting template file: {template_file}") workspace_dir = output_dir / f"workspace_{entities['run']}" workspace_dir.mkdir(parents=True, exist_ok=True) @@ -151,46 +180,44 @@ def brainprep_quasiraw( workspace_dir / "01-reorient", entities, ) - mask_file = interfaces.brainmask( + _, mask_file = interfaces.brainmask( reoriented_anatomical_file, workspace_dir / "02-brainmask", entities, ) - masked_anatomical_file = interfaces.applymask( + bc_anatomical_file, _ = interfaces.biasfield( reoriented_anatomical_file, mask_file, - workspace_dir / "03-applymask", + workspace_dir / "03-biasfield", entities, + quick=quick, ) - scaled_anatomical_file, _ = interfaces.scale( - masked_anatomical_file, - 1, - workspace_dir / "04-scale", - entities, - ) - scaled_mask_file, _ = interfaces.scale( + bc_brain_file = interfaces.applymask( + bc_anatomical_file, mask_file, - 1, - workspace_dir / "05-scale", + workspace_dir / "03-biasfield", entities, ) - bc_anatomical_file, _ = interfaces.biasfield( - scaled_anatomical_file, - scaled_mask_file, - workspace_dir / "06-biasfield", + scaled_anatomical_file, _ = interfaces.scale( + bc_brain_file, + 2 if quick else 1, + workspace_dir / "04-scale", entities, + interpolation="trilinear" if quick else "spline", ) _, affine_transform_file = interfaces.affine( - bc_anatomical_file, - template_file, - workspace_dir / "07-affine", + scaled_anatomical_file, + lowres_template_file if quick else template_file, + workspace_dir / "05-affine", entities, + rigid=rigid, + quick=quick, ) aligned_anatomical_file = interfaces.applyaffine( bc_anatomical_file, template_file, affine_transform_file, - workspace_dir / "08-applyaffine", + workspace_dir / "06-applyaffine", entities, interpolation="spline", ) @@ -198,16 +225,10 @@ def brainprep_quasiraw( mask_file, template_file, affine_transform_file, - workspace_dir / "09-applyaffine", + workspace_dir / "07-applyaffine", entities, interpolation="nearestneighbour", ) - aligned_anatomical_file = interfaces.applymask( - aligned_anatomical_file, - aligned_mask_file, - workspace_dir / "10-applymask", - entities, - ) mod = entities["mod"] basename = "sub-{sub}_ses-{ses}_run-{run}".format(**entities) @@ -215,8 +236,16 @@ def brainprep_quasiraw( output_mask_file = output_dir / f"{basename}_mod-{mod}_brainmask.nii.gz" output_transform_file = output_dir / f"{basename}_mod-{mod}_affine.txt" interfaces.copyfiles( - [aligned_anatomical_file, aligned_mask_file, affine_transform_file], - [output_anatomical_file, output_mask_file, output_transform_file], + [ + aligned_anatomical_file, + aligned_mask_file, + affine_transform_file, + ], + [ + output_anatomical_file, + output_mask_file, + output_transform_file, + ], output_dir, ) @@ -246,6 +275,7 @@ def brainprep_quasiraw( ] ) def brainprep_group_quasiraw( + modality: str, output_dir: Directory, correlation_threshold: float = 0.5, keep_intermediate: bool = False) -> Bunch: @@ -256,22 +286,26 @@ def brainprep_group_quasiraw( This includes: 1) Generate a TSV file containing the mean correlation of each image to - the template. + the template. The optimal scenario is when the correlation is maximized. 2) Apply threshold-based quality checks on the selected quality metrics. 3) Generate a histogram showing the distribution of these quality metrics. - 4) Computing a PCA embedding of the images. - 5) Generating a scatter plot of the first two PCA components with BIDS + 4) Compute a PCA embedding of the images. + 5) Generate a scatter plot of the first two PCA components with BIDS annotations for visual inspection. Parameters ---------- + modality : str + Modality: T1w, T2w or FLAIR. output_dir : Directory Working directory containing all the subjects. correlation_threshold : float - Quality control threshold on the correlation score. Default 0.5. + Quality control threshold on the correlation score. + Default 0.5. keep_intermediate : bool If True, retains intermediate results (i.e., the workspace); useful - for debugging. Default False. + for debugging. + Default False. Returns ------- @@ -279,37 +313,40 @@ def brainprep_group_quasiraw( A dictionary-like object containing: - correlations_file : File - a TSV file containing mean correlation - of each input image to the atlas image quality check (QC) data. - - correlation_histogram_file : File - PNG file containing the + of each input image to the atlas image. + - correlation_histogram_file : File - a PNG file containing the histogram of the computed mean correlations. - pca_file : File - a TSV file containing PCA two first components as two columns named ``pc1`` and ``pc2``, as well as BIDS ``participant_id``, ``session``, and ``run``. - - pca_image_file : File - PNG file containing the two first PCA + - pca_image_file : File - a PNG file containing the two first PCA components with ``participant_id``, ``session``, and ``run`` annotations. + Raises + ------ + ValueError + If the input modality is not supported. + Notes ----- This workflow assumes the subject-level analyses have already been performed. - - A ``qc`` column is added to the TSV QC output table. It contains a - binary flag indicating whether the produced results should be kept: - ``qc = 1`` if the result passes the thresholds, otherwise ``qc = 0``. - + A ``qc`` column is added to the ``correlations_file`` output table. + It contains a binary flag indicating whether the produced results should + be kept: ``qc = 1`` if the result passes the thresholds, otherwise + ``qc = 0``. The associated PNG histograms help verify that the chosen thresholds are neither too restrictive nor too permissive. Examples -------- >>> from brainprep.config import Config - >>> from brainprep.reporting import RSTReport >>> from brainprep.workflow import brainprep_group_quasiraw >>> >>> with Config(dryrun=True, verbose=False): - ... report = RSTReport() ... outputs = brainprep_group_quasiraw( + ... modality="T1w", ... output_dir="/tmp/dataset/derivatives", ... ) >>> outputs @@ -321,30 +358,39 @@ def brainprep_group_quasiraw( ) """ resource_dir = Path(interfaces.__file__).parent.parent / "resources" - template_file = resource_dir / "MNI152_T1_1mm_brain.nii.gz" + if modality not in ("T1w", "T2w", "FLAIR"): + raise ValueError( + f"Modality not supported: {modality}" + ) + modality_ = "T2" if modality == "FLAIR" else modality[:-1] + template_file = resource_dir / f"MNI152_{modality_}_1mm_brain.nii.gz" print_info(f"setting template file: {template_file}") correlations_file = interfaces.mean_correlation( - output_dir / "subjects" / "sub-*" / "ses-*" / "*_T1w.nii.gz", + output_dir / "subjects" / "sub-*" / "ses-*" / f"*_{modality}.nii.gz", template_file, output_dir, correlation_threshold, + suffix=f"_{modality}", ) correlation_histogram_file = interfaces.plot_histogram( correlations_file, "mean_correlation", output_dir, bar_coords=[correlation_threshold], + suffix=f"_{modality}", ) pca_file = interfaces.incremental_pca( - output_dir / "subjects" / "sub-*" / "ses-*" / "*_T1w.nii.gz", + output_dir / "subjects" / "sub-*" / "ses-*" / f"*_{modality}.nii.gz", output_dir, batch_size=50, + suffix=f"_{modality}", ) pca_image_file = interfaces.plot_pca( pca_file, output_dir, + suffix=f"_{modality}", ) return Bunch( diff --git a/brainprep/workflow/sulcirec.py b/brainprep/workflow/sulcirec.py index 74e6b373..20ed3e48 100755 --- a/brainprep/workflow/sulcirec.py +++ b/brainprep/workflow/sulcirec.py @@ -110,8 +110,8 @@ def brainprep_sulcirec( ... ) >>> outputs Bunch( - sulci_graphs_files=[PosixPath('...'), PosixPath('...')], - qc_file=PosixPath('...') + sulci_graphs_files: [PosixPath('...'), PosixPath('...')] + qc_file: PosixPath('...') ) """ entities = kwargs.get("entities", {}) @@ -207,7 +207,7 @@ def brainprep_group_sulcirec( ... ) >>> outputs Bunch( - morphometry_files: [PosixPath('...'), PosixPath('...')], + morphometry_files: [PosixPath('...'), PosixPath('...')] group_stats_file: PosixPath('...') ) """ diff --git a/doc/user_guide/defacing.rst b/doc/user_guide/defacing.rst index a27f1c3d..2641e440 100644 --- a/doc/user_guide/defacing.rst +++ b/doc/user_guide/defacing.rst @@ -30,7 +30,7 @@ Description **Processing Steps** - **Defacing T1w image** - We use the UK-Biobank defacing method that is provided as part of the + We use the defacing method that is provided as part of the standard FSL distribution under the command ``fsl_deface`` :footcite:p:`almagro2018deface`. Similar to other established defacing tools, such as ``mri_deface`` :footcite:p:`bischoff2007deface` @@ -40,30 +40,44 @@ Description A key distinction of ``fsl_deface`` compared with ``mri_deface`` and ``pydeface`` is that it additionally removes the ears, providing a more comprehensive anonymization of head anatomy. - This workflow is applied to T1-weighted (T1w) structural images and can be - propagated to other modalities through rigid alignement. + This workflow is applied to T1-weighted structural images. + +- **Defacing T2w and FLAIR images** + The defacing mask obtained from the T1w image is propagated to other + modalities through rigid alignement. **Quality Control** -- **Overlap score** + +- **Correlation score** + For each aligned image, we compute its correlation with the MNI for T1w or + to T1w for T2w and FLAIR reference image. Images are then sorted in + ascending order of this score, allowing potential outliers to be easily + identified. + +- **Overlap score** For each image, we compute the overlap ratio between the brain mask extracted using FreeSurfer's deep‑learning–based ``mri_synthstrip`` method :footcite:p:`hoopes2022brainmask` and the corresponding defacing mask. Images are then sorted in ascending order of this score, allowing potential outliers to be easily identified. -- **Manual inspection** +- **Manual inspection** Following the overlap-based ranking, a manual quality control step is performed using the generated ``defacemosaic`` figure. This figure allows visual inspection of the defaced image and verification that facial/ear structures have been successfully removed. The most obvious outliers are thus removed. -- **Thresholding** +- **Thresholding** + The correlation score is thresholded at 0.5, meaning that if an image is not + roughly registered to the template, the preprocessing is considered invalid. The overlap score is thresholded at 5%, meaning that if the defacing mask removes a substantial portion of the brain, the preprocessing is considered - invalid. Images with an overlap greater than 5% are flagged as low‑quality. + invalid. + Images with a correlation lower than 0.5 or an overlap greater than 5% are + flagged as low‑quality. Outputs ------- @@ -77,45 +91,68 @@ The structure is organized following the :ref:`brainprep ontology `. defacing/ ├── dataset_description.json ├── figures - │   └── histogram_overlap.png + │   ├── histogram_mean_correlation_.png + │   └── histogram_overlap_.rst + │ ├── report_.rst + │ └── commands_.rst ├── quality_check - │   └── mask_overlap.tsv + │   ├── mask_overlap_.tsv + │ └── mean_correlations_.tsv └── subjects └── sub-01 └── ses-01 ├── figures - │   └── sub-01_ses-01_run-01_mod-T1w_defacemosaic.png + │   └── sub-01_ses-01_run-01_mod-_defacemosaic.png ├── log - │   └── report_.rst - ├── sub-01_ses-01_run-01_mod-T1w_defacemask.nii.gz - ├── sub-01_ses-01_run-01_mod-T1w_defacemask.tsv - └── sub-01_ses-01_run-01_mod-T1w_deface.nii.gz + │ ├── report_.rst + │ └── commands_.rst + ├── quality_check + │ ├── sub-01_ses-01_run-01_mod-_corr.tsv + │ └── sub-01_ses-01_run-01_mod-_maskdiff.tsv + ├── sub-01_ses-01_run-01_mod-_affine.mat + ├── sub-01_ses-01_run-01_mod-_defacemask.nii.gz + ├── sub-01_ses-01_run-01_mod-_maskdiff.tsv + └── sub-01_ses-01_run-01_mod-.nii.gz **Description of contents**: - ``dataset_description.json`` Metadata describing the process, including versioning and processing information. -- ``figures/histogram_overlap.png`` +- ``figures/histogram_mean_correlation_.png`` + Image correlation-to-reference image distribution and applied threshold. +- ``figures/histogram_overlap_.png`` Image overlap distribution and applied threshold. -- ``logs/report_.rst`` +- ``logs/report_.rst`` Contains group-level workflow steps and parameters. -- ``quality_check/mask_overlap.tsv`` +- ``log/commands_.rst`` + Contains group-level executed commands. +- ``quality_check/mask_overlap_.tsv`` Table containing the overlap score for each subject/session/run. The table includes a binary ``qc`` column indicating the quality control result. -- ``subjects/sub-/ses-/figures/sub-01_ses-01_run-01_mod-T1w_defacemosaic.png`` +- ``quality_check/mean_correlations_.tsv`` + Table containing the correlation score for each subject/session/run. The + table includes a binary ``qc`` column indicating the quality control result. +- ``subjects/sub-/ses-/figures/sub-01_ses-01_run-01_mod-_defacemosaic.png`` A visual mosaic showing defacing masks on some slices for quick quality check. -- ``subjects/sub-/ses-/logs/report_.rst`` +- ``subjects/sub-/ses-/logs/report_.rst`` Contains subject-level workflow steps and parameters. -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-T1w_defacemask.nii.gz`` - The binary mask identifying voxels removed during defacing (face and ears). -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-T1w_defacemask.tsv`` +- ``subjects/sub-/ses-/log/commands_.rst`` + Contains subject-level executed commands. +- ``subjects/sub-/ses-/quality_check/sub-01_ses-01_run-01_mod-_corr.tsv`` + Table containing the correlation score. +- ``subjects/sub-/ses-/quality_check/sub-01_ses-01_run-01_mod-_maskdiff.tsv`` A table containing voxel counts and physical volumes (in mm³) for the brain/defacing masks and their intersection. -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-T1w_deface.nii.gz`` - The final defaced T1w image with facial/ear structures removed. +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-_affine.txt`` + Affine transformation parameters (12 DOF) used to align the T1w image to the + MNI template or rigid transformation parameters (6 DOF) used to coregister + the T2w or FLAIR image to the T1w image. +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-_defacemask.nii.gz`` + The binary mask identifying voxels removed during defacing (face and ears). +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-.nii.gz`` + The final defaced T1w, T2w or FLAIR image with facial/ear structures removed. Featured examples ----------------- diff --git a/doc/user_guide/quasiraw.rst b/doc/user_guide/quasiraw.rst index acb9e68c..2383c4cc 100644 --- a/doc/user_guide/quasiraw.rst +++ b/doc/user_guide/quasiraw.rst @@ -78,48 +78,57 @@ The structure is organized following the :ref:`brainprep ontology `. quasiraw/ ├── dataset_description.json ├── figures - │   ├── histogram_mean_correlation.png - │   └── pca.png + │   ├── histogram_mean_correlation_.png + │   └── pca_.png ├── log - │ └── report_.rst + │ ├── report_.rst + │ └── commands_.rst ├── quality_check - │ ├── mean_correlations.tsv - │ └── pca.tsv + │ ├── mean_correlations_.tsv + │ └── pca_.tsv └── subjects └── sub-01 └── ses-01 ├── log - │ └── report_.rst - ├── sub-01_ses-01_run-01_mod-T1w_affine.txt - ├── sub-01_ses-01_run-01_mod-T1w_brainmask.nii.gz - └── sub-01_ses-01_run-01_T1w.nii.gz + │ ├── report_.rst + │ └── commands_.rst + ├── quality_check + │ └── sub-01_ses-01_run-01_mod-_corr.tsv + ├── sub-01_ses-01_run-01_mod-_affine.txt + ├── sub-01_ses-01_run-01_mod-_brainmask.nii.gz + └── sub-01_ses-01_run-01_.nii.gz **Description of contents**: - ``dataset_description.json`` Metadata describing the process, including versioning and processing information. -- ``figures/histogram_mean_correlation.png`` +- ``figures/histogram_mean_correlation_.png`` Image correlation-to-template distribution and applied threshold. -- ``figures/pca.png`` +- ``figures/pca_.png`` Display of the first two PCA components of the generated images. - ``log/report_.rst`` Contains group-level workflow steps and parameters. -- ``quality_check/mask_overlap.tsv`` +- ``log/commands_.rst`` + Contains group-level executed commands. +- ``quality_check/mean_correlations_.tsv`` Table containing the correlation score for each subject/session/run. The table includes a binary ``qc`` column indicating the quality control result. -- ``quality_check/pca.tsv`` +- ``quality_check/pca_.tsv`` Table containing information on the first two PCA components. -- ``subjects/sub-/ses-/log/report_.rst`` +- ``subjects/sub-/ses-/log/report_.rst`` Contains subject-level workflow steps and parameters. -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-T1w_affine.txt`` - Affine transformation parameters (9 DOF) used to align the T1w image to - the MNI template. -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-T1w_brainmask.nii.gz`` +- ``subjects/sub-/ses-/log/commands_.rst`` + Contains subject-level executed commands. +- ``subjects/sub-/ses-/quality_check/sub-01_ses-01_run-01_mod-_corr.tsv`` + Table containing the correlation score. +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-_affine.txt`` + Affine transformation parameters (9 DOF) used to align the T1w, T2w or FLAIR + image to the MNI 1 mm template. +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_mod-_brainmask.nii.gz`` Brain mask generated during skull stripping (e.g., via SynthStrip). -- ``subjects/sub-/ses-/sub-01_ses-01_run-01_T1w.nii.gz`` - The minimally preprocessed T1w image, including skull stripping, bias - correction, and affine alignment. +- ``subjects/sub-/ses-/sub-01_ses-01_run-01_.nii.gz`` + The minimally preprocessed T1w, T2w or FLAIR image in the MNI 1mm space. Featured examples ----------------- diff --git a/examples/workflows/plot_defacing.py b/examples/workflows/plot_defacing.py index e9fbd654..d4eb8c75 100644 --- a/examples/workflows/plot_defacing.py +++ b/examples/workflows/plot_defacing.py @@ -10,20 +10,30 @@ Data ---- -Let's first get some anatomical data. +Let's first get some anatomical data: T1w, T2w and FLAIR.. """ from pathlib import Path +from brainprep.utils import Bunch from brainprep.datasets import OpenMSDataset datadir = Path("/tmp/brainprep-data") datadir.mkdir(parents=True, exist_ok=True) dataset = OpenMSDataset(datadir) -data = dataset.fetch( - subject="01", - modality="T1w", - dtype="cross_sectional", -) +data = Bunch() +for modality in ("T1w", "T2w", "FLAIR"): + data[modality] = Bunch( + sub01=dataset.fetch( + subject="01", + modality=modality, + dtype="cross_sectional", + ), + sub02=dataset.fetch( + subject="02", + modality=modality, + dtype="cross_sectional", + ), + ) print(data) @@ -42,23 +52,25 @@ brainprep_group_defacing, ) from brainprep.config import Config -from brainprep.reporting import RSTReport outdir = Path("/tmp/brainprep-defacing") if outdir.is_dir(): shutil.rmtree(outdir) outdir.mkdir(parents=True, exist_ok=True) with Config(dryrun=True, verbose=True): - report = RSTReport() - brainprep_defacing( - t1_file=data.anat, - output_dir=outdir, - keep_intermediate=True, - ) - print(report) - brainprep_group_defacing( - output_dir=outdir, - ) + for modality, modality_data in data.items(): + for subject_data in modality_data.values(): + outputs = brainprep_defacing( + anatomical_file=subject_data.anat, + output_dir=outdir, + keep_intermediate=True, + ) + outputs.deface_anatomical_file.touch(exist_ok=True) + outputs.mask_file.touch(exist_ok=True) + outputs = brainprep_group_defacing( + modality=modality, + output_dir=outdir, + ) # %% @@ -76,18 +88,33 @@ [ [ "brainprep", "subject-level-defacing", - "--t1_file", str(data.anat), + "--anatomical_file", str(subject_data.anat), + "--output-dir", str(outdir), + "--keep-intermediate", + ] + for subject_data in data["T1w"].values() + ] +) +commands.append( + [ + [ + "brainprep", "subject-level-defacing", + "--anatomical_file", str(subject_data.anat), "--output-dir", str(outdir), "--keep-intermediate", ] + for mod in ("T2w", "FLAIR") + for subject_data in data[mod].values() ] ) commands.append( [ [ "brainprep", "group-level-defacing", + "--modality", modality, "--output-dir", str(outdir), ] + for modality in data.keys() ] ) pprint(commands) diff --git a/examples/workflows/plot_quasiraw.py b/examples/workflows/plot_quasiraw.py index e1a50b75..2f72c975 100644 --- a/examples/workflows/plot_quasiraw.py +++ b/examples/workflows/plot_quasiraw.py @@ -4,13 +4,13 @@ Simple example. -Example on how to run the brain parcellation pre-processing using BrainPrep. +Example on how to run the quasiraw pre-processing using BrainPrep. See :ref:`user guide ` for details. Data ---- -Let's first get some anatomical data. +Let's first get some anatomical data: T1w, T2w and FLAIR. """ from pathlib import Path @@ -20,18 +20,20 @@ datadir = Path("/tmp/brainprep-data") datadir.mkdir(parents=True, exist_ok=True) dataset = OpenMSDataset(datadir) -data = Bunch( - sub01=dataset.fetch( - subject="01", - modality="T1w", - dtype="cross_sectional", - ), - sub02=dataset.fetch( - subject="02", - modality="T1w", - dtype="cross_sectional", - ), -) +data = Bunch() +for modality in ("T1w", "T2w", "FLAIR"): + data[modality] = Bunch( + sub01=dataset.fetch( + subject="01", + modality=modality, + dtype="cross_sectional", + ), + sub02=dataset.fetch( + subject="02", + modality=modality, + dtype="cross_sectional", + ), + ) print(data) @@ -50,24 +52,23 @@ brainprep_group_quasiraw, ) from brainprep.config import Config -from brainprep.reporting import RSTReport outdir = Path("/tmp/brainprep-quasiraw") if outdir.is_dir(): shutil.rmtree(outdir) outdir.mkdir(parents=True, exist_ok=True) with Config(dryrun=True, verbose=True): - for subject_data in data.values(): - report = RSTReport() - brainprep_quasiraw( - anatomical_file=subject_data.anat, + for modality, modality_data in data.items(): + for subject_data in modality_data.values(): + outputs = brainprep_quasiraw( + anatomical_file=subject_data.anat, + output_dir=outdir, + keep_intermediate=True, + ) + outputs = brainprep_group_quasiraw( + modality=modality, output_dir=outdir, - keep_intermediate=True, ) - print(report) - outputs = brainprep_group_quasiraw( - output_dir=outdir, - ) # %% @@ -88,15 +89,19 @@ "--anatomical_file", str(subject_data.anat), "--output-dir", str(outdir), "--keep-intermediate", - ] for subject_data in data.values() + ] + for modality_data in data.values() + for subject_data in modality_data.values() ] ) commands.append( [ [ "brainprep", "group-level-quasiraw", + "--modality", modality, "--output-dir", str(outdir), ] + for modality in data.keys() ] ) pprint(commands)