From e40e9ae554d77ac81f2315ef52126b1c1911f5f2 Mon Sep 17 00:00:00 2001 From: Michael Harms Date: Wed, 22 Nov 2017 15:12:06 -0600 Subject: [PATCH 1/2] Add option for plot; Increase flexibility of input args; Add comments to script --- parse_vNav_Motion.py | 86 ++++++++++++++++++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 15 deletions(-) diff --git a/parse_vNav_Motion.py b/parse_vNav_Motion.py index f92ce34..5fa372f 100644 --- a/parse_vNav_Motion.py +++ b/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 = 'MotionScore (measure=' + measStr + ')' + ax.set(xlabel=xlabel, ylabel=ylabel, title=titleStr) + fig.savefig('MotionScore'+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 'IndividualFrameScoresRMS:', rmsMotionScores + print 'MeanScoreRMSPerMin:', 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 'IndividualFrameScoresMax:', maxMotionScores + print 'MeanScoreMaxPerMin:', np.mean(maxMotionScores) * 60.0 / float(args.tr) + if args.plot: + simplePlot(maxMotionScores, 'Max', args.radius, args.tr) From 247eae3d95c8c28e54b3dc80d23e3f3954b14f05 Mon Sep 17 00:00:00 2001 From: Michael Harms Date: Wed, 22 Nov 2017 15:44:50 -0600 Subject: [PATCH 2/2] Minor naming changes --- parse_vNav_Motion.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/parse_vNav_Motion.py b/parse_vNav_Motion.py index 5fa372f..dd2cb81 100644 --- a/parse_vNav_Motion.py +++ b/parse_vNav_Motion.py @@ -183,9 +183,9 @@ def simplePlot(scores, measStr, radius, TR): ax.plot(scores) xlabel = 'Frame (TR=' + str(TR) + ' s)' ylabel = 'mm (assumed radius=' + str(radius) + ' mm)' - titleStr = 'MotionScore (measure=' + measStr + ')' + titleStr = 'vNavMotionScores (measure=' + measStr + ')' ax.set(xlabel=xlabel, ylabel=ylabel, title=titleStr) - fig.savefig('MotionScore'+measStr+'.png') + fig.savefig('vNavMotionScores'+measStr+'.png') ## Parse arguments @@ -218,15 +218,15 @@ def simplePlot(scores, measStr, radius, TR): if args.rms : if args.verbose: - print 'IndividualFrameScoresRMS:', rmsMotionScores - print 'MeanScoreRMSPerMin:', np.mean(rmsMotionScores) * 60.0 / float(args.tr) + 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 'IndividualFrameScoresMax:', maxMotionScores - print 'MeanScoreMaxPerMin:', np.mean(maxMotionScores) * 60.0 / float(args.tr) + print 'vNavMotionScoresMax:', maxMotionScores + print 'MeanMotionScoreMaxPerMin:', np.mean(maxMotionScores) * 60.0 / float(args.tr) if args.plot: simplePlot(maxMotionScores, 'Max', args.radius, args.tr)