Merge pull request #28 from Hyedryn/sift

Sift

elikopy/core.py

@@ -14,8 +14,10 @@
 import matplotlib
 
 import elikopy.utils
-from elikopy.individual_subject_processing import preproc_solo, dti_solo, white_mask_solo, noddi_solo, diamond_solo, \
-    mf_solo, noddi_amico_solo, ivim_solo, odf_csd_solo, odf_msmtcsd_solo, tracking_solo, verdict_solo, clean_study_solo
+from elikopy.individual_subject_processing import (preproc_solo, dti_solo,
+    white_mask_solo, noddi_solo, diamond_solo, mf_solo, noddi_amico_solo,
+    ivim_solo, odf_csd_solo, odf_msmtcsd_solo, tracking_solo, sift_solo,
+    verdict_solo, clean_study_solo)
 
 from elikopy.utils import submit_job, get_job_state, makedir, tbss_utils, regall_FA, regall, randomise_all
 
@@ -1062,6 +1064,91 @@ def tracking(self, folder_path=None, streamline_number:int=100000, max_angle:int
         f=open(folder_path + "/logs.txt", "a+")
         f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": End of tractography\n")
         f.close()
+
+    def sift(self, folder_path=None, streamline_number: int = 100000,
+             msmtCSD: bool = True, input_filename: str = "tractogram",
+             slurm = None, patient_list_m = None, slurm_email = None,
+             slurm_timeout = None, cpus = None, slurm_mem = None):
+        """Computes the sifted odf using MSMT-CSD for each subject. The outputs are
+        available in the directories <folder_path>/subjects/<subjects_ID>/dMRI/tractography/.
+    
+        example : study.sift_solo()
+    
+        :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 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 20h by a custom timeout.
+        :param slurm_cpus: Replace the default number of slurm cpus of 1 by a custom number of cpus of using slum, or for standard processing, its the number of core available for processing.
+        :param slurm_mem: Replace the default amount of ram allocated to the slurm task (8096MO by cpu) by a custom amount of ram.
+        """
+        log_prefix = "SIFT"
+        folder_path = self._folder_path if folder_path is None else folder_path
+        slurm = self._slurm if slurm is None else slurm
+        slurm_email = self._slurm_email if slurm_email is None else slurm_email
+
+        f = open(folder_path + "/logs.txt", "a+")
+        f.write("["+log_prefix+"] " +
+                datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") +
+                ": Beginning of SIFT with slurm:" + str(slurm) + "\n")
+        f.close()
+
+        dest_success = folder_path + "/subjects/subj_list.json"
+        with open(dest_success, 'r') as f:
+            patient_list = json.load(f)
+
+        if patient_list_m:
+            patient_list = patient_list_m
+
+        core_count = 4 if cpus is None else cpus
+
+        job_list = []
+        f = open(folder_path + "/logs.txt", "a+")
+        for p in patient_list:
+            patient_path = p
+
+            tracking_path = (folder_path + '/subjects/' + patient_path
+                             + "/dMRI/tractography")
+            makedir(tracking_path, folder_path + '/subjects/' + patient_path
+                    + "/dMRI/tractography/sift_logs.txt", log_prefix)
+
+            if slurm:
+                p_job = {
+                        "wrap": "export MKL_NUM_THREADS="+ str(core_count)+" ; export OMP_NUM_THREADS="+ str(core_count)+" ; python -c 'from elikopy.individual_subject_processing import sift_solo; sift_solo(\"" + folder_path + "/\",\"" + p + "\", streamline_number =" + str(streamline_number) + ", msmtCSD=" + str(msmtCSD) + ", core_count=" + str(core_count) + ", input_filename= \"" + input_filename + "\"" + ")'",
+                        "job_name": "SIFT_" + p,
+                        "ntasks": 1,
+                        "cpus_per_task": core_count,
+                        "mem_per_cpu": 4096,
+                        "time": "00:25:00",
+                        "mail_user": slurm_email,
+                        "mail_type": "FAIL",
+                        "output": folder_path + '/subjects/' + patient_path + "/dMRI/tractography/slurm-%j.out",
+                        "error": folder_path + '/subjects/' + patient_path + "/dMRI/tractography/slurm-%j.err",
+                    }
+                p_job["time"] = p_job["time"] if slurm_timeout is None else slurm_timeout
+                p_job["mem_per_cpu"] = p_job["mem_per_cpu"] if slurm_mem is None else slurm_mem
+                p_job_id = {}
+                p_job_id["id"] = submit_job(p_job)
+                p_job_id["name"] = p
+                job_list.append(p_job_id)
+                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:
+                tracking_solo(folder_path + "/", p,  streamline_number=streamline_number, msmtCSD=msmtCSD, input_filename=input_filename, core_count=core_count)
+                matplotlib.pyplot.close(fig='all')
+                f.write("["+log_prefix+"] " + datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") + ": Successfully applied sift_solo on patient %s\n" % p)
+                f.flush()
+        f.close()
+
+        # Wait for all jobs to finish
+        if slurm:
+            elikopy.utils.getJobsState(folder_path, job_list, log_prefix)
+
+        f = open(folder_path + "/logs.txt", "a+")
+        f.write("["+log_prefix+"] " +
+                datetime.datetime.now().strftime("%d.%b %Y %H:%M:%S") +
+                ": End of SIFT\n")
+        f.close()
 
 
     def white_mask(self, maskType, folder_path=None, patient_list_m=None, corr_gibbs=True, debug=False, slurm=None, slurm_email=None, slurm_timeout=None, cpus=None, slurm_mem=None):

elikopy/individual_subject_processing.py

@@ -3533,7 +3533,57 @@ def tracking_solo(folder_path:str, p:str, streamline_number:int=100000,
 
     with open(output_file[:-3]+'json', 'w') as outfile:
         json.dump(params, outfile)
-
+
+
+def sift_solo(folder_path: str, p: str, streamline_number: int = 100000,
+              msmtCSD: bool = True, input_filename: str = 'tractogram',
+              core_count: int = 1):
+    """ Computes the sifted whole brain tractogram of a single patient based on
+    the fod obtained from msmt-CSD.
+
+    :param folder_path: the path to the root directory.
+    :param p: The name of the patient.
+    :param streamline_number: Number of streamlines in the final tractogram. default=100000
+    :param msmtCSD: boolean. If True then uses ODF from msmt-CSD, if False from CSD. default=True
+    :param input_filename: str. Specify input filename for tractogram.
+    """
+
+    from dipy.io.streamline import load_tractogram, save_trk
+
+    patient_path = p
+
+    if msmtCSD:
+        odf_file_path = (folder_path + '/subjects/' + patient_path
+                         + "/dMRI/ODF/MSMT-CSD/"+patient_path
+                         + "_MSMT-CSD_WM_ODF.nii.gz")
+    else:
+        odf_file_path = (folder_path + '/subjects/' + patient_path
+                         + "/dMRI/ODF/CSD/"+patient_path
+                         + "_CSD_SH_ODF_mrtrix.nii.gz")
+
+    tracking_path = folder_path + '/subjects/' + patient_path + "/dMRI/tractography/"
+    mask_path = folder_path + '/subjects/' + patient_path + '/masks/' + patient_path + "_brain_mask_dilated.nii.gz"
+    dwi_path = folder_path + '/subjects/' + patient_path + '/dMRI/preproc/' + patient_path + '_dmri_preproc.nii.gz'
+
+    input_file = tracking_path+patient_path+'_'+input_filename+'.tck'
+    output_file = tracking_path+patient_path+'_'+input_filename+'_sift.tck'
+
+    bashCommand=('tcksift ' + input_file + ' ' + odf_file_path + ' ' +
+                 output_file +
+                 ' -nthreads ' + str(core_count) +
+                 ' -term_number ' + str(streamline_number) +
+                 ' -force')
+
+    sift_log = open(tracking_path+"sift_logs.txt", "a+")
+    process = subprocess.Popen(bashCommand, universal_newlines=True, shell=True,
+                               stdout=sift_log, stderr=subprocess.STDOUT)
+
+    process.communicate()
+    sift_log.close()
+
+    tract = load_tractogram(output_file, dwi_path)
+    save_trk(tract, output_file[:-3]+'trk')
+
 
 def verdict_solo(folder_path, p, core_count=1, small_delta=0.003, big_delta=0.035, G1Ball_1_lambda_iso=0.9e-9, C1Stick_1_lambda_par=[3.05e-9, 10e-9],TumorCells_Dconst=0.9e-9):
     """ Computes the verdict metrics for a single. The outputs are available in the directories <folder_path>/subjects/<subjects_ID>/dMRI/microstructure/verdict/.