SIRF  3.4.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  ISMRMRD::ImageHeader& head()
336  {
337  IMAGE_PROCESSING_SWITCH(type_, return get_head_ref_, ptr_);
338  }
339  const ISMRMRD::ImageHeader& head() const
340  {
341  IMAGE_PROCESSING_SWITCH_CONST(type_, return get_head_ref_const_, ptr_);
342  }
343  std::string attributes() const
344  {
345  std::string attr;
346  IMAGE_PROCESSING_SWITCH_CONST(type_, get_attr_, ptr_, attr);
347  return attr;
348  }
349  void set_attributes(const std::string& attr)
350  {
351  IMAGE_PROCESSING_SWITCH(type_, set_attr_, ptr_, attr);
352  }
353  void show_attributes() const
354  {
355  IMAGE_PROCESSING_SWITCH_CONST(type_, show_attr_, ptr_, 0);
356  }
357  void set_imtype(ISMRMRD::ISMRMRD_ImageTypes imtype)
358  {
359  IMAGE_PROCESSING_SWITCH(type_, set_imtype_, ptr_, imtype);
360  }
361  size_t get_dim(int* dim) const
362  {
363  IMAGE_PROCESSING_SWITCH_CONST(type_, get_dim_, ptr_, dim);
364  size_t n = dim[0];
365  n *= dim[1];
366  n *= dim[2];
367  n *= dim[3];
368  return n;
369  }
370  void get_data(float* data) const
371  {
372  //std::cout << "in get_data\n";
373  //std::cout << "trying new image wrap iterator...\n";
374  ImageWrap::Iterator_const i = begin_const();
375  ImageWrap::Iterator_const stop = end_const();
376  for (; i != stop; ++data, ++i) {
377  *data = *i;
378  }
379  //IMAGE_PROCESSING_SWITCH_CONST(type_, get_data_, ptr_, data);
380  }
381  void set_data(const float* data)
382  {
383  //std::cout << "in set_data\n";
384  for (ImageWrap::Iterator i = begin(); i != end(); ++i, ++data)
385  *i = *data;
386  //IMAGE_PROCESSING_SWITCH(type_, set_data_, ptr_, data);
387  }
388  void fill(float s)
389  {
390  for (ImageWrap::Iterator i = begin(); i != end(); ++i)
391  *i = s;
392  }
393  void scale(float s)
394  {
395  for (ImageWrap::Iterator i = begin(); i != end(); ++i)
396  *i /= s;
397  }
398  void get_complex_data(complex_float_t* data) const
399  {
400  //std::cout << "in get_complex_data\n";
401  //std::cout << "trying new const image wrap iterator...\n";
402  ImageWrap::Iterator_const i = begin_const();
403  ImageWrap::Iterator_const stop = end_const();
404  for (; i != stop; ++data, ++i) {
405  *data = (*i).complex_float();
406  }
407  //IMAGE_PROCESSING_SWITCH_CONST(type_, get_complex_data_, ptr_, data);
408  }
409 
410  void set_complex_data(const complex_float_t* data)
411  {
412  //std::cout << "in set_complex_data\n";
413  //std::cout << "trying new image wrap iterator...\n";
414  ImageWrap::Iterator i = begin();
415  ImageWrap::Iterator stop = end();
416  for (; i != stop; ++i, ++data) {
417  *i = *data;
418  }
419  //IMAGE_PROCESSING_SWITCH(type_, set_complex_data_, ptr_, data);
420  }
421 
422  gadgetron::shared_ptr<ImageWrap> abs() const
423  {
424  return real("abs");
425  }
426 
427  gadgetron::shared_ptr<ImageWrap> real(const std::string& way = "real") const
428  {
429  int dim[4];
430  get_dim(dim);
431  ISMRMRD::Image<float>* ptr_im = new ISMRMRD::Image<float>(dim[0], dim[1], dim[2], dim[3]);
432  ISMRMRD::Image<float>& im = *ptr_im;
433 
434  ISMRMRD::ImageHeader header = head();
435  header.data_type = ISMRMRD::ISMRMRD_FLOAT;
436  header.image_type = ISMRMRD::ISMRMRD_IMTYPE_MAGNITUDE;
437  header.image_series_index = 0;
438  im.setHead(header);
439 
440  im.setAttributeString(attributes());
441 
442  ImageWrap::Iterator_const i = begin_const();
443  ImageWrap::Iterator_const stop = end_const();
444  float* data = im.getDataPtr();
445  if (sirf::iequals(way, "abs"))
446  for (; i != stop; ++data, ++i) {
447  *data = std::abs((*i).complex_float());
448  }
449  else if (sirf::iequals(way, "real"))
450  for (; i != stop; ++data, ++i) {
451  *data = std::real((*i).complex_float());
452  }
453  else
454  THROW("unknown conversion to real specified in ImageWrap::real");
455 
456  ImageWrap* ptr_iw = new ImageWrap(ISMRMRD::ISMRMRD_FLOAT, ptr_im);
457  return gadgetron::shared_ptr<ImageWrap>(ptr_iw);
458  }
459 
461  ISMRMRD::ISMRMRD_DataTypes get_data_type() const
462  {
463  ISMRMRD::ISMRMRD_DataTypes data_type;
464  IMAGE_PROCESSING_SWITCH_CONST(type_, get_data_type_, ptr_, &data_type);
465  return data_type;
466  }
468  bool is_complex() const
469  {
470  ISMRMRD::ISMRMRD_DataTypes data_type = this->get_data_type();
471  if (data_type == ISMRMRD::ISMRMRD_DataTypes::ISMRMRD_CXFLOAT ||
472  data_type == ISMRMRD::ISMRMRD_DataTypes::ISMRMRD_CXDOUBLE) {
473  return true;
474  }
475  else
476  return false;
477  }
478  void write(ISMRMRD::Dataset& dataset) const
479  {
480  IMAGE_PROCESSING_SWITCH_CONST(type_, write_, ptr_, dataset);
481  }
482  void read(ISMRMRD::Dataset& dataset, const char* var, int ind)
483  {
484  IMAGE_PROCESSING_SWITCH(type_, read_, ptr_, dataset, var, ind, &ptr_);
485  }
486  //void axpby(complex_float_t a, const ImageWrap& x, complex_float_t b)
487  //{
488  // IMAGE_PROCESSING_SWITCH(type_, axpby_, x.ptr_image(), a, b);
489  //}
490  void axpby(complex_float_t a, const ImageWrap& x, complex_float_t b,
491  const ImageWrap& y)
492  {
493  xapyb(x, a, y, b);
494  //IMAGE_PROCESSING_SWITCH(type_, axpby_, x.ptr_image(), a, y.ptr_image(), b);
495  }
496  void xapyb(const ImageWrap& x, complex_float_t a,
497  const ImageWrap& y, complex_float_t b)
498  {
499  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), &a,
500  y.ptr_image(), &b, 0, 0);
501  }
502  void xapyb(const ImageWrap& x, complex_float_t a,
503  const ImageWrap& y, const ImageWrap& b)
504  {
505  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), &a,
506  y.ptr_image(), &b, 0, 1);
507  }
508  void xapyb(const ImageWrap& x, const ImageWrap& a,
509  const ImageWrap& y, complex_float_t b)
510  {
511  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), &a,
512  y.ptr_image(), &b, 1, 0);
513  }
514  void xapyb(const ImageWrap& x, const ImageWrap& a,
515  const ImageWrap& y, const ImageWrap& b)
516  {
517  IMAGE_PROCESSING_SWITCH(type_, xapyb_, x.ptr_image(), a.ptr_image(),
518  y.ptr_image(), b.ptr_image(), 1, 1);
519  }
520  void multiply(const ImageWrap& x)
521  {
522  IMAGE_PROCESSING_SWITCH(type_, multiply_, x.ptr_image());
523  }
524  void multiply(const ImageWrap& x, const ImageWrap& y)
525  {
526  IMAGE_PROCESSING_SWITCH(type_, multiply__, x.ptr_image(), y.ptr_image());
527  }
528  void divide(const ImageWrap& x)
529  {
530  IMAGE_PROCESSING_SWITCH(type_, divide_, x.ptr_image());
531  }
532  void divide(const ImageWrap& x, const ImageWrap& y)
533  {
534  IMAGE_PROCESSING_SWITCH(type_, divide__, x.ptr_image(), y.ptr_image());
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  float diff(ImageWrap& iw) const
549  {
550  float s;
551  IMAGE_PROCESSING_SWITCH_CONST(type_, diff_, iw.ptr_image(), &s);
552  return s;
553  }
554  void conjugate()
555  {
556  IMAGE_PROCESSING_SWITCH(type_, conjugate_, ptr_);
557  }
558 
559  private:
560  int type_;
561  void* ptr_;
562  mutable gadgetron::shared_ptr<Iterator> begin_;
563  mutable gadgetron::shared_ptr<Iterator> end_;
564  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
565  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
566 
567  ImageWrap& operator=(const ImageWrap& iw)
568  {
569  //type_ = iw.type();
570  //IMAGE_PROCESSING_SWITCH(type_, copy_, iw.ptr_image());
571  return *this;
572  }
573 
574  template<typename T>
575  void get_data_parameters_
576  (const ISMRMRD::Image<T>* ptr_im, size_t *n, unsigned int *dsize,
577  char** data) const
578  {
579  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
580  unsigned int dim[4];
581  dim[0] = im.getMatrixSizeX();
582  dim[1] = im.getMatrixSizeY();
583  dim[2] = im.getMatrixSizeZ();
584  dim[3] = im.getNumberOfChannels();
585  *data = (char*)im.getDataPtr();
586  *dsize = sizeof(T);
587  *n = dim[0];
588  *n *= dim[1];
589  *n *= dim[2];
590  *n *= dim[3];
591  //*n = ptr_im->getDataSize()/(*dsize);
592  }
593 
594  template<typename T>
595  void copy_(const ISMRMRD::Image<T>* ptr_im)
596  {
597  type_ = ptr_im->getDataType();
598  ptr_ = (void*)new ISMRMRD::Image<T>(*ptr_im);
599  }
600 
601  template<typename T>
602  ISMRMRD::ImageHeader& get_head_ref_(ISMRMRD::Image<T>* ptr_im)
603  {
604  return ptr_im->getHead();
605  }
606 
607  template<typename T>
608  const ISMRMRD::ImageHeader& get_head_ref_const_(const ISMRMRD::Image<T>* ptr_im) const
609  {
610  return ptr_im->getHead();
611  }
612 
613  template<typename T>
614  void set_imtype_(ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_ImageTypes type)
615  {
616  ptr_im->setImageType(type);
617  }
618 
619  template<typename T>
620  void get_size_(const ISMRMRD::Image<T>* ptr_im, size_t& size) const
621  {
622  size = ptr_im->getDataSize();
623  }
624 
625  template<typename T>
626  void get_attr_(const ISMRMRD::Image<T>* ptr_im, std::string& attr) const
627  {
628  ptr_im->getAttributeString(attr);
629  }
630 
631  template<typename T>
632  void set_attr_(ISMRMRD::Image<T>* ptr_im, const std::string& attr)
633  {
634  ptr_im->setAttributeString(attr);
635  }
636 
637  template<typename T>
638  void show_attr_(const ISMRMRD::Image<T>* ptr_im, int i) const
639  {
640  const ISMRMRD::Image<T>& im = *ptr_im;
641  size_t meta_attrib_length = im.getAttributeStringLength();
642  std::cout << "meta_attrib_length: " << meta_attrib_length << '\n';
643  std::string meta_attrib(meta_attrib_length + 1, 0);
644  im.getAttributeString(meta_attrib);
645 
646  //std::cout << "attributes:" << std::endl << meta_attrib << std::endl;
647  if (meta_attrib_length > 0) {
648  ISMRMRD::MetaContainer mc;
649  ISMRMRD::deserialize(meta_attrib.c_str(), mc);
650  for (auto it = mc.begin(); it != mc.end(); ++it) {
651  auto name = it->first.c_str();
652  std::cout << name << ":\n";
653  size_t l = mc.length(name);
654  for (int j = 0; j <l; j++)
655  std::cout << mc.as_str(name, j) << '\n';
656  }
657  std::stringstream meta_xml;
658  ISMRMRD::serialize(mc, meta_xml);
659  std::cout << meta_xml.str().c_str() << '\n';
660  }
661  }
662 
663  template<typename T>
664  void write_
665  (const ISMRMRD::Image<T>* ptr_im, ISMRMRD::Dataset& dataset) const
666  {
667  //std::cout << "appending image..." << std::endl;
668  const ISMRMRD::Image<T>& im = *ptr_im;
669  std::stringstream ss;
670  ss << "image_" << im.getHead().image_series_index;
671  std::string image_varname = ss.str();
672  {
673  Mutex mtx;
674  mtx.lock();
675  dataset.appendImage(image_varname, im);
676  mtx.unlock();
677  }
678  }
679 
680  template<typename T>
681  void read_
682  (const ISMRMRD::Image<T>* ptr,
683  ISMRMRD::Dataset& dataset, const char* var, int index,
684  void** ptr_ptr)
685  {
686  ISMRMRD::Image < T >* ptr_im = new ISMRMRD::Image < T > ;
687  *ptr_ptr = (void*)ptr_im;
688  ISMRMRD::Image<T>& im = *ptr_im;
689  dataset.readImage(var, index, im);
690  //int status = ismrmrd_read_image(&dataset, var, (uint32_t)index, &(im.im));
691  }
692 
693  template<typename T>
694  void get_dim_(const ISMRMRD::Image<T>* ptr_im, int* dim) const
695  {
696  const ISMRMRD::Image<T>& im = *(const ISMRMRD::Image<T>*)ptr_im;
697  dim[0] = im.getMatrixSizeX();
698  dim[1] = im.getMatrixSizeY();
699  dim[2] = im.getMatrixSizeZ();
700  dim[3] = im.getNumberOfChannels();
701  }
702  template<typename T>
703  void get_data_type_(const ISMRMRD::Image<T>* ptr_im, ISMRMRD::ISMRMRD_DataTypes* data_type_ptr) const
704  {
705  const ISMRMRD::Image<T>& im = *static_cast<const ISMRMRD::Image<T>*>(ptr_im);
706  *data_type_ptr = im.getDataType();
707  }
708 
709  template<typename T>
710  void get_data_(const ISMRMRD::Image<T>* ptr_im, float* data) const
711  {
712  const ISMRMRD::Image<T>& im = *ptr_im;
713  const T* ptr = im.getDataPtr();
714  size_t n = im.getNumberOfDataElements();
715  for (size_t i = 0; i < n; i++)
716  data[i] = std::real(ptr[i]);
717  }
718 
719  template<typename T>
720  void set_data_(ISMRMRD::Image<T>* ptr_im, const float* data)
721  {
722  ISMRMRD::Image<T>& im = *ptr_im;
723  T* ptr = im.getDataPtr();
724  size_t n = im.getNumberOfDataElements();
725  for (size_t i = 0; i < n; i++)
726  ptr[i] = (T)data[i];
727  }
728 
729  template<typename T>
730  void get_complex_data_
731  (const ISMRMRD::Image<T>* ptr_im, complex_float_t* data) const
732  {
733  const ISMRMRD::Image<T>& im = *ptr_im;
734  const T* ptr = im.getDataPtr();
735  size_t n = im.getNumberOfDataElements();
736  for (size_t i = 0; i < n; i++)
737  data[i] = complex_float_t(ptr[i]);
738  }
739 
740  template<typename T>
741  void set_complex_data_
742  (ISMRMRD::Image<T>* ptr_im, const complex_float_t* data)
743  {
744  ISMRMRD::Image<T>& im = *ptr_im;
745  ISMRMRD::ISMRMRD_DataTypes type = im.getDataType();
746  T* ptr = im.getDataPtr();
747  size_t n = im.getNumberOfDataElements();
748  if (type == ISMRMRD::ISMRMRD_CXFLOAT || type == ISMRMRD::ISMRMRD_CXDOUBLE)
749  for (size_t i = 0; i < n; i++)
750  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
751  else
752  for (size_t i = 0; i < n; i++)
753  xGadgetronUtilities::convert_complex(data[i], ptr[i]);
754  }
755 
756  //template<typename T>
757  //void axpby_
758  // (const ISMRMRD::Image<T>* ptr_x, complex_float_t a,
759  // complex_float_t b)
760  //{
761  // ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)ptr_;
762  // size_t nx = ptr_x->getNumberOfDataElements();
763  // size_t ny = ptr_y->getNumberOfDataElements();
764  // if (nx != ny)
765  // THROW("sizes mismatch in ImageWrap multiply");
766  // const T* i = ptr_x->getDataPtr();
767  // T* j = ptr_y->getDataPtr();
768  // if (b == complex_float_t(0.0))
769  // for (size_t ii = 0; ii < nx; i++, j++, ii++) {
770  // complex_float_t u = (complex_float_t)*i;
771  // xGadgetronUtilities::convert_complex(a*u, *j);
772  // }
773  // else
774  // for (size_t ii = 0; ii < nx; i++, j++, ii++) {
775  // complex_float_t u = (complex_float_t)*i;
776  // complex_float_t v = (complex_float_t)*j;
777  // xGadgetronUtilities::convert_complex(a*u + b*v, *j);
778  // }
779  //}
780 
781  template<typename T>
782  void axpby_
783  (const ISMRMRD::Image<T>* ptr_x, complex_float_t a,
784  const void* vptr_y, complex_float_t b)
785  {
786  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
787  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
788  size_t n = ptr->getNumberOfDataElements();
789  size_t nx = ptr_x->getNumberOfDataElements();
790  size_t ny = ptr_y->getNumberOfDataElements();
791  if (!(n == nx && n == ny))
792  THROW("sizes mismatch in ImageWrap multiply");
793  const T* i = ptr_x->getDataPtr();
794  const T* j = ptr_y->getDataPtr();
795  T* k = ptr->getDataPtr();
796  if (b == complex_float_t(0.0))
797  for (size_t ii = 0; ii < n; i++, k++, ii++) {
798  complex_float_t u = (complex_float_t)*i;
799  xGadgetronUtilities::convert_complex(a*u, *k);
800  }
801  else
802  for (size_t ii = 0; ii < n; i++, j++, k++, ii++) {
803  complex_float_t u = (complex_float_t)*i;
804  complex_float_t v = (complex_float_t)*j;
805  complex_float_t w = a*u + b*v;
806  xGadgetronUtilities::convert_complex(w, *k);
807  }
808  }
809 
810  template<typename T>
811  void xapyb_
812  (const ISMRMRD::Image<T>* ptr_x, const void* vptr_a,
813  const void* vptr_y, const void* vptr_b, int a_type, int b_type)
814  {
815  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
816  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
817  size_t n = ptr->getNumberOfDataElements();
818  size_t nx = ptr_x->getNumberOfDataElements();
819  if (nx != n)
820  THROW("sizes mismatch in ImageWrap::xapyb: nx != n");
821  size_t ny = ptr_y->getNumberOfDataElements();
822  if (ny != n)
823  THROW("sizes mismatch in ImageWrap::xapyb: ny != n");
824  const T* ix = ptr_x->getDataPtr();
825  const T* iy = ptr_y->getDataPtr();
826  T* i = ptr->getDataPtr();
827  size_t na, ja;
828  size_t nb, jb;
829  const T* ia;
830  const T* ib;
831  if (a_type == 0) {
832  na = 1;
833  ia = (T*)vptr_a;
834  ja = 0;
835  }
836  else {
837  ISMRMRD::Image<T>* ptr_a = (ISMRMRD::Image<T>*)vptr_a;
838  na = ptr_a->getNumberOfDataElements();
839  if (na != n)
840  THROW("sizes mismatch in ImageWrap xapyb: na != n");
841  ia = ptr_a->getDataPtr();
842  ja = 1;
843  }
844  if (b_type == 0) {
845  nb = 1;
846  ib = (T*)vptr_b;
847  jb = 0;
848  }
849  else {
850  ISMRMRD::Image<T>* ptr_b = (ISMRMRD::Image<T>*)vptr_b;
851  nb = ptr_b->getNumberOfDataElements();
852  if (nb != n)
853  THROW("sizes mismatch in ImageWrap xapyb: nb != n");
854  ib = ptr_b->getDataPtr();
855  jb = 1;
856  }
857  for (size_t ii = 0; ii < n; ix++, ia += ja, iy++, ib += jb, i++, ii++) {
858  complex_float_t vx = (complex_float_t)*ix;
859  complex_float_t va = (complex_float_t)*ia;
860  complex_float_t vy = (complex_float_t)*iy;
861  complex_float_t vb = (complex_float_t)*ib;
862  complex_float_t v = vx * va + vy * vb;
863  xGadgetronUtilities::convert_complex(v, *i);
864  }
865  }
866 
867  template<typename T>
868  void multiply_(const ISMRMRD::Image<T>* ptr_x)
869  {
870  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)ptr_;
871  size_t nx = ptr_x->getNumberOfDataElements();
872  size_t ny = ptr_y->getNumberOfDataElements();
873  if (nx != ny)
874  THROW("sizes mismatch in ImageWrap multiply");
875  const T* i = ptr_x->getDataPtr();
876  T* j = ptr_y->getDataPtr();
877  size_t ii = 0;
878  for (; ii < nx; i++, j++, ii++) {
879  complex_float_t u = (complex_float_t)*i;
880  complex_float_t v = (complex_float_t)*j;
881  xGadgetronUtilities::convert_complex(u*v, *j);
882  }
883  }
884 
885  template<typename T>
886  void multiply__(const ISMRMRD::Image<T>* ptr_x, const void* vptr_y)
887  {
888  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
889  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
890  size_t nx = ptr_x->getNumberOfDataElements();
891  size_t ny = ptr_y->getNumberOfDataElements();
892  size_t n = ptr->getNumberOfDataElements();
893  if (!(n == nx && n == ny))
894  THROW("sizes mismatch in ImageWrap multiply");
895  const T* i = ptr_x->getDataPtr();
896  const T* j = ptr_y->getDataPtr();
897  T* k = ptr->getDataPtr();
898  size_t ii = 0;
899  for (; ii < n; i++, j++, k++, ii++) {
900  complex_float_t u = (complex_float_t)*i;
901  complex_float_t v = (complex_float_t)*j;
902  xGadgetronUtilities::convert_complex(u*v, *k);
903  }
904  }
905 
906  template<typename T>
907  void divide_(const ISMRMRD::Image<T>* ptr_x)
908  {
909  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)ptr_;
910  size_t nx = ptr_x->getNumberOfDataElements();
911  size_t ny = ptr_y->getNumberOfDataElements();
912  if (nx != ny)
913  THROW("sizes mismatch in ImageWrap divide 1");
914  const T* i = ptr_x->getDataPtr();
915  T* j = ptr_y->getDataPtr();
916  size_t ii = 0;
917  for (; ii < nx; i++, j++, ii++) {
918  complex_float_t u = (complex_float_t)*i;
919  complex_float_t v = (complex_float_t)*j;
920  xGadgetronUtilities::convert_complex(v / u, *j);
921  }
922  }
923 
924  template<typename T>
925  void divide__(const ISMRMRD::Image<T>* ptr_x, const void* vptr_y)
926  {
927  ISMRMRD::Image<T>* ptr = (ISMRMRD::Image<T>*)ptr_;
928  ISMRMRD::Image<T>* ptr_y = (ISMRMRD::Image<T>*)vptr_y;
929  size_t nx = ptr_x->getNumberOfDataElements();
930  size_t ny = ptr_y->getNumberOfDataElements();
931  size_t n = ptr->getNumberOfDataElements();
932  if (!(n == nx && n == ny))
933  THROW("sizes mismatch in ImageWrap divide 2");
934  const T* i = ptr_x->getDataPtr();
935  const T* j = ptr_y->getDataPtr();
936  T* k = ptr->getDataPtr();
937  size_t ii = 0;
938  for (; ii < n; i++, j++, k++, ii++) {
939  complex_float_t u = (complex_float_t)*i;
940  complex_float_t v = (complex_float_t)*j;
941  xGadgetronUtilities::convert_complex(u / v, *k);
942  }
943  }
944 
945  template<typename T>
946  void dot_(const ISMRMRD::Image<T>* ptr_im, complex_float_t *z) const
947  {
948  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
949  const T* i;
950  const T* j;
951  *z = 0;
952  size_t ii = 0;
953  size_t n = ptr_im->getNumberOfDataElements();
954  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
955  i++, j++, ii++) {
956  complex_float_t u = (complex_float_t)*i;
957  complex_float_t v = (complex_float_t)*j;
958  *z += std::conj(v) * u;
959  }
960  }
961 
962  template<typename T>
963  void norm_(const ISMRMRD::Image<T>* ptr, float *r) const
964  {
965  const T* i;
966  *r = 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 a = (complex_float_t)*i;
971  *r += std::abs(std::conj(a) * a);
972  }
973  *r = std::sqrt(*r);
974  }
975 
976  template<typename T>
977  void diff_(const ISMRMRD::Image<T>* ptr_im, float *s) const
978  {
979  const ISMRMRD::Image<T>* ptr = (const ISMRMRD::Image<T>*)ptr_;
980  const T* i;
981  const T* j;
982  *s = 0;
983  size_t ii = 0;
984  size_t n = ptr_im->getNumberOfDataElements();
985  for (i = ptr->getDataPtr(), j = ptr_im->getDataPtr(); ii < n;
986  i++, j++, ii++) {
987  complex_float_t a = (complex_float_t)*i;
988  complex_float_t b = (complex_float_t)*j;
989  *s += (float)std::abs(b - a);
990  }
991  }
992 
993  template<typename T>
994  void conjugate_(ISMRMRD::Image<T>* ptr)
995  {
996  T* i;
997  size_t ii = 0;
998  size_t n = ptr->getNumberOfDataElements();
999  for (i = ptr->getDataPtr(); ii < n; i++, ii++) {
1000  complex_float_t a = (complex_float_t)*i;
1001  xGadgetronUtilities::convert_complex(std::conj(a), *i);
1002  }
1003  }
1004 
1005  };
1006 }
1007 
1008 #ifdef _MSC_VER
1009 #pragma warning( pop )
1010 #endif
1011 
1012 #endif
sirf::iequals
bool iequals(const std::string &a, const std::string &b)
Case insensitive string comparison, replaces boost::iequals.
Definition: iequals.cpp:7
sirf::Mutex
Definition: xgadgetron_utilities.h:107
iequals.h
Case insensitive string comparison sirf::iequals.
sirf::NumRef
Definition: ANumRef.h:140
sirf
Abstract data container.
Definition: GeometricalInfo.cpp:141
sirf::ImageWrap::Iterator
Definition: gadgetron_image_wrap.h:108
sirf::ImageWrap::is_complex
bool is_complex() const
Is the image wrap complex?
Definition: gadgetron_image_wrap.h:468
sirf::ImageWrap::Iterator_const
Definition: gadgetron_image_wrap.h:179
xgadgetron_utilities.h
Various utilities used by SIRF Gadgetron extensions.
sirf::ImageWrap::get_data_type
ISMRMRD::ISMRMRD_DataTypes get_data_type() const
Get data type.
Definition: gadgetron_image_wrap.h:461
sirf::ImageWrap
Wrapper for ISMRMRD::Image.
Definition: gadgetron_image_wrap.h:106
Generated by   doxygen 1.8.13