Basic unit test for vecreg

unit_tests/fdt/vecreg/feedsInputs

+fsl_course_data/fdt2

unit_tests/fdt/vecreg/feedsRun

+#!/usr/bin/env fslpython
+
+import os
+import sys
+
+sys.path.append('/home/paulmc/Projects/fslpy/')
+
+from pathlib import Path
+
+import numpy as np
+
+from fsl.data.image import Image
+from fsl.wrappers import vecreg, concatxfm, invxfm, convertwarp
+
+from pyfeeds.evaluate import cmpVectorArrays
+
+
+def main():
+    fsldir        = Path(os.environ['FSLDIR'])
+    outdir        = Path(sys.argv[1])
+    indir         = Path(sys.argv[2])/'fsl_course_data'/'fdt2'
+    dyads         = indir/'diffusion.bedpostX'/'dyads1.nii.gz'
+    std           = fsldir/'data'/'standard'/'MNI152_T1_2mm.nii.gz'
+    str2std       = indir/'diffusion.bedpostX'/'xfms'/'str2standard.mat'
+    diff2str      = indir/'structural'/'xfms'/'diff2str.mat'
+    diff2std_warp = indir/'structural'/'xfms'/'diff2standard_warp.nii.gz'
+
+    # Create diff2std/std2diff affines
+    std2diff = outdir/'std2diff.mat'
+    diff2std = outdir/'diff2std.mat'
+    concatxfm(diff2str, str2std, diff2std)
+    invxfm(diff2std, std2diff)
+
+    # create a copy of the diff2std warp
+    # which does not encode the diff2std
+    # affine
+    diff2std_warp_no_premat = outdir/'diff2std_warp_no_premat.nii.gz'
+    convertwarp(diff2std_warp_no_premat, std,
+                warp1=diff2std_warp,
+                premat=std2diff,
+                relout=True)
+
+    # Transform dyads with the full diff2std
+    # warp (which encodes the diff->standard
+    # affine transform) - this will be used
+    # as a benchmark
+    dyads_std = outdir/'dyads1_std.nii.gz'
+    vecreg(dyads, dyads_std, std, warpfield=diff2std_warp)
+
+    # Transform dyads with the raw diff2std
+    # warp, and specify the diff2std affine
+    # as premat. The result should match
+    # that from the first vecreg call above.
+    dyads_std_with_premat = outdir/'dyads1_std_with_premat.nii.gz'
+    vecreg(dyads, dyads_std_with_premat, std,
+           warpfield=diff2std_warp_no_premat,
+           premat=diff2std)
+
+    # compare the result
+    dyads1 = Image(dyads_std).data
+    dyads2 = Image(dyads_std_with_premat).data
+    diff   = cmpVectorArrays(dyads1, dyads2)
+    assert np.isclose(diff, 0, atol=1e-4), diff
+
+    # Create a nonsense warp field with a postmat
+    diff2std_warp_postmat = outdir/'diff2std_warp_postmat.nii.gz'
+    convertwarp(diff2std_warp_postmat, std,
+                warp1=diff2std_warp_no_premat,
+                premat=diff2std,
+                postmat=std2diff)
+
+    # Warp dyads with the nonsense warp
+    dyads_postmat1 = outdir/'dyads_postmat1.nii.gz'
+    vecreg(dyads, dyads_postmat1, std,
+           warpfield=diff2std_warp_postmat)
+
+    # And with the original warp+premat+postmat
+    dyads_postmat2 = outdir/'dyads_postmat2.nii.gz'
+    vecreg(dyads, dyads_postmat2, std,
+           warpfield=diff2std_warp_no_premat,
+           premat=diff2std, postmat=std2diff)
+
+    # And check that the results are equivalent
+    dyads1 = Image(dyads_postmat1).data
+    dyads2 = Image(dyads_postmat2).data
+    diff   = cmpVectorArrays(dyads1, dyads2)
+    assert np.isclose(diff, 0, atol=1e-4), diff
+
+
+if __name__ == '__main__':
+    main()