6.0.2 changes

6.0.2 changes

6.0.2 changes

EddyHigh_b_Test/EddyHigh_b_Feeds.py

+#!/usr/bin/env fslpython
+
+import sys
+import numpy as np
+import nibabel as nib
+import datetime
+import pandas as pd
+
+class EddyHigh_b_FeedsType(object):
+    """ The purpose of this class is to make all the comparisons between
+        newly estimated and precomputed eddy results for plain vanilla eddy
+        on high b-value data.
+    """
+
+    def __init__(self,mask,corr,precomp_corr):
+        # Read corrected images and make sure dimensions are right
+        try:
+            self._mask = nib.load(mask,mmap=False)
+            self._corr = nib.load(corr,mmap=False)
+            self._precomp_corr = nib.load(precomp_corr,mmap=False)
+        except Exception as e:
+            print(str(e))
+            raise Exception('EddyHigh_b_FeedsType:__init__:Error opening corrected image files')
+
+        if not (all(self._mask.header['dim'][1:4] == self._corr.header['dim'][1:4]) and
+                all(self._mask.header['dim'][1:4] == self._precomp_corr.header['dim'][1:4])):
+            raise Exception('EddyHigh_b_FeedsType:__init__:Size mismatch in first three dimensions of corrected images')
+        if not (self._mask.header['dim'][4] == 1 and
+                self._corr.header['dim'][4] == self._precomp_corr.header['dim'][4]):
+            raise Exception('EddyHigh_b_FeedsType:__init__:Size mismatch in fourth dimension of corrected images')
+
+        # Compare images and create statistics of differences
+        try:
+            mask = self._mask.get_data().astype(float)
+            mask = (mask > 0).astype(float)
+            corrdiff = self._corr.get_data().astype(float)
+            corrdiff = abs(corrdiff - self._precomp_corr.get_data().astype(float))
+
+            self._corrdiffmeans = np.zeros(corrdiff.shape[3])
+            for vol in range(0, corrdiff.shape[3]):
+                tmpdiff = np.multiply(mask,corrdiff[:,:,:,vol])
+                self._corrdiffmeans[vol] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+
+        except Exception as e:
+            print(str(e))
+            raise Exception('EddyHigh_b_FeedsType:__init__:Error calculating image statistics')
+
+    def MeanDiffsOfCorrectedImages(self):
+        return self._corrdiffmeans
+
+def main(argv):
+    # This is the main program that tests the output from running
+    # vanilla eddy on high b-value data.
+    try:
+        if len(argv) != 8:
+            print('EddyHigh_b_Feeds.py usage: EddyHigh_b_Feeds.py output_dir prefix mask corrected precomputed_corrected' \
+                  'allowed_mean_diff_corrected allowed_max_diff_corrected ')
+            sys.exit(1)
+        else:
+            output_dir = argv[1]
+            output_prefix = argv[2]
+            mask = argv[3]
+            corrected = argv[4]
+            precomputed_corrected = argv[5]
+            allowed_mean_diff_corrected = float(argv[6])
+            allowed_max_diff_corrected = float(argv[7])
+
+        # Try to create EddyFeedsType object (involves reading all files)
+        try:
+            ef = EddyHigh_b_FeedsType(mask,corrected,precomputed_corrected)
+        except Exception as e:
+            print(str(e))
+            print('main: Error when creating EddyHigh_b_FeedsType object.')
+            sys.exit(1)
+
+        try:
+            passes_test = True
+            # Check pass/fail based on corrected images and fields
+            if ef.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+               ef.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected:
+                passes_test = False
+
+        except Exception as e:
+            print(str(e))
+            print('main: Failed calculating stats for test.')
+            sys.exit(1)
+
+        # Write report
+        try:
+            fp = open(output_dir + '/' + output_prefix + '_EddyHigh_b_Report.txt','w')
+        except Exception as e:
+            print(str(e))
+            print('main: Cannot open report file: ' + output_dir + '/' + output_prefix + '_EddyHigh_b_Report.txt')
+            sys.exit(1)
+        else:
+            try:
+                fp.write('EddyHigh_b_Feeds was run on ' + datetime.datetime.now().strftime("%Y-%m-%d %H:%M") + '\n')
+                fp.write('With the command' + ' '.join(argv) + '\n')
+                if passes_test:
+                    fp.write('\nOverall the test passed\n')
+                else:
+                    fp.write('\nOverall the test failed\n')
+                # Report on differences in corrected images
+                fp.write('\nThe absolute differences, averaged across the mask, for the corrected images were: ' + \
+                         ' '.join(["{0:.4f}".format(elem) for elem in ef.MeanDiffsOfCorrectedImages()]) + '\n')
+                fp.write('That gives mean error: ' + "{0:.4f}".format(ef.MeanDiffsOfCorrectedImages().mean()) + \
+                         ', and a max error: ' + "{0:.4f}".format(ef.MeanDiffsOfCorrectedImages().max()) + '\n')
+                fp.write('The allowed mean of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_mean_diff_corrected) + '\n')
+                fp.write('The allowed maximum of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_max_diff_corrected) + '\n')
+                if ef.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+                   ef.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected:
+                    fp.write('Based on these criteria the test failed\n')
+                else:
+                    fp.write('Based on these criteria the test passed\n')
+                fp.close()
+            except Exception as e:
+                print(str(e))
+                print('main: Problem writing report file: ' + argv[1] + '/' + argv[2] + '_EddyHigh_b_Report.txt')
+                sys.exit(1)
+    except Exception as e:
+        print(str(e))
+        print('main: Unknown problem in body of function')
+        sys.exit(1)
+
+    if passes_test:
+        sys.exit(0)
+    else:
+        sys.exit(1)
+
+
+if __name__ == "__main__":
+    main(sys.argv)
\ No newline at end of file

EddyHigh_b_Test/feedsInputs

+EddyHigh_b_TestData/eddyData

EddyHigh_b_Test/feedsRun

+#! /bin/sh
+#
+# This script runs eddy followed by some
+# tests on the output.
+# It uses some b=7000 data from the pilot phase
+# of the HCP. The "ground truth" is simply the 
+# output from eddy at a time when it seemed to
+# do a good job. The tolerance is quite high. That
+# reflects both that the data has a quite high
+# scale (b=0 values ~2000) and also that the
+# run-to-run variability of eddy is quite high
+# on these data.
+#
+
+odir=$1
+strict=$2
+indir=$3
+#
+# Inputs 1--3 are the once neccessary for feeds to work
+# Additional inputs are optional and intended for testing
+# outside of the feeds context.
+#
+# Input 4 is alternative location of executable
+#
+if [ "$#" -gt 3 ]; then
+    exedir=$4
+   echo "Directory for eddy executables set to ${exedir}"
+else
+   exedir="${FSLDIR}/bin"
+fi 
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    if [ ! -x "$cuda_exe" ]; then
+	echo "No executable ${exedir}/eddy_cuda* exists"
+	exit 1
+    fi
+done
+if [ ! -x "${exedir}/eddy_openmp" ]; then
+   echo "Executable ${exedir}/eddy_openmp does not exist"
+   exit 1
+fi
+ 
+if [ ! -d "$odir" ]; then
+   echo "Output directory ${odir} does not exist" 
+   exit 1;
+fi
+if [ ! -d "$indir" ]; then
+   echo "Input directory ${indir} does not exist" 
+   exit 1;
+fi
+
+# Launch both GPU and CPU versions
+
+cuda_jid=""
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    cuda_jid="$cuda_jid `fsl_sub -q cuda.q ${cuda_exe} --imain=${indir}/EddyHigh_b_TestData/eddyData/testData --mask=${indir}/EddyHigh_b_TestData/eddyData/testMask --bvals=${indir}/EddyHigh_b_TestData/eddyData/testBvals --bvecs=${indir}/EddyHigh_b_TestData/eddyData/testBvecs --index=${indir}/EddyHigh_b_TestData/eddyData/testIndex --acqp=${indir}/EddyHigh_b_TestData/eddyData/testAcqparams --topup=${indir}/EddyHigh_b_TestData/eddyData/testTopup --nvoxhp=5000 --repol --fwhm=10,0,0,0,0 --dont_peas --out=${odir}/eddyCudaOutput${version} --very_verbose`"
+done
+
+omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyHigh_b_TestData/eddyData/testData --mask=${indir}/EddyHigh_b_TestData/eddyData/testMask --bvals=${indir}/EddyHigh_b_TestData/eddyData/testBvals --bvecs=${indir}/EddyHigh_b_TestData/eddyData/testBvecs --index=${indir}/EddyHigh_b_TestData/eddyData/testIndex --acqp=${indir}/EddyHigh_b_TestData/eddyData/testAcqparams --topup=${indir}/EddyHigh_b_TestData/eddyData/testTopup --nvoxhp=5000 --repol --fwhm=10,0,0,0,0 --dont_peas --out=${odir}/eddyOmpOutput --very_verbose`
+# Ensure that slots are being reserved on the queue
+qalter ${omp_jid} -R y
+
+# Wait until they have both finished
+
+while [ : ]; do
+    for jid in ${cuda_jid};
+    do
+	tmp=`qstat -j ${jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -ne 0 ]; then
+	    break
+	fi
+    done
+    if [ $tmp -eq 0 ]; then
+	tmp=`qstat -j ${omp_jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -eq 0 ]; then
+      	    break
+	fi
+    fi
+    sleep 5m
+done
+
+# Check that eddy output exists
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    if [ ! -f ${odir}/eddyCudaOutput${version}.nii* ]; then
+	exit 1
+    fi
+done
+if [ ! -f ${odir}/eddyOmpOutput.nii* ]; then
+   exit 1
+fi
+
+# Define some constants
+
+max_mean_ima_diff=15.0
+max_max_ima_diff=20.0
+
+# Check the results against ground truth
+
+cuda_exit_status=0
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    ./EddyHigh_b_Feeds.py ${odir} Cuda${version} `imglob -extension ${indir}/EddyHigh_b_TestData/eddyData/testMask` `imglob -extension ${odir}/eddyCudaOutput${version}` `imglob -extension ${indir}/EddyHigh_b_TestData/eddyData/Precomputed/eddy_cuda_output` $max_mean_ima_diff $max_max_ima_diff
+    cuda_exit_status=$(($cuda_exit_status + $?))
+done
+
+./EddyHigh_b_Feeds.py ${odir} Omp `imglob -extension ${indir}/EddyHigh_b_TestData/eddyData/testMask` `imglob -extension ${odir}/eddyOmpOutput` `imglob -extension ${indir}/EddyHigh_b_TestData/eddyData/Precomputed/eddy_openmp_output` $max_mean_ima_diff $max_max_ima_diff
+
+omp_exit_status=$?
+
+if [ $cuda_exit_status -gt 0 ] || [ $omp_exit_status -gt 0 ]; then
+    exit 1
+else
+    exit 0
+fi
+
+

EddyHigh_b_Test/nohup.sh

+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyHigh_b blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?

EddyLSRTest/EddyLSRFeeds.py

+#!/usr/bin/env fslpython
+
+import sys
+import numpy as np
+import nibabel as nib
+import datetime
+
+class EddyLSRFeedsType(object):
+    """ The purpose of this class is to make all the comparisons between
+        newly estimated and precomputed eddy results for plain vanilla eddy
+        with LSR-resampling.
+    """
+
+    def __init__(self,mask,corr,precomp_corr):
+        # Read corrected images and make sure dimensions are right
+        try:
+            self._mask = nib.load(mask,mmap=False)
+            self._corr = nib.load(corr,mmap=False)
+            self._precomp_corr = nib.load(precomp_corr,mmap=False)
+        except Exception as e:
+            print(str(e))
+            raise Exception('EddyLSRFeedsType:__init__:Error opening corrected image files')
+
+        if not (all(self._mask.header['dim'][1:4] == self._corr.header['dim'][1:4]) and
+                all(self._mask.header['dim'][1:4] == self._precomp_corr.header['dim'][1:4])):
+            raise Exception('EddyLSRFeedsType:__init__:Size mismatch in first three dimensions of corrected images')
+        if not (self._mask.header['dim'][4] == 1 and
+                self._corr.header['dim'][4] == self._precomp_corr.header['dim'][4]):
+            raise Exception('EddyLSRFeedsType:__init__:Size mismatch in fourth dimension of corrected images')
+
+        # Compare images and create statistics of differences
+        try:
+            mask = self._mask.get_data().astype(float)
+            mask = (mask > 0).astype(float)
+            corrdiff = self._corr.get_data().astype(float)
+            corrdiff = abs(corrdiff - self._precomp_corr.get_data().astype(float))
+
+            self._corrdiffmeans = np.zeros(corrdiff.shape[3])
+            for vol in range(0, corrdiff.shape[3]):
+                tmpdiff = np.multiply(mask,corrdiff[:,:,:,vol])
+                self._corrdiffmeans[vol] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+
+        except Exception as e:
+            print(str(e))
+            raise Exception('EddyLSRFeedsType:__init__:Error calculating image statistics')
+
+    def MeanDiffsOfCorrectedImages(self):
+        return self._corrdiffmeans
+
+def main(argv):
+    # This is the main program that tests the output from running
+    # vanilla eddy with LSR resampling.
+    try:
+        if len(argv) != 8:
+            print('EddyLSRFeeds.py usage: EddyLSRFeeds.py output_dir prefix mask corrected precomputed_corrected' \
+                  'allowed_mean_diff_corrected allowed_max_diff_corrected ')
+            sys.exit(1)
+        else:
+            output_dir = argv[1]
+            output_prefix = argv[2]
+            mask = argv[3]
+            corrected = argv[4]
+            precomputed_corrected = argv[5]
+            allowed_mean_diff_corrected = float(argv[6])
+            allowed_max_diff_corrected = float(argv[7])
+
+        # Try to create EddyFeedsType object (involves reading all files)
+        try:
+            ef = EddyLSRFeedsType(mask,corrected,precomputed_corrected)
+        except Exception as e:
+            print(str(e))
+            print('main: Error when creating EddyLSRFeedsType object.')
+            sys.exit(1)
+
+        try:
+            passes_test = True
+            # Check pass/fail based on corrected images and fields
+            if ef.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+               ef.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected:
+                passes_test = False
+
+        except Exception as e:
+            print(str(e))
+            print('main: Failed calculating stats for test.')
+            sys.exit(1)
+
+        # Write report
+        try:
+            fp = open(output_dir + '/' + output_prefix + '_EddyLSRReport.txt','w')
+        except Exception as e:
+            print(str(e))
+            print('main: Cannot open report file: ' + output_dir + '/' + output_prefix + '_EddyLSRReport.txt')
+            sys.exit(1)
+        else:
+            try:
+                fp.write('EddyLSRFeeds was run on ' + datetime.datetime.now().strftime("%Y-%m-%d %H:%M") + '\n')
+                fp.write('With the command' + ' '.join(argv) + '\n')
+                if passes_test:
+                    fp.write('\nOverall the test passed\n')
+                else:
+                    fp.write('\nOverall the test failed\n')
+                # Report on differences in corrected images
+                fp.write('\nThe absolute differences, averaged across the mask, for the corrected images were: ' + \
+                         ' '.join(["{0:.4f}".format(elem) for elem in ef.MeanDiffsOfCorrectedImages()]) + '\n')
+                fp.write('That gives mean error: ' + "{0:.4f}".format(ef.MeanDiffsOfCorrectedImages().mean()) + \
+                         ', and a max error: ' + "{0:.4f}".format(ef.MeanDiffsOfCorrectedImages().max()) + '\n')
+                fp.write('The allowed mean of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_mean_diff_corrected) + '\n')
+                fp.write('The allowed maximum of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_max_diff_corrected) + '\n')
+                if ef.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+                   ef.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected:
+                    fp.write('Based on these criteria the test failed\n')
+                else:
+                    fp.write('Based on these criteria the test passed\n')
+                fp.close()
+            except Exception as e:
+                print(str(e))
+                print('main: Problem writing report file: ' + argv[1] + '/' + argv[2] + '_EddyLSRReport.txt')
+                sys.exit(1)
+    except Exception as e:
+        print(str(e))
+        print('main: Unknown problem in body of function')
+        sys.exit(1)
+
+    if passes_test:
+        sys.exit(0)
+    else:
+        sys.exit(1)
+
+
+if __name__ == "__main__":
+    main(sys.argv)
\ No newline at end of file

EddyLSRTest/LaunchEddy.sh

+#! /bin/bash
+
+export exedir="/home/fs0/jesper/fsl/src/BuildingEddy_5_0_11_for_PreRelease/eddy/"
+export indir="/vols/Data/fsldev/dataSets/EddyLSRTestData/eddyData"
+export odir="/vols/Data/fsldev/dataSets/EddyLSRTestData/eddyData/Precomputed"
+
+fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/testData --mask=${indir}/testMask --bvals=${indir}/testBvals --bvecs=${indir}/testBvecs --index=${indir}/testIndex --acqp=${indir}/testAcqparams --topup=${indir}/testTopup --out=${odir}/eddy_openmp_output --resamp=lsr --fep --nvoxhp=5000 --dont_sep_offs_move --very_verbose
+
+fsl_sub -q cuda.q ${exedir}/eddy_cuda7.5 --imain=${indir}/testData --mask=${indir}/testMask --bvals=${indir}/testBvals --bvecs=${indir}/testBvecs --index=${indir}/testIndex --acqp=${indir}/testAcqparams --topup=${indir}/testTopup --out=${odir}/eddy_cuda_output --resamp=lsr --fep --nvoxhp=5000 --dont_sep_offs_move --very_verbose
+
+#
+# Run new eddy as well, as a kind of pre-test
+#
+
+export exedir="/home/fs0/jesper/DebugEddyNewNewimage/eddy_new/eddy"
+export indir="/vols/Data/fsldev/dataSets/EddyLSRTestData/eddyData"
+export odir="/vols/Scratch/HCP/Diffusion/jesper/EddyLSRFeedsTest20190129"
+
+fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/testData --mask=${indir}/testMask --bvals=${indir}/testBvals --bvecs=${indir}/testBvecs --index=${indir}/testIndex --acqp=${indir}/testAcqparams --topup=${indir}/testTopup --out=${odir}/eddy_openmp_output --resamp=lsr --fep --nvoxhp=5000 --dont_sep_offs_move --very_verbose
+
+fsl_sub -q cuda.q ${exedir}/eddy_cuda --imain=${indir}/testData --mask=${indir}/testMask --bvals=${indir}/testBvals --bvecs=${indir}/testBvecs --index=${indir}/testIndex --acqp=${indir}/testAcqparams --topup=${indir}/testTopup --out=${odir}/eddy_cuda_output --resamp=lsr --fep --nvoxhp=5000 --dont_sep_offs_move --very_verbose
+
+
+
+

EddyLSRTest/feedsInputs

+EddyLSRTestData/eddyData

EddyLSRTest/feedsRun

+#! /bin/sh
+#
+# This script runs eddy followed by some
+# tests on the output.
+# It uses some data from the pilot phase
+# of the HCP where each diffusion direction was
+# scanned twice with opposing PE-directions. 
+# The "ground truth" is simply the output from eddy 
+# prior to Armadillo and NewNewimage.
+#
+
+odir=$1
+strict=$2
+indir=$3
+#
+# Inputs 1--3 are the once neccessary for feeds to work
+# Additional inputs are optional and intended for testing
+# outside of the feeds context.
+#
+# Input 4 is alternative location of executable
+#
+if [ "$#" -gt 3 ]; then
+    exedir=$4
+   echo "Directory for eddy executables set to ${exedir}"
+else
+   exedir="${FSLDIR}/bin"
+fi 
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    if [ ! -x "$cuda_exe" ]; then
+	echo "No executable ${exedir}/eddy_cuda* exists"
+	exit 1
+    fi
+done
+if [ ! -x "${exedir}/eddy_openmp" ]; then
+   echo "Executable ${exedir}/eddy_openmp does not exist"
+   exit 1
+fi
+ 
+if [ ! -d "$odir" ]; then
+   echo "Output directory ${odir} does not exist" 
+   exit 1;
+fi
+if [ ! -d "$indir" ]; then
+   echo "Input directory ${indir} does not exist" 
+   exit 1;
+fi
+
+# Launch both GPU and CPU versions
+
+cuda_jid=""
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    cuda_jid="$cuda_jid `fsl_sub -q cuda.q ${cuda_exe} --imain=${indir}/EddyLSRTestData/eddyData/testData --mask=${indir}/EddyLSRTestData/eddyData/testMask --bvals=${indir}/EddyLSRTestData/eddyData/testBvals --bvecs=${indir}/EddyLSRTestData/eddyData/testBvecs --index=${indir}/EddyLSRTestData/eddyData/testIndex --acqp=${indir}/EddyLSRTestData/eddyData/testAcqparams --topup=${indir}/EddyLSRTestData/eddyData/testTopup --out=${odir}/eddyCudaOutput${version} --resamp=lsr --fep --nvoxhp=5000 --repol --fwhm=10,0,0,0,0 --dont_peas --very_verbose`"
+done
+
+omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyLSRTestData/eddyData/testData --mask=${indir}/EddyLSRTestData/eddyData/testMask --bvals=${indir}/EddyLSRTestData/eddyData/testBvals --bvecs=${indir}/EddyLSRTestData/eddyData/testBvecs --index=${indir}/EddyLSRTestData/eddyData/testIndex --acqp=${indir}/EddyLSRTestData/eddyData/testAcqparams --topup=${indir}/EddyLSRTestData/eddyData/testTopup --out=${odir}/eddyOmpOutput --resamp=lsr --fep --nvoxhp=5000 --repol --fwhm=10,0,0,0,0 --dont_peas --very_verbose`
+# Ensure that slots are being reserved on the queue
+qalter ${omp_jid} -R y
+
+# Wait until they have both finished
+
+while [ : ]; do
+    for jid in ${cuda_jid};
+    do
+	tmp=`qstat -j ${jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -ne 0 ]; then
+	    break
+	fi
+    done
+    if [ $tmp -eq 0 ]; then
+	tmp=`qstat -j ${omp_jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -eq 0 ]; then
+      	    break
+	fi
+    fi
+    sleep 5m
+done
+
+# Check that eddy output exists
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    if [ ! -f ${odir}/eddyCudaOutput${version}.nii* ]; then
+	exit 1
+    fi
+done
+if [ ! -f ${odir}/eddyOmpOutput.nii* ]; then
+   exit 1
+fi
+
+# Define some constants
+
+max_mean_ima_diff=15.0
+max_max_ima_diff=20.0
+
+# Check the results against ground truth
+
+cuda_exit_status=0
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    ./EddyLSRFeeds.py ${odir} Cuda${version} `imglob -extension ${indir}/EddyLSRTestData/eddyData/testMask` `imglob -extension ${odir}/eddyCudaOutput${version}` `imglob -extension ${indir}/EddyLSRTestData/eddyData/Precomputed/eddy_cuda_output` $max_mean_ima_diff $max_max_ima_diff
+    cuda_exit_status=$(($cuda_exit_status + $?))
+done
+
+./EddyLSRFeeds.py ${odir} Omp `imglob -extension ${indir}/EddyLSRTestData/eddyData/testMask` `imglob -extension ${odir}/eddyOmpOutput` `imglob -extension ${indir}/EddyLSRTestData/eddyData/Precomputed/eddy_openmp_output` $max_mean_ima_diff $max_max_ima_diff 
+
+omp_exit_status=$?
+
+if [ $cuda_exit_status -gt 0 ] || [ $omp_exit_status -gt 0 ]; then
+    exit 1
+else
+    exit 0
+fi
+
+

EddyLSRTest/nohup.sh

+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyLSRTest blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?

EddyMBSTest/MoveBySuscFeeds.py

+#!/usr/bin/env fslpython
+
+import sys
+import numpy as np
+import nibabel as nib
+import datetime
+
+class MoveBySuscFeedsType(object):
+    """ The purpose of this class is to make all the comparisons between
+        newly estimated and precomputed eddy results for the MoveBySusc
+        feeds for eddy.
+    """
+    def __init__(self,mask,corr,deriv_fields,precomp_corr,precomp_deriv_fields):
+        # Read image files and make sure dimensions are right
+        try:
+            self._mask = nib.load(mask,mmap=False)
+            self._corr = nib.load(corr,mmap=False)
+            self._deriv_fields = nib.load(deriv_fields,mmap=False)
+            self._precomp_corr = nib.load(precomp_corr,mmap=False)
+            self._precomp_deriv_fields = nib.load(precomp_deriv_fields,mmap=False)
+        except Exception as e:
+            print(str(e))
+            raise Exception('MoveBySuscFeeds:__init__:Error opening image files')
+
+        if not (all(self._mask.header['dim'][1:4]==self._corr.header['dim'][1:4]) and
+                all(self._mask.header['dim'][1:4]==self._deriv_fields.header['dim'][1:4]) and
+                all(self._mask.header['dim'][1:4]==self._precomp_corr.header['dim'][1:4]) and
+                all(self._mask.header['dim'][1:4]==self._precomp_deriv_fields.header['dim'][1:4])):
+            raise Exception('MoveBySuscFeeds:__init__:Size mismatch in first three dimensions')
+        if not (self._mask.header['dim'][4] == 1 and
+                self._corr.header['dim'][4] == self._precomp_corr.header['dim'][4] and
+                self._deriv_fields.header['dim'][4] == self._precomp_deriv_fields.header['dim'][4]) :
+            raise Exception('MoveBySuscFeeds:__init__:Size mismatch in fourth dimension')
+
+        # Compare images and create statistics of differences
+        try:
+            mask = self._mask.get_data().astype(float)
+            mask = (mask > 0).astype(float)
+            corrdiff = self._corr.get_data().astype(float)
+            corrdiff = abs(corrdiff - self._precomp_corr.get_data().astype(float))
+            fielddiff = self._deriv_fields.get_data().astype(float)
+            fielddiff = abs(fielddiff - self._precomp_deriv_fields.get_data().astype(float))
+
+            self._corrdiffmeans = np.zeros(corrdiff.shape[3])
+            for vol in range(0, corrdiff.shape[3]):
+                tmpdiff = np.multiply(mask,corrdiff[:,:,:,vol])
+                self._corrdiffmeans[vol] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+
+            self._fielddiffmeans = np.zeros(fielddiff.shape[3])
+            for vol in range(0, fielddiff.shape[3]):
+                tmpdiff = np.multiply(mask,fielddiff[:,:,:,vol])
+                self._fielddiffmeans[vol] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+        except Exception as e:
+            print(str(e))
+            raise Exception('MoveBySuscFeeds:__init__:Error calculating image statistics')
+
+    def MeanDiffsOfCorrectedImages(self):
+        return self._corrdiffmeans
+
+    def MeanDiffsOfFieldDerivatives(self):
+        return self._fielddiffmeans
+
+def main(argv):
+    # This is the main program that tests the output from running
+    # the MoveBySuscFeeds tests.
+    try:
+        if len(argv) != 12:
+            print('MoveBySuscFeeds.py usage: MoveBySuscFeeds output_dir prefix mask corrected precomputed_corrected fields precomputed_fields ' \
+                  'allowed_mean_diff_corrected allowed_max_diff_corrected allowed_mean_diff_field allowed_max_diff_field')
+            sys.exit(1)
+        else:
+            output_dir = argv[1]
+            output_prefix = argv[2]
+            mask = argv[3]
+            corrected = argv[4]
+            precomputed_corrected = argv[5]
+            fields = argv[6]
+            precomputed_fields = argv[7]
+            allowed_mean_diff_corrected = float(argv[8])
+            allowed_max_diff_corrected = float(argv[9])
+            allowed_mean_diff_field = float(argv[10])
+            allowed_max_diff_field = float(argv[11])
+
+        # Try to create MoveBySuscFeedsType object (involves reading all files)
+        try:
+            mbsf = MoveBySuscFeedsType(mask,corrected,fields,precomputed_corrected,precomputed_fields)
+        except Exception as e:
+            print(str(e))
+            print('main: Error when creating MoveBySuscFeedsType object.')
+            sys.exit(1)
+
+        # Compare stats to pre-defined limits and make a decision on pass or fail
+        try:
+            passes_test = True
+            if mbsf.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+               mbsf.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected or \
+               mbsf.MeanDiffsOfFieldDerivatives().mean() > allowed_mean_diff_field or \
+               mbsf.MeanDiffsOfFieldDerivatives().max() > allowed_max_diff_field:
+                passes_test = False
+        except Exception as e:
+            print(str(e))
+            print('main: Failed calculating stats for test.')
+            sys.exit(1)
+
+        # Write report
+        try:
+            fp = open(output_dir + '/' + output_prefix + '_MoveBySuscFeedsReport.txt','w')
+        except Exception as e:
+            print(str(e))
+            print('main: Cannot open report file: ' + output_dir + '/' + output_prefix + '_MoveBySuscFeedsReport.txt')
+            sys.exit(1)
+        else:
+            try:
+                fp.write('MoveBySuscFeeds was run on ' + datetime.datetime.now().strftime("%Y-%m-%d %H:%M") + '\n')
+                fp.write('With the command' + ' '.join(argv) + '\n')
+                if passes_test:
+                    fp.write('\nOverall the test passed\n')
+                else:
+                    fp.write('\nOverall the test failed\n')
+                fp.write('\nThe absolute differences, averaged across the mask, for the corrected images were: ' + \
+                         ' '.join(["{0:.4f}".format(elem) for elem in mbsf.MeanDiffsOfCorrectedImages()]) + '\n')
+                fp.write('That gives mean error: ' + "{0:.4f}".format(mbsf.MeanDiffsOfCorrectedImages().mean()) + \
+                         ', and a max error: ' + "{0:.4f}".format(mbsf.MeanDiffsOfCorrectedImages().max()) + '\n')
+                fp.write('The allowed mean of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_mean_diff_corrected) + '\n')
+                fp.write('The allowed maximum of averaged differences for the corrected images was: ' + \
+                         "{0:.4f}".format(allowed_max_diff_corrected) + '\n')
+                if mbsf.MeanDiffsOfCorrectedImages().mean() > allowed_mean_diff_corrected or \
+                   mbsf.MeanDiffsOfCorrectedImages().max() > allowed_max_diff_corrected:
+                    fp.write('Based on these criteria the test failed\n')
+                else:
+                    fp.write('Based on these criteria the test passed\n')
+                fp.write('\nThe absolute differences, averaged across the mask, for the estimated derivative images were: ' + \
+                         ' '.join(["{0:.4f}".format(elem) for elem in mbsf.MeanDiffsOfFieldDerivatives()]) + '\n')
+                fp.write('That gives mean error: ' + "{0:.4f}".format(mbsf.MeanDiffsOfFieldDerivatives().mean()) + \
+                         ', and a max error: ' + "{0:.4f}".format(mbsf.MeanDiffsOfFieldDerivatives().max()) + '\n')
+                fp.write('The allowed mean of averaged differences for the estimated derivative images was: ' + \
+                         "{0:.4f}".format(allowed_mean_diff_field) + '\n')
+                fp.write('The allowed maximum of averaged differences for the estimated derivative images was: ' + \
+                         "{0:.4f}".format(allowed_max_diff_field) + '\n')
+                if mbsf.MeanDiffsOfFieldDerivatives().mean() > allowed_mean_diff_field or \
+                   mbsf.MeanDiffsOfFieldDerivatives().max() > allowed_max_diff_field:
+                    fp.write('Based on these criteria the test failed\n')
+                else:
+                    fp.write('Based on these criteria the test passed\n')
+                fp.close()
+            except Exception as e:
+                print(str(e))
+                print('main: Problem writing report file: ' + output_dir + '/' + output_prefix + '_MoveBySuscFeedsReport.txt')
+                sys.exit(1)
+
+    except Exception as e:
+        print(str(e))
+        print('main: Unknown problem in body of function')
+        sys.exit(1)
+
+    # Exit with success or failure
+    if passes_test:
+        sys.exit(0)
+    else:
+        sys.exit(1)
+
+
+
+if __name__ == "__main__":
+    main(sys.argv)

EddyMBSTest/feedsInputs

+EddyMBSTestData/eddyData

EddyMBSTest/feedsRun

+#! /bin/sh
+#
+# This script runs eddy for testing the 
+# movement-by-susceptibility functionality
+# followed by some tests on the output.
+# It uses simulated data supplied by Mark
+# Graham, UCL.
+# The data is a subset of that used for 
+# the MBS paper. We don't have an "actual"
+# ground truth so instead we compare against
+# the results from the paper (calculated with
+# the Cuda version).
+# The actual data used for the test is 
+# .../SecondSetOfSimulations/ap/ec-volumetric-3x/images_SNR40_1
+#
+
+odir=$1
+strict=$2
+indir=$3
+#
+# Inputs 1--3 are the once neccessary for feeds to work
+# Additional inputs are optional and intended for testing
+# outside of the feeds context.
+#
+# Input 4 is alternative location of executable
+#
+if [ "$#" -gt 3 ]; then
+    exedir=$4
+   echo "Directory for eddy executables set to ${exedir}"
+else
+   exedir="${FSLDIR}/bin"
+fi 
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    if [ ! -x "$cuda_exe" ]; then
+	echo "No executable ${exedir}/eddy_cuda* exists"
+	exit 1
+    fi
+done
+if [ ! -x "${exedir}/eddy_openmp" ]; then
+   echo "Executable ${exedir}/eddy_openmp does not exist"
+   exit 1
+fi
+ 
+if [ ! -d "$odir" ]; then
+   echo "Output directory ${odir} does not exist" 
+   exit 1;
+fi
+if [ ! -d "$indir" ]; then
+   echo "Input directory ${indir} does not exist" 
+   exit 1;
+fi
+
+# Define some constants
+
+max_mean_dfield_diff=0.03
+max_max_dfield_diff=0.06
+max_mean_ima_diff=0.3
+max_max_ima_diff=0.6
+
+# Launch both GPU and CPU versions
+
+cuda_jid=""
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    cuda_jid="$cuda_jid `fsl_sub -q cuda.q ${cuda_exe} --imain=${indir}/EddyMBSTestData/eddyData/testData --acqp=${indir}/EddyMBSTestData/eddyData/testAcqparams --mask=${indir}/EddyMBSTestData/eddyData/testMask --index=${indir}/EddyMBSTestData/eddyData/testIndex --bvecs=${indir}/EddyMBSTestData/eddyData/testBvecs --bvals=${indir}/EddyMBSTestData/eddyData/testBvals --topup=${indir}/EddyMBSTestData/eddyData/testTopup --fwhm=10,5,2,0,0,0,0,0 --niter=8 --nvoxhp=2000 --flm=quadratic --dont_peas --estimate_move_by_susceptibility --mbs_niter=20 --mbs_lambda=10 --mbs_ksp=10 --out=${odir}/eddyCudaOutput${version} --very_verbose`"
+done
+
+omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyMBSTestData/eddyData/testData --acqp=${indir}/EddyMBSTestData/eddyData/testAcqparams --mask=${indir}/EddyMBSTestData/eddyData/testMask --index=${indir}/EddyMBSTestData/eddyData/testIndex --bvecs=${indir}/EddyMBSTestData/eddyData/testBvecs --bvals=${indir}/EddyMBSTestData/eddyData/testBvals --topup=${indir}/EddyMBSTestData/eddyData/testTopup --fwhm=10,5,2,0,0,0,0,0 --niter=8 --nvoxhp=2000 --flm=quadratic --dont_peas --estimate_move_by_susceptibility --mbs_niter=20 --mbs_lambda=10 --mbs_ksp=10 --out=${odir}/eddyOmpOutput --very_verbose`
+# Ensure that slots are being reserved on the queue
+qalter ${omp_jid} -R y
+
+# Wait until they have both finished
+
+while [ : ]; do
+    for jid in ${cuda_jid};
+    do
+	tmp=`qstat -j ${jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -ne 0 ]; then
+	    break
+	fi
+    done
+    if [ $tmp -eq 0 ]; then
+	tmp=`qstat -j ${omp_jid} | wc`
+	tmp=`echo $tmp | awk '{print $1}'`
+	if [ $tmp -eq 0 ]; then
+      	    break
+	fi
+    fi
+    sleep 5m
+done
+
+# Check that eddy output exists
+
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    if [ ! -f ${odir}/eddyCudaOutput${version}.nii* ]; then
+	exit 1
+    fi
+done
+if [ ! -f ${odir}/eddyOmpOutput.nii* ]; then
+   exit 1
+fi
+
+# Check the results against precomputed results
+
+cuda_exit_status=0
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+    tmp=`basename $cuda_exe`
+    version=`echo $tmp | sed 's/eddy_cuda//'`
+    ./MoveBySuscFeeds.py ${odir} Cuda${version} `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/BrainMaskForComparison` `imglob -extension ${odir}/eddyCudaOutput${version}` `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/PrecomputedResults` `imglob -extension ${odir}/eddyCudaOutput${version}.eddy_mbs_first_order_fields` `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/PrecomputedResults.eddy_mbs_first_order_fields` ${max_mean_ima_diff} ${max_max_ima_diff} ${max_mean_dfield_diff} ${max_max_dfield_diff}  
+    cuda_exit_status=$(($cuda_exit_status + $?))
+done
+
+./MoveBySuscFeeds.py ${odir} Omp `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/BrainMaskForComparison` `imglob -extension ${odir}/eddyOmpOutput` `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/PrecomputedResults` `imglob -extension ${odir}/eddyOmpOutput.eddy_mbs_first_order_fields` `imglob -extension ${indir}/EddyMBSTestData/eddyData/Precomputed/PrecomputedResults.eddy_mbs_first_order_fields` ${max_mean_ima_diff} ${max_max_ima_diff} ${max_mean_dfield_diff} ${max_max_dfield_diff}  
+
+omp_exit_status=$?
+
+if [ $cuda_exit_status -gt 0 ] || [ $omp_exit_status -gt 0 ]; then
+    exit 1
+else
+    exit 0
+fi
+

EddyMBSTest/nohup.sh

+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyMBSTest blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?

EddyS2VTest/CompareMovementOverTime.py

-#!/usr/bin/env fslpython
-# CompareMovementOverTime.py - script for comparing two .eddy_movement_over_time files.
-#
-"""
-Usage:
-CompareMovementOverTime.py 'test .eddy_movement_over_time' 'ref .eddy_movement_over_time' max_allowed_max_trans_diff max_allowed_99percentile_trans_diff max_allowed_mean_trans_diff max_allowed_max_rot_diff max_allowed_99percentile_rot_diff max_allowed_mean_rot_diff
-or 
-CompareMovementOverTime.py 'test .eddy_movement_over_time' 'ref .eddy_movement_over_time' max_allowed_max_trans_diff max_allowed_99percentile_trans_diff max_allowed_mean_trans_diff max_allowed_max_rot_diff max_allowed_99percentile_rot_diff max_allowed_mean_rot_diff number_of_slices_in_data first_slice_to_check last_slice_to_check
-"""
-
-from __future__ import print_function # I have no words...
-
-import sys
-import CompareTextOutputDef as ct
-
-if __name__ == '__main__':
-    try:
-        if len(sys.argv) == 1:
-            print(__doc__)
-            exit(1)
-
-        if len(sys.argv) == 10:
-            are_the_same = ct.MovementOverTimeAreSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),float(sys.argv[7]),float(sys.argv[8]),int(sys.argv[9]))
-        elif len(sys.argv) == 13:
-            are_the_same = ct.MovementOverTimeAreSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),float(sys.argv[7]),float(sys.argv[8]),int(sys.argv[9]),int(sys.argv[10]),int(sys.argv[11]),int(sys.argv[12]))
-        else:
-            are_the_same = False
-            print('CompareMovementOverTime: Invalid number of input arguments: {}'.format(len(sys.argv)-1),file=sys.stderr)
-
-        
-        if are_the_same == True:
-            sys.exit(0)
-        else:
-            sys.exit(1)
-    except Exception as e:
-        print('CompareMovementOverTime: Unknown error: {}'.format(e.args[0]),file=sys.stderr)
-        sys.exit(1)
-
-    print('CompareMovementOverTime: Unknown error2',file=sys.stderr)
-    sys.exit(1)
-

EddyS2VTest/CompareMovementOverTime.py~

-#!/usr/bin/env fslpython
-# CompareMovementOverTime.py - script for comparing two .eddy_movement_over_time files.
-#
-"""
-Usage:
-CompareMovementOverTime.py 'test .eddy_movement_over_time' 'ref .eddy_movement_over_time' max_allowed_max_trans_diff max_allowed_99percentile_trans_diff max_allowed_mean_trans_diff max_allowed_max_rot_diff max_allowed_99percentile_rot_diff max_allowed_mean_rot_diff
-or 
-CompareMovementOverTime.py 'test .eddy_movement_over_time' 'ref .eddy_movement_over_time' max_allowed_max_trans_diff max_allowed_99percentile_trans_diff max_allowed_mean_trans_diff max_allowed_max_rot_diff max_allowed_99percentile_rot_diff max_allowed_mean_rot_diff number_of_slices_in_data first_slice_to_check last_slice_to_check
-"""
-
-from __future__ import print_function # I have no words...
-
-import sys
-import CompareTextOutputDef as ct
-
-if __name__ == '__main__':
-    try:
-        if len(sys.argv) == 1:
-            print(__doc__)
-            exit(1)
-
-        if len(sys.argv) == 10:
-            are_the_same = ct.MovementOverTimeAreSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),float(sys.argv[7]),float(sys.argv[8]),int(sys.argv[9]))
-        elif len(sys.argv) == 13:
-            are_the_same = ct.MovementOverTimeAreSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),float(sys.argv[7]),float(sys.argv[8]),int(sys.argv[9]),int(sys.argv[10]),int(sys.argv[11]),int(sys.argv[12]))
-            print('I am here and are_the_same = {}'.format(are_the_same))
-        else:
-            are_the_same = False
-            print('CompareMovementOverTime: Invalid number of input arguments: {}'.format(len(sys.argv)),file=sys.stderr)
-
-        
-        print('I am here and now are_the_same = {}'.format(are_the_same))
-        if are_the_same == True:
-            print('I am about to exit with code 0')
-            sys.exit(0)
-        else:
-            sys.exit(1)
-    except:
-        # raise
-        print('CompareMovementOverTime: Unknown error1',file=sys.stderr)
-        sys.exit(1)
-
-    print('CompareMovementOverTime: Unknown error2',file=sys.stderr)
-    sys.exit(1)
-

EddyS2VTest/CompareOutlierNStdevMaps.py

-#!/usr/bin/env fslpython
-# CompareOutlierNStdevMaps.py - script for comparing two .eddy_outlier_n_stdev_map files.
-#
-"""
-Usage:
-CompareOutlierNStdevMaps.py 'test .eddy_outlier_n_stdev_map' 'ref .eddy_outlier_n_stdev_map' max_allowed_max_relative_z_score_diff_for_z>1 max_allowed_mean_relative_z_score_diff_for_z>1 max_allowed_max_relative_z_score_diff_for_z>2 max_allowed_mean_relative_z_score_diff_for_z>2
-or 
-CompareOutlierNStdevMaps.py 'test .eddy_outlier_n_stdev_map' 'ref .eddy_outlier_n_stdev_map' max_allowed_max_relative_z_score_diff_for_z>1 max_allowed_mean_relative_z_score_diff_for_z>1 max_allowed_max_relative_z_score_diff_for_z>2 max_allowed_mean_relative_z_score_diff_for_z>2 number_of_slices_in_data first_slice_to_check last_slice_to_check
-"""
-
-from __future__ import print_function # I have no words...
-
-import sys
-import CompareTextOutputDef as ct
-
-if __name__ == '__main__':
-    try:
-        if len(sys.argv) == 1:
-            print(__doc__)
-            are_the_same = false
-        if len(sys.argv) == 7:
-            are_the_same = ct.OutlierNStdevMapsAreTheSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]))
-        elif len(sys.argv) == 10:
-            are_the_same = ct.OutlierNStdevMapsAreTheSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),int(sys.argv[7]),int(sys.argv[8]),int(sys.argv[9]))
-        else:
-            print('CompareOutlierNStdevMaps: Invalid number of input arguments: {}'.format(len(sys.argv)-1),file=sys.stderr)
-        if are_the_same == True:
-            exit(0)
-        else:
-            exit(1)
-    except Exception as e:
-        print('CompareOutlierNStdevMaps: Unknown error: {}'.format(e.args[0]),file=sys.stderr)
-        exit(1)
-
-    print('CompareOutlierNStdevMaps: Unknown error',file=sys.stderr)
-    exit(1)
-

EddyS2VTest/CompareOutlierNStdevMaps.py~

-#!/usr/bin/env fslpython
-# CompareOutlierNStdevMaps.py - script for comparing two .eddy_outlier_n_stdev_map files.
-#
-"""
-Usage:
-CompareOutlierNStdevMaps.py 'test .eddy_outlier_n_stdev_map' 'ref .eddy_outlier_n_stdev_map' max_allowed_max_relative_z_score_diff_for_z>1 max_allowed_mean_relative_z_score_diff_for_z>1 max_allowed_max_relative_z_score_diff_for_z>2 max_allowed_mean_relative_z_score_diff_for_z>2
-or 
-CompareOutlierNStdevMaps.py 'test .eddy_outlier_n_stdev_map' 'ref .eddy_outlier_n_stdev_map' max_allowed_max_relative_z_score_diff_for_z>1 max_allowed_mean_relative_z_score_diff_for_z>1 max_allowed_max_relative_z_score_diff_for_z>2 max_allowed_mean_relative_z_score_diff_for_z>2 number_of_slices_in_data first_slice_to_check last_slice_to_check
-"""
-
-from __future__ import print_function # I have no words...
-
-import sys
-import CompareTextOutputDef as ct
-
-if __name__ == '__main__':
-    try:
-        if len(sys.argv) == 1:
-            print(__doc__)
-            are_the_same = false
-        if len(sys.argv) == 7:
-            are_the_same = ct.OutlierNStdevMapsAreTheSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]))
-        elif len(sys.argv) == 10:
-            are_the_same = ct.OutlierNStdevMapsAreTheSame(sys.argv[1],sys.argv[2],float(sys.argv[3]),float(sys.argv[4]),float(sys.argv[5]),float(sys.argv[6]),int(sys.argv[7]),int(sys.argv[8]),int(sys.argv[9]))
-        else:
-            print('CompareOutlierNStdevMaps: Invalid number of input arguments: {}'.format(len(sys.argv)),file=sys.stderr)
-        if are_the_same == True:
-            exit(0)
-        else:
-            exit(1)
-    except:
-        print('CompareOutlierNStdevMaps: Unknown error',file=sys.stderr)
-        exit(1)
-
-    print('CompareOutlierNStdevMaps: Unknown error',file=sys.stderr)
-    exit(1)
-

EddyS2VTest/CompareTextOutputDef.py

-#!/usr/bin/env python
-# CompareTextOutputDef.py - contains defs for comparing textual output from eddy
-# It has code for comparing two .eddy_movement_over_time files and for comparing
-# two .eddy_outlier_n_stdev_map files.
-#
-
-from __future__ import print_function # I have no words...
-
-import            collections
-import            datetime
-import            tempfile
-import            logging
-import            shutil
-import            fnmatch
-import            re
-import            sys
-import            os
-import os.path as op
-import numpy   as np
-
-def MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip, *args):
-    """Compares two .eddy_movement_over_time files and returns true or false.
-    Pass or fail depends on observed differences compared to limits specified by
-    inputs maxtrans, trans99, meantrans, maxrot, rot99 and meanrot. Skip is 0 or 1 if 
-    we are to skip comparison for x-trans (0) or y-trans (1). If not, skip should be -1.
-    USE: MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip)
-    or
-    MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip, nsl, fsl, lsl)
-    """
-
-    if (skip < -1) or (skip > 1):
-        print('MovementOverTimeAreSame: Invalid skip value',file=sys.stderr)
-        raise ValueError('MovementOverTimeAreSame: Invalid skip value')
-
-    if len(args) != 0 and len(args) != 3:
-        print('MovementOverTimeAreSame: Invalid number of optional arguments',file=sys.stderr)
-        raise RuntimeError('MovementOverTimeAreSame: Invalid number of optional arguments')
-
-    if len(args) == 3:
-        nsl = args[0]
-        fsl = args[1]
-        lsl = args[2]
-        if fsl < 0 or lsl > (nsl-1):
-            print('MovementOverTimeAreSame: Invalid fsl or lsl',file=sys.stderr)
-            raise ValueError('MovementOverTimeAreSame: Invalid fsl or lsl')
-
-    try:
-        test_data = np.genfromtxt(testfname)
-        ref_data = np.genfromtxt(reffname)
-    except Exception as e:
-        print('MovementOverTimeAreSame: File read error: {}'.format(e.args[0]),file=sys.stderr)
-        raise e
-
-    if test_data.shape != ref_data.shape or test_data.shape[1] != 6 or test_data.shape[0] % nsl != 0:
-        print('MovementOverTimeAreSame: Shape mismatch between files {} and {}'.format(testfname,reffname),file=sys.stderr)
-        raise ValueError('MovementOverTimeAreSame: Shape mismatch between files {} and {}'.format(testfname,reffname))
-
-    # Remove end slices if nsl, fsl and lsl have been set
-    if len(args) == 3:
-        try:
-            if test_data.shape[0] % nsl:
-                raise ValueError('')
-            for p in range(0,6):
-                tmp = test_data[:,p]
-                tmp = np.reshape(tmp,(nsl,test_data.shape[0]//nsl),'F') 
-                tmp = tmp[fsl:lsl+1,:]
-                if p == 0:
-                    tmp_test_data = np.reshape(tmp,(tmp.size,1),'F')
-                else:
-                    tmp_test_data = np.concatenate((tmp_test_data,np.reshape(tmp,(tmp.size,1),'F')),axis=1)
-                tmp = ref_data[:,p]
-                tmp = np.reshape(tmp,(nsl,test_data.shape[0]//nsl),'F') 
-                tmp = tmp[fsl:lsl+1,:]
-                if p == 0:
-                    tmp_ref_data = np.reshape(tmp,(tmp.size,1),'F')
-                else:
-                    tmp_ref_data = np.concatenate((tmp_ref_data,np.reshape(tmp,(tmp.size,1),'F')),axis=1)
-        except Exception as e:
-            print('MovementOverTimeAreSame: Input nsl={}, fsl={} or lsl={} not suited for files {} and {}'.format(nsl,fsl,lsl,testfname,reffname),file=sys.stderr)
-            raise ValueError('MovementOverTimeAreSame: Input nsl={}, fsl={} or lsl={} not suited for files {} and {}'.format(nsl,fsl,lsl,testfname,reffname))
-
-
-
-    test_data = tmp_test_data
-    ref_data = tmp_ref_data
-
-
-    diff = test_data - ref_data
-    absdiff = np.absolute(diff)
-    if skip == 0:
-        atdiff = absdiff[:,[1, 2]]
-    elif skip == 1:
-        atdiff = absdiff[:,[0, 2]]
-    elif skip == -1:
-        atdiff = absdiff[:,[0, 1, 2]]
-    
-    ardiff = 180.0 * absdiff[:,[3, 4, 5]] / 3.141593    
-    
-    omaxtrans = np.amax(atdiff)                                    # Max translation
-    omaxrot = np.amax(ardiff)                                      # Max rotation
-    omeantrans = np.max(np.mean(atdiff,0))                         # Mean translation along worst direction
-    omeanrot = np.max(np.mean(ardiff,0))                           # Mean rotation around worst axis
-    otrans99 = np.max(np.sort(atdiff,0)[int(0.99*len(atdiff)),:])  # 99th percentile along worst direction
-    orot99 = np.max(np.sort(ardiff,0)[int(0.99*len(atdiff)),:])    # 99th percentile around worst axis
-
-    print('MovementOverTimeAreSame: maxtrans = {}, maxrot = {}'.format(maxtrans,maxrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: omaxtrans = {}, omaxrot = {}'.format(omaxtrans,omaxrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: meantrans = {}, meanrot = {}'.format(meantrans,meanrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: omeantrans = {}, omeanrot = {}'.format(omeantrans,omeanrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: trans99 = {}, rot99 = {}'.format(trans99,rot99),file=sys.stdout)
-    print('MovementOverTimeAreSame: otrans99 = {}, orot99 = {}'.format(otrans99,orot99),file=sys.stdout)
-
-    rval = True
-    if (omaxtrans > maxtrans) or (omaxrot > maxrot) or (otrans99 > trans99) or (orot99 > rot99) or (omeantrans > meantrans) or (omeanrot > meanrot):
-        rval = False
-    
-    return rval
-
-def OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2, *args):
-    """Compares two .eddy_outlier_n_stdev_map files and returns true or false.
-    Pass or fail depends on observed differences compared to limits specified by
-    inputs maxdiff1, meandiff1, maxdiff2 and meandiff2.
-    USE: OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2)
-    or
-    OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2, nsl, fsl, lsl)
-    """
-
-    if len(args) != 0 and len(args) != 3:
-        print('OutlierNStdevMapsAreTheSame: Invalid number of optional arguments',file=sys.stderr)
-        raise ValueError('OutlierNStdevMapsAreTheSame: Invalid number of optional arguments')
-    if len(args) == 3:
-        nsl = args[0]
-        fsl = args[1]
-        lsl = args[2]
-        if fsl < 0 or lsl > (nsl-1):
-            print('OutlierNStdevMapsAreTheSame: Invalid fsl or lsl',file=sys.stderr)
-            raise ValueError('OutlierNStdevMapsAreTheSame: Invalid fsl or lsl')
-
-    try:
-        test_data = np.genfromtxt(testfname,skip_header=1)
-        ref_data = np.genfromtxt(reffname,skip_header=1)
-    except Exception as e:
-        print('OutlierNStdevMapsAreTheSame: File read error: {}'.format(e.args[0]),file=sys.stderr)
-        raise e
-
-    if test_data.shape != ref_data.shape:
-        print('OutlierNStdevMapsAreTheSame: Shape mismatch between files {} and {}'.format(testfname,reffname),file=sys.stderr)
-        raise ValueError('OutlierNStdevMapsAreTheSame: Shape mismatch between files {} and {}'.format(testfname,reffname))
-
-    if len(args) == 3:
-        if nsl != ref_data.shape[1]:
-            print('OutlierNStdevMapsAreTheSame: Shape mismatch between the files {} and {} and nsl={}'.format(testfname,reffname,nsl),file=sys.stderr)
-            raise ValueError('OutlierNStdevMapsAreTheSame: Shape mismatch between the files {} and {} and nsl={}'.format(testfname,reffname,nsl))
-        ref_data = ref_data[:,fsl:lsl+1]
-        test_data = test_data[:,fsl:lsl+1]
-
-    try:
-        # Find vol-slice tuples that are non-zero in both runs
-        nztest = test_data != 0
-        nzref = ref_data != 0
-        if float(np.sum(nzref)-np.sum(nztest)) / np.sum(nzref) > 0.01: # Ensure ~ equal no. of non-zero
-            print('OutlierNStdevMapsAreTheSame: Large mismatch in number of slices considered',file=sys.stdout)
-            return False
-
-        mask = np.multiply(nztest.astype('int'),nzref.astype('int'))
-        indx = np.where(mask > 0)
-        test_data = test_data[indx]
-        ref_data = ref_data[indx]
-
-        # Calculate difference-ratio for all non-zero tuples
-        radiff = np.divide(np.absolute(test_data - ref_data),np.absolute(ref_data))
-
-        # Check for max and mean diff-ratio for voxels with z > 1
-        indx2 = np.where(ref_data > 1.0)
-        omaxdiff1 = np.max(radiff[indx2])
-        omeandiff1 = np.mean(radiff[indx2])
-
-        # Check for max and mean diff-ratio for voxels with z > 2    
-        indx2 = np.where(ref_data > 2.0)
-        omaxdiff2 = np.max(radiff[indx2])
-        omeandiff2 = np.mean(radiff[indx2])
-
-        print('OutlierNStdevMapsAreTheSame: maxdiff1 = {}, omaxdiff1 = {}'.format(maxdiff1,omaxdiff1),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: meandiff1 = {}, omeandiff1 = {}'.format(meandiff1,omeandiff1),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: maxdiff2 = {}, omaxdiff2 = {}'.format(maxdiff2,omaxdiff2),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: meandiff2 = {}, omeandiff2 = {}'.format(meandiff2,omeandiff2),file=sys.stdout)
-
-    except Exception as e:
-        print('OutlierNStdevMapsAreTheSame: Unknown error: {}'.format(e.args[0]),file=sys.stderr)
-        raise e
-
-
-    rval = True
-    if (omaxdiff1 > maxdiff1) or (omeandiff1 > meandiff1) or (omaxdiff2 > maxdiff2) or (omeandiff2 > meandiff2):
-        rval = False
-
-    return rval

EddyS2VTest/CompareTextOutputDef.py~

-#!/usr/bin/env python
-# CompareTextOutputDef.py - contains defs for comparing textual output from eddy
-# It has code for comparing two .eddy_movement_over_time files and for comparing
-# two .eddy_outlier_n_stdev_map files.
-#
-
-from __future__ import print_function # I have no words...
-
-import            collections
-import            datetime
-import            tempfile
-import            logging
-import            shutil
-import            fnmatch
-import            re
-import            sys
-import            os
-import os.path as op
-import numpy   as np
-
-def MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip, *args):
-    """Compares two .eddy_movement_over_time files and returns true or false.
-    Pass or fail depends on observed differences compared to limits specified by
-    inputs maxtrans, trans99, meantrans, maxrot, rot99 and meanrot. Skip is 0 or 1 if 
-    we are to skip comparison for x-trans (0) or y-trans (1). If not, skip should be -1.
-    USE: MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip)
-    or
-    MovementOverTimeAreSame(testfname, reffname, maxtrans, trans99, meantrans, maxrot, rot99, meanrot, skip, nsl, fsl, lsl)
-    """
-
-    if (skip < -1) or (skip > 1):
-        print('MovementOverTimeAreSame: Invalid skip value',file=sys.stderr)
-        return False
-
-    if len(args) != 0 and len(args) != 3:
-        print('MovementOverTimeAreSame: Invalid number of optional arguments',file=sys.stderr)
-        return False
-
-    if len(args) == 3:
-        nsl = args[0]
-        fsl = args[1]
-        lsl = args[2]
-        if fsl < 0 or lsl > (nsl-1):
-            print('MovementOverTimeAreSame: Invalid fsl or lsl',file=sys.stderr)
-            return False
-
-    try:
-        test_data = np.genfromtxt(testfname)
-        ref_data = np.genfromtxt(reffname)
-    except:
-        print('MovementOverTimeAreSame: Unable to read file {} or {}'.format(testfname,reffname),file=sys.stderr)
-        return False
-
-    if test_data.shape != ref_data.shape or test_data.shape[1] != 6 or test_data.shape[0] % nsl != 0:
-        print('MovementOverTimeAreSame: Shape mismatch between files {} and {}'.format(testfname,reffname),file=sys.stderr)
-        return False
-
-    # Remove end slices if nsl, fsl and lsl have been set
-    if len(args) == 3:
-        try:
-            for p in range(0,6):
-                tmp = test_data[:,p]
-                tmp = np.reshape(tmp,(nsl,test_data.shape[0]/nsl),'F') 
-                tmp = tmp[fsl:lsl+1,:]
-                if p == 0:
-                    tmp_test_data = np.reshape(tmp,(tmp.size,1),'F')
-                else:
-                    tmp_test_data = np.concatenate((tmp_test_data,np.reshape(tmp,(tmp.size,1),'F')),axis=1)
-                tmp = ref_data[:,p]
-                tmp = np.reshape(tmp,(nsl,test_data.shape[0]/nsl),'F') 
-                tmp = tmp[fsl:lsl+1,:]
-                if p == 0:
-                    tmp_ref_data = np.reshape(tmp,(tmp.size,1),'F')
-                else:
-                    tmp_ref_data = np.concatenate((tmp_ref_data,np.reshape(tmp,(tmp.size,1),'F')),axis=1)
-        except:
-            print('MovementOverTimeAreSame: Input nsl, fsl or lsl not suited for files {} and {}'.format(testfname,reffname),file=sys.stderr)
-            return False
-
-
-
-    test_data = tmp_test_data
-    ref_data = tmp_ref_data
-
-
-    diff = test_data - ref_data
-    absdiff = np.absolute(diff)
-    if skip == 0:
-        atdiff = absdiff[:,[1, 2]]
-    elif skip == 1:
-        atdiff = absdiff[:,[0, 2]]
-    elif skip == -1:
-        atdiff = absdiff[:,[0, 1, 2]]
-    
-    ardiff = 180.0 * absdiff[:,[3, 4, 5]] / 3.141593    
-    
-    omaxtrans = np.amax(atdiff)                                    # Max translation
-    omaxrot = np.amax(ardiff)                                      # Max rotation
-    omeantrans = np.max(np.mean(atdiff,0))                         # Mean translation along worst direction
-    omeanrot = np.max(np.mean(ardiff,0))                           # Mean rotation around worst axis
-    otrans99 = np.max(np.sort(atdiff,0)[int(0.99*len(atdiff)),:])  # 99th percentile along worst direction
-    orot99 = np.max(np.sort(ardiff,0)[int(0.99*len(atdiff)),:])    # 99th percentile around worst axis
-
-    print('MovementOverTimeAreSame: maxtrans = {}, maxrot = {}'.format(maxtrans,maxrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: omaxtrans = {}, omaxrot = {}'.format(omaxtrans,omaxrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: meantrans = {}, meanrot = {}'.format(meantrans,meanrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: omeantrans = {}, omeanrot = {}'.format(omeantrans,omeanrot),file=sys.stdout)
-    print('MovementOverTimeAreSame: trans99 = {}, rot99 = {}'.format(trans99,rot99),file=sys.stdout)
-    print('MovementOverTimeAreSame: otrans99 = {}, orot99 = {}'.format(otrans99,orot99),file=sys.stdout)
-
-    rval = True
-    if (omaxtrans > maxtrans) or (omaxrot > maxrot) or (otrans99 > trans99) or (orot99 > rot99) or (omeantrans > meantrans) or (omeanrot > meanrot):
-        rval = False
-    
-    print('I am about to return')
-    return rval
-
-def OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2, *args):
-    """Compares two .eddy_outlier_n_stdev_map files and returns true or false.
-    Pass or fail depends on observed differences compared to limits specified by
-    inputs maxdiff1, meandiff1, maxdiff2 and meandiff2.
-    USE: OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2)
-    or
-    OutlierNStdevMapsAreTheSame(testfname, reffname, maxdiff1, meandiff1, maxdiff2, meandiff2, nsl, fsl, lsl)
-    """
-
-    if len(args) != 0 and len(args) != 3:
-        print('OutlierNStdevMapsAreTheSame: Invalid number of optional arguments',file=sys.stderr)
-        return False
-    if len(args) == 3:
-        nsl = args[0]
-        fsl = args[1]
-        lsl = args[2]
-        if fsl < 0 or lsl > (nsl-1):
-            print('OutlierNStdevMapsAreTheSame: Invalid fsl or lsl',file=sys.stderr)
-            return False
-
-    try:
-        test_data = np.genfromtxt(testfname,skip_header=1)
-        ref_data = np.genfromtxt(reffname,skip_header=1)
-    except:
-        print('OutlierNStdevMapsAreTheSame: Unable to read file {} or {}'.format(testfname,reffname),file=sys.stderr)
-        return False
-
-    if test_data.shape != ref_data.shape:
-        print('OutlierNStdevMapsAreTheSame: Shape mismatch between files {} and {}'.format(testfname,reffname),file=sys.stderr)
-        return False
-
-    if len(args) == 3:
-        if nsl != ref_data.shape[1]:
-            print('OutlierNStdevMapsAreTheSame: Shape mismatch between the files {} and {} and nsl'.format(testfname,reffname),file=sys.stderr)
-            return False
-        ref_data = ref_data[:,fsl:lsl+1]
-        test_data = test_data[:,fsl:lsl+1]
-
-    try:
-        # Find vol-slice tuples that are non-zero in both runs
-        nztest = test_data != 0
-        nzref = ref_data != 0
-        if float(np.sum(nzref)-np.sum(nztest)) / np.sum(nzref) > 0.01: # Ensure ~ equal no. of non-zero
-            print('OutlierNStdevMapsAreTheSame: Large mismatch in number of slices considered',file=sys.stdout)
-            return False
-
-        mask = np.multiply(nztest.astype('int'),nzref.astype('int'))
-        indx = np.where(mask > 0)
-        test_data = test_data[indx]
-        ref_data = ref_data[indx]
-
-        # Calculate difference-ratio for all non-zero tuples
-        radiff = np.divide(np.absolute(test_data - ref_data),np.absolute(ref_data))
-
-        # Check for max and mean diff-ratio for voxels with z > 1
-        indx2 = np.where(ref_data > 1.0)
-        omaxdiff1 = np.max(radiff[indx2])
-        omeandiff1 = np.mean(radiff[indx2])
-
-        # Check for max and mean diff-ratio for voxels with z > 2    
-        indx2 = np.where(ref_data > 2.0)
-        omaxdiff2 = np.max(radiff[indx2])
-        omeandiff2 = np.mean(radiff[indx2])
-
-        print('OutlierNStdevMapsAreTheSame: maxdiff1 = {}, omaxdiff1 = {}'.format(maxdiff1,omaxdiff1),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: meandiff1 = {}, omeandiff1 = {}'.format(meandiff1,omeandiff1),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: maxdiff2 = {}, omaxdiff2 = {}'.format(maxdiff2,omaxdiff2),file=sys.stdout)
-        print('OutlierNStdevMapsAreTheSame: meandiff2 = {}, omeandiff2 = {}'.format(meandiff2,omeandiff2),file=sys.stdout)
-
-    except:
-        print('OutlierNStdevMapsAreTheSame: Unknown error',file=sys.stderr)
-        return False
-
-
-    rval = True
-    if (omaxdiff1 > maxdiff1) or (omeandiff1 > meandiff1) or (omaxdiff2 > maxdiff2) or (omeandiff2 > meandiff2):
-        rval = False
-
-    return rval