A more flexible alternative to TO_BE_IGNORED macro (#1207)

* implemented a more flexible alternative to TO_BE_IGNORED macro

* moved IgnoreMask to optional arguments of set_data etc.

* defined ISMRMRD acquisition flags in Gadgetron.py

* implemented acquisition ignoring helper class IgnoreMask in Python

* provided example of user-set acquisition ignoring in acquisition_data.py

* changed default IgnoreMask for reading from file to ~13bffff

* have to use int for ignore mask instead of size_t as SWIG rejects int64

* [ci skip] small amendment in cGT_ISMRMRDAcquisitionsFromFile

* changed default value for all_ in AcquisitionData.__init__ to True

* put debug prints in test_pSynergistic.py

* trying reading MR acquisitions with all_=False

* trying reading MR acquisitions with all_=False in noncartesian_recon.py

* trying literal_eval instead of eval to placate Codacy

* small edit in test_pSynergistic.py to trigger GitHub Actions

* small corrections in MR/acquisition_data.py

* switched to 64-bit acquisition ignore masks

* made some white space adjustments

* switched to employing MRAcquisitionData attribute ignore_mask

* implemented set_ignore_mask method for MR AcquisitionData Python class

* covered all acquisitions ignored cases, fixes #1220

* corrected previous commit

* got rid of object_mask

* attended to the reviewer comments and suggestions

* Update build-test.yml

remove `-e` from the ctest line preventing unit test to complete.

* interfaced MRAcquisitionData::ignore_mask() into Python

* implemented setter for smoothing (convolution) kernel size

* updated CHANGES.md

* replaced C-style casts with reinterpret_cast in cgadgetron.cpp

* if this commit builds, the source of build failures is localised

* corrected AcquisitionData.ignore_mask in Gadgetron.py

* added missing initialisations of MRAcquisitionData.ignore_mask_ in constructors

* switched to defalt value False for all_ in AcquisitionData constructor

* smoothness->smoothing_iterations, smth_kernel_size->conv_kernel_halfsize

* [ci skip] removed some commented-out and unused stuff

* added method size() to C++ ImageData class

---------

Co-authored-by: Edoardo Pasca <[email protected]>

.github/workflows/build-test.yml

@@ -121,7 +121,7 @@ jobs:
     - name: tests
       shell: bash
       run:
-          bash -ev ${GITHUB_WORKSPACE}/SIRF-SuperBuild/docker/ctest_sirf.sh
+          bash -v ${GITHUB_WORKSPACE}/SIRF-SuperBuild/docker/ctest_sirf.sh
 
     - name: Coverage
       shell: bash

CHANGES.md

@@ -22,6 +22,9 @@
   - Several allocate methods in STIR.py, Gadgetron.py and Reg.py are replaced with just one allocate in DataContainer class that does not copy data between Python and C++.
   - `return None` in the method `Datacontainer.shape()` replaced with more Pythonesque `return (0,)`.
 
+* MR
+  - Handling of "irregular" ISMRMRD acquisitions hit what appears to be a bug in recent Gadgetron (HEAD detached at 0670db84). A quick fix applied, Gadgetron issue to be raised.
+
 ## v3.5.0
 
 * GitHub Action: remove temporarily the Ubuntu 20.04 build, #1178

examples/Python/MR/acquisition_data.py

@@ -11,6 +11,8 @@
   -p <path>, --path=<path>    path to data files, defaults to data/examples/MR
                               subfolder of SIRF root folder
   -s <slcs>, --slices=<slcs>  max number of slices to display [default: 8]
+  -i <flag>, --ignore=<flag>  ignore acquisitions with non-zero value of <flag>
+                              (try e.g. --ignore=19 and --file=grappa2_6rep.h5)
   -e <engn>, --engine=<engn>  reconstruction engine [default: Gadgetron]
   --non-interactive           do not show plots
 '''
@@ -38,6 +40,7 @@
 from docopt import docopt
 args = docopt(__doc__, version=__version__)
 
+import ast
 from sirf.Utilities import error, examples_data_path, existing_filepath
 
 # import MR engine module
@@ -50,6 +53,10 @@
 if data_path is None:
     data_path = examples_data_path('MR')
 slcs = int(args['--slices'])
+to_be_ignored = args['--ignore']
+ignore_mask = mr.IgnoreMask(0)
+if to_be_ignored is not None:
+    ignore_mask.ignore(ast.literal_eval(to_be_ignored))
 show_plot = not args['--non-interactive']
 
 def main():
@@ -58,7 +65,7 @@ def main():
     input_file = existing_filepath(data_path, data_file)
 
     # acquisition data will be read from an HDF file input_file
-    acq_data = mr.AcquisitionData(input_file, True)
+    acq_data = mr.AcquisitionData(input_file)
 
     # the raw k-space data is a list of different readouts
     # of different data type (e.g. noise correlation data, navigator data,
@@ -69,11 +76,7 @@ def main():
     ni = acq_data.number_of_readouts()
     print('readouts: total %d, image data %d' % (na, ni))
 
-    # sort acquisition data;
-    # currently performed with respect to (in this order):
-    #    - repetition
-    #    - slice
-    #    - kspace encode step 1
+    # sort acquisition data with respect to the acquisition times
     print('sorting...')
     acq_data.sort()
 
@@ -89,11 +92,13 @@ def main():
     if flags0 & mr.IMAGE_DATA_MASK:
         print('first readout is image data')
     else:
-        # should see this if input data file is test_2D_2x.h5
+        # should see this if input data file is e.g. grappa2_6rep.h5
         print('first readout is not image data')
         a0 = acq_data.as_array(0)
         print('first readout shape: %dx%d' % a0.shape)
 
+    print('Checking acquisitions %d to %d...' % (first, last))
+
     # display flags
     print('Flags')
     print(flags)

examples/Python/MR/coil_sensitivity_maps.py

@@ -11,6 +11,7 @@
   -p <path>, --path=<path>    path to data files, defaults to data/examples/MR
                               subfolder of SIRF root folder
   -i <iter>, --iter=<iter>    number of smoothing iterations [default: 10]
+  -c <size>, --conv=<size>    smoothing convolution kernel halfsize [default: 1]
   -e <engn>, --engine=<engn>  reconstruction engine [default: Gadgetron]
   --non-interactive           do not show plots
 '''
@@ -50,6 +51,7 @@
 if data_path is None:
     data_path = examples_data_path('MR')
 nit = int(args['--iter'])
+cks = int(args['--conv'])
 show_plot = not args['--non-interactive']
 
 def main():
@@ -76,7 +78,8 @@ def main():
     CSMs = mr.CoilSensitivityData()
     #
     # set number of smoothing iterations to suppress noise
-    CSMs.smoothness = nit
+    CSMs.smoothing_iterations = nit
+    CSMs.conv_kernel_halfsize = cks
 
     # calculate coil sensitivity maps directly from the raw k-space data by the
     # Square-Root-of-the-Sum-of-Squares over all coils (SRSS) method
@@ -106,6 +109,7 @@ def main():
     # to the image data prior to the calculation of the coil sensitivity maps
     CIs.calculate(processed_data)
     CSs = mr.CoilSensitivityData()
+    CSs.conv_kernel_halfsize = cks
     print('B) calculating from coil images...')
     CSs.calculate(CIs, method='SRSS(niter=%d)' % nit)
     diff = CSs - CSMs

examples/Python/MR/noncartesian_cg_sense.py

@@ -71,7 +71,7 @@ def ConjugateGradient(rawdata, num_iter = 10, stop_criterion = 1e-7):
 
     print('---\n computing coil sensitivity maps...')
     csms = mr.CoilSensitivityData()
-    csms.smoothness = 10
+    csms.smoothing_iterations = 10
     csms.calculate(rawdata)
 
     # create acquisition model based on the acquisition parameters
@@ -137,7 +137,7 @@ def main():
 
     # acquisition data will be read from an HDF file input_file
     # AcquisitionData.set_storage_scheme('memory')
-    acq_data = mr.AcquisitionData(input_file)
+    acq_data = mr.AcquisitionData(input_file, False)
 
     print('---\n acquisition data norm: %e' % acq_data.norm())
 

examples/Python/MR/noncartesian_recon.py

@@ -71,7 +71,7 @@ def main():
 
     # locate the k-space raw data file adn read
     input_file = existing_filepath(data_path, data_file)
-    acq_data = mr.AcquisitionData(input_file)
+    acq_data = mr.AcquisitionData(input_file, False)
 
     # pre-process acquisition data
     if trajtype != 'radial' and trajtype != 'goldenangle':
@@ -109,7 +109,8 @@ def main():
 
         print('---\n computing coil sensitivity maps...')
         csms = mr.CoilSensitivityData()
-        csms.smoothness = 10
+        csms.smoothing_iterations = 10
+        csms.conv_kernel_halfsize = 3
         csms.calculate(processed_data)
 
         if show_plot:

examples/Python/PETMR/generate_MCIR_data.py

@@ -247,7 +247,7 @@ def main():
     # Create coil sensitivity data
     print("Calculating coil sensitivity map...")
     csm = mr.CoilSensitivityData()
-    csm.smoothness = 500
+    csm.smoothing_iterations = 500
     csm.calculate(template_MR_raw)
 
     # Create interleaved sampling

src/Synergistic/tests/test_pSynergistic.py

@@ -105,7 +105,8 @@ def try_complex_resample(raw_mr_filename):
     sys.stderr.write('# --------------------------------------------------------------------------------- #\n')
     time.sleep(0.5)
 
-    raw_mr = mr.AcquisitionData(raw_mr_filename)
+    raw_mr = mr.AcquisitionData(raw_mr_filename, False)
+    print('%d acquisitions read...' % raw_mr.number())
 
     recon_gadgets = ['NoiseAdjustGadget',
                      'AsymmetricEchoAdjustROGadget',
@@ -120,6 +121,7 @@ def try_complex_resample(raw_mr_filename):
     recon.process()
 
     ismrmrd_im = recon.get_output()
+    print('%d images reconstructed...' % ismrmrd_im.number())
 
     if ismrmrd_im.is_real():
         raise AssertionError("Expected output of reconstruction to be complex")

src/common/include/sirf/common/ImageData.h

@@ -71,6 +71,18 @@ namespace sirf {
 			for (; dst != end; ++dst, ++src)
 				*dst = *src;
 		}
+		size_t size() const
+		{
+			Dimensions dim = dimensions();
+			if (is_empty())
+				return 0;
+			size_t n = 1;
+			for (std::map<std::string, int>::iterator it = dim.begin(); it != dim.end(); ++it) {
+				n *= it->second;
+			}
+			return n;
+		}
+
         void fill(const ImageData& im)
         {
             Iterator_const& src = im.begin();

src/xGadgetron/cGadgetron/cgadgetron.cpp

@@ -320,9 +320,11 @@ cGT_setCSParameter(void* ptr, const char* par, const void* val)
 	CAST_PTR(DataHandle, h_csms, ptr);
 	CoilSensitivitiesVector& csms =
 		objectFromHandle<CoilSensitivitiesVector>(h_csms);
-	if (sirf::iequals(par, "smoothness"))
+	if (sirf::iequals(par, "smoothing_iterations"))
 		csms.set_csm_smoothness(dataFromHandle<int>(val));
 	//csms.set_csm_smoothness(intDataFromHandle(val)); // causes problems with Matlab
+	else if (sirf::iequals(par, "conv_kernel_size"))
+		csms.set_csm_conv_kernel_size(dataFromHandle<int>(val));
 	else
 		return unknownObject("parameter", par, __FILE__, __LINE__);
 	return new DataHandle;
@@ -561,14 +563,26 @@ cGT_sortAcquisitionsByTime(void* ptr_acqs)
 
 extern "C"
 void*
-cGT_ISMRMRDAcquisitionsFromFile(const char* file, int all)
+cGT_ISMRMRDAcquisitionsFromFile(const char* file, int all, size_t ptr)
 {
 	if (!file_exists(file))
 		return fileNotFound(file, __FILE__, __LINE__);
 	try {
+		auto ptr_ignored = reinterpret_cast<unsigned long int*>(ptr);
+		auto ignored = *ptr_ignored;
+		std::cout << "reading from " << file << " using ignore mask ";
+		IgnoreMask mask;
+		if (all)
+			mask.set(0);
+		else
+			mask.set(ignored);
+		std::cout << mask.bits_string() << '\n';
 		shared_ptr<MRAcquisitionData>
 			acquisitions(new AcquisitionsVector);
+		IgnoreMask copy = acquisitions->ignore_mask();
+		acquisitions->set_ignore_mask(mask);
 		acquisitions->read(file, all);
+		acquisitions->set_ignore_mask(copy);
 		return newObjectHandle<MRAcquisitionData>(acquisitions);
 	}
 	CATCH;
@@ -703,6 +717,40 @@ cGT_getAcquisitionDataDimensions(void* ptr_acqs, size_t ptr_dim)
 	CATCH;
 }
 
+extern "C"
+void*
+cGT_setAcquisitionsIgnoreMask(void* ptr_acqs, size_t ptr_im)
+{
+	try {
+		CAST_PTR(DataHandle, h_acqs, ptr_acqs);
+		MRAcquisitionData& acqs =
+			objectFromHandle<MRAcquisitionData>(h_acqs);
+		shared_ptr<ISMRMRD::Acquisition>
+			sptr_acq(new ISMRMRD::Acquisition);
+		auto im = *reinterpret_cast<unsigned long int*>(ptr_im);
+		acqs.set_ignore_mask(im);
+		return new DataHandle;
+	}
+	CATCH;
+}
+
+extern "C"
+void*
+cGT_acquisitionsIgnoreMask(void* ptr_acqs, size_t ptr_im)
+{
+	try {
+		CAST_PTR(DataHandle, h_acqs, ptr_acqs);
+		MRAcquisitionData& acqs =
+			objectFromHandle<MRAcquisitionData>(h_acqs);
+		shared_ptr<ISMRMRD::Acquisition>
+			sptr_acq(new ISMRMRD::Acquisition);
+		auto im = reinterpret_cast<unsigned long int*>(ptr_im);
+		*im = acqs.ignore_mask().bits();
+		return new DataHandle;
+	}
+	CATCH;
+}
+
 extern "C"
 void*
 cGT_acquisitionDataAsArray(void* ptr_acqs, size_t ptr_z, int all)