SIRF  3.7.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  complex_float_t min() const
561  {
562  complex_float_t s;
563  IMAGE_PROCESSING_SWITCH_CONST(type_, min_, ptr_, &s);
564  return s;
565  }
566  float diff(ImageWrap& iw) const
567  {
568  float s;
569  IMAGE_PROCESSING_SWITCH_CONST(type_, diff_, iw.ptr_image(), &s);
570  return s;
571  }
572  void conjugate()
573  {
574  IMAGE_PROCESSING_SWITCH(type_, conjugate_, ptr_);
575  }
576 
577  private:
578  int type_;
579  void* ptr_;
580  mutable gadgetron::shared_ptr<Iterator> begin_;
581  mutable gadgetron::shared_ptr<Iterator> end_;
582  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
583  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
584 
585  ImageWrap& operator=(const ImageWrap& iw)
586  {
587  //type_ = iw.type();
588  //IMAGE_PROCESSING_SWITCH(type_, copy_, iw.ptr_image());
589  return *this;
590  }
591 
592  template<typename T>
593  void get_data_parameters_
594  (const ISMRMRD::Image<T>* ptr_im, size_t *n, unsigned int *dsize,
595  char** data) const
596  {
597  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
598  unsigned int dim[4];
599  dim[0] = im.getMatrixSizeX();
600  dim[1] = im.getMatrixSizeY();
601  dim[2] = im.getMatrixSizeZ();
602  dim[3] = im.getNumberOfChannels();
603  *data = (char*)im.getDataPtr();
604  *dsize = sizeof(T);
605  *n = dim[0];
606  *n *= dim[1];
607  *n *= dim[2];
608  *n *= dim[3];
609  //*n = ptr_im->getDataSize()/(*dsize);
610  }
611 
612  template<typename T>
613  size_t num_data_elm_(const ISMRMRD::Image<T>* ptr) const
614  {
615  return ptr->getNumberOfDataElements();
616  }
617 
618  template<typename T>
619  void copy_(const ISMRMRD::Image<T>* ptr_im)
620  {
621  type_ = ptr_im->getDataType();
622  ptr_ = (void*)new ISMRMRD::Image<T>(*ptr_im);
623  }
624 
625  template<typename T>
626  ISMRMRD::ImageHeader& get_head_ref_(ISMRMRD::Image<T>* ptr_im)
627  {
628  return ptr_im->getHead();
629  }
630 
631  template<typename T>
632  const ISMRMRD::ImageHeader& get_head_ref_const_(const ISMRMRD::Image<T>* ptr_im) const
633  {
634  return ptr_im->getHead();
635  }
636 
637  template<typename T>
638  void set_imtype_(ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_ImageTypes type)
639  {
640  ptr_im->setImageType(type);
641  }
642 
643  template<typename T>
644  void get_size_(const ISMRMRD::Image<T>* ptr_im, size_t& size) const
645  {
646  size = ptr_im->getDataSize();
647  }
648 
649  template<typename T>
650  void get_attr_(const ISMRMRD::Image<T>* ptr_im, std::string& attr) const
651  {
652  ptr_im->getAttributeString(attr);
653  }
654 
655  template<typename T>
656  void set_attr_(ISMRMRD::Image<T>* ptr_im, const std::string& attr)
657  {
658  ptr_im->setAttributeString(attr);
659  }
660 
661  template<typename T>
662  void show_attr_(const ISMRMRD::Image<T>* ptr_im, int i) const
663  {
664  const ISMRMRD::Image<T>& im = *ptr_im;
665  size_t meta_attrib_length = im.getAttributeStringLength();
666  std::cout << "meta_attrib_length: " << meta_attrib_length << '\n';
667  std::string meta_attrib(meta_attrib_length + 1, 0);
668  im.getAttributeString(meta_attrib);
669 
670  //std::cout << "attributes:" << std::endl << meta_attrib << std::endl;
671  if (meta_attrib_length > 0) {
672  ISMRMRD::MetaContainer mc;
673  ISMRMRD::deserialize(meta_attrib.c_str(), mc);
674  for (auto it = mc.begin(); it != mc.end(); ++it) {
675  auto name = it->first.c_str();
676  std::cout << name << ":\n";
677  size_t l = mc.length(name);
678  for (int j = 0; j <l; j++)
679  std::cout << mc.as_str(name, j) << '\n';
680  }
681  std::stringstream meta_xml;
682  ISMRMRD::serialize(mc, meta_xml);
683  std::cout << meta_xml.str().c_str() << '\n';
684  }
685  }
686 
687  template<typename T>
688  void write_
689  (const ISMRMRD::Image<T>* ptr_im, ISMRMRD::Dataset& dataset) const
690  {
691  //std::cout << "appending image..." << std::endl;
692  const ISMRMRD::Image<T>& im = *ptr_im;
693  std::stringstream ss;
694  ss << "image_" << im.getHead().image_series_index;
695  std::string image_varname = ss.str();
696  {
697  Mutex mtx;
698  mtx.lock();
699  dataset.appendImage(image_varname, im);
700  mtx.unlock();
701  }
702  }
703 
704  template<typename T>
705  void read_
706  (const ISMRMRD::Image<T>* ptr,
707  ISMRMRD::Dataset& dataset, const char* var, int index,
708  void** ptr_ptr)
709  {
710  ISMRMRD::Image < T >* ptr_im = new ISMRMRD::Image < T > ;
711  *ptr_ptr = (void*)ptr_im;
712  ISMRMRD::Image<T>& im = *ptr_im;
713  dataset.readImage(var, index, im);
714  //int status = ismrmrd_read_image(&dataset, var, (uint32_t)index, &(im.im));
715  }
716 
717  template<typename T>
718  void get_dim_(const ISMRMRD::Image<T>* ptr_im, int* dim) const
719  {
720  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
721  dim[0] = im.getMatrixSizeX();
722  dim[1] = im.getMatrixSizeY();
723  dim[2] = im.getMatrixSizeZ();
724  dim[3] = im.getNumberOfChannels();
725  }
726  template<typename T>
727  void get_data_type_(const ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_DataTypes* data_type_ptr) const
728  {
729  const ISMRMRD::Image<T>& im = *static_cast<const ISMRMRD::Image<T>*>(ptr_im);
730  *data_type_ptr = im.getDataType();
731  }
732 
733  template<typename T>
734  void get_data_(const ISMRMRD::Image<T>* ptr_im, float* data) const
735  {
736  const ISMRMRD::Image<T>& im = *ptr_im;
737  const T* ptr = im.getDataPtr();
738  size_t n = im.getNumberOfDataElements();
739  for (size_t i = 0; i < n; i++)
740  data[i] = std::real(ptr[i]);
741  }
742 
743  template<typename T>
744  void set_data_(ISMRMRD::Image<T>* ptr_im, const float* data)
745  {
746  ISMRMRD::Image<T>& im = *ptr_im;
747  T* ptr = im.getDataPtr();
748  size_t n = im.getNumberOfDataElements();
749  for (size_t i = 0; i < n; i++)
750  ptr[i] = (T)data[i];
751  }
752 
753  template<typename T>
754  void get_complex_data_
755  (const ISMRMRD::Image<T>* ptr_im, complex_float_t* data) const
756  {
757  const ISMRMRD::Image<T>& im = *ptr_im;
758  const T* ptr = im.getDataPtr();
759  size_t n = im.getNumberOfDataElements();
760  for (size_t i = 0; i < n; i++)
761  data[i] = complex_float_t(ptr[i]);
762  }
763 
764  template<typename T>
765  void set_complex_data_
766  (ISMRMRD::Image<T>* ptr_im, const complex_float_t* data)
767  {
768  ISMRMRD::Image<T>& im = *ptr_im;
769  ISMRMRD::ISMRMRD_DataTypes type = im.getDataType();
770  T* ptr = im.getDataPtr();
771  size_t n = im.getNumberOfDataElements();
772  if (type == ISMRMRD::ISMRMRD_CXFLOAT || type == ISMRMRD::ISMRMRD_CXDOUBLE)
773  for (size_t i = 0; i < n; i++)
774  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
775  else
776  for (size_t i = 0; i < n; i++)
777  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
778  }
779 
780  template<typename T>
781  void axpby_
782  (const ISMRMRD::Image<T>* ptr_x, complex_float_t a,
783  const void* vptr_y, complex_float_t b)
784  {
785  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
786  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
787  size_t n = ptr->getNumberOfDataElements();
788  size_t nx = ptr_x->getNumberOfDataElements();
789  size_t ny = ptr_y->getNumberOfDataElements();
790  if (!(n == nx && n == ny))
791  THROW("sizes mismatch in ImageWrap multiply");
792  const T* i = ptr_x->getDataPtr();
793  const T* j = ptr_y->getDataPtr();
794  T* k = ptr->getDataPtr();
795  if (b == complex_float_t(0.0))
796  for (size_t ii = 0; ii < n; i++, k++, ii++) {
797  complex_float_t u = (complex_float_t)*i;
798  xGadgetronUtilities::convert_complex(a*u, *k);
799  }
800  else
801  for (size_t ii = 0; ii < n; i++, j++, k++, ii++) {
802  complex_float_t u = (complex_float_t)*i;
803  complex_float_t v = (complex_float_t)*j;
804  complex_float_t w = a*u + b*v;
805  xGadgetronUtilities::convert_complex(w, *k);
806  }
807  }
808 
809  template<typename T>
810  void xapyb_
811  (const ISMRMRD::Image<T>* ptr_x, const void* vptr_a,
812  const void* vptr_y, const void* vptr_b, int a_type, int b_type)
813  {
814  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
815  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
816  size_t n = ptr->getNumberOfDataElements();
817  size_t nx = ptr_x->getNumberOfDataElements();
818  if (nx != n)
819  THROW("sizes mismatch in ImageWrap::xapyb: nx != n");
820  size_t ny = ptr_y->getNumberOfDataElements();
821  if (ny != n)
822  THROW("sizes mismatch in ImageWrap::xapyb: ny != n");
823  const T* ix = ptr_x->getDataPtr();
824  const T* iy = ptr_y->getDataPtr();
825  T* i = ptr->getDataPtr();
826  size_t na, ja;
827  size_t nb, jb;
828  const T* ia;
829  const T* ib;
830  if (a_type == 0) {
831  na = 1;
832  ia = (T*)vptr_a;
833  ja = 0;
834  }
835  else {
836  ISMRMRD::Image<T>* ptr_a = (ISMRMRD::Image<T>*)vptr_a;
837  na = ptr_a->getNumberOfDataElements();
838  if (na != n)
839  THROW("sizes mismatch in ImageWrap xapyb: na != n");
840  ia = ptr_a->getDataPtr();
841  ja = 1;
842  }
843  if (b_type == 0) {
844  nb = 1;
845  ib = (T*)vptr_b;
846  jb = 0;
847  }
848  else {
849  ISMRMRD::Image<T>* ptr_b = (ISMRMRD::Image<T>*)vptr_b;
850  nb = ptr_b->getNumberOfDataElements();
851  //if (nb != n)
852  // std::cout << nb << ' ' << n << '\n';
853  if (nb != n)
854  THROW("sizes mismatch in ImageWrap xapyb: nb != n");
855  ib = ptr_b->getDataPtr();
856  jb = 1;
857  }
858  for (size_t ii = 0; ii < n; ix++, ia += ja, iy++, ib += jb, i++, ii++) {
859  complex_float_t vx = (complex_float_t)*ix;
860  complex_float_t va = (complex_float_t)*ia;
861  complex_float_t vy = (complex_float_t)*iy;
862  complex_float_t vb = (complex_float_t)*ib;
863  complex_float_t v = vx * va + vy * vb;
864  xGadgetronUtilities::convert_complex(v, *i);
865  }
866  }
867 
868  template<typename T>
869  void binary_op_(const ISMRMRD::Image<T>* ptr_x, const void* vptr_y,
870  complex_float_t(*f)(complex_float_t, complex_float_t))
871  {
872  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
873  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
874  size_t nx = ptr_x->getNumberOfDataElements();
875  size_t ny = ptr_y->getNumberOfDataElements();
876  size_t n = ptr->getNumberOfDataElements();
877  if (!(n == nx && n == ny))
878  THROW("sizes mismatch in ImageWrap binary_op_");
879  const T* i = ptr_x->getDataPtr();
880  const T* j = ptr_y->getDataPtr();
881  T* k = ptr->getDataPtr();
882  size_t ii = 0;
883  for (; ii < n; i++, j++, k++, ii++) {
884  complex_float_t u = (complex_float_t)*i;
885  complex_float_t v = (complex_float_t)*j;
886  xGadgetronUtilities::convert_complex(f(u, v), *k);
887  }
888  }
889 
890  template<typename T>
891  void semibinary_op_(const ISMRMRD::Image<T>* ptr_x, complex_float_t y,
892  complex_float_t(*f)(complex_float_t, complex_float_t))
893  {
894  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
895  size_t nx = ptr_x->getNumberOfDataElements();
896  size_t n = ptr->getNumberOfDataElements();
897  if (n != nx)
898  THROW("sizes mismatch in ImageWrap semibinary_op_");
899  const T* i = ptr_x->getDataPtr();
900  T* k = ptr->getDataPtr();
901  size_t ii = 0;
902  for (; ii < n; i++, k++, ii++) {
903  complex_float_t x = (complex_float_t)*i;
904  xGadgetronUtilities::convert_complex(f(x, y), *k);
905  }
906  }
907 
908  template<typename T>
909  void unary_op_(const ISMRMRD::Image<T>* ptr_x,
910  complex_float_t(*f)(complex_float_t))
911  {
912  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
913  size_t nx = ptr_x->getNumberOfDataElements();
914  size_t n = ptr->getNumberOfDataElements();
915  if (n != nx)
916  THROW("sizes mismatch in ImageWrap semibinary_op_");
917  const T* i = ptr_x->getDataPtr();
918  T* k = ptr->getDataPtr();
919  size_t ii = 0;
920  for (; ii < n; i++, k++, ii++) {
921  complex_float_t x = (complex_float_t)*i;
922  xGadgetronUtilities::convert_complex(f(x), *k);
923  }
924  }
925 
926  template<typename T>
927  void dot_(const ISMRMRD::Image<T>* ptr_im, complex_float_t *z) const
928  {
929  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
930  const T* i;
931  const T* j;
932  *z = 0;
933  size_t ii = 0;
934  size_t n = ptr_im->getNumberOfDataElements();
935  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
936  i++, j++, ii++) {
937  complex_float_t u = (complex_float_t)*i;
938  complex_float_t v = (complex_float_t)*j;
939  *z += std::conj(v) * u;
940  }
941  }
942 
943  template<typename T>
944  void norm_(const ISMRMRD::Image<T>* ptr, float *r) const
945  {
946  const T* i;
947  *r = 0;
948  size_t ii = 0;
949  size_t n = ptr->getNumberOfDataElements();
950  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
951  complex_float_t a = (complex_float_t)*i;
952  *r += std::abs(std::conj(a) * a);
953  }
954  *r = std::sqrt(*r);
955  }
956 
957  template<typename T>
958  void sum_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
959  {
960  const T* i;
961  *s = 0;
962  size_t ii = 0;
963  size_t n = ptr->getNumberOfDataElements();
964  for (i = ptr->getDataPtr(); ii < n; i++, ii++)
965  *s += (complex_float_t)*i;
966  }
967 
968  template<typename T>
969  void max_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
970  {
971  const T* i;
972  *s = 0;
973  size_t ii = 0;
974  size_t n = ptr->getNumberOfDataElements();
975  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
976  complex_float_t si = (complex_float_t)*i;
977  float r = std::real(*s);
978  float ri = std::real(si);
979  if (ii == 0 || ri > r)
980  *s = si;
981  }
982  }
983 
984  template<typename T>
985  void min_(const ISMRMRD::Image<T>* ptr, complex_float_t* s) const
986  {
987  const T* i;
988  *s = 0;
989  size_t ii = 0;
990  size_t n = ptr->getNumberOfDataElements();
991  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
992  complex_float_t si = (complex_float_t)*i;
993  float r = std::real(*s);
994  float ri = std::real(si);
995  if (ii == 0 || ri < r)
996  *s = si;
997  }
998  }
999 
1000  template<typename T>
1001  void diff_(const ISMRMRD::Image<T>* ptr_im, float *s) const
1002  {
1003  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
1004  const T* i;
1005  const T* j;
1006  *s = 0;
1007  size_t ii = 0;
1008  size_t n = ptr_im->getNumberOfDataElements();
1009  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
1010  i++, j++, ii++) {
1011  complex_float_t a = (complex_float_t)*i;
1012  complex_float_t b = (complex_float_t)*j;
1013  *s += (float)std::abs(b - a);
1014  }
1015  }
1016 
1017  template<typename T>
1018  void conjugate_(ISMRMRD::Image<T>* ptr)
1019  {
1020  T* i;
1021  size_t ii = 0;
1022  size_t n = ptr->getNumberOfDataElements();
1023  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
1024  complex_float_t a = (complex_float_t)*i;
1025  xGadgetronUtilities::convert_complex(std::conj(a), *i);
1026  }
1027  }
1028 
1029  };
1030 }
1031 
1032 #ifdef _MSC_VER
1033 #pragma warning( pop )
1034 #endif
1035 
1036 #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 base class for SIRF image data.
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