From 9f83c6f9f479e29f43decd28bde5d39abc346048 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Wed, 7 Sep 2022 16:26:46 +0100
Subject: [PATCH 1/5] RF: avoid pointless loops depending on the value of
isComplex
---
generalio.cc | 15 ++++++++++-----
1 file changed, 10 insertions(+), 5 deletions(-)
diff --git a/generalio.cc b/generalio.cc
index ad8fc0e..081d1fa 100644
--- a/generalio.cc
+++ b/generalio.cc
@@ -576,11 +576,16 @@ int read_complexvolume(volume<float>& realvols, volume<float>& imagvols,
extensions.clear();
}
- for ( size_t voxel=0;voxel<volsize && !isComplex;voxel++) //Only copy real buffer for non-complex
- rbuffer[voxel]=((float*)buffer)[voxel];
- for ( size_t voxel=0;voxel<volsize && isComplex;voxel++) {
- rbuffer[voxel]=((float *)buffer)[2*voxel];
- ibuffer[voxel]=((float *)buffer)[2*voxel+1];
+ //Only copy real buffer for non-complex
+ if (!isComplex) {
+ for ( size_t voxel=0;voxel<volsize && !isComplex;voxel++)
+ rbuffer[voxel]=((float*)buffer)[voxel];
+ }
+ else {
+ for ( size_t voxel=0;voxel<volsize && isComplex;voxel++) {
+ rbuffer[voxel]=((float *)buffer)[2*voxel];
+ ibuffer[voxel]=((float *)buffer)[2*voxel+1];
+ }
}
realvols.reinitialize(sx,sy,sz,st,rbuffer,true);
imagvols.reinitialize(sx,sy,sz,st,ibuffer,true);
--
GitLab
From 209e5e08bc20ed0e44106cde159f21dfc8d65cfd Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Wed, 7 Sep 2022 16:27:24 +0100
Subject: [PATCH 2/5] BF: Actually copy orientation info from complex input
file to real/imag outputs, instead of just doing L/R flip
---
generalio.cc | 14 ++++----------
1 file changed, 4 insertions(+), 10 deletions(-)
diff --git a/generalio.cc b/generalio.cc
index 081d1fa..c5e770d 100644
--- a/generalio.cc
+++ b/generalio.cc
@@ -593,16 +593,10 @@ int read_complexvolume(volume<float>& realvols, volume<float>& imagvols,
realvols.setdims(niihdr.pixdim[1],niihdr.pixdim[2],niihdr.pixdim[3],niihdr.pixdim[4]);
imagvols.setdims(niihdr.pixdim[1],niihdr.pixdim[2],niihdr.pixdim[3],niihdr.pixdim[4]);
- // swap to Radiological when necessary
- if ( NiftiGetLeftRightOrder(niihdr) != FSL_RADIOLOGICAL ) {
- realvols.RadiologicalFile = false;
- realvols.makeradiological();
- imagvols.RadiologicalFile = false;
- imagvols.makeradiological();
- } else {
- realvols.RadiologicalFile = true;
- imagvols.RadiologicalFile = true;
- }
+
+ set_volume_properties(niihdr,realvols);
+ set_volume_properties(niihdr,imagvols);
+
delete [] buffer;
return 0;
}
--
GitLab
From b932d06df9ba1063ff328fc1d87c2ff4dbdbb4fa Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Wed, 7 Sep 2022 17:14:57 +0100
Subject: [PATCH 3/5] RF: remove isComplex check in for-loop conditionals
---
generalio.cc | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/generalio.cc b/generalio.cc
index c5e770d..7352e31 100644
--- a/generalio.cc
+++ b/generalio.cc
@@ -578,13 +578,13 @@ int read_complexvolume(volume<float>& realvols, volume<float>& imagvols,
//Only copy real buffer for non-complex
if (!isComplex) {
- for ( size_t voxel=0;voxel<volsize && !isComplex;voxel++)
- rbuffer[voxel]=((float*)buffer)[voxel];
+ for ( size_t voxel=0; voxel < volsize; voxel++)
+ rbuffer[voxel] = ((float*)buffer)[voxel];
}
else {
- for ( size_t voxel=0;voxel<volsize && isComplex;voxel++) {
- rbuffer[voxel]=((float *)buffer)[2*voxel];
- ibuffer[voxel]=((float *)buffer)[2*voxel+1];
+ for ( size_t voxel=0; voxel < volsize; voxel++) {
+ rbuffer[voxel] = ((float *)buffer)[2*voxel];
+ ibuffer[voxel] = ((float *)buffer)[2*voxel+1];
}
}
realvols.reinitialize(sx,sy,sz,st,rbuffer,true);
--
GitLab
From 33a28d4f2b0a6662fce8fe21971bb58921793126 Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 8 Sep 2022 15:37:59 +0100
Subject: [PATCH 4/5] BF: Accidentally removed l/r flips in read_complexvolume
---
generalio.cc | 3 +++
test/{test.cc => test_volume_constructor.cc} | 0
2 files changed, 3 insertions(+)
rename test/{test.cc => test_volume_constructor.cc} (100%)
diff --git a/generalio.cc b/generalio.cc
index 7352e31..55cf824 100644
--- a/generalio.cc
+++ b/generalio.cc
@@ -597,6 +597,9 @@ int read_complexvolume(volume<float>& realvols, volume<float>& imagvols,
set_volume_properties(niihdr,realvols);
set_volume_properties(niihdr,imagvols);
+ if (!realvols.RadiologicalFile) realvols.makeradiological();
+ if (!imagvols.RadiologicalFile) imagvols.makeradiological();
+
delete [] buffer;
return 0;
}
diff --git a/test/test.cc b/test/test_volume_constructor.cc
similarity index 100%
rename from test/test.cc
rename to test/test_volume_constructor.cc
--
GitLab
From 397d9a14b60afd8604caa2845258c574550ef6ac Mon Sep 17 00:00:00 2001
From: Paul McCarthy <pauldmccarthy@gmail.com>
Date: Thu, 8 Sep 2022 15:38:56 +0100
Subject: [PATCH 5/5] TEST: Re-arrange test code, and add unit tests for
read/save_complexvolume
---
test/Makefile | 11 +++-
test/feedsRun | 59 ++++++++++++++++--
test/test_complexvolume.cc | 103 ++++++++++++++++++++++++++++++++
test/test_main.cc | 4 ++
test/test_volume_constructor.cc | 9 ++-
5 files changed, 177 insertions(+), 9 deletions(-)
create mode 100644 test/test_complexvolume.cc
create mode 100644 test/test_main.cc
diff --git a/test/Makefile b/test/Makefile
index 621c2b7..7d85699 100644
--- a/test/Makefile
+++ b/test/Makefile
@@ -5,7 +5,7 @@ PROJNAME = test-newimage
TESTXFILES = test-newimage
LIBS = -lfsl-newimage -lfsl-miscmaths -lfsl-NewNifti \
- -lfsl-cprob -lfsl-utils -lfsl-znz
+ -lfsl-cprob -lfsl-utils -lfsl-znz -lboost_unit_test_framework
# The test program can be run against
# an in-source checkout, or against
@@ -30,8 +30,13 @@ endif
all: ${TESTXFILES}
+OBJS = test_main.o test_volume_constructor.o test_complexvolume.o
+
# We add -I.., -L.., -Wl,-rpath so that
# in-source builds take precedence over
# $FSLDEVDIR/$FSLDIR
-test-newimage: test.cc
- $(CXX) -I.. ${CXXFLAGS} -o $@ $< -L.. ${RPATH} ${LDFLAGS}
+%.o: %.cc
+ $(CXX) -I.. ${CXXFLAGS} -c -o $@ $<
+
+test-newimage: ${OBJS}
+ $(CXX) -o $@ $^ -L.. ${RPATH} ${LDFLAGS}
diff --git a/test/feedsRun b/test/feedsRun
index dcf6266..b8d681d 100755
--- a/test/feedsRun
+++ b/test/feedsRun
@@ -1,6 +1,57 @@
-#!/usr/bin/env bash
+#!/usr/bin/env python
-set -e
+import numpy as np
+import nibabel as nib
+import subprocess as sp
+import os.path as op
+import shlex
+import sys
+import traceback
-make
-./test-newimage -l all
+from fsl.utils.tempdir import tempdir
+from fsl.data.image import addExt
+
+
+def sprun(cmd):
+ print(f'RUN {cmd}')
+ sp.run(shlex.split(cmd), check=True)
+
+
+# radio=[True|False] -> whether or not
+# to add a flip on the affine X axis
+def create_test_input_data(seed=1):
+
+ np.random.seed(seed)
+
+ rneuro = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
+ ineuro = np.random.randint(0, 100, (10, 10, 10)).astype(np.float32)
+ rradio = np.flip(rneuro, 0)
+ iradio = np.flip(ineuro, 0)
+
+ # test_complexvolume.cc has these
+ # same values hard-coded, so if
+ # you change one, you must change
+ # the other.
+ affneuro = np.diag([ 3, 3, 3, 1])
+ affradio = np.diag([-3, 3, 3, 1])
+ affneuro[:3, 3] = [10, 20, 30]
+ affradio[:3, 3] = [37, 20, 30]
+
+ cneuro = rneuro + ineuro * 1j
+ cradio = rradio + iradio * 1j
+
+ cneuro = nib.Nifti1Image(cneuro, affneuro)
+ cradio = nib.Nifti1Image(cradio, affradio)
+
+ cneuro.set_qform(*cneuro.get_sform(coded=True))
+ cradio.set_qform(*cradio.get_sform(coded=True))
+
+ cneuro.to_filename(f'complex_neuro.nii.gz')
+ cradio.to_filename(f'complex_radio.nii.gz')
+
+
+if __name__ == '__main__':
+ create_test_input_data()
+ sprun('make clean')
+ sprun('make')
+ sprun('./test-newimage -l test_suite')
diff --git a/test/test_complexvolume.cc b/test/test_complexvolume.cc
new file mode 100644
index 0000000..4b95e94
--- /dev/null
+++ b/test/test_complexvolume.cc
@@ -0,0 +1,103 @@
+#define EXPOSE_TREACHEROUS
+
+#include "newimage/newimageall.h"
+#include "armawrap/newmat.h"
+#include <stdlib.h>
+#include <vector>
+#include <boost/test/unit_test.hpp>
+
+
+BOOST_AUTO_TEST_SUITE(test_complexvolume)
+
+
+using namespace NEWIMAGE;
+using namespace NEWMAT;
+using namespace std;
+
+// Test data used by these tests is generated in the
+// feedsRun script before execution.
+
+// Test that neurologically and radiologically stored images
+// (which are otherwise identical) are loaded correctly
+BOOST_AUTO_TEST_CASE(read_complexvolume_handle_neuro_radio)
+{
+ complexvolume neuro, radio;
+
+ BOOST_CHECK(read_complexvolume(neuro, string("complex_neuro")) == 0);
+ BOOST_CHECK(read_complexvolume(radio, string("complex_radio")) == 0);
+
+ BOOST_CHECK(neuro.xsize() == radio.xsize());
+ BOOST_CHECK(neuro.ysize() == radio.ysize());
+ BOOST_CHECK(neuro.zsize() == radio.zsize());
+
+ for (int x = 0; x < neuro.xsize(); x++) {
+ for (int y = 0; y < neuro.ysize(); y++) {
+ for (int z = 0; z < neuro.zsize(); z++) {
+ BOOST_CHECK(neuro.re(x, y, z) == radio.re(x, y, z));
+ BOOST_CHECK(neuro.im(x, y, z) == radio.im(x, y, z));
+ }}}
+}
+
+
+BOOST_AUTO_TEST_CASE(read_save_complexvolume_round_trip)
+{
+ volume<float> rneuro, ineuro, rradio, iradio,
+ rneuro2, ineuro2, rradio2, iradio2;
+
+ // defined in feedsRun
+ Matrix neuromat(4, 4), radiomat(4, 4);
+ neuromat << 3 << 0 << 0 << 10
+ << 0 << 3 << 0 << 20
+ << 0 << 0 << 3 << 30
+ << 0 << 0 << 0 << 1;
+ radiomat << -3 << 0 << 0 << 37
+ << 0 << 3 << 0 << 20
+ << 0 << 0 << 3 << 30
+ << 0 << 0 << 0 << 1;
+
+ BOOST_CHECK(read_complexvolume(rneuro, ineuro, string("complex_neuro")) == 0);
+ BOOST_CHECK(read_complexvolume(rradio, iradio, string("complex_radio")) == 0);
+ BOOST_CHECK(save_complexvolume(rneuro, ineuro, string("complex_neuro2")) == 0);
+ BOOST_CHECK(save_complexvolume(rradio, iradio, string("complex_radio2")) == 0);
+ BOOST_CHECK(read_complexvolume(rneuro2, ineuro2, string("complex_neuro2")) == 0);
+ BOOST_CHECK(read_complexvolume(rradio2, iradio2, string("complex_radio2")) == 0);
+
+ BOOST_CHECK(rneuro.sform_code() == rneuro2.sform_code());
+ BOOST_CHECK(rneuro.qform_code() == rneuro2.qform_code());
+ BOOST_CHECK(ineuro.sform_code() == ineuro2.sform_code());
+ BOOST_CHECK(ineuro.qform_code() == ineuro2.qform_code());
+ BOOST_CHECK(rradio.sform_code() == rradio2.sform_code());
+ BOOST_CHECK(rradio.qform_code() == rradio2.qform_code());
+
+ // All volumes are presented as radiological - the
+ // L/R flip on neurological volumes is also applied
+ // to the s/qform affines, so we compare all volumes
+ // against the known radiological affine.
+ BOOST_CHECK(rradio .sform_mat() == radiomat);
+ BOOST_CHECK(rradio2.sform_mat() == radiomat);
+ BOOST_CHECK(iradio .sform_mat() == radiomat);
+ BOOST_CHECK(iradio2.sform_mat() == radiomat);
+ BOOST_CHECK(rneuro .sform_mat() == radiomat);
+ BOOST_CHECK(rneuro2.sform_mat() == radiomat);
+ BOOST_CHECK(ineuro .sform_mat() == radiomat);
+ BOOST_CHECK(ineuro2.sform_mat() == radiomat);
+ BOOST_CHECK(rradio .qform_mat() == radiomat);
+ BOOST_CHECK(rradio2.qform_mat() == radiomat);
+ BOOST_CHECK(iradio .qform_mat() == radiomat);
+ BOOST_CHECK(iradio2.qform_mat() == radiomat);
+ BOOST_CHECK(rneuro .qform_mat() == radiomat);
+ BOOST_CHECK(rneuro2.qform_mat() == radiomat);
+ BOOST_CHECK(ineuro .qform_mat() == radiomat);
+ BOOST_CHECK(ineuro2.qform_mat() == radiomat);
+
+ for (int x = 0; x < rneuro.xsize(); x++) {
+ for (int y = 0; y < rneuro.ysize(); y++) {
+ for (int z = 0; z < rneuro.zsize(); z++) {
+ BOOST_CHECK(rneuro(x, y, z) == rneuro2(x, y, z));
+ BOOST_CHECK(ineuro(x, y, z) == ineuro2(x, y, z));
+ BOOST_CHECK(rradio(x, y, z) == rradio2(x, y, z));
+ BOOST_CHECK(iradio(x, y, z) == iradio2(x, y, z));
+ }}}
+}
+
+BOOST_AUTO_TEST_SUITE_END()
diff --git a/test/test_main.cc b/test/test_main.cc
new file mode 100644
index 0000000..58fcf6c
--- /dev/null
+++ b/test/test_main.cc
@@ -0,0 +1,4 @@
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_MODULE newimage
+#include <boost/test/unit_test.hpp>
diff --git a/test/test_volume_constructor.cc b/test/test_volume_constructor.cc
index 08a9395..c06a804 100644
--- a/test/test_volume_constructor.cc
+++ b/test/test_volume_constructor.cc
@@ -1,8 +1,11 @@
#include "newimage/newimageall.h"
#include <stdlib.h>
#include <vector>
-#define BOOST_TEST_MODULE newimage
-#include <boost/test/included/unit_test.hpp>
+#include <boost/test/unit_test.hpp>
+
+
+BOOST_AUTO_TEST_SUITE(test_volume_constructor)
+
using namespace NEWIMAGE;
@@ -76,3 +79,5 @@ BOOST_AUTO_TEST_CASE(volume_4D_reinitialize)
BOOST_TEST(v(xi, yi, zi, ti) == i);
}
}
+
+BOOST_AUTO_TEST_SUITE_END()
--
GitLab