ImageAttributes.h

/*
 * Medical Image Registration ToolKit (MIRTK)
 *
 * Copyright 2008-2015 Imperial College London
 * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel
 * Copyright 2013-2015 Andreas Schuh
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef MIRTK_ImageAttributes_H
#define MIRTK_ImageAttributes_H

#include "mirtk/Indent.h"
#include "mirtk/Matrix.h"
#include "mirtk/Point.h"
#include "mirtk/Array.h"


namespace mirtk {


/**
 * Class which defines the attributes of the imaging geometry
 */
struct ImageAttributes
{
  // ---------------------------------------------------------------------------
  // Attributes
  int _x; ///< Image x-dimension (in voxels)
  int _y; ///< Image y-dimension (in voxels)
  int _z; ///< Image z-dimension (in voxels)
  int _t; ///< Image t-dimension (in voxels)

  double _dx; ///< Voxel x-dimensions (in mm)
  double _dy; ///< Voxel y-dimensions (in mm)
  double _dz; ///< Voxel z-dimensions (in mm)
  double _dt; ///< Voxel t-dimensions (in ms)

  double _xorigin; ///< Image x-origin (in mm)
  double _yorigin; ///< Image y-origin (in mm)
  double _zorigin; ///< Image z-origin (in mm)
  double _torigin; ///< Image t-origin (in ms)

  double _xaxis[3]; ///< Direction of x-axis
  double _yaxis[3]; ///< Direction of y-axis
  double _zaxis[3]; ///< Direction of z-axis

  Matrix _smat; ///< Affine transformation matrix

  const Matrix *_w2i;  ///< Pointer to pre-computed world to image matrix (cf. BaseImage::_matW2I)
  const Matrix *_i2w;  ///< Pointer to pre-computed image to world matrix (cf. BaseImage::_matI2W)

  /// Number of lattice points in image domain
  int NumberOfLatticePoints() const;

  /// Number of spatial lattice points
  int NumberOfSpatialPoints() const;

  /// Number of lattice points
  /// (i.e., NumberOfLatticePoints if _dt != 0, otherwise NumberOfSpatialPoints)
  int NumberOfPoints() const;

  /// Get number of lattice points in x dimension
  int X() const;

  /// Get number of lattice points in y dimension
  int Y() const;

  /// Get number of lattice points in z dimension
  int Z() const;

  /// Get number of lattice points in t dimension
  int T() const;

  /// Get number of lattice points in i-th dimension
  int N(int) const;

  /// Get spacing of lattice points in i-th dimension
  double Spacing(int) const;

  /// Get spacing of lattice points in x dimension
  double XSpacing() const;

  /// Get spacing of lattice points in y dimension
  double YSpacing() const;

  /// Get spacing of lattice points in z dimension
  double ZSpacing() const;

  /// Get spacing of lattice points in t dimension
  double TSpacing() const;

  /// Get voxel size/spacing of lattice points in x dimension
  double XSize() const;

  /// Get voxel size/spacing of lattice points in y dimension
  double YSize() const;

  /// Get voxel size/spacing of lattice points in z dimension
  double ZSize() const;

  /// Get voxel size/spacing of lattice points in t dimension
  double TSize() const;

  // ---------------------------------------------------------------------------
  // Construction/Destruction

  /// Constructor
  ImageAttributes();

  /// Constructor
  ImageAttributes(int, int, double = 1.0, double = 1.0);

  /// Constructor
  ImageAttributes(int, int, int, double = 1.0, double = 1.0, double = 1.0);

  /// Constructor
  ImageAttributes(int, int, int, int, double = 1.0, double = 1.0, double = 1.0, double = 1.0);

  /// Copy constructor
  ImageAttributes(const ImageAttributes &);

  /// Copy operator
  ImageAttributes& operator= (const ImageAttributes &);

  /// Check if attributes are valid
  operator bool() const;

  // ---------------------------------------------------------------------------
  // Comparison

  /// Whether given lattice is fully contained by this image lattice
  bool ContainsInSpace(const ImageAttributes &attr) const;

  /// Whether spatial attributes are equal
  bool EqualInSpace(const ImageAttributes &attr) const;

  /// Whether temporal attributes are equal
  bool EqualInTime(const ImageAttributes &attr) const;

  /// Equality operator
  bool operator==(const ImageAttributes &attr) const;

  /// Inequality operator
  bool operator!=(const ImageAttributes &attr) const;

  // ---------------------------------------------------------------------------
  // Coordinate conversion

  /// Get Index from Lattice
  int LatticeToIndex(int, int, int = 0, int = 0) const;

  /// Get Index from Lattice
  void IndexToLattice(int, int *, int *, int * = NULL, int * = NULL) const;

  /// Get Index from Lattice
  void IndexToLattice(int, int &, int &) const;

  /// Get Index from Lattice
  void IndexToLattice(int, int &, int &, int &) const;

  /// Get Index from Lattice
  void IndexToLattice(int, int &, int &, int &, int &) const;

  /// Get world coordinates (in mm) of lattice point
  void IndexToWorld(int, double &, double &) const;

  /// Get world coordinates (in mm) of lattice point
  void IndexToWorld(int, double &, double &, double &) const;

  /// Get world coordinates (in mm) of lattice point
  void IndexToWorld(int, Point &) const;

  /// Get world coordinates (in mm) of lattice point
  Point IndexToWorld(int) const;

  /// Convert lattice to world coordinate
  void LatticeToWorld(Point &) const;

  /// Convert lattice to world coordinate
  void LatticeToWorld(double &, double &, double &) const;

  /// Compute spatial world coordinates of lattice points (_x * _y * _z)
  void LatticeToWorld(double *, double *, double *) const;

  /// Compute world coordinates of all lattice points (_x * _y * _z * _t)
  void LatticeToWorld(double *, double *, double *, double *) const;

  /// Convert world to lattice coordinate
  void WorldToLattice(double &, double &, double &) const;

  /// Convert lattice to world coordinate
  void WorldToLattice(Point &) const;

  /// Convert lattice to time coordinate
  double LatticeToTime(double) const;

  /// Convert time to lattice coordinate
  double TimeToLattice(double) const;

  /// Put affine world coordinate transformation which is applied
  /// after the image to world coordinate transformation derived from the
  /// imaging geometry when mapping voxel indices to world coordinates.
  /// This transformation can be the inverse of the affine transformation
  /// obtained by an affine registration with this image as source.
  ///
  /// \param[in] m     Homogeneous transformation matrix.
  /// \param[in] apply Whether to apply the translation, rotation, and
  ///                  scaling directly to the image attributes and store
  ///                  only the shearing (if any) as additional transformation.
  void PutAffineMatrix(const Matrix &m, bool apply = false);

  /// Return transformation matrix for lattice to world coordinates
  Matrix GetLatticeToWorldMatrix() const;

  /// Return transformation matrix for world to lattice coordinates
  Matrix GetWorldToLatticeMatrix() const;

  /// Return orientation part of lattice to world coordinate transformation
  Matrix GetLatticeToWorldOrientation() const;

  /// Return orientation part of lattice to world coordinate transformation
  Matrix GetWorldToLatticeOrientation() const;

  /// Alias for GetLatticeToWorldMatrix
  inline Matrix GetImageToWorldMatrix()      const { return GetLatticeToWorldMatrix(); }
  /// Alias for GetWorldToLatticeMatrix
  inline Matrix GetWorldToImageMatrix()      const { return GetWorldToLatticeMatrix(); }
  /// Alias for GetLatticeToWorldOrientation
  inline Matrix GetImageToWorldOrientation() const { return GetLatticeToWorldOrientation(); }
  /// Alias for GetWorldToLatticeOrientation
  inline Matrix GetWorldToImageOrientation() const { return GetWorldToLatticeOrientation(); }

  // ---------------------------------------------------------------------------
  // Voxel index checks

  /// Whether voxel index is within finite image domain
  bool IsInside(int) const;

  /// Whether voxel indices are within finite 2D image domain
  bool IsInside(int, int) const;

  /// Whether voxel indices are within finite 3D image domain
  bool IsInside(int, int, int) const;

  /// Whether voxel indices are within finite 4D image domain
  bool IsInside(int, int, int, int) const;

  /// Whether voxel is index is outside finite image domain
  bool IsOutside(int) const;

  /// Whether voxel indices are outside finite 4D image domain
  bool IsOutside(int, int) const;

  /// Whether voxel indices are outside finite 4D image domain
  bool IsOutside(int, int, int) const;

  /// Whether voxel indices are outside finite 4D image domain
  bool IsOutside(int, int, int, int) const;

  /// Whether voxel index is at boundary of finite image domain
  bool IsBoundary(int) const;

  /// Whether voxel indices are at boundary of finite 2D image domain
  bool IsBoundary(int, int) const;

  /// Whether voxel indices are at boundary of finite 3D image domain
  bool IsBoundary(int, int, int) const;

  /// Whether voxel indices are at boundary of finite 4D image domain
  bool IsBoundary(int, int, int, int) const;

  // ---------------------------------------------------------------------------
  // Output

  /// Image slice area in world space
  double Area() const;

  /// Image volume in world space
  double Volume() const;

  /// Amount of space occupied by the image in n-D world space (excluding time dimension)
  double Space() const;

  /// Print attributes
  void Print(ostream &, Indent = 0) const;

  /// Print attributes
  void Print(Indent = 0) const;

};

////////////////////////////////////////////////////////////////////////////////
// Inline definitions
////////////////////////////////////////////////////////////////////////////////

// -----------------------------------------------------------------------------
inline int ImageAttributes::X() const
{
  return _x;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::Y() const
{
  return _y;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::Z() const
{
  return _z;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::T() const
{
  return _t;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::N(int i) const
{
  if (i == 0) return _x;
  if (i == 1) return _y;
  if (i == 2) return _z;
  if (i == 3) return _t;
  return 0;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::Spacing(int i) const
{
  if (i == 0) return _dx;
  if (i == 1) return _dy;
  if (i == 2) return _dz;
  if (i == 3) return _dt;
  return 0.;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::XSpacing() const
{
  return _dx;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::YSpacing() const
{
  return _dy;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::ZSpacing() const
{
  return _dz;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::TSpacing() const
{
  return _dt;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::XSize() const
{
  return _dx;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::YSize() const
{
  return _dy;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::ZSize() const
{
  return _dz;
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::TSize() const
{
  return _dt;
}

// -----------------------------------------------------------------------------
inline ImageAttributes::operator bool() const
{
  // Note: _dz may be zero for 2D images, _dt may even be negative!
  return _x > 0 && _y > 0 && _z > 0 && _t > 0 && _dx > .0 && _dy > .0 && _dz >= .0;
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::operator==(const ImageAttributes &attr) const
{
  return EqualInSpace(attr) && EqualInTime(attr);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::operator!=(const ImageAttributes &attr) const
{
  return !(*this == attr);
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::LatticeToIndex(int i, int j, int k, int l) const
{
  return i + _x * (j + _y * (k + _z * l));
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToLattice(int index, int *i, int *j, int *k, int *l) const
{
  int n = _x * _y * _z;
  if (l) *l = index / n;
  index = index % n;
  n = _x * _y;
  if (k) *k = index / n;
  *j = index % n / _x;
  *i = index % n % _x;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::NumberOfLatticePoints() const
{
  return _x * _y * _z * _t;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::NumberOfSpatialPoints() const
{
  return _x * _y * _z;
}

// -----------------------------------------------------------------------------
inline int ImageAttributes::NumberOfPoints() const
{
  return _dt ? NumberOfLatticePoints() : NumberOfSpatialPoints();
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToLattice(int index, int &i, int &j) const
{
  IndexToLattice(index, &i, &j);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToLattice(int index, int &i, int &j, int &k) const
{
  IndexToLattice(index, &i, &j, &k);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToLattice(int index, int &i, int &j, int &k, int &l) const
{
  IndexToLattice(index, &i, &j, &k, &l);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToWorld(int idx, double &x, double &y) const
{
  int i, j, k;
  IndexToLattice(idx, i, j, k);
  x = i, y = j;
  double z = k;
  LatticeToWorld(x, y, z);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToWorld(int idx, double &x, double &y, double &z) const
{
  int i, j, k;
  IndexToLattice(idx, i, j, k);
  x = i, y = j, z = k;
  LatticeToWorld(x, y, z);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::IndexToWorld(int idx, Point &p) const
{
  IndexToWorld(idx, p._x, p._y, p._z);
}

// -----------------------------------------------------------------------------
inline Point ImageAttributes::IndexToWorld(int idx) const
{
  Point p;
  IndexToWorld(idx, p);
  return p;
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::LatticeToWorld(double &x, double &y, double &z) const
{
  Transform(_i2w ? *_i2w : GetImageToWorldMatrix(), x, y, z);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::LatticeToWorld(Point &p) const
{
  LatticeToWorld(p._x, p._y, p._z);
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::LatticeToTime(double t) const
{
  return _torigin + t * _dt;
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::WorldToLattice(double &x, double &y, double &z) const
{
  Transform(_w2i ? *_w2i : GetWorldToImageMatrix(), x, y, z);
}

// -----------------------------------------------------------------------------
inline void ImageAttributes::WorldToLattice(Point &p) const
{
  WorldToLattice(p._x, p._y, p._z);
}

// -----------------------------------------------------------------------------
inline double ImageAttributes::TimeToLattice(double t) const
{
  return (_dt ? ((t - _torigin) / _dt) : .0);
}

// -----------------------------------------------------------------------------
inline Matrix ImageAttributes::GetLatticeToWorldOrientation() const
{
  Matrix R(3, 3);
  R(0, 0) = _xaxis[0];
  R(1, 0) = _xaxis[1];
  R(2, 0) = _xaxis[2];
  R(0, 1) = _yaxis[0];
  R(1, 1) = _yaxis[1];
  R(2, 1) = _yaxis[2];
  R(0, 2) = _zaxis[0];
  R(1, 2) = _zaxis[1];
  R(2, 2) = _zaxis[2];
  return R;
}

// -----------------------------------------------------------------------------
inline Matrix ImageAttributes::GetWorldToLatticeOrientation() const
{
  Matrix R(3, 3);
  R(0, 0) = _xaxis[0];
  R(0, 1) = _xaxis[1];
  R(0, 2) = _xaxis[2];
  R(1, 0) = _yaxis[0];
  R(1, 1) = _yaxis[1];
  R(1, 2) = _yaxis[2];
  R(2, 0) = _zaxis[0];
  R(2, 1) = _zaxis[1];
  R(2, 2) = _zaxis[2];
  return R;
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsInside(int idx) const
{
  return (0 <= idx && idx < _x * _y * _z * _t);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsInside(int i, int j) const
{
  return (0 <= i && i < _x) && (0 <= j && j < _y);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsInside(int i, int j, int k) const
{
  return IsInside(i, j) && (0 <= k && k < _z);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsInside(int i, int j, int k, int l) const
{
  return IsInside(i, j, k) && (0 <= l && l < _t);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsOutside(int idx) const
{
  return !IsInside(idx);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsOutside(int i, int j) const
{
  return !IsInside(i, j);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsOutside(int i, int j, int k) const
{
  return !IsInside(i, j, k);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsOutside(int i, int j, int k, int l) const
{
  return !IsInside(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsBoundary(int i, int j) const
{
  return (i == 0 || i == _x - 1) || (j == 0 || j == _y - 1);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsBoundary(int i, int j, int k) const
{
  return IsBoundary(i, j) || (k == 0 || k == _z - 1);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsBoundary(int i, int j, int k, int l) const
{
  return IsBoundary(i, j, k) || (l == 0 || l == _t - 1);
}

// -----------------------------------------------------------------------------
inline bool ImageAttributes::IsBoundary(int idx) const
{
  if (!IsInside(idx)) return false;
  int i, j, k, l;
  IndexToLattice(idx, i, j, k, l);
  if (_t == 1) {
    if (_z == 1) return IsBoundary(i, j);
    else         return IsBoundary(i, j, k);
  }
  return IsBoundary(i, j, k, l);
}

////////////////////////////////////////////////////////////////////////////////
// Image domain helpers
////////////////////////////////////////////////////////////////////////////////

/// Ensure that image domain has orthogonal basis vectors, i.e., no additional
/// affine transformation which contains any shearing
///
/// \returns Image grid which fully contains the input image domain but without
///          additional affine transformation (_smat is identity matrix).
///
/// \remarks The returned image domain need not be an axis-aligned bounding box!
///          When the Geometric Tools Engine (GTEngine) library is available,
///          the minimum-volume bounding box is returned.
ImageAttributes OrthogonalFieldOfView(const ImageAttributes &);

/// This method implements a method to find a common image grid which fully
/// contains all given image grids in world space
///
/// The voxel size of the resulting grid corresponds to the average voxel size
/// of the input images. The orientation and origin are chosen such that the
/// resulting image domain overlaps all input images in world space with a near
/// to minimal image volume. The current implementation only uses an heuristic
/// approach to find such minimal oriented bounding box of the input domains.
/// When the Geometric Tools Engine (GTEngine) library is available, the
/// minimum-volume bounding box is returned.
///
/// Additional information regarding algorithms to find the optimal
/// minimum-volume oriented bounding box (OBB) can be found at:
/// - http://stackoverflow.com/questions/7282805/where-can-i-find-a-c-c-implementation-of-the-minimum-bound-box-algorithm
/// - http://www.geometrictools.com/LibMathematics/Containment/Containment.html
/// - http://link.springer.com/article/10.1007%2FBF00991005?LI=true
///
/// \note The final overall image grid computed by this function must be
///       independent of the order in which the input attributes are given
///       in the input list. Otherwise an inverse consistent registration
///       might still depend on the order of the input images.
ImageAttributes OverallFieldOfView(const Array<ImageAttributes> &);


} // namespace mirtk

#endif // MIRTK_ImageAttributes_H