Merge pull request #37 from Hyedryn/restricted_DTI

Limiting DTI to b=2000

elikopy/core.py

@@ -23,7 +23,8 @@
 
 
 def dicom_to_nifti(folder_path):
-    """ Convert dicom data into compressed nifti. Converted dicoms are then moved to a sub-folder named original_data.
+    """ Convert dicom data into compressed nifti. Converted dicoms are then
+    moved to a sub-folder named original_data.
     The niftis are named patientID_ProtocolName_SequenceName.
 
     :param folder_path: Path to root folder containing all the dicoms
@@ -653,14 +654,18 @@ 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") + ": All the preprocessing operation are finished!\n")
         f.close()
 
-    def dti(self,folder_path=None, patient_list_m=None, maskType="brain_mask_dilated", slurm=None, slurm_email=None, slurm_timeout=None, slurm_cpus=None, slurm_mem=None):
+    def dti(self,folder_path=None, patient_list_m=None, maskType="brain_mask_dilated", use_all_shells: bool = False,
+            slurm=None, slurm_email=None, slurm_timeout=None, slurm_cpus=None, slurm_mem=None):
         """Computes the DTI metrics for each subject using Weighted Least-Squares. The outputs are available in the directories <folder_path>/subjects/<subjects_ID>/dMRI/dti/.
 
         example : study.dti()
 
         :param folder_path: the path to the root directory. default=study_folder
         :param patient_list_m: Define a subset of subjects to process instead of all the available subjects. example : ['patientID1','patientID2','patientID3']. default=None
         :param maskType: Define which mask to use during processing. default="brain_mask_dilated"
+        :param use_all_shells: Boolean. DTI will use all shells available, not just
+        shells <= 2000, this will cause a more defined white matter at the cost of
+        an erronous estimation of the CSF. The default is False.
         :param slurm: Whether to use the Slurm Workload Manager or not (for computer clusters). default=value_during_init
         :param slurm_email: Email adress to send notification if a task fails. default=None
         :param slurm_timeout: Replace the default slurm timeout of 1h by a custom timeout.
@@ -697,7 +702,7 @@ def dti(self,folder_path=None, patient_list_m=None, maskType="brain_mask_dilated
 
             if slurm:
                 p_job = {
-                        "wrap": "python -c 'from elikopy.individual_subject_processing import dti_solo; dti_solo(\"" + folder_path + "/\",\"" + p + "\",maskType=\"" + str(maskType) + "\")'",
+                        "wrap": "python -c 'from elikopy.individual_subject_processing import dti_solo; dti_solo(\"" + folder_path + "/\",\"" + p + "\",maskType=\"" + str(maskType) + "\", use_all_shells=" + str(use_all_shells) +")'",
                         "job_name": "dti_" + p,
                         "ntasks": 1,
                         "cpus_per_task": 1,
@@ -719,7 +724,7 @@ def dti(self,folder_path=None, patient_list_m=None, maskType="brain_mask_dilated
                 f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Patient %s is ready to be processed\n" % p)
                 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:
-                dti_solo(folder_path + "/",p,maskType=maskType)
+                dti_solo(folder_path + "/",p,maskType=maskType,use_all_shells=use_all_shells)
                 matplotlib.pyplot.close(fig='all')
                 f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Successfully applied DTI on patient %s\n" % p)
                 f.flush()

elikopy/individual_subject_processing.py

@@ -629,7 +629,7 @@ def preproc_solo(folder_path, p, reslice=False, reslice_addSlice=False, denoisin
     matplotlib.use('Agg')
     import matplotlib.pyplot as plt
     import dipy.reconst.dti as dti
-    from dipy.align.imaffine import (AffineMap, MutualInformationMetric, AffineRegistration)
+    from dipy.align.imaffine import MutualInformationMetric, AffineRegistration
     from dipy.align.transforms import RigidTransform3D
     from dipy.segment.mask import segment_from_cfa
     from dipy.segment.mask import bounding_box
@@ -1385,13 +1385,19 @@ def print_page(self, images):
     f.close()
 
 
-def dti_solo(folder_path, p, maskType="brain_mask_dilated", report=True):
+def dti_solo(folder_path, p, maskType="brain_mask_dilated",
+             use_all_shells: bool = False, report=True):
     """
-    Computes the DTI metrics for a single subject. The outputs are available in the directories <folder_path>/subjects/<subjects_ID>/dMRI/dti/.
+    Computes the DTI metrics for a single subject. The outputs are available in
+    the directories <folder_path>/subjects/<subjects_ID>/dMRI/dti/.
 
     :param folder_path: the path to the root directory.
     :param p: The name of the patient.
-    :param use_wm_mask: If true a white matter mask is used. The white_matter() function needs to already be applied. default=False
+    :param use_wm_mask: If true a white matter mask is used. The white_matter()
+    function needs to already be applied. default=False
+    :param use_all_shells: Boolean. DTI will use all shells available, not just
+    shells <= 2000, this will cause a more defined white matter at the cost of
+    an erronous estimation of the CSF. The default is False.
     """
     log_prefix = "DTI SOLO"
     print("[" + log_prefix + "] " + datetime.datetime.now().strftime(
@@ -1424,10 +1430,19 @@ def dti_solo(folder_path, p, maskType="brain_mask_dilated", report=True):
     bvals, bvecs = read_bvals_bvecs(
         folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bval",
         folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + "_dmri_preproc.bvec")
+    # Remove shells >2000================================
+    if not use_all_shells:
+        indexes = np.argwhere(bvals < 2000+10)
+        indexes = indexes.squeeze()
+        bvals = bvals[indexes]
+        bvecs = bvecs[indexes]
+        data = data[..., indexes]
+        print('Warning: removing shells above b=2000 for DTI. To disable this, '
+              + 'activate the use_all_shells option.')
     # create the model===================================
     b0_threshold = np.min(bvals)+10
     b0_threshold = max(50, b0_threshold)
-    gtab = gradient_table(bvals, bvecs,b0_threshold=b0_threshold)
+    gtab = gradient_table(bvals, bvecs, b0_threshold=b0_threshold)
     tenmodel = dti.TensorModel(gtab)
     tenfit = tenmodel.fit(data, mask=mask)
     # FA ================================================