From aae7f19fae42ec5c4f8a1d5f3154dc0f4285933e Mon Sep 17 00:00:00 2001
From: Matthew Webster <mwebster@jalapeno.fmrib.ox.ac.uk>
Date: Fri, 16 Aug 2019 14:00:23 +0100
Subject: [PATCH] 6.0.2 changes
6.0.2 changes
6.0.2 changes
---
EddyHigh_b_Test/EddyHigh_b_Feeds.py | 134 ++++++
EddyHigh_b_Test/feedsInputs | 1 +
EddyHigh_b_Test/feedsRun | 128 ++++++
EddyHigh_b_Test/nohup.sh | 2 +
EddyLSRTest/EddyLSRFeeds.py | 133 ++++++
EddyLSRTest/LaunchEddy.sh | 25 ++
EddyLSRTest/feedsInputs | 1 +
EddyLSRTest/feedsRun | 125 ++++++
EddyLSRTest/nohup.sh | 2 +
EddyMBSTest/MoveBySuscFeeds.py | 166 ++++++++
EddyMBSTest/feedsInputs | 1 +
EddyMBSTest/feedsRun | 131 ++++++
EddyMBSTest/nohup.sh | 2 +
EddyS2VTest/CompareMovementOverTime.py | 41 --
EddyS2VTest/CompareMovementOverTime.py~ | 45 ---
EddyS2VTest/CompareOutlierNStdevMaps.py | 37 --
EddyS2VTest/CompareOutlierNStdevMaps.py~ | 37 --
EddyS2VTest/CompareTextOutputDef.py | 197 ---------
EddyS2VTest/CompareTextOutputDef.py~ | 196 ---------
EddyS2VTest/S2VFeeds.py | 381 ++++++++++++++++++
.../CompareTextOutputDef.cpython-35.pyc | Bin 7055 -> 0 bytes
EddyS2VTest/feedsInputs~ | 1 -
EddyS2VTest/feedsRun | 197 ++++-----
EddyS2VTest/feedsRun~ | 190 ---------
EddyS2VTest/nohup.sh | 2 +
EddyS2VTest/rubbestad.txt | 12 +
EddyTest/EddyFeeds.py | 278 +++++++++++++
EddyTest/eddy_cuda.e2284102 | 0
EddyTest/eddy_cuda.e2286236 | 48 ---
EddyTest/eddy_cuda.e2286247 | 0
EddyTest/eddy_cuda.e2294981 | 0
EddyTest/eddy_cuda.e2410534 | 0
EddyTest/eddy_cuda.o2284102 | 88 ----
EddyTest/eddy_cuda.o2286236 | 0
EddyTest/eddy_cuda.o2286247 | 90 -----
EddyTest/eddy_cuda.o2294981 | 90 -----
EddyTest/eddy_cuda.o2410534 | 88 ----
EddyTest/eddy_openmp.e2284103 | 0
EddyTest/eddy_openmp.e2286237 | 48 ---
EddyTest/eddy_openmp.e2286248 | 0
EddyTest/eddy_openmp.e2294982 | 0
EddyTest/eddy_openmp.e2410535 | 0
EddyTest/eddy_openmp.o2284103 | 125 ------
EddyTest/eddy_openmp.o2286237 | 0
EddyTest/eddy_openmp.o2286248 | 127 ------
EddyTest/eddy_openmp.o2294982 | 127 ------
EddyTest/eddy_openmp.o2410535 | 125 ------
EddyTest/eddy_openmp.pe2284103 | 0
EddyTest/eddy_openmp.pe2286237 | 0
EddyTest/eddy_openmp.pe2286248 | 0
EddyTest/eddy_openmp.pe2294982 | 0
EddyTest/eddy_openmp.pe2410535 | 0
EddyTest/eddy_openmp.po2284103 | 0
EddyTest/eddy_openmp.po2286237 | 0
EddyTest/eddy_openmp.po2286248 | 0
EddyTest/eddy_openmp.po2294982 | 0
EddyTest/eddy_openmp.po2410535 | 0
EddyTest/feedsRun | 167 +++-----
EddyTest/feedsRun~ | 171 --------
EddyTest/nohup.sh | 2 +
fsl_course/fdt/bedpostx_gpu/feedsRun | 2 +-
fsl_course/fdt1/test_topup_eddy.sh | 4 +-
fsl_course/fdt2/feedsRun6.sh | 5 +-
fsl_course/feat3/feedsRun3.py | 2 +-
fsl_course/pnm/feedsInputs | 1 +
fsl_course/pnm/feedsRun | 7 +
mist/feedsRun | 1 +
mist/feedsRun~ | 26 --
68 files changed, 1675 insertions(+), 2134 deletions(-)
create mode 100755 EddyHigh_b_Test/EddyHigh_b_Feeds.py
create mode 100755 EddyHigh_b_Test/feedsInputs
create mode 100755 EddyHigh_b_Test/feedsRun
create mode 100755 EddyHigh_b_Test/nohup.sh
create mode 100755 EddyLSRTest/EddyLSRFeeds.py
create mode 100644 EddyLSRTest/LaunchEddy.sh
create mode 100755 EddyLSRTest/feedsInputs
create mode 100755 EddyLSRTest/feedsRun
create mode 100755 EddyLSRTest/nohup.sh
create mode 100755 EddyMBSTest/MoveBySuscFeeds.py
create mode 100755 EddyMBSTest/feedsInputs
create mode 100755 EddyMBSTest/feedsRun
create mode 100755 EddyMBSTest/nohup.sh
delete mode 100755 EddyS2VTest/CompareMovementOverTime.py
delete mode 100755 EddyS2VTest/CompareMovementOverTime.py~
delete mode 100755 EddyS2VTest/CompareOutlierNStdevMaps.py
delete mode 100755 EddyS2VTest/CompareOutlierNStdevMaps.py~
delete mode 100644 EddyS2VTest/CompareTextOutputDef.py
delete mode 100644 EddyS2VTest/CompareTextOutputDef.py~
create mode 100755 EddyS2VTest/S2VFeeds.py
delete mode 100644 EddyS2VTest/__pycache__/CompareTextOutputDef.cpython-35.pyc
delete mode 100755 EddyS2VTest/feedsInputs~
delete mode 100755 EddyS2VTest/feedsRun~
create mode 100755 EddyS2VTest/nohup.sh
create mode 100644 EddyS2VTest/rubbestad.txt
create mode 100755 EddyTest/EddyFeeds.py
delete mode 100755 EddyTest/eddy_cuda.e2284102
delete mode 100755 EddyTest/eddy_cuda.e2286236
delete mode 100755 EddyTest/eddy_cuda.e2286247
delete mode 100755 EddyTest/eddy_cuda.e2294981
delete mode 100755 EddyTest/eddy_cuda.e2410534
delete mode 100755 EddyTest/eddy_cuda.o2284102
delete mode 100755 EddyTest/eddy_cuda.o2286236
delete mode 100755 EddyTest/eddy_cuda.o2286247
delete mode 100755 EddyTest/eddy_cuda.o2294981
delete mode 100755 EddyTest/eddy_cuda.o2410534
delete mode 100755 EddyTest/eddy_openmp.e2284103
delete mode 100755 EddyTest/eddy_openmp.e2286237
delete mode 100755 EddyTest/eddy_openmp.e2286248
delete mode 100755 EddyTest/eddy_openmp.e2294982
delete mode 100755 EddyTest/eddy_openmp.e2410535
delete mode 100755 EddyTest/eddy_openmp.o2284103
delete mode 100755 EddyTest/eddy_openmp.o2286237
delete mode 100755 EddyTest/eddy_openmp.o2286248
delete mode 100755 EddyTest/eddy_openmp.o2294982
delete mode 100755 EddyTest/eddy_openmp.o2410535
delete mode 100755 EddyTest/eddy_openmp.pe2284103
delete mode 100755 EddyTest/eddy_openmp.pe2286237
delete mode 100755 EddyTest/eddy_openmp.pe2286248
delete mode 100755 EddyTest/eddy_openmp.pe2294982
delete mode 100755 EddyTest/eddy_openmp.pe2410535
delete mode 100755 EddyTest/eddy_openmp.po2284103
delete mode 100755 EddyTest/eddy_openmp.po2286237
delete mode 100755 EddyTest/eddy_openmp.po2286248
delete mode 100755 EddyTest/eddy_openmp.po2294982
delete mode 100755 EddyTest/eddy_openmp.po2410535
delete mode 100755 EddyTest/feedsRun~
create mode 100755 EddyTest/nohup.sh
create mode 100644 fsl_course/pnm/feedsInputs
create mode 100755 fsl_course/pnm/feedsRun
delete mode 100755 mist/feedsRun~
diff --git a/EddyHigh_b_Test/EddyHigh_b_Feeds.py b/EddyHigh_b_Test/EddyHigh_b_Feeds.py
new file mode 100755
index 0000000..588dfaf
--- /dev/null
+++ b/EddyHigh_b_Test/EddyHigh_b_Feeds.py
@@ -0,0 +1,134 @@
+#!/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
diff --git a/EddyHigh_b_Test/feedsInputs b/EddyHigh_b_Test/feedsInputs
new file mode 100755
index 0000000..16e9916
--- /dev/null
+++ b/EddyHigh_b_Test/feedsInputs
@@ -0,0 +1 @@
+EddyHigh_b_TestData/eddyData
diff --git a/EddyHigh_b_Test/feedsRun b/EddyHigh_b_Test/feedsRun
new file mode 100755
index 0000000..ee99d99
--- /dev/null
+++ b/EddyHigh_b_Test/feedsRun
@@ -0,0 +1,128 @@
+#! /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
+
+
diff --git a/EddyHigh_b_Test/nohup.sh b/EddyHigh_b_Test/nohup.sh
new file mode 100755
index 0000000..480cedc
--- /dev/null
+++ b/EddyHigh_b_Test/nohup.sh
@@ -0,0 +1,2 @@
+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyHigh_b blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?
diff --git a/EddyLSRTest/EddyLSRFeeds.py b/EddyLSRTest/EddyLSRFeeds.py
new file mode 100755
index 0000000..5fe1589
--- /dev/null
+++ b/EddyLSRTest/EddyLSRFeeds.py
@@ -0,0 +1,133 @@
+#!/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
diff --git a/EddyLSRTest/LaunchEddy.sh b/EddyLSRTest/LaunchEddy.sh
new file mode 100644
index 0000000..3bbd752
--- /dev/null
+++ b/EddyLSRTest/LaunchEddy.sh
@@ -0,0 +1,25 @@
+#! /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
+
+
+
+
diff --git a/EddyLSRTest/feedsInputs b/EddyLSRTest/feedsInputs
new file mode 100755
index 0000000..d2c80fb
--- /dev/null
+++ b/EddyLSRTest/feedsInputs
@@ -0,0 +1 @@
+EddyLSRTestData/eddyData
diff --git a/EddyLSRTest/feedsRun b/EddyLSRTest/feedsRun
new file mode 100755
index 0000000..a2db1ee
--- /dev/null
+++ b/EddyLSRTest/feedsRun
@@ -0,0 +1,125 @@
+#! /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
+
+
diff --git a/EddyLSRTest/nohup.sh b/EddyLSRTest/nohup.sh
new file mode 100755
index 0000000..447f22f
--- /dev/null
+++ b/EddyLSRTest/nohup.sh
@@ -0,0 +1,2 @@
+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyLSRTest blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?
diff --git a/EddyMBSTest/MoveBySuscFeeds.py b/EddyMBSTest/MoveBySuscFeeds.py
new file mode 100755
index 0000000..0c60e72
--- /dev/null
+++ b/EddyMBSTest/MoveBySuscFeeds.py
@@ -0,0 +1,166 @@
+#!/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)
diff --git a/EddyMBSTest/feedsInputs b/EddyMBSTest/feedsInputs
new file mode 100755
index 0000000..0d13a25
--- /dev/null
+++ b/EddyMBSTest/feedsInputs
@@ -0,0 +1 @@
+EddyMBSTestData/eddyData
diff --git a/EddyMBSTest/feedsRun b/EddyMBSTest/feedsRun
new file mode 100755
index 0000000..063b610
--- /dev/null
+++ b/EddyMBSTest/feedsRun
@@ -0,0 +1,131 @@
+#! /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
+
diff --git a/EddyMBSTest/nohup.sh b/EddyMBSTest/nohup.sh
new file mode 100755
index 0000000..c4985c1
--- /dev/null
+++ b/EddyMBSTest/nohup.sh
@@ -0,0 +1,2 @@
+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyMBSTest blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?
diff --git a/EddyS2VTest/CompareMovementOverTime.py b/EddyS2VTest/CompareMovementOverTime.py
deleted file mode 100755
index af50852..0000000
--- a/EddyS2VTest/CompareMovementOverTime.py
+++ /dev/null
@@ -1,41 +0,0 @@
-#!/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)
-
diff --git a/EddyS2VTest/CompareMovementOverTime.py~ b/EddyS2VTest/CompareMovementOverTime.py~
deleted file mode 100755
index e01935a..0000000
--- a/EddyS2VTest/CompareMovementOverTime.py~
+++ /dev/null
@@ -1,45 +0,0 @@
-#!/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)
-
diff --git a/EddyS2VTest/CompareOutlierNStdevMaps.py b/EddyS2VTest/CompareOutlierNStdevMaps.py
deleted file mode 100755
index e5feb39..0000000
--- a/EddyS2VTest/CompareOutlierNStdevMaps.py
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/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)
-
diff --git a/EddyS2VTest/CompareOutlierNStdevMaps.py~ b/EddyS2VTest/CompareOutlierNStdevMaps.py~
deleted file mode 100755
index ee7fed3..0000000
--- a/EddyS2VTest/CompareOutlierNStdevMaps.py~
+++ /dev/null
@@ -1,37 +0,0 @@
-#!/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)
-
diff --git a/EddyS2VTest/CompareTextOutputDef.py b/EddyS2VTest/CompareTextOutputDef.py
deleted file mode 100644
index 0b3cb07..0000000
--- a/EddyS2VTest/CompareTextOutputDef.py
+++ /dev/null
@@ -1,197 +0,0 @@
-#!/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
diff --git a/EddyS2VTest/CompareTextOutputDef.py~ b/EddyS2VTest/CompareTextOutputDef.py~
deleted file mode 100644
index f1b952c..0000000
--- a/EddyS2VTest/CompareTextOutputDef.py~
+++ /dev/null
@@ -1,196 +0,0 @@
-#!/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
diff --git a/EddyS2VTest/S2VFeeds.py b/EddyS2VTest/S2VFeeds.py
new file mode 100755
index 0000000..a51905b
--- /dev/null
+++ b/EddyS2VTest/S2VFeeds.py
@@ -0,0 +1,381 @@
+#!/usr/bin/env fslpython
+
+import sys
+import numpy as np
+import nibabel as nib
+import datetime
+import pandas as pd
+
+class S2VFeedsType(object):
+ """ The purpose of this class is to make all the comparisons between
+ estimated and ground truth eddy results for the Slice-to-vol
+ feeds for eddy. There are some pretty hard assumptions of the
+ data in here, so will need revision if we change the data for the test.
+ """
+
+ def __init__(self,mask,bvals,corr,truth,est_MOT,true_MOT,est_ol,true_ol):
+ # 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._truth = nib.load(truth,mmap=False)
+ if not (all(self._mask.header['dim'][1:4] == self._corr.header['dim'][1:4]) and
+ all(self._mask.header['dim'][1:4] == self._truth.header['dim'][1:4])):
+ raise Exception('S2VFeedsType:__init__:Size mismatch in first three dimensions of images')
+ if not (self._mask.header['dim'][4] == 1 and
+ self._corr.header['dim'][4] == self._truth.header['dim'][4]):
+ raise Exception('S2VFeedsType:__init__:Size mismatch in fourth dimension of images')
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error opening image files')
+
+ # Read bvals file
+ try:
+ self._bvals = pd.read_csv(bvals, delim_whitespace=True, header=None, skiprows=0).values[0]
+ if self._bvals.shape[0] != self._corr.header['dim'][4]:
+ raise Exception('S2VFeedsType:__init__:Size mismatch between bvals and images')
+ b0index = (self._bvals < 100).nonzero()[0]
+ lowbindex = np.logical_and([self._bvals > 100], [self._bvals < 1000]).nonzero()[1]
+ hibindex = (self._bvals > 1000).nonzero()[0]
+ self._indicies = [b0index, lowbindex, hibindex]
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error opening bvals text file')
+
+ # Read outlier text files
+ try:
+ self._est_ol = pd.read_csv(est_ol, delim_whitespace=True, header=None, skiprows=1).values
+ self._true_ol = pd.read_csv(true_ol, sep=',', header=None, decimal='.').values
+ if self._est_ol.shape[0] != self._corr.header['dim'][4] or \
+ self._true_ol.shape[0] != self._corr.header['dim'][4]:
+ raise Exception('S2VFeedsType:__init__:Size mismatch in volume dimension of outlier data')
+ if self._est_ol.shape[1] != self._corr.header['dim'][3] or \
+ self._true_ol.shape[1] != self._corr.header['dim'][3]:
+ raise Exception('S2VFeedsType:__init__:Size mismatch in slice dimension of outlier data')
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error opening outlier text files')
+
+ # Read movement-over-time files
+ try:
+ self._est_MOT = pd.read_csv(est_MOT, delim_whitespace=True, header=None, skiprows=0).values
+ self._true_MOT = pd.read_csv(true_MOT, delim_whitespace=True, header=None, skiprows=0).values
+ if self._est_MOT.shape[1] != 6 or \
+ self._true_MOT.shape[1] != 6:
+ raise Exception('S2VFeedsType:__init__:Size mismatch in parameter dimension of MOT data')
+ if self._est_MOT.shape[0] != np.prod(self._corr.header['dim'][3:5]) or \
+ self._true_MOT.shape[0] != np.prod(self._corr.header['dim'][3:5]):
+ raise Exception('S2VFeedsType:__init__:Size mismatch in time dimension of MOT data')
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error opening movement over time text files')
+
+ # 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._truth.get_data().astype(float))
+
+ self._corrdiffmeans = [np.zeros(len(self._indicies[0])), np.zeros(len(self._indicies[1])), np.zeros(len(self._indicies[2]))]
+ for b in range(0,len(self._corrdiffmeans)):
+ cntr = 0
+ for vol in self._indicies[b]:
+ tmpdiff = np.multiply(mask, corrdiff[:, :, :, vol])
+ self._corrdiffmeans[b][cntr] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+ cntr += 1
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error calculating image statistics')
+
+ # Look extra carefully at the volumes that we know are affected by S2V movement
+ # N.B. this will no longer be valid if we change the data we use for the tests
+ #try:
+ # b0vols = [9, 27, 36, 45, 54, 63, 90, 99]
+ # lowbvols = [2, 4, 5, 6, 7, 8, 10, 11]
+ # hibvols = [41, 44, 57, 58, 79, 80, 89, 98, 100, 102, 103, 104, 106, 107]
+ # self._s2v_indicies = [b0vols, lowbvols, hibvols]
+ #except Exception as e:
+ # print(str(e))
+ # raise Exception('S2VFeedsType:__init__:Error calculating image statistics for S2V affected volumes')
+
+ # Calculate false positives, false negatives and true positives for outliers
+ try:
+ # Find slices with > 500 brain voxels
+ vsl = []
+ for sl in range(0,self._mask.header['dim'][3]):
+ if mask[:,:,sl].sum() > 500:
+ vsl.append(sl)
+ # Go through and count positives and false negatives
+ self._tol = 0; self._tdol = 0; self._tp = 0; self._fp = 0; self._fn = 0
+ for sl in vsl:
+ for vol in range(0,self._true_ol.shape[0]):
+ if not np.isnan(self._true_ol[vol,sl]): # If it is a true outlier
+ self._tol += 1
+ if self._est_ol[vol,sl] == 1:
+ self._tp += 1
+ if self._true_ol[vol,sl] < 0.8: # If it ought to have been detected
+ self._tdol += 1
+ if self._est_ol[vol,sl] != 1:
+ self._fn += 1
+ else:
+ if self._est_ol[vol,sl] == 1:
+ self._fp += 1
+
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error calculating outlier statistics')
+
+ # Compare movement over time to ground truth
+ try:
+ # Map slices with > 500 brain voxels onto time-series
+ index = []
+ for vol in range(0,self._corr.header['dim'][4]):
+ tmpindx = [x+vol*self._corr.header['dim'][3] for x in vsl]
+ index.append(tmpindx)
+ # Calculate differences in mm and degrees
+ self._MOT_diffs = np.zeros((len(index)*len(index[0]),6))
+ cntr = 0
+ mmdeg = np.array([1, 1, 1, 180/np.pi, 180/np.pi, 180/np.pi])
+ for vol in range(0,len(index)):
+ for i in range(0,len(index[vol])):
+ self._MOT_diffs[cntr,:] = np.abs(mmdeg * (self._est_MOT[index[vol][i],:] - self._true_MOT[index[vol][i],:]))
+ cntr += 1
+ except Exception as e:
+ print(str(e))
+ raise Exception('S2VFeedsType:__init__:Error calculating movement over time statistics')
+
+ def BValues(self):
+ return [0, 700, 2000]
+
+ def MeanDiffsOfCorrectedImages(self):
+ return self._corrdiffmeans
+
+ def MeanDiffsOfCorrectedImages_b0(self):
+ return self._corrdiffmeans[0]
+
+ def MeanDiffsOfCorrectedImages_700(self):
+ return self._corrdiffmeans[1]
+
+ def MeanDiffsOfCorrectedImages_2000(self):
+ return self._corrdiffmeans[2]
+
+ def DiffsOfMovementOverTime(self):
+ return self._MOT_diffs
+
+ def MeanDiffsOfMovementOverTime(self):
+ return self._MOT_diffs.mean(axis=0)
+
+ def MaxDiffsOfMovementOverTime(self):
+ return self._MOT_diffs.max(axis=0)
+
+ def OlFalsePositives(self):
+ return self._fp
+
+ def OlFalseNegatives(self):
+ return self._fn
+
+ def OlTruePositives(self):
+ return self._tp
+
+ def OlTotalNumber(self):
+ return self._tol
+
+ def OlDetectableNumber(self):
+ return self._tdol
+
+def main(argv):
+ # This is the main program that tests the output from running
+ # the S2VFeeds tests.
+ try:
+ if len(argv) != 23:
+ print('S2VFeeds.py usage: S2VFeeds output_dir prefix mask bvals corrected truth est_MOT true_MOT ' \
+ 'est_ol true_ol allowed_mean_diff_b0 allowed_max_diff_b0 allowed_mean_diff_700 ' \
+ 'allowed_max_diff_700 allowed_mean_diff_2000 allowed_max_diff_2000 allowed_false_pos ' \
+ 'allowed_false_neg allowed_mean_trans_error allowed_max_trans_error allowed_mean_rot_error ' \
+ 'allowed_max_rot_error')
+ sys.exit(1)
+ else:
+ output_dir = argv[1]
+ output_prefix = argv[2]
+ mask = argv[3]
+ bvals = argv[4]
+ corrected = argv[5]
+ truth = argv[6]
+ est_MOT = argv[7]
+ true_MOT = argv[8]
+ est_ol = argv[9]
+ true_ol = argv[10]
+ allowed_mean_diff_b0 = float(argv[11])
+ allowed_max_diff_b0 = float(argv[12])
+ allowed_mean_diff_700 = float(argv[13])
+ allowed_max_diff_700 = float(argv[14])
+ allowed_mean_diff_2000 = float(argv[15])
+ allowed_max_diff_2000 = float(argv[16])
+ allowed_false_positives = float(argv[17])
+ allowed_false_negatives = float(argv[18])
+ allowed_mean_trans_error = float(argv[19])
+ allowed_max_trans_error = float(argv[20])
+ allowed_mean_rot_error = float(argv[21])
+ allowed_max_rot_error = float(argv[22])
+
+
+ # Try to create MoveBySuscFeedsType object (involves reading all files)
+ try:
+ s2vf = S2VFeedsType(mask,bvals,corrected,truth,est_MOT,true_MOT,est_ol,true_ol)
+ except Exception as e:
+ print(str(e))
+ print('main: Error when creating S2VFeedsType object.')
+ sys.exit(1)
+
+ # Inspect statistics to see if it passes test
+ try:
+ # Image differences
+ passes_test = True
+ if s2vf.MeanDiffsOfCorrectedImages()[0].mean() > allowed_mean_diff_b0 or \
+ s2vf.MeanDiffsOfCorrectedImages()[0].max() > allowed_max_diff_b0 or \
+ s2vf.MeanDiffsOfCorrectedImages()[1].mean() > allowed_mean_diff_700 or \
+ s2vf.MeanDiffsOfCorrectedImages()[1].max() > allowed_max_diff_700 or \
+ s2vf.MeanDiffsOfCorrectedImages()[2].mean() > allowed_mean_diff_2000 or \
+ s2vf.MeanDiffsOfCorrectedImages()[2].max() > allowed_max_diff_2000:
+ passes_test = False
+ # Outliers
+ if s2vf.OlFalsePositives() > allowed_false_positives or \
+ s2vf.OlFalseNegatives() > allowed_false_negatives:
+ passes_test = False
+ # Movement-over-time estimates
+ if s2vf.MeanDiffsOfMovementOverTime()[0:3].mean() > allowed_mean_trans_error or \
+ s2vf.MeanDiffsOfMovementOverTime()[3:6].mean() > allowed_mean_rot_error or \
+ s2vf.MaxDiffsOfMovementOverTime()[0:3].mean() > allowed_max_trans_error or \
+ s2vf.MaxDiffsOfMovementOverTime()[3:6].mean() > allowed_max_rot_error:
+ passes_test = False
+ except Exception as e:
+ print(str(e))
+ print('main: Failed calculating stats for test.')
+ sys.exit(1)
+
+ # Write report file
+ try:
+ fp = open(output_dir + '/' + output_prefix + '_S2VFeedsReport.txt','w')
+ except Exception as e:
+ print(str(e))
+ print('Cannot open report file: ' + output_dir + '/' + output_prefix + '_S2VFeedsReport.txt')
+ sys.exit(1)
+ else:
+ try:
+ fp.write('S2VFeeds 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 corrected image differences separately for b=0, 700 and 2000
+ fp.write('\nThe absolute differences, averaged across the mask, ' \
+ 'for the corrected b = {0:d} images were: '.format(s2vf.BValues()[0]) + \
+ ' '.join(["{0:.4f}".format(elem) for elem in s2vf.MeanDiffsOfCorrectedImages()[0]]) + '\n')
+ fp.write('That gives mean error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[0].mean()) + \
+ ', and a max error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[0].max()) + '\n')
+ fp.write('The allowed mean of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[0],allowed_mean_diff_b0))
+ fp.write('The allowed maximum of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[0],allowed_max_diff_b0))
+ if s2vf.MeanDiffsOfCorrectedImages()[0].mean() > allowed_mean_diff_b0 or \
+ s2vf.MeanDiffsOfCorrectedImages()[0].max() > allowed_max_diff_b0:
+ 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 corrected b = {0:d} images were: '.format(s2vf.BValues()[1]) + \
+ ' '.join(["{0:.4f}".format(elem) for elem in s2vf.MeanDiffsOfCorrectedImages()[1]]) + '\n')
+ fp.write('That gives mean error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[1].mean()) + \
+ ', and a max error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[1].max()) + '\n')
+ fp.write('The allowed mean of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[1],allowed_mean_diff_700))
+ fp.write('The allowed maximum of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[1],allowed_max_diff_700))
+ if s2vf.MeanDiffsOfCorrectedImages()[1].mean() > allowed_mean_diff_700 or \
+ s2vf.MeanDiffsOfCorrectedImages()[1].max() > allowed_max_diff_700:
+ 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 corrected b = {0:d} images were: '.format(s2vf.BValues()[2]) + \
+ ' '.join(["{0:.4f}".format(elem) for elem in s2vf.MeanDiffsOfCorrectedImages()[2]]) + '\n')
+ fp.write('That gives mean error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[2].mean()) + \
+ ', and a max error: ' + "{0:.4f}".format(s2vf.MeanDiffsOfCorrectedImages()[2].max()) + '\n')
+ fp.write('The allowed mean of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[2],allowed_mean_diff_2000))
+ fp.write('The allowed maximum of averaged differences for the corrected b = {0:d} images was: {1:.4f}\n'\
+ .format(s2vf.BValues()[2],allowed_max_diff_2000))
+ if s2vf.MeanDiffsOfCorrectedImages()[2].mean() > allowed_mean_diff_2000 or \
+ s2vf.MeanDiffsOfCorrectedImages()[2].max() > allowed_max_diff_2000:
+ fp.write('Based on these criteria the test failed\n')
+ else:
+ fp.write('Based on these criteria the test passed\n')
+
+ # Report on differences in outlier estimations
+ fp.write('\nThe input data had ' + '{:d}'.format(s2vf.OlTotalNumber()) + ' true outliers')
+ fp.write('\nOf these ' + '{:d}'.format(s2vf.OlDetectableNumber()) + ' had a signal loss greater ' + \
+ 'than 20% and should be detected')
+ fp.write('\nThe analysis yielded ' + '{:d}'.format(s2vf.OlTruePositives()) + ' true positives, ' + \
+ '{:d}'.format(s2vf.OlFalseNegatives()) + ' false negatives and ' + \
+ '{:d}'.format(s2vf.OlFalsePositives()) + ' false positives\n')
+ if s2vf.OlFalsePositives() > allowed_false_positives or \
+ s2vf.OlFalseNegatives() > allowed_false_negatives:
+ fp.write('Based on these criteria the test failed\n')
+ else:
+ fp.write('Based on these criteria the test passed\n')
+ # Report on differences in Movement over time
+ fp.write('\nThe mean absolute errors in movement estimated over time were')
+ fp.write('\nx-trans: {0:.4f} mm, y-trans: {1:.4f} mm, z-trans: {2:.4f} mm'.\
+ format(s2vf.MeanDiffsOfMovementOverTime()[0],s2vf.MeanDiffsOfMovementOverTime()[1],\
+ s2vf.MeanDiffsOfMovementOverTime()[2]))
+ fp.write('\nx-rot: {0:.4f} deg, y-rot: {1:.4f} deg, z-rot: {2:.4f} deg'.\
+ format(s2vf.MeanDiffsOfMovementOverTime()[3],s2vf.MeanDiffsOfMovementOverTime()[4],\
+ s2vf.MeanDiffsOfMovementOverTime()[5]))
+ fp.write('\nThe allowed mean absolute errors in movement estimated over time were')
+ fp.write('Translation: {0:.4f} mm and Rotation: {1:.4f} deg\n'.format(allowed_mean_trans_error,allowed_mean_rot_error))
+ if s2vf.MeanDiffsOfMovementOverTime()[0:3].mean() > allowed_mean_trans_error or \
+ s2vf.MeanDiffsOfMovementOverTime()[3:6].mean() > allowed_mean_rot_error:
+ fp.write('Based on these criteria the test failed\n')
+ else:
+ fp.write('Based on these criteria the test passed\n')
+
+ fp.write('\nThe max absolute errors in movement estimated over time were')
+ fp.write('\nx-trans: {0:.4f} mm, y-trans: {1:.4f} mm, z-trans: {2:.4f} mm'.\
+ format(s2vf.MaxDiffsOfMovementOverTime()[0],s2vf.MaxDiffsOfMovementOverTime()[1],\
+ s2vf.MaxDiffsOfMovementOverTime()[2]))
+ fp.write('\nx-rot: {0:.4f} deg, y-rot: {1:.4f} deg, z-rot: {2:.4f} deg'.\
+ format(s2vf.MaxDiffsOfMovementOverTime()[3],s2vf.MaxDiffsOfMovementOverTime()[4],\
+ s2vf.MaxDiffsOfMovementOverTime()[5]))
+ fp.write('\nThe allowed mean absolute errors in movement estimated over time were')
+ fp.write('Translation: {0:.4f} mm and Rotation: {1:.4f} deg\n'.format(allowed_max_trans_error,allowed_max_rot_error))
+ if s2vf.MaxDiffsOfMovementOverTime()[0:3].mean() > allowed_max_trans_error or \
+ s2vf.MaxDiffsOfMovementOverTime()[3:6].mean() > allowed_max_rot_error:
+ 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 + '_S2VFeedsReport.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)
\ No newline at end of file
diff --git a/EddyS2VTest/__pycache__/CompareTextOutputDef.cpython-35.pyc b/EddyS2VTest/__pycache__/CompareTextOutputDef.cpython-35.pyc
deleted file mode 100644
index 83de7a88f1ed35b3bf712231e20088e428bb0a6a..0000000000000000000000000000000000000000
GIT binary patch
literal 0
HcmV?d00001
literal 7055
zcmc&(O>o=B72X9vQlh9o>fiDY^3%jh6v<Ya$=FG1XB<0Cr*YcIu@lU6G_U|lkVt?9
zKv|^8q?4qVUUJE?ho0JL55DBwLwo3@hh96lbvild)SlbvwBLIRkd!PjN@voN!t(9D
zeLws5Vevj_c67A(+}pqV;FE6&@mJCJ(Li6v6@7um!JiaWA!<@sQb_7+qN)i?V;*Q-
zSb9bq!ZI>CC#+mX=Y^He=z_2c89gMdp^P3D)^J7_g;mVx5wTMe)@arp6V|wJV9^kt
z3SmtM2UhY-Pm*uzON9eCE|z36Ru4Tt-l#PE&Dab4+ndA0z~=(4Cvm-rEBZ}dh&H@9
zFLu$JXIhF~LxeAjn3oG{u@pO+_(Vz>!!%LiV<}=ow1<g?GrPJ7TYV$BYy>K>_CqSz
z<D)7dpo#EH(g<}6ogz~b)Cf}rP$khGVPo+052B4O!X(8!dHW@Mi-jCA1#gE^Z<|kw
zI4^bzqAh5oGn#cJcgMs%|F$sVA+cK&;XNrk<3hGpSnIIZ85Zqvu>*7PjaJ-m1sp`<
z-UklYMFF`F$Hn`;6z{i3#SXk5qZR>*K9_0)J})>&hM*U#nN36aM@g|Sv0dyBTtve4
z+oXLs%sxJ#(rA8>SxrrZ|KeR7P2;lvnT$)cE_ehhZ~W*(*_jY<R4wPHnrKgmiiYq{
za>PbOdy)d5E=4Ilg9*-<_(T_&MSLt9N~M^Cd{<Gb1^Q^2*zTWev?WU8q}UlFiJqpi
z@Lyz;<Dy;ev;Eji62V<T8!0<w3V7G#H|P8N`XJw@#E(aLaqxZW)V{OSy?egjO+;a+
zM7+(3IKc^l*mHQ1$dnnlosMS43G-Bn9G9qNx*y<ZXN>k00i5qlv)9OkX}&*Z_+(Jq
z_CupR(+e@;?i9q9rjle<_;<5I1UuY5JQZ1XXzKGk^#UiuuO6GK6ziE{2YJ<R7Yi>7
zXG(0T8Pm?RSnSNP+n71aZqH}kIqoiSS0&RF%|Rj}>FGQ;bC@;P$6D-TP4}c?I5Wd%
z9^0CxbAgWl<UFN1X+3*bG}C$-(Q~}zv%KZF33GjC0NKM?j3Z)>**LZrj%Cj?8|QT4
zSoQ+5&xrQovFtO<)<rviEPIjJXGQzmV{GTF=$sdA_@&XaB>F`j4q^!@4H-S7&yo6c
z9G)P<)ci$q8g%%}ymOAVoyXJ7v|q;)`bJQz+o2no@q@rzaUG|*Q499mnv3The#4E}
ztGQ;yt73}nJ7(y{jnK#Ru;H3PXjbfM<gOG=e12d@ks9)<rsLLK-$CCugZq&i?zxWX
zc$JD9y8b3~Z7REt83$(7t9fx`Ms;`7t9Y1rzsU+czuv%j&EAhg+mDuE6hBw4ESojg
z_R}$N!ysNZ(OkL0E)o;H6?1*ptD9bAULvn9n_k6awGUhqD<;DcO;DbA5j?G6-u?@$
z*t~dYiNu;&_wv$;`BugB$rj5++d-oWUtRNs%PIixt>3(AzMV$z9R%+VV)uIJuG=;D
zV(dn7#Rmmf-AYP*jisqY^HRzD5-%kP`R9LRbJ_Hx>atm(rW#e13kMRlPwV*;T)1k!
z<?q>5&rvZ0*>Dqs(kUdl1D%9l`$}QIQM>PkW>7JMI#pnH6<fH~pe&6JsLniv18|!*
z6ppu+Pc(Q7o?(A%$8^Im2(OwCJGT!|*NN7(6J^%7?Ye8$yr^c!o7+f;_<`#pnfRW0
z*x{Uf*hxepZ&`-D??tV<C+d2OQl7&6+QSY<8o#O6L?QdksNu!9M<`d1nAn1!Lj0t1
z{lLKQK@QE<N5B33$Nzj{`1;nPH~-R_e)P9%trt%;o@Jx?8kUWG<2#Ext(Q&|3XVuo
zA6fcUV%<A9Gp24@CHC?FYuTNj$_44&t^_ZiD3~4{sgyb_Xqf}T5ge8YdO~+@9sC=g
zu!ION%_TX82Z<gvqa+{24pJv6-lbH(NlBFGRo72S?>2loOVyC(*OTIw>sP{{7VpQ&
z@Xh^Ams2~*V`q^`iGd7{l3YaDmE>^pZMjwfCpGtt$ZNUQ$Y$Vg+Og~7N`~zFQBZBf
zE|Ls;km%GH6rRM0f-p{WU>No<Xv9kwlHu-&lniAjk~Ng=rd3E+Zw+_1%F3&~wQ?K^
ztB}UtGU!AK=MarVNA4y%?lr6$4I9~T?AT@_ijr<3>y{qZ>eeWl4RVAzR*AZBM)eGO
zU|0pLi;27)r@tZp4tG<}Dx_Iw6*(W+b!#Zg4N%!)@@f|oBg9ekYY{#VpZ{@VbuXw!
zt8c)|Re0{Wd#m+k#dV#CU5i$4;<32CcJ~haTumRCcijE>okoni^$oYOQg0@6C#2pC
zEJQy;BZ~8KN|v;zxh=~gf6sDHm*e;?%ctam_Nu%f=k>DG<)kdhS$v<7Ca5R07yIZ}
zF#DqXwvm^4y$I<Nq|!@!RvX92BEGN4SF{V@q+cyBzbo@v_yWAyq*IMP^SI9AVhp>F
zlLN&Nl_@B}fL;McZ}g3z22ms6GIWkLv7@UR3Kc0)jAy$lAb{HS5ov_Fg$_U?s2r#v
zrt+YG2!S?X43YYbAS%igOkptkh!8U>T~wEcO`!~9Q%wNNoY=_|7zILw96^_2HVEPF
z;7AHgQ49f?E@N6r0aHu@Foi>m1IbTXeZNI;`hjUld{|0>X_2rjPy|{E&7LHG(S^1y
z-0ar`sDFfge8}cCL_z7LaFOmmvu^m9>84urm|b{mhOz)<-J%)+r`H2HMwJXEfrGf>
z#uNgokMWD5`6sE~62KS2I8et40;Cg<aJ=yv0mv1vGq53`85hs@`pjE&M!6RKNTWaj
zU#j5HTq%J)*GZVM(++^>1gRs~pV5FAcr8G4fTtL#qOipr0@~Pz38a$(bwnwz0Jbq5
zV3rw~Kj;Ic^XmYD<4PfxG)dl(LO>e~0@?rpl@$a5{KOD&OK%9+t}9v|APRJIRCJb8
zdy2Oc5uf7HzdsrjplN0UH~>2Wn;&J*FdL|3=2$i$4>mH#5blVydPk#fg_i-1d$$6B
z5QM(|t>DZHU`DJXMJm1lr&DwS8n9O(640|Lf(c^-TTW7QKqL$HBNT<lAQC*VsvfY2
zzd=vey^We(|Gz^LdIVioMR`Wr`$VmEK?wKypaGs2w2a&9w1~IudIW56XL~R_$aLc|
zo{_1Up6#_mVBj0Ga%}m}Wqe4bX6y9<mU}p0P$?f$-5L@m2)X87J(b2@VQ-C4A+oWJ
z`rQp%mj+57SLwHI3=}zCq2GG<zflrzXDGXi7&wm<6xez=P4)Yy*IM5jsQmk=8T+7!
z=bOE7M^#S&Z$I)79@n`Q((7vl8pPkPrG&WNx;oGXYCz7!ZUNrkZ@oHLHcR8K@S%az
zSzhZ3XEz)d@~5-h6HYCwf!4}Ezu)tB{osM$1BG{%X2MCz(@|>1s2Qhbf|@C6rm3Og
zFPx=@isO(<tt3~e1}N|Ks8I_|44~$2RAaARZCZIdiko%U%01Y|23dK>+oP>oB-yjn
zoI|rTA3i}e{ojG>b>0nCrYn+)c%s&NgKNrxr%3Q=YVgk*5k5oBv(#KfW99u8mEu;;
zZ=nXajG7(oT88I4`&K@*Nu!mc9`sVZt&FLPThZO&wX_gtFT&@@0#(|oSq7CX4v6;y
zOo(>Sh+!_4QR*HpmS^NbS}d2TSk}rYo}WayYf4kowGs;9moQ>zPs<sS)Gz8Y=ufFO
zEr_3j6v}NqFE3F^d<oW;%1Lnp|K-6OE_Y)ir1uI0v8q*<|Lqf5LjW=^fjDa@c58Kl
z-&Ua-Y;Ae|mX(jT8?jfl3KgGU1SDGMa%Lo25CNe?D|I{GP7G=-EvO^Rjat2_vV|&g
m_MtCM95D{o(5tn;X;j_s6fiA1gJyzPNh|Gr0Rptr{r>>lg!4oI
diff --git a/EddyS2VTest/feedsInputs~ b/EddyS2VTest/feedsInputs~
deleted file mode 100755
index dcfa4da..0000000
--- a/EddyS2VTest/feedsInputs~
+++ /dev/null
@@ -1 +0,0 @@
-EddyTestData/eddyData
diff --git a/EddyS2VTest/feedsRun b/EddyS2VTest/feedsRun
index 2d886ca..f6d8a39 100755
--- a/EddyS2VTest/feedsRun
+++ b/EddyS2VTest/feedsRun
@@ -12,27 +12,30 @@
# if/when there is an openmp implementation of slice-to-vol.
#
-scriptdir=`dirname $0`
odir=$1
-indir=$2
+strict=$2
+indir=$3
#
-# Inputs 1--2 are the ones neccessary for feeds to work
+# Inputs 1--3 are the ones neccessary for feeds to work
# Additional inputs are optional and intended for testing
# outside of the feeds context.
#
-# Input 3 is alternative location of executable
+# Input 4 is alternative location of executable
#
-if [ "$#" -gt 2 ]; then
- exedir=$3
+if [ "$#" -gt 3 ]; then
+ exedir=$4
echo "Directory for eddy executables set to ${exedir}"
else
exedir="${FSLDIR}/bin"
fi
-if [ ! -x "${exedir}/eddy_cuda" ]; then
- echo "Executable ${exedir}/eddy_cuda does not exist"
- exit 1
-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
@@ -47,138 +50,84 @@ if [ ! -d "$indir" ]; then
exit 1;
fi
-# Define some constants
-
-# These are for the comparison of output images divided on b=0, b=700 and b=2000
-max_mean_ima_diff="1.8 0.6 0.5"
-max_max_ima_diff="2.0 0.7 0.6"
-
# Launch both GPU and CPU versions
-cuda_jid=`fsl_sub -q cuda.q ${exedir}/eddy_cuda --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
+cuda_jid=""
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+ tmp=`basename $cuda_exe`
+ version=`echo $tmp | sed 's/eddy_cuda//'`
+ cuda_jid=`fsl_sub -q cuda.q ${cuda_exe} --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --ol_nstd=6 --ol_nvox=500 --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyCudaOutput${version} --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
+done
-# omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
+# omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --ol_nstd=6 --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyOmpOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
+# Ensure that slots are being reserved on the queue
+# qalter ${omp_jid} -R y
# Wait until they have both finished
while [ : ]; do
- tmp=`qstat -j ${cuda_jid} | wc`
- tmp=`echo $tmp | awk '{print $1}'`
-# 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
+ 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
+ break
+ fi
+# 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
-if [ ! -f ${odir}/eddyCudaOutput.nii* ]; then
- exit 1
-fi
+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 ground truth
-
-# First compare output images
+# Define some constants
+# These are for the comparison of output images divided on b=0, b=700 and b=2000
+max_ima_diff="2.0 2.5 0.6 0.7 0.5 0.6"
+# Translations
+max_trans_diff="0.2 1.0"
+# Rotations
+max_rot_diff="0.2 2.0"
+# Outliers
+max_false_pos=10
+max_false_neg=10
-mmidiff=""
-maxidiff=""
-# for oname in eddyCudaOutput eddyOmpOutput; do
-for oname in eddyCudaOutput; do
- fslmaths ${indir}/EddyS2VTestData/eddyData/GroundTruth/ground_truth -sub ${odir}/${oname} -abs ${odir}/${oname}_diff_tmp
- # Divide into different shells
- for bv in 0 700 2000; do
- select_dwi_vols ${odir}/${oname}_diff_tmp ${indir}/EddyS2VTestData/eddyData/testBvals ${odir}/${oname}_diff_tmp_b${bv} ${bv}
- done
- imrm ${odir}/${oname}_diff_tmp
- # Remove top and bottom slides with no data in the eddy output
- for bv in 0 700 2000; do
- fslroi ${odir}/${oname}_diff_tmp_b${bv} ${odir}/${oname}_diff_tmp_b${bv} 2 68 4 80 6 44
- done
- fslroi ${indir}/EddyS2VTestData/eddyData/testMask ${odir}/tmp_mask 2 68 4 80 6 44
+# Check the results against ground truth
- for bv in 0 700 2000; do
- mdiff=$(fslstats -t ${odir}/${oname}_diff_tmp_b${bv} -k ${odir}/tmp_mask -m)
- imrm ${odir}/${oname}_diff_tmp_b${bv}
- echo $mdiff > ${odir}/${oname}.ima_mean_diff_b${bv}.txt
- echo "feedsRun: Mean difference per volume is: ${mdiff} for shell b=${bv}"
- mmidiff="${mmidiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxidiff="${maxidiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
- done
+cuda_exit_status=0
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+ tmp=`basename $cuda_exe`
+ version=`echo $tmp | sed 's/eddy_cuda//'`
+ ./S2VFeeds.py ${odir} Cuda${version} `imglob -extension ${indir}/EddyTestData/eddyData/testMask` ${indir}/EddyS2VTestData/eddyData/testBvals `imglob -extension ${odir}/eddyCudaOutput${version}` `imglob -extension ${indir}/EddyS2VTestData/eddyData/GroundTruth/ground_truth` ${odir}/eddyCudaOutput${version}.eddy_movement_over_time ${indir}/EddyS2VTestData/eddyData/MovementTruth/ground_truth.eddy_movement_over_time ${odir}/eddyCudaOutput${version}.eddy_outlier_map ${indir}/EddyS2VTestData/eddyData/OutlierTruth/dropoutLog_1.txt ${max_ima_diff} ${max_false_pos} ${max_false_neg} ${max_trans_diff} ${max_rot_diff}
+ cuda_exit_status=$(($cuda_exit_status + $?))
done
-echo "Avrage mean abolute differences for shells 0, 700 and 2000 are: ${mmidiff}"
-echo "Acceptable average mean absolute differences for shells 0, 700 and 2000 are: ${max_mean_ima_diff}"
-echo "Maximum mean absolute differences for shells 0, 700 and 2000 are: ${maxidiff}"
-echo "Acceptable maximum mean absolute differences for shells 0, 700 and 2000 are: ${max_max_ima_diff}"
-# Next decide on fail or pass
+omp_exit_status=0
-# Mean diff for b0
-fail=$(echo "$(echo $mmidiff | awk '{print $1}') > $(echo $max_mean_ima_diff | awk '{print $1}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=0 too large"
- exit 1;
-fi
-# Mean diff for b700
-fail=$(echo "$(echo $mmidiff | awk '{print $2}') > $(echo $max_mean_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=700 too large"
- exit 1;
-fi
-# Mean diff for b2000
-fail=$(echo "$(echo $mmidiff | awk '{print $3}') > $(echo $max_mean_ima_diff | awk '{print $3}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=2000 too large"
- exit 1;
-fi
-
-# Max diff for b0
-fail=$(echo "$(echo $maxidiff | awk '{print $1}') > $(echo $max_max_ima_diff | awk '{print $1}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=0 too large"
- exit 1;
-fi
-# Max diff for b700
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > $(echo $max_max_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=700 too large"
- exit 1;
-fi
-# Max diff for b2000
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > $(echo $max_max_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=2000 too large"
- exit 1;
-fi
-
-#
-# Next do checks on output text files
-#
-# First .eddy_movement_over_time
-#
-echo "Now running CompareMovementOverTime.py"
-${scriptdir}/CompareMovementOverTime.py ${odir}/${oname}.eddy_movement_over_time ${indir}/EddyS2VTestData/eddyData/MovementTruth/ground_truth.eddy_movement_over_time 1.0 1.0 0.2 3.0 1.5 0.3 1 55 7 46
-fail=$?
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed in CompareMovementOverTime.py"
- exit 1;
-fi
-
-#
-# Next .eddy_outlier_n_stdev_map
-#
-echo "Now running CompareOutlierNStdevMaps.py"
-${scriptdir}/CompareOutlierNStdevMaps.py ${odir}/${oname}.eddy_outlier_n_stdev_map ${indir}/EddyS2VTestData/eddyData/NStdevMapRef/ref.eddy_outlier_n_stdev_map 0.2 0.03 0.05 0.02 55 7 46
-fail=$?
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed in CompareOutlierNStdevMaps.py"
- exit 1;
+if [ $cuda_exit_status -gt 0 ] || [ $omp_exit_status -gt 0 ]; then
+ exit 1
+else
+ exit 0
fi
-
-exit 0
diff --git a/EddyS2VTest/feedsRun~ b/EddyS2VTest/feedsRun~
deleted file mode 100755
index 482be6c..0000000
--- a/EddyS2VTest/feedsRun~
+++ /dev/null
@@ -1,190 +0,0 @@
-#! /bin/sh
-#
-# This script runs eddy on data with intra-volume movement,
-# followed by some tests on the output.
-# It uses simulated data supplied by Mark
-# Graham, UCL, which allows us to have a ground
-# truth. The data is a subset of that used
-# for the s2v paper.
-# At the time of writing this (May, 2017) the openmp version
-# is not able to do slice-to-vol realignment. Hence only
-# the cuda version is tested. The script should be updated
-# if/when there is an openmp implementation of slice-to-vol.
-#
-
-scriptdir=`dirname $0`
-odir=$1
-indir=$2
-#
-# Inputs 1--2 are the ones neccessary for feeds to work
-# Additional inputs are optional and intended for testing
-# outside of the feeds context.
-#
-# Input 3 is alternative location of executable
-#
-if [ "$#" -gt 2 ]; then
- exedir=$3
- echo "Directory for eddy executables set to ${exedir}"
-else
- exedir="${FSLDIR}/bin"
-fi
-
-if [ ! -x "${exedir}/eddy_cuda" ]; then
- echo "Executable ${exedir}/eddy_cuda does not exist"
- exit 1
-fi
-#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 ${odir} does not exist"
- exit 1;
-fi
-
-# Define some constants
-
-# These are for the comparison of output images divided on b=0, b=700 and b=2000
-max_mean_ima_diff="1.8 0.6 0.5"
-max_max_ima_diff="2.0 0.7 0.6"
-
-# Launch both GPU and CPU versions
-
-if false; then
-
-cuda_jid=`fsl_sub -q cuda.q ${exedir}/eddy_cuda --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
-
-# omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/eddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose`
-
-# Wait until they have both finished
-
-while [ : ]; do
- tmp=`qstat -j ${cuda_jid} | wc`
- tmp=`echo $tmp | awk '{print $1}'`
-# 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
-
-if [ ! -f ${odir}/eddyCudaOutput.nii* ]; then
- exit 1
-fi
-#if [ ! -f ${odir}/eddyOmpOutput.nii* ]; then
-# exit 1
-#fi
-
-# Check the results against ground truth
-
-# First compare output images
-
-mmidiff=""
-maxidiff=""
-# for oname in eddyCudaOutput eddyOmpOutput; do
-for oname in eddyCudaOutput; do
- fslmaths ${indir}/EddyS2VTestData/eddyData/GroundTruth/ground_truth -sub ${odir}/${oname} -abs ${odir}/${oname}_diff_tmp
- # Divide into different shells
- for bv in 0 700 2000; do
- select_dwi_vols ${odir}/${oname}_diff_tmp ${indir}/EddyS2VTestData/eddyData/testBvals ${odir}/${oname}_diff_tmp_b${bv} ${bv}
- done
- imrm ${odir}/${oname}_diff_tmp
- # Remove top and bottom slides with no data in the eddy output
- for bv in 0 700 2000; do
- fslroi ${odir}/${oname}_diff_tmp_b${bv} ${odir}/${oname}_diff_tmp_b${bv} 2 68 4 80 6 44
- done
- fslroi ${indir}/EddyS2VTestData/eddyData/testMask ${odir}/tmp_mask 2 68 4 80 6 44
-
- for bv in 0 700 2000; do
- mdiff=$(fslstats -t ${odir}/${oname}_diff_tmp_b${bv} -k ${odir}/tmp_mask -m)
- imrm ${odir}/${oname}_diff_tmp_b${bv}
- echo $mdiff > ${odir}/${oname}.ima_mean_diff_b${bv}.txt
- echo "feedsRun: Mean difference per volume is: ${mdiff} for shell b=${bv}"
- mmidiff="${mmidiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxidiff="${maxidiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
- done
-done
-echo "Avrage mean abolute differences for shells 0, 700 and 2000 are: ${mmidiff}"
-echo "Acceptable average mean absolute differences for shells 0, 700 and 2000 are: ${max_mean_ima_diff}"
-echo "Maximum mean absolute differences for shells 0, 700 and 2000 are: ${maxidiff}"
-echo "Acceptable maximum mean absolute differences for shells 0, 700 and 2000 are: ${max_max_ima_diff}"
-
-# Next decide on fail or pass
-
-# Mean diff for b0
-fail=$(echo "$(echo $mmidiff | awk '{print $1}') > $(echo $max_mean_ima_diff | awk '{print $1}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=0 too large"
- exit 1;
-fi
-# Mean diff for b700
-fail=$(echo "$(echo $mmidiff | awk '{print $2}') > $(echo $max_mean_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=700 too large"
- exit 1;
-fi
-# Mean diff for b2000
-fail=$(echo "$(echo $mmidiff | awk '{print $3}') > $(echo $max_mean_ima_diff | awk '{print $3}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because mean diff for b=2000 too large"
- exit 1;
-fi
-
-# Max diff for b0
-fail=$(echo "$(echo $maxidiff | awk '{print $1}') > $(echo $max_max_ima_diff | awk '{print $1}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=0 too large"
- exit 1;
-fi
-# Max diff for b700
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > $(echo $max_max_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=700 too large"
- exit 1;
-fi
-# Max diff for b2000
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > $(echo $max_max_ima_diff | awk '{print $2}')" | bc)
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed because max diff for b=2000 too large"
- exit 1;
-fi
-
-fi # End if false
-
-oname="eddyCudaOutput"
-
-#
-# Next do checks on output text files
-#
-# First .eddy_movement_over_time
-#
-echo "Now running CompareMovementOverTime.py"
-${scriptdir}/CompareMovementOverTime.py ${odir}/${oname}.eddy_movement_over_time ${indir}/EddyS2VTestData/eddyData/MovementTruth/ground_truth.eddy_movement_over_time 1.0 1.0 0.2 3.0 1.5 0.3 1 55 7 46
-fail=$?
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed in CompareMovementOverTime.py"
- exit 1;
-fi
-
-#
-# Next .eddy_outlier_n_stdev_map
-#
-echo "Now running CompareOutlierNStdevMaps.py"
-${scriptdir}/CompareOutlierNStdevMaps.py ${odir}/${oname}.eddy_outlier_n_stdev_map ${indir}/EddyS2VTestData/eddyData/NStdevMapRef/ref.eddy_outlier_n_stdev_map 0.2 0.03 0.05 0.02 55 7 46
-fail=$?
-if [ $fail -gt 0 ]; then
- echo "feedsRun: Test failed in CompareOutlierNStdevMaps.py"
- exit 1;
-fi
-
-exit 0
diff --git a/EddyS2VTest/nohup.sh b/EddyS2VTest/nohup.sh
new file mode 100755
index 0000000..87acbcc
--- /dev/null
+++ b/EddyS2VTest/nohup.sh
@@ -0,0 +1,2 @@
+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyS2VTest blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?
diff --git a/EddyS2VTest/rubbestad.txt b/EddyS2VTest/rubbestad.txt
new file mode 100644
index 0000000..3d7a911
--- /dev/null
+++ b/EddyS2VTest/rubbestad.txt
@@ -0,0 +1,12 @@
+export indir=/vols/Data/fsldev/dataSets
+export odir=/vols/Scratch/HCP/Diffusion/jesper/EddyS2VFeedsTest20190129
+exedir=/home/fs0/jesper/DebugEddyNewNewimage/eddy_new/eddy
+
+${exedir}/eddy_cuda --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/NewEddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose
+
+
+export indir=/vols/Data/fsldev/dataSets
+export odir=/vols/Scratch/HCP/Diffusion/jesper/EddyS2VFeedsTest20190129
+exedir=/home/fs0/jesper/fsl/src/BuildingEddyForHCP_2017_09_29/eddy
+
+${exedir}/eddy_cuda8.0 --imain=${indir}/EddyS2VTestData/eddyData/testData --mask=${indir}/EddyS2VTestData/eddyData/testMask --bvals=${indir}/EddyS2VTestData/eddyData/testBvals --bvecs=${indir}/EddyS2VTestData/eddyData/testBvecs --index=${indir}/EddyS2VTestData/eddyData/testIndex --acqp=${indir}/EddyS2VTestData/eddyData/testAcqparams --repol --ol_type=slice --niter=8 --fwhm=10,6,4,2,0,0,0,0 --out=${odir}/OldEddyCudaOutput --nvoxhp=5000 --mporder=16 --s2v_niter=10 --s2v_interp=trilinear --s2v_lambda=1 --slspec=${indir}/EddyS2VTestData/eddyData/testSlspec --very_verbose
diff --git a/EddyTest/EddyFeeds.py b/EddyTest/EddyFeeds.py
new file mode 100755
index 0000000..f427bd1
--- /dev/null
+++ b/EddyTest/EddyFeeds.py
@@ -0,0 +1,278 @@
+#!/usr/bin/env fslpython
+
+import sys
+import numpy as np
+import nibabel as nib
+import datetime
+import pandas as pd
+
+class EddyFeedsType(object):
+ """ The purpose of this class is to make all the comparisons between
+ newly estimated and precomputed eddy results for plain vanilla and
+ outlier eddy.
+ """
+
+ def __init__(self,mask,corr,precomp_corr,disp_field_base,precomp_disp_field_base,ol_txt,precomp_ol_txt):
+ # 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('EddyFeedsType:__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('EddyFeeds:__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('EddyFeedsType:__init__:Size mismatch in fourth dimension of corrected images')
+
+ # Read displacement fields and make sure dimensions are right
+ try:
+ self._fields = []
+ self._precomp_fields = []
+ for vol in range(1,self._corr.header['dim'][4]+1):
+ field_name = disp_field_base + ".{:03d}".format(vol)
+ self._fields.append(self.OpenImageFile(field_name))
+ precomp_field_name = precomp_disp_field_base + str(vol)
+ self._precomp_fields.append(self.OpenImageFile(precomp_field_name))
+ except Exception as e:
+ print(str(e))
+ raise Exception('EddyFeedsType:__init__:Error opening displacement field image files')
+
+ for vol in range (0,self._corr.header['dim'][4]):
+ if not (all(self._mask.header['dim'][1:4] == self._fields[vol].header['dim'][1:4]) and
+ all(self._mask.header['dim'][1:4] == self._precomp_fields[vol].header['dim'][1:4])):
+ raise Exception('EddyFeedsType:__init__:Size mismatch in first three dimensions of displacement field image files')
+ if self._fields[vol].header['dim'][4] != 3 or self._precomp_fields[vol].header['dim'][4] != 3:
+ raise Exception('EddyFeedsType:__init__:Size mismatch in fourth dimension of displacement field image files')
+
+ # Read outlier text files
+ try:
+ self._ol = pd.read_csv(ol_txt, delim_whitespace=True, header=None, skiprows=1).values
+ self._precomp_ol = pd.read_csv(precomp_ol_txt, sep=',', header=None, decimal='.').values
+ if self._ol.shape[0] != self._corr.header['dim'][4] or \
+ self._precomp_ol.shape[0] != self._corr.header['dim'][4]:
+ raise Exception('EddyFeedsType:__init__:Size mismatch in volume dimension of outlier data')
+ if self._ol.shape[1] != self._corr.header['dim'][3] or \
+ self._precomp_ol.shape[1] != self._corr.header['dim'][3]:
+ raise Exception('EddyFeedsType:__init__:Size mismatch in slice dimension of outlier data')
+
+ except Exception as e:
+ print(str(e))
+ raise Exception('EddyFeedsType:__init__:Error opening outlier text files')
+
+ # 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()
+
+ self._fielddiffmeans = np.zeros(len(self._fields))
+ for vol in range(0,len(self._fields)):
+ tmpdiff = self._fields[vol].get_data().astype(float) - self._precomp_fields[vol].get_data().astype(float)
+ tmpdiff = np.square(tmpdiff[:,:,:,0]) + np.square(tmpdiff[:,:,:,1]) + np.square(tmpdiff[:,:,:,2])
+ tmpdiff = np.multiply(mask,np.sqrt(tmpdiff))
+ self._fielddiffmeans[vol] = np.array(mask.shape).prod() * tmpdiff.mean() / mask.sum()
+
+ except Exception as e:
+ print(str(e))
+ raise Exception('EddyFeedsType:__init__:Error calculating image statistics')
+
+ # Calculate false positives, false negatives and true positives for outliers
+ try:
+ # Find slices with > 500 brain voxels
+ vsl = []
+ for sl in range(0,self._mask.header['dim'][3]):
+ if mask[:,:,sl].sum() > 500:
+ vsl.append(sl)
+ # Go through and count positives and false negatives
+ self._tol = 0; self._tdol = 0; self._tp = 0; self._fp = 0; self._fn = 0
+ for sl in vsl:
+ for vol in range(0,self._precomp_ol.shape[0]):
+ if not np.isnan(self._precomp_ol[vol,sl]): # If it is a true outlier
+ self._tol += 1
+ if self._ol[vol,sl] == 1:
+ self._tp += 1
+ if self._precomp_ol[vol,sl] < 0.8: # If it ought to have been detected
+ self._tdol += 1
+ if self._ol[vol,sl] != 1:
+ self._fn += 1
+ else:
+ if self._ol[vol,sl] == 1:
+ self._fp += 1
+
+ except Exception as e:
+ print(str(e))
+ raise Exception('EddyFeedsType:__init__:Error calculating outlier statistics')
+
+ def OpenImageFile(self,fname):
+ # Try to open 'fname', 'fname'.nii and 'fname'.nii.gz in that order
+ try:
+ rval = nib.load(fname,mmap=False)
+ except:
+ try:
+ rval = nib.load(fname+'.nii',mmap=False)
+ except:
+ try:
+ rval = nib.load(fname+'.nii.gz',mmap=False)
+ except:
+ raise Exception('EddyFeedsType:OpenImageFile:Error opening file ' + fname)
+ return rval
+
+ def MeanDiffsOfCorrectedImages(self):
+ return self._corrdiffmeans
+
+ def MeanDiffsOfFields(self):
+ return self._fielddiffmeans
+
+ def OlFalsePositives(self):
+ return self._fp
+
+ def OlFalseNegatives(self):
+ return self._fn
+
+ def OlTruePositives(self):
+ return self._tp
+
+ def OlTotalNumber(self):
+ return self._tol
+
+ def OlDetectableNumber(self):
+ return self._tdol
+
+def main(argv):
+ # This is the main program that tests the output from running
+ # the eddy (vanilla and outliers) tests.
+ try:
+ if len(argv) != 16:
+ print('EddyFeeds.py usage: EddyFeeds output_dir prefix mask corrected precomputed_corrected fields_basename ' \
+ 'precomputed_fields_basename ol precomputed_ol allowed_mean_diff_corrected ' \
+ 'allowed_max_diff_corrected allowed_mean_diff_field allowed_max_diff_field ' \
+ 'allowed_false_positives allowed_false_negatives')
+ sys.exit(1)
+ else:
+ output_dir = argv[1]
+ output_prefix = argv[2]
+ mask = argv[3]
+ corrected = argv[4]
+ precomputed_corrected = argv[5]
+ fields_basename = argv[6]
+ precomputed_fields_basename = argv[7]
+ ol = argv[8]
+ precomputed_ol = argv[9]
+ allowed_mean_diff_corrected = float(argv[10])
+ allowed_max_diff_corrected = float(argv[11])
+ allowed_mean_diff_field = float(argv[12])
+ allowed_max_diff_field = float(argv[13])
+ allowed_false_positives = float(argv[14])
+ allowed_false_negatives = float(argv[15])
+
+ # Try to create EddyFeedsType object (involves reading all files)
+ try:
+ ef = EddyFeedsType(mask,corrected,precomputed_corrected,fields_basename, \
+ precomputed_fields_basename,ol,precomputed_ol)
+ except Exception as e:
+ print(str(e))
+ print('main: Error when creating EddyFeedsType 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 or \
+ ef.MeanDiffsOfFields().mean() > allowed_mean_diff_field or \
+ ef.MeanDiffsOfFields().max() > allowed_max_diff_field:
+ passes_test = False
+
+ # Check pass/fail based on detected outliers
+ if ef.OlFalsePositives() > allowed_false_positives or \
+ ef.OlFalseNegatives() > allowed_false_negatives:
+ 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 + '_EddyReport.txt','w')
+ except Exception as e:
+ print(str(e))
+ print('main: Cannot open report file: ' + output_dir + '/' + output_prefix + '_EddyReport.txt')
+ sys.exit(1)
+ else:
+ try:
+ fp.write('EddyFeeds 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')
+ # Report on differences in estimated displacement fields
+ fp.write('\nThe absolute differences, averaged across the mask, for the estimated displacement ' + \
+ 'images were: ' + ' '.join(["{0:.4f}".format(elem) for elem in ef.MeanDiffsOfFields()]) + '\n')
+ fp.write('That gives mean error: ' + "{0:.4f}".format(ef.MeanDiffsOfFields().mean()) + \
+ ', and a max error: ' + "{0:.4f}".format(ef.MeanDiffsOfFields().max()) + '\n')
+ fp.write('The allowed mean of averaged differences for the estimated displacement images was: ' + \
+ "{0:.4f}".format(allowed_mean_diff_field) + '\n')
+ fp.write('The allowed maximum of averaged differences for the estimated displacement images was: ' + \
+ "{0:.4f}".format(allowed_max_diff_field) + '\n')
+ if ef.MeanDiffsOfFields().mean() > allowed_mean_diff_field or \
+ ef.MeanDiffsOfFields().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')
+ # Report on differences in outlier estimations
+ fp.write('\nThe input data had ' + '{:d}'.format(ef.OlTotalNumber()) + ' true outliers\n')
+ fp.write('\nOf these ' + '{:d}'.format(ef.OlDetectableNumber()) + ' had a signal loss greater ' + \
+ 'than 20% and should be detected\n')
+ fp.write('\nThe analysis yielded ' + '{:d}'.format(ef.OlTruePositives()) + ' true positives, ' + \
+ '{:d}'.format(ef.OlFalseNegatives()) + ' false negatives and ' + \
+ '{:d}'.format(ef.OlFalsePositives()) + ' false positives\n')
+ if ef.OlFalsePositives() > allowed_false_positives or \
+ ef.OlFalseNegatives() > allowed_false_negatives:
+ 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] + '_EddyReport.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
diff --git a/EddyTest/eddy_cuda.e2284102 b/EddyTest/eddy_cuda.e2284102
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_cuda.e2286236 b/EddyTest/eddy_cuda.e2286236
deleted file mode 100755
index 2ecbbab..0000000
--- a/EddyTest/eddy_cuda.e2286236
+++ /dev/null
@@ -1,48 +0,0 @@
-
-Part of FSL (build 508)
-eddy
-Copyright(c) 2015, University of Oxford (Jesper Andersson)
-
-Usage:
-eddy --monsoon
-
-Compulsory arguments (You MUST set one or more of):
- --imain File containing all the images to estimate distortions for
- --mask Mask to indicate brain
- --index File containing indices for all volumes in --imain into --acqp and --topup
- --acqp File containing acquisition parameters
- --bvecs File containing the b-vectors for all volumes in --imain
- --bvals File containing the b-values for all volumes in --imain
- --out Basename for output
-
-Optional arguments (You may optionally specify one or more of):
- --mb Multi-band factor
- --mb_offs Multi-band offset (-1 if bottom slice removed, 1 if top slice removed)
- --topup Base name for output files from topup
- --field Name of file with susceptibility field (in Hz)
- --field_mat Name of rigid body transform for susceptibility field
- --flm First level EC model (linear/quadratic/cubic, default quadratic)
- --slm Second level EC model (none/linear/quadratic, default none)
- --fwhm FWHM for conditioning filter when estimating the parameters (default 0)
- --niter Number of iterations (default 5)
- --fep Fill empty planes in x- or y-directions
- --interp Interpolation model for estimation step (spline/trilinear, default spline)
- --resamp Final resampling method (jac/lsr, default jac)
- --nvoxhp # of voxels used to estimate the hyperparameters (default 1000)
- --initrand Resets rand for when selecting voxels (default false)
- --ff Fudge factor for hyperparameter error variance (default 10.0)
- --repol Detect and replace outlier slices
- --ol_nstd Number of std off to qualify as outlier (default 4)
- --ol_nvox Min # of voxels in a slice for inclusion in outlier detection (default 250)
- --ol_type Type of outliers, slicewise (sw), groupwise (gw) or both (both). (default sw)
- --ol_pos Consider both positive and negative outliers if set (default false)
- --ol_sqr Consider outliers among sums-of-squared differences if set (default false)
- --dont_sep_offs_move Do NOT attempt to separate field offset from subject movement (default false)
- --dont_peas Do NOT perform a post-eddy alignment of shells (default false)
- --data_is_shelled Assume, don't check, that data is shelled (default false)
- -v,--verbose switch on diagnostic messages
- -h,--help display this message
-
-
-
---with_outliers: Option doesn't exist!
diff --git a/EddyTest/eddy_cuda.e2286247 b/EddyTest/eddy_cuda.e2286247
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_cuda.e2294981 b/EddyTest/eddy_cuda.e2294981
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_cuda.e2410534 b/EddyTest/eddy_cuda.e2410534
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_cuda.o2284102 b/EddyTest/eddy_cuda.o2284102
deleted file mode 100755
index 29c8021..0000000
--- a/EddyTest/eddy_cuda.o2284102
+++ /dev/null
@@ -1,88 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-
-...................Allocated GPU # 0...................
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 32.3879
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 69.4732
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6548
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3811
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2441
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Calculating parameter updates
-Iter: 0, Total mss = 0.794029
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 16.7945
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 16.3233
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 16.3347
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 16.1455
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_cuda.o2286236 b/EddyTest/eddy_cuda.o2286236
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_cuda.o2286247 b/EddyTest/eddy_cuda.o2286247
deleted file mode 100755
index 8aeadb5..0000000
--- a/EddyTest/eddy_cuda.o2286247
+++ /dev/null
@@ -1,90 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-
-...................Allocated GPU # 0...................
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 32.3879
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 69.473
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6673
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3908
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2527
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Calculating parameter updates
-Iter: 0, Total mss = 0.776861
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 16.6888
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 16.3153
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 16.2089
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 16.2899
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRegisteredImages
-Running sm.RecycleOutliers
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_cuda.o2294981 b/EddyTest/eddy_cuda.o2294981
deleted file mode 100755
index 2c52400..0000000
--- a/EddyTest/eddy_cuda.o2294981
+++ /dev/null
@@ -1,90 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-
-...................Allocated GPU # 0...................
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 32.3879
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 69.4732
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6548
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3811
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2441
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Calculating parameter updates
-Iter: 0, Total mss = 0.774561
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 16.8246
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 16.3174
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 16.1742
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 16.2312
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRegisteredImages
-Running sm.RecycleOutliers
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_cuda.o2410534 b/EddyTest/eddy_cuda.o2410534
deleted file mode 100755
index e35839a..0000000
--- a/EddyTest/eddy_cuda.o2410534
+++ /dev/null
@@ -1,88 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-
-...................Allocated GPU # 0...................
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 32.3879
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 69.473
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6665
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3911
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.253
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Calculating parameter updates
-Iter: 0, Total mss = 0.771488
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 16.7688
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 16.2489
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 16.1949
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 16.2025
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Entering EddyGpuUtils::LoadPredictionMaker
-Loading prediction maker
-Evaluating prediction maker model
-Checking for outliers
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_openmp.e2284103 b/EddyTest/eddy_openmp.e2284103
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.e2286237 b/EddyTest/eddy_openmp.e2286237
deleted file mode 100755
index 2ecbbab..0000000
--- a/EddyTest/eddy_openmp.e2286237
+++ /dev/null
@@ -1,48 +0,0 @@
-
-Part of FSL (build 508)
-eddy
-Copyright(c) 2015, University of Oxford (Jesper Andersson)
-
-Usage:
-eddy --monsoon
-
-Compulsory arguments (You MUST set one or more of):
- --imain File containing all the images to estimate distortions for
- --mask Mask to indicate brain
- --index File containing indices for all volumes in --imain into --acqp and --topup
- --acqp File containing acquisition parameters
- --bvecs File containing the b-vectors for all volumes in --imain
- --bvals File containing the b-values for all volumes in --imain
- --out Basename for output
-
-Optional arguments (You may optionally specify one or more of):
- --mb Multi-band factor
- --mb_offs Multi-band offset (-1 if bottom slice removed, 1 if top slice removed)
- --topup Base name for output files from topup
- --field Name of file with susceptibility field (in Hz)
- --field_mat Name of rigid body transform for susceptibility field
- --flm First level EC model (linear/quadratic/cubic, default quadratic)
- --slm Second level EC model (none/linear/quadratic, default none)
- --fwhm FWHM for conditioning filter when estimating the parameters (default 0)
- --niter Number of iterations (default 5)
- --fep Fill empty planes in x- or y-directions
- --interp Interpolation model for estimation step (spline/trilinear, default spline)
- --resamp Final resampling method (jac/lsr, default jac)
- --nvoxhp # of voxels used to estimate the hyperparameters (default 1000)
- --initrand Resets rand for when selecting voxels (default false)
- --ff Fudge factor for hyperparameter error variance (default 10.0)
- --repol Detect and replace outlier slices
- --ol_nstd Number of std off to qualify as outlier (default 4)
- --ol_nvox Min # of voxels in a slice for inclusion in outlier detection (default 250)
- --ol_type Type of outliers, slicewise (sw), groupwise (gw) or both (both). (default sw)
- --ol_pos Consider both positive and negative outliers if set (default false)
- --ol_sqr Consider outliers among sums-of-squared differences if set (default false)
- --dont_sep_offs_move Do NOT attempt to separate field offset from subject movement (default false)
- --dont_peas Do NOT perform a post-eddy alignment of shells (default false)
- --data_is_shelled Assume, don't check, that data is shelled (default false)
- -v,--verbose switch on diagnostic messages
- -h,--help display this message
-
-
-
---with_outliers: Option doesn't exist!
diff --git a/EddyTest/eddy_openmp.e2286248 b/EddyTest/eddy_openmp.e2286248
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.e2294982 b/EddyTest/eddy_openmp.e2294982
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.e2410535 b/EddyTest/eddy_openmp.e2410535
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.o2284103 b/EddyTest/eddy_openmp.o2284103
deleted file mode 100755
index bcbdef5..0000000
--- a/EddyTest/eddy_openmp.o2284103
+++ /dev/null
@@ -1,125 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 30.2634
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 71.4462
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6407
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3682
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2258
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.393602
-0.404945
-7.624107
-3.336045
-0.965445
-
-Calculating parameter updates
-Iter: 0, Total mss = 0.799956
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.113134
-0.460072
-7.997223
-3.640580
-2.813176
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.351879
-0.460348
-7.682541
-3.385224
-2.898969
-
-Calculating parameter updates
-Iter: 1, Total mss = 16.7099
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.950022
-0.449921
-8.132596
-5.077930
-4.520842
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.630933
-0.472084
-7.320908
-1.482529
-0.966824
-
-Calculating parameter updates
-Iter: 2, Total mss = 16.2991
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.197911
-0.454846
-7.634029
-3.300519
-2.780603
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.181265
-0.452707
-7.837860
-2.027355
-1.452874
-
-Calculating parameter updates
-Iter: 3, Total mss = 16.3176
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.326802
-0.458960
-7.505953
-4.476644
-3.911405
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.321295
-0.447338
-7.843067
-1.873306
-1.278127
-
-Calculating parameter updates
-Iter: 4, Total mss = 16.2606
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.826346
-0.460438
-7.820621
-1.563020
-1.110227
-
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_openmp.o2286237 b/EddyTest/eddy_openmp.o2286237
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.o2286248 b/EddyTest/eddy_openmp.o2286248
deleted file mode 100755
index cdd1763..0000000
--- a/EddyTest/eddy_openmp.o2286248
+++ /dev/null
@@ -1,127 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 30.2634
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 71.3976
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.654
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3726
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2296
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.781577
-0.397896
-7.600491
-4.419513
-2.102231
-
-Calculating parameter updates
-Iter: 0, Total mss = 0.791013
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.079796
-0.457485
-7.445464
-3.597443
-2.688557
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.993701
-0.475485
-7.568944
-1.113893
-0.638403
-
-Calculating parameter updates
-Iter: 1, Total mss = 16.8324
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.019007
-0.458953
-8.417402
-2.078327
-1.577257
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.282552
-0.468170
-7.376894
-3.346710
-2.906947
-
-Calculating parameter updates
-Iter: 2, Total mss = 16.3073
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.337659
-0.445033
-7.231357
-4.483216
-3.916252
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.506793
-0.465650
-7.748341
-3.707770
-3.215642
-
-Calculating parameter updates
-Iter: 3, Total mss = 16.4035
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.212213
-0.476378
-7.490193
-4.340426
-3.826962
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.258636
-0.467637
-8.267385
-4.409565
-3.911504
-
-Calculating parameter updates
-Iter: 4, Total mss = 16.313
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.753859
-0.430826
-7.275254
-1.513822
-1.081081
-
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRegisteredImages
-Running sm.RecycleOutliers
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_openmp.o2294982 b/EddyTest/eddy_openmp.o2294982
deleted file mode 100755
index c048395..0000000
--- a/EddyTest/eddy_openmp.o2294982
+++ /dev/null
@@ -1,127 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 30.2634
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 71.4505
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6378
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3696
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2274
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 3.695336
-0.419596
-8.511405
-7.360930
-5.019484
-
-Calculating parameter updates
-Iter: 0, Total mss = 0.793321
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.373124
-0.464419
-6.346941
-4.019980
-3.100530
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.631008
-0.449926
-8.586970
-3.740501
-3.259632
-
-Calculating parameter updates
-Iter: 1, Total mss = 16.8156
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.787431
-0.438354
-8.679474
-4.898026
-4.361649
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.247011
-0.445063
-7.402384
-4.405869
-3.862211
-
-Calculating parameter updates
-Iter: 2, Total mss = 16.3174
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 2.872214
-0.449623
-8.215931
-5.965500
-5.448123
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.660907
-0.457384
-7.893076
-3.807008
-3.319137
-
-Calculating parameter updates
-Iter: 3, Total mss = 16.3301
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.566298
-0.461746
-8.001916
-2.550802
-2.030200
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.293848
-0.455347
-7.411443
-2.850402
-2.316608
-
-Calculating parameter updates
-Iter: 4, Total mss = 16.2609
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.322931
-0.457442
-7.576957
-1.178943
-0.638636
-
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRegisteredImages
-Running sm.RecycleOutliers
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_openmp.o2410535 b/EddyTest/eddy_openmp.o2410535
deleted file mode 100755
index 665d35a..0000000
--- a/EddyTest/eddy_openmp.o2410535
+++ /dev/null
@@ -1,125 +0,0 @@
-Reading images
-Performing volume-to-volume registration
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 0, Total mss = 30.2634
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 1, Total mss = 71.4428
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 2, Total mss = 52.6283
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 3, Total mss = 52.3574
-Loading prediction maker
-Evaluating prediction maker model
-Calculating parameter updates
-Iter: 4, Total mss = 52.2153
-Running sm.Setb0Reference
-Running sm.PolateB0MovPar
-Running Register
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.108949
-0.405655
-7.539981
-3.458005
-1.114089
-
-Calculating parameter updates
-Iter: 0, Total mss = 0.773125
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.979832
-0.447969
-8.008780
-2.571105
-1.650663
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.553886
-0.468308
-7.954551
-4.676186
-4.190508
-
-Calculating parameter updates
-Iter: 1, Total mss = 16.8301
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.102304
-0.442443
-8.047528
-3.229302
-2.705240
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.546262
-0.460540
-8.821884
-4.674658
-4.185523
-
-Calculating parameter updates
-Iter: 2, Total mss = 16.3353
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.874422
-0.457946
-8.299983
-2.281360
-1.774702
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 1.326475
-0.452787
-7.433080
-4.470536
-3.925454
-
-Calculating parameter updates
-Iter: 3, Total mss = 16.2504
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: 0.770232
-0.479239
-7.027165
-3.833927
-3.346439
-
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -1.255352
-0.446087
-7.977991
-1.917732
-1.403634
-
-Calculating parameter updates
-Iter: 4, Total mss = 16.3182
-Running sm.SetDWIReference
-Aligning shells (running PostEddyAlignShells)
-Performing final outlier check
-Loading prediction maker
-Evaluating prediction maker model
-Estimated hyperparameters: -0.080947
-0.463076
-7.670441
-0.720256
-0.258951
-
-Running sm.WriteOutlierFreeData
-Running sm.WriteParameterFile
-Running sm.WriteRegisteredImages
-Running sm.WriteRotatedBVecs
-Running sm.WriteMovementRMS
-Running sm.WriteDisplacementFields
diff --git a/EddyTest/eddy_openmp.pe2284103 b/EddyTest/eddy_openmp.pe2284103
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.pe2286237 b/EddyTest/eddy_openmp.pe2286237
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.pe2286248 b/EddyTest/eddy_openmp.pe2286248
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.pe2294982 b/EddyTest/eddy_openmp.pe2294982
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.pe2410535 b/EddyTest/eddy_openmp.pe2410535
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.po2284103 b/EddyTest/eddy_openmp.po2284103
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.po2286237 b/EddyTest/eddy_openmp.po2286237
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.po2286248 b/EddyTest/eddy_openmp.po2286248
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.po2294982 b/EddyTest/eddy_openmp.po2294982
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/eddy_openmp.po2410535 b/EddyTest/eddy_openmp.po2410535
deleted file mode 100755
index e69de29..0000000
diff --git a/EddyTest/feedsRun b/EddyTest/feedsRun
index 410e4eb..d3a7306 100755
--- a/EddyTest/feedsRun
+++ b/EddyTest/feedsRun
@@ -8,26 +8,6 @@
# for the outlier paper.
#
-mean_dfield_difference() {
- local mdd_ii=$1
- local mdd_odir=$2
- local mdd_obase=$3
- local mdd_tfname=$4
- local mdd_mask=$5
- if [ ${#mdd_ii} -eq 1 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.00${mdd_ii}.nii.gz"
- elif [ ${#mdd_ii} -eq 2 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.0${mdd_ii}.nii.gz"
- elif [ ${#mdd_ii} -eq 3 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.${mdd_ii}.nii.gz"
- fi
- fslmaths ${mdd_tfname} -sub ${mdd_ofname} -mul ${mdd_mask} -sqr -Tmean -mul 3 -sqrt ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp
- local mdd_mdiff=`fslstats ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp -k ${mdd_mask} -M`
- # imrm ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp
- echo ${mdd_mdiff}
- return 0
-}
-
odir=$1
strict=$2
indir=$3
@@ -45,10 +25,13 @@ else
exedir="${FSLDIR}/bin"
fi
-if [ ! -x "${exedir}/eddy_cuda" ]; then
- echo "Executable ${exedir}/eddy_cuda does not exist"
- exit 1
-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
@@ -59,113 +42,87 @@ if [ ! -d "$odir" ]; then
exit 1;
fi
if [ ! -d "$indir" ]; then
- echo "Input directory ${odir} does not exist"
+ echo "Input directory ${indir} does not exist"
exit 1;
fi
# Define some constants
-max_mean_dfield_diff=0.543
-max_max_dfield_diff=0.877
-max_mean_ima_diff=1.49
-max_max_ima_diff=3.18
+max_mean_dfield_diff=0.5
+max_max_dfield_diff=1.0
+max_mean_ima_diff=1.5
+max_max_ima_diff=3.0
+allowed_false_positives=1
+allowed_false_negatives=1
# Launch both GPU and CPU versions
-cuda_jid=`fsl_sub -q cuda.q ${exedir}/eddy_cuda --imain=${indir}/EddyTestData/eddyData/testData --mask=${indir}/EddyTestData/eddyData/testMask --bvals=${indir}/EddyTestData/eddyData/testBvals --bvecs=${indir}/EddyTestData/eddyData/testBvecs --index=${indir}/EddyTestData/eddyData/testIndex --acqp=${indir}/EddyTestData/eddyData/testAcqparams --repol --fwhm=10,0,0,0,0 --out=${odir}/eddyCudaOutput --dfields -v`
+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}/EddyTestData/eddyData/testData --mask=${indir}/EddyTestData/eddyData/testMask --bvals=${indir}/EddyTestData/eddyData/testBvals --bvecs=${indir}/EddyTestData/eddyData/testBvecs --index=${indir}/EddyTestData/eddyData/testIndex --acqp=${indir}/EddyTestData/eddyData/testAcqparams --repol --fwhm=10,0,0,0,0 --out=${odir}/eddyCudaOutput${version} --dfields -v`"
+done
omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyTestData/eddyData/testData --mask=${indir}/EddyTestData/eddyData/testMask --bvals=${indir}/EddyTestData/eddyData/testBvals --bvecs=${indir}/EddyTestData/eddyData/testBvecs --index=${indir}/EddyTestData/eddyData/testIndex --acqp=${indir}/EddyTestData/eddyData/testAcqparams --repol --fwhm=10,0,0,0,0 --out=${odir}/eddyOmpOutput --dfields -v`
+# Ensure that slots are being reserved on the queue
+qalter ${omp_jid} -R y
# Wait until they have both finished
while [ : ]; do
- tmp=`qstat -j ${cuda_jid} | wc`
- tmp=`echo $tmp | awk '{print $1}'`
- 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
+ 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
-if [ ! -f ${odir}/eddyCudaOutput.nii* ]; then
- exit 1
-fi
+
+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 ground truth
-# First compare displacement fields
-
-mmddiff=""
-maxddiff=""
-for obase in eddyCudaOutput eddyOmpOutput; do
- mdiff=""
- for ((ii=1; ii<=108; ii+=1)); do
- mdiff="${mdiff} $(mean_dfield_difference ${ii} ${odir} ${obase} ${indir}/EddyTestData/eddyData/DisplacementFields/DisplacementField${ii} ${indir}/EddyTestData/eddyData/testMask)"
- done
- echo $mdiff > ${odir}/${obase}.dfield_mean_diff.txt
- mmddiff="${mmddiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxddiff="${maxddiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
+cuda_exit_status=0
+for cuda_exe in ${exedir}/eddy_cuda*;
+do
+ tmp=`basename $cuda_exe`
+ version=`echo $tmp | sed 's/eddy_cuda//'`
+ ./EddyFeeds.py ${odir} Cuda${version} `imglob -extension ${indir}/EddyTestData/eddyData/testMask` `imglob -extension ${odir}/eddyCudaOutput${version}` `imglob -extension ${indir}/EddyTestData/eddyData/GroundTruth/ground_truth` ${odir}/eddyCudaOutput${version}.eddy_displacement_fields ${indir}/EddyTestData/eddyData/DisplacementFields/DisplacementField ${odir}/eddyCudaOutput${version}.eddy_outlier_map ${indir}/EddyTestData/eddyData/OutlierTruth/dropoutLog_1.txt ${max_mean_ima_diff} ${max_max_ima_diff} ${max_mean_dfield_diff} ${max_max_dfield_diff} ${allowed_false_positives} ${allowed_false_negatives}
+ cuda_exit_status=$(($cuda_exit_status + $?))
done
-echo "mmddiff = ${mmddiff}"
-echo "maxddiff = ${maxddiff}"
-
-# Next compare output images
-
-mmidiff=""
-maxidiff=""
-for oname in eddyCudaOutput eddyOmpOutput; do
- fslmaths ${indir}/EddyTestData/eddyData/GroundTruth/ground_truth -sub ${odir}/${oname} -abs ${odir}/${oname}_diff_tmp
- mdiff=$(fslstats -t ${odir}/${oname}_diff_tmp -k ${indir}/EddyTestData/eddyData/testMask -m)
- imrm ${odir}/${oname}_diff_tmp
- echo $mdiff > ${odir}/${oname}.ima_mean_diff.txt
- mmidiff="${mmidiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxidiff="${maxidiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
-done
-echo "mmidiff = ${mmidiff}"
-echo "maxidiff = ${maxidiff}"
-# Next decide on fail or pass
+./EddyFeeds.py ${odir} Omp `imglob -extension ${indir}/EddyTestData/eddyData/testMask` `imglob -extension ${odir}/eddyOmpOutput` `imglob -extension ${indir}/EddyTestData/eddyData/GroundTruth/ground_truth` ${odir}/eddyOmpOutput.eddy_displacement_fields ${indir}/EddyTestData/eddyData/DisplacementFields/DisplacementField ${odir}/eddyOmpOutput.eddy_outlier_map ${indir}/EddyTestData/eddyData/OutlierTruth/dropoutLog_1.txt ${max_mean_ima_diff} ${max_max_ima_diff} ${max_mean_dfield_diff} ${max_max_dfield_diff} ${allowed_false_positives} ${allowed_false_negatives}
+omp_exit_status=$?
-fail=$(echo "$(echo $mmddiff | awk '{print $1}') > ${max_mean_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmddiff | awk '{print $2}') > ${max_mean_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxddiff | awk '{print $1}') > ${max_max_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxddiff | awk '{print $2}') > ${max_max_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmidiff | awk '{print $1}') > ${max_mean_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmidiff | awk '{print $2}') > ${max_mean_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxidiff | awk '{print $1}') > ${max_max_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > ${max_max_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
+if [ $cuda_exit_status -gt 0 ] || [ $omp_exit_status -gt 0 ]; then
+ exit 1
+else
+ exit 0
fi
-exit 0
+
diff --git a/EddyTest/feedsRun~ b/EddyTest/feedsRun~
deleted file mode 100755
index 410e4eb..0000000
--- a/EddyTest/feedsRun~
+++ /dev/null
@@ -1,171 +0,0 @@
-#! /bin/sh
-#
-# This script runs eddy followed by some
-# tests on the output.
-# It uses simulated data supplied by Mark
-# Graham, UCL, which allows us to have a ground
-# truth. The data is a subset of that used
-# for the outlier paper.
-#
-
-mean_dfield_difference() {
- local mdd_ii=$1
- local mdd_odir=$2
- local mdd_obase=$3
- local mdd_tfname=$4
- local mdd_mask=$5
- if [ ${#mdd_ii} -eq 1 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.00${mdd_ii}.nii.gz"
- elif [ ${#mdd_ii} -eq 2 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.0${mdd_ii}.nii.gz"
- elif [ ${#mdd_ii} -eq 3 ]; then
- local mdd_ofname="${mdd_odir}/${mdd_obase}.eddy_displacement_fields.${mdd_ii}.nii.gz"
- fi
- fslmaths ${mdd_tfname} -sub ${mdd_ofname} -mul ${mdd_mask} -sqr -Tmean -mul 3 -sqrt ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp
- local mdd_mdiff=`fslstats ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp -k ${mdd_mask} -M`
- # imrm ${mdd_odir}/${mdd_obase}_${mdd_ii}_tmp
- echo ${mdd_mdiff}
- return 0
-}
-
-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
-
-if [ ! -x "${exedir}/eddy_cuda" ]; then
- echo "Executable ${exedir}/eddy_cuda does not exist"
- exit 1
-fi
-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 ${odir} does not exist"
- exit 1;
-fi
-
-# Define some constants
-
-max_mean_dfield_diff=0.543
-max_max_dfield_diff=0.877
-max_mean_ima_diff=1.49
-max_max_ima_diff=3.18
-
-# Launch both GPU and CPU versions
-
-cuda_jid=`fsl_sub -q cuda.q ${exedir}/eddy_cuda --imain=${indir}/EddyTestData/eddyData/testData --mask=${indir}/EddyTestData/eddyData/testMask --bvals=${indir}/EddyTestData/eddyData/testBvals --bvecs=${indir}/EddyTestData/eddyData/testBvecs --index=${indir}/EddyTestData/eddyData/testIndex --acqp=${indir}/EddyTestData/eddyData/testAcqparams --repol --fwhm=10,0,0,0,0 --out=${odir}/eddyCudaOutput --dfields -v`
-
-omp_jid=`fsl_sub -q long.q -s openmp,6 ${exedir}/eddy_openmp --imain=${indir}/EddyTestData/eddyData/testData --mask=${indir}/EddyTestData/eddyData/testMask --bvals=${indir}/EddyTestData/eddyData/testBvals --bvecs=${indir}/EddyTestData/eddyData/testBvecs --index=${indir}/EddyTestData/eddyData/testIndex --acqp=${indir}/EddyTestData/eddyData/testAcqparams --repol --fwhm=10,0,0,0,0 --out=${odir}/eddyOmpOutput --dfields -v`
-
-# Wait until they have both finished
-
-while [ : ]; do
- tmp=`qstat -j ${cuda_jid} | wc`
- tmp=`echo $tmp | awk '{print $1}'`
- 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
-
-if [ ! -f ${odir}/eddyCudaOutput.nii* ]; then
- exit 1
-fi
-if [ ! -f ${odir}/eddyOmpOutput.nii* ]; then
- exit 1
-fi
-
-# Check the results against ground truth
-
-# First compare displacement fields
-
-mmddiff=""
-maxddiff=""
-for obase in eddyCudaOutput eddyOmpOutput; do
- mdiff=""
- for ((ii=1; ii<=108; ii+=1)); do
- mdiff="${mdiff} $(mean_dfield_difference ${ii} ${odir} ${obase} ${indir}/EddyTestData/eddyData/DisplacementFields/DisplacementField${ii} ${indir}/EddyTestData/eddyData/testMask)"
- done
- echo $mdiff > ${odir}/${obase}.dfield_mean_diff.txt
- mmddiff="${mmddiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxddiff="${maxddiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
-done
-echo "mmddiff = ${mmddiff}"
-echo "maxddiff = ${maxddiff}"
-
-# Next compare output images
-
-mmidiff=""
-maxidiff=""
-for oname in eddyCudaOutput eddyOmpOutput; do
- fslmaths ${indir}/EddyTestData/eddyData/GroundTruth/ground_truth -sub ${odir}/${oname} -abs ${odir}/${oname}_diff_tmp
- mdiff=$(fslstats -t ${odir}/${oname}_diff_tmp -k ${indir}/EddyTestData/eddyData/testMask -m)
- imrm ${odir}/${oname}_diff_tmp
- echo $mdiff > ${odir}/${oname}.ima_mean_diff.txt
- mmidiff="${mmidiff} $(echo $mdiff | awk '{sum=0; for (i=1; i<=NF; i++) { sum+= $i } print sum/NF}')"
- maxidiff="${maxidiff} $(tr ' ' '\n' <<<$mdiff | sort -r | head -n1)"
-done
-echo "mmidiff = ${mmidiff}"
-echo "maxidiff = ${maxidiff}"
-
-# Next decide on fail or pass
-
-fail=$(echo "$(echo $mmddiff | awk '{print $1}') > ${max_mean_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmddiff | awk '{print $2}') > ${max_mean_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxddiff | awk '{print $1}') > ${max_max_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxddiff | awk '{print $2}') > ${max_max_dfield_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmidiff | awk '{print $1}') > ${max_mean_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $mmidiff | awk '{print $2}') > ${max_mean_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxidiff | awk '{print $1}') > ${max_max_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-fail=$(echo "$(echo $maxidiff | awk '{print $2}') > ${max_max_ima_diff}" | bc)
-if [ $fail -gt 0 ]; then
- exit 1;
-fi
-
-exit 0
diff --git a/EddyTest/nohup.sh b/EddyTest/nohup.sh
new file mode 100755
index 0000000..a7b6766
--- /dev/null
+++ b/EddyTest/nohup.sh
@@ -0,0 +1,2 @@
+./feedsRun /vols/Data/fsldev/pyfeeds_results/EddyTestTaylor blah /vols/Data/fsldev/dataSets /opt/fmrib/fsltmp/fsl_3b4d64ae/bin
+echo Final Status $?
diff --git a/fsl_course/fdt/bedpostx_gpu/feedsRun b/fsl_course/fdt/bedpostx_gpu/feedsRun
index 2c1b1e6..a5aed2b 100755
--- a/fsl_course/fdt/bedpostx_gpu/feedsRun
+++ b/fsl_course/fdt/bedpostx_gpu/feedsRun
@@ -10,7 +10,7 @@ unset SGE_ROOT
outdir=$1
datadir=$2
-cp -r $datadir/fsl_course_data/fdt2/subj1 $outdir/
+cp -Lr $datadir/fsl_course_data/fdt2/subj1 $outdir/
bedpostx_gpu $outdir/subj1 --nf=2 --fudge=1 --bi=1000 --seed=1234
diff --git a/fsl_course/fdt1/test_topup_eddy.sh b/fsl_course/fdt1/test_topup_eddy.sh
index 8e1a4f0..8777b26 100755
--- a/fsl_course/fdt1/test_topup_eddy.sh
+++ b/fsl_course/fdt1/test_topup_eddy.sh
@@ -35,13 +35,13 @@ applytopup --imain=$outdir/nodif,$subdatadir/nodif_PA --topup=$outdir/topup_AP_P
bet $outdir/hifi_nodif $outdir/hifi_nodif_brain -m -f 0.2
# Correct for eddy currents, actually just make sure eddy ( both pre- and post- 5.0.9 naming scheme) runs and then simulate CTRL+C ( as per course instructions )
-if [ `bash -c "test -e eddy_openmp; echo \$?; exit 0"` -eq "0" ]; then
+if [ `which -s eddy_openmp; echo \$?; exit 0` -eq "0" ]; then
eddy_openmp --imain=$subdatadir/dwidata --mask=$outdir/hifi_nodif_brain_mask \
--index=$subdatadir/index.txt --acqp=$subdatadir/acqparams.txt \
--bvecs=$subdatadir/bvecs --bvals=$subdatadir/bvals \
--fwhm=0 --topup=$outdir/topup_AP_PA_b0 --flm=quadratic \
--out=$outdir/eddy_unwarped_images & sleep 10; kill -TERM $!
-elif [ `bash -c "test -e eddy; echo \$?; exit 0"` -eq "0" ]; then
+elif [ `which -s eddy; echo \$?; exit 0` -eq "0" ]; then
eddy --imain=$subdatadir/dwidata --mask=$outdir/hifi_nodif_brain_mask \
--index=$subdatadir/index.txt --acqp=$subdatadir/acqparams.txt \
--bvecs=$subdatadir/bvecs --bvals=$subdatadir/bvals \
diff --git a/fsl_course/fdt2/feedsRun6.sh b/fsl_course/fdt2/feedsRun6.sh
index 08f48a9..2af5726 100755
--- a/fsl_course/fdt2/feedsRun6.sh
+++ b/fsl_course/fdt2/feedsRun6.sh
@@ -9,11 +9,8 @@ outdir=$1
datadir=$2
subdatadir=$datadir/fsl_course_data/fdt2
-cat ${subdatadir}/LGNTracts/waypoints.txt | sed -e "s|/vols/Data/fsldev/dataSets|$datadir|g" > ${outdir}/waypoints.txt
-
-
# Tracking from voxels in MNI152 standard space
-probtrackx2 -x ${subdatadir}/subj1_2fibres.bedpostX/SEED_LGN_LEFT.nii.gz -l --onewaycondition -c 0.2 -S 2000 --steplength=0.5 -P 5000 --fibthresh=0.01 --distthresh=0.0 --sampvox=0.0 --xfm=${subdatadir}/subj1_2fibres.bedpostX/xfms/standard2diff_warp.nii.gz --invxfm=${subdatadir}/subj1_2fibres.bedpostX/xfms/diff2standard_warp.nii.gz --avoid=${subdatadir}/subj1_2fibres.bedpostX/EXCLUSION.nii.gz --stop=${subdatadir}/subj1_2fibres.bedpostX/TARGET_V1_LEFT.nii.gz --forcedir --opd -s ${subdatadir}/subj1_2fibres.bedpostX/merged -m ${subdatadir}/subj1_2fibres.bedpostX/nodif_brain_mask --dir=${outdir}/LGNTracts --waypoints=${outdir}/waypoints.txt --waycond=AND --rseed=10
+probtrackx2 -x ${subdatadir}/subj1_2fibres.bedpostX/SEED_LGN_LEFT.nii.gz -l --onewaycondition -c 0.2 -S 2000 --steplength=0.5 -P 5000 --fibthresh=0.01 --distthresh=0.0 --sampvox=0.0 --xfm=${subdatadir}/subj1_2fibres.bedpostX/xfms/standard2diff_warp.nii.gz --invxfm=${subdatadir}/subj1_2fibres.bedpostX/xfms/diff2standard_warp.nii.gz --avoid=${subdatadir}/subj1_2fibres.bedpostX/EXCLUSION.nii.gz --stop=${subdatadir}/subj1_2fibres.bedpostX/TARGET_V1_LEFT.nii.gz --forcedir --opd -s ${subdatadir}/subj1_2fibres.bedpostX/merged -m ${subdatadir}/subj1_2fibres.bedpostX/nodif_brain_mask --dir=${outdir}/LGNTracts --waypoints=${subdatadir}/subj1_2fibres.bedpostX/TARGET_V1_LEFT.nii.gz --waycond=AND --rseed=10
exit 0
diff --git a/fsl_course/feat3/feedsRun3.py b/fsl_course/feat3/feedsRun3.py
index 56e3b0e..019d6d9 100755
--- a/fsl_course/feat3/feedsRun3.py
+++ b/fsl_course/feat3/feedsRun3.py
@@ -25,7 +25,7 @@ fslDir = os.environ['FSLDIR']
# copy FSF to test location and update paths
origFeatDir = op.join(dataDir, "fsl_course_data/fmri/fluency2/fmri.feat")
-origFSF = op.join(origFeatDir, "design.fsf")
+origFSF = op.join(origFeatDir, "design_new.fsf")
origFSLDir = "/opt/fmrib/fsl-alpha"
origDataDir = "/vols/Data/fsldev/dataSets"
diff --git a/fsl_course/pnm/feedsInputs b/fsl_course/pnm/feedsInputs
new file mode 100644
index 0000000..d300b17
--- /dev/null
+++ b/fsl_course/pnm/feedsInputs
@@ -0,0 +1 @@
+fsl_course_data/fmri3/motion
\ No newline at end of file
diff --git a/fsl_course/pnm/feedsRun b/fsl_course/pnm/feedsRun
new file mode 100755
index 0000000..3d4eaff
--- /dev/null
+++ b/fsl_course/pnm/feedsRun
@@ -0,0 +1,7 @@
+#!/bin/bash -e
+
+iDir=$1
+oDir=$2
+
+pnm_evs -i ${iDir}/naughty.nii.gz -c ${iDir}/pnm/prebaked/naughty_card.txt -r ${iDir}/pnm/prebaked/naughty_resp.txt -o ${oDir}/naughty --tr=0.933 --oc=4 --or=4 --multc=2 --multr=2 --rvt=${iDir}/pnm/prebaked/naughty_rvt.txt --rvtsmooth=10 --slicetiming=${iDir}/pnm/slice_timings.txt
+pnm_evs -i ${iDir}/nice.nii.gz -c ${iDir}/pnm/prebaked/nice_card.txt -r ${iDir}/pnm/prebaked/nice_resp.txt -o ${oDir}/nice --tr=0.933 --oc=4 --or=4 --multc=2 --multr=2 --rvt=${iDir}/pnm/prebaked/nice_rvt.txt --rvtsmooth=10 --slicetiming=${iDir}/pnm/slice_timings.txt
diff --git a/mist/feedsRun b/mist/feedsRun
index 31052d5..a2f712a 100755
--- a/mist/feedsRun
+++ b/mist/feedsRun
@@ -11,6 +11,7 @@ datadir=$2/mistData/
pushd $outdir > /dev/null
cp -LR $datadir/* .
+chmod -R u+w .
mist_1_train hippocampus thalamus red_nucleus
mist_2_fit
diff --git a/mist/feedsRun~ b/mist/feedsRun~
deleted file mode 100755
index 2db5b76..0000000
--- a/mist/feedsRun~
+++ /dev/null
@@ -1,26 +0,0 @@
-#!/bin/bash
-
-
-unset SGE_ROOT
-
-set -e
-
-outdir=$1
-datadir=$2/mist/
-
-
-pushd $outdir > /dev/null
-
-cp -Lr $datadir/* .
-
-mist_1_train
-mist_2_fit
-
-# save a bit of space - remove input data
-for subj in ./???; do
- rm $subj/T1.nii.gz
- rm $subj/T1_brain.nii.gz
- rm $subj/FA.nii.gz
-done
-
-popd > /dev/null
--
GitLab