diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index c896b7b8a..58c627dc6 100644 --- a/.github/workflows/build-test.yml +++ b/.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 diff --git a/CHANGES.md b/CHANGES.md index e662b6aa3..f65d44892 100644 --- a/CHANGES.md +++ b/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 diff --git a/examples/Python/MR/acquisition_data.py b/examples/Python/MR/acquisition_data.py index 3f8457570..31077c43f 100644 --- a/examples/Python/MR/acquisition_data.py +++ b/examples/Python/MR/acquisition_data.py @@ -11,6 +11,8 @@ -p , --path= path to data files, defaults to data/examples/MR subfolder of SIRF root folder -s , --slices= max number of slices to display [default: 8] + -i , --ignore= ignore acquisitions with non-zero value of + (try e.g. --ignore=19 and --file=grappa2_6rep.h5) -e , --engine= 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) diff --git a/examples/Python/MR/coil_sensitivity_maps.py b/examples/Python/MR/coil_sensitivity_maps.py index 5b0c4fe4c..96459fc9d 100644 --- a/examples/Python/MR/coil_sensitivity_maps.py +++ b/examples/Python/MR/coil_sensitivity_maps.py @@ -11,6 +11,7 @@ -p , --path= path to data files, defaults to data/examples/MR subfolder of SIRF root folder -i , --iter= number of smoothing iterations [default: 10] + -c , --conv= smoothing convolution kernel halfsize [default: 1] -e , --engine= 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 diff --git a/examples/Python/MR/noncartesian_cg_sense.py b/examples/Python/MR/noncartesian_cg_sense.py index e4b5e86e9..764c2a1ab 100644 --- a/examples/Python/MR/noncartesian_cg_sense.py +++ b/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()) diff --git a/examples/Python/MR/noncartesian_recon.py b/examples/Python/MR/noncartesian_recon.py index 2d18031c3..3bfeee77a 100644 --- a/examples/Python/MR/noncartesian_recon.py +++ b/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: diff --git a/examples/Python/PETMR/generate_MCIR_data.py b/examples/Python/PETMR/generate_MCIR_data.py index d8db5f1fc..eb46d77c0 100644 --- a/examples/Python/PETMR/generate_MCIR_data.py +++ b/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 diff --git a/src/Synergistic/tests/test_pSynergistic.py b/src/Synergistic/tests/test_pSynergistic.py index 5322b035e..8164ee84b 100644 --- a/src/Synergistic/tests/test_pSynergistic.py +++ b/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") diff --git a/src/common/include/sirf/common/ImageData.h b/src/common/include/sirf/common/ImageData.h index f7bc11d14..efd0f1d10 100644 --- a/src/common/include/sirf/common/ImageData.h +++ b/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::iterator it = dim.begin(); it != dim.end(); ++it) { + n *= it->second; + } + return n; + } + void fill(const ImageData& im) { Iterator_const& src = im.begin(); diff --git a/src/xGadgetron/cGadgetron/cgadgetron.cpp b/src/xGadgetron/cGadgetron/cgadgetron.cpp index f6990b351..35d332c54 100644 --- a/src/xGadgetron/cGadgetron/cgadgetron.cpp +++ b/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(h_csms); - if (sirf::iequals(par, "smoothness")) + if (sirf::iequals(par, "smoothing_iterations")) csms.set_csm_smoothness(dataFromHandle(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(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(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 acquisitions(new AcquisitionsVector); + IgnoreMask copy = acquisitions->ignore_mask(); + acquisitions->set_ignore_mask(mask); acquisitions->read(file, all); + acquisitions->set_ignore_mask(copy); return newObjectHandle(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(h_acqs); + shared_ptr + sptr_acq(new ISMRMRD::Acquisition); + auto im = *reinterpret_cast(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(h_acqs); + shared_ptr + sptr_acq(new ISMRMRD::Acquisition); + auto im = reinterpret_cast(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) diff --git a/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp b/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp index aca8f722a..6b66f7eca 100644 --- a/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp +++ b/src/xGadgetron/cGadgetron/gadgetron_data_containers.cpp @@ -98,8 +98,7 @@ MRAcquisitionData::write(const std::string &filename) const void MRAcquisitionData::read(const std::string& filename_ismrmrd_with_ext, int all) { - - bool const verbose = true; + bool const verbose = true; if( verbose ) std::cout<< "Started reading acquisitions from " << filename_ismrmrd_with_ext << std::endl; @@ -109,10 +108,8 @@ MRAcquisitionData::read(const std::string& filename_ismrmrd_with_ext, int all) mtx.lock(); ISMRMRD::Dataset d(filename_ismrmrd_with_ext.c_str(),"dataset", false); d.readHeader(this->acqs_info_); - - ISMRMRD::IsmrmrdHeader hdr = acqs_info_.get_IsmrmrdHeader(); - - uint32_t num_acquis = d.getNumberOfAcquisitions(); + ISMRMRD::IsmrmrdHeader hdr = acqs_info_.get_IsmrmrdHeader(); + uint32_t num_acquis = d.getNumberOfAcquisitions(); mtx.unlock(); std::stringstream str; @@ -144,7 +141,8 @@ MRAcquisitionData::read(const std::string& filename_ismrmrd_with_ext, int all) if( verbose ) { if( i_acqu%( num_acquis/10 ) == 0 ) - std::cout << std::ceil(float(i_acqu) / num_acquis * 100) << "%.." << std::flush; + std::cout << std::ceil(float(i_acqu) / num_acquis * 100) + << "%.." << std::flush; } ISMRMRD::Acquisition acq; @@ -152,10 +150,11 @@ MRAcquisitionData::read(const std::string& filename_ismrmrd_with_ext, int all) d.readAcquisition( i_acqu, acq); mtx.unlock(); - if(all || !TO_BE_IGNORED(acq)) - this->append_acquisition( acq ); + IgnoreMask ignore_mask = this->ignore_mask(); + if (all || !ignore_mask.ignored(acq.flags())) + this->append_acquisition(acq); } - this->sort_by_time(); + this->sort_by_time(); if( verbose ) std::cout<< "\nFinished reading acquisitions from " << filename_ismrmrd_with_ext << std::endl; } @@ -202,6 +201,7 @@ MRAcquisitionData::get_acquisitions_dimensions(size_t ptr_dim) const } num_acq++; } + ASSERT(num_acq > 0, "All acquisitions ignored, dimensions undefined,"); int const num_dims = 3; dim[0] = ns; @@ -217,16 +217,19 @@ uint16_t MRAcquisitionData::get_trajectory_dimensions(void) const ASSERT(na > 0, "You are asking for dimensions on an empty acquisition container. Please don't..."); ISMRMRD::Acquisition acq; - uint16_t traj_dims = 65535; + uint16_t traj_dims; + int num_acq = 0; for (int i = 0; i < na; ++i) { if (!get_acquisition(i, acq)) continue; - if (traj_dims == 65535) + if (num_acq == 0) traj_dims = acq.trajectory_dimensions(); else if (acq.trajectory_dimensions() != traj_dims) throw LocalisedException("Not every acquisition in your container has the same trajectory dimension." , __FILE__, __LINE__); + num_acq++; } + ASSERT(num_acq > 0, "All acquisitions ignored, trajectory dimensions undefined."); return traj_dims; } @@ -236,13 +239,14 @@ void MRAcquisitionData::get_kspace_dimensions(std::vector& dims) const ASSERT(na > 0, "You are asking for dimensions on an empty acquisition container. Please don't..."); ISMRMRD::Acquisition acq; - int nro = -1; + int nro; int nc; + int num_acq = 0; for (int i = 0; i < na; ++i) { if (!get_acquisition(i, acq)) continue; - if (nro == -1) { + if (num_acq == 0) { nro = acq.number_of_samples(); nc = acq.active_channels(); } @@ -252,7 +256,9 @@ void MRAcquisitionData::get_kspace_dimensions(std::vector& dims) const if (acq.number_of_samples() != nro) throw std::runtime_error("The number of readout points is not consistent within this container."); } + num_acq++; } + ASSERT(num_acq > 0, "All acquisitions ignored, some k-space dimensions undefined"); ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader(); ISMRMRD::Encoding e = hdr.encoding[0]; @@ -1289,7 +1295,7 @@ AcquisitionsVector* AcquisitionsVector::clone_impl() const { AcquisitionsVector* ptr_ad = - new AcquisitionsVector(this->acqs_info_); + new AcquisitionsVector(this->acqs_info_, this->ignore_mask_); for (int i = 0; i < number(); i++) { ISMRMRD::Acquisition acq; get_acquisition(i, acq); @@ -1330,7 +1336,8 @@ AcquisitionsVector::set_data(const complex_float_t* z, int all) for (int a = 0, i = 0; a < na; a++) { int ia = index(a); ISMRMRD::Acquisition& acq = *acqs_[ia]; - if (!all && TO_BE_IGNORED(acq)) { + IgnoreMask ignore_mask = this->ignore_mask(); + if (!all && ignore_mask.ignored(acq.flags())) { std::cout << "ignoring acquisition " << ia << '\n'; continue; } @@ -2667,7 +2674,6 @@ void CoilSensitivitiesVector::combine_images_with_coilmaps(GadgetronImageData& c if( img_dims != csm_dims) throw LocalisedException("The data dimensions of the image don't match the sensitivity maps.", __FILE__, __LINE__); - CFImage* ptr_dst_img = new CFImage(Nx, Ny, Nz, 1); //urgh this is so horrible sirf::ImageWrap iw_dst(ISMRMRD::ISMRMRD_CXFLOAT, ptr_dst_img ); @@ -2697,7 +2703,6 @@ void CoilSensitivitiesVector::combine_images_with_coilmaps(GadgetronImageData& c } } - void CoilSensitivitiesVector::calculate(CoilImagesVector& iv) { @@ -2758,10 +2763,6 @@ void CoilSensitivitiesVector::calculate_csm } } - int* object_mask = new int[nx*ny*nz]; - memset(object_mask, 0, nx*ny*nz * sizeof(int)); - - ISMRMRD::NDArray v(cm0); ISMRMRD::NDArray w(cm0); float* ptr_img = img.getDataPtr(); @@ -2778,17 +2779,9 @@ void CoilSensitivitiesVector::calculate_csm } } - float max_im = max_(nx, ny, nz, ptr_img); - float small_grad = max_im * 2 / (nx + ny + 0.0f); - for (int i = 0; i < 3; i++) - smoothen_(nx, ny, nz, nc, v.getDataPtr(), w.getDataPtr(), 0, 1); - float noise = max_diff_(nx, ny, nz, nc, small_grad, - v.getDataPtr(), cm0.getDataPtr()); - mask_noise_(nx, ny, nz, ptr_img, noise, object_mask); - for (int i = 0; i < csm_smoothness_; i++) - smoothen_(nx, ny, nz, nc, cm0.getDataPtr(), w.getDataPtr(), //0, 1); - object_mask, 1); + smoothen_(nx, ny, nz, nc, cm0.getDataPtr(), w.getDataPtr(), + csm_conv_kernel_halfsize_); for (unsigned int z = 0; z < nz; z++) { for (unsigned int y = 0; y < ny; y++) { @@ -2808,7 +2801,7 @@ void CoilSensitivitiesVector::calculate_csm for (unsigned int x = 0; x < nx; x++, i++) { float r = img(x, y, z); float s; - if (r != 0.0 && object_mask[i]) + if (r != 0.0) s = (float)(1.0 / r); else s = 0.0; @@ -2820,28 +2813,12 @@ void CoilSensitivitiesVector::calculate_csm } } - delete[] object_mask; -} - - - -void CoilSensitivitiesVector::mask_noise_ -(int nx, int ny, int nz, float* u, float noise, int* mask) -{ - int i = 0; - for (int iz = 0; iz < nz; iz++) - for (int iy = 0; iy < ny; iy++) - for (int ix = 0; ix < nx; ix++, i++) { - float t = fabs(u[i]); - mask[i] = (t > noise); - } } void CoilSensitivitiesVector::smoothen_ (int nx, int ny, int nz, int nc, - complex_float_t* u, complex_float_t* v, - int* obj_mask, int w) + complex_float_t* u, complex_float_t* v, int w) { const complex_float_t ONE(1.0, 0.0); const complex_float_t TWO(2.0, 0.0); @@ -2849,10 +2826,6 @@ CoilSensitivitiesVector::smoothen_ for (int iz = 0, k = 0; iz < nz; iz++) for (int iy = 0; iy < ny; iy++) for (int ix = 0; ix < nx; ix++, i++, k++) { - if (obj_mask && !obj_mask[k]) { - v[i] = u[i]; - continue; - } int n = 0; complex_float_t r(0.0, 0.0); complex_float_t s(0.0, 0.0); @@ -2864,7 +2837,7 @@ CoilSensitivitiesVector::smoothen_ continue; int j = i + jx + jy*nx; int l = k + jx + jy*nx; - if (i != j && (!obj_mask || obj_mask[l])) { + if (i != j) { n++; r += ONE; s += u[j]; @@ -2878,44 +2851,3 @@ CoilSensitivitiesVector::smoothen_ memcpy(u, v, nx*ny*nz*nc * sizeof(complex_float_t)); } -float -CoilSensitivitiesVector::max_(int nx, int ny, int nz, float* u) -{ - float r = 0.0; - int i = 0; - for (int iz = 0; iz < nz; iz++) - for (int iy = 0; iy < ny; iy++) - for (int ix = 0; ix < nx; ix++, i++) { - float t = fabs(u[i]); - if (t > r) - r = t; - } - return r; -} - -float -CoilSensitivitiesVector::max_diff_ -(int nx, int ny, int nz, int nc, float small_grad, - complex_float_t* u, complex_float_t* v) -{ - int nxy = nx*ny; - int nxyz = nxy*nz; - float s = 0.0f; - for (int ic = 0; ic < nc; ic++) { - for (int iz = 0; iz < nz; iz++) { - for (int iy = 1; iy < ny - 1; iy++) { - for (int ix = 1; ix < nx - 1; ix++) { - int i = ix + nx*iy + nxy*iz + nxyz*ic; - float gx = std::abs(u[i + 1] - u[i - 1]) / 2.0f; - float gy = std::abs(u[i + nx] - u[i - nx]) / 2.0f; - float g = (float)std::sqrt(gx*gx + gy*gy); - float si = std::abs(u[i] - v[i]); - if (g <= small_grad && si > s) - s = si; - } - } - } - } - return s; -} - diff --git a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/cgadgetron.h b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/cgadgetron.h index 396cfbea5..2751b3aed 100644 --- a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/cgadgetron.h +++ b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/cgadgetron.h @@ -64,8 +64,10 @@ extern "C" { void* cGT_AcquisitionModelBackward(void* ptr_am, const void* ptr_acqs); // acquisition data methods - void* cGT_ISMRMRDAcquisitionsFromFile(const char* file, int all); + void* cGT_ISMRMRDAcquisitionsFromFile(const char* file, int all, size_t ptr); void* cGT_ISMRMRDAcquisitionsFile(const char* file); + void* cGT_setAcquisitionsIgnoreMask(void* ptr_acqs, size_t ptr_im); + void* cGT_acquisitionsIgnoreMask(void* ptr_acqs, size_t ptr_im); void* cGT_processAcquisitions(void* ptr_proc, void* ptr_input); void* cGT_acquisitionFromContainer(void* ptr_acqs, unsigned int acq_num); void* cGT_appendAcquisition(void* ptr_acqs, void* ptr_acq); diff --git a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h index 80e30bc66..39443028e 100644 --- a/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h +++ b/src/xGadgetron/cGadgetron/include/sirf/Gadgetron/gadgetron_data_containers.h @@ -53,20 +53,6 @@ limitations under the License. //#define SIRF_DYNAMIC_CAST(T, X, Y) T& X = (T&)Y #define SIRF_DYNAMIC_CAST(T, X, Y) T& X = dynamic_cast(Y) -/*! -\ingroup MR -\brief Acquisitions filter. - -Some acquisitions do not participate directly in the reconstruction process -(e.g. noise calibration acquisitions). -*/ -#define TO_BE_IGNORED(acq) \ - (!(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_MEASUREMENT) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_REVERSE) && \ - (acq).flags() >= (1 << (ISMRMRD::ISMRMRD_ACQ_IS_NOISE_MEASUREMENT - 1))) - /*! \ingroup MR \brief Serialized ISMRMRD acquisition header (cf. ismrmrd.h). @@ -129,6 +115,79 @@ namespace sirf { mutable bool have_header_; }; + /*! + \ingroup Gadgetron Data Containers + \brief Class enabling ignoring acquisitions with certain ISMRMRD acquisition flags. + + Most of the acquisitiondata can be represented by a multidimensional array. + However, the input raw data may also contain acquisitions serving some special + purposes different from the rest, e.g. acquisitions with noise calibration data + (ISMRMRD flag 19). These may have incompatible shape, and must be ignored e.g. + when calculating multidimensional data dimensions or performing algebraic operations + on acquisition data by some reconstruction algorithms. Class IgnoreMask encapsulates + a 64-bit integer with non-zero bits representing 'forbidden' ISMRMRD flags, i.e. any + acquisition acq that has at least one of these flags in its bit field returned by + its acq.flags() method will be ignored. + */ + class IgnoreMask { + public: + IgnoreMask(unsigned long int mask = (1 << 18)) : + ignore_(mask), max_(8*sizeof(mask)) {} + void set(unsigned long int mask) + { + ignore_ = mask; + } + void ignore(unsigned int i) + { + if (i < 1 || i > max_) + return; + unsigned long int one = 1; + ignore_ = ignore_ | (one << (i - 1)); + } + void ignore_not(unsigned int i) + { + if (i < 1 || i > max_) + return; + unsigned long int one = 1; + ignore_ = ignore_ & ~(one << (i - 1)); + } + bool bit(unsigned int i) const + { + if (i < 1 || i > max_) + return true; + unsigned long int one = 1; + return ignore_ & (one << (i - 1)); + } + unsigned long int bits() const + { + return ignore_; + } + bool ignored(unsigned long int bits) const + { + return bits & ignore_; + } + std::string bits_string() const + { + unsigned int size = max_; + unsigned long int one = 1; + unsigned long int bitmask = (one << (size - 1)); + std::stringstream str; + for (unsigned int i = 0; i < size; i++) { + if (ignore_ & (bitmask >> i)) + str << '1'; + else + str << '0'; + if ((i + 1) % 4 == 0) + str << ' '; + } + str << '\n'; + return str.str(); + } + private: + unsigned long int ignore_; + unsigned int max_; + }; + /*! \ingroup Gadgetron Data Containers \brief Class to keep track of order in k-space @@ -434,8 +493,8 @@ namespace sirf { name, name of k-space dimension to modify min_max_ctr = minimium, maximum and center value of sampled region */ - void set_encoding_limits(const std::string& name, \ - const std::tuple min_max_ctr) + void set_encoding_limits(const std::string& name, + const std::tuple min_max_ctr) { ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader(); ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits; @@ -514,8 +573,10 @@ namespace sirf { virtual gadgetron::shared_ptr get_acquisition_sptr(unsigned int num) = 0; - virtual int get_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) const = 0; - virtual void set_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) = 0; + virtual int get_acquisition(unsigned int, + ISMRMRD::Acquisition&) const = 0; + virtual void set_acquisition(unsigned int, + ISMRMRD::Acquisition&) = 0; virtual void append_acquisition(ISMRMRD::Acquisition& acq) = 0; virtual void copy_acquisitions_info(const MRAcquisitionData& ac) = 0; @@ -529,7 +590,7 @@ namespace sirf { virtual void set_data(const complex_float_t* z, int all = 1) = 0; virtual void get_data(complex_float_t* z, int all = 1); - virtual void set_user_floats(float const * const z, int const idx); + virtual void set_user_floats(float const * const z, int const idx); virtual bool is_complex() const { @@ -591,9 +652,11 @@ namespace sirf { AcquisitionsInfo acquisitions_info() const { return acqs_info_; } void set_acquisitions_info(std::string info) { acqs_info_ = info; } - void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;} + void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;} + IgnoreMask ignore_mask() const { return ignore_mask_; } + void set_ignore_mask(IgnoreMask ignore_mask) const { ignore_mask_= ignore_mask; } - ISMRMRD::TrajectoryType get_trajectory_type() const; + ISMRMRD::TrajectoryType get_trajectory_type() const; void set_trajectory_type(const ISMRMRD::TrajectoryType type); @@ -618,8 +681,8 @@ namespace sirf { } bool undersampled() const; - int get_acquisitions_dimensions(size_t ptr_dim) const; - void get_kspace_dimensions(std::vector& dims) const; + int get_acquisitions_dimensions(size_t ptr_dim) const; + void get_kspace_dimensions(std::vector& dims) const; uint16_t get_trajectory_dimensions(void) const; void sort(); @@ -627,31 +690,32 @@ namespace sirf { bool sorted() const { return sorted_; } void set_sorted(bool sorted) { sorted_ = sorted; } - //! Function to get the indices of the acquisitions belonging to different dimensions of k-space - /*! - * All acquisitions belong to only one subset in a multi-dimensional k-space. - * This function returns a vector of sets of indices belonging to the acquisitions of the individual subsets. - */ - std::vector get_kspace_order() const; - - //! Function to get the all KSpaceSubset's of the MRAcquisitionData - std::vector get_kspace_sorting() const { return this->sorting_; } + //! Function to get the indices of the acquisitions belonging to different dimensions of k-space + /*! + * All acquisitions belong to only one subset in a multi-dimensional k-space. + * This function returns a vector of sets of indices belonging to the acquisitions of the individual subsets. + */ + std::vector get_kspace_order() const; - //! Function to go through Acquisitions and assign them to their k-space dimension - /*! - * All acquisitions belong to only one subset in a multi-dimensional k-space. This function goes through - * all acquisitions in the container, extracts their subset (i.e. which slice contrast etc.) and stores - * this information s.t. consisten subsets (i.e. all acquisitions belonging to the same slice) can be - * extracted. - */ - void organise_kspace(); + //! Function to get the all KSpaceSubset's of the MRAcquisitionData + std::vector get_kspace_sorting() const { return this->sorting_; } - virtual std::vector get_flagged_acquisitions_index(const std::vector flags) const; - virtual std::vector get_slice_encoding_index(const unsigned kspace_encode_step_2) const; + //! Function to go through Acquisitions and assign them to their k-space dimension + /*! + * All acquisitions belong to only one subset in a multi-dimensional k-space. This function goes through + * all acquisitions in the container, extracts their subset (i.e. which slice contrast etc.) and stores + * this information s.t. consisten subsets (i.e. all acquisitions belonging to the same slice) can be + * extracted. + */ + void organise_kspace(); + virtual std::vector get_flagged_acquisitions_index + (const std::vector flags) const; + virtual std::vector get_slice_encoding_index + (const unsigned kspace_encode_step_2) const; - virtual void get_subset(MRAcquisitionData& subset, const std::vector subset_idx) const; - virtual void set_subset(const MRAcquisitionData &subset, const std::vector subset_idx); + virtual void get_subset(MRAcquisitionData& subset, const std::vector subset_idx) const; + virtual void set_subset(const MRAcquisitionData &subset, const std::vector subset_idx); std::vector index() { return index_; } const std::vector& index() const { return index_; } @@ -659,7 +723,8 @@ namespace sirf { int index(int i) const { const std::size_t ni = index_.size(); - if (i < 0 || (ni > 0 && static_cast(i) >= ni) || static_cast(i) >= number()) + if (i < 0 || (ni > 0 && static_cast(i) >= ni) + || static_cast(i) >= number()) THROW("Aquisition number is out of range"); if (ni > 0) return index_[i]; @@ -667,23 +732,26 @@ namespace sirf { return i; } - /*! - \brief Reader for ISMRMRD::Acquisition from ISMRMRD file. - * filename_ismrmrd_with_ext: filename of ISMRMRD rawdata file with .h5 extension. - * - * In case the ISMRMRD::Dataset constructor throws an std::runtime_error the reader catches it, - * displays the message and throws it again. - * To avoid reading noise samples and other calibration data, the TO_BE_IGNORED macro is employed - * to exclude potentially incompatible input. - */ + /*! + \brief Reader for ISMRMRD::Acquisition from ISMRMRD file. + * In case the ISMRMRD::Dataset constructor throws an std::runtime_error the reader catches it, + * displays the message and throws it again. + * To avoid reading noise samples and other calibration data, IgnoreMask may be set + * (see method set_ignore_mask) to exclude potentially incompatible input. + + * filename_ismrmrd_with_ext: filename of ISMRMRD rawdata file with .h5 extension. + * all: overrider of the ignore mask - non-zero value forces reading all acquisitions. + */ void read(const std::string& filename_ismrmrd_with_ext, int all = 0); protected: bool sorted_ = false; std::vector index_; - std::vector sorting_; + std::vector sorting_; AcquisitionsInfo acqs_info_; + mutable IgnoreMask ignore_mask_; //= IgnoreMask(); + // new MRAcquisitionData objects will be created from this template // using same_acquisitions_container() static gadgetron::shared_ptr acqs_templ_; @@ -704,13 +772,15 @@ namespace sirf { */ class AcquisitionsVector : public MRAcquisitionData { public: - AcquisitionsVector(const std::string& filename_with_ext, int all = 0) - { - this->read(filename_with_ext, all); - } + AcquisitionsVector(const std::string& filename_with_ext, int all = 0, IgnoreMask ignore_mask = IgnoreMask()) + { + this->set_ignore_mask(ignore_mask); + this->read(filename_with_ext, all); + } - AcquisitionsVector(const AcquisitionsInfo& info = AcquisitionsInfo()) + AcquisitionsVector(const AcquisitionsInfo& info = AcquisitionsInfo(), IgnoreMask ignore_mask = IgnoreMask()) { + this->set_ignore_mask(ignore_mask); acqs_info_ = info; } virtual void empty(); @@ -728,15 +798,13 @@ namespace sirf { int ind = index(num); return acqs_[ind]; } - virtual int get_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) const + virtual int get_acquisition(unsigned int num, + ISMRMRD::Acquisition& acq) const { + IgnoreMask ignore_mask = this->ignore_mask(); int ind = index(num); acq = *acqs_[ind]; - if (!(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_MEASUREMENT) && \ - !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_REVERSE) && \ - (acq).flags() >= (1 << (ISMRMRD::ISMRMRD_ACQ_IS_NOISE_MEASUREMENT - 1))) + if (ignore_mask.ignored(acq.flags())) return 0; return 1; } @@ -755,11 +823,11 @@ namespace sirf { virtual AcquisitionsVector* same_acquisitions_container (const AcquisitionsInfo& info) const { - return new AcquisitionsVector(info); + return new AcquisitionsVector(info, ignore_mask_); } virtual ObjectHandle* new_data_container_handle() const { - DataContainer* ptr = new AcquisitionsVector(acqs_info_); + DataContainer* ptr = new AcquisitionsVector(acqs_info_, ignore_mask_); return new ObjectHandle (gadgetron::shared_ptr(ptr)); } @@ -767,7 +835,7 @@ namespace sirf { new_acquisitions_container() { return gadgetron::unique_ptr - (new AcquisitionsVector(acqs_info_)); + (new AcquisitionsVector(acqs_info_, ignore_mask_)); } private: @@ -1449,7 +1517,14 @@ namespace sirf { throw std::runtime_error("This has not been implemented yet."); } - void set_csm_smoothness(int s){csm_smoothness_ = s;} + void set_csm_smoothness(int s) + { + csm_smoothness_ = s; + } + void set_csm_conv_kernel_size(int w) + { + csm_conv_kernel_halfsize_ = w; + } void calculate(CoilImagesVector& iv); void calculate(const MRAcquisitionData& acq) @@ -1480,7 +1555,10 @@ namespace sirf { private: int csm_smoothness_ = 0; - void smoothen_(int nx, int ny, int nz, int nc, complex_float_t* u, complex_float_t* v, int* obj_mask, int w); + int csm_conv_kernel_halfsize_ = 1; + void smoothen_(int nx, int ny, int nz, int nc, complex_float_t* u, complex_float_t* v, + //int* obj_mask, + int w); void mask_noise_(int nx, int ny, int nz, float* u, float noise, int* mask); float max_diff_(int nx, int ny, int nz, int nc, float small_grad, complex_float_t* u, complex_float_t* v); float max_(int nx, int ny, int nz, float* u); diff --git a/src/xGadgetron/pGadgetron/Gadgetron.py b/src/xGadgetron/pGadgetron/Gadgetron.py index e93329647..111c52f18 100644 --- a/src/xGadgetron/pGadgetron/Gadgetron.py +++ b/src/xGadgetron/pGadgetron/Gadgetron.py @@ -58,7 +58,7 @@ MAX_ACQ_DIMENSIONS = 16 # mask for image-related acquisitions -IMAGE_DATA_MASK = 0x11BFFFF +IMAGE_DATA_MASK = 0x13BFFFF # image type ISMRMRD_IMTYPE_MAGNITUDE = 1 @@ -66,7 +66,7 @@ ISMRMRD_IMTYPE_REAL = 3 ISMRMRD_IMTYPE_IMAG = 4 -#image data type +# image data type ISMRMRD_USHORT = 1 ##, /**< corresponds to uint16_t */ ISMRMRD_SHORT = 2 ##, /**< corresponds to int16_t */ ISMRMRD_UINT = 3 ##, /**< corresponds to uint32_t */ @@ -76,6 +76,68 @@ ISMRMRD_CXFLOAT = 7 ##, /**< corresponds to complex float */ ISMRMRD_CXDOUBLE = 8 ## /**< corresponds to complex double */ +# ISMRMRD acquisition flags +ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP1 = 1 +ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1 = 2 +ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP2 = 3 +ISMRMRD_ACQ_LAST_IN_ENCODE_STEP2 = 4 +ISMRMRD_ACQ_FIRST_IN_AVERAGE = 5 +ISMRMRD_ACQ_LAST_IN_AVERAGE = 6 +ISMRMRD_ACQ_FIRST_IN_SLICE = 7 +ISMRMRD_ACQ_LAST_IN_SLICE = 8 +ISMRMRD_ACQ_FIRST_IN_CONTRAST = 9 +ISMRMRD_ACQ_LAST_IN_CONTRAST = 10 +ISMRMRD_ACQ_FIRST_IN_PHASE = 11 +ISMRMRD_ACQ_LAST_IN_PHASE = 12 +ISMRMRD_ACQ_FIRST_IN_REPETITION = 13 +ISMRMRD_ACQ_LAST_IN_REPETITION = 14 +ISMRMRD_ACQ_FIRST_IN_SET = 15 +ISMRMRD_ACQ_LAST_IN_SET = 16 +ISMRMRD_ACQ_FIRST_IN_SEGMENT = 17 +ISMRMRD_ACQ_LAST_IN_SEGMENT = 18 +ISMRMRD_ACQ_IS_NOISE_MEASUREMENT = 19 +ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION = 20 +ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING = 21 +ISMRMRD_ACQ_IS_REVERSE = 22 +ISMRMRD_ACQ_IS_NAVIGATION_DATA = 23 +ISMRMRD_ACQ_IS_PHASECORR_DATA = 24 +ISMRMRD_ACQ_LAST_IN_MEASUREMENT = 25 +ISMRMRD_ACQ_IS_HPFEEDBACK_DATA = 26 +ISMRMRD_ACQ_IS_DUMMYSCAN_DATA = 27 +ISMRMRD_ACQ_IS_RTFEEDBACK_DATA = 28 +ISMRMRD_ACQ_IS_SURFACECOILCORRECTIONSCAN_DATA = 29 + +ISMRMRD_ACQ_COMPRESSION1 = 53 +ISMRMRD_ACQ_COMPRESSION2 = 54 +ISMRMRD_ACQ_COMPRESSION3 = 55 +ISMRMRD_ACQ_COMPRESSION4 = 56 +ISMRMRD_ACQ_USER1 = 57 +ISMRMRD_ACQ_USER2 = 58 +ISMRMRD_ACQ_USER3 = 59 +ISMRMRD_ACQ_USER4 = 60 +ISMRMRD_ACQ_USER5 = 61 +ISMRMRD_ACQ_USER6 = 62 +ISMRMRD_ACQ_USER7 = 63 +ISMRMRD_ACQ_USER8 = 64 + + +# acquisition ignoring helper class +class IgnoreMask(object): + def __init__(self, mask=None): + if mask is None: + self.mask = 1 << 18 + else: + self.mask = mask + def set(self, mask): + self.mask = mask + def ignore(self, bit): + self.mask = self.mask | 1 << (bit - 1) + def ignore_not(self, bit): + self.mask = self.mask & ~(1 << (bit - 1)) + def ignored(self, bits): + return bits & self.mask + + # data path finding helper functions @deprecated( deprecated_in="2.0.0", removed_in="4.0", current_version=sirf.__version__, @@ -547,7 +609,8 @@ class CoilSensitivityData(ImageData): ''' def __init__(self): self.handle = None - self.smoothness = 0 + self.smoothing_iterations = 0 + self.conv_kernel_halfsize = 1 def __del__(self): if self.handle is not None: pyiutil.deleteDataHandle(self.handle) @@ -573,7 +636,8 @@ def calculate(self, data, method=None): pyiutil.deleteDataHandle(self.handle) self.handle = pygadgetron.cGT_CoilSensitivities('') check_status(self.handle) - nit = self.smoothness + nit = self.smoothing_iterations + w = self.conv_kernel_halfsize # convolution kernel size is (2w+1)-by-(2w+1) pixels if method is not None: method_name, parm_list = name_and_parameters(method) @@ -584,7 +648,8 @@ def calculate(self, data, method=None): method_name = 'SRSS' parm = {} - parms.set_int_par(self.handle, 'coil_sensitivity', 'smoothness', nit) + parms.set_int_par(self.handle, 'coil_sensitivity', 'smoothing_iterations', nit) + parms.set_int_par(self.handle, 'coil_sensitivity', 'conv_kernel_size', w) if isinstance(data, AcquisitionData): self.__calc_from_acquisitions(data, method_name) @@ -823,12 +888,14 @@ class AcquisitionData(DataContainer): Class for an MR acquisitions container. Each item is a 2D complex array of acquisition samples for each coil. ''' - def __init__(self, file=None, all_=False): + def __init__(self, file=None, all_=False, ignored=IgnoreMask()): self.handle = None self.sorted = False self.info = None if file is not None: - self.handle = pygadgetron.cGT_ISMRMRDAcquisitionsFromFile(file, 1*all_) + mask = numpy.ndarray((1,), dtype=numpy.int64) + mask[0] = ignored.mask + self.handle = pygadgetron.cGT_ISMRMRDAcquisitionsFromFile(file, 1*all_, mask.ctypes.data) check_status(self.handle) def __del__(self): @@ -856,6 +923,14 @@ def new_acquisition_data(self, empty=True): new_ad.handle = pygadgetron.cGT_cloneAcquisitions(self.handle) check_status(new_ad.handle) return new_ad + def set_ignore_mask(self, ignored): + mask = numpy.ndarray((1,), dtype=numpy.int64) + mask[0] = ignored.mask + try_calling(pygadgetron.cGT_setAcquisitionsIgnoreMask(self.handle, mask.ctypes.data)) + def ignore_mask(self): + mask = numpy.ndarray((1,), dtype=numpy.int64) + try_calling(pygadgetron.cGT_acquisitionsIgnoreMask(self.handle, mask.ctypes.data)) + return mask[0] def number_of_readouts(self, select='image'): if select == 'image': dim = self.dimensions()