dwi preprocessing using metadata

Dockerfile

@@ -55,7 +55,7 @@ RUN apt-get update && \
     curl -sSL http://neuro.debian.net/lists/trusty.us-ca.full >> /etc/apt/sources.list.d/neurodebian.sources.list && \
     apt-key adv --recv-keys --keyserver hkp://pgp.mit.edu:80 0xA5D32F012649A5A9 && \
     apt-get update && \
-    apt-get install -y fsl-core=5.0.9-3~nd14.04+1 && \
+    apt-get install -y fsl-core=5.0.9-1~nd+1+nd16.04+1 && \
     rm -rf /var/lib/apt/lists/* /tmp/* /var/tmp/*
 
 # Configure environment
@@ -72,7 +72,7 @@ ENV FSLOUTPUTTYPE=NIFTI_GZ
 RUN echo "cHJpbnRmICJrcnp5c3p0b2YuZ29yZ29sZXdza2lAZ21haWwuY29tXG41MTcyXG4gKkN2dW12RVYzelRmZ1xuRlM1Si8yYzFhZ2c0RVxuIiA+IC9vcHQvZnJlZXN1cmZlci9saWNlbnNlLnR4dAo=" | base64 -d | sh
 
 # Install Connectome Workbench
-RUN apt-get update && apt-get -y install connectome-workbench=1.2.3-1~nd14.04+1
+RUN apt-get update && apt-get -y install connectome-workbench=1.2.3-1~nd16.04+1
 
 ENV CARET7DIR=/usr/bin
 

run.py

@@ -162,6 +162,7 @@ def run_diffusion_processsing(**args):
       '--echospacing="{echospacing}" '+ \
       '--PEdir={PEdir} ' + \
       '--gdcoeffs="NONE" ' + \
+      '--dwiname="{dwiname}" ' + \
       '--printcom=""'
     cmd = cmd.format(**args)
     run(cmd, cwd=args["path"], env={"OMP_NUM_THREADS": str(args["n_cpus"])})
@@ -343,14 +344,15 @@ def run_diffusion_processsing(**args):
                                                 type='bold',
                                                 extensions=["nii.gz", "nii"])]
         for fmritcs in bolds:
-            fmriname = "_".join(fmritcs.split("sub-")[-1].split("_")[1:]).split(".")[0]
+            fmriname = "_".join(fmritcs.split("sub-")[-1].split("_")[1:-1]).split(".")[0]
             assert fmriname
 
             fmriscout = fmritcs.replace("_bold", "_sbref")
             if not os.path.exists(fmriscout):
                 fmriscout = "NONE"
 
             fieldmap_set = layout.get_fieldmap(fmritcs)
+            print(fieldmap_set)
             if fieldmap_set and fieldmap_set["type"] == "epi":
                 SEPhaseNeg = None
                 SEPhasePos = None
@@ -409,31 +411,86 @@ def run_diffusion_processsing(**args):
                 if stage in args.stages:
                     stage_func()
 
-        dwis = layout.get(subject=subject_label, type='dwi',
-                                                 extensions=["nii.gz", "nii"])
-
-        # print(dwis)
-        # acqs = set(layout.get(target='acquisition', return_type='id',
-        #                       subject=subject_label, type='dwi',
-        #                       extensions=["nii.gz", "nii"]))
-        # print(acqs)
-        # posData = []
-        # negData = []
-        # for acq in acqs:
-        #     pos = "EMPTY"
-        #     neg = "EMPTY"
-        #     dwis = layout.get(subject=subject_label,
-        #                       type='dwi', acquisition=acq,
-        #                       extensions=["nii.gz", "nii"])
-        #     assert len(dwis) <= 2
-        #     for dwi in dwis:
-        #         dwi = dwi.filename
-        #         if "-" in layout.get_metadata(dwi)["PhaseEncodingDirection"]:
-        #             neg = dwi
-        #         else:
-        #             pos = dwi
-        #     posData.append(pos)
-        #     negData.append(neg)
-        #
-        # print(negData)
-        # print(posData)
+        posData = []
+        negData = []
+        PEdir = "None"
+        dwiname = "Diffusion"
+        dirnums = []
+        onerun = False
+
+        numruns = set(layout.get(target='run', return_type='id',
+                                 subject=subject_label, type='dwi',
+                                 extensions=["nii.gz", "nii"]))
+        # accounts for multiple runs, number of directions, and phase encoding directions
+        if not numruns:
+            onerun= True
+            numruns = {'run-01'}
+        if numruns:
+            for numrun in numruns:
+                if not onerun:
+                    bvals = [f.filename for f in layout.get(subject=subject_label,
+                                                            type='dwi', run=numrun,
+                                                            extensions=["bval"])]
+                else:
+                    bvals = [f.filename for f in layout.get(subject=subject_label,
+                                                            type='dwi', extensions=["bval"])]
+                ## find number of directions by reading bval files, then create dictionary with corresponding
+                # bval file name, number of directions, dwi image file name, and phase encoding direction (i or j).
+                dwi_dict = {'bvalFile':[], 'bval':[], 'dwiFile':[], 'direction':[]}
+                for bvalfile in bvals: # find number of directions
+                    with open(bvalfile) as f:
+                        bvalues = [bvalue for line in f for bvalue in line.split()]
+                    # fill in the rest of dictionary
+                    dwi_dict['bvalFile'].append(bvalfile)
+                    dwi_dict['bval'].append(len(bvalues) - 1)
+                    dwiFile = glob(os.path.join(os.path.dirname(bvalfile),'{0}.nii*'.format(os.path.basename(bvalfile).split('.')[0]))) # ensures bval file has same name as dwi file
+                    assert len(dwiFile) == 1
+                    dwi_dict['dwiFile'].append(dwiFile[0])
+                    dwi_dict['direction'].append(layout.get_metadata(dwiFile[0])["PhaseEncodingDirection"][0])
+
+                # check if length of lists in dictionary are the same
+                n = len(dwi_dict['bvalFile'])
+                assert all(len(dwi_dict[k]) == n for k,v in dwi_dict.items())
+
+                for dirnum in set(dwi_dict['bval']):
+                    ## the following statement extracts index values in dwi_dict['bval'] if the value matches
+                    # "dirnum", which is the number of directions (i.e. 98 or 99). These index values are used
+                    # to find the corresponding PE directions, dwi file names, etc. in the dictionary
+                    idxs = { i for k,v in dwi_dict.iteritems() for i in range(0,len(dwi_dict['bval'])) if v[i] == dirnum }
+                    PEdirNums = set([dwi_dict['direction'][i] for i in idxs])
+                    for PEdirNum in PEdirNums:
+                        dwis = [ dwi_dict['dwiFile'][i] for i in idxs if dwi_dict['direction'][i] == PEdirNum ]
+                        assert len(dwis) <= 2
+                        dwiname = "Diffusion" + "_dir-" + str(dirnum) + "_" + numrun + "_corr_" + str(PEdirNum)
+                        if "j" in PEdirNum:
+                            PEdir = 2
+                        elif "i" in PEdirNum:
+                            PEdir = 1
+                        else:
+                            RuntimeError("Phase encoding direction not specified for diffusion data.")
+                        pos = "EMPTY"
+                        neg = "EMPTY"
+                        gdcoeffs = "None"
+                        for dwi in dwis:
+                            if "-" in layout.get_metadata(dwi)["PhaseEncodingDirection"]:
+                                neg = dwi
+                            else:
+                                pos = dwi
+                        try:
+                            echospacing = layout.get_metadata(pos)["EffectiveEchoSpacing"] * 1000
+                            dwi_stage_dict = OrderedDict([("DiffusionPreprocessing", partial(run_diffusion_processsing,
+                                                                                             posData=pos,
+                                                                                             negData=neg,
+                                                                                             path=args.output_dir,
+                                                                                             subject="sub-%s" % subject_label,
+                                                                                             echospacing=echospacing,
+                                                                                             PEdir=PEdir,
+                                                                                             gdcoeffs="NONE",
+                                                                                             dwiname=dwiname,
+                                                                                             n_cpus=args.n_cpus))])
+                            for stage, stage_func in dwi_stage_dict.iteritems():
+                                if stage in args.stages:
+                                    stage_func()
+                        except NameError:
+                            print("You may have missing diffusion data in the positive phase encoding direction.")
+