SIRF  3.6.0
gadgetron_image_wrap.h
Go to the documentation of this file.
1 /*
2 SyneRBI Synergistic Image Reconstruction Framework (SIRF)
3 Copyright 2015 - 2020 Rutherford Appleton Laboratory STFC
4 Copyright 2020, 2021 University College London
5 
6 This is software developed for the Collaborative Computational
7 Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
8 (http://www.ccpsynerbi.ac.uk/).
9 
10 Licensed under the Apache License, Version 2.0 (the "License");
11 you may not use this file except in compliance with the License.
12 You may obtain a copy of the License at
13 http://www.apache.org/licenses/LICENSE-2.0
14 Unless required by applicable law or agreed to in writing, software
15 distributed under the License is distributed on an "AS IS" BASIS,
16 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 See the License for the specific language governing permissions and
18 limitations under the License.
19 
20 */
21 
31 #ifndef GADGETRON_IMAGE_WRAP_TYPE
32 #define GADGETRON_IMAGE_WRAP_TYPE
33 
34 #include <ismrmrd/ismrmrd.h>
35 #include <ismrmrd/dataset.h>
36 #include <ismrmrd/meta.h>
37 #include <ismrmrd/xml.h>
38 
39 #include "sirf/common/ANumRef.h"
40 #include "sirf/common/iequals.h"
41 #include "sirf/Gadgetron/cgadgetron_shared_ptr.h"
42 #include "sirf/Gadgetron/xgadgetron_utilities.h"
43 
44 #define IMAGE_PROCESSING_SWITCH(Type, Operation, Arguments, ...)\
45  if (Type == ISMRMRD::ISMRMRD_USHORT)\
46  Operation ((ISMRMRD::Image<unsigned short>*) Arguments, ##__VA_ARGS__);\
47  else if (Type == ISMRMRD::ISMRMRD_SHORT)\
48  Operation ((ISMRMRD::Image<short>*) Arguments, ##__VA_ARGS__);\
49  else if (Type == ISMRMRD::ISMRMRD_UINT)\
50  Operation ((ISMRMRD::Image<unsigned int>*) Arguments, ##__VA_ARGS__);\
51  else if (Type == ISMRMRD::ISMRMRD_INT)\
52  Operation ((ISMRMRD::Image<int>*) Arguments, ##__VA_ARGS__);\
53  else if (Type == ISMRMRD::ISMRMRD_FLOAT)\
54  Operation ((ISMRMRD::Image<float>*) Arguments, ##__VA_ARGS__);\
55  else if (Type == ISMRMRD::ISMRMRD_DOUBLE)\
56  Operation ((ISMRMRD::Image<double>*) Arguments, ##__VA_ARGS__);\
57  else if (Type == ISMRMRD::ISMRMRD_CXFLOAT)\
58  Operation ((ISMRMRD::Image< std::complex<float> >*) Arguments, \
59  ##__VA_ARGS__);\
60  else if (Type == ISMRMRD::ISMRMRD_CXDOUBLE)\
61  Operation ((ISMRMRD::Image< std::complex<double> >*) Arguments, \
62  ##__VA_ARGS__);\
63  else\
64  throw std::domain_error("unknown data type in IMAGE_PROCESSING_SWITCH");
65 
66 #define IMAGE_PROCESSING_SWITCH_CONST(Type, Operation, Arguments, ...)\
67  if (Type == ISMRMRD::ISMRMRD_USHORT)\
68  Operation ((const ISMRMRD::Image<unsigned short>*) Arguments, \
69  ##__VA_ARGS__);\
70  else if (Type == ISMRMRD::ISMRMRD_SHORT)\
71  Operation ((const ISMRMRD::Image<short>*) Arguments, ##__VA_ARGS__);\
72  else if (Type == ISMRMRD::ISMRMRD_UINT)\
73  Operation ((const ISMRMRD::Image<unsigned int>*) Arguments, ##__VA_ARGS__);\
74  else if (Type == ISMRMRD::ISMRMRD_INT)\
75  Operation ((const ISMRMRD::Image<int>*) Arguments, ##__VA_ARGS__);\
76  else if (Type == ISMRMRD::ISMRMRD_FLOAT)\
77  Operation ((const ISMRMRD::Image<float>*) Arguments, ##__VA_ARGS__);\
78  else if (Type == ISMRMRD::ISMRMRD_DOUBLE)\
79  Operation ((const ISMRMRD::Image<double>*) Arguments, ##__VA_ARGS__);\
80  else if (Type == ISMRMRD::ISMRMRD_CXFLOAT)\
81  Operation ((const ISMRMRD::Image< std::complex<float> >*) Arguments, \
82  ##__VA_ARGS__);\
83  else if (Type == ISMRMRD::ISMRMRD_CXDOUBLE)\
84  Operation ((const ISMRMRD::Image< std::complex<double> >*) Arguments, \
85  ##__VA_ARGS__);\
86  else\
87  throw std::domain_error("unknown data type in IMAGE_PROCESSING_SWITCH_CONST");
88 
89 
90 typedef ISMRMRD::Image<complex_float_t> CFImage;
91 typedef ISMRMRD::Image<complex_double_t> CDImage;
92 
93 #ifdef _MSC_VER
94 #pragma warning( push )
95 // disable warning about loss of precision. we convert to float in various places
96 #pragma warning( disable : 4244 )
97 #endif
98 
99 namespace sirf {
100 
106  class ImageWrap {
107  public:
108  class Iterator {
109  public:
110  Iterator(int type, void* data, unsigned int dsize, size_t n) :
111  type_(type), ptr_((char*)data), dsize_(dsize), n_(n), i_(0),
112  ref_(data, type)
113  {}
114  Iterator(const Iterator& iter)
115  {
116  type_ = iter.type_;
117  ptr_ = iter.ptr_;
118  dsize_ = iter.dsize_;
119  n_ = iter.n_;
120  i_ = iter.i_;
121  sptr_iter_ = iter.sptr_iter_;
122  ref_.copy(iter.ref_);
123  }
124  bool operator!=(const Iterator& i) const
125  {
126  return ptr_ != i.ptr_;
127  }
128  bool operator==(const Iterator& i) const
129  {
130  return ptr_ == i.ptr_;
131  }
132  Iterator& operator++()
133  {
134  if (i_ >= n_)
135  throw std::out_of_range("cannot advance out-of-range iterator");
136  i_++;
137  ptr_ += dsize_;
138  return *this;
139  }
140  // too inefficient
141  //virtual Iterator& operator++(int)
142  //{
143  // //sptr_iter_.reset(new Iterator(type_, ptr_, dsize_, n_));
144  // sptr_iter_.reset(new Iterator(*this));
145  // if (i_ >= n_)
146  // throw std::out_of_range("cannot advance out-of-range iterator");
147  // i_++;
148  // ptr_ += dsize_;
149  // return *sptr_iter_;
150  //}
151  NumRef& operator*()
152  {
153  if (i_ >= n_)
154  throw std::out_of_range
155  ("cannot dereference out-of-range iterator");
156  ref_.set_ptr(ptr_);
157  return ref_;
158  }
159  Iterator& operator=(const Iterator& iter)
160  {
161  type_ = iter.type_;
162  ptr_ = iter.ptr_;
163  dsize_ = iter.dsize_;
164  n_ = iter.n_;
165  i_ = iter.i_;
166  sptr_iter_ = iter.sptr_iter_;
167  ref_.copy(iter.ref_);
168  return *this;
169  }
170  private:
171  int type_;
172  char* ptr_;
173  unsigned int dsize_;
174  size_t n_;
175  size_t i_;
176  gadgetron::shared_ptr<Iterator> sptr_iter_;
177  NumRef ref_;
178  };
179  class Iterator_const {
180  public:
181  Iterator_const(int type, void* data, unsigned int dsize, size_t n) :
182  type_(type), ptr_((char*)data), dsize_(dsize), n_(n), i_(0),
183  ref_(data, type)
184  {}
185  Iterator_const(const Iterator_const& iter)
186  {
187  type_ = iter.type_;
188  ptr_ = iter.ptr_;
189  dsize_ = iter.dsize_;
190  n_ = iter.n_;
191  i_ = iter.i_;
192  sptr_iter_ = iter.sptr_iter_;
193  ref_.copy(iter.ref_);
194  }
195  bool operator!=(const Iterator_const& i) const
196  {
197  return ptr_ != i.ptr_;
198  }
199  bool operator==(const Iterator_const& i) const
200  {
201  return ptr_ == i.ptr_;
202  }
203  Iterator_const& operator++()
204  {
205  if (i_ >= n_)
206  throw std::out_of_range("cannot advance out-of-range iterator");
207  i_++;
208  ptr_ += dsize_;
209  return *this;
210  }
211  //virtual Iterator_const& operator++(int)
212  //{
213  // //sptr_iter_.reset(new Iterator_const(type_, ptr_, dsize_, n_));
214  // sptr_iter_.reset(new Iterator_const(*this));
215  // if (i_ >= n_)
216  // throw std::out_of_range("cannot advance out-of-range iterator");
217  // i_++;
218  // ptr_ += dsize_;
219  // return *sptr_iter_;
220  //}
221  const NumRef& operator*() const
222  {
223  if (i_ >= n_) {
224  std::cout << i_ << ' ' << n_ << '\n';
225  throw std::out_of_range
226  ("cannot dereference out-of-range iterator");
227  }
228  ref_.set_ptr(ptr_);
229  return ref_;
230  }
231  Iterator_const& operator=(const Iterator_const& iter)
232  {
233  type_ = iter.type_;
234  ptr_ = iter.ptr_;
235  dsize_ = iter.dsize_;
236  n_ = iter.n_;
237  i_ = iter.i_;
238  sptr_iter_ = iter.sptr_iter_;
239  ref_.copy(iter.ref_);
240  return *this;
241  }
242  private:
243  int type_;
244  char* ptr_;
245  unsigned int dsize_;
246  size_t n_;
247  size_t i_;
248  gadgetron::shared_ptr<Iterator_const> sptr_iter_;
249  mutable NumRef ref_;
250  };
251  ImageWrap(uint16_t type, void* ptr_im)
252  {
253  type_ = type;
254  ptr_ = ptr_im;
255  }
256  ImageWrap(uint16_t type, ISMRMRD::Dataset& dataset, const char* var, int index) noexcept(false)
257  {
258  type_ = type;
259  IMAGE_PROCESSING_SWITCH(type_, read_, ptr_, dataset, var, index, &ptr_);
260  }
261  ImageWrap(const ImageWrap& iw) noexcept(false)
262  {
263  type_ = iw.type();
264  IMAGE_PROCESSING_SWITCH(type_, copy_, iw.ptr_image());
265  }
266  ~ImageWrap() noexcept
267  {
268  try {
269  IMAGE_PROCESSING_SWITCH(type_, delete, ptr_);
270  }
271  catch(...)
272  {
273  // catch exceptions in destructor to prevent crashes
274  }
275  }
276  int type() const
277  {
278  return type_;
279  }
280  void* ptr_image()
281  {
282  return ptr_;
283  }
284  const void* ptr_image() const
285  {
286  return ptr_;
287  }
288  Iterator& begin()
289  {
290  size_t n;
291  unsigned int dsize;
292  char* ptr;
293  IMAGE_PROCESSING_SWITCH
294  (type_, get_data_parameters_, ptr_, &n, &dsize, &ptr);
295  begin_.reset(new Iterator(type_, ptr, dsize, n));
296  return *begin_;
297  }
298  Iterator_const& begin_const() const
299  {
300  size_t n;
301  unsigned int dsize;
302  char* ptr;
303  IMAGE_PROCESSING_SWITCH_CONST
304  (type_, get_data_parameters_, ptr_, &n, &dsize, &ptr);
305  begin_const_.reset(new Iterator_const(type_, ptr, dsize, n));
306  return *begin_const_;
307  }
308  Iterator& end()
309  {
310  size_t n;
311  unsigned int dsize;
312  char* ptr;
313  IMAGE_PROCESSING_SWITCH
314  (type_, get_data_parameters_, ptr_, &n, &dsize, &ptr);
315  end_.reset(new Iterator(type_, ptr + n*dsize, dsize, n));
316  return *end_;
317  }
318  Iterator_const& end_const() const
319  {
320  size_t n;
321  unsigned int dsize;
322  char* ptr;
323  IMAGE_PROCESSING_SWITCH_CONST
324  (type_, get_data_parameters_, ptr_, &n, &dsize, &ptr);
325  end_const_.reset(new Iterator_const(type_, ptr + n*dsize, dsize, n));
326  //std::cout << type_ << ' ' << n << ' ' << dsize << '\n';
327  return *end_const_;
328  }
329  size_t size() const
330  {
331  size_t s;
332  IMAGE_PROCESSING_SWITCH_CONST(type_, get_size_, ptr_, s);
333  return s;
334  }
335  size_t num_data_elm() const
336  {
337  IMAGE_PROCESSING_SWITCH_CONST(type_, return num_data_elm_, ptr_);
338  }
339  ISMRMRD::ImageHeader& head()
340  {
341  IMAGE_PROCESSING_SWITCH(type_, return get_head_ref_, ptr_);
342  }
343  const ISMRMRD::ImageHeader& head() const
344  {
345  IMAGE_PROCESSING_SWITCH_CONST(type_, return get_head_ref_const_, ptr_);
346  }
347  std::string attributes() const
348  {
349  std::string attr;
350  IMAGE_PROCESSING_SWITCH_CONST(type_, get_attr_, ptr_, attr);
351  return attr;
352  }
353  void set_attributes(const std::string& attr)
354  {
355  IMAGE_PROCESSING_SWITCH(type_, set_attr_, ptr_, attr);
356  }
357  void show_attributes() const
358  {
359  IMAGE_PROCESSING_SWITCH_CONST(type_, show_attr_, ptr_, 0);
360  }
361  void set_imtype(ISMRMRD::ISMRMRD_ImageTypes imtype)
362  {
363  IMAGE_PROCESSING_SWITCH(type_, set_imtype_, ptr_, imtype);
364  }
365  size_t get_dim(int* dim) const
366  {
367  IMAGE_PROCESSING_SWITCH_CONST(type_, get_dim_, ptr_, dim);
368  size_t n = dim[0];
369  n *= dim[1];
370  n *= dim[2];
371  n *= dim[3];
372  return n;
373  }
374  void get_data(float* data) const
375  {
376  //std::cout << "in get_data\n";
377  //std::cout << "trying new image wrap iterator...\n";
378  ImageWrap::Iterator_const i = begin_const();
379  ImageWrap::Iterator_const stop = end_const();
380  for (; i != stop; ++data, ++i) {
381  *data = *i;
382  }
383  //IMAGE_PROCESSING_SWITCH_CONST(type_, get_data_, ptr_, data);
384  }
385  void set_data(const float* data)
386  {
387  //std::cout << "in set_data\n";
388  for (ImageWrap::Iterator i = begin(); i != end(); ++i, ++data)
389  *i = *data;
390  //IMAGE_PROCESSING_SWITCH(type_, set_data_, ptr_, data);
391  }
392  void fill(float s)
393  {
394  for (ImageWrap::Iterator i = begin(); i != end(); ++i)
395  *i = s;
396  }
397  void scale(float s)
398  {
399  for (ImageWrap::Iterator i = begin(); i != end(); ++i)
400  *i /= s;
401  }
402  void get_complex_data(complex_float_t* data) const
403  {
404  //std::cout << "in get_complex_data\n";
405  //std::cout << "trying new const image wrap iterator...\n";
406  ImageWrap::Iterator_const i = begin_const();
407  ImageWrap::Iterator_const stop = end_const();
408  for (; i != stop; ++data, ++i) {
409  *data = (*i).complex_float();
410  }
411  //IMAGE_PROCESSING_SWITCH_CONST(type_, get_complex_data_, ptr_, data);
412  }
413 
414  void set_complex_data(const complex_float_t* data)
415  {
416  //std::cout << "in set_complex_data\n";
417  //std::cout << "trying new image wrap iterator...\n";
418  ImageWrap::Iterator i = begin();
419  ImageWrap::Iterator stop = end();
420  for (; i != stop; ++i, ++data) {
421  *i = *data;
422  }
423  //IMAGE_PROCESSING_SWITCH(type_, set_complex_data_, ptr_, data);
424  }
425 
426  gadgetron::shared_ptr<ImageWrap> abs() const
427  {
428  return real("abs");
429  }
430 
431  gadgetron::shared_ptr<ImageWrap> real(const std::string& way = "real") const
432  {
433  int dim[4];
434  get_dim(dim);
435  ISMRMRD::Image<float>* ptr_im = new ISMRMRD::Image<float>(dim[0], dim[1], dim[2], dim[3]);
436  ISMRMRD::Image<float>& im = *ptr_im;
437 
438  ISMRMRD::ImageHeader header = head();
439  header.data_type = ISMRMRD::ISMRMRD_FLOAT;
440  header.image_type = ISMRMRD::ISMRMRD_IMTYPE_MAGNITUDE;
441  header.image_series_index = 0;
442  im.setHead(header);
443 
444  im.setAttributeString(attributes());
445 
446  ImageWrap::Iterator_const i = begin_const();
447  ImageWrap::Iterator_const stop = end_const();
448  float* data = im.getDataPtr();
449  if (sirf::iequals(way, "abs"))
450  for (; i != stop; ++data, ++i) {
451  *data = std::abs((*i).complex_float());
452  }
453  else if (sirf::iequals(way, "real"))
454  for (; i != stop; ++data, ++i) {
455  *data = std::real((*i).complex_float());
456  }
457  else
458  THROW("unknown conversion to real specified in ImageWrap::real");
459 
460  ImageWrap* ptr_iw = new ImageWrap(ISMRMRD::ISMRMRD_FLOAT, ptr_im);
461  return gadgetron::shared_ptr<ImageWrap>(ptr_iw);
462  }
463 
465  ISMRMRD::ISMRMRD_DataTypes get_data_type() const
466  {
467  ISMRMRD::ISMRMRD_DataTypes data_type;
468  IMAGE_PROCESSING_SWITCH_CONST(type_, get_data_type_, ptr_, &data_type);
469  return data_type;
470  }
472  bool is_complex() const
473  {
474  ISMRMRD::ISMRMRD_DataTypes data_type = this->get_data_type();
475  if (data_type == ISMRMRD::ISMRMRD_DataTypes::ISMRMRD_CXFLOAT ||
476  data_type == ISMRMRD::ISMRMRD_DataTypes::ISMRMRD_CXDOUBLE) {
477  return true;
478  }
479  else
480  return false;
481  }
482  void write(ISMRMRD::Dataset& dataset) const
483  {
484  IMAGE_PROCESSING_SWITCH_CONST(type_, write_, ptr_, dataset);
485  }
486  void read(ISMRMRD::Dataset& dataset, const char* var, int ind)
487  {
488  IMAGE_PROCESSING_SWITCH(type_, read_, ptr_, dataset, var, ind, &ptr_);
489  }
490  //void axpby(complex_float_t a, const ImageWrap& x, complex_float_t b)
491  //{
492  // IMAGE_PROCESSING_SWITCH(type_, axpby_, x.ptr_image(), a, b);
493  //}
494  void axpby(complex_float_t a, const ImageWrap& x, complex_float_t b,
495  const ImageWrap& y)
496  {
497  xapyb(x, a, y, b);
498  //IMAGE_PROCESSING_SWITCH(type_, axpby_, x.ptr_image(), a, y.ptr_image(), b);
499  }
500  void xapyb(const ImageWrap& x, complex_float_t a,
501  const ImageWrap& y, complex_float_t b)
502  {
503  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), &a,
504  y.ptr_image(), &b, 0, 0);
505  }
506  void xapyb(const ImageWrap& x, complex_float_t a,
507  const ImageWrap& y, const ImageWrap& b)
508  {
509  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), &a,
510  y.ptr_image(), b.ptr_image(), 0, 1);
511  }
512  void xapyb(const ImageWrap& x, const ImageWrap& a,
513  const ImageWrap& y, complex_float_t b)
514  {
515  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), a.ptr_image(),
516  y.ptr_image(), &b, 1, 0);
517  }
518  void xapyb(const ImageWrap& x, const ImageWrap& a,
519  const ImageWrap& y, const ImageWrap& b)
520  {
521  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), a.ptr_image(),
522  y.ptr_image(), b.ptr_image(), 1, 1);
523  }
524  void binary_op(const ImageWrap& x, const ImageWrap& y, complex_float_t(*f)(complex_float_t, complex_float_t))
525  {
526  IMAGE_PROCESSING_SWITCH(type_, binary_op_, x.ptr_image(), y.ptr_image(), f);
527  }
528  void semibinary_op(const ImageWrap& x, complex_float_t y, complex_float_t(*f)(complex_float_t, complex_float_t))
529  {
530  IMAGE_PROCESSING_SWITCH(type_, semibinary_op_, x.ptr_image(), y, f);
531  }
532  void unary_op(const ImageWrap& x, complex_float_t(*f)(complex_float_t))
533  {
534  IMAGE_PROCESSING_SWITCH(type_, unary_op_, x.ptr_image(), f);
535  }
536  complex_float_t dot(const ImageWrap& iw) const
537  {
538  complex_float_t z;
539  IMAGE_PROCESSING_SWITCH_CONST(type_, dot_, iw.ptr_image(), &z);
540  return z;
541  }
542  float norm() const
543  {
544  float r;
545  IMAGE_PROCESSING_SWITCH_CONST(type_, norm_, ptr_, &r);
546  return r;
547  }
548  complex_float_t sum() const
549  {
550  complex_float_t s;
551  IMAGE_PROCESSING_SWITCH_CONST(type_, sum_, ptr_, &s);
552  return s;
553  }
554  complex_float_t max() const
555  {
556  complex_float_t s;
557  IMAGE_PROCESSING_SWITCH_CONST(type_, max_, ptr_, &s);
558  return s;
559  }
560  float diff(ImageWrap& iw) const
561  {
562  float s;
563  IMAGE_PROCESSING_SWITCH_CONST(type_, diff_, iw.ptr_image(), &s);
564  return s;
565  }
566  void conjugate()
567  {
568  IMAGE_PROCESSING_SWITCH(type_, conjugate_, ptr_);
569  }
570 
571  private:
572  int type_;
573  void* ptr_;
574  mutable gadgetron::shared_ptr<Iterator> begin_;
575  mutable gadgetron::shared_ptr<Iterator> end_;
576  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
577  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
578 
579  ImageWrap& operator=(const ImageWrap& iw)
580  {
581  //type_ = iw.type();
582  //IMAGE_PROCESSING_SWITCH(type_, copy_, iw.ptr_image());
583  return *this;
584  }
585 
586  template<typename T>
587  void get_data_parameters_
588  (const ISMRMRD::Image<T>* ptr_im, size_t *n, unsigned int *dsize,
589  char** data) const
590  {
591  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
592  unsigned int dim[4];
593  dim[0] = im.getMatrixSizeX();
594  dim[1] = im.getMatrixSizeY();
595  dim[2] = im.getMatrixSizeZ();
596  dim[3] = im.getNumberOfChannels();
597  *data = (char*)im.getDataPtr();
598  *dsize = sizeof(T);
599  *n = dim[0];
600  *n *= dim[1];
601  *n *= dim[2];
602  *n *= dim[3];
603  //*n = ptr_im->getDataSize()/(*dsize);
604  }
605 
606  template<typename T>
607  size_t num_data_elm_(const ISMRMRD::Image<T>* ptr) const
608  {
609  return ptr->getNumberOfDataElements();
610  }
611 
612  template<typename T>
613  void copy_(const ISMRMRD::Image<T>* ptr_im)
614  {
615  type_ = ptr_im->getDataType();
616  ptr_ = (void*)new ISMRMRD::Image<T>(*ptr_im);
617  }
618 
619  template<typename T>
620  ISMRMRD::ImageHeader& get_head_ref_(ISMRMRD::Image<T>* ptr_im)
621  {
622  return ptr_im->getHead();
623  }
624 
625  template<typename T>
626  const ISMRMRD::ImageHeader& get_head_ref_const_(const ISMRMRD::Image<T>* ptr_im) const
627  {
628  return ptr_im->getHead();
629  }
630 
631  template<typename T>
632  void set_imtype_(ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_ImageTypes type)
633  {
634  ptr_im->setImageType(type);
635  }
636 
637  template<typename T>
638  void get_size_(const ISMRMRD::Image<T>* ptr_im, size_t& size) const
639  {
640  size = ptr_im->getDataSize();
641  }
642 
643  template<typename T>
644  void get_attr_(const ISMRMRD::Image<T>* ptr_im, std::string& attr) const
645  {
646  ptr_im->getAttributeString(attr);
647  }
648 
649  template<typename T>
650  void set_attr_(ISMRMRD::Image<T>* ptr_im, const std::string& attr)
651  {
652  ptr_im->setAttributeString(attr);
653  }
654 
655  template<typename T>
656  void show_attr_(const ISMRMRD::Image<T>* ptr_im, int i) const
657  {
658  const ISMRMRD::Image<T>& im = *ptr_im;
659  size_t meta_attrib_length = im.getAttributeStringLength();
660  std::cout << "meta_attrib_length: " << meta_attrib_length << '\n';
661  std::string meta_attrib(meta_attrib_length + 1, 0);
662  im.getAttributeString(meta_attrib);
663 
664  //std::cout << "attributes:" << std::endl << meta_attrib << std::endl;
665  if (meta_attrib_length > 0) {
666  ISMRMRD::MetaContainer mc;
667  ISMRMRD::deserialize(meta_attrib.c_str(), mc);
668  for (auto it = mc.begin(); it != mc.end(); ++it) {
669  auto name = it->first.c_str();
670  std::cout << name << ":\n";
671  size_t l = mc.length(name);
672  for (int j = 0; j <l; j++)
673  std::cout << mc.as_str(name, j) << '\n';
674  }
675  std::stringstream meta_xml;
676  ISMRMRD::serialize(mc, meta_xml);
677  std::cout << meta_xml.str().c_str() << '\n';
678  }
679  }
680 
681  template<typename T>
682  void write_
683  (const ISMRMRD::Image<T>* ptr_im, ISMRMRD::Dataset& dataset) const
684  {
685  //std::cout << "appending image..." << std::endl;
686  const ISMRMRD::Image<T>& im = *ptr_im;
687  std::stringstream ss;
688  ss << "image_" << im.getHead().image_series_index;
689  std::string image_varname = ss.str();
690  {
691  Mutex mtx;
692  mtx.lock();
693  dataset.appendImage(image_varname, im);
694  mtx.unlock();
695  }
696  }
697 
698  template<typename T>
699  void read_
700  (const ISMRMRD::Image<T>* ptr,
701  ISMRMRD::Dataset& dataset, const char* var, int index,
702  void** ptr_ptr)
703  {
704  ISMRMRD::Image < T >* ptr_im = new ISMRMRD::Image < T > ;
705  *ptr_ptr = (void*)ptr_im;
706  ISMRMRD::Image<T>& im = *ptr_im;
707  dataset.readImage(var, index, im);
708  //int status = ismrmrd_read_image(&dataset, var, (uint32_t)index, &(im.im));
709  }
710 
711  template<typename T>
712  void get_dim_(const ISMRMRD::Image<T>* ptr_im, int* dim) const
713  {
714  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
715  dim[0] = im.getMatrixSizeX();
716  dim[1] = im.getMatrixSizeY();
717  dim[2] = im.getMatrixSizeZ();
718  dim[3] = im.getNumberOfChannels();
719  }
720  template<typename T>
721  void get_data_type_(const ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_DataTypes* data_type_ptr) const
722  {
723  const ISMRMRD::Image<T>& im = *static_cast<const ISMRMRD::Image<T>*>(ptr_im);
724  *data_type_ptr = im.getDataType();
725  }
726 
727  template<typename T>
728  void get_data_(const ISMRMRD::Image<T>* ptr_im, float* data) const
729  {
730  const ISMRMRD::Image<T>& im = *ptr_im;
731  const T* ptr = im.getDataPtr();
732  size_t n = im.getNumberOfDataElements();
733  for (size_t i = 0; i < n; i++)
734  data[i] = std::real(ptr[i]);
735  }
736 
737  template<typename T>
738  void set_data_(ISMRMRD::Image<T>* ptr_im, const float* data)
739  {
740  ISMRMRD::Image<T>& im = *ptr_im;
741  T* ptr = im.getDataPtr();
742  size_t n = im.getNumberOfDataElements();
743  for (size_t i = 0; i < n; i++)
744  ptr[i] = (T)data[i];
745  }
746 
747  template<typename T>
748  void get_complex_data_
749  (const ISMRMRD::Image<T>* ptr_im, complex_float_t* data) const
750  {
751  const ISMRMRD::Image<T>& im = *ptr_im;
752  const T* ptr = im.getDataPtr();
753  size_t n = im.getNumberOfDataElements();
754  for (size_t i = 0; i < n; i++)
755  data[i] = complex_float_t(ptr[i]);
756  }
757 
758  template<typename T>
759  void set_complex_data_
760  (ISMRMRD::Image<T>* ptr_im, const complex_float_t* data)
761  {
762  ISMRMRD::Image<T>& im = *ptr_im;
763  ISMRMRD::ISMRMRD_DataTypes type = im.getDataType();
764  T* ptr = im.getDataPtr();
765  size_t n = im.getNumberOfDataElements();
766  if (type == ISMRMRD::ISMRMRD_CXFLOAT || type == ISMRMRD::ISMRMRD_CXDOUBLE)
767  for (size_t i = 0; i < n; i++)
768  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
769  else
770  for (size_t i = 0; i < n; i++)
771  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
772  }
773 
774  template<typename T>
775  void axpby_
776  (const ISMRMRD::Image<T>* ptr_x, complex_float_t a,
777  const void* vptr_y, complex_float_t b)
778  {
779  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
780  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
781  size_t n = ptr->getNumberOfDataElements();
782  size_t nx = ptr_x->getNumberOfDataElements();
783  size_t ny = ptr_y->getNumberOfDataElements();
784  if (!(n == nx && n == ny))
785  THROW("sizes mismatch in ImageWrap multiply");
786  const T* i = ptr_x->getDataPtr();
787  const T* j = ptr_y->getDataPtr();
788  T* k = ptr->getDataPtr();
789  if (b == complex_float_t(0.0))
790  for (size_t ii = 0; ii < n; i++, k++, ii++) {
791  complex_float_t u = (complex_float_t)*i;
792  xGadgetronUtilities::convert_complex(a*u, *k);
793  }
794  else
795  for (size_t ii = 0; ii < n; i++, j++, k++, ii++) {
796  complex_float_t u = (complex_float_t)*i;
797  complex_float_t v = (complex_float_t)*j;
798  complex_float_t w = a*u + b*v;
799  xGadgetronUtilities::convert_complex(w, *k);
800  }
801  }
802 
803  template<typename T>
804  void xapyb_
805  (const ISMRMRD::Image<T>* ptr_x, const void* vptr_a,
806  const void* vptr_y, const void* vptr_b, int a_type, int b_type)
807  {
808  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
809  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
810  size_t n = ptr->getNumberOfDataElements();
811  size_t nx = ptr_x->getNumberOfDataElements();
812  if (nx != n)
813  THROW("sizes mismatch in ImageWrap::xapyb: nx != n");
814  size_t ny = ptr_y->getNumberOfDataElements();
815  if (ny != n)
816  THROW("sizes mismatch in ImageWrap::xapyb: ny != n");
817  const T* ix = ptr_x->getDataPtr();
818  const T* iy = ptr_y->getDataPtr();
819  T* i = ptr->getDataPtr();
820  size_t na, ja;
821  size_t nb, jb;
822  const T* ia;
823  const T* ib;
824  if (a_type == 0) {
825  na = 1;
826  ia = (T*)vptr_a;
827  ja = 0;
828  }
829  else {
830  ISMRMRD::Image<T>* ptr_a = (ISMRMRD::Image<T>*)vptr_a;
831  na = ptr_a->getNumberOfDataElements();
832  if (na != n)
833  THROW("sizes mismatch in ImageWrap xapyb: na != n");
834  ia = ptr_a->getDataPtr();
835  ja = 1;
836  }
837  if (b_type == 0) {
838  nb = 1;
839  ib = (T*)vptr_b;
840  jb = 0;
841  }
842  else {
843  ISMRMRD::Image<T>* ptr_b = (ISMRMRD::Image<T>*)vptr_b;
844  nb = ptr_b->getNumberOfDataElements();
845  //if (nb != n)
846  // std::cout << nb << ' ' << n << '\n';
847  if (nb != n)
848  THROW("sizes mismatch in ImageWrap xapyb: nb != n");
849  ib = ptr_b->getDataPtr();
850  jb = 1;
851  }
852  for (size_t ii = 0; ii < n; ix++, ia += ja, iy++, ib += jb, i++, ii++) {
853  complex_float_t vx = (complex_float_t)*ix;
854  complex_float_t va = (complex_float_t)*ia;
855  complex_float_t vy = (complex_float_t)*iy;
856  complex_float_t vb = (complex_float_t)*ib;
857  complex_float_t v = vx * va + vy * vb;
858  xGadgetronUtilities::convert_complex(v, *i);
859  }
860  }
861 
862  template<typename T>
863  void binary_op_(const ISMRMRD::Image<T>* ptr_x, const void* vptr_y,
864  complex_float_t(*f)(complex_float_t, complex_float_t))
865  {
866  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
867  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
868  size_t nx = ptr_x->getNumberOfDataElements();
869  size_t ny = ptr_y->getNumberOfDataElements();
870  size_t n = ptr->getNumberOfDataElements();
871  if (!(n == nx && n == ny))
872  THROW("sizes mismatch in ImageWrap binary_op_");
873  const T* i = ptr_x->getDataPtr();
874  const T* j = ptr_y->getDataPtr();
875  T* k = ptr->getDataPtr();
876  size_t ii = 0;
877  for (; ii < n; i++, j++, k++, ii++) {
878  complex_float_t u = (complex_float_t)*i;
879  complex_float_t v = (complex_float_t)*j;
880  xGadgetronUtilities::convert_complex(f(u, v), *k);
881  }
882  }
883 
884  template<typename T>
885  void semibinary_op_(const ISMRMRD::Image<T>* ptr_x, complex_float_t y,
886  complex_float_t(*f)(complex_float_t, complex_float_t))
887  {
888  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
889  size_t nx = ptr_x->getNumberOfDataElements();
890  size_t n = ptr->getNumberOfDataElements();
891  if (n != nx)
892  THROW("sizes mismatch in ImageWrap semibinary_op_");
893  const T* i = ptr_x->getDataPtr();
894  T* k = ptr->getDataPtr();
895  size_t ii = 0;
896  for (; ii < n; i++, k++, ii++) {
897  complex_float_t x = (complex_float_t)*i;
898  xGadgetronUtilities::convert_complex(f(x, y), *k);
899  }
900  }
901 
902  template<typename T>
903  void unary_op_(const ISMRMRD::Image<T>* ptr_x,
904  complex_float_t(*f)(complex_float_t))
905  {
906  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
907  size_t nx = ptr_x->getNumberOfDataElements();
908  size_t n = ptr->getNumberOfDataElements();
909  if (n != nx)
910  THROW("sizes mismatch in ImageWrap semibinary_op_");
911  const T* i = ptr_x->getDataPtr();
912  T* k = ptr->getDataPtr();
913  size_t ii = 0;
914  for (; ii < n; i++, k++, ii++) {
915  complex_float_t x = (complex_float_t)*i;
916  xGadgetronUtilities::convert_complex(f(x), *k);
917  }
918  }
919 
920  template<typename T>
921  void dot_(const ISMRMRD::Image<T>* ptr_im, complex_float_t *z) const
922  {
923  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
924  const T* i;
925  const T* j;
926  *z = 0;
927  size_t ii = 0;
928  size_t n = ptr_im->getNumberOfDataElements();
929  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
930  i++, j++, ii++) {
931  complex_float_t u = (complex_float_t)*i;
932  complex_float_t v = (complex_float_t)*j;
933  *z += std::conj(v) * u;
934  }
935  }
936 
937  template<typename T>
938  void norm_(const ISMRMRD::Image<T>* ptr, float *r) const
939  {
940  const T* i;
941  *r = 0;
942  size_t ii = 0;
943  size_t n = ptr->getNumberOfDataElements();
944  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
945  complex_float_t a = (complex_float_t)*i;
946  *r += std::abs(std::conj(a) * a);
947  }
948  *r = std::sqrt(*r);
949  }
950 
951  template<typename T>
952  void sum_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
953  {
954  const T* i;
955  *s = 0;
956  size_t ii = 0;
957  size_t n = ptr->getNumberOfDataElements();
958  for (i = ptr->getDataPtr(); ii < n; i++, ii++)
959  *s += (complex_float_t)*i;
960  }
961 
962  template<typename T>
963  void max_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
964  {
965  const T* i;
966  *s = 0;
967  size_t ii = 0;
968  size_t n = ptr->getNumberOfDataElements();
969  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
970  complex_float_t si = (complex_float_t)*i;
971  float r = std::real(*s);
972  float ri = std::real(si);
973  if (ri > r)
974  *s = si;
975  }
976  }
977 
978  template<typename T>
979  void diff_(const ISMRMRD::Image<T>* ptr_im, float *s) const
980  {
981  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
982  const T* i;
983  const T* j;
984  *s = 0;
985  size_t ii = 0;
986  size_t n = ptr_im->getNumberOfDataElements();
987  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
988  i++, j++, ii++) {
989  complex_float_t a = (complex_float_t)*i;
990  complex_float_t b = (complex_float_t)*j;
991  *s += (float)std::abs(b - a);
992  }
993  }
994 
995  template<typename T>
996  void conjugate_(ISMRMRD::Image<T>* ptr)
997  {
998  T* i;
999  size_t ii = 0;
1000  size_t n = ptr->getNumberOfDataElements();
1001  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
1002  complex_float_t a = (complex_float_t)*i;
1003  xGadgetronUtilities::convert_complex(std::conj(a), *i);
1004  }
1005  }
1006 
1007  };
1008 }
1009 
1010 #ifdef _MSC_VER
1011 #pragma warning( pop )
1012 #endif
1013 
1014 #endif
sirf::ImageWrap::Iterator_const
Definition: gadgetron_image_wrap.h:179
sirf::ImageWrap::Iterator
Definition: gadgetron_image_wrap.h:108
sirf::ImageWrap
Wrapper for ISMRMRD::Image.
Definition: gadgetron_image_wrap.h:106
sirf::ImageWrap::is_complex
bool is_complex() const
Is the image wrap complex?
Definition: gadgetron_image_wrap.h:472
sirf::ImageWrap::get_data_type
ISMRMRD::ISMRMRD_DataTypes get_data_type() const
Get data type.
Definition: gadgetron_image_wrap.h:465
sirf::NumRef
Definition: ANumRef.h:140
iequals.h
Defines sirf::iequals.
sirf
Abstract data container.
Definition: GeometricalInfo.cpp:141
sirf::iequals
bool iequals(const std::string &a, const std::string &b)
Case insensitive string comparison, replaces boost::iequals.
Definition: iequals.cpp:7
xgadgetron_utilities.h
Various utilities used by SIRF Gadgetron extensions.
Generated by doxygen 1.9.1