From c6b561a8c97c563c3056c5f3c6a6621ba667b664 Mon Sep 17 00:00:00 2001 From: Quentin Dessain <5383293+Hyedryn@users.noreply.github.com> Date: Sun, 7 Apr 2024 11:32:03 +0200 Subject: [PATCH] New masking for T1 and DWI, new registration technique for DWI->T1 transform (#40) * New masking strategy for preproc * Changed connectivity in clean_mask * Edited slurm command for better lemaitre4 HCP cluster compatibility * Improved T1 mask for synb0disco * Added transform from DWI2T1 using fsl epi_reg * Improved T1 mask with mri_synth_strip --- elikopy/core.py | 8 +- elikopy/individual_subject_processing.py | 403 +++++++++++++---------- elikopy/registration.py | 94 +++++- elikopy/utils.py | 8 +- 4 files changed, 323 insertions(+), 190 deletions(-) diff --git a/elikopy/core.py b/elikopy/core.py index 2d1e540..5db2492 100644 --- a/elikopy/core.py +++ b/elikopy/core.py @@ -459,7 +459,7 @@ def patient_list(self, folder_path=None, bids_path=None, reverseEncoding=True): f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Patient list generated\n") f.close() - def preproc(self, folder_path=None, reslice=False, reslice_addSlice=False, denoising=False, denoising_algorithm="mppca_mrtrix", gibbs=False, topup=False, topupConfig=None, forceSynb0DisCo=False, useGPUsynb0DisCo=False, eddy=False, biasfield=False, biasfield_bsplineFitting=[100,3], biasfield_convergence=[1000,0.001], patient_list_m=None, starting_state=None, bet_median_radius=2, bet_numpass=1, bet_dilate=2, static_files_path=None, cuda=None, cuda_name="eddy_cuda10.1", s2v=[0,5,1,'trilinear'], olrep=[False, 4, 250, 'sw'], eddy_additional_arg="", slurm=None, slurm_email=None, slurm_timeout=None, cpus=None, slurm_mem=None, qc_reg=True, niter=5, slspec_gc_path=None, report=True): + def preproc(self, folder_path=None, reslice=False, reslice_addSlice=False, denoising=False, gibbs=False, topup=False, topupConfig=None, forceSynb0DisCo=False, useGPUsynb0DisCo=False, eddy=False, biasfield=False, biasfield_bsplineFitting=[100,3], biasfield_convergence=[1000,0.001], patient_list_m=None, starting_state=None, bet_median_radius=2, bet_numpass=1, bet_dilate=2, static_files_path=None, cuda=None, cuda_name="eddy_cuda10.1", s2v=[0,5,1,'trilinear'], olrep=[False, 4, 250, 'sw'], eddy_additional_arg="", slurm=None, slurm_email=None, slurm_timeout=None, cpus=None, slurm_mem=None, qc_reg=True, niter=5, slspec_gc_path=None, report=True): """ Performs data preprocessing. By default only the brain extraction is enabled. Optional preprocessing steps include : reslicing, denoising, gibbs ringing correction, susceptibility field estimation, EC-induced distortions and motion correction, bias field correction. The results are stored in the preprocessing subfolder of each study subject /subjects//dMRI/preproc. @@ -499,8 +499,6 @@ def preproc(self, folder_path=None, reslice=False, reslice_addSlice=False, denoi """ assert starting_state in (None, "denoising", "gibbs", "topup", "eddy", "biasfield", "report", "post_report", "topup_synb0DisCo_Registration", "topup_synb0DisCo_Inference", "topup_synb0DisCo_Apply", "topup_synb0DisCo_topup"), 'invalid starting state!' - if denoising==True: - assert denoising_algorithm in ["patch2self", "mppca_mrtrix", "mppca_dipy"], 'invalid denoising algorithm!' if starting_state=="denoising": assert denoising == True, 'if starting_state is denoising, denoising must be True!' if starting_state=="gibbs": @@ -548,7 +546,7 @@ def preproc(self, folder_path=None, reslice=False, reslice_addSlice=False, denoi p_job = { "wrap": "export OMP_NUM_THREADS="+str(tot_cpu)+" ; export FSLPARALLEL="+str(tot_cpu)+" ; python -c 'from elikopy.individual_subject_processing import preproc_solo; preproc_solo(\"" + folder_path + "/\",\"" + p + "\",eddy=" + str( eddy) + ",biasfield=" + str(biasfield) + ",biasfield_convergence=[" + str(biasfield_convergence[0]) + "," + str(biasfield_convergence[1]) + "],biasfield_bsplineFitting=[" + str(biasfield_bsplineFitting[0]) + "," + str(biasfield_bsplineFitting[1]) + "],denoising=" + str( - denoising) + ",denoising_algorithm=\"" + str(denoising_algorithm) +"\",reslice=" + str(reslice) + ",reslice_addSlice=" + str(reslice_addSlice) + ",gibbs=" + str( + denoising) + ",reslice=" + str(reslice) + ",reslice_addSlice=" + str(reslice_addSlice) + ",gibbs=" + str( gibbs) + ",topup=" + str(topup) + ",forceSynb0DisCo=" + str(forceSynb0DisCo) + ",useGPUsynb0DisCo=" + str(useGPUsynb0DisCo) + ",topupConfig=\"" + str(topupConfig) + "\",starting_state=\"" + str(starting_state) + "\",static_files_path=\""+ static_files_path +"\" ,bet_median_radius=" + str( bet_median_radius) + ",bet_dilate=" + str(bet_dilate) + ", qc_reg=" + str(qc_reg) + ", report=" + str(report) + ", slspec_gc_path=" + str(slspec_gc_path) + ", core_count=" + str(core_count)+ ", niter=" + str(niter)+",bet_numpass=" + str( bet_numpass) + ",cuda=" + str(cuda) + ",cuda_name=\"" + str(cuda_name) + "\",s2v=[" + str(s2v[0]) + "," + str(s2v[1]) + "," + str(s2v[2]) + ",\"" + str(s2v[3]) + "\"],olrep=[" + str(olrep[0]) + "," + str(olrep[1]) + "," + str(olrep[2]) + ",\"" + str(olrep[3]) + "\"], " + "eddy_additional_arg=\"" + str(eddy_additional_arg) + "\" )'", @@ -595,7 +593,7 @@ def preproc(self, folder_path=None, reslice=False, reslice_addSlice=False, denoi f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Successfully submited job %s using slurm\n" % p_job_id) else: core_count = 1 if cpus is None else cpus - preproc_solo(folder_path + "/",p,reslice=reslice,reslice_addSlice=reslice_addSlice,denoising=denoising, denoising_algorithm=denoising_algorithm, gibbs=gibbs, + preproc_solo(folder_path + "/",p,reslice=reslice,reslice_addSlice=reslice_addSlice,denoising=denoising, gibbs=gibbs, topup=topup, topupConfig=topupConfig, forceSynb0DisCo=forceSynb0DisCo, useGPUsynb0DisCo=useGPUsynb0DisCo, eddy=eddy,biasfield=biasfield, biasfield_bsplineFitting=biasfield_bsplineFitting, biasfield_convergence=biasfield_convergence, starting_state=starting_state, diff --git a/elikopy/individual_subject_processing.py b/elikopy/individual_subject_processing.py index e18261b..be93292 100644 --- a/elikopy/individual_subject_processing.py +++ b/elikopy/individual_subject_processing.py @@ -5,6 +5,7 @@ import sys import numpy as np import math +from scipy.ndimage.morphology import binary_dilation import subprocess from elikopy.utils import makedir @@ -13,7 +14,7 @@ print = functools.partial(print, flush=True) -def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoising=False, denoising_algorithm="mppca_mrtrix", gibbs=False, topup=False, topupConfig=None, forceSynb0DisCo=False, useGPUsynb0DisCo=False, eddy=False, biasfield=False, biasfield_bsplineFitting=[100,3], biasfield_convergence=[1000,0.001], static_files_path=None, starting_state=None, bet_median_radius=2, bet_numpass=1, bet_dilate=2, cuda=False, cuda_name="eddy_cuda10.1", s2v=[0,5,1,'trilinear'], olrep=[False, 4, 250, 'sw'], eddy_additional_arg="", qc_reg=True, core_count=1, niter=5, report=True, slspec_gc_path=None): +def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoising=False, gibbs=False, topup=False, topupConfig=None, forceSynb0DisCo=False, useGPUsynb0DisCo=False, eddy=False, biasfield=False, biasfield_bsplineFitting=[100,3], biasfield_convergence=[1000,0.001], static_files_path=None, starting_state=None, bet_median_radius=2, bet_numpass=1, bet_dilate=2, cuda=False, cuda_name="eddy_cuda10.1", s2v=[0,5,1,'trilinear'], olrep=[False, 4, 250, 'sw'], eddy_additional_arg="", qc_reg=True, core_count=1, niter=5, report=True, slspec_gc_path=None): """ Performs data preprocessing on a single subject. By default only the brain extraction is enabled. Optional preprocessing steps include : reslicing, denoising, gibbs ringing correction, susceptibility field estimation, EC-induced distortions and motion correction, bias field correction. The results are stored in the preprocessing subfolder of the study subject /subjects//dMRI/preproc. @@ -52,7 +53,6 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin assert starting_state in (None,"None", "denoising", "gibbs", "topup", "eddy", "biasfield", "report", "topup_synb0DisCo_Registration", "topup_synb0DisCo_Inference", "topup_synb0DisCo_Apply", "topup_synb0DisCo_topup"), 'invalid starting state!' if starting_state == "denoising": assert denoising == True, 'if starting_state is denoising, denoising must be True!' - assert denoising_algorithm in ["patch2self", "mppca_mrtrix", "mppca_dipy"], 'invalid denoising algorithm!' if starting_state == "gibbs": assert gibbs == True, 'if starting_state is gibbs, gibbs must be True!' if starting_state in ("topup", "topup_synb0DisCo_Registration", "topup_synb0DisCo_Inference", "topup_synb0DisCo_Apply", "topup_synb0DisCo_topup"): @@ -105,12 +105,11 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin from dipy.io.image import load_nifti, save_nifti from dipy.segment.mask import median_otsu - from dipy.denoise.localpca import mppca nifti_path = folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.nii.gz' if (starting_state == None): data, affine, voxel_size = load_nifti(nifti_path, return_voxsize=True) - + curr_dmri = data reslice_path = folder_path + '/subjects/' + patient_path + "/dMRI/preproc/reslice" if reslice and starting_state == None: makedir(reslice_path, folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", log_prefix) @@ -122,6 +121,8 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin if reslice_addSlice: data = np.insert(data, np.size(data,2), 0, axis=2) + curr_dmri = data + nifti_path = reslice_path + '/' + patient_path + '_reslice.nii.gz' save_nifti(reslice_path + '/' + patient_path + '_reslice.nii.gz', data, affine) print("[" + log_prefix + "] " + datetime.datetime.now().strftime( "%d.%b %Y %H:%M:%S") + ": Reslice completed for patient %s \n" % p) @@ -129,25 +130,16 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( "%d.%b %Y %H:%M:%S") + ": Reslice completed for patient %s \n" % p) f.close() + if reslice: + nifti_path = reslice_path + '/' + patient_path + '_reslice.nii.gz' if starting_state == None: from elikopy.utils import clean_mask - b0_mask, mask = median_otsu(data, median_radius=bet_median_radius, numpass=bet_numpass, vol_idx=range(0, np.shape(data)[3]), dilate=bet_dilate) + _, mask = median_otsu(curr_dmri, median_radius=bet_median_radius, numpass=bet_numpass, vol_idx=range(0, np.shape(curr_dmri)[3]), dilate=bet_dilate) mask = clean_mask(mask) - # Apply cleaned mask to b0_mask - b0_mask = b0_mask * mask[..., None] save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz',mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz',b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) - - _, mask_nodilate = median_otsu(data, median_radius=bet_median_radius, numpass=bet_numpass, - vol_idx=range(0, np.shape(data)[3]), dilate=None) - mask_nodilate = clean_mask(mask_nodilate) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', - mask_nodilate.astype(np.float32), affine) print("[" + log_prefix + "] " + datetime.datetime.now().strftime( "%d.%b %Y %H:%M:%S") + ": Brain extraction completed for patient %s \n" % p) @@ -157,21 +149,15 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin f.close() if not denoising and not eddy and not gibbs and not topup and not biasfield: - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', mask.astype(np.float32), affine) + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz', curr_dmri.astype(np.float32), affine) shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval") shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bvec", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bvec") - if denoising_algorithm in ['mppca_mrtrix', 'mppca_dipy']: - denoising_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/mppca' - denoising_ext = '_mppca.nii.gz' - elif denoising_algorithm == "patch2self": - from dipy.denoise.patch2self import patch2self + denoising_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/mppca' + denoising_ext = '_mppca.nii.gz' - denoising_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/patch2self' - denoising_ext = '_patch2self.nii.gz' if denoising and starting_state!="gibbs" and starting_state!="eddy" and (starting_state not in ("topup", "topup_synb0DisCo_Registration", "topup_synb0DisCo_Inference", "topup_synb0DisCo_Apply", "topup_synb0DisCo_topup")) and starting_state!="biasfield" and starting_state!="report": print("[" + log_prefix + "] " + datetime.datetime.now().strftime( @@ -184,40 +170,17 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( "%d.%b %Y %H:%M:%S") + ": Denoising launched for patient %s \n" % p) - if (starting_state == "denoising"): - mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz' - b0_mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz' - b0_mask, affine, voxel_size = load_nifti(b0_mask_path, return_voxsize=True) - mask, _ = load_nifti(mask_path) - - if denoising_algorithm == "mppca_mrtrix": - import subprocess - bet_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz' - bashCommand = "dwidenoise -nthreads " + str(core_count) + " " + bet_path + \ - " " + denoising_path + '/' + patient_path + '_mppca.nii.gz' +\ - " -noise " + denoising_path + '/' + patient_path + '_sigmaNoise.nii.gz -force' - - process = subprocess.Popen(bashCommand, universal_newlines=True, shell=True, stdout=f, - stderr=subprocess.STDOUT) - - output, error = process.communicate() + import subprocess + bashCommand = "dwidenoise -nthreads " + str(core_count) + " " + nifti_path + \ + " " + denoising_path + '/' + patient_path + '_mppca.nii.gz' +\ + " -noise " + denoising_path + '/' + patient_path + '_sigmaNoise.nii.gz -force' - denoised, _ = load_nifti(denoising_path + '/' + patient_path + '_mppca.nii.gz') + process = subprocess.Popen(bashCommand, universal_newlines=True, shell=True, stdout=f, + stderr=subprocess.STDOUT) - elif denoising_algorithm == "mppca_dipy": - pr = math.ceil((np.shape(b0_mask)[3] ** (1 / 3) - 1) / 2) - denoised, sigma = mppca(b0_mask, patch_radius=pr, return_sigma=True, mask = mask) - save_nifti(denoising_path + '/' + patient_path + '_sigmaNoise.nii.gz', sigma.astype(np.float32), affine) - save_nifti(denoising_path + '/' + patient_path + '_mppca.nii.gz', denoised.astype(np.float32), affine) - else: - bvals_path = folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval' - bvals = np.loadtxt(bvals_path) - denoised = patch2self(b0_mask, bvals) - save_nifti(denoising_path + '/' + patient_path + '_patch2self.nii.gz', denoised.astype(np.float32), affine) + output, error = process.communicate() + denoised, _ = load_nifti(denoising_path + '/' + patient_path + '_mppca.nii.gz') - # computes and save the residuals - rms_diff = np.sqrt((data - denoised) ** 2) - save_nifti(denoising_path + '/' + patient_path + '_patch2self_residuals.nii.gz', rms_diff.astype(np.float32), affine) f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Denoising finished for patient %s \n" % p) print("[" + log_prefix + "] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Denoising finished for patient %s \n" % p) @@ -227,13 +190,11 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin print("[" + log_prefix + "] " + datetime.datetime.now().strftime( "%d.%b %Y %H:%M:%S") + ": End of denoising for patient %s \n" % p) - b0_mask = denoised + curr_dmri = denoised if not eddy and not gibbs and not topup and not biasfield: - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', - b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz', + curr_dmri.astype(np.float32), affine) shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval") shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bvec", @@ -245,23 +206,20 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin makedir(gibbs_path, folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", log_prefix) if (starting_state == "gibbs"): - mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz' if not denoising: - b0_mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz' + curr_dmri_path = nifti_path else: - b0_mask_path = denoising_path + '/' + patient_path + denoising_ext - b0_mask, affine, voxel_size = load_nifti(b0_mask_path, return_voxsize=True) - mask, _ = load_nifti(mask_path) + curr_dmri_path = denoising_path + '/' + patient_path + denoising_ext + curr_dmri, affine, voxel_size = load_nifti(curr_dmri_path, return_voxsize=True) - data = gibbs_removal(b0_mask, num_processes=core_count) + data = gibbs_removal(curr_dmri, num_processes=core_count) corrected_path = folder_path + '/subjects/' + patient_path + "/dMRI/preproc/gibbs/" + patient_path + '_gibbscorrected.nii.gz' save_nifti(corrected_path, data.astype(np.float32), affine) + curr_dmri = data if not eddy and not topup and not biasfield: - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', - data.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz', + curr_dmri.astype(np.float32), affine) shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval") shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bvec", @@ -270,7 +228,6 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin # Explicitly freeing memory import gc denoised = None - b0_mask = None mask = None data = None affine = None @@ -293,7 +250,7 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin elif denoising: imain_tot = denoising_path + '/' + patient_path + denoising_ext else: - imain_tot = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz' + imain_tot = nifti_path multiple_encoding=False topup_log = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/topup/topup_logs.txt", "a+") @@ -371,21 +328,20 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin # wait until topup finish output, error = process.communicate() - if not eddy: - inindex="" - first=True - for r in roi: - if first: - inindex = str(topup_index[r-1]) - else: - inindex = inindex + "," + str(topup_index[r-1]) + inindex="" + first=True + for r in roi: + if first: + inindex = str(topup_index[r-1]) + else: + inindex = inindex + "," + str(topup_index[r-1]) - bashCommand2 = 'applytopup --imain="' + imain_tot + '" --inindex='+inindex+' --datain="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --topup="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_estimate" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr"' + bashCommand2 = 'applytopup --imain="' + imain_tot + '" --inindex='+inindex+' --datain="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --topup="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_estimate" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr"' - process2 = subprocess.Popen(bashCommand2, universal_newlines=True, shell=True, stdout=topup_log, - stderr=subprocess.STDOUT) - # wait until apply topup finish - output, error = process2.communicate() + process2 = subprocess.Popen(bashCommand2, universal_newlines=True, shell=True, stdout=topup_log, + stderr=subprocess.STDOUT) + # wait until apply topup finish + output, error = process2.communicate() else: f = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", "a+") @@ -414,13 +370,12 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin synb0DisCo_starting_step = "topup" synb0DisCo(folder_path,topup_path,patient_path,starting_step=synb0DisCo_starting_step,topup=True,gpu=useGPUsynb0DisCo, static_files_path=static_files_path) - if not eddy: - bashCommand2 = 'export OMP_NUM_THREADS='+str(core_count)+' ; export FSLPARALLEL='+str(core_count)+' ; applytopup --imain="' + imain_tot + '" --inindex=1 --datain="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --topup="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_estimate" --method=jac --interp=spline --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr"' + bashCommand2 = 'export OMP_NUM_THREADS='+str(core_count)+' ; export FSLPARALLEL='+str(core_count)+' ; applytopup --imain="' + imain_tot + '" --inindex=1 --datain="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --topup="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_estimate" --method=jac --interp=spline --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr"' - process2 = subprocess.Popen(bashCommand2, universal_newlines=True, shell=True, stdout=topup_log, - stderr=subprocess.STDOUT) - # wait until apply topup finish - output, error = process2.communicate() + process2 = subprocess.Popen(bashCommand2, universal_newlines=True, shell=True, stdout=topup_log, + stderr=subprocess.STDOUT) + # wait until apply topup finish + output, error = process2.communicate() topup_log.close() @@ -432,28 +387,90 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin "%d.%b %Y %H:%M:%S") + ": End of topup for patient %s \n" % p) f.close() + + ## Compute pre eddy/biasfield mask + gc.collect() + from elikopy.utils import clean_mask + + topup_corr_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr.nii.gz' + + topup_corr, affine = load_nifti(topup_corr_path) + topup_corr_b0_ref = topup_corr[..., 0] + dwiref_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_corr_dwiref.nii.gz' + save_nifti(dwiref_path, topup_corr_b0_ref.astype(np.float32), affine) + topup_corr_b0_ref = None + topup_corr = None + gc.collect() + + + # Step 1 : dwi2mask + bvec_path = folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec' + bval_path = folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval' + dwi2mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_space-dwi_type-dwi2mask_brainmask.nii.gz' + cmd = f"dwi2mask -fslgrad {bvec_path} {bval_path} {topup_corr_path} {dwi2mask_path} -force " + import subprocess + print("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": dwi2mask launched for patient %s \n" % p + " with bash command " + cmd) + f = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", "a+") + f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": dwi2mask launched for patient %s \n" % p + " with bash command " + cmd) + f.flush() + process = subprocess.Popen(cmd, universal_newlines=True, shell=True, stdout=f, stderr=subprocess.STDOUT) + # wait until dwi2mask finish + output, error = process.communicate() + f.flush() + f.close() + + gc.collect() + # Step 2 : mri_synth_strip on dwiref + mrisynthstrip_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_space-dwi_type-mrisynthstrip_brainmask.nii.gz' + cmd = f"mri_synth_strip -i {dwiref_path} -m {mrisynthstrip_path}" + print("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": mri_synth_strip launched for patient %s \n" % p + " with bash command " + cmd) + f = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", "a+") + f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": mri_synth_strip launched for patient %s \n" % p + " with bash command " + cmd) + f.flush() + process = subprocess.Popen(cmd, universal_newlines=True, shell=True, stdout=f, stderr=subprocess.STDOUT) + output, error = process.communicate() + f.flush() + f.close() + + # Step 3 : median otsu on preprocess data + topup_corr, affine = load_nifti(topup_corr_path) + _, mask = median_otsu(topup_corr, median_radius=2, numpass=1, vol_idx=range(0, np.shape(topup_corr)[3]), + dilate=2) + mask = clean_mask(mask) + save_nifti( + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_type-otsu_dilate-2_brainmask.nii.gz', + mask.astype(np.float32), affine) + + # Step 4: Apply all masks to preprocess data + dwi2mask_mask, _ = load_nifti(dwi2mask_path) + mrisynthstrip_mask, _ = load_nifti(mrisynthstrip_path) + + full_mask = np.logical_or(mask, np.logical_or(dwi2mask_mask, mrisynthstrip_mask)) + full_mask = binary_dilation(full_mask, iterations=1) + + topup_corr_full_mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_brain_mask.nii.gz' + save_nifti(topup_corr_full_mask_path, + full_mask.astype(np.float32), affine) + if not eddy and not biasfield: data, affine = load_nifti( folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + "_unwarped.nii.gz") - b0_mask, mask = median_otsu(data, median_radius=2, numpass=1, vol_idx=range(0, np.shape(data)[3]), dilate=2) - mask = clean_mask(mask) - # Apply cleaned mask to b0_mask - b0_mask = b0_mask * mask[..., None] - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', - b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) - - _, mask_nodilate = median_otsu(data, median_radius=2, numpass=1, - vol_idx=range(0, np.shape(data)[3]), dilate=None) - mask_nodilate = clean_mask(mask_nodilate) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', - mask_nodilate.astype(np.float32), affine) - + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz', data.astype(np.float32), affine) shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval") shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bvec", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bvec") + save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', full_mask.astype(np.float32), affine) + + if topup: + topup_corr_full_mask_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_brain_mask.nii.gz' + processing_mask = topup_corr_full_mask_path + else: + processing_mask = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz' if eddy and starting_state!="biasfield" and starting_state!="report": print("[" + log_prefix + "] " + datetime.datetime.now().strftime( @@ -475,32 +492,32 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin slspec_path = folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'slspec.txt' if slspec_gc_path is not None and os.path.isdir(slspec_gc_path): if gibbs: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' elif denoising: - bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + nifti_path + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --ge_slspecs="' + slspec_gc_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' elif os.path.isfile(slspec_path): if gibbs: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' elif denoising: - bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + nifti_path + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --slspec="' + slspec_path + '" --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: if gibbs: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' elif denoising: - bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + nifti_path + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --mporder=' + str(s2v[0]) + ' --s2v_niter=' + str(s2v[1]) + ' --s2v_lambda=' + str(s2v[2]) + ' --s2v_interp=' + s2v[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: if gibbs: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/gibbs/' + patient_path + '_gibbscorrected.nii.gz" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' elif denoising: - bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + denoising_path + '/' + patient_path + denoising_ext + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' else: - bashCommand = eddycmd + ' --imain="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz" --mask="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_binary_mask.nii.gz" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' + bashCommand = eddycmd + ' --imain="' + nifti_path + '" --mask="' + processing_mask + '" --acqp="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'acqparams.txt" --index="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + 'index.txt" --bvecs="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bvec" --bvals="' + folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + '_raw_dmri.bval" --out="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + '_eddy_corr" --verbose --cnr_maps --residuals --repol=' + str(olrep[0]) + ' --ol_nstd=' + str(olrep[1]) + ' --ol_nvox=' + str(olrep[2]) + ' --ol_type=' + olrep[3] + ' --fwhm=' + fwhm + ' --niter=' + str(niter) + ' --slm=linear' if topup: bashCommand = bashCommand + ' --topup="' + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/topup/' + patient_path + '_topup_estimate"' @@ -530,24 +547,12 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin "%d.%b %Y %H:%M:%S") + ": End of eddy for patient %s \n" % p) f.close() - data, affine = load_nifti( - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + "_eddy_corr.nii.gz") - b0_mask, mask = median_otsu(data, median_radius=2, numpass=1, vol_idx=range(0, np.shape(data)[3]), dilate=2) - mask = clean_mask(mask) - # Apply cleaned mask to b0_mask - b0_mask = b0_mask * mask[..., None] - - _, mask_nodilate = median_otsu(data, median_radius=2, numpass=1, - vol_idx=range(0, np.shape(data)[3]), dilate=None) - mask_nodilate = clean_mask(mask_nodilate) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', - mask_nodilate.astype(np.float32), affine) - if not biasfield: - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', - b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) + data, affine = load_nifti( + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/eddy/' + patient_path + "_eddy_corr.nii.gz") + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz', + data.astype(np.float32), affine) + shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval") shutil.copyfile( @@ -569,33 +574,16 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin elif denoising: inputImage = denoising_path + '/' + patient_path + denoising_ext else: - inputImage = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/bet/' + patient_path + '_mask.nii.gz' - - #bashCommand= "N4BiasFieldCorrection -i " + inputImage + " -o [" + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_corr.nii.gz, " + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_est.nii.gz] -d 4" - + inputImage = nifti_path - bashCommand = 'export OMP_NUM_THREADS='+str(core_count)+' ; dwibiascorrect ants {} {} -fslgrad {} {} -mask {} -bias {} -scratch {} -force -info -nthreads {}'.format( + bashCommand = 'export OMP_NUM_THREADS='+str(core_count)+' ; dwibiascorrect ants {} {} -fslgrad {} {} -bias {} -scratch {} -force -info -nthreads {}'.format( inputImage, folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_corr.nii.gz", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bvec", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval", - folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_est.nii.gz", folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/tmp', core_count) - '''bashCommand = 'export OMP_NUM_THREADS=' + str( - core_count) + ' ; dwibiascorrect ants {} {} -fslgrad {} {} -mask {} -bias {} -scratch {} -force -info -nthreads {} -ants.b {} -ants.c {} '.format( - inputImage, - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_corr.nii.gz", - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bvec", - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval", - folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + "_biasfield_est.nii.gz", - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/tmp', - core_count, - biasfield_bsplineFitting, - biasfield_convergence)''' - import subprocess bashcmd = bashCommand.split() print("[" + log_prefix + "] " + datetime.datetime.now().strftime( @@ -613,27 +601,9 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin output, error = process.communicate() biasfield_log.close() - #shutil.rmtree(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/tmp', ignore_errors=True) - shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + '_biasfield_corr.nii.gz', - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz') + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz') - data, affine = load_nifti( - folder_path + '/subjects/' + patient_path + '/dMRI/preproc/biasfield/' + patient_path + '_biasfield_corr.nii.gz') - b0_mask, mask = median_otsu(data, median_radius=2, numpass=1, vol_idx=range(0, np.shape(data)[3]), dilate=2) - mask = clean_mask(mask) - # Apply cleaned mask to b0_mask - b0_mask = b0_mask * mask[..., None] - save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', - b0_mask.astype(np.float32), affine) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', - mask.astype(np.float32), affine) - - _, mask_nodilate = median_otsu(data, median_radius=2, numpass=1, - vol_idx=range(0, np.shape(data)[3]), dilate=None) - mask_nodilate = clean_mask(mask_nodilate) - save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', - mask_nodilate.astype(np.float32), affine) if not eddy: shutil.copyfile(folder_path + '/subjects/' + patient_path + '/dMRI/raw/' + patient_path + "_raw_dmri.bval", @@ -657,6 +627,83 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin f.write(output) f.truncate() + # Generate b0 ref: + preproc, affine = load_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz') + b0_ref = preproc[..., 0] + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dwiref.nii.gz', b0_ref.astype(np.float32), affine) + + preproc = None + b0_ref = None + gc.collect() + #### Generate final mask #### + from elikopy.utils import clean_mask + + # Step 1 : dwi2mask + preproc_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz' + bvec_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.bvec' + bval_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.bval' + dwi2mask_path = folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_space-dwi_type-dwi2mask_brainmask.nii.gz' + cmd = f"dwi2mask -fslgrad {bvec_path} {bval_path} {preproc_path} {dwi2mask_path} -force" + import subprocess + print("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": dwi2mask launched for patient %s \n" % p + " with bash command " + cmd) + f = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", "a+") + f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": dwi2mask launched for patient %s \n" % p + " with bash command " + cmd) + f.flush() + process = subprocess.Popen(cmd, universal_newlines=True, shell=True, stdout=f, stderr=subprocess.STDOUT) + # wait until dwi2mask finish + output, error = process.communicate() + f.flush() + f.close() + + # Step 2 : mri_synth_strip on dwiref + dwiref_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dwiref.nii.gz' + mrisynthstrip_path = folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_space-dwi_type-mrisynthstrip_brainmask.nii.gz' + cmd = f"mri_synth_strip -i {dwiref_path} -m {mrisynthstrip_path}" + print("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": mri_synth_strip launched for patient %s \n" % p + " with bash command " + cmd) + f = open(folder_path + '/subjects/' + patient_path + "/dMRI/preproc/preproc_logs.txt", "a+") + f.write("[" + log_prefix + "] " + datetime.datetime.now().strftime( + "%d.%b %Y %H:%M:%S") + ": mri_synth_strip launched for patient %s \n" % p + " with bash command " + cmd) + f.flush() + process = subprocess.Popen(cmd, universal_newlines=True, shell=True, stdout=f, stderr=subprocess.STDOUT) + output, error = process.communicate() + f.flush() + f.close() + + # Step 3 : median otsu on preprocess data + preproc, affine = load_nifti( + folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc_nomask.nii.gz') + preproc_masked, mask = median_otsu(preproc, median_radius=2, numpass=1, vol_idx=range(0, np.shape(preproc)[3]), dilate=2) + mask = clean_mask(mask) + save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_type-otsu_dilate-2_brainmask.nii.gz', + mask.astype(np.float32), affine) + + # Step 4: Apply all masks to preprocess data + dwi2mask_mask, _ = load_nifti(dwi2mask_path) + mrisynthstrip_mask, _ = load_nifti(mrisynthstrip_path) + + all_mask = mask + dwi2mask_mask + mrisynthstrip_mask + full_mask = np.zeros_like(all_mask) + full_mask[all_mask >= 2] = 1 + full_mask = clean_mask(full_mask) + + full_mask_inclusive = np.logical_or(mask, np.logical_or(dwi2mask_mask, mrisynthstrip_mask)) + full_mask_inclusive = clean_mask(full_mask_inclusive) + # dilate mask + full_mask_inclusive = binary_dilation(full_mask_inclusive, iterations=1) + preproc_masked = preproc * full_mask_inclusive[..., np.newaxis] + + save_nifti(folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz', + preproc_masked.astype(np.float32), affine) + save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask.nii.gz', + full_mask.astype(np.float32), affine) + save_nifti(folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + '_brain_mask_dilated.nii.gz', + full_mask_inclusive.astype(np.float32), affine) + + + if not report: return @@ -680,7 +727,6 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin from dipy.align.transforms import RigidTransform3D from dipy.segment.mask import segment_from_cfa from dipy.segment.mask import bounding_box - from scipy.ndimage.morphology import binary_dilation from os.path import isdir from skimage import measure from fpdf import FPDF @@ -708,8 +754,13 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin reslice_data, reslice_affine = load_nifti(preproc_path + "reslice/" + patient_path + "_reslice.nii.gz") # bet data (stage to compare with final) - bet_data, bet_affine = load_nifti(preproc_path + "bet/" + patient_path + "_mask.nii.gz") mask_raw, mask_raw_affine = load_nifti(preproc_path + "bet/" + patient_path + "_binary_mask.nii.gz") + if bool_reslice: + bet_data = reslice_data * mask_raw[..., np.newaxis] + bet_affine = reslice_affine + else: + bet_data = raw_data * mask_raw[..., np.newaxis] + bet_affine = raw_affine # mppca data bool_mppca = isdir(os.path.join(preproc_path, "mppca")) diff --git a/elikopy/registration.py b/elikopy/registration.py index 72b5b04..a63b528 100644 --- a/elikopy/registration.py +++ b/elikopy/registration.py @@ -22,7 +22,7 @@ def getTransform(static_volume_file, moving_volume_file, mask_file=None, onlyAffine=False, - diffeomorph=True, sanity_check=False, DWI=False): + diffeomorph=True, sanity_check=False, DWI=False, affine_map=None): ''' @@ -56,8 +56,9 @@ def getTransform(static_volume_file, moving_volume_file, mask_file=None, onlyAff if sanity_check or onlyAffine: - identity = np.eye(4) - affine_map = AffineMap(identity, + if affine_map is None: + affine_map = np.eye(4) + affine_map = AffineMap(affine_map, static.shape, static_grid2world, moving.shape, moving_grid2world) @@ -243,6 +244,84 @@ def applyTransformToAllMapsInFolder(input_folder, output_folder, mapping, mappin continue + +def regToT1fromB0FSL(reg_path, T1_subject, DWI_subject, mask_file, metrics_dic, folderpath, p, mapping_T1_to_T1MNI, T1_MNI, + mask_static, FA_MNI, longitudinal_transform=None): + if os.path.exists(reg_path + 'mapping_DWI_B0FSL_to_T1.p'): + with open(reg_path + 'mapping_DWI_B0FSL_to_T1.p', 'rb') as handle: + mapping_DWI_to_T1 = pickle.load(handle) + else: + if not (os.path.exists(reg_path)): + try: + os.makedirs(reg_path) + except OSError: + print("Creation of the directory %s failed" % reg_path) + + DWI_B0_subject = os.path.join(folderpath, 'subjects', p, 'dMRI', 'preproc', p + '_dwiref.nii.gz') + T1_subject_raw = os.path.join(folderpath, 'subjects', p, 'T1', p + '_T1.nii.gz') + T1_brain_subject = os.path.join(folderpath, 'subjects', p, 'T1', p + '_T1_brain.nii.gz') + b0fsl_reg_path = os.path.join(folderpath, 'subjects', p, 'reg', 'B0FSL_to_T1') + os.makedirs(b0fsl_reg_path, exist_ok=True) + cmd = f"epi_reg --epi={DWI_B0_subject} --t1={T1_subject_raw} --t1brain={T1_brain_subject} --out={b0fsl_reg_path}/{p}_B0toT1 " + cmd_part2 = f"c3d_affine_tool -ref {T1_subject_raw} -src {DWI_B0_subject} {b0fsl_reg_path}/{p}_B0toT1.mat -fsl2ras -oitk {b0fsl_reg_path}/{p}_B0toT1_ANTS.txt -o {b0fsl_reg_path}/{p}_B0toT1_ANTS.mat" + cmd += " && " + cmd_part2 + import subprocess + process = subprocess.Popen(cmd, shell=True, universal_newlines=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + out, err = process.communicate() + + affine_map_fslb0 = np.loadtxt(f"{b0fsl_reg_path}/{p}_B0toT1_ANTS.mat") + mapping_DWI_to_T1 = getTransform(T1_subject, DWI_subject, mask_file=mask_file, onlyAffine=True, diffeomorph=False, sanity_check=False, DWI=True, affine_map=affine_map_fslb0) + with open(reg_path + 'mapping_DWI_B0FSL_to_T1.p', 'wb') as handle: + pickle.dump(mapping_DWI_to_T1, handle, protocol=pickle.HIGHEST_PROTOCOL) + + if not (os.path.exists(folderpath + "/subjects/" + p + "/masks/reg/")): + try: + os.makedirs(folderpath + "/subjects/" + p + "/masks/reg/") + except OSError: + print("Creation of the directory %s failed" % folderpath + "/subjects/" + p + "/masks/reg/") + + for maskType in ["brain_mask_dilated","brain_mask", "wm_mask_MSMT", "wm_mask_AP", "wm_mask_FSL_T1", + "wm_mask_Freesurfer_T1"]: + in_mask_path = folderpath + "/subjects/" + p + "/masks/" + p + "_" + maskType + ".nii.gz" + if os.path.exists(in_mask_path): + if longitudinal_transform is not None: + reg_mask_path = folderpath + "/subjects/" + p + "/masks/reg/" + p + "_B0FSL_" + maskType + "_longitudinal.nii.gz" + applyTransform(in_mask_path, mapping_DWI_to_T1, mapping_2=longitudinal_transform, mapping_3=mapping_T1_to_T1MNI, static_file=T1_MNI, + output_path=reg_mask_path, mask_file=None, binary=False, inverse=False, + mask_static=mask_static, static_fa_file=FA_MNI) + else: + reg_mask_path = folderpath + "/subjects/" + p + "/masks/reg/" + p + "_B0FSL_" + maskType + ".nii.gz" + applyTransform(in_mask_path, mapping_DWI_to_T1, mapping_2=mapping_T1_to_T1MNI, static_file=T1_MNI, + output_path=reg_mask_path, mask_file=None, binary=False, inverse=False, + mask_static=mask_static, static_fa_file=FA_MNI) + + for key, value in metrics_dic.items(): + input_folder = folderpath + '/subjects/' + p + '/dMRI/microstructure/' + value + '/' + if longitudinal_transform is not None: + output_folder = folderpath + '/subjects/' + p + '/dMRI/microstructure/' + value + '_CommonSpace_T1_B0FSL_longitudinal/' + else: + output_folder = folderpath + '/subjects/' + p + '/dMRI/microstructure/' + value + '_CommonSpace_T1_B0FSL/' + + if not (os.path.exists(output_folder)): + try: + os.makedirs(output_folder) + except OSError: + print("Creation of the directory %s failed" % output_folder) + + print("Start of applyTransformToAllMapsInFolder for metrics ", value, ":", key) + if longitudinal_transform is not None: + applyTransformToAllMapsInFolder(input_folder, output_folder, mapping_DWI_to_T1, mapping_2=longitudinal_transform, + mapping_3=mapping_T1_to_T1MNI, static_file=T1_MNI, mask_file=mask_file, + keywordList=[p, key], inverse=False, mask_static=mask_static, + static_fa_file=FA_MNI) + else: + applyTransformToAllMapsInFolder(input_folder, output_folder, mapping_DWI_to_T1, mapping_2=mapping_T1_to_T1MNI, + static_file=T1_MNI, mask_file=mask_file, keywordList=[p, key], inverse=False, + mask_static=mask_static, static_fa_file=FA_MNI) + + + + def regToT1fromB0(reg_path, T1_subject, DWI_subject, mask_file, metrics_dic, folderpath, p, mapping_T1_to_T1MNI, T1_MNI, mask_static, FA_MNI, longitudinal_transform=None): if os.path.exists(reg_path + 'mapping_DWI_B0_to_T1.p'): @@ -432,7 +511,7 @@ def regToT1fromAP(reg_path, T1_subject, AP_subject, mask_file, metrics_dic, fold mask_static=mask_static, static_fa_file=FA_MNI) -def regallDWIToT1wToT1wCommonSpace(folder_path, p, DWI_type="AP", maskType=None, T1_filepath=None, T1wCommonSpace_filepath="${FSLDIR}/data/standard/MNI152_T1_1mm_brain.nii.gz", T1wCommonSpaceMask_filepath="${FSLDIR}/data/standard/MNI152_T1_1mm_brain_mask.nii.gz", metrics_dic={'_FA': 'dti', 'RD': 'dti', 'AD': 'dti', 'MD': 'dti'}, longitudinal=False): +def regallDWIToT1wToT1wCommonSpace(folder_path, p, DWI_type="B0FSL", maskType="brain_mask", T1_filepath=None, T1wCommonSpace_filepath="${FSLDIR}/data/standard/MNI152_T1_1mm_brain.nii.gz", T1wCommonSpaceMask_filepath="${FSLDIR}/data/standard/MNI152_T1_1mm_brain_mask.nii.gz", metrics_dic={'_FA': 'dti', 'RD': 'dti', 'AD': 'dti', 'MD': 'dti'}, longitudinal=False): preproc_folder = folder_path + '/subjects/' + p + '/dMRI/preproc/' T1_CommonSpace = os.path.expandvars(T1wCommonSpace_filepath) FA_MNI = os.path.expandvars('${FSLDIR}/data/standard/FSL_HCP1065_FA_1mm.nii.gz') @@ -440,7 +519,7 @@ def regallDWIToT1wToT1wCommonSpace(folder_path, p, DWI_type="AP", maskType=None, assert maskType in ["brain_mask_dilated","brain_mask", "wm_mask_MSMT", "wm_mask_AP", "wm_mask_FSL_T1", "wm_mask_Freesurfer_T1", None], "The mask parameter must be one of the following : brain_mask_dilated, brain_mask, wm_mask_MSMT, wm_mask_AP, wm_mask_FSL_T1, wm_mask_Freesurfer_T1, None" - assert DWI_type in ["AP", "WMFOD", "BO"], "The DWI_type parameter must be one of the following : AP, WMFOD, BO" + assert DWI_type in ["AP", "WMFOD", "B0", "B0FSL"], "The DWI_type parameter must be one of the following : AP, WMFOD, B0, B0FSL" mask_path = "" if maskType is not None and os.path.isfile(folder_path + '/subjects/' + p + "/masks/" + p + '_' + maskType + '.nii.gz'): @@ -546,6 +625,11 @@ def regallDWIToT1wToT1wCommonSpace(folder_path, p, DWI_type="AP", maskType=None, regToT1fromWMFOD(reg_path, T1_subject, WM_FOD_subject, mask_path, metrics_dic, folder_path, p, mapping_T1w_to_T1wCommonSpace, T1_CommonSpace, mask_static, FA_MNI, longitudinal_transform=mapping_T1w_to_T1wRef) elif DWI_type == "AP": regToT1fromAP(reg_path, T1_subject, AP_subject, mask_path, metrics_dic, folder_path, p, mapping_T1w_to_T1wCommonSpace, T1_CommonSpace, mask_static, FA_MNI, longitudinal_transform=mapping_T1w_to_T1wRef) + elif DWI_type == "B0FSL": + regToT1fromB0FSL(reg_path, T1_subject, AP_subject, mask_path, metrics_dic, folder_path, p, + mapping_T1w_to_T1wCommonSpace, T1_CommonSpace, mask_static, FA_MNI, + longitudinal_transform=mapping_T1w_to_T1wRef) + else: print("DWI_type not recognized") diff --git a/elikopy/utils.py b/elikopy/utils.py index e584e32..04a2a43 100644 --- a/elikopy/utils.py +++ b/elikopy/utils.py @@ -47,6 +47,7 @@ def submit_job(job_info): slurm_cmd.append("--%s=%s" % (key, value)) if script: slurm_cmd.append(job_info["script"]) + slurm_cmd.append("--hint=multithread") print("[INFO] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Generated slurm batch command: '%s'" % slurm_cmd) @@ -678,8 +679,7 @@ def synb0DisCo(folder_path, topuppath, patient_path, static_files_path=None, sta # Skull strip T1 - bet = "bet " + synb0path + "/T1.nii.gz " + synb0path + \ - "/T1_mask.nii.gz -R" # + " -f 0.4 -g -0.2" + bet = shutil.copyfile(os.path.join(folder_path,"subjects",patient_path,"T1",f"{patient_path}_T1_brain.nii.gz"), synb0path + "/T1_mask.nii.gz") # epi_reg distorted b0 to T1; wont be perfect since B0 is distorted @@ -692,7 +692,7 @@ def synb0DisCo(folder_path, topuppath, patient_path, static_files_path=None, sta "/b0.nii.gz " + synb0path + "/epi_reg_d.mat -fsl2ras -oitk " + \ synb0path + "/epi_reg_d_ANTS.txt" - # ANTs register T1 to atla + # ANTs register T1 to atlas antsRegistrationSyNQuick = "antsRegistrationSyNQuick.sh -d 3 -f " + static_files_path + \ "/atlases/mni_icbm152_t1_tal_nlin_asym_09c.nii.gz -m " + \ synb0path + "/T1.nii.gz -o " + synb0path + "/ANTS" @@ -1695,7 +1695,7 @@ def clean_mask(mask): center = tuple([np.average(indices) for indices in np.where(mask == 1)]) center = tuple([int(point) for point in center]) - mask = flood(mask, center) + mask = flood(mask, center, connectivity=1) mask_cleaned = np.zeros((mask.shape)) mask_cleaned[mask] = 1