Add option for plot; Increase flexibility of input args; Add comments…

parse_vNav_Motion.py

@@ -1,21 +1,41 @@
+# Author: Dylan Tisdall
+# https://github.com/MRIMotionCorrection/parse_vNav_Motion
+#
+# Additions 11/22/2017 by Michael Harms
+#   Options for plots and intermediate outputs
+#   Allowed both --rms and --max simultaneously as input args
+#   Added 'radius' as optional input arg
+#   Added comments to help navigate script
+
 import dicom
 import os
 import numpy as np
 import itertools
 import glob
 import argparse
+import matplotlib.pyplot as plt
+from pprint import pprint
+
 
 def normalize(x):
   return x / np.sqrt(np.dot(x,x))
 
-def readRotAndTrans(paths):
+def readRotAndTrans(paths,verbose):
   files = itertools.chain.from_iterable([glob.glob(path) for path in paths])
 
+  # Sort DICOMs by AcquisitionNumber (important!)
   ds = sorted([dicom.read_file(x) for x in files], key=lambda dcm: dcm.AcquisitionNumber)
- 
+
   head = [(np.array([1,0,0,0]),np.array([0,0,0]))]
 
-  return list(itertools.chain.from_iterable([head, [(np.array(map(float, y[1:5])), map(float, y[6:9])) for y in [str.split(x.ImageComments) for x in ds[1:]]]]))
+  # Return the "raw" quaternion values stored in the DICOM ImageComments field
+  quaternions = list(itertools.chain.from_iterable([head, [(np.array(map(float, y[1:5])), map(float, y[6:9])) for y in [str.split(x.ImageComments) for x in ds[1:]]]]))
+
+  if verbose:
+    print 'Quaternions:'
+    pprint(quaternions)
+
+  return quaternions
 
 def angleAxisToQuaternion(a):
   w = np.cos(a[0] / 2.0)
@@ -127,6 +147,10 @@ def motionEntryToHomogeneousTransform(e) :
   t[3,:] = [0,0,0,1]
   return np.matrix(t)
 
+# MaxMotion measure is similar to the motion score in Tisdall et al. 2012 (MRM, 68:389-399),
+# but not identical. In particular, the definition of t_rotmax here is the same as the 
+# computation in Eq. (1) of the paper, but Eq. (3) from the paper (the final score) 
+# would be simply: score = t_rotmax + np.linalg.norm(trans)
 def diffTransformToMaxMotion(t, radius):
   angleAxis = quaternionToAxisAngle(rotationMatrixToQuaternion(t[0:3, 0:3]))
   angle = angleAxis[0]
@@ -141,6 +165,10 @@ def diffTransformToMaxMotion(t, radius):
     np.linalg.norm(trans) * np.linalg.norm(trans)
     )
 
+# RMSmotion measure based on:
+# Jenkinson, M., 1999. Measuring transformation error by RMS deviation. 
+# Technical Report No. TR99MJ1. FMRIB, Oxford
+# http://www.fmrib.ox.ac.uk/datasets/techrep/tr99mj1/tr99mj1/tr99mj1.html
 def diffTransformToRMSMotion(t, radius):
   rotMatMinusIdentity = t[0:3,0:3] - np.array([[1,0,0],[0,1,0],[0,0,1]])
   trans = np.ravel(t[0:3,3])
@@ -150,28 +178,56 @@ def diffTransformToRMSMotion(t, radius):
     np.dot(trans, trans)    
     )
 
+def simplePlot(scores, measStr, radius, TR):
+  fig, ax = plt.subplots()
+  ax.plot(scores)
+  xlabel = 'Frame (TR=' + str(TR) + ' s)'
+  ylabel = 'mm (assumed radius=' + str(radius) + ' mm)'
+  titleStr = 'vNavMotionScores (measure=' + measStr + ')'
+  ax.set(xlabel=xlabel, ylabel=ylabel, title=titleStr)
+  fig.savefig('vNavMotionScores'+measStr+'.png')
+
+
+## Parse arguments
 parser = argparse.ArgumentParser()
 
-parser.add_argument('--tr', required=True)
-parser.add_argument('--input', nargs='+', required=True)
-output_type = parser.add_mutually_exclusive_group(required=True)
-output_type.add_argument('--rms', action='store_true')
-output_type.add_argument('--max', action='store_true')
+group1 = parser.add_argument_group('required')
+group1.add_argument('--tr', type=float, required=True, help='TR (sec)')
+group1.add_argument('--input', nargs='+', required=True, help='List of input DICOMs')
+group2 = parser.add_argument_group('measure (at least one required)')
+group2.add_argument('--rms', action='store_true')
+group2.add_argument('--max', action='store_true')
+# Optional arguments
+parser.add_argument('--radius', type=float, default=100, help='radius (mm) of assumed sphere (default: %(default)s)')
+parser.add_argument('--plot', help='output plot of chosen measure across frames', action='store_true')
+parser.add_argument('-v','--verbose', help='increase output verbosity', action='store_true')
 
 args = parser.parse_args()
 
-transforms = [motionEntryToHomogeneousTransform(e) for e in readRotAndTrans(args.input)]
+if args.rms is False and args.max is False:
+  parser.error('At least one of --rms and --max is required.')
 
-diffTransforms = [ts[1] * np.linalg.inv(ts[0]) for ts in zip(transforms[0:], transforms[1:])]
+## Perform calcs and generate output
+transforms = [motionEntryToHomogeneousTransform(e) for e in readRotAndTrans(args.input,args.verbose)]
 
-rmsMotionScores = [diffTransformToRMSMotion(t, 100) for t in diffTransforms]
+diffTransforms = [ts[1] * np.linalg.inv(ts[0]) for ts in zip(transforms[0:], transforms[1:])]
 
-maxMotionScores = [diffTransformToMaxMotion(t, 100) for t in diffTransforms]
+rmsMotionScores = [diffTransformToRMSMotion(t, args.radius) for t in diffTransforms]
 
+maxMotionScores = [diffTransformToMaxMotion(t, args.radius) for t in diffTransforms]
 
 if args.rms :
-  print np.mean(rmsMotionScores) * 60.0 / float(args.tr)
-elif args.max :
-  print np.mean(maxMotionScores) * 60.0 / float(args.tr)
+  if args.verbose:
+    print 'vNavMotionScoresRMS:', rmsMotionScores
+  print 'MeanMotionScoreRMSPerMin:', np.mean(rmsMotionScores) * 60.0 / float(args.tr)
+  if args.plot:
+    simplePlot(rmsMotionScores, 'RMS', args.radius, args.tr)
+
+if args.max :
+  if args.verbose:
+    print 'vNavMotionScoresMax:', maxMotionScores
+  print 'MeanMotionScoreMaxPerMin:', np.mean(maxMotionScores) * 60.0 / float(args.tr)
+  if args.plot:
+    simplePlot(maxMotionScores, 'Max', args.radius, args.tr)