libodin-dev 1.8.8-1 / usr / include / odinpara / geometry.h

/***************************************************************************
                          geometry.h  -  description
                             -------------------
    begin                : Wed Apr 17 2002
    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 GEOMETRY_H
#define GEOMETRY_H


#include <odinpara/jdxblock.h>
#include <odinpara/jdxnumbers.h>
#include <odinpara/jdxarrays.h>

#include <odinpara/odinpara.h>


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



/**
  * Enum to specify which boundary of a slice to take for a given direction
  */
enum sliceBoundary {lowerBound=0,upperBound,n_boundaries};


/**
  * The principal slice direction in the anatomical system of reference
  * - sagittal: Slice is perpendicular to left-right direction
  * - coronal: Slice is perpendicular to anterior-posterior direction
  * - axial: Slice is perpendicular to head-foot direction
  */
enum sliceOrientation {sagittal=0,coronal,axial,n_orientations};


/**
  * The geometry mode
  * - slicepack: A stack of equidistant parallel slices
  * - voxel_3d: A single 3D volume
  */
enum geometryMode {slicepack=0,voxel_3d,n_geometry_modes};


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


/**
  * \brief  Rotation Matrix
  *
  * This class represents 3 by 3 rotation matrix
  */
class RotMatrix : public virtual Labeled {
 public:

/**
  * default constructor for a 3 by 3  rotation matrix
  */
  RotMatrix(const STD_string& object_label = "unnamedRotMatrix" );

/**
  * constructs a rotation matrix from an existing rotation matrix
  */
   RotMatrix(const RotMatrix& sct);

/**
  * returns a row of rotation matrix
  */
   dvector& operator [] (unsigned int index) {return matrix[index%3];} // make sure index is lower than 3

/**
  * returns a row of rotation matrix
  */
   const dvector& operator [] (unsigned int index) const {return matrix[index%3];} // make sure index is lower than 3

/**
  * assignment operator
  */
  RotMatrix& operator = (const RotMatrix& sct);

/**
  * returns true if srm is element-wise equal to this. Two elements are
  * considered equal if their difference does not exceed GEO_CHECK_LIMIT.
  *
  */
  bool operator == (const RotMatrix& srm) const;

  // dummy comparison operator for lists
  bool operator <  (const RotMatrix& srm) const;

/**
  * returns the product this * vec
  */
  dvector operator * (const dvector& vec) const;

/**
  * returns the product this * matrix
  */
  RotMatrix operator * (const RotMatrix& matrix) const;

/**
  * Returns the matrix as a 2-dim farray
  */
  operator farray () const;

/**
  * Create anti-clockwise in-plane rotation matrix with angle 'phi' (in rad).
  */
  RotMatrix& set_inplane_rotation(float phi);


/**
  * returns a string describing the rotation matrix
  */
  STD_string print() const;

 private:
  friend class RotMatrixVector;

  class rowVec : public dvector {
   public:
    rowVec() : dvector(3) {};
    ~rowVec() {};
  };

//  bool check_and_correct();

  rowVec matrix[3];

//  static dvector returndummy;
};


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


/**
  * \brief  Geometry Settings
  *
  * Depending on the selected geometryMode, this class describes a pack
  * of parallel and congruent slices, or a 3D volume, .i.e. a voxel.
  *
  * The patients x,y,z coordinate system used for the gradients is as follows:
  * If standing in front of the magnet with the feet of the patient
  * pointing towards you, the axes are oriented
  *
  * - x: from right to left (in the final image left/right will then have the correct polarity)
  * - y: from bottom to top (floor to ceiling)
  * - z: pointing away from you to the patients head
  *
  *
  * The three orientation angles 'heightAngle', 'azimutAngle', and 'inplaneAngle'
  * (in degree) describe the orientation of the slicepack/voxel relative to
  * the x/y/z-coordinate system of the gradients.
  * If all these angles are zero, the read/phase/slice-system
  * of the slicepack/voxel is the same as the x/y/z-system of
  * the gradients. Otherwise, the system is rotated by the following
  * angles (in the given order):
  *
  * - azimutAngle: Anti-clockwise rotation around y-axis.
  * - heightAngle: Angle between slice vector (which is perp. to the slice) and x-z-plane
  * - inplaneAngle: Clockwise rotation around slice vector.
  *
  */
class Geometry : public JcampDxBlock {

public:

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

/**
  * Copy constructor
  */
  Geometry(const Geometry& ia);

/**
  * Assignment operator
  */
  Geometry& operator = (const Geometry& ia);

/**
  * Sets the geometry selection mode
  */
  Geometry& set_Mode(geometryMode mode);

/**
  * Returns the geometry selection mode
  */
  geometryMode get_Mode() const {return geometryMode(int(Mode));}


/**
  * Sets the Field Of View in the specified direction
  */
  Geometry& set_FOV(direction dir,double fov);

/**
  * Returns the Field Of View in the specified direction
  */
  double get_FOV(direction dir) const;

/**
  * Sets the spatial offset in the specified direction
  */
  Geometry& set_offset(direction dir,double offset);

/**
  * Returns the spatial offset in the specified direction
  */
  double get_offset(direction dir) const;

/**
  * Sets the number of slices of the current slice pack
  */
  Geometry& set_nSlices(unsigned int nslices);

/**
  * Returns the number of slices of the current slice pack
  */
  unsigned int get_nSlices() const {return nSlices;};

/**
  * Sets the slice thickness of the current slice pack
  */
  Geometry& set_sliceThickness(double thick);

/**
  * Sets the slice interslice distance
  */
  Geometry& set_sliceDistance(double dist);

/**
  * Returns the slice interslice distance
  */
  double get_sliceDistance() const {return sliceDistance;}

/**
  * Returns the normalised vector for the read direction in the laboratory system of reference
  */
  dvector get_readVector() const;

/**
  * Returns the normalised vector for the phase direction in the laboratory system of reference
  */
  dvector get_phaseVector() const;

/**
  * Returns the normalised vector for the slice direction in the laboratory system of reference
  */
  dvector get_sliceVector() const;

/**
  * Returns the vector of spatial offsets in slice direction for the different slices of the slicepack.
  */
  dvector get_sliceOffsetVector() const;

/**
  * Returns the slice thickness for the slices
  */
  double get_sliceThickness() const {return sliceThickness;}


/**
  * Returns a 5-dimensional array that contains the corner points of the slices/voxels where indexing is as follows:
  * (unsigned int slice,sliceBoundary boundary1,sliceBoundary boundary2,sliceBoundary boundary3, axis direction) with
  * - slice: The slice for the point
  * - boundary1: The boundary in the read direction
  * - boundary2: The boundary in the phase direction
  * - boundary3: The boundary in the slice direction, the size is 1 in slicepack mode
  * - direction: The x,y,z index
  *
  * The points are transformed to slice 'backgrslice' the coordinate system of 'background'.
  *
  */
  darray get_cornerPoints(const Geometry& background, unsigned int backgrslice) const;


/**
  * Set the angle parameters to match the specified anatomical slice orientation
  */
  Geometry& set_orientation(sliceOrientation orientation);


/**
  * Returns the principal anatomical slice orientation of slice normal 'svec'.
  */
  static sliceOrientation get_slice_orientation(const dvector& svec);

/**
  * Returns the principal anatomical slice orientation
  */
  sliceOrientation get_orientation() const {return get_slice_orientation(get_sliceVector());}


/**
  * Returns the orientation of the slice, i.e. the three orientation angles
  * 'heightAng', 'azimutAng' and 'inplaneAng'.
  * 'revSlice' returns whether the slice direction (handness) is reversed.
  */
  void get_orientation(double& heightAng, double& azimutAng, double& inplaneAng, bool& revSlice) const;

/**
  * Specifies the orientation of the slice pack by giving the three orientation angles
  * 'heightAng', 'azimutAng' and 'inplaneAng'.
  * Optionally, the slice direction can be reversed by setting 'revSlice' to true.
  */
  Geometry& set_orientation(double heightAng, double azimutAng, double inplaneAng, bool revSlice=false);

/**
  * Specifies the orientation of the slicepack/voxel by giving the three
  * vectors 'readvec', 'phasevec', and 'slicevec' which point in read,
  * phase, and slice direction. The vector 'offset' must contain the spatial
  * offset of the center of the slicepack/voxel.
  */
  Geometry& set_orientation_and_offset(const dvector& readvec, const dvector& phasevec, const dvector& slicevec, const dvector& centervec);



/**
  * Returns a vector pointing to the center of the slicepack/voxel.
  */
  dvector get_center() const;

/**
  * Returns the rotation matrix to convert the gradient strengths in the
  * sequence coordinate system (read,phase,slice) to the laboratory system (x,y,z).
  * If 'transpose' is set to true, the transposed matrix is returned.
  */
  RotMatrix get_gradrotmatrix(bool transpose=false) const;


/**
  * Coordinate transformation from sequence coordinate system (read,phase,slice)
  * of slice 'slice' to laboratory system (x,y,z), or vice versa if
  * 'inverse' is set to 'true'.
  */
  dvector transform(const dvector& rpsvec, bool inverse=false) const;

/**
  * Transpose geometry in-plane, 'reverse_read' and 'reverse_phase'
  * can be used to reverse read/phase direction before transposing.
  */
  Geometry& transpose_inplane(bool reverse_read=false, bool reverse_phase=false);


/**
  * Resets the slicepacks geometry to its inintial state
  */
  Geometry& reset();

/**
  * Updates internal relations of the geometry parameters
  */
  Geometry& update();



 private:

  dvector get_readVector_inplane() const;
  dvector get_phaseVector_inplane() const;

  void append_all_members();

  double deg2rad(double degree) const {return degree/180.0*PII;}

  JDXenum Mode;

  JDXdouble FOVread;
  JDXdouble offsetRead;
  JDXdouble FOVphase;
  JDXdouble offsetPhase;
  JDXdouble FOVslice;
  JDXdouble offsetSlice;

  JDXdouble heightAngle;
  JDXdouble azimutAngle;
  JDXdouble inplaneAngle;
  JDXbool   reverseSlice;

  JDXint    nSlices;
  JDXdouble sliceDistance;
  JDXdouble sliceThickness;
  JDXintArr sliceOrder;

  JDXaction Reset;
  JDXaction Transpose;


  // cache to speed up repetitive calls to transform()
  mutable bool cache_up2date;
  mutable bool inv_trans_cache;
  mutable double rotmat_cache[3][3];
  mutable double offset_cache[3];
};


/** @}
  */


#endif