dwicat: Support 3D images

bin/dwicat

@@ -17,7 +17,7 @@
 
 
 
-import json
+import json, shutil
 
 
 
@@ -36,7 +36,7 @@ def usage(cmdline): #pylint: disable=unused-variable
 
 
 def execute(): #pylint: disable=unused-variable
-  from mrtrix3 import MRtrixError #pylint: disable=no-name-in-module, import-outside-toplevel
+  from mrtrix3 import CONFIG, MRtrixError #pylint: disable=no-name-in-module, import-outside-toplevel
   from mrtrix3 import app, image, path, run #pylint: disable=no-name-in-module, import-outside-toplevel
 
   num_inputs = len(app.ARGS.inputs)
@@ -45,13 +45,26 @@ def execute(): #pylint: disable=unused-variable
 
   # check input data
   def check_header(header):
-    if len(header.size()) != 4:
-      raise MRtrixError('Image "' + header.name() + '" is not a 4D image series')
+    if len(header.size()) > 4:
+      raise MRtrixError('Image "' + header.name() + '" contains more than 4 dimensions')
     if not 'dw_scheme' in header.keyval():
       raise MRtrixError('Image "' + header.name() + '" does not contain a gradient table')
-    num_grad_lines = len(header.keyval()['dw_scheme'])
-    if num_grad_lines != header.size()[3]:
-      raise MRtrixError('Number of lines in gradient table for image "' + header.name() + '" (' + str(num_grad_lines) + ') does not match number of volumes (' + str(header.size()[3]) + ')')
+    dw_scheme = header.keyval()['dw_scheme']
+    try:
+      if isinstance(dw_scheme[0], list):
+        num_grad_lines = len(dw_scheme)
+      elif (isinstance(dw_scheme[0], ( int, float))) and len(dw_scheme) >= 4:
+        num_grad_lines = 1
+      else:
+        raise MRtrixError
+    except (IndexError, MRtrixError):
+      raise MRtrixError('Image "' + header.name() + '" contains gradient table of unknown format')
+    if len(header.size()) == 4:
+      num_volumes = header.size()[3]
+      if num_grad_lines != num_volumes:
+        raise MRtrixError('Number of lines in gradient table for image "' + header.name() + '" (' + str(num_grad_lines) + ') does not match number of volumes (' + str(num_volumes) + ')')
+    elif not (num_grad_lines == 1 and len(dw_scheme) >= 4 and dw_scheme[3] <= float(CONFIG.get('BZeroThreshold', 10.0))):
+      raise MRtrixError('Image "' + header.name() + '" is 3D, and cannot be validated as a b=0 volume')
 
   first_header = image.Header(path.from_user(app.ARGS.inputs[0], False))
   check_header(first_header)
@@ -87,7 +100,12 @@ def execute(): #pylint: disable=unused-variable
 
   # extract b=0 volumes within each input series
   for index in range(0, num_inputs):
-    run.command('dwiextract ' + str(index) + 'in.mif ' + str(index) + 'b0.mif -bzero')
+    infile = str(index) + 'in.mif'
+    outfile = str(index) + 'b0.mif'
+    if len(image.Header(infile).size()) > 3:
+      run.command('dwiextract ' + infile + ' ' + outfile + ' -bzero')
+    else:
+      run.function(shutil.copyfile, infile, outfile)
 
   mask_option = ' -mask_input mask.mif -mask_target mask.mif' if app.ARGS.mask else ''
 

cmd/dwiextract.cpp

@@ -67,6 +67,8 @@ void usage ()
 void run()
 {
   auto input_image = Image<float>::open (argument[0]);
+  if (input_image.ndim() < 4)
+    throw Exception ("Epected input image to contain more than three dimensions");
   auto grad = DWI::get_DW_scheme (input_image);
 
   // Want to support non-shell-like data if it's just a straight extraction

testing/scripts/tests/dwicat

@@ -2,4 +2,4 @@ mkdir -p ../tmp/dwicat && mrconvert BIDS/sub-01/dwi/sub-01_dwi.nii.gz -fslgrad B
 dwicat tmp01_b1000.mif tmp03_b2000.mif tmp03_b3000.mif ../tmp/dwicat/sharedb0_unmasked.mif -force && testing_diff_image ../tmp/dwicat/sharedb0_unmasked.mif tmp02.mif -frac 1e-6
 dwiextract tmp.mif tmp11_b0.mif -shell 0 -force && dwiextract tmp.mif tmp11_b1000.mif -shell 1000 -force && dwiextract tmp.mif tmp11_b2000.mif -shell 2000 -force && dwiextract tmp.mif tmp11_b3000.mif -shell 3000 -force && mrconvert tmp11_b0.mif -coord 3 0,1 - | mrcat - tmp11_b1000.mif tmp12_b1000.mif -axis 3 -force && mrconvert tmp11_b0.mif -coord 3 2,3 - | mrcat - tmp11_b2000.mif tmp12_b2000.mif -axis 3 -force && mrconvert tmp11_b0.mif -coord 3 4,5 - | mrcat - tmp11_b3000.mif tmp12_b3000.mif -axis 3 -force && mrcalc tmp12_b2000.mif 0.2 -mult tmp13_b2000.mif -force && mrcalc tmp12_b3000.mif 5.0 -mult tmp13_b3000.mif -force && mrcat tmp12_b1000.mif tmp12_b2000.mif tmp12_b3000.mif tmp14.mif -axis 3 -force && dwicat tmp12_b1000.mif tmp13_b2000.mif tmp13_b3000.mif ../tmp/dwicat/ownb0_masked.mif -mask BIDS/sub-01/dwi/sub-01_brainmask.nii.gz -force && testing_diff_image ../tmp/dwicat/ownb0_masked.mif tmp14.mif -frac 0.01
 dwicat tmp12_b1000.mif tmp13_b2000.mif tmp13_b3000.mif ../tmp/dwicat/ownb0_unmasked.mif -force && testing_diff_image ../tmp/dwicat/ownb0_unmasked.mif tmp14.mif -frac 0.02
-
+dwiextract tmp12_b1000.mif -bzero - | mrconvert - tmp12_b1000_1b0.mif -coord 3 0 -axes 0,1,2 -force && dwicat tmp12_b1000.mif tmp12_b1000_1b0.mif ../tmp/dwicat/3dinput.mif -force