SIRF  3.6.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 int mask = (1 << 18)) :
135  ignore_(mask), max_(8*sizeof(mask)) {}
136  void set(unsigned 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 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 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 int one = 1;
159  return ignore_ & (one << (i - 1));
160  }
161  unsigned long int bits() const
162  {
163  return ignore_;
164  }
165  bool ignored(unsigned 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 int one = 1;
173  unsigned 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 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  // elementwise multiplication
313  // y := x .* y
314  static void multiply
315  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
316  // multiply by scalar
317  // y := x * y
318  static void multiply
319  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
320  // add scalar
321  // y := x + y
322  static void add
323  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
324  // elementwise division
325  // y := x ./ y
326  static void divide
327  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
328  // elementwise maximum
329  // y := std::real(x) > std::real(y) ? x : y
330  static void maximum
331  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
332  static void maximum
333  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
334  // elementwise minimum
335  // y := std::real(x) < std::real(y) ? x : y
336  static void minimum
337  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
338  static void minimum
339  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
340  // y := pow(x, y)
341  static void power
342  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
343  static void power
344  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
345  // y := exp(x)
346  static void exp
347  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
348  // y := log(x)
349  static void log
350  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
351  // y := sqrt(x)
352  static void sqrt
353  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
354  // y := sign(x) (x < 0: -1, x == 0: 0, x > 0: 1)
355  static void sign
356  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
357  // l2 norm of x
358  // y := abs(x)
359  static void abs
360  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
361  static float norm(const ISMRMRD::Acquisition& acq_x);
362 
363  // type and dimension of an ISMRMRD::Acquisition parameter
364  static void ismrmrd_par_info(const char* par, int* output)
365  {
366  // default type: int, dimension: 1
367  output[0] = 0;
368  output[1] = 1;
369 
370  // check parameter type
371  if (sirf::iequals(par, "sample_time_us") ||
372  sirf::iequals(par, "position") ||
373  sirf::iequals(par, "read_dir") ||
374  sirf::iequals(par, "phase_dir") ||
375  sirf::iequals(par, "slice_dir") ||
376  sirf::iequals(par, "patient_table_position") ||
377  sirf::iequals(par, "user_float"))
378  output[0] = 1; // float
379 
380  // check parameter dimension
381  if (sirf::iequals(par, "physiology_time_stamp"))
382  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
383  else if (sirf::iequals(par, "channel_mask"))
384  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
385  else if (sirf::iequals(par, "position") ||
386  sirf::iequals(par, "read_dir") ||
387  sirf::iequals(par, "phase_dir") ||
388  sirf::iequals(par, "slice_dir") ||
389  sirf::iequals(par, "patient_table_position"))
390  output[1] = 3;
391  else if (sirf::iequals(par, "user_int") ||
392  sirf::iequals(par, "idx_user"))
393  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_INTS;
394  else if (sirf::iequals(par, "user_float"))
395  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_FLOATS;
396  }
397  // value of an ISMRMRD::Acquisition int parameter
398  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
399  const char* name, unsigned long long int* v)
400  {
401  if (sirf::iequals(name, "version"))
402  *v = ((unsigned int)acq.version());
403  else if (sirf::iequals(name, "flags"))
404  *v = ((unsigned long long int)acq.flags());
405  else if (sirf::iequals(name, "measurement_uid"))
406  *v = ((unsigned int)acq.measurement_uid());
407  else if (sirf::iequals(name, "scan_counter"))
408  *v = ((unsigned int)acq.scan_counter());
409  else if (sirf::iequals(name, "acquisition_time_stamp"))
410  *v = ((unsigned int)acq.acquisition_time_stamp());
411  else if (sirf::iequals(name, "number_of_samples"))
412  *v = ((unsigned int)acq.number_of_samples());
413  else if (sirf::iequals(name, "available_channels"))
414  *v = ((unsigned int)acq.available_channels());
415  else if (sirf::iequals(name, "active_channels"))
416  *v = ((unsigned int)acq.active_channels());
417  else if (sirf::iequals(name, "discard_pre"))
418  *v = ((unsigned int)acq.discard_pre());
419  else if (sirf::iequals(name, "discard_post"))
420  *v = ((unsigned int)acq.discard_post());
421  else if (sirf::iequals(name, "center_sample"))
422  *v = ((unsigned int)acq.center_sample());
423  else if (sirf::iequals(name, "encoding_space_ref"))
424  *v = ((unsigned int)acq.encoding_space_ref());
425  else if (sirf::iequals(name, "trajectory_dimensions"))
426  *v = ((unsigned int)acq.trajectory_dimensions());
427  else if (sirf::iequals(name, "kspace_encode_step_1"))
428  *v = ((unsigned int)acq.idx().kspace_encode_step_1);
429  else if (sirf::iequals(name, "kspace_encode_step_2"))
430  *v = ((unsigned int)acq.idx().kspace_encode_step_2);
431  else if (sirf::iequals(name, "average"))
432  *v = ((unsigned int)acq.idx().average);
433  else if (sirf::iequals(name, "slice"))
434  *v = ((unsigned int)acq.idx().slice);
435  else if (sirf::iequals(name, "contrast"))
436  *v = ((unsigned int)acq.idx().contrast);
437  else if (sirf::iequals(name, "phase"))
438  *v = ((unsigned int)acq.idx().phase);
439  else if (sirf::iequals(name, "repetition"))
440  *v = ((unsigned int)acq.idx().repetition);
441  else if (sirf::iequals(name, "set"))
442  *v = ((unsigned int)acq.idx().set);
443  else if (sirf::iequals(name, "segment"))
444  *v = ((unsigned int)acq.idx().segment);
445  else if (sirf::iequals(name, "physiology_time_stamp")) {
446  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
447  const uint32_t* pts = acq.physiology_time_stamp();
448  for (int i = 0; i < n; i++)
449  v[i] = (unsigned int)pts[i];
450  }
451  else if (sirf::iequals(name, "channel_mask")) {
452  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
453  const uint64_t* pts = acq.channel_mask();
454  for (int i = 0; i < n; i++)
455  v[i] = (unsigned long long int)pts[i];
456  }
457  }
458  // value of an ISMRMRD::Acquisition float parameter
459  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
460  const char* name, float* v)
461  {
462  if (sirf::iequals(name, "sample_time_us"))
463  *v = acq.sample_time_us();
464  else if (sirf::iequals(name, "position")) {
465  float* u = acq.position();
466  for (int i = 0; i < 3; i++)
467  v[i] = u[i];
468  }
469  else if (sirf::iequals(name, "read_dir")) {
470  float* u = acq.read_dir();
471  for (int i = 0; i < 3; i++)
472  v[i] = u[i];
473  }
474  else if (sirf::iequals(name, "phase_dir")) {
475  float* u = acq.phase_dir();
476  for (int i = 0; i < 3; i++)
477  v[i] = u[i];
478  }
479  else if (sirf::iequals(name, "slice_dir")) {
480  float* u = acq.slice_dir();
481  for (int i = 0; i < 3; i++)
482  v[i] = u[i];
483  }
484  else if (sirf::iequals(name, "patient_table_position")) {
485  float* u = acq.patient_table_position();
486  for (int i = 0; i < 3; i++)
487  v[i] = u[i];
488  }
489  }
496  void set_encoding_limits(const std::string& name,
497  const std::tuple<unsigned short,unsigned short,unsigned short> min_max_ctr)
498  {
499  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
500  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
501 
502  ISMRMRD::Limit limit;
503  limit.minimum = std::get<0>(min_max_ctr);
504  limit.maximum = std::get<1>(min_max_ctr);
505  limit.center = std::get<2>(min_max_ctr);
506 
507  if(sirf::iequals(name, "kspace_encoding_step_1"))
508  enc_limits.kspace_encoding_step_1.get() = limit;
509  else if(sirf::iequals(name, "kspace_encoding_step_2"))
510  enc_limits.kspace_encoding_step_2.get() = limit;
511  else if(sirf::iequals(name, "average"))
512  enc_limits.average.get() = limit;
513  else if(sirf::iequals(name, "slice"))
514  enc_limits.slice.get() = limit;
515  else if(sirf::iequals(name, "contrast"))
516  enc_limits.contrast.get() = limit;
517  else if(sirf::iequals(name, "phase"))
518  enc_limits.phase.get() = limit;
519  else if(sirf::iequals(name, "repetition"))
520  enc_limits.repetition.get() = limit;
521  else if(sirf::iequals(name, "set"))
522  enc_limits.set.get() = limit;
523  else if(sirf::iequals(name, "segment"))
524  enc_limits.segment.get() = limit;
525  else
526  throw std::runtime_error("You passed a name that is not an encoding limit.");
527 
528  hdr.encoding[0].encodingLimits = enc_limits;
529  std::stringstream serialised_hdr;
530  ISMRMRD::serialize(hdr, serialised_hdr);
531  this->set_acquisitions_info(AcquisitionsInfo(serialised_hdr.str()));
532  }
533 
534  std::tuple<unsigned short,unsigned short,unsigned short>
535  get_encoding_limits(const std::string& name) const
536  {
537  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
538  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
539 
540  ISMRMRD::Limit limit;
541 
542  if(sirf::iequals(name, "kspace_encoding_step_1"))
543  limit = enc_limits.kspace_encoding_step_1.get();
544  else if(sirf::iequals(name, "kspace_encoding_step_2"))
545  limit = enc_limits.kspace_encoding_step_2.get();
546  else if(sirf::iequals(name, "average"))
547  limit = enc_limits.average.get();
548  else if(sirf::iequals(name, "slice"))
549  limit = enc_limits.slice.get();
550  else if(sirf::iequals(name, "contrast"))
551  limit = enc_limits.contrast.get();
552  else if(sirf::iequals(name, "phase"))
553  limit = enc_limits.phase.get();
554  else if(sirf::iequals(name, "repetition"))
555  limit = enc_limits.repetition.get();
556  else if(sirf::iequals(name, "set"))
557  limit = enc_limits.set.get();
558  else if(sirf::iequals(name, "segment"))
559  limit = enc_limits.segment.get();
560  else
561  throw std::runtime_error("You passed a name that is not an encoding limit.");
562 
563  return std::make_tuple(limit.minimum, limit.maximum, limit.center);
564  }
565 
566  // abstract methods
567 
568  virtual void empty() = 0;
569  virtual void take_over(MRAcquisitionData&) = 0;
570 
571  // the number of acquisitions in the container
572  virtual unsigned int number() const = 0;
573 
574  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
575  get_acquisition_sptr(unsigned int num) = 0;
576  virtual int get_acquisition(unsigned int,
577  ISMRMRD::Acquisition&) const = 0;
578  virtual void set_acquisition(unsigned int,
579  ISMRMRD::Acquisition&) = 0;
580  virtual void append_acquisition(ISMRMRD::Acquisition& acq) = 0;
581 
582  virtual void copy_acquisitions_info(const MRAcquisitionData& ac) = 0;
583  virtual void copy_acquisitions_data(const MRAcquisitionData& ac) = 0;
584 
585  // 'export' constructors: workaround for creating 'ABC' objects
586  virtual gadgetron::unique_ptr<MRAcquisitionData> new_acquisitions_container() = 0;
587  virtual MRAcquisitionData*
588  same_acquisitions_container(const AcquisitionsInfo& info) const = 0;
589 
590  virtual void set_data(const complex_float_t* z, int all = 1) = 0;
591  virtual void get_data(complex_float_t* z, int all = 1);
592 
593  virtual void set_user_floats(float const * const z, int const idx);
594 
595  virtual bool is_complex() const
596  {
597  return true;
598  }
599 
600  // acquisition data algebra
602  virtual void sum(void* ptr) const;
603  virtual void max(void* ptr) const;
604  virtual void dot(const DataContainer& dc, void* ptr) const;
605  complex_float_t dot(const DataContainer& a_x)
606  {
607  complex_float_t z;
608  dot(a_x, &z);
609  return z;
610  }
611  virtual void axpby(
612  const void* ptr_a, const DataContainer& a_x,
613  const void* ptr_b, const DataContainer& a_y);
614  virtual void xapyb(
615  const DataContainer& a_x, const DataContainer& a_a,
616  const DataContainer& a_y, const DataContainer& a_b);
617  virtual void xapyb(
618  const DataContainer& a_x, const void* ptr_a,
619  const DataContainer& a_y, const void* ptr_b)
620  {
621  axpby(ptr_a, a_x, ptr_b, a_y);
622  }
623  virtual void xapyb(
624  const DataContainer& a_x, const void* ptr_a,
625  const DataContainer& a_y, const DataContainer& a_b);
626  virtual void multiply(const DataContainer& x, const DataContainer& y);
627  virtual void divide(const DataContainer& x, const DataContainer& y);
628  virtual void maximum(const DataContainer& x, const DataContainer& y);
629  virtual void minimum(const DataContainer& x, const DataContainer& y);
630  virtual void power(const DataContainer& x, const DataContainer& y);
631  virtual void multiply(const DataContainer& x, const void* y);
632  virtual void add(const DataContainer& x, const void* ptr_y);
633  virtual void maximum(const DataContainer& x, const void* y);
634  virtual void minimum(const DataContainer& x, const void* y);
635  virtual void power(const DataContainer& x, const void* y);
636  virtual void exp(const DataContainer& x);
637  virtual void log(const DataContainer& x);
638  virtual void sqrt(const DataContainer& x);
639  virtual void sign(const DataContainer& x);
640  virtual void abs(const DataContainer& x);
641  virtual float norm() const;
642 
643  virtual void write(const std::string &filename) const;
644 
645  // regular methods
646  void binary_op(const DataContainer& a_x, const DataContainer& a_y,
647  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
648  void semibinary_op(const DataContainer& a_x, complex_float_t y,
649  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&, complex_float_t));
650  void unary_op(const DataContainer& a_x,
651  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
652 
653  AcquisitionsInfo acquisitions_info() const { return acqs_info_; }
654  void set_acquisitions_info(std::string info) { acqs_info_ = info; }
655  void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;}
656  IgnoreMask ignore_mask() const { return ignore_mask_; }
657  void set_ignore_mask(IgnoreMask ignore_mask) const { ignore_mask_= ignore_mask; }
658 
659  ISMRMRD::TrajectoryType get_trajectory_type() const;
660 
661  void set_trajectory_type(const ISMRMRD::TrajectoryType type);
662 
663  void set_trajectory(const uint16_t traj_dim, float* traj)
664  {
665  ISMRMRD::Acquisition acq;
666  for(int i=0; i<number(); ++i)
667  {
668  get_acquisition(i, acq);
669  const uint16_t num_samples = acq.number_of_samples();
670  const uint16_t num_channels = acq.active_channels();
671  acq.resize(num_samples,num_channels,traj_dim);
672  int const offset = i*traj_dim*num_samples;
673  acq.setTraj(traj + offset);
674  set_acquisition(i, acq);
675  }
676  }
677 
678  gadgetron::unique_ptr<MRAcquisitionData> clone() const
679  {
680  return gadgetron::unique_ptr<MRAcquisitionData>(this->clone_impl());
681  }
682 
683  bool undersampled() const;
684  int get_acquisitions_dimensions(size_t ptr_dim) const;
685  void get_kspace_dimensions(std::vector<size_t>& dims) const;
686  uint16_t get_trajectory_dimensions(void) const;
687 
688  void sort();
689  void sort_by_time();
690  bool sorted() const { return sorted_; }
691  void set_sorted(bool sorted) { sorted_ = sorted; }
692 
694 
698  std::vector<KSpaceSubset::SetType > get_kspace_order() const;
699 
701  std::vector<KSpaceSubset> get_kspace_sorting() const { return this->sorting_; }
702 
704 
710  void organise_kspace();
711 
712  virtual std::vector<int> get_flagged_acquisitions_index
713  (const std::vector<ISMRMRD::ISMRMRD_AcquisitionFlags> flags) const;
714  virtual std::vector<int> get_slice_encoding_index
715  (const unsigned kspace_encode_step_2) const;
716 
717  virtual void get_subset(MRAcquisitionData& subset, const std::vector<int> subset_idx) const;
718  virtual void set_subset(const MRAcquisitionData &subset, const std::vector<int> subset_idx);
719 
720  std::vector<int> index() { return index_; }
721  const std::vector<int>& index() const { return index_; }
722 
723  int index(int i) const
724  {
725  const std::size_t ni = index_.size();
726  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni)
727  || static_cast<unsigned>(i) >= number())
728  THROW("Aquisition number is out of range");
729  if (ni > 0)
730  return index_[i];
731  else
732  return i;
733  }
734 
745  void read(const std::string& filename_ismrmrd_with_ext, int all = 0);
746 
747  protected:
748  bool sorted_ = false;
749  std::vector<int> index_;
750  std::vector<KSpaceSubset> sorting_;
751  AcquisitionsInfo acqs_info_;
752 
753  mutable IgnoreMask ignore_mask_; //= IgnoreMask();
754 
755  // new MRAcquisitionData objects will be created from this template
756  // using same_acquisitions_container()
757  static gadgetron::shared_ptr<MRAcquisitionData> acqs_templ_;
758 
759  virtual MRAcquisitionData* clone_impl() const = 0;
760 
761  private:
762 
763  };
764 
773  class AcquisitionsVector : public MRAcquisitionData {
774  public:
775  AcquisitionsVector(const std::string& filename_with_ext, int all = 0, IgnoreMask ignore_mask = IgnoreMask())
776  {
777  this->set_ignore_mask(ignore_mask);
778  this->read(filename_with_ext, all);
779  }
780 
781  AcquisitionsVector(const AcquisitionsInfo& info = AcquisitionsInfo(), IgnoreMask ignore_mask = IgnoreMask())
782  {
783  this->set_ignore_mask(ignore_mask);
784  acqs_info_ = info;
785  }
786  virtual void empty();
787  virtual void take_over(MRAcquisitionData& ad) {}
788  virtual unsigned int number() const { return (unsigned int)acqs_.size(); }
789  virtual unsigned int items() const { return (unsigned int)acqs_.size(); }
790  virtual void append_acquisition(ISMRMRD::Acquisition& acq)
791  {
792  acqs_.push_back(gadgetron::shared_ptr<ISMRMRD::Acquisition>
793  (new ISMRMRD::Acquisition(acq)));
794  }
795  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
796  get_acquisition_sptr(unsigned int num)
797  {
798  int ind = index(num);
799  return acqs_[ind];
800  }
801  virtual int get_acquisition(unsigned int num,
802  ISMRMRD::Acquisition& acq) const
803  {
804  IgnoreMask ignore_mask = this->ignore_mask();
805  int ind = index(num);
806  acq = *acqs_[ind];
807  if (ignore_mask.ignored(acq.flags()))
808  return 0;
809  return 1;
810  }
811  virtual void set_acquisition(unsigned int num, ISMRMRD::Acquisition& acq)
812  {
813  int ind = index(num);
814  *acqs_[ind] = acq;
815  }
816  virtual void copy_acquisitions_info(const MRAcquisitionData& ac)
817  {
818  acqs_info_ = ac.acquisitions_info();
819  }
820  virtual void copy_acquisitions_data(const MRAcquisitionData& ac);
821  virtual void set_data(const complex_float_t* z, int all = 1);
822 
823  virtual AcquisitionsVector* same_acquisitions_container
824  (const AcquisitionsInfo& info) const
825  {
826  return new AcquisitionsVector(info, ignore_mask_);
827  }
828  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
829  {
830  DataContainer* ptr = new AcquisitionsVector(acqs_info_, ignore_mask_);
831  return new ObjectHandle<DataContainer>
832  (gadgetron::shared_ptr<DataContainer>(ptr));
833  }
834  virtual gadgetron::unique_ptr<MRAcquisitionData>
835  new_acquisitions_container()
836  {
837  return gadgetron::unique_ptr<MRAcquisitionData>
838  (new AcquisitionsVector(acqs_info_, ignore_mask_));
839  }
840 
841  private:
842  std::vector<gadgetron::shared_ptr<ISMRMRD::Acquisition> > acqs_;
843  virtual AcquisitionsVector* clone_impl() const;
844  virtual void conjugate_impl();
845  };
846 
853  class ISMRMRDImageData : public ImageData {
854  public:
855  //ISMRMRDImageData(ISMRMRDImageData& id, const char* attr,
856  //const char* target); //does not build, have to be in the derived class
857 
858  virtual void empty() = 0;
859  virtual unsigned int number() const = 0;
860  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
861  (unsigned int im_num) = 0;
862  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
863  (unsigned int im_num) const = 0;
864 // virtual ImageWrap& image_wrap(unsigned int im_num) = 0;
865 // virtual const ImageWrap& image_wrap(unsigned int im_num) const = 0;
866  virtual void append(int image_data_type, void* ptr_image) = 0;
867  virtual void append(const ImageWrap& iw) = 0;
868  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw) = 0;
869  virtual gadgetron::shared_ptr<ISMRMRDImageData> abs() const = 0;
870  virtual gadgetron::shared_ptr<ISMRMRDImageData> real() const = 0;
871  virtual void clear_data() = 0;
872  virtual void set_image_type(int imtype) = 0;
873  virtual void get_data(complex_float_t* data) const;
874  virtual void set_data(const complex_float_t* data);
875  virtual void get_real_data(float* data) const;
876  virtual void set_real_data(const float* data);
877  virtual int read(std::string filename, std::string variable = "", int iv = -1);
878  virtual void write(const std::string &filename, const std::string &groupname, const bool dicom) const;
879  virtual void write(const std::string &filename) const
880  {
881  size_t size = filename.size();
882  std::string suff = filename.substr(size - 4, 4);
883  if (suff == std::string(".dcm")) {
884  std::string prefix = filename.substr(0, size - 4);
885  this->write(prefix, "", true);
886  }
887  else {
888  auto found = filename.find_last_of("/\\");
889  auto slash_found = found;
890  if (found == std::string::npos)
891  found = filename.find_last_of(".");
892  else
893  found = filename.substr(found + 1).find_last_of(".");
894  if (found == std::string::npos)
895  this->write(filename + ".h5", "", false);
896  else {
897  std::string ext = filename.substr(slash_found + found + 1);
898  if (ext == std::string(".h5"))
899  this->write(filename, "", false);
900  else
901  std::cerr << "WARNING: writing ISMRMRD images to "
902  << ext << "-files not implemented, "
903  << "please convert to Nifti images\n";
904  }
905  }
906  }
907  virtual Dimensions dimensions() const
908  {
909  Dimensions dim;
910  const ImageWrap& iw = image_wrap(0);
911  int d[5];
912  iw.get_dim(d);
913  dim["x"] = d[0];
914  dim["y"] = d[1];
915  dim["z"] = d[2];
916  dim["c"] = d[3];
917  dim["n"] = number();
918  return dim;
919  }
920  virtual void get_image_dimensions(unsigned int im_num, int* dim) const
921  {
922  if (im_num >= number())
923  dim[0] = dim[1] = dim[2] = dim[3] = 0;
924  const ImageWrap& iw = image_wrap(im_num);
925  iw.get_dim(dim);
926  }
927  bool check_dimension_consistency() const
928  {
929  size_t const num_dims = 4;
930  std::vector<int> first_img_dims(num_dims), temp_img_dims(num_dims);
931 
932  this->get_image_dimensions(0, &first_img_dims[0]);
933 
934  bool dims_match = true;
935  for(int i=1; i<number(); ++i)
936  {
937  this->get_image_dimensions(0, &temp_img_dims[0]);
938  dims_match *= (first_img_dims == temp_img_dims);
939  }
940  return dims_match;
941  }
942  virtual gadgetron::shared_ptr<ISMRMRDImageData>
943  new_images_container() const = 0;
944  virtual gadgetron::shared_ptr<ISMRMRDImageData>
945  clone(const char* attr, const char* target) = 0;
946  virtual int image_data_type(unsigned int im_num) const
947  {
948  return image_wrap(im_num).type();
949  }
950  virtual size_t num_data_elm() const
951  {
952  return image_wrap(0).num_data_elm();
953  }
954 
955  virtual float norm() const;
957  virtual void sum(void* ptr) const;
958  virtual void max(void* ptr) const;
959  virtual void dot(const DataContainer& dc, void* ptr) const;
960  virtual void axpby(
961  const void* ptr_a, const DataContainer& a_x,
962  const void* ptr_b, const DataContainer& a_y);
963  virtual void xapyb(
964  const DataContainer& a_x, const void* ptr_a,
965  const DataContainer& a_y, const void* ptr_b)
966  {
967  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
968  ComplexFloat_ b(*static_cast<const complex_float_t*>(ptr_b));
969  xapyb_(a_x, a, a_y, b);
970  }
971  virtual void xapyb(
972  const DataContainer& a_x, const void* ptr_a,
973  const DataContainer& a_y, const DataContainer& a_b)
974  {
975  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
976  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
977  xapyb_(a_x, a, a_y, b);
978  }
979  //virtual void xapyb(
980  // const DataContainer& a_x, const DataContainer& a_a,
981  // const DataContainer& a_y, const void* ptr_b)
982  //{
983  // SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
984  // ComplexFloat_ b(*(complex_float_t*)ptr_b);
985  // xapyb_(a_x, a, a_y, b);
986  //}
987  virtual void xapyb(
988  const DataContainer& a_x, const DataContainer& a_a,
989  const DataContainer& a_y, const DataContainer& a_b)
990  {
991  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
992  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
993  xapyb_(a_x, a, a_y, b);
994  }
995  virtual void multiply(const DataContainer& x, const DataContainer& y);
996  virtual void divide(const DataContainer& x, const DataContainer& y);
997  virtual void maximum(const DataContainer& x, const DataContainer& y);
998  virtual void minimum(const DataContainer& x, const DataContainer& y);
999  virtual void power(const DataContainer& x, const DataContainer& y);
1000  virtual void multiply(const DataContainer& x, const void* ptr_y);
1001  virtual void add(const DataContainer& x, const void* ptr_y);
1002  virtual void maximum(const DataContainer& x, const void* ptr_y);
1003  virtual void minimum(const DataContainer& x, const void* ptr_y);
1004  virtual void power(const DataContainer& x, const void* ptr_y);
1005  virtual void exp(const DataContainer& x);
1006  virtual void log(const DataContainer& x);
1007  virtual void sqrt(const DataContainer& x);
1008  virtual void sign(const DataContainer& x);
1009  virtual void abs(const DataContainer& x);
1010 
1011  void binary_op(
1012  const DataContainer& a_x, const DataContainer& a_y,
1013  complex_float_t(*f)(complex_float_t, complex_float_t));
1014  void semibinary_op(
1015  const DataContainer& a_x, complex_float_t y,
1016  complex_float_t(*f)(complex_float_t, complex_float_t));
1017  void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t));
1018 
1019  void fill(float s);
1020  void scale(float s);
1021  complex_float_t dot(const DataContainer& a_x)
1022  {
1023  complex_float_t z;
1024  dot(a_x, &z);
1025  return z;
1026  }
1027  void axpby(
1028  complex_float_t a, const DataContainer& a_x,
1029  complex_float_t b, const DataContainer& a_y)
1030  {
1031  axpby(&a, a_x, &b, a_y);
1032  }
1033  void xapyb(
1034  const DataContainer& a_x, complex_float_t a,
1035  const DataContainer& a_y, complex_float_t b)
1036  {
1037  xapyb(a_x, &a, a_y, &b);
1038  }
1039  gadgetron::unique_ptr<ISMRMRDImageData> clone() const
1040  {
1041  return gadgetron::unique_ptr<ISMRMRDImageData>(this->clone_impl());
1042  }
1043 
1044  virtual void sort() = 0;
1045  bool sorted() const { return sorted_; }
1046  void set_sorted(bool sorted) { sorted_ = sorted; }
1047  std::vector<int> index() { return index_; }
1048  const std::vector<int>& index() const { return index_; }
1049  int index(int i) const
1050  {
1051  const std::size_t ni = index_.size();
1052  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni) || static_cast<unsigned>(i) >= number())
1053  THROW("Image number is out of range. You tried to look up an image number that is not inside the container.");
1054  if (ni > 0)
1055  return index_[i];
1056  else
1057  return i;
1058  }
1059  ImageWrap& image_wrap(unsigned int im_num)
1060  {
1061  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
1062  return *sptr_iw;
1063  }
1064  const ImageWrap& image_wrap(unsigned int im_num) const
1065  {
1066  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
1067  sptr_image_wrap(im_num);
1068  return *sptr_iw;
1069  }
1071  void set_meta_data(const AcquisitionsInfo &acqs_info);
1073  const AcquisitionsInfo &get_meta_data() const { return acqs_info_; }
1074 
1075  protected:
1076  bool sorted_=false;
1077  std::vector<int> index_;
1078  AcquisitionsInfo acqs_info_;
1080  virtual ISMRMRDImageData* clone_impl() const = 0;
1081  virtual void conjugate_impl();
1082 
1083  private:
1084  class ComplexFloat_ {
1085  public:
1086  ComplexFloat_(complex_float_t v) : v_(v) {}
1087  unsigned int number() const
1088  {
1089  return 0;
1090  }
1091  complex_float_t image_wrap(unsigned int i)
1092  {
1093  return v_;
1094  }
1095  size_t num_data_elm()
1096  {
1097  return 1;
1098  }
1099  private:
1100  complex_float_t v_;
1101  };
1102 
1103  template<class A, class B>
1104  void xapyb_(const DataContainer& a_x, A& a, const DataContainer& a_y, B& b)
1105  {
1106  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, x, a_x);
1107  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, y, a_y);
1108  unsigned int nx = x.number();
1109  unsigned int na = a.number();
1110  unsigned int ny = y.number();
1111  unsigned int nb = b.number();
1112  //std::cout << nx << ' ' << ny << '\n';
1113  //std::cout << na << ' ' << nb << '\n';
1114  if (nx != ny)
1115  THROW("ImageData sizes mismatch in axpby");
1116  if (na > 0 && na != nx)
1117  THROW("ImageData sizes mismatch in axpby");
1118  if (nb > 0 && nb != nx)
1119  THROW("ImageData sizes mismatch in axpby");
1120  unsigned int n = number();
1121  if (n > 0) {
1122  if (n != nx)
1123  THROW("ImageData sizes mismatch in axpby");
1124  for (unsigned int i = 0; i < nx; i++)
1125  image_wrap(i).xapyb(x.image_wrap(i), a.image_wrap(i),
1126  y.image_wrap(i), b.image_wrap(i));
1127  }
1128  else {
1129  for (unsigned int i = 0; i < nx; i++) {
1130  const ImageWrap& u = x.image_wrap(i);
1131  const ImageWrap& v = y.image_wrap(i);
1132  ImageWrap w(u);
1133  w.xapyb(u, a.image_wrap(i), v, b.image_wrap(i));
1134  append(w);
1135  }
1136  }
1137  this->set_meta_data(x.get_meta_data());
1138  }
1139 
1140  };
1141 
1142  typedef ISMRMRDImageData GadgetronImageData;
1143 
1152  class GadgetronImagesVector : public GadgetronImageData {
1153  public:
1154  typedef ImageData::Iterator BaseIter;
1155  typedef ImageData::Iterator_const BaseIter_const;
1156  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator
1157  ImageWrapIter;
1158  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator
1159  ImageWrapIter_const;
1160  class Iterator : public ImageData::Iterator {
1161  public:
1162  Iterator(ImageWrapIter iw, int n, int i, const ImageWrap::Iterator& it) :
1163  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end())
1164  {}
1165  Iterator(const Iterator& iter) : iw_(iter.iw_), n_(iter.n_), i_(iter.i_),
1166  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1167  {}
1168  Iterator& operator=(const Iterator& iter)
1169  {
1170  iw_ = iter.iw_;
1171  n_ = iter.n_;
1172  i_ = iter.i_;
1173  iter_ = iter.iter_;
1174  end_ = iter.end_;
1175  sptr_iter_ = iter.sptr_iter_;
1176  return *this;
1177  }
1178  virtual bool operator==(const BaseIter& ai) const
1179  {
1180  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1181  return iter_ == i.iter_;
1182  }
1183  virtual bool operator!=(const BaseIter& ai) const
1184  {
1185  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1186  return iter_ != i.iter_;
1187  }
1188  Iterator& operator++()
1189  {
1190  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1191  throw std::out_of_range("cannot advance out-of-range iterator");
1192  ++iter_;
1193  if (iter_ == end_ && i_ < n_ - 1) {
1194  ++i_;
1195  ++iw_;
1196  iter_ = (**iw_).begin();
1197  end_ = (**iw_).end();
1198  }
1199  return *this;
1200  }
1201  Iterator& operator++(int)
1202  {
1203  sptr_iter_.reset(new Iterator(*this));
1204  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1205  throw std::out_of_range("cannot advance out-of-range iterator");
1206  ++iter_;
1207  if (iter_ == end_ && i_ < n_ - 1) {
1208  ++i_;
1209  ++iw_;
1210  iter_ = (**iw_).begin();
1211  end_ = (**iw_).end();
1212  }
1213  return *sptr_iter_;
1214  }
1215  NumRef& operator*()
1216  {
1217  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1218  throw std::out_of_range
1219  ("cannot dereference out-of-range iterator");
1220  return *iter_;
1221  }
1222  private:
1223  //std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator iw_;
1224  ImageWrapIter iw_;
1225  int n_;
1226  int i_;
1227  ImageWrap::Iterator iter_;
1228  ImageWrap::Iterator end_;
1229  gadgetron::shared_ptr<Iterator> sptr_iter_;
1230  };
1231 
1232  class Iterator_const : public ImageData::Iterator_const {
1233  public:
1234  Iterator_const(ImageWrapIter_const iw, int n, int i,
1235  const ImageWrap::Iterator_const& it) :
1236  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end_const())
1237  {}
1238  Iterator_const(const Iterator_const& iter) : iw_(iter.iw_),
1239  n_(iter.n_), i_(iter.i_),
1240  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1241  {}
1242  Iterator_const& operator=(const Iterator_const& iter)
1243  {
1244  iw_ = iter.iw_;
1245  n_ = iter.n_;
1246  i_ = iter.i_;
1247  iter_ = iter.iter_;
1248  end_ = iter.end_;
1249  sptr_iter_ = iter.sptr_iter_;
1250  return *this;
1251  }
1252  bool operator==(const BaseIter_const& ai) const
1253  {
1254  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1255  return iter_ == i.iter_;
1256  }
1257  bool operator!=(const BaseIter_const& ai) const
1258  {
1259  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1260  return iter_ != i.iter_;
1261  }
1262  Iterator_const& operator++()
1263  {
1264  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1265  throw std::out_of_range("cannot advance out-of-range iterator");
1266  ++iter_;
1267  if (iter_ == end_ && i_ < n_ - 1) {
1268  ++i_;
1269  ++iw_;
1270  iter_ = (**iw_).begin_const();
1271  end_ = (**iw_).end_const();
1272  }
1273  return *this;
1274  }
1275  // causes crashes and very inefficient anyway
1276  //Iterator_const& operator++(int)
1277  //{
1278  // sptr_iter_.reset(new Iterator_const(*this));
1279  // if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1280  // throw std::out_of_range("cannot advance out-of-range iterator");
1281  // ++iter_;
1282  // if (iter_ == end_ && i_ < n_ - 1) {
1283  // ++i_;
1284  // ++iw_;
1285  // iter_ = (**iw_).begin_const();
1286  // end_ = (**iw_).end_const();
1287  // }
1288  // return *sptr_iter_;
1289  //}
1290  const NumRef& operator*() const
1291  {
1292  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1293  throw std::out_of_range
1294  ("cannot dereference out-of-range iterator");
1295  ref_.copy(*iter_);
1296  return ref_;
1297  //return *iter_;
1298  }
1299  private:
1300  //std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator iw_;
1301  ImageWrapIter_const iw_;
1302  int n_;
1303  int i_;
1304  ImageWrap::Iterator_const iter_;
1305  ImageWrap::Iterator_const end_;
1306  mutable NumRef ref_;
1307  gadgetron::shared_ptr<Iterator_const> sptr_iter_;
1308  };
1309 
1310  GadgetronImagesVector() : images_()
1311  {}
1312 
1323  GadgetronImagesVector(const MRAcquisitionData& ad, const bool coil_resolved=false);
1324  GadgetronImagesVector(const GadgetronImagesVector& images);
1325  GadgetronImagesVector(GadgetronImagesVector& images, const char* attr,
1326  const char* target);
1327  virtual void empty()
1328  {
1329  images_.clear();
1330  }
1331  virtual unsigned int items() const
1332  {
1333  return (unsigned int)images_.size();
1334  }
1335  virtual unsigned int number() const
1336  {
1337  return (unsigned int)images_.size();
1338  }
1339  virtual void append(int image_data_type, void* ptr_image)
1340  {
1341  images_.push_back(gadgetron::shared_ptr<ImageWrap>
1342  (new ImageWrap(image_data_type, ptr_image)));
1343  }
1344  virtual void append(CFImage& img)
1345  {
1346  void* vptr_img = new CFImage(img);
1347  this->append(7, vptr_img);
1348  }
1349  virtual void append(const ImageWrap& iw)
1350  {
1351  images_.push_back(gadgetron::shared_ptr<ImageWrap>(new ImageWrap(iw)));
1352  }
1353  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw)
1354  {
1355  images_.push_back(sptr_iw);
1356  }
1357  virtual gadgetron::shared_ptr<GadgetronImageData> abs() const;
1358  virtual gadgetron::shared_ptr<GadgetronImageData> real() const;
1359  virtual void clear_data()
1360  {
1361  std::vector<gadgetron::shared_ptr<ImageWrap> > empty_data;
1362  images_.swap(empty_data);
1363  }
1364  virtual void sort();
1365  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
1366  (unsigned int im_num)
1367  {
1368  int i = index(im_num);
1369  return images_.at(i);
1370  }
1371  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
1372  (unsigned int im_num) const
1373  {
1374  int i = index(im_num);
1375  return images_.at(i);
1376  }
1377 /* virtual ImageWrap& image_wrap(unsigned int im_num)
1378  {
1379  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
1380  return *sptr_iw;
1381  }
1382  virtual const ImageWrap& image_wrap(unsigned int im_num) const
1383  {
1384  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
1385  sptr_image_wrap(im_num);
1386  return *sptr_iw;
1387  }
1388 */
1389  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
1390  {
1391  return new ObjectHandle<DataContainer>
1392  (gadgetron::shared_ptr<DataContainer>(new_images_container()));
1393  }
1394  virtual gadgetron::shared_ptr<GadgetronImageData> new_images_container() const
1395  {
1396  gadgetron::shared_ptr<GadgetronImageData> sptr_img
1397  ((GadgetronImageData*)new GadgetronImagesVector());
1398  sptr_img->set_meta_data(get_meta_data());
1399  return sptr_img;
1400  }
1401  virtual gadgetron::shared_ptr<GadgetronImageData>
1402  clone(const char* attr, const char* target)
1403  {
1404  return gadgetron::shared_ptr<GadgetronImageData>
1405  (new GadgetronImagesVector(*this, attr, target));
1406  }
1407 
1408  virtual Iterator& begin()
1409  {
1410  ImageWrapIter iw = images_.begin();
1411  begin_.reset(new Iterator(iw, images_.size(), 0, (**iw).begin()));
1412  return *begin_;
1413  }
1414  virtual Iterator& end()
1415  {
1416  ImageWrapIter iw = images_.begin();
1417  int n = images_.size();
1418  for (int i = 0; i < n - 1; i++)
1419  ++iw;
1420  end_.reset(new Iterator(iw, n, n - 1, (**iw).end()));
1421  return *end_;
1422  }
1423  virtual Iterator_const& begin() const
1424  {
1425  ImageWrapIter_const iw = images_.begin();
1426  begin_const_.reset
1427  (new Iterator_const(iw, images_.size(), 0, (**iw).begin_const()));
1428  return *begin_const_;
1429  }
1430  virtual Iterator_const& end() const
1431  {
1432  ImageWrapIter_const iw = images_.begin();
1433  int n = images_.size();
1434  for (int i = 0; i < n - 1; i++)
1435  ++iw;
1436  end_const_.reset
1437  (new Iterator_const(iw, n, n - 1, (**iw).end_const()));
1438  return *end_const_;
1439  }
1440  virtual void set_image_type(int image_type);
1441  virtual void get_data(complex_float_t* data) const;
1442  virtual void set_data(const complex_float_t* data);
1443  virtual void get_real_data(float* data) const;
1444  virtual void set_real_data(const float* data);
1445 
1447  std::unique_ptr<GadgetronImagesVector> clone() const
1448  {
1449  return std::unique_ptr<GadgetronImagesVector>(this->clone_impl());
1450  }
1451 
1453  void print_header(const unsigned im_num);
1454 
1456  virtual bool is_complex() const;
1457 
1459  virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out);
1460 
1462  virtual void set_up_geom_info();
1463 
1464  private:
1466  virtual GadgetronImagesVector* clone_impl() const
1467  {
1468  return new GadgetronImagesVector(*this);
1469  }
1470 
1471  std::vector<gadgetron::shared_ptr<ImageWrap> > images_;
1472  mutable gadgetron::shared_ptr<Iterator> begin_;
1473  mutable gadgetron::shared_ptr<Iterator> end_;
1474  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
1475  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
1476  };
1477 
1483  class CoilImagesVector : public GadgetronImagesVector
1484  {
1485  public:
1486  CoilImagesVector() : GadgetronImagesVector(){}
1487  void calculate(const MRAcquisitionData& ad);
1488  protected:
1489  std::unique_ptr<MRAcquisitionData> extract_calibration_data(const MRAcquisitionData& ad) const;
1490  gadgetron::shared_ptr<FourierEncoding> sptr_enc_;
1491  };
1492 
1507  class CoilSensitivitiesVector : public GadgetronImagesVector
1508  {
1509 
1510  public:
1511 
1512  CoilSensitivitiesVector() : GadgetronImagesVector(){}
1513  CoilSensitivitiesVector(MRAcquisitionData& ad) :
1514  GadgetronImagesVector(ad, true){}
1515  CoilSensitivitiesVector(const char * file)
1516  {
1517  throw std::runtime_error("This has not been implemented yet.");
1518  }
1519 
1520  void set_csm_smoothness(int s)
1521  {
1522  csm_smoothness_ = s;
1523  }
1524  void set_csm_conv_kernel_size(int w)
1525  {
1526  csm_conv_kernel_halfsize_ = w;
1527  }
1528 
1529  void calculate(CoilImagesVector& iv);
1530  void calculate(const MRAcquisitionData& acq)
1531  {
1532  CoilImagesVector ci;
1533  ci.calculate(acq);
1534  calculate(ci);
1535  }
1536 
1537  CFImage get_csm_as_cfimage(size_t const i) const;
1538  CFImage get_csm_as_cfimage(const KSpaceSubset::TagType tag, const int offset) const;
1539 
1540 
1541  void get_dim(size_t const num_csm, int* dim) const
1542  {
1543  GadgetronImagesVector::get_image_dimensions(num_csm, dim);
1544  }
1545 
1546  void forward(GadgetronImageData& img, const GadgetronImageData& combined_img)const;
1547  void backward(GadgetronImageData& combined_img, const GadgetronImageData& img)const;
1548 
1549  protected:
1550 
1551  void coilchannels_from_combined_image(GadgetronImageData& img, const GadgetronImageData& combined_img) const;
1552  void combine_images_with_coilmaps(GadgetronImageData& combined_img, const GadgetronImageData& img) const;
1553 
1554  void calculate_csm(ISMRMRD::NDArray<complex_float_t>& cm, ISMRMRD::NDArray<float>& img, ISMRMRD::NDArray<complex_float_t>& csm);
1555 
1556  private:
1557  int csm_smoothness_ = 0;
1558  int csm_conv_kernel_halfsize_ = 1;
1559  void smoothen_(int nx, int ny, int nz, int nc, complex_float_t* u, complex_float_t* v,
1560  //int* obj_mask,
1561  int w);
1562  void mask_noise_(int nx, int ny, int nz, float* u, float noise, int* mask);
1563  float max_diff_(int nx, int ny, int nz, int nc, float small_grad, complex_float_t* u, complex_float_t* v);
1564  float max_(int nx, int ny, int nz, float* u);
1565 
1566  };
1567 
1568 void match_img_header_to_acquisition(CFImage& img, const ISMRMRD::Acquisition& acq);
1569 
1570 }
1571 
1572 #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:773
sirf::CoilImagesVector
A coil images container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1484
sirf::CoilSensitivitiesVector
A coil sensitivities container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1508
sirf::DataContainer
Definition: DataContainer.h:42
sirf::GadgetronImagesVector::Iterator_const
Definition: gadgetron_data_containers.h:1232
sirf::GadgetronImagesVector::Iterator
Definition: gadgetron_data_containers.h:1160
sirf::GadgetronImagesVector
A vector implementation of the abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:1152
sirf::GadgetronImagesVector::is_complex
virtual bool is_complex() const
Is complex?
Definition: gadgetron_data_containers.cpp:2251
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:2352
sirf::GadgetronImagesVector::reorient
virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out)
Reorient image. Requires that dimensions match.
Definition: gadgetron_data_containers.cpp:2259
sirf::GadgetronImagesVector::clone
std::unique_ptr< GadgetronImagesVector > clone() const
Clone and return as unique pointer.
Definition: gadgetron_data_containers.h:1447
sirf::GadgetronImagesVector::print_header
void print_header(const unsigned im_num)
Print header info.
Definition: gadgetron_data_containers.cpp:2214
sirf::ISMRMRDImageData
Abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:853
sirf::ISMRMRDImageData::set_meta_data
void set_meta_data(const AcquisitionsInfo &acqs_info)
Set the meta data.
Definition: gadgetron_data_containers.cpp:1990
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:971
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:1935
sirf::ISMRMRDImageData::multiply
virtual void multiply(const DataContainer &x, const DataContainer &y)
*this = the elementwise product x*y
Definition: gadgetron_data_containers.cpp:1577
sirf::ISMRMRDImageData::exp
virtual void exp(const DataContainer &x)
*this = the elementwise exp(x)
Definition: gadgetron_data_containers.cpp:1663
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:963
sirf::ISMRMRDImageData::minimum
virtual void minimum(const DataContainer &x, const DataContainer &y)
*this = the elementwise min(x, y)
Definition: gadgetron_data_containers.cpp:1627
sirf::ISMRMRDImageData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1698
sirf::ISMRMRDImageData::maximum
virtual void maximum(const DataContainer &x, const DataContainer &y)
*this = the elementwise max(x, y)
Definition: gadgetron_data_containers.cpp:1609
sirf::ISMRMRDImageData::sqrt
virtual void sqrt(const DataContainer &x)
*this = the elementwise sqrt(x)
Definition: gadgetron_data_containers.cpp:1677
sirf::ISMRMRDImageData::log
virtual void log(const DataContainer &x)
*this = the elementwise log(x)
Definition: gadgetron_data_containers.cpp:1670
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:1454
sirf::ISMRMRDImageData::divide
virtual void divide(const DataContainer &x, const DataContainer &y)
*this = the elementwise ratio x / y
Definition: gadgetron_data_containers.cpp:1601
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:987
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:1427
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:1593
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:1470
sirf::ISMRMRDImageData::power
virtual void power(const DataContainer &x, const DataContainer &y)
*this = the elementwise pow(x, y)
Definition: gadgetron_data_containers.cpp:1645
sirf::ISMRMRDImageData::sign
virtual void sign(const DataContainer &x)
*this = the elementwise sign(x)
Definition: gadgetron_data_containers.cpp:1684
sirf::ISMRMRDImageData::get_meta_data
const AcquisitionsInfo & get_meta_data() const
Get the meta data.
Definition: gadgetron_data_containers.h:1073
sirf::ISMRMRDImageData::sum
virtual void sum(void *ptr) const
below all void* are actually complex_float_t*
Definition: gadgetron_data_containers.cpp:1441
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:1374
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:1393
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:617
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:1142
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:701
sirf::MRAcquisitionData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1056
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:1172
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:496
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 data container.
Definition: GeometricalInfo.cpp:141
sirf::iequals
bool iequals(const std::string &a, const std::string &b)
Case insensitive string comparison, replaces boost::iequals.
Definition: iequals.cpp:7
Generated by doxygen 1.9.1