BaseImage.h

/*
 * Medical Image Registration ToolKit (MIRTK)
 *
 * Copyright 2008-2017 Imperial College London
 * Copyright 2008-2013 Daniel Rueckert, Julia Schnabel
 * Copyright 2013-2017 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_BaseImage_H
#define MIRTK_BaseImage_H

#include "mirtk/Object.h"
#include "mirtk/Voxel.h"
#include "mirtk/Point.h"
#include "mirtk/PointSet.h"
#include "mirtk/Vector.h"
#include "mirtk/Vector3.h"
#include "mirtk/Matrix.h"
#include "mirtk/ImageAttributes.h"


// When VTK support is enabled during the build of the Image module,
// the vtable of the image classes contains additional entries compared to
// a build without VTK! Therefore, MIRTK_Image_WITH_VTK must reflect whether
// the interface includes those VTK support functions or not.
#include "mirtk/ImageConfig.h"
#if MIRTK_Image_WITH_VTK
class vtkStructuredPoints;
#endif // MIRTK_Image_WITH_VTK


namespace mirtk {


////////////////////////////////////////////////////////////////////////////////
// Forward declaration of used specialized image types
////////////////////////////////////////////////////////////////////////////////

// Forward declaration of generic image type
template <class VoxelType> class GenericImage;

/// Binary image as used for masks (0: off, otherwise: on)
typedef GenericImage<BinaryPixel> BinaryImage;

/// FIXME: Use double3 as voxel type instead
typedef GenericImage<double> WorldCoordsImage;

////////////////////////////////////////////////////////////////////////////////
// Definition of base image class
////////////////////////////////////////////////////////////////////////////////

/**
 * Abstract generic class for 2D or 3D images
 *
 * This abstract base class provides generic support for 2D and 3D image
 * classes. It provides functions for accessing image dimension and origin
 * as well as voxel dimensions. It also provides functions for conversion
 * between image and world coordinates.
 */
class BaseImage : public Object
{
  mirtkAbstractMacro(BaseImage);

  // ---------------------------------------------------------------------------
  // Types

public:

  /// Orientation codes (same as defined in nifti1_io.h, e.g., NIFTI_L2R)
  enum OrientationCode
  {
    Orientation_Unknown = 0,
    L2R = 1, ///< Left to Right
    R2L = 2, ///< Right to Left
    P2A = 3, ///< Posterior to Anterior
    A2P = 4, ///< Anterior to Posterior
    I2S = 5, ///< Inferior to Superior
    S2I = 6  ///< Superior to Inferior
  };

  // ---------------------------------------------------------------------------
  // Attributes

protected:

  /// Image attributes
  ImageAttributes _attr;

  /// Total number of voxels
  int _NumberOfVoxels;

  /// Transformation matrix from image coordinates to (transformed) world coordinates
  Matrix _matI2W;

  /// Transformation matrix from (transformed) world coordinates to image coordinates
  Matrix _matW2I;

  /// Foreground mask
  BinaryImage *_mask;

  /// Whether this image owns the mask image
  ///
  /// Set to a value above 1 when mask was generated by InitializeMask,
  /// in which case the mask will be deleted by the n-th ClearMask call,
  /// where n is the number of preceeding (nested) InitializeMask calls.
  /// A value of 1 denotes a mask image set by the owner which this image
  /// should delete upon destruction.
  int _maskOwner;

  /// Background value - may also be NaN for floating point images
  double _bg;

  /// Whether a background value was set
  bool _bgSet;

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

  /// Default constructor
  BaseImage();

  /// Constructor
  BaseImage(const ImageAttributes &, int = -1);

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

public:

  /// Destructor
  virtual ~BaseImage();

  // ---------------------------------------------------------------------------
  // Initialization

  /// Read file and construct image
  static BaseImage *New(const char *);
  
  /// Construct image copy of same type
  static BaseImage *New(const BaseImage *);

  /// Construct image with given voxel type
  static BaseImage *New(int);

  /// Create copy of this image
  virtual BaseImage *Copy() const;

protected:

  /// Update coordinate transformation
  void UpdateMatrix();

  /// Puts attributes of image
  void PutAttributes(const ImageAttributes &);

public:

  /// Initialize image
  virtual void Initialize(const ImageAttributes &, int = -1) = 0;

  /// Assignment operator
  BaseImage &operator =(const BaseImage &);

  /// Clear image
  virtual void Clear() = 0;

  // ---------------------------------------------------------------------------
  // Lattice

  /// Gets the image attributes
  const ImageAttributes &Attributes() const;

  /// Returns the total number of voxels
  int NumberOfVoxels() const;

  /// Returns the total number of spatial voxels
  int NumberOfSpatialVoxels() const;

  /// Returns the number of voxels in the x-direction
  int X() const;
  
  /// Returns the number of voxels in the y-direction
  int Y() const;
  
  /// Returns the number of voxels in the z-direction
  int Z() const;
  
  /// Returns the number of voxels in the t-direction
  int T() const;

  /// Returns the number of vector components (i.e., 1 for scalar images)
  virtual int N() const;

  /// Returns the number of voxels in the x-direction
  int GetX() const;

  /// Returns the number of voxels in the y-direction
  int GetY() const;

  /// Returns the number of voxels in the z-direction
  int GetZ() const;

  /// Returns the number of voxels in the t-direction
  int GetT() const;

  /// Returns the size of a voxel in the x-direction
  double XSize() const;

  /// Returns the size of a voxel in the y-direction
  double YSize() const;

  /// Returns the size of a voxel in the z-direction
  double ZSize() const;

  /// Returns the size of a voxel in the t-direction
  double TSize() const;

  /// Returns the size of a voxel in the x-direction
  double GetXSize() const;

  /// Returns the size of a voxel in the y-direction
  double GetYSize() const;

  /// Returns the size of a voxel in the z-direction
  double GetZSize() const;

  /// Returns the size of a voxel in the t-direction
  double GetTSize() const;

  /// Set temporal voxel size, i.e., to zero for vector field and non-zero for temporal sequence
  void PutTSize(double);

  /// Voxel dimensions get access
  void GetPixelSize(double &, double &) const;

  /// Voxel dimensions get access
  void GetPixelSize(double &, double &, double &) const;

  /// Voxel dimensions get access
  void GetPixelSize(double &, double &, double &, double &) const;

  /// Voxel dimensions put access
  void PutPixelSize(double, double, double);

  /// Voxel dimensions put access
  void PutPixelSize(double, double, double, double);

  /// Image origin get access
  Point GetOrigin() const;

  /// Image origin get access
  void GetOrigin(double &, double &, double &) const;

  /// Image origin get access
  void GetOrigin(double &, double &, double &, double &) const;

  /// Image origin put access
  void PutOrigin(const Point &);

  /// Image origin put access
  void PutOrigin(double, double, double);

  /// Image origin put access
  void PutOrigin(double, double, double, double);

  /// Put temporal origin
  void PutTOrigin(double);
  
  /// Get temporal origin
  double GetTOrigin() const;

  /// Put image x- and y-axis and z-axis
  void PutOrientation(double *, double *, double * = NULL);

  /// Get image x- and y-axis and z-axis
  void GetOrientation(double *, double *, double * = NULL) const;

  /// Get orientation of axis relative to patient
  void Orientation(OrientationCode &, OrientationCode &, OrientationCode &) 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.
  void PutAffineMatrix(const Matrix &, bool = false);

  /// Reset 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.
  void ResetAffineMatrix();

  /// Get 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
  const Matrix &GetAffineMatrix() const;

  /// Function to convert pixel to index
  int VoxelToIndex(int, int, int = 0, int = 0) const;
  
  /// Function to convert index to pixel coordinates
  void IndexToVoxel(int, int &, int &) const;
  
  /// Function to convert index to pixel coordinates
  void IndexToVoxel(int, int &, int &, int &) const;
  
  /// Function to convert index to pixel coordinates
  void IndexToVoxel(int, int &, int &, int &, int &) const;

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

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

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

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

  /// Image to world coordinate conversion with two doubles
  void ImageToWorld(double &, double &) const;

  /// Image to world coordinate conversion with three doubles
  void ImageToWorld(double &, double &, double &) const;

  /// Image to world coordinate conversion with a given point
  void ImageToWorld(Point &) const;

  /// Convert vector w.r.t. image axes to vector w.r.t. world axes
  void ImageToWorld(Vector3 &) const;

  /// Image to world coordinate map for image domain
  ///
  /// \note Stores x, y, and z components as vector image,
  ///       i.e., in t dimension of image.
  ///
  /// \param[out] i2w Vector image with world coordinates of voxels.
  ///                 Need not be initialized before. No reallocation
  ///                 takes place if the image has already the right size.
  /// \param[in] _3D  If \c false, only the x and y coordinates are
  ///                 stored if the image is two dimensional. Otherwise,
  ///                 all three x, y, and z components are stored.
  void ImageToWorld(WorldCoordsImage &i2w, bool _3D = true) const;

  /// Image to world coordinate map for image domain
  ///
  /// \note Stores x, y, and z components at consecutive memory
  ///       locations for faster access.
  ///
  /// \param[out] i2w Memory for pre-computed world coordinates.
  /// \param[in] _3D  If \c false, only the x and y coordinates are
  ///                 stored if the image is two dimensional. Otherwise,
  ///                 all three x, y, and z components are stored.
  void ImageToWorld(double *i2w, bool _3D = true) const;

  /// Adds world coordinates of each voxel to given point set
  void ImageToWorld(PointSet &) const;

  /// World to image coordinate conversion with two doubles
  void WorldToImage(double &, double &) const;

  /// World to image coordinate conversion with three doubles
  void WorldToImage(double &, double &, double &) const;
  
  /// World to image coordinate conversion with a given point
  void WorldToImage(Point &) const;

  /// Convert vector w.r.t. world axes to vector w.r.t. image axes
  void WorldToImage(Vector3 &) const;

  /// Return transformation matrix for image to world coordinates
  const Matrix &GetImageToWorldMatrix() const;
  
  /// Return transformation matrix for world to image coordinates
  const Matrix &GetWorldToImageMatrix() const;

  /// Image to time coordinate conversion
  double ImageToTime(double) const;

  /// Time to image coordinate conversion
  double TimeToImage(double) const;

  /// Checks if this image shares the same spatial attributes with another image
  bool HasSpatialAttributesOf(const BaseImage *) const;

  /// Returns true if point is within the field of view of image
  bool IsInFOV(double, double, double);
  
  /// Whether image is uninitialized
  bool IsEmpty() const;

  // ---------------------------------------------------------------------------
  // Type independent access to image data

  /// Function for pixel get access as double
  virtual double GetAsDouble(int) const;
  
  /// Function for pixel get access as double
  virtual double GetAsDouble(int, int, int = 0, int = 0) const = 0;
  
  /// Function for pixel put access
  virtual void PutAsDouble(int, double);
  
  /// Function for pixel put access
  virtual void PutAsDouble(int, int, double);
  
  /// Function for pixel put access
  virtual void PutAsDouble(int, int, int, double);

  /// Function for pixel put access
  virtual void PutAsDouble(int, int, int, int, double) = 0;

  /// Function for pixel get access as double
  virtual void GetAsVector(Vector &, int) const;

  /// Function for pixel get access as double
  virtual void GetAsVector(Vector &, int, int, int = 0, int = 0) const = 0;

  /// Function for pixel get access as double
  virtual Vector GetAsVector(int) const;

  /// Function for pixel get access as double
  virtual Vector GetAsVector(int, int, int = 0, int = 0) const;

  /// Function for pixel put access
  virtual void PutAsVector(int, const Vector &);

  /// Function for pixel put access
  virtual void PutAsVector(int, int, const Vector &);

  /// Function for pixel put access
  virtual void PutAsVector(int, int, int, const Vector &);

  /// Function for pixel put access
  virtual void PutAsVector(int, int, int, int, const Vector &) = 0;

  // ---------------------------------------------------------------------------
  // Foreground region

  /// Set foreground mask
  void PutMask(BinaryImage *, bool = false);

  /// Get foreground mask (optionally, take over ownership)
  BinaryImage *GetMask(bool = false);

  /// Get foreground mask
  const BinaryImage *GetMask() const;

  /// Whether this image has a foreground mask which has been
  /// either set via PutMask or created upon InitializeMask
  bool HasMask() const;

  /// Whether this image has a foreground mask
  /// either set via PutMask with ownership transver or created upon InitializeMask
  bool OwnsMask() const;

  /// Initialize mask if not done yet or none is set
  void InitializeMask(int t = -1, bool = false);

  /// Clear mask upon n-th call after n preceeding InitializeMask calls
  void ClearMask(bool = false);

  /// Put background value
  void PutBackgroundValueAsDouble(double);

  /// Put background value
  virtual void PutBackgroundValueAsDouble(double, bool);

  /// Get background value
  double GetBackgroundValueAsDouble() const;

  /// Change background value
  ///
  /// When a background value has been set before, this function replaces all image
  /// values equal to this previous background value by the new background value.
  /// Otherwise, it is equivalent to PutBackgroundValueAsDouble.
  void ResetBackgroundValueAsDouble(double);

  /// Clear background value
  void ClearBackgroundValue();

  /// Whether a background value has been set
  bool HasBackgroundValue() const;

  /// Whether voxel is within foreground without index-out-of-bounds check
  bool IsForeground(int) const;

  /// Whether voxel is within foreground without index-out-of-bounds check
  bool IsForeground(int, int, int = 0, int = 0) const;

  /// Whether voxel is within background without index-out-of-bounds check
  bool IsBackground(int) const;

  /// Whether voxel is within background without index-out-of-bounds check
  bool IsBackground(int, int, int = 0, int = 0) const;

  /// 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;

  /// Whether voxel is index is within finite image domain and part of foreground region
  bool IsInsideForeground(int) const;

  /// Whether voxel indices are within finite image domain and part of foreground region
  bool IsInsideForeground(int, int, int = 0, int = 0) const;

  /// Whether voxel is index is outside finite image domain or part of background region
  bool IsOutsideForeground(int) const;

  /// Whether voxel indices are outside finite image domain or part of background region
  bool IsOutsideForeground(int, int, int = 0, int = 0) const;

  /// Whether all voxels within a 2D bounding region are inside foreground region
  bool IsBoundingBoxInsideForeground(int, int, int, int) const;

  /// Whether all voxels within a 3D bounding region are inside foreground region
  bool IsBoundingBoxInsideForeground(int, int, int, int, int, int) const;

  /// Whether all voxels within a 4D bounding region are inside foreground region
  bool IsBoundingBoxInsideForeground(int, int, int, int, int, int, int, int) const;

  /// Whether at least one neighboring voxel is outside the finite foreground region
  bool IsNextToBackground(int) const;

  /// Whether at least one neighboring voxel is outside the finite foreground region
  bool IsNextToBackground(int, int, int = 0, int = 0) const;

  /// Whether at least one neighboring voxel is inside the finite foreground region
  bool IsNextToForeground(int) const;

  /// Whether at least one neighboring voxel is inside the finite foreground region
  bool IsNextToForeground(int, int, int = 0, int = 0) const;

  /// Whether any voxel is within background
  bool HasBackground() const;

  /// Get 2D bounding box of image foreground
  void BoundingBox(int &, int &, int &, int &) const;

  /// Get 3D bounding box of image foreground
  void BoundingBox(int &, int &, int &, int &, int &, int &) const;

  /// Get 3D+t bounding box of image foreground
  void BoundingBox(int &, int &, int &, int &, int &, int &, int &, int &) const;

  /// Determine center of mass of image foreground
  ///
  /// \param[out] center Centroid of image foreground.
  ///
  /// \returns Number of foreground intensities encountered or zero if the image does
  ///          not contain any foreground, in which case the \p center is invalid.
  int CenterOfForeground(Point &center) const;

  /// Determine center of mass of image foreground
  ///
  /// \param[out] center  Centroid of image foreground.
  /// \param[in]  padding Padding value. Voxels with a value less or equal are considered background.
  ///                     When this value is NaN, all non-NaN image values are considered.
  ///
  /// \returns Number of foreground intensities encountered or zero if the image does
  ///          not contain any foreground, in which case the \p center is invalid.
  int CenterOfForeground(Point &center, double padding) const;

  /// Get cropped image foreground bounding region
  ///
  /// \param[in] i1         Lower voxel index in x direction.
  /// \param[in] j1         Lower voxel index in y direction.
  /// \param[in] k1         Lower voxel index in z direction.
  /// \param[in] i2         Upper voxel index in x direction.
  /// \param[in] j2         Upper voxel index in y direction.
  /// \param[in] k2         Upper voxel index in z direction.
  /// \param[in] orthogonal Whether to orthogonalize image axes. Useful only in case
  ///                       of input images where a previous affine (12 DoFs) alignment
  ///                       has been applied to the attributes.
  ///
  /// \returns Attributes of foreground region.
  ImageAttributes ForegroundDomain(int i1, int j1, int k1,
                                   int i2, int j2, int k2,
                                   bool orthogonal = true) const;

  /// Determine minimal axes-aligned foreground bounding region
  ///
  /// \param[in] orthogonal Whether to orthogonalize image axes. Useful only in case
  ///                       of input images where a previous affine (12 DoFs) alignment
  ///                       has been applied to the attributes.
  ///
  /// \returns Attributes of foreground region.
  ImageAttributes ForegroundDomain(bool orthogonal = true) const;

  /// Determine minimal axes-aligned foreground bounding region
  ///
  /// \param[in] padding    Padding value. Voxels with a value less or equal the specified
  ///                       value are considered background. When the padding value is NaN,
  ///                       all NaN image values are considered background.
  /// \param[in] orthogonal Whether to orthogonalize image axes. Useful only in case
  ///                       of input images where a previous affine (12 DoFs) alignment
  ///                       has been applied to the attributes.
  ///
  /// \returns Attributes of foreground region.
  ImageAttributes ForegroundDomain(double padding, bool orthogonal = true) const;

  // ---------------------------------------------------------------------------
  // Region-of-interest extraction

  /// Get image consisting of specified 2D slice
  virtual void GetRegion(BaseImage *&, int, int) const = 0;

  /// Get image consisting of specified 3D subregion
  virtual void GetRegion(BaseImage *&, int, int, int, int, int, int) const = 0;

  /// Get image consisting of specified 4D subregion
  virtual void GetRegion(BaseImage *&, int, int, int, int, int, int, int, int) const = 0;

  /// Get time instance (i.e., frame) or channel of image
  virtual void GetFrame(BaseImage *&, int, int = -1) const = 0;

  // ---------------------------------------------------------------------------
  // Common image statistics

  /// Minimum and maximum pixel values get accessor
  virtual void GetMinMaxAsDouble(double &, double &) const;

  /// Minimum and maximum pixel values put accessor
  virtual void PutMinMaxAsDouble(double, double);
  
  // ---------------------------------------------------------------------------
  // Access to raw image data

  /// Function for pixel access via pointers
  virtual void *GetDataPointer(int = 0) = 0;

  /// Function for pixel access via pointers
  virtual const void *GetDataPointer(int = 0) const = 0;

  /// Function for pixel access via pointers
  virtual void *GetDataPointer(int, int, int = 0, int = 0) = 0;
  
  /// Function for pixel access via pointers
  virtual const void *GetDataPointer(int, int, int = 0, int = 0) const = 0;

  /// Function which returns pixel scalar type
  virtual int GetDataType() const = 0;

  /// Function which returns size of pixel scalar type in bytes
  virtual int GetDataTypeSize() const = 0;

  /// Function which returns the minimum value the pixel can hold without overflowing
  virtual double GetDataTypeMin() const = 0;

  /// Function which returns the minimum value the pixel can hold without overflowing
  virtual double GetDataTypeMax() const = 0;

  // ---------------------------------------------------------------------------
  // Common image manipulations
  
  virtual void ReflectX(bool modify_axes = false) = 0;  ///< Reflect image along x
  virtual void ReflectY(bool modify_axes = false) = 0;  ///< Reflect image along y
  virtual void ReflectZ(bool modify_axes = false) = 0;  ///< Reflect image along z
  virtual void ReflectT(bool modify_axes = false) = 0;  ///< Reflect image along t

  virtual void FlipXY(bool modify_origin = false) = 0; ///< Flip x and y axis, always also swaps voxel size
  virtual void FlipXZ(bool modify_origin = false) = 0; ///< Flip x and z axis, always also swaps voxel size
  virtual void FlipYZ(bool modify_origin = false) = 0; ///< Flip y and z axis, always also swaps voxel size
  virtual void FlipXT(bool modify_origin = false) = 0; ///< Flip x and t axis, always also swaps voxel size
  virtual void FlipYT(bool modify_origin = false) = 0; ///< Flip y and t axis, always also swaps voxel size
  virtual void FlipZT(bool modify_origin = false) = 0; ///< Flip z and t axis, always also swaps voxel size

  virtual void SwapXY(bool modify_axes = true) = 0; ///< Swap x and y axis
  virtual void SwapXZ(bool modify_axes = true) = 0; ///< Swap x and z axis
  virtual void SwapYZ(bool modify_axes = true) = 0; ///< Swap y and z axis
  virtual void SwapXT(bool modify_axes = true) = 0; ///< Swap x and t axis
  virtual void SwapYT(bool modify_axes = true) = 0; ///< Swap y and t axis
  virtual void SwapZT(bool modify_axes = true) = 0; ///< Swap z and t axis

  // ---------------------------------------------------------------------------
  // VTK interface
  #if MIRTK_Image_WITH_VTK

  /// Return the corresponding VTK scalar type
  ///
  /// \note Use only when MIRTK_Image_WITH_VTK is 1.
  int ImageToVTKScalarType() const;

  /// Convert image to VTK structured points
  ///
  /// \note Use only when MIRTK_Image_WITH_VTK is 1.
  virtual void ImageToVTK(vtkStructuredPoints *) const;

  /// Convert VTK structured points to image
  ///
  /// \note Use only when MIRTK_Image_WITH_VTK is 1.
  virtual void VTKToImage(vtkStructuredPoints *);

  #endif // MIRTK_Image_WITH_VTK
  // ---------------------------------------------------------------------------
  // I/O

  /// Read image from file
  virtual void Read(const char *) = 0;

  /// Write image to file
  virtual void Write(const char *) const = 0;

  /// Print image information
  virtual void Print(Indent = 0) const;

  // ---------------------------------------------------------------------------
  // Emulation of GenericImage<VoxelType>

  // The following type definitions and "generic" voxel type specific member
  // function are required by image functions/filters which are templated
  // over the image type. If BaseImage is provided as template argument,
  // these functions and filters access the image data (indirectly) through the
  // abstract interface methods with implicit conversion to Vector.
  // For improved performance and no data conversion, these image functions
  // and filters can, however, be instantiated with a specific image type
  // such as GreyImage in which case no data conversion takes place.
  //
  // \sa GenericExtrapolateImageFunction, GenericInterpolateImageFunction

  /// Default voxel type used by generic interpolate/extrapolate image functions
  /// when instantiated with BaseImage as template argument as done by the
  /// general interpolate image functions which thus can interpolate any scalar image.
  typedef double                             VoxelType;
  typedef voxel_info<VoxelType>::ScalarType  ScalarType;
  typedef voxel_info<VoxelType>::RealType    RealType;
  typedef voxel_info<RealType>::ScalarType   RealScalarType;

  /// Get pixel value at voxel with given index
  VoxelType Get(int) const;

  /// Get pixel value at voxel with given lattice coordinates
  VoxelType Get(int, int, int = 0, int = 0) const;

  // ---------------------------------------------------------------------------
  // Deprecated

  /// Gets the image attributes
  /// \deprecated Use Attributes instead.
  const ImageAttributes &GetImageAttributes() const;

  /// Returns the total number of voxels
  /// \deprecated Use NumberOfVoxels instead.
  int GetNumberOfVoxels() const;

  /// Voxel dimensions get access
  void GetPixelSize(double *, double *, double *) const;

  /// Voxel dimensions get access
  void GetPixelSize(double *, double *, double *, double *) const;

  /// Minimum and maximum pixel values get accessor
  void GetMinMaxAsDouble(double *, double *) const;

  /// \returns Raw pointer to contiguous image data.
  /// \deprecated Use GetDataPointer instead.
  void *GetScalarPointer(int = 0, int = 0, int = 0, int = 0);

  /// \returns Raw pointer to contiguous image data.
  /// \deprecated Use GetDataPointer instead.
  const void *GetScalarPointer(int = 0, int = 0, int = 0, int = 0) const;

  /// \returns Enumeration value corresponding to pixel scalar type.
  /// \deprecated Use GetVoxelType instead.
  int GetScalarType() const;

  /// \returns Size of pixel scalar type in bytes.
  /// \deprecated Use GetVoxelTypeSize instead.
  int GetScalarTypeSize() const;

  /// \returns Minimum value a pixel can hold without overflowing.
  /// \deprecated Use GetVoxelTypeMin instead.
  double GetScalarTypeMin() const;

  /// \returns Maximum value a pixel can hold without overflowing.
  /// \deprecated Use GetVoxelTypeMax instead.
  double GetScalarTypeMax() const;

};

/// Alternative/backwards compatible type name
typedef BaseImage Image;


} // namespace mirtk

#endif // MIRTK_BaseImage_H

////////////////////////////////////////////////////////////////////////////////
// Inline definitions of BaseImage
////////////////////////////////////////////////////////////////////////////////

#ifndef MIRTK_BaseImage_HH
#define MIRTK_BaseImage_HH

// BinaryImage must be completely defined here already such that it can
// be used in the inline definitions. The methods using the mask must be inline
// in order for them to be efficient enough as they will be called for each voxel
#include "mirtk/GenericImage.h"
#include "mirtk/Math.h"


namespace mirtk {

// =============================================================================
// Lattice
// =============================================================================

// -----------------------------------------------------------------------------
inline const ImageAttributes &BaseImage::Attributes() const
{
  return _attr;
}

// -----------------------------------------------------------------------------
inline int BaseImage::NumberOfVoxels() const
{
  return _NumberOfVoxels;
}

// -----------------------------------------------------------------------------
inline int BaseImage::NumberOfSpatialVoxels() const
{
  return _attr.NumberOfSpatialPoints();
}

// -----------------------------------------------------------------------------
inline int BaseImage::X() const
{
  return _attr._x;
}

// -----------------------------------------------------------------------------
inline int BaseImage::Y() const
{
  return _attr._y;
}

// -----------------------------------------------------------------------------
inline int BaseImage::Z() const
{
  return _attr._z;
}

// -----------------------------------------------------------------------------
inline int BaseImage::T() const
{
  return _attr._t;
}

// -----------------------------------------------------------------------------
inline int BaseImage::N() const
{
  return 1;
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetX() const
{
  return _attr._x;
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetY() const
{
  return _attr._y;
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetZ() const
{
  return _attr._z;
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetT() const
{
  return _attr._t;
}

// -----------------------------------------------------------------------------
inline double BaseImage::XSize() const
{
  return _attr._dx;
}

// -----------------------------------------------------------------------------
inline double BaseImage::YSize() const
{
  return _attr._dy;
}

// -----------------------------------------------------------------------------
inline double BaseImage::ZSize() const
{
  return _attr._dz;
}

// -----------------------------------------------------------------------------
inline double BaseImage::TSize() const
{
  return _attr._dt;
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetXSize() const
{
  return _attr._dx;
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetYSize() const
{
  return _attr._dy;
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetZSize() const
{
  return _attr._dz;
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetTSize() const
{
  return _attr._dt;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutTSize(double dt)
{
  _attr._dt = dt;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutPixelSize(double dx, double dy, double dz)
{
  _attr._dx = dx;
  _attr._dy = dy;
  _attr._dz = dz;
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutPixelSize(double dx, double dy, double dz, double dt)
{
  _attr._dx = dx;
  _attr._dy = dy;
  _attr._dz = dz;
  _attr._dt = dt;
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetPixelSize(double &dx, double &dy) const
{
  dx = _attr._dx;
  dy = _attr._dy;
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetPixelSize(double &dx, double &dy, double &dz) const
{
  dx = _attr._dx;
  dy = _attr._dy;
  dz = _attr._dz;
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetPixelSize(double &dx, double &dy, double &dz, double &dt) const
{
  dx = _attr._dx;
  dy = _attr._dy;
  dz = _attr._dz;
  dt = _attr._dt;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutOrigin(const Point &p)
{
  _attr._xorigin = p._x;
  _attr._yorigin = p._y;
  _attr._zorigin = p._z;
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutOrigin(double x, double y, double z)
{
  _attr._xorigin = x;
  _attr._yorigin = y;
  _attr._zorigin = z;
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutOrigin(double x, double y, double z, double t)
{
  _attr._xorigin = x;
  _attr._yorigin = y;
  _attr._zorigin = z;
  _attr._torigin = t;
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetOrigin(double &x, double &y, double &z) const
{
  x = _attr._xorigin;
  y = _attr._yorigin;
  z = _attr._zorigin;
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetOrigin(double &x, double &y, double &z, double &t) const
{
  x = _attr._xorigin;
  y = _attr._yorigin;
  z = _attr._zorigin;
  t = _attr._torigin;
}

// -----------------------------------------------------------------------------
inline Point BaseImage::GetOrigin() const
{
  return Point(_attr._xorigin, _attr._yorigin, _attr._zorigin);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutTOrigin(double t)
{
  _attr._torigin = t;
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetTOrigin() const
{
  return _attr._torigin;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutOrientation(double *xaxis, double *yaxis, double *zaxis)
{
  const size_t n = 3 * sizeof(double);
  memcpy(_attr._xaxis, xaxis, n);
  memcpy(_attr._yaxis, yaxis, n);
  if (zaxis) memcpy(_attr._zaxis, zaxis, n);
  else {
    _attr._zaxis[0] = _attr._xaxis[1] * _attr._yaxis[2] - _attr._xaxis[2] * _attr._yaxis[1];
    _attr._zaxis[1] = _attr._xaxis[2] * _attr._yaxis[0] - _attr._xaxis[0] * _attr._yaxis[2];
    _attr._zaxis[2] = _attr._xaxis[0] * _attr._yaxis[1] - _attr._xaxis[1] * _attr._yaxis[0];
  }
  UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetOrientation(double *xaxis, double *yaxis, double *zaxis) const
{
  const size_t n = 3 * sizeof(double);
  memcpy(xaxis, _attr._xaxis, n);
  memcpy(yaxis, _attr._yaxis, n);
  if (zaxis) memcpy(zaxis, _attr._zaxis, n);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAffineMatrix(const Matrix &mat, bool apply)
{
  _attr.PutAffineMatrix(mat, apply);
  this->UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline const Matrix &BaseImage::GetAffineMatrix() const
{
  return _attr._smat;
}

// -----------------------------------------------------------------------------
inline void BaseImage::ResetAffineMatrix()
{
  _attr._smat.Ident();
  this->UpdateMatrix();
}

// -----------------------------------------------------------------------------
inline int BaseImage::VoxelToIndex(int i, int j, int k, int l) const
{
  return _attr.LatticeToIndex(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToVoxel(int idx, int &i, int &j) const
{
  _attr.IndexToLattice(idx, i, j);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToVoxel(int idx, int &i, int &j, int &k) const
{
  _attr.IndexToLattice(idx, i, j, k);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToVoxel(int idx, int &i, int &j, int &k, int &l) const
{
  _attr.IndexToLattice(idx, i, j, k, l);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToWorld(int idx, double &x, double &y) const
{
  _attr.IndexToWorld(idx, x, y);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToWorld(int idx, double &x, double &y, double &z) const
{
  _attr.IndexToWorld(idx, x, y, z);
}

// -----------------------------------------------------------------------------
inline void BaseImage::IndexToWorld(int idx, Point &p) const
{
  _attr.IndexToWorld(idx, p);
}

// -----------------------------------------------------------------------------
inline Point BaseImage::IndexToWorld(int idx) const
{
  return _attr.IndexToWorld(idx);
}

// -----------------------------------------------------------------------------
inline void BaseImage::ImageToWorld(double &x, double &y) const
{
  double a = _matI2W(0, 0) * x + _matI2W(0, 1) * y + _matI2W(0, 3);
  double b = _matI2W(1, 0) * x + _matI2W(1, 1) * y + _matI2W(1, 3);
  x = a, y = b;
}

// -----------------------------------------------------------------------------
inline void BaseImage::ImageToWorld(double &x, double &y, double &z) const
{
  double a = _matI2W(0, 0) * x + _matI2W(0, 1) * y + _matI2W(0, 2) * z + _matI2W(0, 3);
  double b = _matI2W(1, 0) * x + _matI2W(1, 1) * y + _matI2W(1, 2) * z + _matI2W(1, 3);
  double c = _matI2W(2, 0) * x + _matI2W(2, 1) * y + _matI2W(2, 2) * z + _matI2W(2, 3);
  x = a, y = b, z = c;
}

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

// -----------------------------------------------------------------------------
inline void BaseImage::ImageToWorld(Vector3 &v) const
{
  double a = _matI2W(0, 0) * v._x + _matI2W(0, 1) * v._y + _matI2W(0, 2) * v._z;
  double b = _matI2W(1, 0) * v._x + _matI2W(1, 1) * v._y + _matI2W(1, 2) * v._z;
  double c = _matI2W(2, 0) * v._x + _matI2W(2, 1) * v._y + _matI2W(2, 2) * v._z;
  v._x = a, v._y = b, v._z = c;
}

// -----------------------------------------------------------------------------
inline void BaseImage::WorldToImage(double &x, double &y) const
{
  double a = _matW2I(0, 0) * x + _matW2I(0, 1) * y + _matW2I(0, 3);
  double b = _matW2I(1, 0) * x + _matW2I(1, 1) * y + _matW2I(1, 3);
  x = a, y = b;
}

// -----------------------------------------------------------------------------
inline void BaseImage::WorldToImage(double &x, double &y, double &z) const
{
  double a = _matW2I(0, 0) * x + _matW2I(0, 1) * y + _matW2I(0, 2) * z + _matW2I(0, 3);
  double b = _matW2I(1, 0) * x + _matW2I(1, 1) * y + _matW2I(1, 2) * z + _matW2I(1, 3);
  double c = _matW2I(2, 0) * x + _matW2I(2, 1) * y + _matW2I(2, 2) * z + _matW2I(2, 3);
  x = a, y = b, z = c;
}

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

// -----------------------------------------------------------------------------
inline void BaseImage::WorldToImage(Vector3 &v) const
{
  double a = _matW2I(0, 0) * v._x + _matW2I(0, 1) * v._y + _matW2I(0, 2) * v._z;
  double b = _matW2I(1, 0) * v._x + _matW2I(1, 1) * v._y + _matW2I(1, 2) * v._z;
  double c = _matW2I(2, 0) * v._x + _matW2I(2, 1) * v._y + _matW2I(2, 2) * v._z;
  v._x = a, v._y = b, v._z = c;
}

// -----------------------------------------------------------------------------
inline const Matrix &BaseImage::GetImageToWorldMatrix() const
{
  return _matI2W;
}

// -----------------------------------------------------------------------------
inline const Matrix &BaseImage::GetWorldToImageMatrix() const
{
  return _matW2I;
}

// -----------------------------------------------------------------------------
inline double BaseImage::ImageToTime(double t) const
{
  return _attr._torigin + t * _attr._dt;
}

// -----------------------------------------------------------------------------
inline double BaseImage::TimeToImage(double t) const
{
  return (_attr._dt ? ((t - _attr._torigin) / _attr._dt) : .0);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::HasSpatialAttributesOf(const BaseImage *other) const
{
  return _attr.EqualInSpace(other->GetImageAttributes());
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInFOV(double x, double y, double z)
{
  this->WorldToImage(x, y, z);
  return (-0.5 <= x && x < _attr._x - 0.5) &&
         (-0.5 <= y && y < _attr._y - 0.5) &&
         (-0.5 <= z && z < _attr._z - 0.5);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsEmpty() const
{
  return ((_attr._x < 1) || (_attr._y < 1) || (_attr._z < 1) || (_attr._t < 1));
}

// =============================================================================
// Type independent access to image data
// =============================================================================

// -----------------------------------------------------------------------------
inline double BaseImage::GetAsDouble(int idx) const
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  return this->GetAsDouble(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsDouble(int idx, double v)
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  this->PutAsDouble(i, j, k, l, v);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsDouble(int i, int j, double v)
{
  this->PutAsDouble(i, j, 0, 0, v);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsDouble(int i, int j, int k, double v)
{
  this->PutAsDouble(i, j, k, 0, v);
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetAsVector(Vector &v, int idx) const
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  this->GetAsVector(v, i, j, k, l);
}

// -----------------------------------------------------------------------------
inline Vector BaseImage::GetAsVector(int idx) const
{
  Vector v;
  this->GetAsVector(v, idx);
  return v;
}

// -----------------------------------------------------------------------------
inline Vector BaseImage::GetAsVector(int i, int j, int k, int l) const
{
  Vector v;
  this->GetAsVector(v, i, j, k, l);
  return v;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsVector(int idx, const Vector &v)
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  this->PutAsVector(i, j, k, l, v);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsVector(int i, int j, const Vector &v)
{
  this->PutAsVector(i, j, 0, 0, v);
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutAsVector(int i, int j, int k, const Vector &v)
{
  this->PutAsVector(i, j, k, 0, v);
}

// =============================================================================
// Foreground region
// =============================================================================

// -----------------------------------------------------------------------------
inline BinaryImage *BaseImage::GetMask(bool surrender_ownership)
{
  if (surrender_ownership) _maskOwner = false;
  return _mask;
}

// -----------------------------------------------------------------------------
inline const BinaryImage *BaseImage::GetMask() const
{
  return _mask;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::HasMask() const
{
  return _mask != NULL;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::OwnsMask() const
{
  return _mask != NULL && _maskOwner;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutBackgroundValueAsDouble(double value)
{
  _bg    = value;
  _bgSet = true;
}

// -----------------------------------------------------------------------------
inline void BaseImage::PutBackgroundValueAsDouble(double value, bool threshold)
{
  _bg    = value;
  _bgSet = true;
  if (threshold && !IsNaN(_bg)) {
    for (int idx = 0; idx < _NumberOfVoxels; idx++) {
      if (this->GetAsDouble(idx) < _bg) this->PutAsDouble(idx, _bg);
    }
  }
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetBackgroundValueAsDouble() const
{
  return _bg;
}

// -----------------------------------------------------------------------------
inline void BaseImage::ClearBackgroundValue()
{
  _bgSet = false;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::HasBackgroundValue() const
{
  return _bgSet;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBackground(int idx) const
{
  if (_mask) {
    if (_mask->T() != _attr._t) idx = idx % (_attr._x * _attr._y * _attr._z);
    return (_mask->Get(idx) == BinaryPixel(0));
  } else if (_bgSet) {
    return AreEqualOrNaN(this->GetAsDouble(idx), _bg, 1e-6);
  }
  return false;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBackground(int i, int j, int k, int l) const
{
  if (_mask) {
    if (_mask->T() != _attr._t) l = 0;
    return (_mask->Get(i, j, k, l) == BinaryPixel(0));
  } else if (_bgSet) {
    return AreEqualOrNaN(this->GetAsDouble(i, j, k, l), _bg, 1e-6);
  }
  return false;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsForeground(int idx) const
{
  return !IsBackground(idx);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsForeground(int i, int j, int k, int l) const
{
  return !IsBackground(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInside(int idx) const
{
  return (0 <= idx && idx < _NumberOfVoxels);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInside(int i, int j) const
{
  return _attr.IsInside(i, j);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInside(int i, int j, int k) const
{
  return _attr.IsInside(i, j, k);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInside(int i, int j, int k, int l) const
{
  return _attr.IsInside(i, j, k, l);
}

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

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

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

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

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundary(int idx) const
{
  return _attr.IsBoundary(idx);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundary(int i, int j) const
{
  return _attr.IsBoundary(i, j);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundary(int i, int j, int k) const
{
  return _attr.IsBoundary(i, j, k);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundary(int i, int j, int k, int l) const
{
  return _attr.IsBoundary(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInsideForeground(int idx) const
{
  return IsInside(idx) && IsForeground(idx);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsInsideForeground(int i, int j, int k, int l) const
{
  return IsInside(i, j, k, l) && IsForeground(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsOutsideForeground(int idx) const
{
  return !IsInsideForeground(idx);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsOutsideForeground(int i, int j, int k, int l) const
{
  return !IsInsideForeground(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundingBoxInsideForeground(int i1, int j1,
                                                     int i2, int j2) const
{
  if (i1 < 0 || i2 >= _attr._x ||
      j1 < 0 || j2 >= _attr._y) return false;
  if (_mask) {
    for (int j = j1; j != j2; j++) {
      for (int i = i1; i != i2; i++) {
        if (!_mask->Get(i, j)) return false;
      }
    }
  } else if (_bgSet) {
    for (int j = j1; j != j2; j++) {
      for (int i = i1; i != i2; i++) {
        if (AreEqualOrNaN(this->GetAsDouble(i, j), _bg, 1e-6)) return false;
      }
    }
  }
  return true;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundingBoxInsideForeground(int i1, int j1, int k1,
                                                     int i2, int j2, int k2) const
{
  if (i1 < 0 || i2 >= _attr._x ||
      j1 < 0 || j2 >= _attr._y ||
      k1 < 0 || k2 >= _attr._z) return false;
  if (_mask) {
    for (int k = k1; k != k2; k++) {
      for (int j = j1; j != j2; j++) {
        for (int i = i1; i != i2; i++) {
          if (!_mask->Get(i, j, k)) return false;
        }
      }
    }
  } else if (_bgSet) {
    for (int k = k1; k != k2; k++) {
      for (int j = j1; j != j2; j++) {
        for (int i = i1; i != i2; i++) {
          if (AreEqualOrNaN(this->GetAsDouble(i, j, k), _bg, 1e-6)) return false;
        }
      }
    }
  }
  return true;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsBoundingBoxInsideForeground(int i1, int j1, int k1, int l1,
                                                     int i2, int j2, int k2, int l2) const
{
  if (i1 < 0 || i2 >= _attr._x ||
      j1 < 0 || j2 >= _attr._y ||
      k1 < 0 || k2 >= _attr._z ||
      l1 < 0 || l2 >= _attr._t) return false;
  if (_mask) {
    if (_mask->T() == 1) {
      for (int k = k1; k != k2; k++) {
        for (int j = j1; j != j2; j++) {
          for (int i = i1; i != i2; i++) {
            if (!_mask->Get(i, j, k)) return false;
          }
        }
      }
    } else {
      for (int l = l1; l != l2; l++) {
        for (int k = k1; k != k2; k++) {
          for (int j = j1; j != j2; j++) {
            for (int i = i1; i != i2; i++) {
              if (!_mask->Get(i, j, k, l)) return false;
            }
          }
        }
      }
    }
  } else if (_bgSet) {
    double value;
    for (int l = l1; l != l2; l++) {
      for (int k = k1; k != k2; k++) {
        for (int j = j1; j != j2; j++) {
          for (int i = i1; i != i2; i++) {
            value = this->GetAsDouble(i, j, k, l);
            if (AreEqualOrNaN(this->GetAsDouble(i, j, k, l), _bg, 1e-6)) return false;
          }
        }
      }
    }
  }
  return true;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsNextToBackground(int i, int j, int k, int l) const
{
  for (int nl = l - 1; nl <= l + 1; ++nl)
  for (int nk = k - 1; nk <= k + 1; ++nk)
  for (int nj = j - 1; nj <= j + 1; ++nj)
  for (int ni = i - 1; ni <= i + 1; ++ni) {
    if (ni != 0 || nj != 0 || nk != 0 || nl != 0) {
      if (IsOutsideForeground(ni, nj, nk, nl)) {
        return true;
      }
    }
  }
  return false;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsNextToBackground(int idx) const
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  return IsNextToBackground(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsNextToForeground(int i, int j, int k, int l) const
{
  for (int nl = l - 1; nl <= l + 1; ++nl)
  for (int nk = k - 1; nk <= k + 1; ++nk)
  for (int nj = j - 1; nj <= j + 1; ++nj)
  for (int ni = i - 1; ni <= i + 1; ++ni) {
    if (ni != 0 || nj != 0 || nk != 0 || nl != 0) {
      if (IsInsideForeground(ni, nj, nk, nl)) {
        return true;
      }
    }
  }
  return false;
}

// -----------------------------------------------------------------------------
inline bool BaseImage::IsNextToForeground(int idx) const
{
  int i, j, k, l;
  IndexToVoxel(idx, i, j, k, l);
  return IsNextToForeground(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline bool BaseImage::HasBackground() const
{
  if (HasBackgroundValue() || GetMask()) {
    for (int i = 0; i < _NumberOfVoxels; i++) {
      if (IsBackground(i)) return true;
    }
  }
  return false;
}

// =============================================================================
// Emulation of GenericImage<VoxelType>
// =============================================================================

// -----------------------------------------------------------------------------
inline double BaseImage::Get(int idx) const
{
  return this->GetAsDouble(idx);
}

// -----------------------------------------------------------------------------
inline double BaseImage::Get(int i, int j, int k, int l) const
{
  return this->GetAsDouble(i, j, k, l);
}

// =============================================================================
// Deprecated
// =============================================================================

// -----------------------------------------------------------------------------
inline const ImageAttributes &BaseImage::GetImageAttributes() const
{
  return Attributes();
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetNumberOfVoxels() const
{
  return NumberOfVoxels();
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetPixelSize(double *dx, double *dy, double *dz) const
{
  *dx = _attr._dx;
  *dy = _attr._dy;
  *dz = _attr._dz;
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetPixelSize(double *dx, double *dy, double *dz, double *dt) const
{
  *dx = _attr._dx;
  *dy = _attr._dy;
  *dz = _attr._dz;
  *dt = _attr._dt;
}

// -----------------------------------------------------------------------------
inline void BaseImage::GetMinMaxAsDouble(double *min, double *max) const
{
  this->GetMinMaxAsDouble(*min, *max);
}

// -----------------------------------------------------------------------------
inline void *BaseImage::GetScalarPointer(int i, int j, int k, int l)
{
  return this->GetDataPointer(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline const void *BaseImage::GetScalarPointer(int i, int j, int k, int l) const
{
  return this->GetDataPointer(i, j, k, l);
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetScalarType() const
{
  return this->GetDataType();
}

// -----------------------------------------------------------------------------
inline int BaseImage::GetScalarTypeSize() const
{
  return this->GetDataTypeSize();
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetScalarTypeMin() const
{
  return this->GetDataTypeMin();
}

// -----------------------------------------------------------------------------
inline double BaseImage::GetScalarTypeMax() const
{
  return this->GetDataTypeMax();
}


} // namespace mirtk

#endif // MIRTK_BaseImage_HH