SIRF  3.8.0
gadgetron_data_containers.h
Go to the documentation of this file.
1 /*
2 SyneRBI Synergistic Image Reconstruction Framework (SIRF)
3 Copyright 2015 - 2023 Rutherford Appleton Laboratory STFC
4 Copyright 2020 - 2023 University College London
5 Copyright 2020 - 2023 Physikalisch-Technische Bundesanstalt (PTB)
6 
7 This is software developed for the Collaborative Computational
8 Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
9 (http://www.ccpsynerbi.ac.uk/).
10 
11 Licensed under the Apache License, Version 2.0 (the "License");
12 you may not use this file except in compliance with the License.
13 You may obtain a copy of the License at
14 http://www.apache.org/licenses/LICENSE-2.0
15 Unless required by applicable law or agreed to in writing, software
16 distributed under the License is distributed on an "AS IS" BASIS,
17 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 See the License for the specific language governing permissions and
19 limitations under the License.
20 
21 */
22 
33 #ifndef GADGETRON_DATA_CONTAINERS
34 #define GADGETRON_DATA_CONTAINERS
35 
36 #include <string>
37 #include <vector>
38 #include <tuple>
39 
40 //#include <boost/algorithm/string.hpp>
41 
42 #include <ismrmrd/ismrmrd.h>
43 #include <ismrmrd/dataset.h>
44 
45 #include "sirf/common/DataContainer.h"
46 #include "sirf/common/ImageData.h"
47 #include "sirf/common/multisort.h"
48 #include "sirf/Gadgetron/cgadgetron_shared_ptr.h"
49 #include "sirf/Gadgetron/gadgetron_image_wrap.h"
50 
51 #include "sirf/iUtilities/LocalisedException.h"
52 
53 //#define SIRF_DYNAMIC_CAST(T, X, Y) T& X = (T&)Y
54 #define SIRF_DYNAMIC_CAST(T, X, Y) T& X = dynamic_cast<T&>(Y)
55 
62 namespace sirf {
63 
64  class FourierEncoding;
65  class CartesianFourierEncoding;
66 #if GADGETRON_TOOLBOXES_AVAILABLE
67  class RPEFourierEncoding;
68 #endif
69 
70  class AcquisitionsInfo {
71  public:
72  AcquisitionsInfo(std::string data = "") : data_(data)
73  {
74  if (data.empty())
75  have_header_ = false;
76  else {
77  deserialize();
78  have_header_ = true;
79  }
80  }
81  AcquisitionsInfo& operator=(std::string data)
82  {
83  data_ = data;
84  if (data.empty())
85  have_header_ = false;
86  else {
87  deserialize();
88  have_header_ = true;
89  }
90 
91  return *this;
92  }
93  const char* c_str() const { return data_.c_str(); }
94  operator std::string&() { return data_; }
95  operator const std::string&() const { return data_; }
96  bool empty() const { return data_.empty(); }
97  const ISMRMRD::IsmrmrdHeader& get_IsmrmrdHeader() const
98  {
99  if (!have_header_)
100  deserialize();
101  return header_;
102  }
103 
104  private:
105  void deserialize() const
106  {
107  if (!this->empty())
108  { header_ = ISMRMRD::IsmrmrdHeader();
109  ISMRMRD::deserialize(data_.c_str(), header_);
110  }
111  have_header_ = true;
112  }
113  std::string data_;
114  mutable ISMRMRD::IsmrmrdHeader header_;
115  mutable bool have_header_;
116  };
117 
132  class IgnoreMask {
133  public:
134  IgnoreMask(unsigned long long int mask = (1 << 18)) :
135  ignore_(mask), max_(8*sizeof(mask)) {}
136  void set(unsigned long long int mask)
137  {
138  ignore_ = mask;
139  }
140  void ignore(unsigned int i)
141  {
142  if (i < 1 || i > max_)
143  return;
144  unsigned long long int one = 1;
145  ignore_ = ignore_ | (one << (i - 1));
146  }
147  void ignore_not(unsigned int i)
148  {
149  if (i < 1 || i > max_)
150  return;
151  unsigned long long int one = 1;
152  ignore_ = ignore_ & ~(one << (i - 1));
153  }
154  bool bit(unsigned int i) const
155  {
156  if (i < 1 || i > max_)
157  return true;
158  unsigned long long int one = 1;
159  return ignore_ & (one << (i - 1));
160  }
161  unsigned long long int bits() const
162  {
163  return ignore_;
164  }
165  bool ignored(unsigned long long int bits) const
166  {
167  return bits & ignore_;
168  }
169  std::string bits_string() const
170  {
171  unsigned int size = max_;
172  unsigned long long int one = 1;
173  unsigned long long int bitmask = (one << (size - 1));
174  std::stringstream str;
175  for (unsigned int i = 0; i < size; i++) {
176  if (ignore_ & (bitmask >> i))
177  str << '1';
178  else
179  str << '0';
180  if ((i + 1) % 4 == 0)
181  str << ' ';
182  }
183  str << '\n';
184  return str.str();
185  }
186  private:
187  unsigned long long int ignore_;
188  unsigned int max_;
189  };
190 
202  class KSpaceSubset
203  {
204  static int const num_kspace_dims_ = 7 + ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_INTS;
205 
206  public:
207 
208  typedef std::array<int, num_kspace_dims_> TagType;
209  typedef std::vector<int> SetType;
210 
211  KSpaceSubset(){
212  for(int i=0; i<num_kspace_dims_; ++i)
213  this->tag_[i] = -1;
214  }
215 
216  KSpaceSubset(TagType tag){
217  this->tag_ = tag;
218  this->idx_set_ = {};
219  }
220 
221  KSpaceSubset(TagType tag, SetType idx_set){
222  this->tag_ = tag;
223  this->idx_set_ = idx_set;
224  }
225 
226  TagType get_tag(void) const {return tag_;}
227  SetType get_idx_set(void) const {return idx_set_;}
228  void add_idx_to_set(size_t const idx){this->idx_set_.push_back(idx);}
229 
230  bool is_first_set() const {
231  bool is_first= (tag_[0] == 0);
232  if(is_first)
233  {
234  for(int dim=2; dim<num_kspace_dims_; ++dim)
235  is_first *= (tag_[dim] == 0);
236  }
237  return is_first;
238  }
239 
240  static void print_tag(const TagType& tag);
241  static void print_acquisition_tag(ISMRMRD::Acquisition acq);
242 
244 
247  static TagType get_tag_from_acquisition(ISMRMRD::Acquisition acq);
248 
250 
253  static TagType get_tag_from_img(const CFImage& img);
254 
255  private:
257 
261  TagType tag_;
262 
264 
268  SetType idx_set_;
269  };
270 
276  class MRAcquisitionData : public DataContainer {
277  public:
278  // static methods
279 
280  // ISMRMRD acquisitions algebra: acquisitions viewed as vectors of
281  // acquisition data
282  static void binary_op
283  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y,
284  complex_float_t (*f)(complex_float_t, complex_float_t));
285  static void semibinary_op
286  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y,
287  complex_float_t(*f)(complex_float_t, complex_float_t));
288  static void unary_op
289  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y,
290  complex_float_t(*f)(complex_float_t));
291  // y := a x + b y
292  static void axpby
293  (complex_float_t a, const ISMRMRD::Acquisition& acq_x,
294  complex_float_t b, ISMRMRD::Acquisition& acq_y);
295  static void xapyb
296  (const ISMRMRD::Acquisition& acq_x, complex_float_t a,
297  ISMRMRD::Acquisition& acq_y, complex_float_t b);
298  static void xapyb
299  (const ISMRMRD::Acquisition& acq_x, complex_float_t a,
300  ISMRMRD::Acquisition& acq_y, const ISMRMRD::Acquisition& acq_b);
301  static void xapyb
302  (const ISMRMRD::Acquisition& acq_x, const ISMRMRD::Acquisition& acq_a,
303  ISMRMRD::Acquisition& acq_y, const ISMRMRD::Acquisition& acq_b);
304 
305  // the inner (l2) product of x and y
306  static complex_float_t dot
307  (const ISMRMRD::Acquisition& acq_x, const ISMRMRD::Acquisition& acq_y);
308  // the sum of the elements of x
309  static complex_float_t sum(const ISMRMRD::Acquisition& acq_x);
310  // the value of the element of x with the largest real part
311  static complex_float_t max(const ISMRMRD::Acquisition& acq_x);
312  static complex_float_t min(const ISMRMRD::Acquisition& acq_x);
313  // elementwise multiplication
314  // y := x .* y
315  static void multiply
316  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
317  // multiply by scalar
318  // y := x * y
319  static void multiply
320  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
321  // add scalar
322  // y := x + y
323  static void add
324  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
325  // elementwise division
326  // y := x ./ y
327  static void divide
328  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
329  // elementwise maximum
330  // y := std::real(x) > std::real(y) ? x : y
331  static void maximum
332  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
333  static void maximum
334  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
335  // elementwise minimum
336  // y := std::real(x) < std::real(y) ? x : y
337  static void minimum
338  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
339  static void minimum
340  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
341  // y := pow(x, y)
342  static void power
343  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
344  static void power
345  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
346  // y := exp(x)
347  static void exp
348  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
349  // y := log(x)
350  static void log
351  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
352  // y := sqrt(x)
353  static void sqrt
354  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
355  // y := sign(x) (x < 0: -1, x == 0: 0, x > 0: 1)
356  static void sign
357  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
358  // l2 norm of x
359  // y := abs(x)
360  static void abs
361  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
362  static float norm(const ISMRMRD::Acquisition& acq_x);
363 
364  // type and dimension of an ISMRMRD::Acquisition parameter
365  static void ismrmrd_par_info(const char* par, int* output)
366  {
367  // default type: int, dimension: 1
368  output[0] = 0;
369  output[1] = 1;
370 
371  // check parameter type
372  if (sirf::iequals(par, "sample_time_us") ||
373  sirf::iequals(par, "position") ||
374  sirf::iequals(par, "read_dir") ||
375  sirf::iequals(par, "phase_dir") ||
376  sirf::iequals(par, "slice_dir") ||
377  sirf::iequals(par, "patient_table_position") ||
378  sirf::iequals(par, "user_float"))
379  output[0] = 1; // float
380 
381  // check parameter dimension
382  if (sirf::iequals(par, "physiology_time_stamp"))
383  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
384  else if (sirf::iequals(par, "channel_mask"))
385  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
386  else if (sirf::iequals(par, "position") ||
387  sirf::iequals(par, "read_dir") ||
388  sirf::iequals(par, "phase_dir") ||
389  sirf::iequals(par, "slice_dir") ||
390  sirf::iequals(par, "patient_table_position"))
391  output[1] = 3;
392  else if (sirf::iequals(par, "user_int") ||
393  sirf::iequals(par, "idx_user"))
394  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_INTS;
395  else if (sirf::iequals(par, "user_float"))
396  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_FLOATS;
397  }
398  // value of an ISMRMRD::Acquisition int parameter
399  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
400  const char* name, unsigned long long int* v)
401  {
402  if (sirf::iequals(name, "version"))
403  *v = ((unsigned int)acq.version());
404  else if (sirf::iequals(name, "flags"))
405  *v = ((unsigned long long int)acq.flags());
406  else if (sirf::iequals(name, "measurement_uid"))
407  *v = ((unsigned int)acq.measurement_uid());
408  else if (sirf::iequals(name, "scan_counter"))
409  *v = ((unsigned int)acq.scan_counter());
410  else if (sirf::iequals(name, "acquisition_time_stamp"))
411  *v = ((unsigned int)acq.acquisition_time_stamp());
412  else if (sirf::iequals(name, "number_of_samples"))
413  *v = ((unsigned int)acq.number_of_samples());
414  else if (sirf::iequals(name, "available_channels"))
415  *v = ((unsigned int)acq.available_channels());
416  else if (sirf::iequals(name, "active_channels"))
417  *v = ((unsigned int)acq.active_channels());
418  else if (sirf::iequals(name, "discard_pre"))
419  *v = ((unsigned int)acq.discard_pre());
420  else if (sirf::iequals(name, "discard_post"))
421  *v = ((unsigned int)acq.discard_post());
422  else if (sirf::iequals(name, "center_sample"))
423  *v = ((unsigned int)acq.center_sample());
424  else if (sirf::iequals(name, "encoding_space_ref"))
425  *v = ((unsigned int)acq.encoding_space_ref());
426  else if (sirf::iequals(name, "trajectory_dimensions"))
427  *v = ((unsigned int)acq.trajectory_dimensions());
428  else if (sirf::iequals(name, "kspace_encode_step_1"))
429  *v = ((unsigned int)acq.idx().kspace_encode_step_1);
430  else if (sirf::iequals(name, "kspace_encode_step_2"))
431  *v = ((unsigned int)acq.idx().kspace_encode_step_2);
432  else if (sirf::iequals(name, "average"))
433  *v = ((unsigned int)acq.idx().average);
434  else if (sirf::iequals(name, "slice"))
435  *v = ((unsigned int)acq.idx().slice);
436  else if (sirf::iequals(name, "contrast"))
437  *v = ((unsigned int)acq.idx().contrast);
438  else if (sirf::iequals(name, "phase"))
439  *v = ((unsigned int)acq.idx().phase);
440  else if (sirf::iequals(name, "repetition"))
441  *v = ((unsigned int)acq.idx().repetition);
442  else if (sirf::iequals(name, "set"))
443  *v = ((unsigned int)acq.idx().set);
444  else if (sirf::iequals(name, "segment"))
445  *v = ((unsigned int)acq.idx().segment);
446  else if (sirf::iequals(name, "physiology_time_stamp")) {
447  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
448  const uint32_t* pts = acq.physiology_time_stamp();
449  for (int i = 0; i < n; i++)
450  v[i] = (unsigned int)pts[i];
451  }
452  else if (sirf::iequals(name, "channel_mask")) {
453  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
454  const uint64_t* pts = acq.channel_mask();
455  for (int i = 0; i < n; i++)
456  v[i] = (unsigned long long int)pts[i];
457  }
458  }
459  // value of an ISMRMRD::Acquisition float parameter
460  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
461  const char* name, float* v)
462  {
463  if (sirf::iequals(name, "sample_time_us"))
464  *v = acq.sample_time_us();
465  else if (sirf::iequals(name, "position")) {
466  float* u = acq.position();
467  for (int i = 0; i < 3; i++)
468  v[i] = u[i];
469  }
470  else if (sirf::iequals(name, "read_dir")) {
471  float* u = acq.read_dir();
472  for (int i = 0; i < 3; i++)
473  v[i] = u[i];
474  }
475  else if (sirf::iequals(name, "phase_dir")) {
476  float* u = acq.phase_dir();
477  for (int i = 0; i < 3; i++)
478  v[i] = u[i];
479  }
480  else if (sirf::iequals(name, "slice_dir")) {
481  float* u = acq.slice_dir();
482  for (int i = 0; i < 3; i++)
483  v[i] = u[i];
484  }
485  else if (sirf::iequals(name, "patient_table_position")) {
486  float* u = acq.patient_table_position();
487  for (int i = 0; i < 3; i++)
488  v[i] = u[i];
489  }
490  }
497  void set_encoding_limits(const std::string& name,
498  const std::tuple<unsigned short,unsigned short,unsigned short> min_max_ctr)
499  {
500  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
501  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
502 
503  ISMRMRD::Limit limit;
504  limit.minimum = std::get<0>(min_max_ctr);
505  limit.maximum = std::get<1>(min_max_ctr);
506  limit.center = std::get<2>(min_max_ctr);
507 
508  if(sirf::iequals(name, "kspace_encoding_step_1"))
509  enc_limits.kspace_encoding_step_1.get() = limit;
510  else if(sirf::iequals(name, "kspace_encoding_step_2"))
511  enc_limits.kspace_encoding_step_2.get() = limit;
512  else if(sirf::iequals(name, "average"))
513  enc_limits.average.get() = limit;
514  else if(sirf::iequals(name, "slice"))
515  enc_limits.slice.get() = limit;
516  else if(sirf::iequals(name, "contrast"))
517  enc_limits.contrast.get() = limit;
518  else if(sirf::iequals(name, "phase"))
519  enc_limits.phase.get() = limit;
520  else if(sirf::iequals(name, "repetition"))
521  enc_limits.repetition.get() = limit;
522  else if(sirf::iequals(name, "set"))
523  enc_limits.set.get() = limit;
524  else if(sirf::iequals(name, "segment"))
525  enc_limits.segment.get() = limit;
526  else
527  throw std::runtime_error("You passed a name that is not an encoding limit.");
528 
529  hdr.encoding[0].encodingLimits = enc_limits;
530  std::stringstream serialised_hdr;
531  ISMRMRD::serialize(hdr, serialised_hdr);
532  this->set_acquisitions_info(AcquisitionsInfo(serialised_hdr.str()));
533  }
534 
535  std::tuple<unsigned short,unsigned short,unsigned short>
536  get_encoding_limits(const std::string& name) const
537  {
538  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
539  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
540 
541  ISMRMRD::Limit limit;
542 
543  if(sirf::iequals(name, "kspace_encoding_step_1"))
544  limit = enc_limits.kspace_encoding_step_1.get();
545  else if(sirf::iequals(name, "kspace_encoding_step_2"))
546  limit = enc_limits.kspace_encoding_step_2.get();
547  else if(sirf::iequals(name, "average"))
548  limit = enc_limits.average.get();
549  else if(sirf::iequals(name, "slice"))
550  limit = enc_limits.slice.get();
551  else if(sirf::iequals(name, "contrast"))
552  limit = enc_limits.contrast.get();
553  else if(sirf::iequals(name, "phase"))
554  limit = enc_limits.phase.get();
555  else if(sirf::iequals(name, "repetition"))
556  limit = enc_limits.repetition.get();
557  else if(sirf::iequals(name, "set"))
558  limit = enc_limits.set.get();
559  else if(sirf::iequals(name, "segment"))
560  limit = enc_limits.segment.get();
561  else
562  throw std::runtime_error("You passed a name that is not an encoding limit.");
563 
564  return std::make_tuple(limit.minimum, limit.maximum, limit.center);
565  }
566 
567  // abstract methods
568 
569  virtual void empty() = 0;
570  virtual void take_over(MRAcquisitionData&) = 0;
571 
572  // the number of acquisitions in the container
573  virtual unsigned int number() const = 0;
574 
575  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
576  get_acquisition_sptr(unsigned int num) = 0;
577  virtual int get_acquisition(unsigned int,
578  ISMRMRD::Acquisition&) const = 0;
579  virtual void set_acquisition(unsigned int,
580  ISMRMRD::Acquisition&) = 0;
581  virtual void append_acquisition(ISMRMRD::Acquisition& acq) = 0;
582 
583  virtual void copy_acquisitions_info(const MRAcquisitionData& ac) = 0;
584  virtual void copy_acquisitions_data(const MRAcquisitionData& ac) = 0;
585 
586  // 'export' constructors: workaround for creating 'ABC' objects
587  virtual gadgetron::unique_ptr<MRAcquisitionData> new_acquisitions_container() = 0;
588  virtual MRAcquisitionData*
589  same_acquisitions_container(const AcquisitionsInfo& info) const = 0;
590 
591  virtual void set_data(const complex_float_t* z, int all = 1) = 0;
592  virtual void get_data(complex_float_t* z, int all = 1);
593 
594  virtual void set_user_floats(float const * const z, int const idx);
595 
596  virtual bool is_complex() const
597  {
598  return true;
599  }
600 
601  // acquisition data algebra
603  virtual void sum(void* ptr) const;
604  virtual void max(void* ptr) const;
605  virtual void min(void* ptr) const;
606  virtual void dot(const DataContainer& dc, void* ptr) const;
607  complex_float_t dot(const DataContainer& a_x)
608  {
609  complex_float_t z;
610  dot(a_x, &z);
611  return z;
612  }
613  virtual void axpby(
614  const void* ptr_a, const DataContainer& a_x,
615  const void* ptr_b, const DataContainer& a_y);
616  virtual void xapyb(
617  const DataContainer& a_x, const DataContainer& a_a,
618  const DataContainer& a_y, const DataContainer& a_b);
619  virtual void xapyb(
620  const DataContainer& a_x, const void* ptr_a,
621  const DataContainer& a_y, const void* ptr_b)
622  {
623  axpby(ptr_a, a_x, ptr_b, a_y);
624  }
625  virtual void xapyb(
626  const DataContainer& a_x, const void* ptr_a,
627  const DataContainer& a_y, const DataContainer& a_b);
628  virtual void multiply(const DataContainer& x, const DataContainer& y);
629  virtual void divide(const DataContainer& x, const DataContainer& y);
630  virtual void maximum(const DataContainer& x, const DataContainer& y);
631  virtual void minimum(const DataContainer& x, const DataContainer& y);
632  virtual void power(const DataContainer& x, const DataContainer& y);
633  virtual void multiply(const DataContainer& x, const void* y);
634  virtual void add(const DataContainer& x, const void* ptr_y);
635  virtual void maximum(const DataContainer& x, const void* y);
636  virtual void minimum(const DataContainer& x, const void* y);
637  virtual void power(const DataContainer& x, const void* y);
638  virtual void exp(const DataContainer& x);
639  virtual void log(const DataContainer& x);
640  virtual void sqrt(const DataContainer& x);
641  virtual void sign(const DataContainer& x);
642  virtual void abs(const DataContainer& x);
643  virtual float norm() const;
644 
645  virtual void write(const std::string &filename) const;
646 
647  // regular methods
648  void binary_op(const DataContainer& a_x, const DataContainer& a_y,
649  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
650  void semibinary_op(const DataContainer& a_x, complex_float_t y,
651  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&, complex_float_t));
652  void unary_op(const DataContainer& a_x,
653  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
654 
655  AcquisitionsInfo acquisitions_info() const { return acqs_info_; }
656  void set_acquisitions_info(std::string info) { acqs_info_ = info; }
657  void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;}
658  IgnoreMask ignore_mask() const { return ignore_mask_; }
659  void set_ignore_mask(IgnoreMask ignore_mask) const { ignore_mask_= ignore_mask; }
660 
661  ISMRMRD::TrajectoryType get_trajectory_type() const;
662 
663  void set_trajectory_type(const ISMRMRD::TrajectoryType type);
664 
665  void set_trajectory(const uint16_t traj_dim, float* traj)
666  {
667  ISMRMRD::Acquisition acq;
668  for(int i=0; i<number(); ++i)
669  {
670  get_acquisition(i, acq);
671  const uint16_t num_samples = acq.number_of_samples();
672  const uint16_t num_channels = acq.active_channels();
673  acq.resize(num_samples,num_channels,traj_dim);
674  int const offset = i*traj_dim*num_samples;
675  acq.setTraj(traj + offset);
676  set_acquisition(i, acq);
677  }
678  }
679 
680  gadgetron::unique_ptr<MRAcquisitionData> clone() const
681  {
682  return gadgetron::unique_ptr<MRAcquisitionData>(this->clone_impl());
683  }
684 
685  bool undersampled() const;
686  int get_acquisitions_dimensions(size_t ptr_dim) const;
687  void get_kspace_dimensions(std::vector<size_t>& dims) const;
688  uint16_t get_trajectory_dimensions(void) const;
689 
690  void sort();
691  void sort_by_time();
692  bool sorted() const { return sorted_; }
693  void set_sorted(bool sorted) { sorted_ = sorted; }
694 
696 
700  std::vector<KSpaceSubset::SetType > get_kspace_order() const;
701 
703  std::vector<KSpaceSubset> get_kspace_sorting() const { return this->sorting_; }
704 
706 
712  void organise_kspace();
713 
714  virtual std::vector<int> get_flagged_acquisitions_index
715  (const std::vector<ISMRMRD::ISMRMRD_AcquisitionFlags> flags) const;
716  virtual std::vector<int> get_slice_encoding_index
717  (const unsigned kspace_encode_step_2) const;
718 
719  virtual void get_subset(MRAcquisitionData& subset, const std::vector<int> subset_idx) const;
720  virtual void set_subset(const MRAcquisitionData &subset, const std::vector<int> subset_idx);
721 
722  std::vector<int> index() { return index_; }
723  const std::vector<int>& index() const { return index_; }
724 
725  int index(int i) const
726  {
727  const std::size_t ni = index_.size();
728  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni)
729  || static_cast<unsigned>(i) >= number())
730  THROW("Aquisition number is out of range");
731  if (ni > 0)
732  return index_[i];
733  else
734  return i;
735  }
736 
747  void read(const std::string& filename_ismrmrd_with_ext, int all = 0);
748 
749  protected:
750  bool sorted_ = false;
751  std::vector<int> index_;
752  std::vector<KSpaceSubset> sorting_;
753  AcquisitionsInfo acqs_info_;
754 
755  mutable IgnoreMask ignore_mask_; //= IgnoreMask();
756 
757  // new MRAcquisitionData objects will be created from this template
758  // using same_acquisitions_container()
759  static gadgetron::shared_ptr<MRAcquisitionData> acqs_templ_;
760 
761  virtual MRAcquisitionData* clone_impl() const = 0;
762 
763  private:
764 
765  };
766 
775  class AcquisitionsVector : public MRAcquisitionData {
776  public:
777  AcquisitionsVector(const std::string& filename_with_ext, int all = 0, IgnoreMask ignore_mask = IgnoreMask())
778  {
779  this->set_ignore_mask(ignore_mask);
780  this->read(filename_with_ext, all);
781  }
782 
783  AcquisitionsVector(const AcquisitionsInfo& info = AcquisitionsInfo(), IgnoreMask ignore_mask = IgnoreMask())
784  {
785  this->set_ignore_mask(ignore_mask);
786  acqs_info_ = info;
787  }
788  virtual void empty();
789  virtual void take_over(MRAcquisitionData& ad) {}
790  virtual unsigned int number() const { return (unsigned int)acqs_.size(); }
791  virtual unsigned int items() const { return (unsigned int)acqs_.size(); }
792  virtual void append_acquisition(ISMRMRD::Acquisition& acq)
793  {
794  acqs_.push_back(gadgetron::shared_ptr<ISMRMRD::Acquisition>
795  (new ISMRMRD::Acquisition(acq)));
796  }
797  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
798  get_acquisition_sptr(unsigned int num)
799  {
800  int ind = index(num);
801  return acqs_[ind];
802  }
803  virtual int get_acquisition(unsigned int num,
804  ISMRMRD::Acquisition& acq) const
805  {
806  IgnoreMask ignore_mask = this->ignore_mask();
807  int ind = index(num);
808  acq = *acqs_[ind];
809  if (ignore_mask.ignored(acq.flags()))
810  return 0;
811  return 1;
812  }
813  virtual void set_acquisition(unsigned int num, ISMRMRD::Acquisition& acq)
814  {
815  int ind = index(num);
816  *acqs_[ind] = acq;
817  }
818  virtual void copy_acquisitions_info(const MRAcquisitionData& ac)
819  {
820  acqs_info_ = ac.acquisitions_info();
821  }
822  virtual void copy_acquisitions_data(const MRAcquisitionData& ac);
823  virtual void set_data(const complex_float_t* z, int all = 1);
824 
825  virtual AcquisitionsVector* same_acquisitions_container
826  (const AcquisitionsInfo& info) const
827  {
828  return new AcquisitionsVector(info, ignore_mask_);
829  }
830  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
831  {
832  DataContainer* ptr = new AcquisitionsVector(acqs_info_, ignore_mask_);
833  return new ObjectHandle<DataContainer>
834  (gadgetron::shared_ptr<DataContainer>(ptr));
835  }
836  virtual gadgetron::unique_ptr<MRAcquisitionData>
837  new_acquisitions_container()
838  {
839  return gadgetron::unique_ptr<MRAcquisitionData>
840  (new AcquisitionsVector(acqs_info_, ignore_mask_));
841  }
842 
843  private:
844  std::vector<gadgetron::shared_ptr<ISMRMRD::Acquisition> > acqs_;
845  virtual AcquisitionsVector* clone_impl() const;
846  virtual void conjugate_impl();
847  };
848 
855  class ISMRMRDImageData : public ImageData {
856  public:
857  //ISMRMRDImageData(ISMRMRDImageData& id, const char* attr,
858  //const char* target); //does not build, have to be in the derived class
859 
860  virtual void empty() = 0;
861  virtual unsigned int number() const = 0;
862  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
863  (unsigned int im_num) = 0;
864  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
865  (unsigned int im_num) const = 0;
866 // virtual ImageWrap& image_wrap(unsigned int im_num) = 0;
867 // virtual const ImageWrap& image_wrap(unsigned int im_num) const = 0;
868  virtual void append(int image_data_type, void* ptr_image) = 0;
869  virtual void append(const ImageWrap& iw) = 0;
870  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw) = 0;
871  virtual gadgetron::shared_ptr<ISMRMRDImageData> abs() const = 0;
872  virtual gadgetron::shared_ptr<ISMRMRDImageData> real() const = 0;
873  virtual void clear_data() = 0;
874  virtual void set_image_type(int imtype) = 0;
875  virtual void get_data(complex_float_t* data) const;
876  virtual void set_data(const complex_float_t* data);
877  virtual void get_real_data(float* data) const;
878  virtual void set_real_data(const float* data);
879  virtual int read(std::string filename, std::string variable = "", int iv = -1);
880  virtual void write(const std::string &filename, const std::string &groupname, const bool dicom) const;
881  virtual void write(const std::string &filename) const
882  {
883  size_t size = filename.size();
884  std::string suff = filename.substr(size - 4, 4);
885  if (suff == std::string(".dcm")) {
886  std::string prefix = filename.substr(0, size - 4);
887  this->write(prefix, "", true);
888  }
889  else {
890  auto found = filename.find_last_of("/\\");
891  auto slash_found = found;
892  if (found == std::string::npos)
893  found = filename.find_last_of(".");
894  else
895  found = filename.substr(found + 1).find_last_of(".");
896  if (found == std::string::npos)
897  this->write(filename + ".h5", "", false);
898  else {
899  std::string ext = filename.substr(slash_found + found + 1);
900  if (ext == std::string(".h5"))
901  this->write(filename, "", false);
902  else
903  std::cerr << "WARNING: writing ISMRMRD images to "
904  << ext << "-files not implemented, "
905  << "please convert to Nifti images\n";
906  }
907  }
908  }
909  virtual Dimensions dimensions() const
910  {
911  Dimensions dim;
912  const ImageWrap& iw = image_wrap(0);
913  int d[5];
914  iw.get_dim(d);
915  dim["x"] = d[0];
916  dim["y"] = d[1];
917  dim["z"] = d[2];
918  dim["c"] = d[3];
919  dim["n"] = number();
920  return dim;
921  }
922  virtual void get_image_dimensions(unsigned int im_num, int* dim) const
923  {
924  if (im_num >= number())
925  dim[0] = dim[1] = dim[2] = dim[3] = 0;
926  const ImageWrap& iw = image_wrap(im_num);
927  iw.get_dim(dim);
928  }
929  bool check_dimension_consistency() const
930  {
931  size_t const num_dims = 4;
932  std::vector<int> first_img_dims(num_dims), temp_img_dims(num_dims);
933 
934  this->get_image_dimensions(0, &first_img_dims[0]);
935 
936  bool dims_match = true;
937  for(int i=1; i<number(); ++i)
938  {
939  this->get_image_dimensions(0, &temp_img_dims[0]);
940  dims_match *= (first_img_dims == temp_img_dims);
941  }
942  return dims_match;
943  }
944  virtual gadgetron::shared_ptr<ISMRMRDImageData>
945  new_images_container() const = 0;
946  virtual gadgetron::shared_ptr<ISMRMRDImageData>
947  clone(const char* attr, const char* target) = 0;
948  virtual int image_data_type(unsigned int im_num) const
949  {
950  return image_wrap(im_num).type();
951  }
952  virtual size_t num_data_elm() const
953  {
954  return image_wrap(0).num_data_elm();
955  }
956 
957  virtual float norm() const;
959  virtual void sum(void* ptr) const;
960  virtual void max(void* ptr) const;
961  virtual void min(void* ptr) const;
962  virtual void dot(const DataContainer& dc, void* ptr) const;
963  virtual void axpby(
964  const void* ptr_a, const DataContainer& a_x,
965  const void* ptr_b, const DataContainer& a_y);
966  virtual void xapyb(
967  const DataContainer& a_x, const void* ptr_a,
968  const DataContainer& a_y, const void* ptr_b)
969  {
970  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
971  ComplexFloat_ b(*static_cast<const complex_float_t*>(ptr_b));
972  xapyb_(a_x, a, a_y, b);
973  }
974  virtual void xapyb(
975  const DataContainer& a_x, const void* ptr_a,
976  const DataContainer& a_y, const DataContainer& a_b)
977  {
978  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
979  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
980  xapyb_(a_x, a, a_y, b);
981  }
982  //virtual void xapyb(
983  // const DataContainer& a_x, const DataContainer& a_a,
984  // const DataContainer& a_y, const void* ptr_b)
985  //{
986  // SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
987  // ComplexFloat_ b(*(complex_float_t*)ptr_b);
988  // xapyb_(a_x, a, a_y, b);
989  //}
990  virtual void xapyb(
991  const DataContainer& a_x, const DataContainer& a_a,
992  const DataContainer& a_y, const DataContainer& a_b)
993  {
994  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
995  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
996  xapyb_(a_x, a, a_y, b);
997  }
998  virtual void multiply(const DataContainer& x, const DataContainer& y);
999  virtual void divide(const DataContainer& x, const DataContainer& y);
1000  virtual void maximum(const DataContainer& x, const DataContainer& y);
1001  virtual void minimum(const DataContainer& x, const DataContainer& y);
1002  virtual void power(const DataContainer& x, const DataContainer& y);
1003  virtual void multiply(const DataContainer& x, const void* ptr_y);
1004  virtual void add(const DataContainer& x, const void* ptr_y);
1005  virtual void maximum(const DataContainer& x, const void* ptr_y);
1006  virtual void minimum(const DataContainer& x, const void* ptr_y);
1007  virtual void power(const DataContainer& x, const void* ptr_y);
1008  virtual void exp(const DataContainer& x);
1009  virtual void log(const DataContainer& x);
1010  virtual void sqrt(const DataContainer& x);
1011  virtual void sign(const DataContainer& x);
1012  virtual void abs(const DataContainer& x);
1013 
1014  void binary_op(
1015  const DataContainer& a_x, const DataContainer& a_y,
1016  complex_float_t(*f)(complex_float_t, complex_float_t));
1017  void semibinary_op(
1018  const DataContainer& a_x, complex_float_t y,
1019  complex_float_t(*f)(complex_float_t, complex_float_t));
1020  void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t));
1021 
1022  void fill(float s);
1023  void scale(float s);
1024  complex_float_t dot(const DataContainer& a_x)
1025  {
1026  complex_float_t z;
1027  dot(a_x, &z);
1028  return z;
1029  }
1030  void axpby(
1031  complex_float_t a, const DataContainer& a_x,
1032  complex_float_t b, const DataContainer& a_y)
1033  {
1034  axpby(&a, a_x, &b, a_y);
1035  }
1036  void xapyb(
1037  const DataContainer& a_x, complex_float_t a,
1038  const DataContainer& a_y, complex_float_t b)
1039  {
1040  xapyb(a_x, &a, a_y, &b);
1041  }
1042  gadgetron::unique_ptr<ISMRMRDImageData> clone() const
1043  {
1044  return gadgetron::unique_ptr<ISMRMRDImageData>(this->clone_impl());
1045  }
1046 
1047  virtual void sort() = 0;
1048  bool sorted() const { return sorted_; }
1049  void set_sorted(bool sorted) { sorted_ = sorted; }
1050  std::vector<int> index() { return index_; }
1051  const std::vector<int>& index() const { return index_; }
1052  int index(int i) const
1053  {
1054  const std::size_t ni = index_.size();
1055  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni) || static_cast<unsigned>(i) >= number())
1056  THROW("Image number is out of range. You tried to look up an image number that is not inside the container.");
1057  if (ni > 0)
1058  return index_[i];
1059  else
1060  return i;
1061  }
1062  ImageWrap& image_wrap(unsigned int im_num)
1063  {
1064  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
1065  return *sptr_iw;
1066  }
1067  const ImageWrap& image_wrap(unsigned int im_num) const
1068  {
1069  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
1070  sptr_image_wrap(im_num);
1071  return *sptr_iw;
1072  }
1074  void set_meta_data(const AcquisitionsInfo &acqs_info);
1076  const AcquisitionsInfo &get_meta_data() const { return acqs_info_; }
1077 
1078  protected:
1079  bool sorted_=false;
1080  std::vector<int> index_;
1081  AcquisitionsInfo acqs_info_;
1083  virtual ISMRMRDImageData* clone_impl() const = 0;
1084  virtual void conjugate_impl();
1085 
1086  private:
1087  class ComplexFloat_ {
1088  public:
1089  ComplexFloat_(complex_float_t v) : v_(v) {}
1090  unsigned int number() const
1091  {
1092  return 0;
1093  }
1094  complex_float_t image_wrap(unsigned int i)
1095  {
1096  return v_;
1097  }
1098  size_t num_data_elm()
1099  {
1100  return 1;
1101  }
1102  private:
1103  complex_float_t v_;
1104  };
1105 
1106  template<class A, class B>
1107  void xapyb_(const DataContainer& a_x, A& a, const DataContainer& a_y, B& b)
1108  {
1109  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, x, a_x);
1110  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, y, a_y);
1111  unsigned int nx = x.number();
1112  unsigned int na = a.number();
1113  unsigned int ny = y.number();
1114  unsigned int nb = b.number();
1115  //std::cout << nx << ' ' << ny << '\n';
1116  //std::cout << na << ' ' << nb << '\n';
1117  if (nx != ny)
1118  THROW("ImageData sizes mismatch in axpby");
1119  if (na > 0 && na != nx)
1120  THROW("ImageData sizes mismatch in axpby");
1121  if (nb > 0 && nb != nx)
1122  THROW("ImageData sizes mismatch in axpby");
1123  unsigned int n = number();
1124  if (n > 0) {
1125  if (n != nx)
1126  THROW("ImageData sizes mismatch in axpby");
1127  for (unsigned int i = 0; i < nx; i++)
1128  image_wrap(i).xapyb(x.image_wrap(i), a.image_wrap(i),
1129  y.image_wrap(i), b.image_wrap(i));
1130  }
1131  else {
1132  for (unsigned int i = 0; i < nx; i++) {
1133  const ImageWrap& u = x.image_wrap(i);
1134  const ImageWrap& v = y.image_wrap(i);
1135  ImageWrap w(u);
1136  w.xapyb(u, a.image_wrap(i), v, b.image_wrap(i));
1137  append(w);
1138  }
1139  }
1140  this->set_meta_data(x.get_meta_data());
1141  }
1142 
1143  };
1144 
1145  typedef ISMRMRDImageData GadgetronImageData;
1146 
1155  class GadgetronImagesVector : public GadgetronImageData {
1156  public:
1157  typedef ImageData::Iterator BaseIter;
1158  typedef ImageData::Iterator_const BaseIter_const;
1159  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator
1160  ImageWrapIter;
1161  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator
1162  ImageWrapIter_const;
1163  class Iterator : public ImageData::Iterator {
1164  public:
1165  Iterator(ImageWrapIter iw, int n, int i, const ImageWrap::Iterator& it) :
1166  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end())
1167  {}
1168  Iterator(const Iterator& iter) : iw_(iter.iw_), n_(iter.n_), i_(iter.i_),
1169  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1170  {}
1171  Iterator& operator=(const Iterator& iter)
1172  {
1173  iw_ = iter.iw_;
1174  n_ = iter.n_;
1175  i_ = iter.i_;
1176  iter_ = iter.iter_;
1177  end_ = iter.end_;
1178  sptr_iter_ = iter.sptr_iter_;
1179  return *this;
1180  }
1181  virtual bool operator==(const BaseIter& ai) const
1182  {
1183  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1184  return iter_ == i.iter_;
1185  }
1186  virtual bool operator!=(const BaseIter& ai) const
1187  {
1188  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1189  return iter_ != i.iter_;
1190  }
1191  Iterator& operator++()
1192  {
1193  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1194  throw std::out_of_range("cannot advance out-of-range iterator");
1195  ++iter_;
1196  if (iter_ == end_ && i_ < n_ - 1) {
1197  ++i_;
1198  ++iw_;
1199  iter_ = (**iw_).begin();
1200  end_ = (**iw_).end();
1201  }
1202  return *this;
1203  }
1204  Iterator& operator++(int)
1205  {
1206  sptr_iter_.reset(new Iterator(*this));
1207  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1208  throw std::out_of_range("cannot advance out-of-range iterator");
1209  ++iter_;
1210  if (iter_ == end_ && i_ < n_ - 1) {
1211  ++i_;
1212  ++iw_;
1213  iter_ = (**iw_).begin();
1214  end_ = (**iw_).end();
1215  }
1216  return *sptr_iter_;
1217  }
1218  NumRef& operator*()
1219  {
1220  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1221  throw std::out_of_range
1222  ("cannot dereference out-of-range iterator");
1223  return *iter_;
1224  }
1225  private:
1226  //std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator iw_;
1227  ImageWrapIter iw_;
1228  int n_;
1229  int i_;
1230  ImageWrap::Iterator iter_;
1231  ImageWrap::Iterator end_;
1232  gadgetron::shared_ptr<Iterator> sptr_iter_;
1233  };
1234 
1235  class Iterator_const : public ImageData::Iterator_const {
1236  public:
1237  Iterator_const(ImageWrapIter_const iw, int n, int i,
1238  const ImageWrap::Iterator_const& it) :
1239  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end_const())
1240  {}
1241  Iterator_const(const Iterator_const& iter) : iw_(iter.iw_),
1242  n_(iter.n_), i_(iter.i_),
1243  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1244  {}
1245  Iterator_const& operator=(const Iterator_const& iter)
1246  {
1247  iw_ = iter.iw_;
1248  n_ = iter.n_;
1249  i_ = iter.i_;
1250  iter_ = iter.iter_;
1251  end_ = iter.end_;
1252  sptr_iter_ = iter.sptr_iter_;
1253  return *this;
1254  }
1255  bool operator==(const BaseIter_const& ai) const
1256  {
1257  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1258  return iter_ == i.iter_;
1259  }
1260  bool operator!=(const BaseIter_const& ai) const
1261  {
1262  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1263  return iter_ != i.iter_;
1264  }
1265  Iterator_const& operator++()
1266  {
1267  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1268  throw std::out_of_range("cannot advance out-of-range iterator");
1269  ++iter_;
1270  if (iter_ == end_ && i_ < n_ - 1) {
1271  ++i_;
1272  ++iw_;
1273  iter_ = (**iw_).begin_const();
1274  end_ = (**iw_).end_const();
1275  }
1276  return *this;
1277  }
1278  // causes crashes and very inefficient anyway
1279  //Iterator_const& operator++(int)
1280  //{
1281  // sptr_iter_.reset(new Iterator_const(*this));
1282  // if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1283  // throw std::out_of_range("cannot advance out-of-range iterator");
1284  // ++iter_;
1285  // if (iter_ == end_ && i_ < n_ - 1) {
1286  // ++i_;
1287  // ++iw_;
1288  // iter_ = (**iw_).begin_const();
1289  // end_ = (**iw_).end_const();
1290  // }
1291  // return *sptr_iter_;
1292  //}
1293  const NumRef& operator*() const
1294  {
1295  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1296  throw std::out_of_range
1297  ("cannot dereference out-of-range iterator");
1298  ref_.copy(*iter_);
1299  return ref_;
1300  //return *iter_;
1301  }
1302  private:
1303  //std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator iw_;
1304  ImageWrapIter_const iw_;
1305  int n_;
1306  int i_;
1307  ImageWrap::Iterator_const iter_;
1308  ImageWrap::Iterator_const end_;
1309  mutable NumRef ref_;
1310  gadgetron::shared_ptr<Iterator_const> sptr_iter_;
1311  };
1312 
1313  GadgetronImagesVector() : images_()
1314  {}
1315 
1326  GadgetronImagesVector(const MRAcquisitionData& ad, const bool coil_resolved=false);
1327  GadgetronImagesVector(const GadgetronImagesVector& images);
1328  GadgetronImagesVector(GadgetronImagesVector& images, const char* attr,
1329  const char* target);
1330  virtual void empty()
1331  {
1332  images_.clear();
1333  }
1334  virtual unsigned int items() const
1335  {
1336  return (unsigned int)images_.size();
1337  }
1338  virtual unsigned int number() const
1339  {
1340  return (unsigned int)images_.size();
1341  }
1342  virtual void append(int image_data_type, void* ptr_image)
1343  {
1344  images_.push_back(gadgetron::shared_ptr<ImageWrap>
1345  (new ImageWrap(image_data_type, ptr_image)));
1346  }
1347  virtual void append(CFImage& img)
1348  {
1349  void* vptr_img = new CFImage(img);
1350  this->append(7, vptr_img);
1351  }
1352  virtual void append(const ImageWrap& iw)
1353  {
1354  images_.push_back(gadgetron::shared_ptr<ImageWrap>(new ImageWrap(iw)));
1355  }
1356  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw)
1357  {
1358  images_.push_back(sptr_iw);
1359  }
1360  virtual gadgetron::shared_ptr<GadgetronImageData> abs() const;
1361  virtual gadgetron::shared_ptr<GadgetronImageData> real() const;
1362  virtual void clear_data()
1363  {
1364  std::vector<gadgetron::shared_ptr<ImageWrap> > empty_data;
1365  images_.swap(empty_data);
1366  }
1367  virtual void sort();
1368  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
1369  (unsigned int im_num)
1370  {
1371  int i = index(im_num);
1372  return images_.at(i);
1373  }
1374  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
1375  (unsigned int im_num) const
1376  {
1377  int i = index(im_num);
1378  return images_.at(i);
1379  }
1380 /* virtual ImageWrap& image_wrap(unsigned int im_num)
1381  {
1382  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
1383  return *sptr_iw;
1384  }
1385  virtual const ImageWrap& image_wrap(unsigned int im_num) const
1386  {
1387  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
1388  sptr_image_wrap(im_num);
1389  return *sptr_iw;
1390  }
1391 */
1392  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
1393  {
1394  return new ObjectHandle<DataContainer>
1395  (gadgetron::shared_ptr<DataContainer>(new_images_container()));
1396  }
1397  virtual gadgetron::shared_ptr<GadgetronImageData> new_images_container() const
1398  {
1399  gadgetron::shared_ptr<GadgetronImageData> sptr_img
1400  ((GadgetronImageData*)new GadgetronImagesVector());
1401  sptr_img->set_meta_data(get_meta_data());
1402  return sptr_img;
1403  }
1404  virtual gadgetron::shared_ptr<GadgetronImageData>
1405  clone(const char* attr, const char* target)
1406  {
1407  return gadgetron::shared_ptr<GadgetronImageData>
1408  (new GadgetronImagesVector(*this, attr, target));
1409  }
1410 
1411  virtual Iterator& begin()
1412  {
1413  ImageWrapIter iw = images_.begin();
1414  begin_.reset(new Iterator(iw, images_.size(), 0, (**iw).begin()));
1415  return *begin_;
1416  }
1417  virtual Iterator& end()
1418  {
1419  ImageWrapIter iw = images_.begin();
1420  int n = images_.size();
1421  for (int i = 0; i < n - 1; i++)
1422  ++iw;
1423  end_.reset(new Iterator(iw, n, n - 1, (**iw).end()));
1424  return *end_;
1425  }
1426  virtual Iterator_const& begin() const
1427  {
1428  ImageWrapIter_const iw = images_.begin();
1429  begin_const_.reset
1430  (new Iterator_const(iw, images_.size(), 0, (**iw).begin_const()));
1431  return *begin_const_;
1432  }
1433  virtual Iterator_const& end() const
1434  {
1435  ImageWrapIter_const iw = images_.begin();
1436  int n = images_.size();
1437  for (int i = 0; i < n - 1; i++)
1438  ++iw;
1439  end_const_.reset
1440  (new Iterator_const(iw, n, n - 1, (**iw).end_const()));
1441  return *end_const_;
1442  }
1443  virtual void set_image_type(int image_type);
1444  virtual void get_data(complex_float_t* data) const;
1445  virtual void set_data(const complex_float_t* data);
1446  virtual void get_real_data(float* data) const;
1447  virtual void set_real_data(const float* data);
1448 
1450  std::unique_ptr<GadgetronImagesVector> clone() const
1451  {
1452  return std::unique_ptr<GadgetronImagesVector>(this->clone_impl());
1453  }
1454 
1456  void print_header(const unsigned im_num);
1457 
1459  virtual bool is_complex() const;
1460 
1462  virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out);
1463 
1465  virtual void set_up_geom_info();
1466 
1467  private:
1469  virtual GadgetronImagesVector* clone_impl() const
1470  {
1471  return new GadgetronImagesVector(*this);
1472  }
1473 
1474  std::vector<gadgetron::shared_ptr<ImageWrap> > images_;
1475  mutable gadgetron::shared_ptr<Iterator> begin_;
1476  mutable gadgetron::shared_ptr<Iterator> end_;
1477  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
1478  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
1479  };
1480 
1486  class CoilImagesVector : public GadgetronImagesVector
1487  {
1488  public:
1489  CoilImagesVector() : GadgetronImagesVector(){}
1490  void calculate(const MRAcquisitionData& ad);
1491  protected:
1492  std::unique_ptr<MRAcquisitionData> extract_calibration_data(const MRAcquisitionData& ad) const;
1493  gadgetron::shared_ptr<FourierEncoding> sptr_enc_;
1494  };
1495 
1510  class CoilSensitivitiesVector : public GadgetronImagesVector
1511  {
1512 
1513  public:
1514 
1515  CoilSensitivitiesVector() : GadgetronImagesVector(){}
1516  CoilSensitivitiesVector(MRAcquisitionData& ad) :
1517  GadgetronImagesVector(ad, true){}
1518  CoilSensitivitiesVector(const char * file)
1519  {
1520  throw std::runtime_error("This has not been implemented yet.");
1521  }
1522 
1523  void set_csm_smoothness(int s)
1524  {
1525  csm_smoothness_ = s;
1526  }
1527  void set_csm_conv_kernel_size(int w)
1528  {
1529  csm_conv_kernel_halfsize_ = w;
1530  }
1531 
1532  void calculate(CoilImagesVector& iv);
1533  void calculate(const MRAcquisitionData& acq)
1534  {
1535  CoilImagesVector ci;
1536  ci.calculate(acq);
1537  calculate(ci);
1538  }
1539 
1540  CFImage get_csm_as_cfimage(size_t const i) const;
1541  CFImage get_csm_as_cfimage(const KSpaceSubset::TagType tag, const int offset) const;
1542 
1543 
1544  void get_dim(size_t const num_csm, int* dim) const
1545  {
1546  GadgetronImagesVector::get_image_dimensions(num_csm, dim);
1547  }
1548 
1549  void forward(GadgetronImageData& img, const GadgetronImageData& combined_img)const;
1550  void backward(GadgetronImageData& combined_img, const GadgetronImageData& img)const;
1551 
1552  protected:
1553 
1554  void coilchannels_from_combined_image(GadgetronImageData& img, const GadgetronImageData& combined_img) const;
1555  void combine_images_with_coilmaps(GadgetronImageData& combined_img, const GadgetronImageData& img) const;
1556 
1557  void calculate_csm(ISMRMRD::NDArray<complex_float_t>& cm, ISMRMRD::NDArray<float>& img, ISMRMRD::NDArray<complex_float_t>& csm);
1558 
1559  private:
1560  int csm_smoothness_ = 0;
1561  int csm_conv_kernel_halfsize_ = 1;
1562  void smoothen_(int nx, int ny, int nz, int nc, complex_float_t* u, complex_float_t* v,
1563  //int* obj_mask,
1564  int w);
1565  void mask_noise_(int nx, int ny, int nz, float* u, float noise, int* mask);
1566  float max_diff_(int nx, int ny, int nz, int nc, float small_grad, complex_float_t* u, complex_float_t* v);
1567  float max_(int nx, int ny, int nz, float* u);
1568 
1569  };
1570 
1571 void match_img_header_to_acquisition(CFImage& img, const ISMRMRD::Acquisition& acq);
1572 
1573 }
1574 
1575 #endif
ObjectHandle
Definition: DataHandle.h:159
sirf::AcquisitionsInfo
Definition: gadgetron_data_containers.h:70
sirf::AcquisitionsVector
A vector implementation of the abstract MR acquisition data container class.
Definition: gadgetron_data_containers.h:775
sirf::CoilImagesVector
A coil images container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1487
sirf::CoilSensitivitiesVector
A coil sensitivities container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1511
sirf::DataContainer
Abstract data container with numerical operations.
Definition: DataContainer.h:58
sirf::GadgetronImagesVector::Iterator_const
Definition: gadgetron_data_containers.h:1235
sirf::GadgetronImagesVector::Iterator
Definition: gadgetron_data_containers.h:1163
sirf::GadgetronImagesVector
A vector implementation of the abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:1155
sirf::GadgetronImagesVector::is_complex
virtual bool is_complex() const
Is complex?
Definition: gadgetron_data_containers.cpp:2316
sirf::GadgetronImagesVector::set_up_geom_info
virtual void set_up_geom_info()
Populate the geometrical info metadata (from the image's own metadata)
Definition: gadgetron_data_containers.cpp:2417
sirf::GadgetronImagesVector::reorient
virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out)
Reorient image. Requires that dimensions match.
Definition: gadgetron_data_containers.cpp:2324
sirf::GadgetronImagesVector::clone
std::unique_ptr< GadgetronImagesVector > clone() const
Clone and return as unique pointer.
Definition: gadgetron_data_containers.h:1450
sirf::GadgetronImagesVector::print_header
void print_header(const unsigned im_num)
Print header info.
Definition: gadgetron_data_containers.cpp:2279
sirf::ISMRMRDImageData
Abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:855
sirf::ISMRMRDImageData::min
virtual void min(void *ptr) const
calculates the value of this container's element with the smallest real part
Definition: gadgetron_data_containers.cpp:1519
sirf::ISMRMRDImageData::set_meta_data
void set_meta_data(const AcquisitionsInfo &acqs_info)
Set the meta data.
Definition: gadgetron_data_containers.cpp:2055
sirf::ISMRMRDImageData::xapyb
virtual void xapyb(const DataContainer &a_x, const void *ptr_a, const DataContainer &a_y, const DataContainer &a_b)
*this = elementwise sum of x*a and elementwise y*b
Definition: gadgetron_data_containers.h:974
sirf::ISMRMRDImageData::conjugate_impl
virtual void conjugate_impl()
we assume data to be real, complex data containers must override this
Definition: gadgetron_data_containers.cpp:2000
sirf::ISMRMRDImageData::multiply
virtual void multiply(const DataContainer &x, const DataContainer &y)
*this = the elementwise product x*y
Definition: gadgetron_data_containers.cpp:1642
sirf::ISMRMRDImageData::exp
virtual void exp(const DataContainer &x)
*this = the elementwise exp(x)
Definition: gadgetron_data_containers.cpp:1728
sirf::ISMRMRDImageData::xapyb
virtual void xapyb(const DataContainer &a_x, const void *ptr_a, const DataContainer &a_y, const void *ptr_b)
alternative interface to the above
Definition: gadgetron_data_containers.h:966
sirf::ISMRMRDImageData::minimum
virtual void minimum(const DataContainer &x, const DataContainer &y)
*this = the elementwise min(x, y)
Definition: gadgetron_data_containers.cpp:1692
sirf::ISMRMRDImageData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1763
sirf::ISMRMRDImageData::maximum
virtual void maximum(const DataContainer &x, const DataContainer &y)
*this = the elementwise max(x, y)
Definition: gadgetron_data_containers.cpp:1674
sirf::ISMRMRDImageData::sqrt
virtual void sqrt(const DataContainer &x)
*this = the elementwise sqrt(x)
Definition: gadgetron_data_containers.cpp:1742
sirf::ISMRMRDImageData::log
virtual void log(const DataContainer &x)
*this = the elementwise log(x)
Definition: gadgetron_data_containers.cpp:1735
sirf::ISMRMRDImageData::max
virtual void max(void *ptr) const
calculates the value of this container's element with the largest real part
Definition: gadgetron_data_containers.cpp:1503
sirf::ISMRMRDImageData::divide
virtual void divide(const DataContainer &x, const DataContainer &y)
*this = the elementwise ratio x / y
Definition: gadgetron_data_containers.cpp:1666
sirf::ISMRMRDImageData::xapyb
virtual void xapyb(const DataContainer &a_x, const DataContainer &a_a, const DataContainer &a_y, const DataContainer &a_b)
*this = elementwise sum of two elementwise products x*a and y*b
Definition: gadgetron_data_containers.h:990
sirf::ISMRMRDImageData::clone_impl
virtual ISMRMRDImageData * clone_impl() const =0
Clone helper function. Don't use.
sirf::ISMRMRDImageData::dot
virtual void dot(const DataContainer &dc, void *ptr) const
calculates the dot product of this container with another one
Definition: gadgetron_data_containers.cpp:1476
sirf::ISMRMRDImageData::add
virtual void add(const DataContainer &x, const void *ptr_y)
*this = the sum x + y with scalar y
Definition: gadgetron_data_containers.cpp:1658
sirf::ISMRMRDImageData::axpby
virtual void axpby(const void *ptr_a, const DataContainer &a_x, const void *ptr_b, const DataContainer &a_y)
*this = the linear combination of x and y
Definition: gadgetron_data_containers.cpp:1535
sirf::ISMRMRDImageData::power
virtual void power(const DataContainer &x, const DataContainer &y)
*this = the elementwise pow(x, y)
Definition: gadgetron_data_containers.cpp:1710
sirf::ISMRMRDImageData::sign
virtual void sign(const DataContainer &x)
*this = the elementwise sign(x)
Definition: gadgetron_data_containers.cpp:1749
sirf::ISMRMRDImageData::get_meta_data
const AcquisitionsInfo & get_meta_data() const
Get the meta data.
Definition: gadgetron_data_containers.h:1076
sirf::ISMRMRDImageData::sum
virtual void sum(void *ptr) const
below all void* are actually complex_float_t*
Definition: gadgetron_data_containers.cpp:1490
sirf::IgnoreMask
Class enabling ignoring acquisitions with certain ISMRMRD acquisition flags.
Definition: gadgetron_data_containers.h:132
sirf::ImageData::Iterator_const
Definition: ImageData.h:52
sirf::ImageData::Iterator
Definition: ImageData.h:44
sirf::ImageData
Definition: ImageData.h:38
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::KSpaceSubset
Class to keep track of order in k-space.
Definition: gadgetron_data_containers.h:203
sirf::KSpaceSubset::get_tag_from_img
static TagType get_tag_from_img(const CFImage &img)
Function to get k-space dimension tag from an ISMRMRD::Image.
Definition: gadgetron_data_containers.cpp:1423
sirf::KSpaceSubset::get_tag_from_acquisition
static TagType get_tag_from_acquisition(ISMRMRD::Acquisition acq)
Function to get k-space dimension tag from an ISMRMRD::Acquisition.
Definition: gadgetron_data_containers.cpp:1442
sirf::MRAcquisitionData
Abstract MR acquisition data container class.
Definition: gadgetron_data_containers.h:276
sirf::MRAcquisitionData::xapyb
virtual void xapyb(const DataContainer &a_x, const void *ptr_a, const DataContainer &a_y, const void *ptr_b)
alternative interface to the above
Definition: gadgetron_data_containers.h:619
sirf::MRAcquisitionData::get_kspace_order
std::vector< KSpaceSubset::SetType > get_kspace_order() const
Function to get the indices of the acquisitions belonging to different dimensions of k-space.
Definition: gadgetron_data_containers.cpp:1191
sirf::MRAcquisitionData::get_kspace_sorting
std::vector< KSpaceSubset > get_kspace_sorting() const
Function to get the all KSpaceSubset's of the MRAcquisitionData.
Definition: gadgetron_data_containers.h:703
sirf::MRAcquisitionData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1105
sirf::MRAcquisitionData::read
void read(const std::string &filename_ismrmrd_with_ext, int all=0)
Reader for ISMRMRD::Acquisition from ISMRMRD file. In case the ISMRMRD::Dataset constructor throws an...
Definition: gadgetron_data_containers.cpp:99
sirf::MRAcquisitionData::organise_kspace
void organise_kspace()
Function to go through Acquisitions and assign them to their k-space dimension.
Definition: gadgetron_data_containers.cpp:1221
sirf::MRAcquisitionData::set_encoding_limits
void set_encoding_limits(const std::string &name, const std::tuple< unsigned short, unsigned short, unsigned short > min_max_ctr)
Setter for the encoding limits in the header of the acquisition. inputs: name, name of k-space dimens...
Definition: gadgetron_data_containers.h:497
sirf::NumRef
Definition: ANumRef.h:140
sirf::VoxelisedGeometricalInfo
Definition: GeometricalInfo.h:50
gadgetron_image_wrap.h
Specification file for a wrapper class for ISMRMRD::Image.
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
Generated by doxygen 1.9.1