From a9c3e1767bf44417788a753c50f5d9c86b2ae863 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 30 Jul 2020 11:55:12 +0100
Subject: [PATCH] ENH: Rewrote fslcreatehd test in python, expanded tests to
 handle 4th dim of (0, 1), and to handle modification of existing image

---
 unit_tests/fslcreatehd/feedsRun | 265 +++++++++++++++++++++-----------
 1 file changed, 176 insertions(+), 89 deletions(-)

diff --git a/unit_tests/fslcreatehd/feedsRun b/unit_tests/fslcreatehd/feedsRun
index 8563ff9..136e138 100755
--- a/unit_tests/fslcreatehd/feedsRun
+++ b/unit_tests/fslcreatehd/feedsRun
@@ -1,96 +1,183 @@
-#!/bin/bash
+#!/usr/bin/env python
 
-set -e
+import os
+import sys
+import subprocess as sp
 
+import numpy   as np
+import nibabel as nib
 
-function neq() {
-    n1=$1
-    n2=$2
-    result=`echo "$n1 == $n2" | bc`
 
-    if [[ "$result" == "0" ]]; then
-        echo "$n1 != $n2"
-        return 0;
-    fi
-
-    return 1;
-}
-
-function fslcreatehd_test() {
-
-    dim0=$1;    shift
-    dim1=$1;    shift
-    dim2=$1;    shift
-    dim3=$1;    shift
-    dim4=$1;    shift
-    pixdim1=$1; shift
-    pixdim2=$1; shift
-    pixdim3=$1; shift
-    pixdim4=$1; shift
-    originx=$1; shift
-    originy=$1; shift
-    originz=$1; shift
-    dtype=$1;   shift
-    fname=$1;   shift
-
-    args="$dim0 $dim1 $dim2 $dim3 $pixdim1 $pixdim2 $pixdim3 $pixdim4 $originx  $originy  $originz $dtype $fname"
-
-    originx=`echo " $originx * $pixdim1" | bc`
-    originy=`echo "-$originy * $pixdim2"  | bc`
-    originz=`echo "-$originz * $pixdim3"  | bc`
-
-    got_dim0=`   fslhd $fname | egrep "^dim0"      | tr -s '\t' | cut -f 2`
-    got_dim1=`   fslhd $fname | egrep "^dim1"      | tr -s '\t' | cut -f 2`
-    got_dim2=`   fslhd $fname | egrep "^dim2"      | tr -s '\t' | cut -f 2`
-    got_dim3=`   fslhd $fname | egrep "^dim3"      | tr -s '\t' | cut -f 2`
-    got_dim4=`   fslhd $fname | egrep "^dim4"      | tr -s '\t' | cut -f 2`
-    got_pixdim1=`fslhd $fname | egrep "^pixdim1"   | tr -s '\t' | cut -f 2`
-    got_pixdim2=`fslhd $fname | egrep "^pixdim2"   | tr -s '\t' | cut -f 2`
-    got_pixdim3=`fslhd $fname | egrep "^pixdim3"   | tr -s '\t' | cut -f 2`
-    got_pixdim4=`fslhd $fname | egrep "^pixdim4"   | tr -s '\t' | cut -f 2`
-    got_dtype=`  fslhd $fname | egrep "^datatype"  | tr -s '\t' | cut -f 2`
-    got_originx=`fslhd $fname | egrep "^sto_xyz:1"              | cut -d ' ' -f 4`
-    got_originy=`fslhd $fname | egrep "^sto_xyz:2"              | cut -d ' ' -f 4`
-    got_originz=`fslhd $fname | egrep "^sto_xyz:3"              | cut -d ' ' -f 4`
-
-    if neq $dim0    $got_dim0;    then echo $args; return 1; fi
-    if neq $dim1    $got_dim1;    then echo $args; return 1; fi
-    if neq $dim2    $got_dim2;    then echo $args; return 1; fi
-    if neq $dim3    $got_dim3;    then echo $args; return 1; fi
-    if neq $dim4    $got_dim4;    then echo $args; return 1; fi
-    if neq $pixdim1 $got_pixdim1; then echo $args; return 1; fi
-    if neq $pixdim2 $got_pixdim2; then echo $args; return 1; fi
-    if neq $pixdim3 $got_pixdim3; then echo $args; return 1; fi
-    if neq $pixdim4 $got_pixdim4; then echo $args; return 1; fi
-    if neq $dtype   $got_dtype;   then echo $args; return 1; fi
-    if neq $originx $got_originx; then echo $args; return 1; fi
-    if neq $originy $got_originy; then echo $args; return 1; fi
-    if neq $originz $got_originz; then echo $args; return 1; fi
-
-    return 0
+# 2=char, 4=short, 8=int, 16=float, 64=double
+DTYPE_MAPPING = {
+    2  : np.uint8,
+    4  : np.int16,
+    8  : np.int32,
+    16 : np.float32,
+    64 : np.float64
 }
 
 
-tests="
-5 5 5 1     2   2   2    1     0 0 0     2
-5 5 5 1     2   2   2    1     0 0 0     4
-5 5 5 1     2   2   2    1     0 0 0     8
-5 5 5 1     2   2   2    1     0 0 0     16
-5 5 5 1     2   2   2    1     1 2 3     64
-5 5 5 1     0.5 1.5 1.25 1     0 0 0     2
-5 5 5 1     0.5 1.5 1.25 1.5   0 0 0     2
-5 5 5 1     0.5 1.5 1.25 0.5   0 0 0     2
-
-30 30 30 5   5 10 3 5   10 20 10  2
-"
-
-echo "$tests" | while IFS= read -r test; do
-  if [[ $test =~ ^\ *$ ]]; then
-    continue
-  fi
-
-  rm `imglob -extension fslcreatehd_test_img` &> /dev/null || true
-  fslcreatehd        $test fslcreatehd_test_img
-  fslcreatehd_test 4 $test fslcreatehd_test_img
-
-done
+def validate_result(imgfile, expshape, exppixdim, exporigin, expdtype):
+
+    img       = nib.load(imgfile)
+    hdr       = img.header
+    expshape  = list(expshape)
+    exppixdim = list(exppixdim)
+    exporigin = list(exporigin)
+
+    if (len(expshape) == 3) or \
+       expshape[3] in (0, 1):
+        expndims = 3
+    else:
+        expndims = 4
+
+    expaffine         = np.diag(exppixdim[:3] + [1])
+    expaffine[0,  0] *= -1
+    expaffine[:3, 3]  =  exporigin
+    expaffine[0,  3] *=  exppixdim[0]
+    expaffine[1,  3] *= -exppixdim[1]
+    expaffine[2,  3] *= -exppixdim[2]
+
+    assert len(img.shape)        == expndims
+    assert list(img.shape)       == expshape[ :expndims]
+    assert list(hdr.get_zooms()) == exppixdim[:expndims]
+    assert img.get_data_dtype()  == expdtype
+    assert np.all(np.isclose(img.affine, expaffine))
+
+
+def create_image(shape, pixdim, origin, dtype):
+    pixdim        = list(pixdim)
+    affine        = np.diag(pixdim[:3] + [1])
+    affine[:3, 3] = origin
+    data          = np.random.randint(1, 100, shape).astype(dtype)
+    hdr           = nib.Nifti1Header()
+
+    hdr.set_data_dtype(dtype)
+    hdr.set_data_shape(shape)
+    hdr.set_zooms(pixdim[:len(shape)])
+
+    return nib.Nifti1Image(data, affine, hdr)
+
+
+def test_new_file():
+    tests = [
+
+        # 4th dim of size 0 or 1 should
+        # result in a 3D image (this
+        # is coded in validate_result)
+        ' 5  5  5 0     2    2   2    0      0  0  0    2',
+        ' 5  5  5 0     2    2   2    1      0  0  0    2',
+        ' 5  5  5 1     2    2   2    1      0  0  0    2',
+        ' 5  5  5 1     2    2   2    1      0  0  0    4',
+        ' 5  5  5 1     2    2   2    1      0  0  0    8',
+        ' 5  5  5 1     2    2   2    1      0  0  0    16',
+        ' 5  5  5 1     2    2   2    1      1  2  3    64',
+        ' 5  5  5 1     0.5  1.5 1.25 1      0  0  0    2',
+        ' 5  5  5 1     0.5  1.5 1.25 1.5    0  0  0    2',
+        ' 5  5  5 1     0.5  1.5 1.25 0.5    0  0  0    2',
+        '30 30 30 5     5   10   3    5     10 20 10    2',
+    ]
+
+    for test in tests:
+        args    = list(test.split())
+        imgfile = 'image.nii.gz'
+
+        sp.run(['fslcreatehd'] + args + [imgfile])
+
+        args      = [float(a) for a in args]
+        expshape  = args[:4]
+        exppixdim = args[4:8]
+        exporigin = args[8:11]
+        expdtype  = DTYPE_MAPPING[args[11]]
+
+        try:
+            validate_result(
+                imgfile, expshape, exppixdim, exporigin, expdtype)
+        finally:
+            os.remove(imgfile)
+
+
+def test_existing_file():
+
+    # each test is a tuple containing:
+    #  - dtype_of_existing_image
+    #  - shape_of_existing_image
+    #  - exp_shape (set to None if image should not be overwritten)
+    #  - fslcreatehd_args
+    tests = [
+        # 4th dim of 0 or 1 should result in a 3D image
+        (2, (5, 5, 5), None, '5 5 5 0  2 2 2 1  0 0 0 2'),
+        (2, (5, 5, 5), None, '5 5 5 0  2 2 2 1  1 2 3 2'),
+        (2, (5, 5, 5), None, '5 5 5 1  2 2 2 2  0 0 0 2'),
+
+        # dtype arg should be ignored for existing images
+        (2, (5, 5, 5), None, '5 5 5 0  2 2 2 1  0 0 0 4'),
+        (4, (5, 5, 5), None, '5 5 5 0  2 2 2 1  0 0 0 2'),
+
+        # data should be overwritten if any
+        # of the first 3 dims are different,
+        # or if a 4th dim is specified
+        (2, (5, 5, 5), (5, 3, 5),    '5 3 5 0  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5), (5, 3, 5),    '5 3 5 0  2 2 2 2  0 0 0 4'),
+        (2, (5, 5, 5), (5, 3, 5),    '5 3 5 1  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5), (5, 3, 5),    '5 3 5 1  2 2 2 2  1 2 3 2'),
+        (2, (5, 5, 5), (5, 3, 5, 5), '5 3 5 5  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5), (5, 3, 5, 5), '5 3 5 5  2 2 2 2  1 2 3 2'),
+
+        # 4D - same rules apply
+        (2, (5, 5, 5, 5), None, '5 5 5 5  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), None, '5 5 5 5  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), None, '5 5 5 5  2 2 2 2  0 0 0 4'),
+        (4, (5, 5, 5, 5), None, '5 5 5 5  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), None, '5 5 5 5  2 2 2 2  1 2 3 2'),
+
+        # data overwritten if nelements
+        # change
+        (2, (5, 5, 5, 5), (5, 5, 5),    '5 5 5 0  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 5, 5),    '5 5 5 1  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 3, 5),    '5 3 5 0  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 3, 5),    '5 3 5 1  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 3, 5),    '5 3 5 1  2 2 2 2  1 2 3 2'),
+        (2, (5, 5, 5, 5), (5, 5, 5, 2), '5 5 5 2  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 3, 5, 5), '5 3 5 5  2 2 2 2  0 0 0 2'),
+        (2, (5, 5, 5, 5), (5, 3, 5, 5), '5 3 5 5  2 2 2 2  1 2 3 2'),
+    ]
+
+    for (dtype, shape, expshape, args) in tests:
+        imgfile = 'image.nii.gz'
+        dtype   = DTYPE_MAPPING[dtype]
+        args    = list(args.split())
+        img     = create_image(shape, (1, 1, 1, 1), (0, 0, 0), dtype)
+
+        img.to_filename(imgfile)
+
+        sp.run(['fslcreatehd'] + args + [imgfile])
+
+        checkdata = expshape is None
+
+        if expshape is None:
+            expshape = shape
+
+        args      = [float(a) for a in args]
+        exppixdim = args[4:8]
+        exporigin = args[8:11]
+
+        try:
+            validate_result(imgfile, expshape, exppixdim, exporigin, dtype)
+
+            if checkdata:
+                data     = np.asanyarray(img.dataobj)
+                datacopy = np.asanyarray(nib.load(imgfile).dataobj)
+                assert np.all(data == datacopy)
+
+        finally:
+            os.remove(imgfile)
+
+
+def main():
+    test_new_file()
+    test_existing_file()
+
+if __name__ == '__main__':
+    sys.exit(main())
-- 
GitLab