21672251336

Committed 04 Feb 2026 12:54PM UTC coverage: 78.442% (-0.07%) from 78.514%

Build # 21672251336

Build Type

Pull #35

github

Committed by

web-flow

Commit Message

Merge 63641075b into 796543878

Pull Request Pull Request #35: Group level vbm

Coverage Stats

3 of 14 new or added lines in 2 files covered. (21.43%)

2 existing lines in 2 files now uncovered.

1470 of 1874 relevant lines covered (78.44%)

0.78 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

46.28

/brainprep/interfaces/utils.py

##########################################################################
# NSAp - Copyright (C) CEA, 2021 - 2025
# Distributed under the terms of the CeCILL-B license, as published by
# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
# for details.
##########################################################################

"""
Tools.
"""


import glob
import gzip
import shutil
from pathlib import Path

import nibabel
import numpy as np
import pandas as pd
from scipy.stats import pearsonr
from sklearn.decomposition import IncrementalPCA

from ..reporting import log_runtime
from ..typing import (
    Directory,
    File,
)
from ..utils import (
    coerceparams,
    coerce_to_path,
    outputdir,
    parse_bids_keys,
)
from ..wrappers import pywrapper
from .plotting import plot_pca


@coerceparams
@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def mask_diff(
        im1_file: File,
        im2_file: File,
        output_dir: Directory,
        entities: dict,
        threshold: float = 0.6,
        dryrun: bool = False) -> tuple[File]:
    """
    Computes a defacing mask by thresholding the intensity difference between
    two input images.

    Parameters
    ----------
    im1_file : File
        Path to the first image.
    im2_file : File
        Path to the second image.
    output_dir : Directory
        Directory where the defacing mask will be saved.
    entities : dict
        A dictionary of parsed BIDS entities including modality.
    threshold : float
        Threshold for intensity difference to define the mask. Default 0.6.
    dryrun : bool
        If True, skip actual computation and file writing. Default False.

    Returns
    -------
    mask_file : File
        Path to the saved mask image.
    """
    basename = "sub-{sub}_ses-{ses}_run-{run}_mod-T1w_defacemask".format(
        **entities)
    mask_file = output_dir / f"{basename}.nii.gz"

    if not dryrun:
        im1 = nibabel.load(im1_file)
        im2 = nibabel.load(im2_file)
        mask = np.abs(im1.get_fdata() - im2.get_fdata())
        indices = np.where(mask > threshold)
        mask[...] = 0
        mask[indices] = 1
        im_mask = nibabel.Nifti1Image(mask, im1.affine)
        nibabel.save(im_mask, mask_file)

    return (mask_file, )


@coerceparams
@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def copyfiles(
        source_image_files: list[File],
        destination_image_files: list[File],
        output_dir: Directory,
        dryrun: bool = False) -> None:
    """
    Copy input image files.

    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.
    output_dir : Directory
        Directory where the images are copied.
    dryrun : bool
        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,
                                       strict=True):
            shutil.copy(src_path, dest_path)


@coerceparams
@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def movedir(
        source_dir: Directory,
        output_dir: Directory,
        content: bool = False,
        dryrun: bool = False) -> tuple[Directory]:
    """
    Move input directory.

    Parameters
    ----------
    source_dir : Directory
        Path to the directory to be moved.
    output_dir : Directory
        Directory where the folder is moved.
    content : bool
        If True, copy only the content of the source directory. Default False.
    dryrun : bool
        If True, skip actual computation and file writing. Default False.

    Returns
    -------
    target_directory : Directory
        Path to the moved directory.

    Raises
    ------
    ValueError
        If `source_dir` is not a directory.
    """
    target_directory = output_dir
    if not content:
        target_directory /= source_dir.name
    if not dryrun:
        if not source_dir.is_dir():
            raise ValueError(
                f"Source '{source_dir}' is not a directory."
            )
        if not content:
            shutil.copytree(source_dir, output_dir)
        else:
            for item in source_dir.iterdir():
                target = output_dir / item.name
                if item.is_dir():
                    shutil.copytree(item, target, dirs_exist_ok=True)
                else:
                    shutil.copy2(item, target)
            if not any(source_dir.iterdir()):
                source_dir.rmdir()
    return (target_directory, )


@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def ungzfile(
        input_file: File,
        output_file: File,
        output_dir: Directory,
        dryrun: bool = False) -> tuple[File]:
    """
    Ungzip input file.

    Parameters
    ----------
    input_file : File
        Path to the file to ungzip.
    output_file : File
        Path to the ungzip file.
    output_dir : Directory
        Directory where the unzip file is created.
    dryrun : bool
        If True, skip actual computation and file writing. Default False.

    Returns
    -------
    output_file : File
        Path to the ungzip file.

    Raises
    ------
    ValueError
        If the input file is not compressed.
    """
    if input_file.suffix != ".gz":
        raise ValueError(
            f"The input file is not compressed: {input_file}"
        )

    if not dryrun:
        with gzip.open(input_file, "rb") as gzfobj:
            output_file.write_bytes(gzfobj.read())

    return (output_file, )


@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def mean_correlation(
        image_files_regex: str,
        atlas_file: File,
        output_dir: Directory,
        correlation_threshold: float = 0.5,
        dryrun: bool = False) -> tuple[File]:
    """
    Compute the mean Pearson correlation between a reference image and a list
    of other images.

    Parameters
    ----------
    image_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.
    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.
    dryrun : bool
        If True, skip actual computation and file writing. Default False.

    Returns
    -------
    correlations_file : File
        A TSV file containing the Pearson correlation coefficient between the
        atlas image and each image pointed by the input REGEX.

    Raises
    ------
    ValueError
        If the atlas and an image have incompatible shape or geometry.
    """
    correlations_file = output_dir / "mean_correlations.tsv"

    if not dryrun:

        image_files = coerce_to_path(glob.glob(
            str(image_files_regex)),
            expected_type=list[File]
        )
        atlas_im = nibabel.load(atlas_file)
        atlas_arr = atlas_im.get_fdata()

        scores = pd.DataFrame(
            columns=(
                "participant_id",
                "session",
                "run",
                "mean_correlation",
            )
        )
        for path in image_files:
            entities = parse_bids_keys(Path(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}"
                )
            corr, _ = pearsonr(
                atlas_arr.flatten(),
                arr.flatten(),
            )
            scores.loc[len(scores)] = [
                entities["sub"],
                entities["ses"],
                entities["run"],
                corr,
            ]

        scores["qc"] = (
            scores["mean_correlation"] > correlation_threshold
        ).astype(int)
        scores = scores.sort_values(by=["participant_id","session","run"])
        scores.to_csv(
            correlations_file,
            index=False,
            sep="\t",
        )

    return (correlations_file, )


@outputdir
@log_runtime(
    bunched=False)
@pywrapper
def incremental_pca(
        image_files_regex: str,
        output_dir: Directory,
        batch_size: int = 10,
        dryrun: bool = False) -> tuple[File]:
    """
    Perform an Incremental PCA with 2 components on a collection of images
    matched by a regex pattern, processing them in batches.

    The function loads all images matching the provided regex, splits them
    into batches, and incrementally fits a PCA model using scikit-learn's
    ``IncrementalPCA``. Each image is flattened into a 1D vector before
    processing. After fitting, the function transforms all batches to obtain
    the first two principal components for each image. These components are
    saved in a TSV file as two columns named ``pc1`` and ``pc2``. BIDS
    entities (``participant_id``, ``session``, ``run``) are extracted from
    filenames using ``parse_bids_keys`` and included in the output table.

    Parameters
    ----------
    image_files_regex : str
        A REGEX to image files, each representing an image,
        all images must have the same size.
    output_dir : Directory
        Directory where a TSV file containing the values of the first two
        components created by the PCA ill be saved, a Directory containing
        all the graph of all batch.
    batch_size : int
        Number of images to use in each batch. If None, a single batch is used.
        Default is 10.
    dryrun : bool
        If True, skip actual computation and file writing. Default False.

    Returns
    -------
    pca_file : File
        Path to the generated ``pca.tsv`` file containing the PCA results.

    Raises
    ------
    ValueError
        If no image matches the regex pattern.
        If the dataset contains fewer than 2 images, which prevents PCA
        computation.
    """
    pca_file = output_dir / "pca.tsv"

    if not dryrun:

        image_files = glob.glob(str(image_files_regex))
        n_images = len(image_files)
        if n_images == 0:
            raise ValueError(
                f"No image matches the regex pattern: {image_files_regex}"
            )
        if n_images < 2:
            raise ValueError(
                f"The dataset contains fewer than 2 images: {n_images}"
            )
        batches = [
            image_files[idx:idx + batch_size]
            for idx in range(0, len(image_files), batch_size)
        ]

        pca = IncrementalPCA(n_components=2)
        for batch_files in batches:
            data = [
                nibabel.load(_file).get_fdata().flatten()
                for _file in batch_files
            ]
            pca.partial_fit(data)

        df = []
        for batch_files in batches:
            data = [
                nibabel.load(_file).get_fdata().flatten()
                for _file in batch_files
            ]
            components = pca.transform(data)
            info = [
                parse_bids_keys(Path(_file))
                for _file in batch_files
            ]
            partial_df = pd.DataFrame({
                "participant_id": [item["sub"] for item in info],
                "session": [item["ses"] for item in info],
                "run": [item["run"] for item in info],
                "pc1": components[:, 0],
                "pc2": components[:, 1],
                "explained_variance_ratio_pc1": [
                    pca.explained_variance_ratio_[0],
                ] * len(info),
                "explained_variance_ratio_pc2": [
                    pca.explained_variance_ratio_[1],
                ] * len(info),
            })
            df.append(partial_df)
        df = pd.concat(df, ignore_index=True)
        df.to_csv(pca_file, index=False, sep="\t")

    return (pca_file, )

1	##########################################################################
2	# NSAp - Copyright (C) CEA, 2021 - 2025
3	# Distributed under the terms of the CeCILL-B license, as published by
4	# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
5	# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
6	# for details.
7	##########################################################################
8
9	"""
10	Tools.
11	"""
12
13
14	import glob	1✔
15	import gzip	1✔
16	import shutil	1✔
17	from pathlib import Path	1✔
18
19	import nibabel	1✔
20	import numpy as np	1✔
21	import pandas as pd	1✔
22	from scipy.stats import pearsonr	1✔
23	from sklearn.decomposition import IncrementalPCA	1✔
24
25	from ..reporting import log_runtime	1✔
26	from ..typing import (	1✔
27	Directory,
28	File,
29	)
30	from ..utils import (	1✔
31	coerceparams,
32	coerce_to_path,
33	outputdir,
34	parse_bids_keys,
35	)
36	from ..wrappers import pywrapper	1✔
37	from .plotting import plot_pca	1✔
38
39
40	@coerceparams	1✔
41	@outputdir	1✔
42	@log_runtime(	1✔
43	bunched=False)
44	@pywrapper	1✔
45	def mask_diff(	1✔
46	im1_file: File,
47	im2_file: File,
48	output_dir: Directory,
49	entities: dict,
50	threshold: float = 0.6,
51	dryrun: bool = False) -> tuple[File]:
52	"""
53	Computes a defacing mask by thresholding the intensity difference between
54	two input images.
55
56	Parameters
57	----------
58	im1_file : File
59	Path to the first image.
60	im2_file : File
61	Path to the second image.
62	output_dir : Directory
63	Directory where the defacing mask will be saved.
64	entities : dict
65	A dictionary of parsed BIDS entities including modality.
66	threshold : float
67	Threshold for intensity difference to define the mask. Default 0.6.
68	dryrun : bool
69	If True, skip actual computation and file writing. Default False.
70
71	Returns
72	-------
73	mask_file : File
74	Path to the saved mask image.
75	"""
76	basename = "sub-{sub}_ses-{ses}_run-{run}_mod-T1w_defacemask".format(	×
77	**entities)
78	mask_file = output_dir / f"{basename}.nii.gz"	×
79
80	if not dryrun:	×
81	im1 = nibabel.load(im1_file)	×
82	im2 = nibabel.load(im2_file)	×
83	mask = np.abs(im1.get_fdata() - im2.get_fdata())	×
84	indices = np.where(mask > threshold)	×
85	mask[...] = 0	×
86	mask[indices] = 1	×
87	im_mask = nibabel.Nifti1Image(mask, im1.affine)	×
88	nibabel.save(im_mask, mask_file)	×
89
90	return (mask_file, )	×
91
92
93	@coerceparams	1✔
94	@outputdir	1✔
95	@log_runtime(	1✔
96	bunched=False)
97	@pywrapper	1✔
98	def copyfiles(	1✔
99	source_image_files: list[File],
100	destination_image_files: list[File],
101	output_dir: Directory,
102	dryrun: bool = False) -> None:
103	"""
104	Copy input image files.
105
106	Parameters
107	----------
108	source_image_files : list[File]
109	Path to the image to be copied.
110	destination_image_files : list[File]
111	Path to the locations where images will be copied.
112	output_dir : Directory
113	Directory where the images are copied.
114	dryrun : bool
115	If True, skip actual computation and file writing. Default False.
116	"""
117	if not dryrun:	1✔
118	for src_path, dest_path in zip(source_image_files,	×
119	destination_image_files,
120	strict=True):
121	shutil.copy(src_path, dest_path)	×
122
123
124	@coerceparams	1✔
125	@outputdir	1✔
126	@log_runtime(	1✔
127	bunched=False)
128	@pywrapper	1✔
129	def movedir(	1✔
130	source_dir: Directory,
131	output_dir: Directory,
132	content: bool = False,
133	dryrun: bool = False) -> tuple[Directory]:
134	"""
135	Move input directory.
136
137	Parameters
138	----------
139	source_dir : Directory
140	Path to the directory to be moved.
141	output_dir : Directory
142	Directory where the folder is moved.
143	content : bool
144	If True, copy only the content of the source directory. Default False.
145	dryrun : bool
146	If True, skip actual computation and file writing. Default False.
147
148	Returns
149	-------
150	target_directory : Directory
151	Path to the moved directory.
152
153	Raises
154	------
155	ValueError
156	If `source_dir` is not a directory.
157	"""
158	target_directory = output_dir	1✔
159	if not content:	1✔
160	target_directory /= source_dir.name	1✔
161	if not dryrun:	1✔
162	if not source_dir.is_dir():	×
163	raise ValueError(	×
164	f"Source '{source_dir}' is not a directory."
165	)
166	if not content:	×
167	shutil.copytree(source_dir, output_dir)	×
168	else:
169	for item in source_dir.iterdir():	×
170	target = output_dir / item.name	×
171	if item.is_dir():	×
172	shutil.copytree(item, target, dirs_exist_ok=True)	×
173	else:
174	shutil.copy2(item, target)	×
175	if not any(source_dir.iterdir()):	×
176	source_dir.rmdir()	×
177	return (target_directory, )	1✔
178
179
180	@outputdir	1✔
181	@log_runtime(	1✔
182	bunched=False)
183	@pywrapper	1✔
184	def ungzfile(	1✔
185	input_file: File,
186	output_file: File,
187	output_dir: Directory,
188	dryrun: bool = False) -> tuple[File]:
189	"""
190	Ungzip input file.
191
192	Parameters
193	----------
194	input_file : File
195	Path to the file to ungzip.
196	output_file : File
197	Path to the ungzip file.
198	output_dir : Directory
199	Directory where the unzip file is created.
200	dryrun : bool
201	If True, skip actual computation and file writing. Default False.
202
203	Returns
204	-------
205	output_file : File
206	Path to the ungzip file.
207
208	Raises
209	------
210	ValueError
211	If the input file is not compressed.
212	"""
213	if input_file.suffix != ".gz":	1✔
214	raise ValueError(	×
215	f"The input file is not compressed: {input_file}"
216	)
217
218	if not dryrun:	1✔
219	with gzip.open(input_file, "rb") as gzfobj:	×
220	output_file.write_bytes(gzfobj.read())	×
221
222	return (output_file, )	1✔
223
224
225	@outputdir	1✔
226	@log_runtime(	1✔
227	bunched=False)
228	@pywrapper	1✔
229	def mean_correlation(	1✔
230	image_files_regex: str,
231	atlas_file: File,
232	output_dir: Directory,
233	correlation_threshold: float = 0.5,
234	dryrun: bool = False) -> tuple[File]:
235	"""
236	Compute the mean Pearson correlation between a reference image and a list
237	of other images.
238
239	Parameters
240	----------
241	image_files_regex : str
242	A REGEX to image files, each representing an image of the same shape
243	and geometry as `atlas_file`.
244	atlas_file : File
245	An file representing the reference image.
246	output_dir : Directory
247	Directory where a TSV file containing the mean correlation values is
248	created.
249	correlation_threshold : float
250	Quality control threshold on the correlation score. Default 0.5.
251	dryrun : bool
252	If True, skip actual computation and file writing. Default False.
253
254	Returns
255	-------
256	correlations_file : File
257	A TSV file containing the Pearson correlation coefficient between the
258	atlas image and each image pointed by the input REGEX.
259
260	Raises
261	------
262	ValueError
263	If the atlas and an image have incompatible shape or geometry.
264	"""
265	correlations_file = output_dir / "mean_correlations.tsv"	1✔
266
267	if not dryrun:	1✔
268
NEW 269	image_files = coerce_to_path(glob.glob(	×
270	str(image_files_regex)),
271	expected_type=list[File]
272	)
273	atlas_im = nibabel.load(atlas_file)	×
274	atlas_arr = atlas_im.get_fdata()	×
275
276	scores = pd.DataFrame(	×
277	columns=(
278	"participant_id",
279	"session",
280	"run",
281	"mean_correlation",
282	)
283	)
284	for path in image_files:	×
285	entities = parse_bids_keys(Path(path))	×
286	im = nibabel.load(path)	×
287	arr = atlas_im.get_fdata()	×
288	if atlas_arr.shape != arr.shape:	×
289	raise ValueError(	×
290	f"Atlas and image have incompatible shape: {path}"
291	)
292	if not np.allclose(atlas_im.affine, im.affine):	×
293	raise ValueError(	×
294	f"Atlas and image have incompatible orientation: {path}"
295	)
296	corr, _ = pearsonr(	×
297	atlas_arr.flatten(),
298	arr.flatten(),
299	)
300	scores.loc[len(scores)] = [	×
301	entities["sub"],
302	entities["ses"],
303	entities["run"],
304	corr,
305	]
306
307	scores["qc"] = (	×
308	scores["mean_correlation"] > correlation_threshold
309	).astype(int)
NEW 310	scores = scores.sort_values(by=["participant_id","session","run"])	×
UNCOV 311	scores.to_csv(	×
312	correlations_file,
313	index=False,
314	sep="\t",
315	)
316
317	return (correlations_file, )	1✔
318
319
320	@outputdir	1✔
321	@log_runtime(	1✔
322	bunched=False)
323	@pywrapper	1✔
324	def incremental_pca(	1✔
325	image_files_regex: str,
326	output_dir: Directory,
327	batch_size: int = 10,
328	dryrun: bool = False) -> tuple[File]:
329	"""
330	Perform an Incremental PCA with 2 components on a collection of images
331	matched by a regex pattern, processing them in batches.
332
333	The function loads all images matching the provided regex, splits them
334	into batches, and incrementally fits a PCA model using scikit-learn's
335	``IncrementalPCA``. Each image is flattened into a 1D vector before
336	processing. After fitting, the function transforms all batches to obtain
337	the first two principal components for each image. These components are
338	saved in a TSV file as two columns named ``pc1`` and ``pc2``. BIDS
339	entities (``participant_id``, ``session``, ``run``) are extracted from
340	filenames using ``parse_bids_keys`` and included in the output table.
341
342	Parameters
343	----------
344	image_files_regex : str
345	A REGEX to image files, each representing an image,
346	all images must have the same size.
347	output_dir : Directory
348	Directory where a TSV file containing the values of the first two
349	components created by the PCA ill be saved, a Directory containing
350	all the graph of all batch.
351	batch_size : int
352	Number of images to use in each batch. If None, a single batch is used.
353	Default is 10.
354	dryrun : bool
355	If True, skip actual computation and file writing. Default False.
356
357	Returns
358	-------
359	pca_file : File
360	Path to the generated ``pca.tsv`` file containing the PCA results.
361
362	Raises
363	------
364	ValueError
365	If no image matches the regex pattern.
366	If the dataset contains fewer than 2 images, which prevents PCA
367	computation.
368	"""
369	pca_file = output_dir / "pca.tsv"	1✔
370
371	if not dryrun:	1✔
372
373	image_files = glob.glob(str(image_files_regex))	×
374	n_images = len(image_files)	×
375	if n_images == 0:	×
376	raise ValueError(	×
377	f"No image matches the regex pattern: {image_files_regex}"
378	)
379	if n_images < 2:	×
380	raise ValueError(	×
381	f"The dataset contains fewer than 2 images: {n_images}"
382	)
383	batches = [	×
384	image_files[idx:idx + batch_size]
385	for idx in range(0, len(image_files), batch_size)
386	]
387
388	pca = IncrementalPCA(n_components=2)	×
389	for batch_files in batches:	×
390	data = [	×
391	nibabel.load(_file).get_fdata().flatten()
392	for _file in batch_files
393	]
394	pca.partial_fit(data)	×
395
396	df = []	×
397	for batch_files in batches:	×
398	data = [	×
399	nibabel.load(_file).get_fdata().flatten()
400	for _file in batch_files
401	]
402	components = pca.transform(data)	×
403	info = [	×
404	parse_bids_keys(Path(_file))
405	for _file in batch_files
406	]
407	partial_df = pd.DataFrame({	×
408	"participant_id": [item["sub"] for item in info],
409	"session": [item["ses"] for item in info],
410	"run": [item["run"] for item in info],
411	"pc1": components[:, 0],
412	"pc2": components[:, 1],
413	"explained_variance_ratio_pc1": [
414	pca.explained_variance_ratio_[0],
415	] * len(info),
416	"explained_variance_ratio_pc2": [
417	pca.explained_variance_ratio_[1],
418	] * len(info),
419	})
420	df.append(partial_df)	×
421	df = pd.concat(df, ignore_index=True)	×
422	df.to_csv(pca_file, index=False, sep="\t")	×
423
424	return (pca_file, )	1✔

neurospin-deepinsight / brainprep / 21672251336

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous