SIRF  3.5.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 
63 #define TO_BE_IGNORED(acq) \
64  (!(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION) && \
65  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING) && \
66  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_MEASUREMENT) && \
67  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_REVERSE) && \
68  (acq).flags() >= (1 << (ISMRMRD::ISMRMRD_ACQ_IS_NOISE_MEASUREMENT - 1)))
69 
76 namespace sirf {
77 
78  class FourierEncoding;
79  class CartesianFourierEncoding;
80 #if GADGETRON_TOOLBOXES_AVAILABLE
81  class RPEFourierEncoding;
82 #endif
83 
84  class AcquisitionsInfo {
85  public:
86  AcquisitionsInfo(std::string data = "") : data_(data)
87  {
88  if (data.empty())
89  have_header_ = false;
90  else {
91  deserialize();
92  have_header_ = true;
93  }
94  }
95  AcquisitionsInfo& operator=(std::string data)
96  {
97  data_ = data;
98  if (data.empty())
99  have_header_ = false;
100  else {
101  deserialize();
102  have_header_ = true;
103  }
104 
105  return *this;
106  }
107  const char* c_str() const { return data_.c_str(); }
108  operator std::string&() { return data_; }
109  operator const std::string&() const { return data_; }
110  bool empty() const { return data_.empty(); }
111  const ISMRMRD::IsmrmrdHeader& get_IsmrmrdHeader() const
112  {
113  if (!have_header_)
114  deserialize();
115  return header_;
116  }
117 
118  private:
119  void deserialize() const
120  {
121  if (!this->empty())
122  { header_ = ISMRMRD::IsmrmrdHeader();
123  ISMRMRD::deserialize(data_.c_str(), header_);
124  }
125  have_header_ = true;
126  }
127  std::string data_;
128  mutable ISMRMRD::IsmrmrdHeader header_;
129  mutable bool have_header_;
130  };
131 
143  class KSpaceSubset
144  {
145  static int const num_kspace_dims_ = 7 + ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_INTS;
146 
147  public:
148 
149  typedef std::array<int, num_kspace_dims_> TagType;
150  typedef std::vector<int> SetType;
151 
152  KSpaceSubset(){
153  for(int i=0; i<num_kspace_dims_; ++i)
154  this->tag_[i] = -1;
155  }
156 
157  KSpaceSubset(TagType tag){
158  this->tag_ = tag;
159  this->idx_set_ = {};
160  }
161 
162  KSpaceSubset(TagType tag, SetType idx_set){
163  this->tag_ = tag;
164  this->idx_set_ = idx_set;
165  }
166 
167  TagType get_tag(void) const {return tag_;}
168  SetType get_idx_set(void) const {return idx_set_;}
169  void add_idx_to_set(size_t const idx){this->idx_set_.push_back(idx);}
170 
171  bool is_first_set() const {
172  bool is_first= (tag_[0] == 0);
173  if(is_first)
174  {
175  for(int dim=2; dim<num_kspace_dims_; ++dim)
176  is_first *= (tag_[dim] == 0);
177  }
178  return is_first;
179  }
180 
181  static void print_tag(const TagType& tag);
182  static void print_acquisition_tag(ISMRMRD::Acquisition acq);
183 
185 
188  static TagType get_tag_from_acquisition(ISMRMRD::Acquisition acq);
189 
191 
194  static TagType get_tag_from_img(const CFImage& img);
195 
196  private:
198 
202  TagType tag_;
203 
205 
209  SetType idx_set_;
210  };
211 
217  class MRAcquisitionData : public DataContainer {
218  public:
219  // static methods
220 
221  // ISMRMRD acquisitions algebra: acquisitions viewed as vectors of
222  // acquisition data
223  static void binary_op
224  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y,
225  complex_float_t (*f)(complex_float_t, complex_float_t));
226  static void semibinary_op
227  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y,
228  complex_float_t(*f)(complex_float_t, complex_float_t));
229  static void unary_op
230  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y,
231  complex_float_t(*f)(complex_float_t));
232  // y := a x + b y
233  static void axpby
234  (complex_float_t a, const ISMRMRD::Acquisition& acq_x,
235  complex_float_t b, ISMRMRD::Acquisition& acq_y);
236  static void xapyb
237  (const ISMRMRD::Acquisition& acq_x, complex_float_t a,
238  ISMRMRD::Acquisition& acq_y, complex_float_t b);
239  static void xapyb
240  (const ISMRMRD::Acquisition& acq_x, complex_float_t a,
241  ISMRMRD::Acquisition& acq_y, const ISMRMRD::Acquisition& acq_b);
242  static void xapyb
243  (const ISMRMRD::Acquisition& acq_x, const ISMRMRD::Acquisition& acq_a,
244  ISMRMRD::Acquisition& acq_y, const ISMRMRD::Acquisition& acq_b);
245 
246  // the inner (l2) product of x and y
247  static complex_float_t dot
248  (const ISMRMRD::Acquisition& acq_x, const ISMRMRD::Acquisition& acq_y);
249  // the sum of the elements of x
250  static complex_float_t sum(const ISMRMRD::Acquisition& acq_x);
251  // the value of the element of x with the largest real part
252  static complex_float_t max(const ISMRMRD::Acquisition& acq_x);
253  // elementwise multiplication
254  // y := x .* y
255  static void multiply
256  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
257  // multiply by scalar
258  // y := x * y
259  static void multiply
260  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
261  // add scalar
262  // y := x + y
263  static void add
264  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
265  // elementwise division
266  // y := x ./ y
267  static void divide
268  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
269  // elementwise maximum
270  // y := std::real(x) > std::real(y) ? x : y
271  static void maximum
272  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
273  static void maximum
274  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
275  // elementwise minimum
276  // y := std::real(x) < std::real(y) ? x : y
277  static void minimum
278  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
279  static void minimum
280  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
281  // y := pow(x, y)
282  static void power
283  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
284  static void power
285  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y, complex_float_t y);
286  // y := exp(x)
287  static void exp
288  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
289  // y := log(x)
290  static void log
291  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
292  // y := sqrt(x)
293  static void sqrt
294  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
295  // y := sign(x) (x < 0: -1, x == 0: 0, x > 0: 1)
296  static void sign
297  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
298  // l2 norm of x
299  // y := abs(x)
300  static void abs
301  (const ISMRMRD::Acquisition& acq_x, ISMRMRD::Acquisition& acq_y);
302  static float norm(const ISMRMRD::Acquisition& acq_x);
303 
304  // type and dimension of an ISMRMRD::Acquisition parameter
305  static void ismrmrd_par_info(const char* par, int* output)
306  {
307  // default type: int, dimension: 1
308  output[0] = 0;
309  output[1] = 1;
310 
311  // check parameter type
312  if (sirf::iequals(par, "sample_time_us") ||
313  sirf::iequals(par, "position") ||
314  sirf::iequals(par, "read_dir") ||
315  sirf::iequals(par, "phase_dir") ||
316  sirf::iequals(par, "slice_dir") ||
317  sirf::iequals(par, "patient_table_position") ||
318  sirf::iequals(par, "user_float"))
319  output[0] = 1; // float
320 
321  // check parameter dimension
322  if (sirf::iequals(par, "physiology_time_stamp"))
323  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
324  else if (sirf::iequals(par, "channel_mask"))
325  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
326  else if (sirf::iequals(par, "position") ||
327  sirf::iequals(par, "read_dir") ||
328  sirf::iequals(par, "phase_dir") ||
329  sirf::iequals(par, "slice_dir") ||
330  sirf::iequals(par, "patient_table_position"))
331  output[1] = 3;
332  else if (sirf::iequals(par, "user_int") ||
333  sirf::iequals(par, "idx_user"))
334  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_INTS;
335  else if (sirf::iequals(par, "user_float"))
336  output[1] = ISMRMRD::ISMRMRD_Constants::ISMRMRD_USER_FLOATS;
337  }
338  // value of an ISMRMRD::Acquisition int parameter
339  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
340  const char* name, unsigned long long int* v)
341  {
342  if (sirf::iequals(name, "version"))
343  *v = ((unsigned int)acq.version());
344  else if (sirf::iequals(name, "flags"))
345  *v = ((unsigned long long int)acq.flags());
346  else if (sirf::iequals(name, "measurement_uid"))
347  *v = ((unsigned int)acq.measurement_uid());
348  else if (sirf::iequals(name, "scan_counter"))
349  *v = ((unsigned int)acq.scan_counter());
350  else if (sirf::iequals(name, "acquisition_time_stamp"))
351  *v = ((unsigned int)acq.acquisition_time_stamp());
352  else if (sirf::iequals(name, "number_of_samples"))
353  *v = ((unsigned int)acq.number_of_samples());
354  else if (sirf::iequals(name, "available_channels"))
355  *v = ((unsigned int)acq.available_channels());
356  else if (sirf::iequals(name, "active_channels"))
357  *v = ((unsigned int)acq.active_channels());
358  else if (sirf::iequals(name, "discard_pre"))
359  *v = ((unsigned int)acq.discard_pre());
360  else if (sirf::iequals(name, "discard_post"))
361  *v = ((unsigned int)acq.discard_post());
362  else if (sirf::iequals(name, "center_sample"))
363  *v = ((unsigned int)acq.center_sample());
364  else if (sirf::iequals(name, "encoding_space_ref"))
365  *v = ((unsigned int)acq.encoding_space_ref());
366  else if (sirf::iequals(name, "trajectory_dimensions"))
367  *v = ((unsigned int)acq.trajectory_dimensions());
368  else if (sirf::iequals(name, "kspace_encode_step_1"))
369  *v = ((unsigned int)acq.idx().kspace_encode_step_1);
370  else if (sirf::iequals(name, "kspace_encode_step_2"))
371  *v = ((unsigned int)acq.idx().kspace_encode_step_2);
372  else if (sirf::iequals(name, "average"))
373  *v = ((unsigned int)acq.idx().average);
374  else if (sirf::iequals(name, "slice"))
375  *v = ((unsigned int)acq.idx().slice);
376  else if (sirf::iequals(name, "contrast"))
377  *v = ((unsigned int)acq.idx().contrast);
378  else if (sirf::iequals(name, "phase"))
379  *v = ((unsigned int)acq.idx().phase);
380  else if (sirf::iequals(name, "repetition"))
381  *v = ((unsigned int)acq.idx().repetition);
382  else if (sirf::iequals(name, "set"))
383  *v = ((unsigned int)acq.idx().set);
384  else if (sirf::iequals(name, "segment"))
385  *v = ((unsigned int)acq.idx().segment);
386  else if (sirf::iequals(name, "physiology_time_stamp")) {
387  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_PHYS_STAMPS;
388  const uint32_t* pts = acq.physiology_time_stamp();
389  for (int i = 0; i < n; i++)
390  v[i] = (unsigned int)pts[i];
391  }
392  else if (sirf::iequals(name, "channel_mask")) {
393  int n = ISMRMRD::ISMRMRD_Constants::ISMRMRD_CHANNEL_MASKS;
394  const uint64_t* pts = acq.channel_mask();
395  for (int i = 0; i < n; i++)
396  v[i] = (unsigned long long int)pts[i];
397  }
398  }
399  // value of an ISMRMRD::Acquisition float parameter
400  static void ismrmrd_par_value(ISMRMRD::Acquisition& acq,
401  const char* name, float* v)
402  {
403  if (sirf::iequals(name, "sample_time_us"))
404  *v = acq.sample_time_us();
405  else if (sirf::iequals(name, "position")) {
406  float* u = acq.position();
407  for (int i = 0; i < 3; i++)
408  v[i] = u[i];
409  }
410  else if (sirf::iequals(name, "read_dir")) {
411  float* u = acq.read_dir();
412  for (int i = 0; i < 3; i++)
413  v[i] = u[i];
414  }
415  else if (sirf::iequals(name, "phase_dir")) {
416  float* u = acq.phase_dir();
417  for (int i = 0; i < 3; i++)
418  v[i] = u[i];
419  }
420  else if (sirf::iequals(name, "slice_dir")) {
421  float* u = acq.slice_dir();
422  for (int i = 0; i < 3; i++)
423  v[i] = u[i];
424  }
425  else if (sirf::iequals(name, "patient_table_position")) {
426  float* u = acq.patient_table_position();
427  for (int i = 0; i < 3; i++)
428  v[i] = u[i];
429  }
430  }
437  void set_encoding_limits(const std::string& name, \
438  const std::tuple<unsigned short,unsigned short,unsigned short> min_max_ctr)
439  {
440  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
441  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
442 
443  ISMRMRD::Limit limit;
444  limit.minimum = std::get<0>(min_max_ctr);
445  limit.maximum = std::get<1>(min_max_ctr);
446  limit.center = std::get<2>(min_max_ctr);
447 
448  if(sirf::iequals(name, "kspace_encoding_step_1"))
449  enc_limits.kspace_encoding_step_1.get() = limit;
450  else if(sirf::iequals(name, "kspace_encoding_step_2"))
451  enc_limits.kspace_encoding_step_2.get() = limit;
452  else if(sirf::iequals(name, "average"))
453  enc_limits.average.get() = limit;
454  else if(sirf::iequals(name, "slice"))
455  enc_limits.slice.get() = limit;
456  else if(sirf::iequals(name, "contrast"))
457  enc_limits.contrast.get() = limit;
458  else if(sirf::iequals(name, "phase"))
459  enc_limits.phase.get() = limit;
460  else if(sirf::iequals(name, "repetition"))
461  enc_limits.repetition.get() = limit;
462  else if(sirf::iequals(name, "set"))
463  enc_limits.set.get() = limit;
464  else if(sirf::iequals(name, "segment"))
465  enc_limits.segment.get() = limit;
466  else
467  throw std::runtime_error("You passed a name that is not an encoding limit.");
468 
469  hdr.encoding[0].encodingLimits = enc_limits;
470  std::stringstream serialised_hdr;
471  ISMRMRD::serialize(hdr, serialised_hdr);
472  this->set_acquisitions_info(AcquisitionsInfo(serialised_hdr.str()));
473  }
474 
475  std::tuple<unsigned short,unsigned short,unsigned short>
476  get_encoding_limits(const std::string& name) const
477  {
478  ISMRMRD::IsmrmrdHeader hdr = this->acquisitions_info().get_IsmrmrdHeader();
479  ISMRMRD::EncodingLimits enc_limits = hdr.encoding[0].encodingLimits;
480 
481  ISMRMRD::Limit limit;
482 
483  if(sirf::iequals(name, "kspace_encoding_step_1"))
484  limit = enc_limits.kspace_encoding_step_1.get();
485  else if(sirf::iequals(name, "kspace_encoding_step_2"))
486  limit = enc_limits.kspace_encoding_step_2.get();
487  else if(sirf::iequals(name, "average"))
488  limit = enc_limits.average.get();
489  else if(sirf::iequals(name, "slice"))
490  limit = enc_limits.slice.get();
491  else if(sirf::iequals(name, "contrast"))
492  limit = enc_limits.contrast.get();
493  else if(sirf::iequals(name, "phase"))
494  limit = enc_limits.phase.get();
495  else if(sirf::iequals(name, "repetition"))
496  limit = enc_limits.repetition.get();
497  else if(sirf::iequals(name, "set"))
498  limit = enc_limits.set.get();
499  else if(sirf::iequals(name, "segment"))
500  limit = enc_limits.segment.get();
501  else
502  throw std::runtime_error("You passed a name that is not an encoding limit.");
503 
504  return std::make_tuple(limit.minimum, limit.maximum, limit.center);
505  }
506 
507  // abstract methods
508 
509  virtual void empty() = 0;
510  virtual void take_over(MRAcquisitionData&) = 0;
511 
512  // the number of acquisitions in the container
513  virtual unsigned int number() const = 0;
514 
515  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
516  get_acquisition_sptr(unsigned int num) = 0;
517  virtual int get_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) const = 0;
518  virtual void set_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) = 0;
519  virtual void append_acquisition(ISMRMRD::Acquisition& acq) = 0;
520 
521  virtual void copy_acquisitions_info(const MRAcquisitionData& ac) = 0;
522  virtual void copy_acquisitions_data(const MRAcquisitionData& ac) = 0;
523 
524  // 'export' constructors: workaround for creating 'ABC' objects
525  virtual gadgetron::unique_ptr<MRAcquisitionData> new_acquisitions_container() = 0;
526  virtual MRAcquisitionData*
527  same_acquisitions_container(const AcquisitionsInfo& info) const = 0;
528 
529  virtual void set_data(const complex_float_t* z, int all = 1) = 0;
530  virtual void get_data(complex_float_t* z, int all = 1);
531 
532  virtual void set_user_floats(float const * const z, int const idx);
533 
534  virtual bool is_complex() const
535  {
536  return true;
537  }
538 
539  // acquisition data algebra
541  virtual void sum(void* ptr) const;
542  virtual void max(void* ptr) const;
543  virtual void dot(const DataContainer& dc, void* ptr) const;
544  complex_float_t dot(const DataContainer& a_x)
545  {
546  complex_float_t z;
547  dot(a_x, &z);
548  return z;
549  }
550  virtual void axpby(
551  const void* ptr_a, const DataContainer& a_x,
552  const void* ptr_b, const DataContainer& a_y);
553  virtual void xapyb(
554  const DataContainer& a_x, const DataContainer& a_a,
555  const DataContainer& a_y, const DataContainer& a_b);
556  virtual void xapyb(
557  const DataContainer& a_x, const void* ptr_a,
558  const DataContainer& a_y, const void* ptr_b)
559  {
560  axpby(ptr_a, a_x, ptr_b, a_y);
561  }
562  virtual void xapyb(
563  const DataContainer& a_x, const void* ptr_a,
564  const DataContainer& a_y, const DataContainer& a_b);
565  virtual void multiply(const DataContainer& x, const DataContainer& y);
566  virtual void divide(const DataContainer& x, const DataContainer& y);
567  virtual void maximum(const DataContainer& x, const DataContainer& y);
568  virtual void minimum(const DataContainer& x, const DataContainer& y);
569  virtual void power(const DataContainer& x, const DataContainer& y);
570  virtual void multiply(const DataContainer& x, const void* y);
571  virtual void add(const DataContainer& x, const void* ptr_y);
572  virtual void maximum(const DataContainer& x, const void* y);
573  virtual void minimum(const DataContainer& x, const void* y);
574  virtual void power(const DataContainer& x, const void* y);
575  virtual void exp(const DataContainer& x);
576  virtual void log(const DataContainer& x);
577  virtual void sqrt(const DataContainer& x);
578  virtual void sign(const DataContainer& x);
579  virtual void abs(const DataContainer& x);
580  virtual float norm() const;
581 
582  virtual void write(const std::string &filename) const;
583 
584  // regular methods
585  void binary_op(const DataContainer& a_x, const DataContainer& a_y,
586  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
587  void semibinary_op(const DataContainer& a_x, complex_float_t y,
588  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&, complex_float_t));
589  void unary_op(const DataContainer& a_x,
590  void(*f)(const ISMRMRD::Acquisition&, ISMRMRD::Acquisition&));
591 
592  AcquisitionsInfo acquisitions_info() const { return acqs_info_; }
593  void set_acquisitions_info(std::string info) { acqs_info_ = info; }
594  void set_acquisitions_info(const AcquisitionsInfo info) { acqs_info_ = info;}
595 
596  ISMRMRD::TrajectoryType get_trajectory_type() const;
597 
598  void set_trajectory_type(const ISMRMRD::TrajectoryType type);
599 
600  void set_trajectory(const uint16_t traj_dim, float* traj)
601  {
602  ISMRMRD::Acquisition acq;
603  for(int i=0; i<number(); ++i)
604  {
605  get_acquisition(i, acq);
606  const uint16_t num_samples = acq.number_of_samples();
607  const uint16_t num_channels = acq.active_channels();
608  acq.resize(num_samples,num_channels,traj_dim);
609  int const offset = i*traj_dim*num_samples;
610  acq.setTraj(traj + offset);
611  set_acquisition(i, acq);
612  }
613  }
614 
615  gadgetron::unique_ptr<MRAcquisitionData> clone() const
616  {
617  return gadgetron::unique_ptr<MRAcquisitionData>(this->clone_impl());
618  }
619 
620  bool undersampled() const;
621  int get_acquisitions_dimensions(size_t ptr_dim) const;
622  void get_kspace_dimensions(std::vector<size_t>& dims) const;
623  uint16_t get_trajectory_dimensions(void) const;
624 
625  void sort();
626  void sort_by_time();
627  bool sorted() const { return sorted_; }
628  void set_sorted(bool sorted) { sorted_ = sorted; }
629 
631 
635  std::vector<KSpaceSubset::SetType > get_kspace_order() const;
636 
638  std::vector<KSpaceSubset> get_kspace_sorting() const { return this->sorting_; }
639 
641 
647  void organise_kspace();
648 
649  virtual std::vector<int> get_flagged_acquisitions_index(const std::vector<ISMRMRD::ISMRMRD_AcquisitionFlags> flags) const;
650  virtual std::vector<int> get_slice_encoding_index(const unsigned kspace_encode_step_2) const;
651 
652 
653  virtual void get_subset(MRAcquisitionData& subset, const std::vector<int> subset_idx) const;
654  virtual void set_subset(const MRAcquisitionData &subset, const std::vector<int> subset_idx);
655 
656  std::vector<int> index() { return index_; }
657  const std::vector<int>& index() const { return index_; }
658 
659  int index(int i) const
660  {
661  const std::size_t ni = index_.size();
662  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni) || static_cast<unsigned>(i) >= number())
663  THROW("Aquisition number is out of range");
664  if (ni > 0)
665  return index_[i];
666  else
667  return i;
668  }
669 
679  void read(const std::string& filename_ismrmrd_with_ext, int all = 0);
680 
681  protected:
682  bool sorted_ = false;
683  std::vector<int> index_;
684  std::vector<KSpaceSubset> sorting_;
685  AcquisitionsInfo acqs_info_;
686 
687  // new MRAcquisitionData objects will be created from this template
688  // using same_acquisitions_container()
689  static gadgetron::shared_ptr<MRAcquisitionData> acqs_templ_;
690 
691  virtual MRAcquisitionData* clone_impl() const = 0;
692 
693  private:
694 
695  };
696 
705  class AcquisitionsVector : public MRAcquisitionData {
706  public:
707  AcquisitionsVector(const std::string& filename_with_ext, int all = 0)
708  {
709  this->read(filename_with_ext, all);
710  }
711 
712  AcquisitionsVector(const AcquisitionsInfo& info = AcquisitionsInfo())
713  {
714  acqs_info_ = info;
715  }
716  virtual void empty();
717  virtual void take_over(MRAcquisitionData& ad) {}
718  virtual unsigned int number() const { return (unsigned int)acqs_.size(); }
719  virtual unsigned int items() const { return (unsigned int)acqs_.size(); }
720  virtual void append_acquisition(ISMRMRD::Acquisition& acq)
721  {
722  acqs_.push_back(gadgetron::shared_ptr<ISMRMRD::Acquisition>
723  (new ISMRMRD::Acquisition(acq)));
724  }
725  virtual gadgetron::shared_ptr<ISMRMRD::Acquisition>
726  get_acquisition_sptr(unsigned int num)
727  {
728  int ind = index(num);
729  return acqs_[ind];
730  }
731  virtual int get_acquisition(unsigned int num, ISMRMRD::Acquisition& acq) const
732  {
733  int ind = index(num);
734  acq = *acqs_[ind];
735  if (!(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION) && \
736  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING) && \
737  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_LAST_IN_MEASUREMENT) && \
738  !(acq).isFlagSet(ISMRMRD::ISMRMRD_ACQ_IS_REVERSE) && \
739  (acq).flags() >= (1 << (ISMRMRD::ISMRMRD_ACQ_IS_NOISE_MEASUREMENT - 1)))
740  return 0;
741  return 1;
742  }
743  virtual void set_acquisition(unsigned int num, ISMRMRD::Acquisition& acq)
744  {
745  int ind = index(num);
746  *acqs_[ind] = acq;
747  }
748  virtual void copy_acquisitions_info(const MRAcquisitionData& ac)
749  {
750  acqs_info_ = ac.acquisitions_info();
751  }
752  virtual void copy_acquisitions_data(const MRAcquisitionData& ac);
753  virtual void set_data(const complex_float_t* z, int all = 1);
754 
755  virtual AcquisitionsVector* same_acquisitions_container
756  (const AcquisitionsInfo& info) const
757  {
758  return new AcquisitionsVector(info);
759  }
760  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
761  {
762  DataContainer* ptr = new AcquisitionsVector(acqs_info_);
763  return new ObjectHandle<DataContainer>
764  (gadgetron::shared_ptr<DataContainer>(ptr));
765  }
766  virtual gadgetron::unique_ptr<MRAcquisitionData>
767  new_acquisitions_container()
768  {
769  return gadgetron::unique_ptr<MRAcquisitionData>
770  (new AcquisitionsVector(acqs_info_));
771  }
772 
773  private:
774  std::vector<gadgetron::shared_ptr<ISMRMRD::Acquisition> > acqs_;
775  virtual AcquisitionsVector* clone_impl() const;
776  virtual void conjugate_impl();
777  };
778 
785  class ISMRMRDImageData : public ImageData {
786  public:
787  //ISMRMRDImageData(ISMRMRDImageData& id, const char* attr,
788  //const char* target); //does not build, have to be in the derived class
789 
790  virtual void empty() = 0;
791  virtual unsigned int number() const = 0;
792  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
793  (unsigned int im_num) = 0;
794  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
795  (unsigned int im_num) const = 0;
796 // virtual ImageWrap& image_wrap(unsigned int im_num) = 0;
797 // virtual const ImageWrap& image_wrap(unsigned int im_num) const = 0;
798  virtual void append(int image_data_type, void* ptr_image) = 0;
799  virtual void append(const ImageWrap& iw) = 0;
800  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw) = 0;
801  virtual gadgetron::shared_ptr<ISMRMRDImageData> abs() const = 0;
802  virtual gadgetron::shared_ptr<ISMRMRDImageData> real() const = 0;
803  virtual void clear_data() = 0;
804  virtual void set_image_type(int imtype) = 0;
805  virtual void get_data(complex_float_t* data) const;
806  virtual void set_data(const complex_float_t* data);
807  virtual void get_real_data(float* data) const;
808  virtual void set_real_data(const float* data);
809  virtual int read(std::string filename, std::string variable = "", int iv = -1);
810  virtual void write(const std::string &filename, const std::string &groupname, const bool dicom) const;
811  virtual void write(const std::string &filename) const
812  {
813  size_t size = filename.size();
814  std::string suff = filename.substr(size - 4, 4);
815  if (suff == std::string(".dcm")) {
816  std::string prefix = filename.substr(0, size - 4);
817  this->write(prefix, "", true);
818  }
819  else {
820  auto found = filename.find_last_of("/\\");
821  auto slash_found = found;
822  if (found == std::string::npos)
823  found = filename.find_last_of(".");
824  else
825  found = filename.substr(found + 1).find_last_of(".");
826  if (found == std::string::npos)
827  this->write(filename + ".h5", "", false);
828  else {
829  std::string ext = filename.substr(slash_found + found + 1);
830  if (ext == std::string(".h5"))
831  this->write(filename, "", false);
832  else
833  std::cerr << "WARNING: writing ISMRMRD images to "
834  << ext << "-files not implemented, "
835  << "please convert to Nifti images\n";
836  }
837  }
838  }
839  virtual Dimensions dimensions() const
840  {
841  Dimensions dim;
842  const ImageWrap& iw = image_wrap(0);
843  int d[5];
844  iw.get_dim(d);
845  dim["x"] = d[0];
846  dim["y"] = d[1];
847  dim["z"] = d[2];
848  dim["c"] = d[3];
849  dim["n"] = number();
850  return dim;
851  }
852  virtual void get_image_dimensions(unsigned int im_num, int* dim) const
853  {
854  if (im_num >= number())
855  dim[0] = dim[1] = dim[2] = dim[3] = 0;
856  const ImageWrap& iw = image_wrap(im_num);
857  iw.get_dim(dim);
858  }
859  bool check_dimension_consistency() const
860  {
861  size_t const num_dims = 4;
862  std::vector<int> first_img_dims(num_dims), temp_img_dims(num_dims);
863 
864  this->get_image_dimensions(0, &first_img_dims[0]);
865 
866  bool dims_match = true;
867  for(int i=1; i<number(); ++i)
868  {
869  this->get_image_dimensions(0, &temp_img_dims[0]);
870  dims_match *= (first_img_dims == temp_img_dims);
871  }
872  return dims_match;
873  }
874  virtual gadgetron::shared_ptr<ISMRMRDImageData>
875  new_images_container() const = 0;
876  virtual gadgetron::shared_ptr<ISMRMRDImageData>
877  clone(const char* attr, const char* target) = 0;
878  virtual int image_data_type(unsigned int im_num) const
879  {
880  return image_wrap(im_num).type();
881  }
882  virtual size_t num_data_elm() const
883  {
884  return image_wrap(0).num_data_elm();
885  }
886 
887  virtual float norm() const;
889  virtual void sum(void* ptr) const;
890  virtual void max(void* ptr) const;
891  virtual void dot(const DataContainer& dc, void* ptr) const;
892  virtual void axpby(
893  const void* ptr_a, const DataContainer& a_x,
894  const void* ptr_b, const DataContainer& a_y);
895  virtual void xapyb(
896  const DataContainer& a_x, const void* ptr_a,
897  const DataContainer& a_y, const void* ptr_b)
898  {
899  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
900  ComplexFloat_ b(*static_cast<const complex_float_t*>(ptr_b));
901  xapyb_(a_x, a, a_y, b);
902  }
903  virtual void xapyb(
904  const DataContainer& a_x, const void* ptr_a,
905  const DataContainer& a_y, const DataContainer& a_b)
906  {
907  ComplexFloat_ a(*static_cast<const complex_float_t*>(ptr_a));
908  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
909  xapyb_(a_x, a, a_y, b);
910  }
911  //virtual void xapyb(
912  // const DataContainer& a_x, const DataContainer& a_a,
913  // const DataContainer& a_y, const void* ptr_b)
914  //{
915  // SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
916  // ComplexFloat_ b(*(complex_float_t*)ptr_b);
917  // xapyb_(a_x, a, a_y, b);
918  //}
919  virtual void xapyb(
920  const DataContainer& a_x, const DataContainer& a_a,
921  const DataContainer& a_y, const DataContainer& a_b)
922  {
923  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, a, a_a);
924  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, b, a_b);
925  xapyb_(a_x, a, a_y, b);
926  }
927  virtual void multiply(const DataContainer& x, const DataContainer& y);
928  virtual void divide(const DataContainer& x, const DataContainer& y);
929  virtual void maximum(const DataContainer& x, const DataContainer& y);
930  virtual void minimum(const DataContainer& x, const DataContainer& y);
931  virtual void power(const DataContainer& x, const DataContainer& y);
932  virtual void multiply(const DataContainer& x, const void* ptr_y);
933  virtual void add(const DataContainer& x, const void* ptr_y);
934  virtual void maximum(const DataContainer& x, const void* ptr_y);
935  virtual void minimum(const DataContainer& x, const void* ptr_y);
936  virtual void power(const DataContainer& x, const void* ptr_y);
937  virtual void exp(const DataContainer& x);
938  virtual void log(const DataContainer& x);
939  virtual void sqrt(const DataContainer& x);
940  virtual void sign(const DataContainer& x);
941  virtual void abs(const DataContainer& x);
942 
943  void binary_op(
944  const DataContainer& a_x, const DataContainer& a_y,
945  complex_float_t(*f)(complex_float_t, complex_float_t));
946  void semibinary_op(
947  const DataContainer& a_x, complex_float_t y,
948  complex_float_t(*f)(complex_float_t, complex_float_t));
949  void unary_op(const DataContainer& a_x, complex_float_t(*f)(complex_float_t));
950 
951  void fill(float s);
952  void scale(float s);
953  complex_float_t dot(const DataContainer& a_x)
954  {
955  complex_float_t z;
956  dot(a_x, &z);
957  return z;
958  }
959  void axpby(
960  complex_float_t a, const DataContainer& a_x,
961  complex_float_t b, const DataContainer& a_y)
962  {
963  axpby(&a, a_x, &b, a_y);
964  }
965  void xapyb(
966  const DataContainer& a_x, complex_float_t a,
967  const DataContainer& a_y, complex_float_t b)
968  {
969  xapyb(a_x, &a, a_y, &b);
970  }
971  gadgetron::unique_ptr<ISMRMRDImageData> clone() const
972  {
973  return gadgetron::unique_ptr<ISMRMRDImageData>(this->clone_impl());
974  }
975 
976  virtual void sort() = 0;
977  bool sorted() const { return sorted_; }
978  void set_sorted(bool sorted) { sorted_ = sorted; }
979  std::vector<int> index() { return index_; }
980  const std::vector<int>& index() const { return index_; }
981  int index(int i) const
982  {
983  const std::size_t ni = index_.size();
984  if (i < 0 || (ni > 0 && static_cast<std::size_t>(i) >= ni) || static_cast<unsigned>(i) >= number())
985  THROW("Image number is out of range. You tried to look up an image number that is not inside the container.");
986  if (ni > 0)
987  return index_[i];
988  else
989  return i;
990  }
991  ImageWrap& image_wrap(unsigned int im_num)
992  {
993  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
994  return *sptr_iw;
995  }
996  const ImageWrap& image_wrap(unsigned int im_num) const
997  {
998  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
999  sptr_image_wrap(im_num);
1000  return *sptr_iw;
1001  }
1003  void set_meta_data(const AcquisitionsInfo &acqs_info);
1005  const AcquisitionsInfo &get_meta_data() const { return acqs_info_; }
1006 
1007  protected:
1008  bool sorted_=false;
1009  std::vector<int> index_;
1010  AcquisitionsInfo acqs_info_;
1012  virtual ISMRMRDImageData* clone_impl() const = 0;
1013  virtual void conjugate_impl();
1014 
1015  private:
1016  class ComplexFloat_ {
1017  public:
1018  ComplexFloat_(complex_float_t v) : v_(v) {}
1019  unsigned int number() const
1020  {
1021  return 0;
1022  }
1023  complex_float_t image_wrap(unsigned int i)
1024  {
1025  return v_;
1026  }
1027  size_t num_data_elm()
1028  {
1029  return 1;
1030  }
1031  private:
1032  complex_float_t v_;
1033  };
1034 
1035  template<class A, class B>
1036  void xapyb_(const DataContainer& a_x, A& a, const DataContainer& a_y, B& b)
1037  {
1038  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, x, a_x);
1039  SIRF_DYNAMIC_CAST(const ISMRMRDImageData, y, a_y);
1040  unsigned int nx = x.number();
1041  unsigned int na = a.number();
1042  unsigned int ny = y.number();
1043  unsigned int nb = b.number();
1044  //std::cout << nx << ' ' << ny << '\n';
1045  //std::cout << na << ' ' << nb << '\n';
1046  if (nx != ny)
1047  THROW("ImageData sizes mismatch in axpby");
1048  if (na > 0 && na != nx)
1049  THROW("ImageData sizes mismatch in axpby");
1050  if (nb > 0 && nb != nx)
1051  THROW("ImageData sizes mismatch in axpby");
1052  unsigned int n = number();
1053  if (n > 0) {
1054  if (n != nx)
1055  THROW("ImageData sizes mismatch in axpby");
1056  for (unsigned int i = 0; i < nx; i++)
1057  image_wrap(i).xapyb(x.image_wrap(i), a.image_wrap(i),
1058  y.image_wrap(i), b.image_wrap(i));
1059  }
1060  else {
1061  for (unsigned int i = 0; i < nx; i++) {
1062  const ImageWrap& u = x.image_wrap(i);
1063  const ImageWrap& v = y.image_wrap(i);
1064  ImageWrap w(u);
1065  w.xapyb(u, a.image_wrap(i), v, b.image_wrap(i));
1066  append(w);
1067  }
1068  }
1069  this->set_meta_data(x.get_meta_data());
1070  }
1071 
1072  };
1073 
1074  typedef ISMRMRDImageData GadgetronImageData;
1075 
1084  class GadgetronImagesVector : public GadgetronImageData {
1085  public:
1086  typedef ImageData::Iterator BaseIter;
1087  typedef ImageData::Iterator_const BaseIter_const;
1088  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator
1089  ImageWrapIter;
1090  typedef std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator
1091  ImageWrapIter_const;
1092  class Iterator : public ImageData::Iterator {
1093  public:
1094  Iterator(ImageWrapIter iw, int n, int i, const ImageWrap::Iterator& it) :
1095  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end())
1096  {}
1097  Iterator(const Iterator& iter) : iw_(iter.iw_), n_(iter.n_), i_(iter.i_),
1098  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1099  {}
1100  Iterator& operator=(const Iterator& iter)
1101  {
1102  iw_ = iter.iw_;
1103  n_ = iter.n_;
1104  i_ = iter.i_;
1105  iter_ = iter.iter_;
1106  end_ = iter.end_;
1107  sptr_iter_ = iter.sptr_iter_;
1108  return *this;
1109  }
1110  virtual bool operator==(const BaseIter& ai) const
1111  {
1112  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1113  return iter_ == i.iter_;
1114  }
1115  virtual bool operator!=(const BaseIter& ai) const
1116  {
1117  SIRF_DYNAMIC_CAST(const Iterator, i, ai);
1118  return iter_ != i.iter_;
1119  }
1120  Iterator& operator++()
1121  {
1122  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1123  throw std::out_of_range("cannot advance out-of-range iterator");
1124  ++iter_;
1125  if (iter_ == end_ && i_ < n_ - 1) {
1126  ++i_;
1127  ++iw_;
1128  iter_ = (**iw_).begin();
1129  end_ = (**iw_).end();
1130  }
1131  return *this;
1132  }
1133  Iterator& operator++(int)
1134  {
1135  sptr_iter_.reset(new Iterator(*this));
1136  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1137  throw std::out_of_range("cannot advance out-of-range iterator");
1138  ++iter_;
1139  if (iter_ == end_ && i_ < n_ - 1) {
1140  ++i_;
1141  ++iw_;
1142  iter_ = (**iw_).begin();
1143  end_ = (**iw_).end();
1144  }
1145  return *sptr_iter_;
1146  }
1147  NumRef& operator*()
1148  {
1149  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1150  throw std::out_of_range
1151  ("cannot dereference out-of-range iterator");
1152  return *iter_;
1153  }
1154  private:
1155  //std::vector<gadgetron::shared_ptr<ImageWrap> >::iterator iw_;
1156  ImageWrapIter iw_;
1157  int n_;
1158  int i_;
1159  ImageWrap::Iterator iter_;
1160  ImageWrap::Iterator end_;
1161  gadgetron::shared_ptr<Iterator> sptr_iter_;
1162  };
1163 
1164  class Iterator_const : public ImageData::Iterator_const {
1165  public:
1166  Iterator_const(ImageWrapIter_const iw, int n, int i,
1167  const ImageWrap::Iterator_const& it) :
1168  iw_(iw), n_(n), i_(i), iter_(it), end_((**iw).end_const())
1169  {}
1170  Iterator_const(const Iterator_const& iter) : iw_(iter.iw_),
1171  n_(iter.n_), i_(iter.i_),
1172  iter_(iter.iter_), end_(iter.end_), sptr_iter_(iter.sptr_iter_)
1173  {}
1174  Iterator_const& operator=(const Iterator_const& iter)
1175  {
1176  iw_ = iter.iw_;
1177  n_ = iter.n_;
1178  i_ = iter.i_;
1179  iter_ = iter.iter_;
1180  end_ = iter.end_;
1181  sptr_iter_ = iter.sptr_iter_;
1182  return *this;
1183  }
1184  bool operator==(const BaseIter_const& ai) const
1185  {
1186  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1187  return iter_ == i.iter_;
1188  }
1189  bool operator!=(const BaseIter_const& ai) const
1190  {
1191  SIRF_DYNAMIC_CAST(const Iterator_const, i, ai);
1192  return iter_ != i.iter_;
1193  }
1194  Iterator_const& operator++()
1195  {
1196  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1197  throw std::out_of_range("cannot advance out-of-range iterator");
1198  ++iter_;
1199  if (iter_ == end_ && i_ < n_ - 1) {
1200  ++i_;
1201  ++iw_;
1202  iter_ = (**iw_).begin_const();
1203  end_ = (**iw_).end_const();
1204  }
1205  return *this;
1206  }
1207  // causes crashes and very inefficient anyway
1208  //Iterator_const& operator++(int)
1209  //{
1210  // sptr_iter_.reset(new Iterator_const(*this));
1211  // if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1212  // throw std::out_of_range("cannot advance out-of-range iterator");
1213  // ++iter_;
1214  // if (iter_ == end_ && i_ < n_ - 1) {
1215  // ++i_;
1216  // ++iw_;
1217  // iter_ = (**iw_).begin_const();
1218  // end_ = (**iw_).end_const();
1219  // }
1220  // return *sptr_iter_;
1221  //}
1222  const NumRef& operator*() const
1223  {
1224  if (i_ >= n_ || (i_ == n_ - 1 && iter_ == end_))
1225  throw std::out_of_range
1226  ("cannot dereference out-of-range iterator");
1227  ref_.copy(*iter_);
1228  return ref_;
1229  //return *iter_;
1230  }
1231  private:
1232  //std::vector<gadgetron::shared_ptr<ImageWrap> >::const_iterator iw_;
1233  ImageWrapIter_const iw_;
1234  int n_;
1235  int i_;
1236  ImageWrap::Iterator_const iter_;
1237  ImageWrap::Iterator_const end_;
1238  mutable NumRef ref_;
1239  gadgetron::shared_ptr<Iterator_const> sptr_iter_;
1240  };
1241 
1242  GadgetronImagesVector() : images_()
1243  {}
1244 
1255  GadgetronImagesVector(const MRAcquisitionData& ad, const bool coil_resolved=false);
1256  GadgetronImagesVector(const GadgetronImagesVector& images);
1257  GadgetronImagesVector(GadgetronImagesVector& images, const char* attr,
1258  const char* target);
1259  virtual void empty()
1260  {
1261  images_.clear();
1262  }
1263  virtual unsigned int items() const
1264  {
1265  return (unsigned int)images_.size();
1266  }
1267  virtual unsigned int number() const
1268  {
1269  return (unsigned int)images_.size();
1270  }
1271  virtual void append(int image_data_type, void* ptr_image)
1272  {
1273  images_.push_back(gadgetron::shared_ptr<ImageWrap>
1274  (new ImageWrap(image_data_type, ptr_image)));
1275  }
1276  virtual void append(CFImage& img)
1277  {
1278  void* vptr_img = new CFImage(img);
1279  this->append(7, vptr_img);
1280  }
1281  virtual void append(const ImageWrap& iw)
1282  {
1283  images_.push_back(gadgetron::shared_ptr<ImageWrap>(new ImageWrap(iw)));
1284  }
1285  virtual void append(gadgetron::shared_ptr<ImageWrap> sptr_iw)
1286  {
1287  images_.push_back(sptr_iw);
1288  }
1289  virtual gadgetron::shared_ptr<GadgetronImageData> abs() const;
1290  virtual gadgetron::shared_ptr<GadgetronImageData> real() const;
1291  virtual void clear_data()
1292  {
1293  std::vector<gadgetron::shared_ptr<ImageWrap> > empty_data;
1294  images_.swap(empty_data);
1295  }
1296  virtual void sort();
1297  virtual gadgetron::shared_ptr<ImageWrap> sptr_image_wrap
1298  (unsigned int im_num)
1299  {
1300  int i = index(im_num);
1301  return images_.at(i);
1302  }
1303  virtual gadgetron::shared_ptr<const ImageWrap> sptr_image_wrap
1304  (unsigned int im_num) const
1305  {
1306  int i = index(im_num);
1307  return images_.at(i);
1308  }
1309 /* virtual ImageWrap& image_wrap(unsigned int im_num)
1310  {
1311  gadgetron::shared_ptr<ImageWrap> sptr_iw = sptr_image_wrap(im_num);
1312  return *sptr_iw;
1313  }
1314  virtual const ImageWrap& image_wrap(unsigned int im_num) const
1315  {
1316  const gadgetron::shared_ptr<const ImageWrap>& sptr_iw =
1317  sptr_image_wrap(im_num);
1318  return *sptr_iw;
1319  }
1320 */
1321  virtual ObjectHandle<DataContainer>* new_data_container_handle() const
1322  {
1323  return new ObjectHandle<DataContainer>
1324  (gadgetron::shared_ptr<DataContainer>(new_images_container()));
1325  }
1326  virtual gadgetron::shared_ptr<GadgetronImageData> new_images_container() const
1327  {
1328  gadgetron::shared_ptr<GadgetronImageData> sptr_img
1329  ((GadgetronImageData*)new GadgetronImagesVector());
1330  sptr_img->set_meta_data(get_meta_data());
1331  return sptr_img;
1332  }
1333  virtual gadgetron::shared_ptr<GadgetronImageData>
1334  clone(const char* attr, const char* target)
1335  {
1336  return gadgetron::shared_ptr<GadgetronImageData>
1337  (new GadgetronImagesVector(*this, attr, target));
1338  }
1339 
1340  virtual Iterator& begin()
1341  {
1342  ImageWrapIter iw = images_.begin();
1343  begin_.reset(new Iterator(iw, images_.size(), 0, (**iw).begin()));
1344  return *begin_;
1345  }
1346  virtual Iterator& end()
1347  {
1348  ImageWrapIter iw = images_.begin();
1349  int n = images_.size();
1350  for (int i = 0; i < n - 1; i++)
1351  ++iw;
1352  end_.reset(new Iterator(iw, n, n - 1, (**iw).end()));
1353  return *end_;
1354  }
1355  virtual Iterator_const& begin() const
1356  {
1357  ImageWrapIter_const iw = images_.begin();
1358  begin_const_.reset
1359  (new Iterator_const(iw, images_.size(), 0, (**iw).begin_const()));
1360  return *begin_const_;
1361  }
1362  virtual Iterator_const& end() const
1363  {
1364  ImageWrapIter_const iw = images_.begin();
1365  int n = images_.size();
1366  for (int i = 0; i < n - 1; i++)
1367  ++iw;
1368  end_const_.reset
1369  (new Iterator_const(iw, n, n - 1, (**iw).end_const()));
1370  return *end_const_;
1371  }
1372  virtual void set_image_type(int image_type);
1373  virtual void get_data(complex_float_t* data) const;
1374  virtual void set_data(const complex_float_t* data);
1375  virtual void get_real_data(float* data) const;
1376  virtual void set_real_data(const float* data);
1377 
1379  std::unique_ptr<GadgetronImagesVector> clone() const
1380  {
1381  return std::unique_ptr<GadgetronImagesVector>(this->clone_impl());
1382  }
1383 
1385  void print_header(const unsigned im_num);
1386 
1388  virtual bool is_complex() const;
1389 
1391  virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out);
1392 
1394  virtual void set_up_geom_info();
1395 
1396  private:
1398  virtual GadgetronImagesVector* clone_impl() const
1399  {
1400  return new GadgetronImagesVector(*this);
1401  }
1402 
1403  std::vector<gadgetron::shared_ptr<ImageWrap> > images_;
1404  mutable gadgetron::shared_ptr<Iterator> begin_;
1405  mutable gadgetron::shared_ptr<Iterator> end_;
1406  mutable gadgetron::shared_ptr<Iterator_const> begin_const_;
1407  mutable gadgetron::shared_ptr<Iterator_const> end_const_;
1408  };
1409 
1415  class CoilImagesVector : public GadgetronImagesVector
1416  {
1417  public:
1418  CoilImagesVector() : GadgetronImagesVector(){}
1419  void calculate(const MRAcquisitionData& ad);
1420  protected:
1421  std::unique_ptr<MRAcquisitionData> extract_calibration_data(const MRAcquisitionData& ad) const;
1422  gadgetron::shared_ptr<FourierEncoding> sptr_enc_;
1423  };
1424 
1439  class CoilSensitivitiesVector : public GadgetronImagesVector
1440  {
1441 
1442  public:
1443 
1444  CoilSensitivitiesVector() : GadgetronImagesVector(){}
1445  CoilSensitivitiesVector(MRAcquisitionData& ad) :
1446  GadgetronImagesVector(ad, true){}
1447  CoilSensitivitiesVector(const char * file)
1448  {
1449  throw std::runtime_error("This has not been implemented yet.");
1450  }
1451 
1452  void set_csm_smoothness(int s){csm_smoothness_ = s;}
1453 
1454  void calculate(CoilImagesVector& iv);
1455  void calculate(const MRAcquisitionData& acq)
1456  {
1457  CoilImagesVector ci;
1458  ci.calculate(acq);
1459  calculate(ci);
1460  }
1461 
1462  CFImage get_csm_as_cfimage(size_t const i) const;
1463  CFImage get_csm_as_cfimage(const KSpaceSubset::TagType tag, const int offset) const;
1464 
1465 
1466  void get_dim(size_t const num_csm, int* dim) const
1467  {
1468  GadgetronImagesVector::get_image_dimensions(num_csm, dim);
1469  }
1470 
1471  void forward(GadgetronImageData& img, const GadgetronImageData& combined_img)const;
1472  void backward(GadgetronImageData& combined_img, const GadgetronImageData& img)const;
1473 
1474  protected:
1475 
1476  void coilchannels_from_combined_image(GadgetronImageData& img, const GadgetronImageData& combined_img) const;
1477  void combine_images_with_coilmaps(GadgetronImageData& combined_img, const GadgetronImageData& img) const;
1478 
1479  void calculate_csm(ISMRMRD::NDArray<complex_float_t>& cm, ISMRMRD::NDArray<float>& img, ISMRMRD::NDArray<complex_float_t>& csm);
1480 
1481  private:
1482  int csm_smoothness_ = 0;
1483  void smoothen_(int nx, int ny, int nz, int nc, complex_float_t* u, complex_float_t* v, int* obj_mask, int w);
1484  void mask_noise_(int nx, int ny, int nz, float* u, float noise, int* mask);
1485  float max_diff_(int nx, int ny, int nz, int nc, float small_grad, complex_float_t* u, complex_float_t* v);
1486  float max_(int nx, int ny, int nz, float* u);
1487 
1488  };
1489 
1490 void match_img_header_to_acquisition(CFImage& img, const ISMRMRD::Acquisition& acq);
1491 
1492 }
1493 
1494 #endif
ObjectHandle
Definition: DataHandle.h:159
sirf::AcquisitionsInfo
Definition: gadgetron_data_containers.h:84
sirf::AcquisitionsVector
A vector implementation of the abstract MR acquisition data container class.
Definition: gadgetron_data_containers.h:705
sirf::CoilImagesVector
A coil images container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1416
sirf::CoilSensitivitiesVector
A coil sensitivities container based on the GadgetronImagesVector class.
Definition: gadgetron_data_containers.h:1440
sirf::DataContainer
Definition: DataContainer.h:42
sirf::GadgetronImagesVector::Iterator_const
Definition: gadgetron_data_containers.h:1164
sirf::GadgetronImagesVector::Iterator
Definition: gadgetron_data_containers.h:1092
sirf::GadgetronImagesVector
A vector implementation of the abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:1084
sirf::GadgetronImagesVector::is_complex
virtual bool is_complex() const
Is complex?
Definition: gadgetron_data_containers.cpp:2244
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:2345
sirf::GadgetronImagesVector::reorient
virtual void reorient(const VoxelisedGeometricalInfo3D &geom_info_out)
Reorient image. Requires that dimensions match.
Definition: gadgetron_data_containers.cpp:2252
sirf::GadgetronImagesVector::clone
std::unique_ptr< GadgetronImagesVector > clone() const
Clone and return as unique pointer.
Definition: gadgetron_data_containers.h:1379
sirf::GadgetronImagesVector::print_header
void print_header(const unsigned im_num)
Print header info.
Definition: gadgetron_data_containers.cpp:2207
sirf::ISMRMRDImageData
Abstract Gadgetron image data container class.
Definition: gadgetron_data_containers.h:785
sirf::ISMRMRDImageData::set_meta_data
void set_meta_data(const AcquisitionsInfo &acqs_info)
Set the meta data.
Definition: gadgetron_data_containers.cpp:1983
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:903
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:1928
sirf::ISMRMRDImageData::multiply
virtual void multiply(const DataContainer &x, const DataContainer &y)
*this = the elementwise product x*y
Definition: gadgetron_data_containers.cpp:1570
sirf::ISMRMRDImageData::exp
virtual void exp(const DataContainer &x)
*this = the elementwise exp(x)
Definition: gadgetron_data_containers.cpp:1656
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:895
sirf::ISMRMRDImageData::minimum
virtual void minimum(const DataContainer &x, const DataContainer &y)
*this = the elementwise min(x, y)
Definition: gadgetron_data_containers.cpp:1620
sirf::ISMRMRDImageData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1691
sirf::ISMRMRDImageData::maximum
virtual void maximum(const DataContainer &x, const DataContainer &y)
*this = the elementwise max(x, y)
Definition: gadgetron_data_containers.cpp:1602
sirf::ISMRMRDImageData::sqrt
virtual void sqrt(const DataContainer &x)
*this = the elementwise sqrt(x)
Definition: gadgetron_data_containers.cpp:1670
sirf::ISMRMRDImageData::log
virtual void log(const DataContainer &x)
*this = the elementwise log(x)
Definition: gadgetron_data_containers.cpp:1663
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:1447
sirf::ISMRMRDImageData::divide
virtual void divide(const DataContainer &x, const DataContainer &y)
*this = the elementwise ratio x / y
Definition: gadgetron_data_containers.cpp:1594
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:919
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:1420
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:1586
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:1463
sirf::ISMRMRDImageData::power
virtual void power(const DataContainer &x, const DataContainer &y)
*this = the elementwise pow(x, y)
Definition: gadgetron_data_containers.cpp:1638
sirf::ISMRMRDImageData::sign
virtual void sign(const DataContainer &x)
*this = the elementwise sign(x)
Definition: gadgetron_data_containers.cpp:1677
sirf::ISMRMRDImageData::get_meta_data
const AcquisitionsInfo & get_meta_data() const
Get the meta data.
Definition: gadgetron_data_containers.h:1005
sirf::ISMRMRDImageData::sum
virtual void sum(void *ptr) const
below all void* are actually complex_float_t*
Definition: gadgetron_data_containers.cpp:1434
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:144
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:1367
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:1386
sirf::MRAcquisitionData
Abstract MR acquisition data container class.
Definition: gadgetron_data_containers.h:217
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:556
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:1136
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:638
sirf::MRAcquisitionData::norm
virtual float norm() const
returns the norm of this container viewed as a vector
Definition: gadgetron_data_containers.cpp:1050
sirf::MRAcquisitionData::read
void read(const std::string &filename_ismrmrd_with_ext, int all=0)
Reader for ISMRMRD::Acquisition from ISMRMRD file. filename_ismrmrd_with_ext: filename of ISMRMRD raw...
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:1166
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:437
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