diff --git a/docs/harmonisation.md b/docs/harmonisation.md index 4e2e9a61..a7e7d964 100644 --- a/docs/harmonisation.md +++ b/docs/harmonisation.md @@ -90,6 +90,7 @@ You can tune this command using additional variables and flags as detailed bello |```-demos ```| The name of the csv file containing the demographic information as detailled in the [guidelines](https://meld-graph.readthedocs.io/en/latest/prepare_data.html#prepare-the-demographic-information-required-only-to-compute-the-harmonisation-parameters). An example 'demographics_file.csv' is provided in the .| | **Optional variables** | |```--parallelise``` | use this flag to speed up the segmentation by running Freesurfer/FastSurfer on multiple subjects in parallel. | +|```--threads ``` | use this variable to speed up the segmentation of a single subject by running Freesurfer with several threads (OpenMP). Contrary to ```--parallelise```, the subjects are still processed one after another. Only applies to Freesurfer, and is ignored when ```--parallelise``` is used. | |```--fastsurfer``` | use this flag to use FastSurfer instead of Freesurfer. Requires FastSurfer installed. | |```--harmo_only``` | Use this flag to do all the processes up to the harmonisation. Usefull if you want to harmonise on some subjects but do not wish to predict on them | diff --git a/docs/interpret_results.md b/docs/interpret_results.md index 2193a4eb..82b1d3ee 100644 --- a/docs/interpret_results.md +++ b/docs/interpret_results.md @@ -114,6 +114,8 @@ After viewing the MELD PDF report, it is then important to visualise the predict The predictions are saved as NIFTI files in the folder: /output/predictions_reports//predictions +The file to open is `prediction.nii.gz` (both hemispheres combined). The subfolders of that directory contain the intermediate files created while registering the prediction back into the native volume, see [the description of script 3](https://meld-graph.readthedocs.io/en/latest/run_prediction_pipeline.html#script-3-lesions-prediction-meld-reports). + **It is important to check that the clusters detected are not due to obvious FreeSurfer reconstruction errors, scan artifacts etc.** To quality control the predicted lesion and the surfaces, you can use the HTML viewer mentioned above or use Freeview (you will need a standalone version of freeview) diff --git a/docs/run_prediction_pipeline.md b/docs/run_prediction_pipeline.md index 5107f59f..f2c06b99 100644 --- a/docs/run_prediction_pipeline.md +++ b/docs/run_prediction_pipeline.md @@ -84,8 +84,11 @@ You can tune the MELD pipeline command using additional variables and flags as d | **Optional variables** | | ```-harmo_code ``` | provide the harmonisation code if you want to harmonise your data before prediction. This requires to have [computed the harmonisation parameters](https://meld-graph.readthedocs.io/en/latest/harmonisation.html) beforehand. The harmonisation code should start with H, e.g. H1. | |```--parallelise``` | use this flag to speed up the segmentation by running Freesurfer/FastSurfer on multiple subjects in parallel. | +|```--threads ``` | use this variable to speed up the segmentation of a single subject by running Freesurfer with several threads (OpenMP). Contrary to ```--parallelise```, the subjects are still processed one after another. Only applies to Freesurfer, and is ignored when ```--parallelise``` is used. | |```--fastsurfer``` | use this flag to use FastSurfer instead of Freesurfer. (Requires FastSurfer installed for native installation). | |```--skip_feature_extraction``` | use this flag to skips the segmentation and features extraction (processes from script1). Usefull if you already have these outputs and you just want to run the preprocessing and the predictions (e.g: after harmonisation) | +|```--skip_prediction``` | use this flag to skip the predictions (processes from script3). The segmentation and the preprocessing are still run. This is different from ```--harmo_only```, which also skips the feature extraction steps that are only needed for the predictions. | +|```--skip_feature_plotting``` | use this flag to skip the plotting of the features for QC during the preprocessing (processes from script2), as this can take a long time. | |**More advanced variables** | |```--no_nifti```| use this flag to run to all the processes up saving the predictions as surface vectors in the hdf5 file. Does not produce produce nifti and pdf outputs.| |```--no_report``` | use this flag to do all the processes up to creating the prediction as a nifti file. Does not produce the pdf reports. | @@ -172,6 +175,30 @@ singularity exec meld_graph.sif /bin/bash -c "cd /app && source \$FREESURFER_HOM ::: :::: +To run a large number of subjects in parallel, e.g. on a multi-core machine or on a cluster, you can call the 3 scripts of the pipeline (see [below](#additional-information-about-the-pipeline)) yourself and parallelise the two steps that work on one subject at a time. Compared to ```--parallelise```, this lets you give several threads to each Freesurfer reconstruction and lets your job scheduler handle the queueing. + +**Step 1 - segmentation of each subject, one job per subject.** Each job uses ```--threads``` cores, so keep (number of jobs) x (number of threads) below the number of cores available: + +```bash +python scripts/new_patient_pipeline/run_script_segmentation.py -id --threads +``` + +**Step 2 - feature extraction and harmonisation, once for all the subjects.** The harmonisation parameters are computed from all the subjects of a scanner/site together, so this step cannot be split per subject: + +```bash +python scripts/new_patient_pipeline/new_pt_pipeline.py -ids input/subjects_list.txt -harmo_code -demos input/demographics_file.csv --skip_prediction --skip_feature_plotting +``` + +The Freesurfer reconstructions of step 1 are detected and skipped, the features of each subject are extracted into the feature matrix of the scanner, and the predictions and the QC feature plots are left out. + +**Step 3 - predictions and reports, one job per subject**, parallelised in the same way as step 1: + +```bash +python scripts/new_patient_pipeline/run_script_prediction.py -id -harmo_code -demos input/demographics_file.csv +``` + +Adapt this for containerization (usually apptainer/singularity for HPC) and your job scheduler (e.g., SLURM) as appropriate for your system. + ## Additional information about the pipeline The pipeline is split into 3 main scripts as illustrated below and detailed in the next sub-sections. @@ -316,9 +343,33 @@ singularity exec meld_graph.sif /bin/bash -c "cd /app && source \$FREESURFER_HOM This script : 1. Run the MELD classifier and predict lesion on new subject 2. Register the prediction back into the native nifti MRI. Results are stored in output/predictions_reports//predictions. -3. Create MELD reports with predicted lesion location on inflated brain, on native MRI and associated saliencies. Reports are stored in output/predictions_reports//predictions/reports. +3. Create MELD reports with predicted lesion location on inflated brain, on native MRI and associated saliencies. Reports are stored in output/predictions_reports//reports. + +The predictions folder contains the final prediction in a nifti volume (aligned with the input T1w) as well as intermediate files (prediction on the template fsaverage_sym surface space, the subjects native surface space, etc.): + +``` +output/predictions_reports// +├── predictions/ +│ ├── prediction.nii.gz <- final prediction (both hemispheres), in the native T1w space +│ ├── lh.prediction.nii.gz <- left hemisphere only +│ ├── rh.prediction.nii.gz <- right hemisphere only +│ ├── fsaverage_sym/ <- intermediate: prediction on the template surface +│ │ ├── lh.prediction.mgh +│ │ └── rh.prediction.mgh +│ ├── surf_native/ <- intermediate: prediction on the native surface +│ │ ├── lh.prediction.mgh +│ │ └── rh.prediction.mgh +│ ├── vol_freesurfer/ <- intermediate: prediction in the Freesurfer conformed volume +│ │ ├── lh.prediction.mgz +│ │ └── rh.prediction.mgz +│ └── vol_native/ <- intermediate: prediction in the native T1 volume, mgz format +│ ├── lh.prediction.mgz +│ └── rh.prediction.mgz +└── reports/ <- MELD PDF report and the figures it is built from +``` Notes: +- The in previous versions, these intermediate files were stored in the Freesurfer subject directory, so there was a risk of multiple predictions using the same Freesurfer outputs (e.g. with different harmonisation parameters) overwriting each other. - Features need to have been processed using script 2 and Freesurfer outputs need to be available for each subject Example to use it on one patient without harmonisation: diff --git a/meld_graph/data_preprocessing.py b/meld_graph/data_preprocessing.py index a9b2f483..218d924f 100644 --- a/meld_graph/data_preprocessing.py +++ b/meld_graph/data_preprocessing.py @@ -23,6 +23,7 @@ import meld_graph.distributedCombat as dc import meld_graph.mesh_tools as mt import meld_graph.meld_plotting as mpt +from meld_graph.hdf5_utils import open_hdf5_file class Preprocess: """ @@ -599,9 +600,9 @@ def load_covars(self, subject_ids=None, demographic_file=DEMOGRAPHIC_FEATURES_FI def save_norm_combat_parameters(self, feature, estimates, hdf5_file): """Save estimates from combat and normalisation parameters on hdf5""" if not os.path.isfile(hdf5_file): - hdf5_file_context = h5py.File(hdf5_file, "a") + hdf5_file_context = open_hdf5_file(hdf5_file, mode="a", create_parent=True) else: - hdf5_file_context = h5py.File(hdf5_file, "r+") + hdf5_file_context = open_hdf5_file(hdf5_file, mode="r+") with hdf5_file_context as f: list_params = list(set(estimates)) @@ -625,7 +626,7 @@ def save_norm_combat_parameters(self, feature, estimates, hdf5_file): def read_norm_combat_parameters(self, feature, hdf5_file): """reconstruct estimates dictionnary from the combat parameters hdf5 file""" - hdf5_file_context = h5py.File(hdf5_file, "r") + hdf5_file_context = open_hdf5_file(hdf5_file, mode="r") estimates = {} with hdf5_file_context as f: feat_dir = f[feature] diff --git a/meld_graph/evaluation.py b/meld_graph/evaluation.py index 8162fdee..b3394c62 100644 --- a/meld_graph/evaluation.py +++ b/meld_graph/evaluation.py @@ -1,5 +1,6 @@ import logging import os +import tempfile import torch import torch_geometric.data from meld_graph.dataset import GraphDataset @@ -18,6 +19,7 @@ import itertools import seaborn as sns from meld_graph.paths import MELD_DATA_PATH +from meld_graph.hdf5_utils import open_hdf5_file # for saliency - do not force people to have this try: @@ -906,27 +908,17 @@ def save_prediction(self, subject, prediction, dataset_str="prediction", dtype=N dtype = prediction.dtype filename = os.path.join(self.results_dir, f"predictions{suffix}.hdf5") - if not os.path.isfile(filename): - mode = "a" - else: - mode = "r+" - done = False - while not done: - try: - with h5py.File(filename, mode=mode) as f: - self.log.info(f"saving {dataset_str} for {subject}") - for i, hemi in enumerate(["lh", "rh"]): - shape = tuple([nvert_hemi] + list(prediction.shape[1:])) - # create dataset - dset = f.require_dataset(f"{subject}/{hemi}/{dataset_str}", shape=shape, dtype=dtype) - # save prediction in dataset - dset[:] = prediction[i * nvert_hemi : (i + 1) * nvert_hemi] - # if dataset_str == "prediction": - # save threshold as attribute in dataset - # dset.attrs["threshold"] = self.threshold - done = True - except OSError: - done = False + with open_hdf5_file(filename, mode="a", create_parent=True) as f: + self.log.info(f"saving {dataset_str} for {subject}") + for i, hemi in enumerate(["lh", "rh"]): + shape = tuple([nvert_hemi] + list(prediction.shape[1:])) + # create dataset + dset = f.require_dataset(f"{subject}/{hemi}/{dataset_str}", shape=shape, dtype=dtype) + # save prediction in dataset + dset[:] = prediction[i * nvert_hemi : (i + 1) * nvert_hemi] + # if dataset_str == "prediction": + # save threshold as attribute in dataset + # dset.attrs["threshold"] = self.threshold def load_prediction(self, subject, dataset_str="prediction", suffix=""): """ @@ -937,7 +929,7 @@ def load_prediction(self, subject, dataset_str="prediction", suffix=""): # cannot load data self.log.debug(f'file {filename} does not exist') return None - with h5py.File(filename, mode='r') as f: + with open_hdf5_file(filename, mode="r") as f: prediction = [] try: for i, hemi in enumerate(["lh", "rh"]): @@ -1155,22 +1147,23 @@ def create_surface_plots(coords, faces, overlay, flat_map=True, limits=None): else: vmin = limits[0] vmax = limits[1] - tmp_file = os.path.join(MELD_DATA_PATH,'tmp.png') - msp.plot_surf( - coords, - faces, - overlay, - flat_map=flat_map, - rotate=[90, 270], - filename=tmp_file, - vmin=vmin, - vmax=vmax, - ) - im = Image.open(tmp_file) - im = trim(im) - im = im.convert("RGBA") - im1 = np.array(im) - os.remove(tmp_file) + with tempfile.TemporaryDirectory() as tmpdir: + tmp_file = os.path.join(tmpdir,'tmp.png') + msp.plot_surf( + coords, + faces, + overlay, + flat_map=flat_map, + rotate=[90, 270], + filename=tmp_file, + vmin=vmin, + vmax=vmax, + ) + im = Image.open(tmp_file) + im = trim(im) + im = im.convert("RGBA") + im1 = np.array(im) + os.remove(tmp_file) return im1 @@ -1195,7 +1188,7 @@ def plot_roc_multiple(roc_dictionary, roc_curves_thresholds): def load_prediction(subject, hdf5, dset="prediction"): """load network predictions""" results = {} - with h5py.File(hdf5, "r") as f: + with open_hdf5_file(hdf5, mode="r") as f: for hemi in ["lh", "rh"]: results[hemi] = f[subject][hemi][dset][:] return results diff --git a/meld_graph/hdf5_utils.py b/meld_graph/hdf5_utils.py new file mode 100644 index 00000000..814bdcc3 --- /dev/null +++ b/meld_graph/hdf5_utils.py @@ -0,0 +1,57 @@ +import os +import random +import time + +import h5py + + +# Modes that require the file to already exist. h5py raises FileNotFoundError +# (a subclass of OSError) for these, which must not be retried. +_READ_MODES = ("r", "r+") + + +def open_hdf5_file(file_path, mode="r", max_retries=100, sleep_max_seconds=1.0, create_parent=False): + """Open an HDF5 file with bounded retries for transient OSErrors. + + The retries exist to survive lock contention when several processes work on + the same file. A missing file is not transient, so it fails immediately + instead of being retried for the full timeout. + + Args: + file_path: HDF5 path. + mode: Mode passed to h5py.File. + max_retries: Maximum number of open attempts. + sleep_max_seconds: Max random sleep between retries. + create_parent: If True, create parent directory before opening. + + Returns: + Open h5py.File handle. + + Raises: + FileNotFoundError: if mode requires an existing file and there is none. + OSError: if opening fails after all retries. + """ + if max_retries < 1: + raise ValueError("max_retries must be >= 1") + + if mode in _READ_MODES and not os.path.isfile(file_path): + raise FileNotFoundError(f"No such HDF5 file: {file_path} (mode={mode})") + + if create_parent: + parent_dir = os.path.dirname(file_path) + if parent_dir: + os.makedirs(parent_dir, exist_ok=True) + + last_error = None + for attempt in range(max_retries): + try: + return h5py.File(file_path, mode) + except OSError as err: + last_error = err + if attempt >= max_retries - 1: + break + time.sleep(random.random() * sleep_max_seconds) + + raise OSError( + f"Could not open HDF5 file after {max_retries} retries: {file_path} (mode={mode})" + ) from last_error diff --git a/meld_graph/meld_cohort.py b/meld_graph/meld_cohort.py index d0f8cb80..2ef6afca 100644 --- a/meld_graph/meld_cohort.py +++ b/meld_graph/meld_cohort.py @@ -21,6 +21,7 @@ import logging import meld_graph.mesh_tools as mt import scipy +from meld_graph.hdf5_utils import open_hdf5_file class MeldCohort: @@ -191,13 +192,13 @@ def _site_hdf5(self, site_code, group, write=False, hdf5_file_root=None): hdf5_file_root = self.hdf5_file_root p = os.path.join(self.data_dir, f"MELD_{site_code}", hdf5_file_root.format(site_code=site_code, group=group)) - # open existing file or create new one + if os.path.isfile(p) and not write: - f = h5py.File(p, "r") + f = open_hdf5_file(p, mode="r") elif os.path.isfile(p) and write: - f = h5py.File(p, "r+") + f = open_hdf5_file(p, mode="r+") elif not os.path.isfile(p) and write: - f = h5py.File(p, "a") + f = open_hdf5_file(p, mode="a", create_parent=True) else: f = None try: @@ -648,10 +649,7 @@ def write_feature_values(self, feature, feature_values, hemis=["lh", "rh"], hdf5 n_vert_cortex = sum(self.cohort.cortex_mask) # open hdf5 file if hdf5_file is not None: - if not os.path.isfile(hdf5_file): - hdf5_file_context = h5py.File(hdf5_file, "a") - else: - hdf5_file_context = h5py.File(hdf5_file, "r+") + hdf5_file_context = open_hdf5_file(hdf5_file, mode="a", create_parent=True) else: hdf5_file_context = self.cohort._site_hdf5( self.site_code, self.group, write=True, hdf5_file_root=hdf5_file_root diff --git a/meld_graph/mesh_tools.py b/meld_graph/mesh_tools.py index baaa522b..2c813f93 100644 --- a/meld_graph/mesh_tools.py +++ b/meld_graph/mesh_tools.py @@ -1,5 +1,6 @@ #Tools for mesh-based operations #Smoothing, +import tempfile import nibabel as nb import numpy as np import subprocess @@ -99,36 +100,34 @@ def smooth_array(input_array, neighbours, n_iter=70, cortex_mask=None): def save_mgh(filename, array, demo): """save mgh file using nibabel and imported demo mgh file""" - rand_int = np.random.randint(100) - mmap = np.memmap("/tmp/tmp" + str(rand_int), dtype="float32", mode="w+", shape=demo.get_data().shape) - mmap[:, 0, 0] = array[:] - output = nb.MGHImage(mmap, demo.affine, demo.header) - nb.save(output, filename) + with tempfile.NamedTemporaryFile() as mmap_file: + mmap = np.memmap(mmap_file.name, dtype="float32", mode="w+", shape=demo.get_data().shape) + mmap[:, 0, 0] = array[:] + output = nb.MGHImage(mmap, demo.affine, demo.header) + nb.save(output, filename) def smoothing_fs(overlay, fwhm, subject="fsaverage_sym", hemi="lh", subjects_dir=MELD_PARAMS_PATH): os.environ["SUBJECTS_DIR"] = subjects_dir """smooth surface overlay on fsaverage_sym""" - tmpdir = "/tmp/" + str(np.random.randint(1000000)) - os.mkdir(tmpdir) - dum = nb.load(os.path.join(subjects_dir, subject, "surf", hemi + ".white.avg.area.mgh")) - save_mgh(os.path.join(tmpdir, hemi + ".tmp.mgh"), overlay, dum) - subprocess.call( - "mris_fwhm --s " - + subject - + " --hemi " - + hemi - + " --cortex --smooth-only --fwhm " - + str(fwhm) - + " --i " - + os.path.join(tmpdir, hemi + ".tmp.mgh") - + " --o " - + os.path.join(tmpdir, hemi + ".sm_tmp.mgh"), - shell=True, - ) - overlay_smoothed = load_mgh(os.path.join(tmpdir, hemi + ".sm_tmp.mgh")) - subprocess.call("rm -r " + tmpdir, shell=True) + with tempfile.TemporaryDirectory() as tmpdir: + dum = nb.load(os.path.join(subjects_dir, subject, "surf", hemi + ".white.avg.area.mgh")) + save_mgh(os.path.join(tmpdir, hemi + ".tmp.mgh"), overlay, dum) + subprocess.call( + "mris_fwhm --s " + + subject + + " --hemi " + + hemi + + " --cortex --smooth-only --fwhm " + + str(fwhm) + + " --i " + + os.path.join(tmpdir, hemi + ".tmp.mgh") + + " --o " + + os.path.join(tmpdir, hemi + ".sm_tmp.mgh"), + shell=True, + ) + overlay_smoothed = load_mgh(os.path.join(tmpdir, hemi + ".sm_tmp.mgh")) return overlay_smoothed diff --git a/meld_graph/test/test_predict_newsubject.py b/meld_graph/test/test_predict_newsubject.py index d77838c6..5e0fb67f 100644 --- a/meld_graph/test/test_predict_newsubject.py +++ b/meld_graph/test/test_predict_newsubject.py @@ -17,6 +17,7 @@ from meld_graph.test.utils import create_test_demos from meld_graph.paths import MELD_DATA_PATH from meld_graph.download_data import get_test_data +from meld_graph.hdf5_utils import open_hdf5_file def get_data_parameters(): data_parameters = { @@ -33,7 +34,7 @@ def get_data_parameters(): def load_prediction(subject,hdf5): results={} - with h5py.File(hdf5, "r") as f: + with open_hdf5_file(hdf5, mode="r") as f: for hemi in ['lh','rh']: results[hemi] = f[subject][hemi]['prediction_clustered'][:] return results diff --git a/scripts/classifier/calculate_dropout_confidence.py b/scripts/classifier/calculate_dropout_confidence.py index 6b949780..9b201219 100644 --- a/scripts/classifier/calculate_dropout_confidence.py +++ b/scripts/classifier/calculate_dropout_confidence.py @@ -13,6 +13,7 @@ from meld_classifier.meld_cohort import MeldCohort from meld_graph.evaluation import Evaluator +from meld_graph.hdf5_utils import open_hdf5_file import json import os import h5py @@ -94,6 +95,6 @@ # save dropout parameters in hdf5 filename = os.path.join(eva.save_dir, "results", f"predictions{suffix}{eva.dropout_suffix}.hdf5") - with h5py.File(filename, mode='r+') as f: + with open_hdf5_file(filename, mode="r+") as f: f.attrs['dropout_p'] = args.p f.attrs['dropout_n'] = args.n \ No newline at end of file diff --git a/scripts/classifier/calculate_thresholds.py b/scripts/classifier/calculate_thresholds.py index b4c6cf1a..3943d476 100644 --- a/scripts/classifier/calculate_thresholds.py +++ b/scripts/classifier/calculate_thresholds.py @@ -5,6 +5,7 @@ from meld_classifier.meld_cohort import MeldCohort,MeldSubject import sklearn.metrics as metrics from meld_graph.evaluation import load_prediction, sens_spec_curves, roc_curves, plot_roc_multiple +from meld_graph.hdf5_utils import open_hdf5_file import pandas as pd import itertools import seaborn as sns @@ -31,7 +32,7 @@ def calculate_roc(model_path, pred_fname,cohort,thresholds): # save_dir = os.path.join(model_path, 'results_best_model') # get list of subjects - with h5py.File(os.path.join(save_dir, pred_fname), "r") as f: + with open_hdf5_file(os.path.join(save_dir, pred_fname), mode="r") as f: subjects = list(f.keys()) # load individual subject predictions over folds for subj in subjects: diff --git a/scripts/data_preparation/extract_features/io_meld.py b/scripts/data_preparation/extract_features/io_meld.py index 85a6ba08..c21239ec 100644 --- a/scripts/data_preparation/extract_features/io_meld.py +++ b/scripts/data_preparation/extract_features/io_meld.py @@ -1,8 +1,10 @@ +import tempfile import numpy as np import nibabel as nb import pandas as pd import os import h5py +from meld_graph.hdf5_utils import open_hdf5_file def load_mgh(filename): """ import mgh file using nibabel. returns flattened data array""" @@ -19,12 +21,13 @@ def import_mgh(filename): array_data=np.ndarray.flatten(mmap_data) return array_data; -def save_mgh(filename,array, demo): - """ save mgh file using nibabel and imported demo mgh file""" - mmap=np.memmap('/tmp/tmp', dtype='float32', mode='w+', shape=demo.get_data().shape) - mmap[:,0,0]=array[:] - output=nb.MGHImage(mmap, demo.affine, demo.header) - nb.save(output, filename) +def save_mgh(filename, array, demo): + """save mgh file using nibabel and imported demo mgh file""" + with tempfile.NamedTemporaryFile() as mmap_file: + mmap = np.memmap(mmap_file.name, dtype="float32", mode="w+", shape=demo.get_data().shape) + mmap[:, 0, 0] = array[:] + output = nb.MGHImage(mmap, demo.affine, demo.header) + nb.save(output, filename) #function to load subject features @@ -150,9 +153,9 @@ def save_subject(fs_id,features,medial_wall,subject_dir, demographic_file, outp hdf5_file = os.path.join(output_dir,site_code+"_"+c_p+"_featurematrix.hdf5") if hdf5_file is not None: if not os.path.isfile(hdf5_file): - f = h5py.File(hdf5_file, "a") + f = open_hdf5_file(hdf5_file, mode="a", create_parent=True) else: - f = h5py.File(hdf5_file, "r+") + f = open_hdf5_file(hdf5_file, mode="r+") for h in hemis: group=f.require_group(os.path.join(site_code,scanner,c_p,fs_id,h)) for f_name in features: diff --git a/scripts/manage_results/move_predictions_to_mgh.py b/scripts/manage_results/move_predictions_to_mgh.py index 32beb811..6108ff84 100644 --- a/scripts/manage_results/move_predictions_to_mgh.py +++ b/scripts/manage_results/move_predictions_to_mgh.py @@ -1,4 +1,5 @@ import os +import tempfile import numpy as np import h5py from meld_graph.meld_cohort import MeldCohort @@ -6,23 +7,25 @@ import nibabel as nb import argparse from meld_graph.tools_pipeline import get_m +from meld_graph.hdf5_utils import open_hdf5_file def load_prediction(subject, hdf5, prediction_name="prediction_clustered"): results = {} - with h5py.File(hdf5, "r") as f: + with open_hdf5_file(hdf5, mode="r") as f: for hemi in ["lh", "rh"]: results[hemi] = f[subject][hemi][prediction_name][:] return results def save_mgh(filename, array, demo): """save mgh file using nibabel and imported demo mgh file""" - mmap = np.memmap("/tmp/tmp", dtype="float32", mode="w+", shape=demo.get_data().shape) - mmap[:, 0, 0] = array[:] - output = nb.MGHImage(mmap, demo.affine, demo.header) - nb.save(output, filename) + with tempfile.NamedTemporaryFile() as mmap_file: + mmap = np.memmap(mmap_file.name, dtype="float32", mode="w+", shape=demo.get_data().shape) + mmap[:, 0, 0] = array[:] + output = nb.MGHImage(mmap, demo.affine, demo.header) + nb.save(output, filename) -def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, verbose=False): +def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, output_dir, verbose=False): ''' move meld predictions from hdf to mgh freesurfer volume. Outputs are saved into freesurfer subject directory inputs: subject_ids : subjects ID in an array @@ -30,10 +33,6 @@ def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, verbose=F prediction_file : hdf5 file containing the MELD predictions ''' c = MeldCohort() - # create classifier directory if not exist - classifier_dir = os.path.join(subjects_dir, subject_id, "xhemi", "classifier") - if not os.path.isdir(classifier_dir): - os.mkdir(classifier_dir) predictions = load_prediction(subject_id, prediction_file, prediction_name="cluster_thresholded_salient") for hemi in ["lh", "rh"]: prediction_h = predictions[hemi] @@ -44,7 +43,9 @@ def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, verbose=F except: print(get_m(f'Could not load {os.path.join(subjects_dir, subject_id, "xhemi", "surf_meld", f"{hemi}.on_lh.thickness.mgh")} ', subject_id, 'ERROR')) return False - filename = os.path.join(subjects_dir, subject_id, "xhemi", "classifier", f"{hemi}.prediction.mgh") + save_dir=os.path.join(output_dir, subject_id, 'predictions') + filename = os.path.join(save_dir, 'fsaverage_sym', f"{hemi}.prediction.mgh") + os.makedirs(os.path.dirname(filename), exist_ok=True) save_mgh(filename, overlay, demo) @@ -66,6 +67,7 @@ def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, verbose=F "--subjects_dir", default="", help="folder containing freesurfer outputs. It will store predictions there" ) parser.add_argument("--list_ids", default=None, help="texte file containing list of ids to process") + parser.add_argument("--output_dir", default="", help="folder to save predictions in mgh format") args = parser.parse_args() @@ -87,4 +89,4 @@ def move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, verbose=F if os.path.isfile(prediction_file): for subject_id in subjids: - move_predictions_to_mgh(subject_id, subjects_dir, prediction_file) \ No newline at end of file + move_predictions_to_mgh(subject_id, subjects_dir, prediction_file, args.output_dir) \ No newline at end of file diff --git a/scripts/manage_results/plot_prediction_report.py b/scripts/manage_results/plot_prediction_report.py index 43ea81dd..55f5106f 100644 --- a/scripts/manage_results/plot_prediction_report.py +++ b/scripts/manage_results/plot_prediction_report.py @@ -14,6 +14,7 @@ SCRIPTS_DIR, ) import os +import tempfile import json import glob import h5py @@ -36,6 +37,7 @@ from datetime import date from fpdf import FPDF from meld_graph.tools_pipeline import get_m, get_anat_files +from meld_graph.hdf5_utils import open_hdf5_file class PDF(FPDF): def lines(self): @@ -131,7 +133,7 @@ def imagey(self,im, y): def load_prediction(subject,hdf5): results={} - with h5py.File(hdf5, "r") as f: + with open_hdf5_file(hdf5, mode="r") as f: for hemi in ['lh','rh']: results[hemi] = f[subject][hemi]['prediction'][:] return results @@ -139,31 +141,32 @@ def load_prediction(subject,hdf5): def create_surface_plots(surf,prediction,c, base_size=20): """plot and reload surface images""" cmap, colors = load_cmap() - tmp_file = os.path.join(MELD_DATA_PATH,'tmp.png') - msp.plot_surf(surf['coords'], - surf['faces'],prediction, - rotate=[90], - mask=prediction==0,pvals=np.ones_like(c.cortex_mask), - colorbar=False,vmin=1,vmax=len(colors) ,cmap=cmap, - base_size=base_size, - filename=tmp_file) - im = Image.open(tmp_file) - im = trim(im) - im = im.convert("RGBA") - im1 = np.array(im) - msp.plot_surf(surf['coords'], - surf['faces'],prediction, - rotate=[270], - mask=prediction==0,pvals=np.ones_like(c.cortex_mask), - colorbar=False,vmin=1,vmax=len(colors),cmap=cmap, - base_size=base_size, - filename=tmp_file) - im = Image.open(tmp_file) - im = trim(im) - im = im.convert("RGBA") - im2 = np.array(im) - plt.close('all') - os.remove(tmp_file) + with tempfile.TemporaryDirectory() as tmpdir: + tmp_file = os.path.join(tmpdir,'tmp.png') + msp.plot_surf(surf['coords'], + surf['faces'],prediction, + rotate=[90], + mask=prediction==0,pvals=np.ones_like(c.cortex_mask), + colorbar=False,vmin=1,vmax=len(colors) ,cmap=cmap, + base_size=base_size, + filename=tmp_file) + im = Image.open(tmp_file) + im = trim(im) + im = im.convert("RGBA") + im1 = np.array(im) + msp.plot_surf(surf['coords'], + surf['faces'],prediction, + rotate=[270], + mask=prediction==0,pvals=np.ones_like(c.cortex_mask), + colorbar=False,vmin=1,vmax=len(colors),cmap=cmap, + base_size=base_size, + filename=tmp_file) + im = Image.open(tmp_file) + im = trim(im) + im = im.convert("RGBA") + im2 = np.array(im) + plt.close('all') + os.remove(tmp_file) return im1,im2 def load_cluster(file, subject): @@ -205,10 +208,11 @@ def get_cluster_location(cluster_array): def save_mgh(filename, array, demo): """save mgh file using nibabel and imported demo mgh file""" - mmap = np.memmap("/tmp/tmp", dtype="float32", mode="w+", shape=demo.get_data().shape) - mmap[:, 0, 0] = array[:] - output = nb.MGHImage(mmap, demo.affine, demo.header) - nb.save(output, filename) + with tempfile.NamedTemporaryFile() as mmap_file: + mmap = np.memmap(mmap_file.name, dtype="float32", mode="w+", shape=demo.get_data().shape) + mmap[:, 0, 0] = array[:] + output = nb.MGHImage(mmap, demo.affine, demo.header) + nb.save(output, filename) def load_cmap(): """ create the colors dictionarry for the clusters""" @@ -458,6 +462,10 @@ def generate_prediction_report( "anat": nb.load(t1_file), "pred": nb.load(prediction_file), } + if len(imgs["anat"].shape) > 3: + # if the input image has a 4th (time/frame) dimension of length 1, + # remove it (otherwise resampling and possibly other steps will fail) + imgs["anat"] = nb.funcs.squeeze_image(imgs["anat"]) # # Resample and move to same shape and affine than t1 imgs["pred"] = image.resample_img( imgs["pred"], diff --git a/scripts/manage_results/register_back_to_xhemi.py b/scripts/manage_results/register_back_to_xhemi.py index ac017a51..b8d564c0 100644 --- a/scripts/manage_results/register_back_to_xhemi.py +++ b/scripts/manage_results/register_back_to_xhemi.py @@ -18,11 +18,17 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' if not os.path.isdir(opj(subjects_dir,template)): shutil.copytree(opj(os.environ['FREESURFER_HOME'],'subjects',template), opj(subjects_dir, os.path.basename(template))) + save_dir=opj(output_dir,subject_id,'predictions') + os.makedirs(save_dir, exist_ok=True) + # Moves left hemi from fsaverage to native space # --src is the source image i.e. the map you want to move back so change to the name of the cluster map in fsaverage_sym that you want to move back # --trg is the target image i.e. the name of the map you want to create in the subject's native space # the rest is the registration files - command = f'SUBJECTS_DIR={subjects_dir} mris_apply_reg --src {subjects_dir}/{subject_id}/xhemi/classifier/lh.prediction.mgh --trg {subjects_dir}/{subject_id}/surf/lh.prediction.mgh --streg {subjects_dir}/fsaverage_sym/surf/lh.sphere.reg {subjects_dir}/{subject_id}/surf/lh.sphere.reg --nnf ' + + surf_native_dir = opj(save_dir, 'surf_native') + os.makedirs(surf_native_dir, exist_ok=True) + command = f'SUBJECTS_DIR={subjects_dir} mris_apply_reg --src {save_dir}/fsaverage_sym/lh.prediction.mgh --trg {surf_native_dir}/lh.prediction.mgh --streg {subjects_dir}/fsaverage_sym/surf/lh.sphere.reg {subjects_dir}/{subject_id}/surf/lh.sphere.reg --nnf ' proc = Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf-8') stdout, stderr= proc.communicate() if verbose: @@ -31,8 +37,8 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' print(get_m(f'COMMAND failing : {command} with error {stderr}', subject_id, 'ERROR')) return False - # Moves the right hemi back from fsaverage to native. There are 2 steps - command = f'SUBJECTS_DIR={subjects_dir} mris_apply_reg --src {subjects_dir}/{subject_id}/xhemi/classifier/rh.prediction.mgh --trg {subjects_dir}/{subject_id}/surf/rh.prediction.mgh --streg {subjects_dir}/fsaverage_sym/surf/lh.sphere.reg {subjects_dir}/{subject_id}/xhemi/surf/lh.fsaverage_sym.sphere.reg --nnf' + # Moves the right hemi back from fsaverage to native. + command = f'SUBJECTS_DIR={subjects_dir} mris_apply_reg --src {save_dir}/fsaverage_sym/rh.prediction.mgh --trg {surf_native_dir}/rh.prediction.mgh --streg {subjects_dir}/fsaverage_sym/surf/lh.sphere.reg {subjects_dir}/{subject_id}/xhemi/surf/lh.fsaverage_sym.sphere.reg --nnf' proc = Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf-8') stdout, stderr= proc.communicate() if verbose: @@ -41,11 +47,16 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' print(get_m(f'COMMAND failing : {command} with error {stderr}', subject_id, 'ERROR')) return False + vol_freesurfer_dir = opj(save_dir, 'vol_freesurfer') + vol_native_dir = opj(save_dir, 'vol_native') + os.makedirs(vol_freesurfer_dir, exist_ok=True) + os.makedirs(vol_native_dir, exist_ok=True) + #correct from interpolation error for hemi in ['lh','rh']: #map from surface back to vol - command = f'SUBJECTS_DIR={subjects_dir} mri_surf2vol --identity {subject_id} --template {subjects_dir}/{subject_id}/mri/T1.mgz --o {subjects_dir}/{subject_id}/mri/{hemi}.prediction.mgz --hemi {hemi} --surfval {subjects_dir}/{subject_id}/surf/{hemi}.prediction.mgh --fillribbon' + command = f'SUBJECTS_DIR={subjects_dir} mri_surf2vol --identity {subject_id} --template {subjects_dir}/{subject_id}/mri/T1.mgz --o {vol_freesurfer_dir}/{hemi}.prediction.mgz --hemi {hemi} --surfval {surf_native_dir}/{hemi}.prediction.mgh --fillribbon' proc = Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf-8') stdout, stderr= proc.communicate() if verbose: @@ -55,7 +66,7 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' return False #register back to original volume - command = f'SUBJECTS_DIR={subjects_dir} mri_vol2vol --mov {subjects_dir}/{subject_id}/mri/{hemi}.prediction.mgz --targ {subjects_dir}/{subject_id}/mri/orig/001.mgz --regheader --o {subjects_dir}/{subject_id}/mri/{hemi}.prediction.mgz --nearest' + command = f'SUBJECTS_DIR={subjects_dir} mri_vol2vol --mov {vol_freesurfer_dir}/{hemi}.prediction.mgz --targ {subjects_dir}/{subject_id}/mri/orig/001.mgz --regheader --o {vol_native_dir}/{hemi}.prediction.mgz --nearest' proc = Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf-8') stdout, stderr= proc.communicate() if verbose: @@ -65,7 +76,7 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' return False #convert to nifti - command = f'SUBJECTS_DIR={subjects_dir} mri_convert {subjects_dir}/{subject_id}/mri/{hemi}.prediction.mgz {subjects_dir}/{subject_id}/mri/{hemi}.prediction.nii.gz -rt nearest' + command = f'SUBJECTS_DIR={subjects_dir} mri_convert {vol_native_dir}/{hemi}.prediction.mgz {vol_native_dir}/{hemi}.prediction.nii.gz -rt nearest' proc = Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE, encoding='utf-8') stdout, stderr= proc.communicate() if verbose: @@ -74,12 +85,10 @@ def register_subject_to_xhemi(subject_id, subjects_dir, output_dir, template = ' print(get_m(f'COMMAND failing : {command} with error {stderr}', subject_id, 'ERROR')) return False - #move files - save_dir=opj(output_dir,subject_id,'predictions') - os.makedirs(save_dir, exist_ok=True) + - shutil.move(f'{subjects_dir}/{subject_id}/mri/lh.prediction.nii.gz', f'{save_dir}/lh.prediction.nii.gz') - shutil.move(f'{subjects_dir}/{subject_id}/mri/rh.prediction.nii.gz', f'{save_dir}/rh.prediction.nii.gz') + shutil.move(f'{vol_native_dir}/lh.prediction.nii.gz', f'{save_dir}/lh.prediction.nii.gz') + shutil.move(f'{vol_native_dir}/rh.prediction.nii.gz', f'{save_dir}/rh.prediction.nii.gz') #combine vols from left and right hemis command=f'mri_concat --i {save_dir}/lh.prediction.nii.gz --i {save_dir}/rh.prediction.nii.gz --o {save_dir}/prediction.nii.gz --combine' diff --git a/scripts/new_patient_pipeline/new_pt_pipeline.py b/scripts/new_patient_pipeline/new_pt_pipeline.py index e4492794..63cd63ae 100644 --- a/scripts/new_patient_pipeline/new_pt_pipeline.py +++ b/scripts/new_patient_pipeline/new_pt_pipeline.py @@ -52,11 +52,22 @@ def flush(self): action="store_true", ) parser.add_argument("--parallelise", - help="parallelise segmentation", + help="run the freesurfer segmentation of several subjects at the same time, " + "as one single-threaded recon-all process per subject. " + "--threads is ignored in this mode.", required=False, default=False, action="store_true", ) + parser.add_argument("--threads", + help="number of threads (openmp) to use within a single freesurfer segmentation. " + "This runs one recon-all call at a time and speeds up that one call; " + "it does not process several subjects at once (see --parallelise), " + "and it is ignored when --parallelise is given.", + required=False, + default=1, + type=int, + ) parser.add_argument('-demos', '--demographic_file', type=str, help='provide the demographic files for the harmonisation', @@ -73,6 +84,14 @@ def flush(self): action="store_true", help='Skip the segmentation and extraction of the MELD features', ) + parser.add_argument('--skip_prediction', + action="store_true", + help='Skip the prediction step. This is different from --harmo_only (which skips some feature extraction steps necessary for predictions).', + ) + parser.add_argument('--skip_feature_plotting', + action="store_true", + help='Skip the plotting of the features for QC during preprocessing (because this can take a long time)', + ) parser.add_argument('--no_nifti', action="store_true", default=False, @@ -143,6 +162,7 @@ def flush(self): sub_id=args.id, use_parallel=args.parallelise, use_fastsurfer=args.fastsurfer, + threads=args.threads, verbose = args.debug_mode ) if result == False: @@ -159,11 +179,12 @@ def flush(self): list_ids=args.list_ids, sub_id=args.id, harmonisation_only = args.harmo_only, + skip_feature_plotting=args.skip_feature_plotting ) #--------------------------------------------------------------------------------- ### PREDICTION ### - if not args.harmo_only: + if not args.harmo_only and not args.skip_prediction: print(get_m(f'Call script prediction', None, 'SCRIPT 3')) result = run_script_prediction( harmo_code = args.harmo_code, diff --git a/scripts/new_patient_pipeline/run_script_prediction.py b/scripts/new_patient_pipeline/run_script_prediction.py index bd184253..6b1b1c99 100644 --- a/scripts/new_patient_pipeline/run_script_prediction.py +++ b/scripts/new_patient_pipeline/run_script_prediction.py @@ -139,6 +139,7 @@ def run_script_prediction(list_ids=None, sub_id=None, harmo_code='noHarmo', no_p result = move_predictions_to_mgh(subject_id=subject_id, subjects_dir=subjects_dir, prediction_file=prediction_file, + output_dir=predictions_output_dir, verbose=verbose) if result == False: print(get_m(f'One step of the pipeline has failed. Process has been aborted for this subject', subject_id, 'ERROR')) diff --git a/scripts/new_patient_pipeline/run_script_preprocessing.py b/scripts/new_patient_pipeline/run_script_preprocessing.py index 63444ea7..053a9586 100644 --- a/scripts/new_patient_pipeline/run_script_preprocessing.py +++ b/scripts/new_patient_pipeline/run_script_preprocessing.py @@ -58,7 +58,14 @@ def check_demographic_file(demographic_file, subject_ids): if len(np.unique(ages))<=1: sys.exit(get_m(f'There is no variance in the ages provided. Harmonisation will fail', None, 'ERROR')) -def run_data_processing_new_subjects(subject_ids, harmo_code, compute_harmonisation = False, harmonisation_only = False, demographic_file=None, output_dir=BASE_PATH, withoutflair=False): +def run_data_processing_new_subjects(subject_ids, + harmo_code, + compute_harmonisation = False, + harmonisation_only = False, + demographic_file=None, + output_dir=BASE_PATH, + withoutflair=False, + skip_feature_plotting=False): # Set features and smoothed values if withoutflair: @@ -204,19 +211,28 @@ def run_data_processing_new_subjects(subject_ids, harmo_code, compute_harmonisat ### PLOT FEATURES FOR QC #----------------------------------------------------------------- - features_to_plot = [ ".inter_z.asym.intra_z" + feature for feature in features_combat] - c_norm = MeldCohort(hdf5_file_root="{site_code}_{group}_featurematrix_combat.hdf5", dataset=tmp.name, data_dir=BASE_PATH) - plot = Preprocess(c_norm, - site_codes=[harmo_code], - write_output_file=None, - data_dir=output_dir) - - print(get_m(f'Plot features to QC', None, 'STEP')) - plot.plot_subject_features(features_to_plot) + if not skip_feature_plotting: + features_to_plot = [ ".inter_z.asym.intra_z" + feature for feature in features_combat] + c_norm = MeldCohort(hdf5_file_root="{site_code}_{group}_featurematrix_combat.hdf5", dataset=tmp.name, data_dir=BASE_PATH) + plot = Preprocess(c_norm, + site_codes=[harmo_code], + write_output_file=None, + data_dir=output_dir) + + print(get_m(f'Plot features to QC', None, 'STEP')) + plot.plot_subject_features(features_to_plot) tmp.close() -def run_script_preprocessing(list_ids=None, sub_id=None, harmo_code='noHarmo', output_dir=BASE_PATH, demographic_file=None, harmonisation_only=False, withoutflair=False, verbose=False): +def run_script_preprocessing(list_ids=None, + sub_id=None, + harmo_code='noHarmo', + output_dir=BASE_PATH, + demographic_file=None, + harmonisation_only=False, + withoutflair=False, + skip_feature_plotting=False, + verbose=False): harmo_code = str(harmo_code) subject_id=None subject_ids=None @@ -268,7 +284,8 @@ def run_script_preprocessing(list_ids=None, sub_id=None, harmo_code='noHarmo', o demographic_file=demographic_file, harmonisation_only = harmonisation_only, output_dir=output_dir, - withoutflair=withoutflair) + withoutflair=withoutflair, + skip_feature_plotting=skip_feature_plotting) if __name__ == '__main__': @@ -308,6 +325,11 @@ def run_script_preprocessing(list_ids=None, sub_id=None, harmo_code='noHarmo', o default=False, help="do not use flair information", ) + parser.add_argument("--skip_feature_plotting", + action="store_true", + default=False, + help="skip the plotting of the features for QC during preprocessing", + ) parser.add_argument("--debug_mode", help="mode to debug error", required=False, @@ -359,5 +381,6 @@ def run_script_preprocessing(list_ids=None, sub_id=None, harmo_code='noHarmo', o demographic_file=args.demographic_file, harmonisation_only = args.harmo_only, withoutflair=args.withoutflair, + skip_feature_plotting=args.skip_feature_plotting, verbose = args.debug_mode, ) \ No newline at end of file diff --git a/scripts/new_patient_pipeline/run_script_segmentation.py b/scripts/new_patient_pipeline/run_script_segmentation.py index d812eb6a..42d0519e 100644 --- a/scripts/new_patient_pipeline/run_script_segmentation.py +++ b/scripts/new_patient_pipeline/run_script_segmentation.py @@ -52,6 +52,23 @@ def check_xhemi_outputs(): #TODO pass +# workaround for the FLAIR pial reconstruction bug in freesurfer +PLACE_MM_PIAL_SURF_OPTS = "PlaceMMPialSurf --mm_min_inside 50 --mm_max_inside 200 --mm_min_outside 10 --mm_max_outside 50" + +def write_expert_opts(subject_id, lines): + """write a per-subject recon-all expert options file and return its path + + A per-subject file is used rather than a global-expert-options.txt in the + SUBJECTS_DIR, because the latter is shared by every subject: parallel runs + would race on it, and it would silently apply to subjects processed later + that do not need it. + """ + expert_opts_path = opj(MELD_DATA_PATH, 'output', 'fs_outputs', 'expertopts', f'{subject_id}.opts') + os.makedirs(os.path.dirname(expert_opts_path), exist_ok=True) + with open(expert_opts_path, 'w') as f: + f.write('\n'.join(lines) + '\n') + return expert_opts_path + def fastsurfer_subject(subject, fs_folder, verbose=False): # run fastsurfer segmentation on 1 subject subject_id = subject['id'] @@ -128,7 +145,7 @@ def fastsurfer_flair(subject, fs_folder, verbose=False): print(get_m(f'COMMAND failing : {command} with error {stderr}', None, 'ERROR')) return False -def freesurfer_subject(subject, fs_folder, verbose=False): +def freesurfer_subject(subject, fs_folder, threads=1, freesurfer_args=None, verbose=False): #run freesurfer recon-all segmentation on 1 subject subject_id = subject['id'] @@ -151,17 +168,25 @@ def freesurfer_subject(subject, fs_folder, verbose=False): # setup cortical segmentation command if subject_flair_path != None: print(get_m('Segmentation using T1 and FLAIR with Freesurfer', subject_id, 'STEP 1')) + + # apply the "PlaceMMPialSurf" bugfix, needed for the FLAIRpial step + expert_opts_path = write_expert_opts(subject_id, [PLACE_MM_PIAL_SURF_OPTS]) + command = format( - "$FREESURFER_HOME/bin/recon-all -sd {} -s {} -i {} -FLAIR {} -FLAIRpial -all".format( - fs_folder, subject_id, subject_t1_path, subject_flair_path + "$FREESURFER_HOME/bin/recon-all -sd {} -s {} -i {} -FLAIR {} -FLAIRpial -all -threads {} -expert {}".format( + fs_folder, subject_id, subject_t1_path, subject_flair_path, threads, expert_opts_path ) ) else: print(get_m('Segmentation using T1 only with Freesurfer', subject_id, 'STEP 1')) command = format( - "$FREESURFER_HOME/bin/recon-all -sd {} -s {} -i {} -all".format(fs_folder, subject_id, subject_t1_path) + "$FREESURFER_HOME/bin/recon-all -sd {} -s {} -i {} -all -threads {}".format(fs_folder, subject_id, subject_t1_path, threads) ) + if freesurfer_args is not None: + command += ' ' + ' '.join(freesurfer_args) + print(get_m(f'Additional freesurfer arguments added: {freesurfer_args}', subject_id, 'INFO')) + # call Freesurfer print(get_m('Start cortical parcellation (up to 6h). Please wait', subject_id, 'INFO')) print(get_m(f'Results will be stored in {fs_folder}', subject_id, 'INFO')) @@ -225,7 +250,7 @@ def extract_features(subject_id, fs_folder, output_dir, verbose=False): if result == False: return False -def run_subjects_segmentation_parallel(subject_ids, num_procs=10, harmo_code="noHarmo", use_fastsurfer=False, verbose=False): +def run_subjects_segmentation_parallel(subject_ids, num_procs=10, harmo_code="noHarmo", use_fastsurfer=False, freesurfer_args=None, verbose=False): # parallel version of the pipeline, finish each stage for all subjects first ### SEGMENTATION ### @@ -278,7 +303,7 @@ def run_subjects_segmentation_parallel(subject_ids, num_procs=10, harmo_code="no ## processing with freesurfer: segmentation pool = multiprocessing.Pool(processes=num_procs, initializer=init, initargs=[multiprocessing.Lock()]) subject_ids_failed=[] - for i,result in enumerate(pool.imap(partial(freesurfer_subject, fs_folder=fs_folder, verbose=verbose), subjects_dict)): + for i,result in enumerate(pool.imap(partial(freesurfer_subject, fs_folder=fs_folder, freesurfer_args=freesurfer_args, verbose=verbose), subjects_dict)): if result==False: print(get_m(f'Subject removed from futur process because a step in the pipeline failed', subject_ids[i], 'ERROR')) subject_ids_failed.append(subject_ids[i]) @@ -307,7 +332,7 @@ def run_subjects_segmentation_parallel(subject_ids, num_procs=10, harmo_code="no return subject_ids -def run_subject_segmentation(subject_id, harmo_code="noHarmo", use_fastsurfer=False, verbose=False): +def run_subject_segmentation(subject_id, harmo_code="noHarmo", use_fastsurfer=False, threads=1, freesurfer_args=None, verbose=False): # pipeline to segment the brain, exract surface-based features for 1 subject ### SEGMENTATION ### @@ -341,7 +366,7 @@ def run_subject_segmentation(subject_id, harmo_code="noHarmo", use_fastsurfer=Fa else: ## processing with freesurfer: segmentation init(multiprocessing.Lock()) - result = freesurfer_subject(subject_dict,fs_folder, verbose=verbose) + result = freesurfer_subject(subject_dict, fs_folder, threads=threads, freesurfer_args=freesurfer_args, verbose=verbose) if result == False: return False @@ -352,7 +377,7 @@ def run_subject_segmentation(subject_id, harmo_code="noHarmo", use_fastsurfer=Fa return False -def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', use_parallel=False, use_fastsurfer=False, verbose=False ): +def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', use_parallel=False, use_fastsurfer=False, threads=1, freesurfer_args=None, verbose=False ): harmo_code = str(harmo_code) subject_id=None subject_ids=None @@ -375,7 +400,7 @@ def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', us if subject_id != None: #launch segmentation and feature extraction for 1 subject - result = run_subject_segmentation(subject_id, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, verbose=verbose) + result = run_subject_segmentation(subject_id, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, threads=threads, freesurfer_args=freesurfer_args, verbose=verbose) if result == False: print(get_m(f'One step of the pipeline has failed. Process has been aborted for this subject', subject_id, 'ERROR')) return False @@ -383,7 +408,7 @@ def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', us if use_parallel: #launch segmentation and feature extraction in parallel print(get_m(f'Run subjects in parallel', None, 'INFO')) - subject_ids_succeed = run_subjects_segmentation_parallel(subject_ids, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, verbose=verbose) + subject_ids_succeed = run_subjects_segmentation_parallel(subject_ids, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, freesurfer_args=freesurfer_args, verbose=verbose) subject_ids_failed= list(set(subject_ids).difference(subject_ids_succeed)) if len(subject_ids_failed): print(get_m(f'One step of the pipeline has failed. Process has been aborted for subjects {subject_ids_failed}', None, 'ERROR')) @@ -393,7 +418,7 @@ def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', us print(get_m(f'Run subjects one after another', None, 'INFO')) subject_ids_failed=[] for subj in subject_ids: - result = run_subject_segmentation(subj, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, verbose=verbose) + result = run_subject_segmentation(subj, harmo_code = harmo_code, use_fastsurfer = use_fastsurfer, threads=threads, freesurfer_args=freesurfer_args, verbose=verbose) if result == False: print(get_m(f'One step of the pipeline has failed. Process has been aborted for this subject', subj, 'ERROR')) subject_ids_failed.append(subj) @@ -432,11 +457,30 @@ def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', us action="store_true", ) parser.add_argument("--parallelise", - help="parallelise segmentation", + help="run the freesurfer segmentation of several subjects at the same time, " + "as one single-threaded recon-all process per subject. " + "--threads is ignored in this mode.", required=False, default=False, action="store_true", ) + parser.add_argument("--threads", + help="number of threads (openmp) to use within a single freesurfer segmentation. " + "This runs one recon-all call at a time and speeds up that one call; " + "it does not process several subjects at once (see --parallelise), " + "and it is ignored when --parallelise is given.", + required=False, + default=1, + type=int, + ) + parser.add_argument("--freesurfer_args", + help="additional arguments to pass to freesurfer recon-all. " + "Quote each one so it is not parsed as an option of this script, " + "e.g. --freesurfer_args '-cw256' '-notal-check'", + required=False, + default=None, + type=str, + nargs="*") parser.add_argument("--debug_mode", help="mode to debug error", required=False, @@ -483,8 +527,10 @@ def run_script_segmentation(list_ids=None, sub_id=None, harmo_code='noHarmo', us harmo_code = args.harmo_code, list_ids=args.list_ids, sub_id=args.id, - use_parallel=args.parallelise, + use_parallel=args.parallelise, use_fastsurfer=args.fastsurfer, + threads=args.threads, + freesurfer_args=args.freesurfer_args, verbose = args.debug_mode )