dev branch for testing mrtrix bids inputs

designer2/designer.py

@@ -114,7 +114,7 @@ def execute(): #pylint: disable=unused-variable
         app.ARGS.input, app.ARGS.fslbval, app.ARGS.fslbvec, app.ARGS.bids)
 
     # ensure inputs are reasonable for subsequent dwi processing
-    assert_inputs(dwi_metadata, app.ARGS.pe_dir, app.ARGS.pf)
+    # assert_inputs(dwi_metadata, app.ARGS.pe_dir, app.ARGS.pf)
 
     # convert input data to .mif format and concatenate
     convert_input_data(dwi_metadata)
@@ -126,6 +126,8 @@ def execute(): #pylint: disable=unused-variable
                   index = shell_rows)
     print('input DWI data has properties:')
     print(shell_df)
+
+    import pdb; pdb.set_trace()
 
     app.goto_scratch_dir()
 

designer2/tmi.py

@@ -153,7 +153,7 @@ def execute(): #pylint: disable=unused-variable
         app.ARGS.input, app.ARGS.fslbval, app.ARGS.fslbvec, app.ARGS.bids)
     logger.info("Input information obtained.", extra={"dwi_metadata": dwi_metadata})
 
-    assert_inputs(dwi_metadata, None, None)
+    #assert_inputs(dwi_metadata, None, None)
     convert_input_data(dwi_metadata)
     logger.info("Input data converted successfully.")
 

lib/designer_input_utils.py

@@ -164,12 +164,18 @@ def convert_pe_dir_to_ijk(args_pe_dir):
     return pe_dir_app
 
 def assert_inputs(dwi_metadata, args_pe_dir, args_pf):
+    """
+    This function has been depreciated. Please use the convert_input_data function instead.
+    This function will be deleted after testing.
+    """
+
     import json
     from mrtrix3 import app, MRtrixError
     import numpy as np
 
     bidslist = dwi_metadata['bidslist']
-    try:
+    # if bids files are provided, check that the user inputs match the bids files
+    try:        
         bids = [json.load(open(i)) for i in bidslist]
 
         TE_bids = []
@@ -248,8 +254,13 @@ def assert_inputs(dwi_metadata, args_pe_dir, args_pf):
         pf = pf_bids
     else:
         pf = args_pf
-
-    def convert_to_float(frac_str):
+
+    dwi_metadata['pe_dir'] = pe_dir
+    dwi_metadata['pf'] = convert_to_float(pf)
+    dwi_metadata['TE'] = TE
+    dwi_metadata['bshape'] = bshape
+
+def convert_to_float(frac_str):
         try:
             return float(frac_str)
         except ValueError:
@@ -262,11 +273,6 @@ def convert_to_float(frac_str):
             frac = float(num) / float(denom)
             return whole - frac if whole < 0 else whole + frac
 
-    dwi_metadata['pe_dir'] = pe_dir
-    dwi_metadata['pf'] = convert_to_float(pf)
-    dwi_metadata['TE'] = TE
-    dwi_metadata['bshape'] = bshape
-
 def convert_input_data(dwi_metadata):
     """
     convert input data to .mif mrtrix format
@@ -276,6 +282,7 @@ def convert_input_data(dwi_metadata):
     from mrtrix3 import run, image, MRtrixError, app
     import numpy as np
     import os
+    import warnings
 
     miflist = []
     phaselist = []
@@ -284,28 +291,43 @@ def convert_input_data(dwi_metadata):
     bshapelist = []
     dwi_ind_size = [[0,0,0,0]]
 
+    te_per_series = []
+    pf_per_series = []
+    ped_per_series = []
+
     dwi_n_list = dwi_metadata['dwi_list']
     isdicom = dwi_metadata['isdicom']
     bveclist = dwi_metadata['bvecs']
     bvallist = dwi_metadata['bvals']
     dwi_ext = dwi_metadata['dwi_ext']
-    te_per_series = dwi_metadata['TE']
-    bshape_per_series = dwi_metadata['bshape']
     phase_n_list = dwi_metadata['phase_list']
     phase_ext = dwi_metadata['phase_ext']
     bshapelist_input = dwi_metadata['bshapelist']
     telist_input = dwi_metadata['telist']
+    bidslist = dwi_metadata['bidslist']
 
     if len(dwi_n_list) == 1:
         if not isdicom:
             if os.path.exists(bvallist[0]):
-                cmd = ('mrconvert -fslgrad %s %s %s%s %s/dwi.mif' % 
-                    (bveclist[0], bvallist[0], ''.join(dwi_n_list), ''.join(dwi_ext), app.SCRATCH_DIR))
-                run.command(cmd)
+                if os.path.exists(bidslist[0]):
+                    cmd = ('mrconvert -fslgrad %s %s -json_import %s %s%s %s/dwi.mif' % 
+                        (bveclist[0], bvallist[0], bidslist[0],''.join(dwi_n_list), ''.join(dwi_ext), app.SCRATCH_DIR))
+                    run.command(cmd)
+                else:
+                    print('... no bids files identified')
+                    cmd = ('mrconvert -fslgrad %s %s %s%s %s/dwi.mif' % 
+                        (bveclist[0], bvallist[0], ''.join(dwi_n_list), ''.join(dwi_ext), app.SCRATCH_DIR))
+                    run.command(cmd)
             elif ~os.path.exists(bvallist[0]) and ('.mif' in ''.join(dwi_ext)):
-                cmd = ('mrconvert %s%s %s/dwi.mif' % 
-                    (''.join(dwi_n_list), ''.join(dwi_ext), app.SCRATCH_DIR))
-                run.command(cmd)
+                if os.path.exists(bidslist[0]):
+                    cmd = ('mrconvert %s%s -json_import %s %s/dwi.mif' % 
+                        (''.join(dwi_n_list), ''.join(dwi_ext), bidslist[0], app.SCRATCH_DIR))
+                    run.command(cmd)
+                else:
+                    print('... no bids files identified')
+                    cmd = ('mrconvert %s%s %s/dwi.mif' % 
+                        (''.join(dwi_n_list), ''.join(dwi_ext), app.SCRATCH_DIR))
+                    run.command(cmd)
             else:
                 raise MRtrixError('please make sure that inputs are either .nii files accompanied by corresponding .bval and .bvec files, or a .mif file with embedded gradient information')
         else:
@@ -314,18 +336,77 @@ def convert_input_data(dwi_metadata):
             run.command(cmd)
         dwi_header = image.Header('%s/dwi.mif' % (app.SCRATCH_DIR))
         dwi_ind_size.append([ int(s) for s in dwi_header.size() ])
+
+        if app.ARGS.pf:
+            pf_per_series_app = convert_to_float(app.ARGS.pf)
+            try:
+                pf_per_series_bids = dwi_header.keyval()['PartialFourier']
+                if pf_per_series_app != pf_per_series_bids:
+                    warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value')
+            except:
+                pass
+            pf_per_series.append(pf_per_series_app)
+        else:
+            try:
+                pf_per_series.append(dwi_header.keyval()['PartialFourier'])
+            except:
+                warnings.warn('No partial fourier factor found in header, assuming full sampling')
+                pf_per_series.append(1)
+
+        if app.ARGS.pe_dir:
+            pe_dir_app = app.ARGS.pe_dir
+            pe_dir = convert_pe_dir_to_ijk(pe_dir_app)
+            try:
+                ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection']
+                if ped_per_series_bids != pe_dir_app:
+                    warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value')
+            except:
+                pass
+            ped_per_series.append(pe_dir)
+        else:
+            try:
+                ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection'])
+            except:
+                raise MRtrixError('No phase encoding direction found in header, please specify manually')
+
+        if app.ARGS.echo_time:
+            te_per_series.append(app.ARGS.echo_time)
+            try:
+                te_per_series_bids = dwi_header.keyval()['EchoTime']
+                if te_per_series_bids != app.ARGS.echo_time:
+                    warnings.warn('User defined echo time does not match that found in bids json. Using user defined value')
+            except:
+                pass
+        else:   
+            try:
+                te_per_series.append(dwi_header.keyval()['EchoTime'])
+            except:
+                warnings.warn('No echo time found in header, assuming 0')
+                te_per_series.append(0)
 
     else:
         for idx,i in enumerate(dwi_n_list):
             if not isdicom:
                 if os.path.exists(bvallist[idx]):
-                    cmd = ('mrconvert -fslgrad %s %s %s%s %s/dwi%s.mif' % 
-                        (bveclist[idx], bvallist[idx], i, dwi_ext[idx], app.SCRATCH_DIR, str(idx)))
-                    run.command(cmd)
+                    if os.path.exists(bidslist[idx]):
+                        cmd = ('mrconvert -fslgrad %s %s -json_import %s %s%s %s/dwi%s.mif' % 
+                            (bveclist[idx], bvallist[idx], bidslist[idx], i, dwi_ext[idx], app.SCRATCH_DIR, str(idx)))
+                        run.command(cmd)
+                    else:
+                        print('... no bids files identified for %s' % i)
+                        cmd = ('mrconvert -fslgrad %s %s %s%s %s/dwi%s.mif' % 
+                            (bveclist[idx], bvallist[idx], i, dwi_ext[idx], app.SCRATCH_DIR, str(idx)))
+                        run.command(cmd)
                 elif ~os.path.exists(bvallist[idx]) and ('.mif' in dwi_ext[idx]):
-                    cmd = ('mrconvert %s%s %s/dwi%s.mif' % 
-                        (i, dwi_ext[idx], app.SCRATCH_DIR, str(idx)))
-                    run.command(cmd)
+                    if os.path.exists(bidslist[idx]):
+                        cmd = ('mrconvert %s%s -json_import %s %s/dwi%s.mif' % 
+                            (i, dwi_ext[idx], bidslist[idx], app.SCRATCH_DIR, str(idx)))
+                        run.command(cmd)
+                    else:
+                        print('... no bids files identified for %s' % i)
+                        cmd = ('mrconvert %s%s %s/dwi%s.mif' % 
+                            (i, dwi_ext[idx], app.SCRATCH_DIR, str(idx)))
+                        run.command(cmd)
                 else:
                     raise MRtrixError('please make sure that inputs are either .nii files accompanied by corresponding .bval and .bvec files, or a .mif file with embedded gradient information')
             else:
@@ -335,6 +416,53 @@ def convert_input_data(dwi_metadata):
             dwi_header = image.Header('%s/dwi%s.mif' % (app.SCRATCH_DIR, str(idx)))
             dwi_ind_size.append([ int(s) for s in dwi_header.size() ])
             miflist.append('%s/dwi%s.mif' % (app.SCRATCH_DIR, str(idx)))
+
+            if app.ARGS.pf:
+                pf_per_series_app = convert_to_float(app.ARGS.pf)
+                try:
+                    pf_per_series_bids = dwi_header.keyval()['PartialFourier']
+                    if pf_per_series_app != pf_per_series_bids:
+                        warnings.warn('User defined partial fourier factor does not match that found in bids json. Using user defined value')
+                except:
+                    pass
+                pf_per_series.append(pf_per_series_app)
+            else:
+                try:
+                    pf_per_series.append(dwi_header.keyval()['PartialFourier'])
+                except:
+                    warnings.warn('No partial fourier factor found in header, assuming full sampling')
+                    pf_per_series.append(1)
+
+            if app.ARGS.pe_dir:
+                pe_dir_app = app.ARGS.pe_dir
+                pe_dir = convert_pe_dir_to_ijk(pe_dir_app)
+                try:
+                    ped_per_series_bids = dwi_header.keyval()['PhaseEncodingDirection']
+                    if ped_per_series_bids != pe_dir_app:
+                        warnings.warn('User defined phase encoding direction does not match that found in bids json. Using user defined value')
+                except:
+                    pass
+                ped_per_series.append(pe_dir)
+            else:
+                try:
+                    ped_per_series.append(dwi_header.keyval()['PhaseEncodingDirection'])
+                except:
+                    raise MRtrixError('No phase encoding direction found in header, please specify manually')
+
+            if app.ARGS.echo_time:
+                te_per_series.append(app.ARGS.echo_time)
+                try:
+                    te_per_series_bids = dwi_header.keyval()['EchoTime']
+                    if te_per_series_bids != app.ARGS.echo_time:
+                        warnings.warn('User defined echo time does not match that found in bids json. Using user defined value')
+                except:
+                    pass
+            else:   
+                try:
+                    te_per_series.append(dwi_header.keyval()['EchoTime'])
+                except:
+                    warnings.warn('No echo time found in header, assuming 0')
+                    te_per_series.append(0)
 
         DWImif = ' '.join(miflist)
         cmd = ('mrcat -axis 3 %s %s/dwi.mif' % (DWImif, app.SCRATCH_DIR))
@@ -355,6 +483,30 @@ def convert_input_data(dwi_metadata):
         raise MRtrixError('No phase files found or phases are in an incompatible format')
 
 
+    dwi_metadata['TE'] = te_per_series
+    dwi_metadata['pf'] = pf_per_series
+    dwi_metadata['pe_dir'] = ped_per_series
+
+    if app.ARGS.bshape:
+        bshape = [float(i) for i in app.ARGS.bshape.rsplit(',')]
+    else:
+        bshape = [1] * len(dwi_metadata['dwi_list'])
+    dwi_metadata['bshape'] = bshape
+    bshape_per_series = dwi_metadata['bshape']
+
+    try:
+        if (len(set(te_per_series)) > 1) and (not app.ARGS.rpe_te):
+            raise MRtrixError('If data has variable echo time and no RPE TE is specified, please use the -rpe_te flag to specify the RPE TE')
+    except:
+        pass
+
+    if not all(x == ped_per_series[0] for x in ped_per_series):
+        raise MRtrixError('input series have different phase encoding directions, series should be processed separately')
+
+    if not all(x == pf_per_series[0] for x in pf_per_series):
+        raise MRtrixError('input series have different partial fourier factors,', +
+                            'series should be processed separately')
+
     # get diffusion header info - check to make sure all values are valid for processing
     dwi_header = image.Header('%s/dwi.mif' % (app.SCRATCH_DIR))
     dwi_size = [ int(s) for s in dwi_header.size() ]
@@ -372,6 +524,8 @@ def convert_input_data(dwi_metadata):
 
     dwi_ind_size = [i + [1] if len(i)==3 else i for i in dwi_ind_size]
 
+    import pdb; pdb.set_trace()
+
     nvols = [i[3] for i in dwi_ind_size]
     for idx,i in enumerate(dwi_n_list):
         if len(dwi_n_list) == 1:

output_tmp.bval

@@ -0,0 +1 @@
+1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000

output_tmp.bvec

@@ -0,0 +1,3 @@
+0.838977888593843 -0.838633701381505 -0.884006280555446 0.88337499857346 0.410882897894761 0.379082928117134 0.334181013394585 0.302920140678514 -0.302524859964087 -0.334104134807124 -0.379005890832917 -0.411239143750211
+0.3572989525 0.3570328728 0.3570321133 0.3572989994 0.3574089112 0.8625878365 0.8630940346 0.3577861662 0.3574088346 0.8625883481 0.8630937514 0.3577861251
+0.410430945515036 -0.411364853483884 0.30176309582217 -0.303291999465761 -0.83870979155511 -0.335019936503161 0.378671015209572 0.883305410161624 -0.883593591050831 -0.379889153280358 0.333801903804951 0.838374293074628

rpg_cpp/fftw-3.3.10/FFTW3Config.cmake

@@ -0,0 +1,17 @@
+# defined since 2.8.3
+if (CMAKE_VERSION VERSION_LESS 2.8.3)
+  get_filename_component (CMAKE_CURRENT_LIST_DIR ${CMAKE_CURRENT_LIST_FILE} PATH)
+endif ()
+
+# Allows loading FFTW3 settings from another project
+set (FFTW3_CONFIG_FILE "${CMAKE_CURRENT_LIST_FILE}")
+
+set (FFTW3_LIBRARIES fftw3)
+set (FFTW3_LIBRARY_DIRS /Users/benaron/Documents/DESIGNER-v2/rpg_cpp/fftw-3.3.10/build/lib)
+set (FFTW3_INCLUDE_DIRS /Users/benaron/Documents/DESIGNER-v2/rpg_cpp/fftw-3.3.10/build/include)
+
+include ("${CMAKE_CURRENT_LIST_DIR}/FFTW3LibraryDepends.cmake")
+
+if (CMAKE_VERSION VERSION_LESS 2.8.3)
+  set (CMAKE_CURRENT_LIST_DIR)
+endif ()

rpg_cpp/fftw-3.3.10/FFTW3ConfigVersion.cmake

@@ -0,0 +1,12 @@
+
+set (PACKAGE_VERSION "3.3.10")
+
+# Check whether the requested PACKAGE_FIND_VERSION is compatible
+if ("${PACKAGE_VERSION}" VERSION_LESS "${PACKAGE_FIND_VERSION}")
+  set (PACKAGE_VERSION_COMPATIBLE FALSE)
+else ()
+  set (PACKAGE_VERSION_COMPATIBLE TRUE)
+  if ("${PACKAGE_VERSION}" VERSION_EQUAL "${PACKAGE_FIND_VERSION}")
+    set (PACKAGE_VERSION_EXACT TRUE)
+  endif ()
+endif ()