libodin-dev 1.8.8-2ubuntu1 / usr / include / odinpara / reco.h

/***************************************************************************
                          reco.h  -  Parameters for image reconstruction
                             -------------------
    begin                : Thu May 06 2004
    copyright            : (C) 2000-2014 by Thies H. Jochimsen
    email                : thies@jochimsen.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

#ifndef RECO_H
#define RECO_H

#include <tjutils/tjvallist.h>

#include <odinpara/protocol.h>


/**
  * @addtogroup odinpara
  * @{
  */

// some limits
#define MAX_NUMOF_READOUT_SHAPES        10
#define MAX_NUMOF_KSPACE_TRAJS          10
#define MAX_NUMOF_ADC_WEIGHTING_VECTORS 10

#define PROTOCOL_BLOCK_LABEL "Protocol"

/**
  * Flags(bits) to mark ADCs:
  * - recoLastInChunkBit:    last ADC in a collection of consecutive ADCs generated by the system, an undefined number of zeroes may follow
  * - recoReflectBit:        ADC is reflected
  * - recoEpiNavPolarityBit: Indicates that first echo in EPI navigator readout is reflected
  */
enum recoFlags {
  recoLastInChunkBit    =1,
  recoReflectBit        =2,
  recoEpiNavPolarityBit =4,
  recoAllBits          =255
};


///////////////////////////////////////////////////////////

/**
  * Dimensions of acquired data, order is reversed for easy indexing of C-style arrays
  * - userdef:    Extra dimension which can be used by the sequence programmer
  * - te:         Echo-time index in a multi-echo acquistion, e.g. for field-map calculation
  * - dti:        Direction index for diffusion-tensor imaging
  * - average:    Average index to add up ADCs
  * - cycle:      Cycle index (spiral imaging, PROPELLER blades)
  * - slice:      Slice index
  * - line3d:     3D phase encoding index
  * - line:       2D phase encoding index
  * - echo:       Echo index in a multi-echo readout, e.g. an CPMG or EPI echo train
  * - epi:        Index of EPI readout train, e.g. in a 3D sequence
  * - templtype:  Type of template prescan
  * - navigator:  Type of navigator scan
  * - freq:       Frequency offset (will be used by multi-frequency reconstruction in odinreco)
  * - channel:    Channel index of phased-array coil
  * - readout:    Readout index
  * - repetition: Repetition dimension
  *
  */
enum recoDim                      { userdef=0, te,  dti,  average,  cycle,  slice,  line3d,  line,  echo,  epi,  templtype,  navigator,  freq,  channel,  readout,  repetition, n_recoDims};
static const char* recoDimLabel[]={"userdef", "te","dti","average","cycle","slice","line3d","line","echo","epi","templtype","navigator","freq","channel","readout","repetition"};
static const int n_recoIndexDims=templtype+1; // only the first dims are required for indexing of ADCs, the remaining are used only by odinreco


///////////////////////////////////////////////////////////


/**
  *
  * Enum to distinguish between different types of template prescan ADCs:
  * - no_template:        Not a template ADC
  * - phasecorr_template: Template ADC for 1D phase correction
  * - fieldmap_template:  Template to calculate fieldmap, the echo dimension will index the different TEs
  * - grappa_template:    Training data to calculate GRAPPA weights
  */
enum templateType                    {no_template=0, phasecorr_template, fieldmap_template, grappa_template, n_templateTypes};
static const char templateTypeChar[]={'N',           'P',                'F',               'G'};

///////////////////////////////////////////////////////////


/**
  *
  * Enum to distinguish between different types of navigator ADCs:
  * - no_navigator:  Not a navigator ADC
  * - epi_navigator: Siemens EPI phase correction navigator
  */
enum navigatorType                    {no_navigator=0, epi_navigator, n_navigatorTypes};
static const char navigatorTypeChar[]={'n',            'e'};


///////////////////////////////////////////////////////////


/**
  *
  * This struct holds the k-space position of one ADC
  */
struct kSpaceCoord {

  kSpaceCoord() {reset2defaults();}

  // for STD_list
  bool operator == (const kSpaceCoord& ri) const;
  bool operator <  (const kSpaceCoord& ri) const;


  static int string2index(const STD_string& str, recoDim dim);
  static STD_string index2string(int index, recoDim dim, int numof=0);

  STD_string printcoord(const unsigned short* numof_cache) const;
//  friend STD_ostream& operator << (STD_ostream& s,const kSpaceCoord& kcoord) { return s << kcoord.printcoord();}
  bool parsecoord(const STD_string& str);
  void reset2defaults();

  static STD_string print_header(const unsigned short* numof_cache);


  mutable unsigned int number;  // unique index of the ADC
  unsigned int reps;    // number of times this ADC is repeated
  unsigned short adcSize; // size with oversampling
  unsigned char channels; // coil channels
  unsigned short preDiscard;  // points to discard at the beginning of an ADC
  unsigned short postDiscard; // points to discard at the end of an ADC
  unsigned short concat; // number of concatenated acquisition this ADC is composed of
  float oversampling; // oversampling factor
  float relcenter; // relative position of k-space center

  short readoutIndex; // index of the readout shape to be used for ramp-sampling correction, -1 if none
  short trajIndex;   // index of the used k-space trajectory, -1 if none
  short weightIndex;   // index of the used ADC weighting vector, -1 if none
  short dtIndex; // index of the dwell time used

  unsigned short index[n_recoIndexDims]; // multi-dimensional reco index to locate ADC in k-space

  unsigned char flags; // recoFlags



 private:
  friend class JDXkSpaceCoords;

  static void assign_parsepos(const STD_string& header);


  static int parsepos_number;
  static int parsepos_reps;
  static int parsepos_adcSize;
  static int parsepos_channels;
  static int parsepos_preDiscard;
  static int parsepos_postDiscard;
  static int parsepos_concat;
  static int parsepos_oversampling;
  static int parsepos_relcenter;
  static int parsepos_readoutIndex;
  static int parsepos_trajIndex;
  static int parsepos_weightIndex;
  static int parsepos_dtIndex;
  static int parsepos_index[n_recoIndexDims];
  static int parsepos_lastinchunk;
  static int parsepos_reflect;
  static int max_parsepos;
  
};

//////////////////////////////////////////////////////////////////////


/**
  *
  * This class holds the assignment of consecutive ADCs to k-space positions
  */
class JDXkSpaceCoords : public JcampDxClass {

  enum JDXkSpaceCoords_state {coords_in_list, has_vec_cache, has_vec_alloc};

 public:
  JDXkSpaceCoords();
  ~JDXkSpaceCoords() {clear();}

  unsigned int size() const {create_vec_cache(); return vec_cache.size();}

  const kSpaceCoord& operator [] (unsigned int i) const {create_vec_cache(); return *(vec_cache[i]);}

  void clear(); // clear list before retrieving k-space ordering from sequence

  JDXkSpaceCoords& append_coord(const kSpaceCoord& coord);

 private:
  friend class RecoPars;

  void create_vec_cache() const;

  // overwriting virtual functions from JcampDxClass
  bool parsevalstring (const STD_string& parstring);
  STD_string printvalstring() const;
  STD_ostream& print2stream(STD_ostream& os) const;
  const char* get_typeInfo() const {return "kSpaceCoords";}
  JcampDxClass* create_copy() const {return new JDXkSpaceCoords(*this);}

  mutable STD_list<kSpaceCoord> coordlist;

  mutable JDXkSpaceCoords_state state;
  mutable STD_vector<kSpaceCoord*> vec_cache; // use pointers for better performance
  mutable unsigned short numof_cache[n_recoIndexDims];
};

//////////////////////////////////////////////////////////////////////

/**
  *
  * Structure to return ordering of k-space data while traversing the sequence tree
  */
struct RecoValList : public ValList<int> {
  RecoValList(const STD_string& object_label="unnamedRecoValList", unsigned int repetitions=1) : ValList<int>(object_label,repetitions) {}
};

//////////////////////////////////////////////////////////////////////

/**
  *
  * This class holds the assignment of consecutive ADCs to k-space positions
  */
class JDXrecoValList : public RecoValList, public JcampDxClass {

 public:
  JDXrecoValList(const STD_string& ldrlabel="unnamedJDXrecoValList");
  JDXrecoValList(const JDXrecoValList& jdrvl) {JDXrecoValList::operator = (jdrvl);}
  JDXrecoValList& operator = (const JDXrecoValList& jdrvl);

  JDXrecoValList& operator = (const RecoValList& rvl); // for assignment from get_recovallist in SeqMethod

  // overwriting virtual functions from JcampDxClass
  bool parsevalstring (const STD_string& parstring);
  STD_string printvalstring() const;
  STD_ostream& print2stream(STD_ostream& os) const;
  const char* get_typeInfo() const {return "RecoValList";}
  JcampDxClass* create_copy() const {return new JDXrecoValList(*this);}

};


//////////////////////////////////////////////////////////////////////



/**
  * This class is used to hold all information about the reconstruction
  */
class RecoPars : public JcampDxBlock {

 public:

/**
  * Constructs a RecoPars with the given label
  */
  RecoPars(const STD_string& label="unnamedRecoPars");

/**
  * Constructs a copy of 'sr'
  */
  RecoPars(const RecoPars& sr);

/**
  * Assignment operator
  */
  RecoPars& operator = (const RecoPars& sr);

/**
  * Returns the scan protocol.
  */
  const Protocol& get_protocol() const {return prot;}

/**
  * Returns the const assignment of ADC indices to k-space lines (for odinreco)
  */
  const JDXkSpaceCoords& get_kSpaceCoords() const {return kSpaceCoords;}

/**
  * Returns the data type of the acquired ADCs
  */
  const STD_string& get_DataFormat() const {return DataFormat;}

/**
  * Returns whether the of the acquired ADCs is little endian
  */
  bool is_LittleEndian() const {return LittleEndian;}

/**
  * Returns the file name of raw data
  */
  const STD_string& get_RawFile() const {return RawFile;}

/**
  * Returns the header size of raw data
  */
  unsigned int get_RawHeaderSize() const {return RawHeaderSize;}

/**
  * Returns the relative offset which should be added by the reco
  */
  const fvector& get_RelativeOffset() const {return RelativeOffset;}

/**
  * Extra command-line arguments for processing of final images
  */
  const STD_string& get_ImageProc() const {return ImageProc;}

/**
  * Returns the value vector attached to reco dimension 'dim'
  */
  const dvector& get_DimValues(recoDim dim) const {return DimValues[dim];}

/**
  * Returns the scaling factors for ADC channels
  */
  fvector get_ChannelScales() const {return ChannelScaling;}

/**
  * Returns the i'th dwell time (including oversampling)
  */
  double get_DwellTime(unsigned int i) const;


/**
  * Returns the number of ADC chunks, a chunk is a collection of
  * one or more ADCs stored by the system in one block. It can be
  * followed by an arbitrary number of zeroes.
  */
  int get_NumOfAdcChunks() const;


/**
  * Attaches a vector 'vals' of values to each step in the given dimension 'dim'
  */
  RecoPars& set_DimValues(recoDim dim, const dvector& vals);


/**
  * Appends a new ADC readout-gradient shape and its target grid size. Returns its index
  */
  int append_readout_shape(const fvector& shape, unsigned int dstsize);

/**
  * Appends a new k-space trajectory and returns its index
  */
  int append_kspace_traj(const farray& kspace_traj);

/**
  * Appends a new ADC weighting vector and returns its index
  */
  int append_adc_weight_vec(const cvector& weightvec);

/**
  * Appends a new ADC dwell time (including oversampling) and returns its index
  */
  int append_dwell_time(double dt);


/**
  * Returns the i'th ADC readout-gradient and its target grid size
  */
  void get_ReadoutShape(unsigned int i, fvector& shape, unsigned int& dstsize) const;


/**
  * Returns the number of k-space trajectories
  */
  unsigned int numof_kSpaceTraj() const;

/**
  * Returns the i'th k-space trajectory
  */
  const farray& get_kSpaceTraj(unsigned int i) const;


/**
  * Returns the i'th complex weighting vector with which the corresponding ADCs will be multiplied
  */
  const cvector& get_AdcWeightVector(unsigned int i) const;


/**
  * Returns the total number of ADC samples, if 'discard' is set to 'true'
  * only the effective ADC size is returned, i.e. without 'preDiscard' and
  * 'postDiscard'
  */
  LONGEST_INT get_TotalNumOfSamples(bool discard=false) const;

/**
  * Returns the total number of ADCs
  */
  unsigned int numof_adcs() const {return kSpaceOrdering.size();}

/**
  * Returns the i'th ADC of the sequence
  */
  const kSpaceCoord& get_kSpaceCoord(unsigned int i) const;

/**
  * Set a string describing the reconstruction pipeline (chain of steps/functors).
  */
  RecoPars& set_Recipe(const STD_string& recipe) {Recipe=recipe; return *this;}

/**
  * Returns a string describing the reconstruction pipeline (chain of steps/functors).
  */
  STD_string get_Recipe() const {return Recipe;}

/**
  * The extra chain of reconstruction steps (functors) described by 'recipe'
  * wil be performed after 3D reconstruction (including channel combination) but
  * before data is stored on disk.
  */
  RecoPars& set_PostProc3D(const STD_string& recipe) {PostProc3D=recipe; return *this;}

/**
  * Returns the extra chain of reconstruction steps (functors) which will be
  * performed after 3D reconstruction (including channel combination) but before
  * data is stored on disk.
  */
  STD_string get_PostProc3D() const {return PostProc3D;}

/**
  * The extra chain of reconstruction steps (functors) described by 'recipe'
  * wil be performed after 3D FFT but before channel combination.
  */
  RecoPars& set_PreProc3D(const STD_string& recipe) {PreProc3D=recipe; return *this;}

/**
  * Returns the extra chain of reconstruction steps (functors) described by 'recipe'
  * wil be performed after 3D FFT but before channel combination.
  */
  STD_string get_PreProc3D() const {return PreProc3D;}

/**
  * Extra command-line options for the reco.
  */
  RecoPars& set_CmdLineOpts(const STD_string& opts) {CmdLineOpts=opts; return *this;}

/**
  * Returns the extra chain of reconstruction steps (functors) which will be
  * performed after 3D reconstruction before data is stored on disk.
  */
  STD_string get_CmdLineOpts() const {return CmdLineOpts;}


 private:
  friend class SeqAcq;
  friend class SeqEpiDriverParavision;
  friend class SeqMethod;

  void common_init();

  void reset();
  void create_cache() const;

  void append_all_members();

  Protocol  prot;

  JDXstring DataFormat;
  JDXbool   LittleEndian;
  JDXstring RawFile;
  JDXint    RawHeaderSize;
  JDXtriple RelativeOffset;
  JDXstring ImageProc;
  JDXfloatArr ChannelScaling;
  JDXdoubleArr DwellTime;
  JDXfloatArr   ReadoutShape[MAX_NUMOF_READOUT_SHAPES];
  JDXintArr     ReadoutDstSize;
  JDXfloatArr   kSpaceTraj[MAX_NUMOF_KSPACE_TRAJS];
  JDXcomplexArr AdcWeightVector[MAX_NUMOF_ADC_WEIGHTING_VECTORS];
  JDXdoubleArr DimValues[n_recoIndexDims];
  JDXstring Recipe;
  JDXstring PreProc3D;
  JDXstring PostProc3D;
  JDXstring CmdLineOpts;

  JDXkSpaceCoords kSpaceCoords;
  JDXrecoValList kSpaceOrdering;

  mutable STD_vector<int> kSpaceOrdering_cache; // fast cache of ADC indices
  mutable bool cache_is_up2date;

};



//////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////



/**
  * Class to hold complex sensitivity
  */
class CoilSensitivity : public JcampDxBlock {

 public:

/**
  * Constructs a CoilSensitivity with the given label
  */
  CoilSensitivity(const STD_string& label="unnamedCoilSensitivity");

/**
  * Constructs a copy of 'cs'
  */
  CoilSensitivity(const CoilSensitivity& cs);

/**
  * Assignment operator
  */
  CoilSensitivity& operator = (const CoilSensitivity& cs);


/**
  * Accepts a 4-dim complex array 'sens_map' with the dimensions
  * (coil,z,y,x) and the corresponding FOVs
  */
  CoilSensitivity& set_sensitivity_map(const carray& sens_map, float FOVx, float FOVy, float FOVz);

/**
  * Returns complex sensitivity value for the given channel and spatial position,
  * uses trilinear interpolation
  */
  STD_complex get_sensitivity_value(unsigned int channel, float x, float y, float z) const;

/**
  * Returns the number of channels of the coil
  */
  unsigned int get_numof_channels() const {return SensitivityMap.get_extent()[0];}


 private:
  void append_all_members();

  JDXtriple     FOV;
  JDXcomplexArr SensitivityMap;
};


//////////////////////////////////////////////////////////////////////



/** @}
  */



#endif