RealVoxel.H

Go to the documentation of this file.

/*!@file PointCloud/RealVoxel.H Open Scene Graph Utils */

// //////////////////////////////////////////////////////////////////// //
// The iLab Neuromorphic Vision C++ Toolkit - Copyright (C) 2000-2003   //
// by the University of Southern California (USC) and the iLab at USC.  //
// See http://iLab.usc.edu for information about this project.          //
// //////////////////////////////////////////////////////////////////// //
// Major portions of the iLab Neuromorphic Vision Toolkit are protected //
// under the U.S. patent ``Computation of Intrinsic Perceptual Saliency //
// in Visual Environments, and Applications'' by Christof Koch and      //
// Laurent Itti, California Institute of Technology, 2001 (patent       //
// pending; application number 09/912,225 filed July 23, 2001; see      //
// http://pair.uspto.gov/cgi-bin/final/home.pl for current status).     //
// //////////////////////////////////////////////////////////////////// //
// This file is part of the iLab Neuromorphic Vision C++ Toolkit.       //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit 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.     //
//                                                                      //
// The iLab Neuromorphic Vision C++ Toolkit is distributed in the hope  //
// that it will be useful, but WITHOUT ANY WARRANTY; without even the   //
// implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR      //
// PURPOSE.  See the GNU General Public License for more details.       //
//                                                                      //
// You should have received a copy of the GNU General Public License    //
// along with the iLab Neuromorphic Vision C++ Toolkit; if not, write   //
// to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,   //
// Boston, MA 02111-1307 USA.                                           //
// //////////////////////////////////////////////////////////////////// //
//
// Primary maintainer for this file: T. Nathan Mundhenk <mundhenk@usc.edu>
// $HeadURL: svn://isvn.usc.edu/software/invt/trunk/saliency/src/PointCloud/RealVoxel.H $
// $Id: RealVoxel.H 10794 2009-02-08 06:21:09Z itti $
//

#ifndef REALVOXEL_H_DEFINED
#define REALVOXEL_H_DEFINED

#include "Image/ArrayData.H"
#include "Image/Pixels.H"
#include "PointCloud/Point.H"
#include "PointCloud/VoxelImage.H"

#include <cstdlib>
#include <iostream>
#include <string>

// ######################################################################
// ######################################################################
// ######################################################################
//! Create, store and get real voxels (non-sparse voxel images), Base Class
template <class T>
class RealVoxel : public
VoxelImage<T,RealVoxel<T>,ArrayHandle<T> >
{
public:
  //typedef T VOX;
  typedef ArrayHandle<T> BASE_IMAGE;

  //! Create a new real voxel image which is blank
  RealVoxel();

  //! Copy a new real voxel image of size x,y,z
  RealVoxel(const T* inarray,
            const unsigned int size_x,
            const unsigned int size_y,
            const unsigned int size_z);

  //! Create a new real voxel image of size x,y,z
  RealVoxel(const unsigned int size_x,
            const unsigned int size_y,
            const unsigned int size_z,
            const InitPolicy init = NO_INIT);

  //! Create a new real voxel image of coord size
  RealVoxel(const Point<int,3>& coord_size,
            const InitPolicy init = NO_INIT);

  //! Create a new real voxel image from another
  RealVoxel(const RealVoxel<T>& A);

  virtual ~RealVoxel();

  //! Free memory and switch to uninitialized state.
  /*! Note that it is \b NOT necessary to call this function to ensure
      proper cleanup, that will be done in the destructor by
      default. Rather, freeMem() is offered just as a performance
      optimization, to allow you to release a potentially large chunk
      of memory when you are finished using it. */
  void freeMem();

  //! Swap the contents of two images
  void swap(RealVoxel<T>& other);

  //! Use existing memory.
  /*! This is potentially dangerous and should really be avoided. The
      only case where this really is useful is to attach an image to
      an existing memory segment that is shared or into which data is
      streaming via DMA. */
  void attach(T* array,
                     const unsigned int width,
                     const unsigned int height,
                     const unsigned int depth);

  //! Detach previously attach()'ed image.
  /*! The main purpose of detach() is to make sure that this Image
      object does not continue to point at attach()'ed memory after
      that memory has been freed. Nevertheless, this function is not
      strictly necessary to ensure correct memory handling, since the
      Image destructor will only try to free memory if that memory is
      owned (i.e., not attach()'ed). All that it does is release any
      association with a previously attach()'ed memory block, by
      setting the current image to a new empty (zero-by-zero)
      ArrayData object. We assume that attached memory will be
      destroyed later. */
  void detach();

  //! Return a new image object with a deep copy of the underlying data
  RealVoxel<T> deepcopy() const;

  //! Free mem and realloc new array (array contents are lost).
  void resize(const unsigned int size_x,
              const unsigned int size_y,
              const unsigned int size_z,
              const bool do_clear = false);

  //! Free mem and realloc new array (array contents are lost).
  void resize(const Point<int,3>& coord_size,
                     const bool do_clear = false);

  //! clear contents (or set to given value)
  void clear(const T& val = T());

  //! inefficent access to get x,y,z if possible use iterator
  /*! if the virtual voxel does not exist (nothing is there) then it will
      return false and vox will be null.
  */
  virtual bool getVal(const unsigned int x,
                      const unsigned int y,
                      const unsigned int z,
                      T& vox);

  //! inefficent access to get x,y,z if possible use iterator
  bool getVal(const Point<int,3>& coord,T& vox);

  //! Return the virtual voxels. Bool will be false if vox is empty.
  virtual bool get(BASE_IMAGE& ivox) const;

  //! store a voxel at x,y,z
  virtual void setVal(const unsigned int x,
                      const unsigned int y,
                      const unsigned int z,
                      const T& vox);

  //! store a voxel at x,y,z using a point
  void setVal(const Point<int,3>& coord,const T& vox);

  //! Sets the virtual voxel
  virtual void set(const BASE_IMAGE& ivox);

  //! Returns read-only (const) pointer to internal image array
  const T* getArrayPtr() const;

  //! Returns read/write (non-const) pointer to internal image array
  T* getArrayPtr();

  //! Test whether point falls inside array boundaries
  bool coordsOk(const unsigned int i,
                const unsigned int j,
                const unsigned int k) const;

  //! Test whether point falls inside array boundaries
  bool coordsOk(const Point<int,3>& coord) const;

  //! Check whether image is non-empty (i.e., non-zero height and width).
  virtual bool initialized() const;

#ifndef INVT_MEM_DEBUG
  // standard iterators are just pointers to the data:
  //! Read/write iterator.
  typedef T* iterator;
  //! Read-only iterator.
  typedef const T* const_iterator;
#else
  // range-checked iterators can detect attempts to access past image bounds:
  //! Read/write iterator.
  typedef CheckedIterator<T> iterator;
  //! Read-only iterator.
  typedef CheckedIterator<const T> const_iterator;
#endif

  //! Returns the constant head iterator
  const_iterator begin() const;

  //! Returns the constant end iterator
  const_iterator end() const;

  //! Returns the head iterator
  iterator beginw();

  //! Returns the end iterator
  iterator endw();

  //! get the range of the voxels in real space
  //virtual void range(VoxRange& range) const;

  //! number of points we now have (same as getSize)
  virtual int size() const;

  //! Return an Image Slice, set the normal axis' to -1
  /*! This will return a single INVT image slice of the full voxel
      image. To use this, set the values in the Point object to -1
      for the axis that lay on the image plane.

      Example: I want to return an image at a depth of 12. That is, z=12

      In this case, the returned image will be of size X (width) by Y (height)
      make a Point like Point<int,3>(-1,-1,12) and pass it with the
      method call.

      The bool is for compatability with sparse voxel images
  */
  virtual bool getSlice(const Point<int,3>& coord, Image<T>& img) const;

  //! Return a vector column, set the normal axis to -1
  /*! This will return a column in a std::vector from the voxel. To use
      it, set the ortho axis to -1. So if you want to get a column at
      X = 12 and Y = 18 then set Z to -1 since its the orthogonal plane.

      Example: I want to get a column of voxels from offset X = 12 and Y = 18

      The returned vector will be the of size Z (depth). Creat a point like
      Point<int,3>(12,18,-1). Since the Z is -1, the method will know
      we want the column othogonal to the z plane.
  */
  virtual bool getColumn(const Point<int,3>& coord, std::vector<T>& vec) const;

  //! Return a single value (voxel)
  //virtual bool getVal(const Point<int,3>& coord, T& val);

  //! Set an Image Slice, set the normal axis' to -1
  /*! This will set a single voxel slice to the value of an INVT image
      To use this, set the values in the Point object to -1
      for the axis that lay on the image plane. The image must be the same
      size dimensions parallel to the image slice. So if setting a
      slice at Z=12, then the X (width) and Y (height) must be the same size
      as the voxel image.

      Example: I want to set a slice at a row of 22. That is, y=22

      In this case, the image must be of size X (width) by Z (depth)
      make a Point like Point<int,3>(-1,22,-1) and pass it with the
      method call.
  */
  virtual void setSlice(const Point<int,3>& coord, const Image<T>& img);

  //! Set a vector column, set the normal axis to -1
  /*! This will set a column with a std::vector of voxels. To use
      it, set the ortho axis to -1. So if you want to set a column at
      X = 12 and Y = 18 then set Z to -1 since its the orthogonal plane. The
      vector must be the same size as the axis requested. So if the image
      size Z (depth) is 256, then the vector needs to be that size in order
      to set along the z axis.

      Example: I want to set a column of voxels from offset X = 12 and Y = 18

      The the vector must be the of size Z (depth). Creat a point like
      Point<int,3>(12,18,-1). Since the Z is -1, the method will know
      we want the column othogonal to the z plane.
  */
  virtual void setColumn(const Point<int,3>& coord, const std::vector<T>& vec);

  //! Set a single value
  //virtual void setVal(const Point<int,3>& coord, const T& val);

  //! get our height - Y (itsSizeY)
  virtual unsigned int getHeight() const;

  //! get our width - X (itsSizeX)
  virtual unsigned int getWidth() const;

  //! get out depth - Z (itsSizeZ)
  virtual unsigned int getDepth() const;

  //! get total size (number of voxels)
  unsigned int getSize() const;

  //! Check if *this is the same size as the other thing
  template <class C>
  bool isSameSize(const C& other) const;

  //! Check if the image is 1D, i.e., width == 1 or height == 1
  bool is1D() const;

  //! Check if the image is square, i.e., width == height == depth
  /*! Ok, it's really a cube, but we keep the name from the Image class */
  bool isSquare() const;

  //! Access voxel image elements through C array index interface
  /*! To get an element at x,y,z use the following computation:

    x + y * width + z * width * height

    Example: get a voxel at position (x,y,z) = (13,23,56) in RealVoxel
    image ivox...

    T voxel = ivox[13 + 23 * ivox.getWidth() +
    56 * ivox.getWidth() * ivox.getHeight()];
  */
  T& operator[](const unsigned int index);

  //! Access voxel image elements through C array index interface
  /*! To get an element at x,y,z use the following computation:

    x + y * width + z * width * height

    Example: get a voxel at position (x,y,z) = (13,23,56) in RealVoxel
    image ivox...

    T voxel = ivox[13 + 23 * ivox.getWidth() +
    56 * ivox.getWidth() * ivox.getHeight()];
  */
  const T& operator[](const unsigned int index) const;

  //! Access voxel image elements through C array index interface
  T& operator[](const Point<int,3>& p);

  //! Access voxel image elements through C array index interface
  const T& operator[](const Point<int,3>& p) const;

  //! Equality Operator
  bool operator==(const RealVoxel<T>& ivox) const;

  //! copy operator
  virtual RealVoxel<T> & operator=(const RealVoxel<T>& ivox);

  //! RealVoxel + RealVoxel operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator+(const RealVoxel<T2>& ivox) const;

  //! RealVoxel - RealVoxel operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator-(const RealVoxel<T2>& ivox) const;

  //! RealVoxel * RealVoxel operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator*(const RealVoxel<T2>& ivox) const;

  //! RealVoxel / RealVoxel operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator/(const RealVoxel<T2>& ivox) const;

  //! RealVoxel operator which does the work
  template <class T2, class MATH_OP>
  RealVoxel<typename promote_trait<T,T2>::TP>
  basicOp(const RealVoxel<T2>& ivox) const;

  //! RealVoxel += RealVoxel operator
  template <class T2>
  RealVoxel<T> & operator+=(const RealVoxel<T2>& ivox);

  //! RealVoxel -= RealVoxel operator
  template <class T2>
  RealVoxel<T> & operator-=(const RealVoxel<T2>& ivox);

  //! RealVoxel *= RealVoxel operator
  template <class T2>
  RealVoxel<T> & operator*=(const RealVoxel<T2>& ivox);

  //! RealVoxel /= RealVoxel operator
  template <class T2>
  RealVoxel<T> & operator/=(const RealVoxel<T2>& ivox);

  //! copy operator
  virtual RealVoxel<T> & operator=(const T& val);

  //! RealVoxel + value operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator+(const T2& val) const;

  //! RealVoxel - value operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator-(const T2& val) const;

  //! RealVoxel * value operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator*(const T2& val) const;

  //! RealVoxel / value operator
  template <class T2>
  RealVoxel<typename promote_trait<T,T2>::TP>
  operator/(const T2& val) const;

  //! RealVoxel value operator which does the work
  template <class T2, class MATH_OP>
  RealVoxel<typename promote_trait<T,T2>::TP>
  valueOp(const T2& val) const;

  //! RealVoxel += value operator
  template <class T2>
  RealVoxel<T> & operator+=(const T2& val);

  //! RealVoxel -= value operator
  template <class T2>
  RealVoxel<T> & operator-=(const T2& val);

  //! RealVoxel *= value operator
  template <class T2>
  RealVoxel<T> & operator*=(const T2& val);

  //! RealVoxel /= value operator
  template <class T2>
  RealVoxel<T> & operator/=(const T2& val);

  // ############################################################

  //! Add an image to this voxel image at each slice along an axis
  template <class T2>
  void addImage(const imgvox::axis axis, const Image<T2>& img);

  //! Subract an image from this voxel image at each slice along an axis
  /*! Here we have 2DImage - VoxelImage */
  template <class T2>
  void subtractImage(const imgvox::axis axis, const Image<T2>& img);

  //! Subract an image from this voxel image at each slice along an axis
  /*! Here we have VoxelImage - 2DImage */
  template <class T2>
  void subtractImageFrom(const imgvox::axis axis, const Image<T2>& img);

  //! Multiply an image by this voxel image at each slice along an axis
  template <class T2>
  void multiplyImage(const imgvox::axis axis, const Image<T2>& img);

  //! Divide an image by this voxel image at each slice along an axis
  /*! Here we have 2DImage / VoxelImage */
  template <class T2>
  void divideImage(const imgvox::axis axis, const Image<T2>& img);

  //! Divide an image by this voxel image at each slice along an axis
  /*! Here we have VoxelImage / 2DImage */
  template <class T2>
  void divideImageFrom(const imgvox::axis axis, const Image<T2>& img);

  // This method does the work for all the methods listed above
  template <class T2, class MATH_OP>
  void sliceOp(const imgvox::axis axis, const Image<T2>& img);

  // ############################################################
  /*! @name Functions for testing/debugging only */
  //@{

  //! For testing/debugging only.
  /*! See if we are pointing to the same ArrayData<T> as is the other
      Image<T>. */
  bool hasSameData(const RealVoxel<T>& b) const;

  //! For testing/debugging only.
  /*! Returns the current reference count. */
  long refCount() const throw();

  //! For testing/debugging only.
  /*! Check if the ArrayHandle is shared. */
  bool isShared() const throw();

  //@}
protected:

  //! Returns a newish semi-copy (shallow copy) object for = overload
  virtual RealVoxel<T> getNew(const RealVoxel<T>& vvox);

  unsigned int itsSizeX;
  unsigned int itsSizeY;
  unsigned int itsSizeZ;
  unsigned int itsSize;

  //! Handel to the actual storage container for all the voxels
  BASE_IMAGE itsVox;

  const ArrayData<T>& impl() const;
  ArrayData<T>& uniq();


};
// ######################################################################
// ######################################################################
// ##### Special Helper Classes
// ######################################################################
// ######################################################################

template <class T> class RealVoxelAdd
{
public:
  RealVoxelAdd();
  ~RealVoxelAdd();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
template <class T> class RealVoxelSub
{
public:
  RealVoxelSub();
  ~RealVoxelSub();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
template <class T> class RealVoxelSubInv
{
public:
  RealVoxelSubInv();
  ~RealVoxelSubInv();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
template <class T> class RealVoxelMult
{
public:
  RealVoxelMult();
  ~RealVoxelMult();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
template <class T> class RealVoxelDiv
{
public:
  RealVoxelDiv();
  ~RealVoxelDiv();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
template <class T> class RealVoxelDivInv
{
public:
  RealVoxelDivInv();
  ~RealVoxelDivInv();
  template <class T2>
  typename promote_trait<T,T2>::TP op(const T& a, const T2& b) const;
};

// ######################################################################
// ######################################################################
// INLINE Methods
// ######################################################################
// ######################################################################
template <class T> inline void
RealVoxel<T>::freeMem()
{
  RealVoxel<T> empty;
  this->swap(empty);
  itsSizeX = 0; itsSizeY = 0; itsSizeZ = 0; itsSize = 0;
}

// ######################################################################
template <class T> inline void
RealVoxel<T>::swap(RealVoxel<T>& other)
{
  itsVox.swap(other.itsVox);

  // We swap the array memory , but we need to copy the
  // size values since we maintain them locally

  const unsigned int x = other.itsSizeX;
  other.itsSizeX = itsSizeX; itsSizeX = x;

  const unsigned int y = other.itsSizeY;
  other.itsSizeY = itsSizeY; itsSizeY = y;

  const unsigned int z = other.itsSizeZ;
  other.itsSizeZ = itsSizeZ; itsSizeZ = z;

  const unsigned int s = other.itsSize;
  other.itsSize  = itsSize; itsSize  = s;
}

// ######################################################################
template <class T> inline void
RealVoxel<T>::attach(T* array,
                     const unsigned int width,
                     const unsigned int height,
                     const unsigned int depth)
{
  itsSizeX = width; itsSizeY = height; itsSizeZ = depth;
  itsSize = width * height * depth;
  ArrayHandle<T> attached(new ArrayData<T>(Dims(1, itsSize),
                                           array, WRITE_THRU));
  itsVox.swap(attached);
}

// ######################################################################
template <class T> inline void
RealVoxel<T>::detach()
{
  ArrayHandle<T> emptyHdl;
  itsVox.swap(emptyHdl);
  itsSizeX = 0; itsSizeY = 0; itsSizeZ = 0; itsSize = 0;
}

// ######################################################################
template <class T> inline RealVoxel<T>
RealVoxel<T>::deepcopy() const
{
  return RealVoxel<T>(this->getArrayPtr(), itsSizeX, itsSizeY, itsSizeZ);
}

// ######################################################################
template <class T> inline void
RealVoxel<T>::resize(const unsigned int size_x,
                     const unsigned int size_y,
                     const unsigned int size_z,
                     const bool do_clear)
{

  if((size_x != itsSizeX) && (size_y != itsSizeY) && (size_z != itsSizeZ))
  {
    itsSize = size_x * size_y * size_z;
    ArrayHandle<T> resized(new ArrayData<T>(Dims(1,itsSize),
                                            do_clear ? ZEROS:NO_INIT));
    itsVox.swap(resized);
    itsSizeX = size_x;
    itsSizeY = size_y;
    itsSizeZ = size_z;
  }
  else // we're keeping the same size, so just clear() if necessary
  {
    if (do_clear) clear( T() );
  }
}

//######################################################################
template <class T> inline void
RealVoxel<T>::resize(const Point<int,3>& coord_size,
                     const bool do_clear)
{
  resize(coord_size.P[0],coord_size.P[1],coord_size.P[2],do_clear);
}

// ######################################################################
template <class T> inline
void RealVoxel<T>::clear(const T& val)
{
  // Check if we have a shared implementation... if yes, then for
  // efficiency we should release our copy before doing the iterative
  // loop, which would otherwise unnecessarily make a unique copy of
  // the shared data for us, which we would then ceremoniously ignore
  // as we clear it to a new value.
  if (itsVox.isShared())
    *this = RealVoxel<T>(itsSizeX, itsSizeY, itsSizeZ, NO_INIT);

  for (iterator itr = beginw(), stop = endw(); itr != stop; ++itr)
    *itr = val;
}

// ######################################################################
template <class T> inline bool
RealVoxel<T>::getVal(const Point<int,3>& coord,T& vox)
{
  return getVal(coord.P[0],coord.P[1],coord.P[2],vox);
}

// ######################################################################
template <class T> inline void
RealVoxel<T>::setVal(const Point<int,3>& coord,const T& vox)
{
  setVal(coord.P[0],coord.P[1],coord.P[2],vox);
}

// ######################################################################
template <class T> inline
const T* RealVoxel<T>::getArrayPtr() const
{
  return impl().data();
}

// ######################################################################
template <class T> inline
T* RealVoxel<T>::getArrayPtr()
{
  return uniq().dataw();
}

// ######################################################################
template <class T> inline bool
RealVoxel<T>::coordsOk(const unsigned int i,
                       const unsigned int j,
                       const unsigned int k) const
{
  // Since this is an unsigned int, we do  not need to check for
  // coords less than zero
  return (i < getWidth() && j < getHeight() && k < getDepth());
}

// ######################################################################
template <class T> inline bool
RealVoxel<T>::coordsOk(const Point<int,3>& coord) const
{
  // Since this is an unsigned int, we do  not need to check for
  // coords less than zero
  return (coord.P[0] < (int)itsSizeX &&
          coord.P[1] < (int)itsSizeY &&
          coord.P[2] < (int)itsSizeZ);
}

// ######################################################################
// ##### Iterators:
// ######################################################################

#ifndef INVT_MEM_DEBUG

template<class T> inline
typename RealVoxel<T>::const_iterator RealVoxel<T>::begin() const
{ return impl().data(); }

template<class T> inline
typename RealVoxel<T>::const_iterator RealVoxel<T>::end() const
{ return impl().end(); }

template<class T> inline
typename RealVoxel<T>::iterator RealVoxel<T>::beginw()
{ return uniq().dataw(); }

template<class T> inline
typename RealVoxel<T>::iterator RealVoxel<T>::endw()
{ return uniq().endw(); }

#else

template <class T> inline
typename RealVoxel<T>::const_iterator RealVoxel<T>::begin() const
{ return const_iterator(impl().data(), impl().data(), impl().end()); }

template <class T> inline
typename RealVoxel<T>::const_iterator RealVoxel<T>::end() const
{ return const_iterator(impl().end(), impl().data(), impl().end()); }

template <class T> inline
typename RealVoxel<T>::iterator RealVoxel<T>::beginw()
{ return iterator(uniq().dataw(), uniq().dataw(), uniq().endw()); }

template <class T> inline
typename RealVoxel<T>::iterator RealVoxel<T>::endw()
{ return iterator(uniq().endw(), uniq().dataw(), uniq().endw()); }

#endif

// ######################################################################
template <class T> inline unsigned int
RealVoxel<T>::getSize() const
{ return itsSize; }

// ######################################################################
template <class T> template <class C> inline
bool RealVoxel<T>::isSameSize(const C& other) const
{
  return itsSizeX == other.itsSizeX &&
         itsSizeY == other.itsSizeY &&
         itsSizeZ == other.itsSizeZ;
}

// ######################################################################
template <class T> inline
bool RealVoxel<T>::is1D() const
{ return (itsSizeX == 1) || (itsSizeY == 1) || (itsSizeZ == 1); }

// ######################################################################
template <class T> inline
bool RealVoxel<T>::isSquare() const
{ return ((itsSizeX == itsSizeY) && (itsSizeX == itsSizeZ)); }

// ######################################################################
// ##### General Overloaded RealVoxel Operators
// ######################################################################

template <class T> inline
T& RealVoxel<T>::operator[](unsigned int index)
{
  ASSERT(index >= 0 && index < this->getSize());
  return beginw()[index];
}

// ######################################################################
template <class T> inline
const T& RealVoxel<T>::operator[](unsigned int index) const
{
  ASSERT(index >= 0 && index < this->getSize());
  return begin()[index];
}

// ######################################################################
template <class T> inline
T& RealVoxel<T>::operator[](const Point<int,3>& p)
{
  ASSERT(this->coordsOk(p));
  return this->beginw()[p.P[0] +
                        p.P[1] * this->getWidth() +
                        p.P[2] * this->getWidth() * this->getHeight()];
}

// ######################################################################
template <class T> inline
const T& RealVoxel<T>::operator[](const Point<int,3>& p) const
{
  ASSERT(this->coordsOk(p));
  return this->begin()[p.P[0] +
                       p.P[1] * this->getWidth() +
                       p.P[2] * this->getWidth() * this->getHeight()];
}

// ######################################################################
template <class T> inline bool
RealVoxel<T>::operator==(const RealVoxel<T>& ivox) const
{
  if (!this->isSameSize(ivox)) return false;

  return std::equal(this->begin(), this->end(), ivox.begin());
}

// ######################################################################
template <class T> RealVoxel<T> &
RealVoxel<T>::operator=(const RealVoxel<T>& ivox)
{
  RealVoxel<T> v_copy( ivox );
  this->swap(v_copy);
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator+(const RealVoxel<T2>& ivox) const
{ return basicOp<RealVoxelAdd>(ivox); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator-(const RealVoxel<T2>& ivox) const
{ return basicOp<RealVoxelSub>(ivox); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator*(const RealVoxel<T2>& ivox) const
{ return basicOp<RealVoxelMult>(ivox); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator/(const RealVoxel<T2>& ivox) const
{ return basicOp<RealVoxelDiv>(ivox); }

// ######################################################################
// Operate on the VoxelImage using an Image Slice
template <class T> template <class T2, class MATH_OP>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::basicOp(const RealVoxel<T2>& ivox) const
{
  MATH_OP math_op;
  ASSERT(isSameSize(ivox));
  typedef typename promote_trait<T,T2>::TP TPRO;
  RealVoxel<TPRO> result(itsSizeX,itsSizeY,itsSizeZ,NO_INIT);
  typename RealVoxel<TPRO>::iterator ritr     = result.beginw();
  typename RealVoxel<T2>::const_iterator vitr = ivox.begin();
  for(typename RealVoxel<T>::const_iterator itr = begin(), stop = end();
      itr != stop; ++itr, ++vitr, ++ritr)
    *ritr = math_op.op((*itr),(*vitr));
  return result;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator+=(const RealVoxel<T2>& ivox)
{
  ASSERT(isSameSize(ivox));
  typename RealVoxel<T2>::const_iterator vitr = ivox.begin();
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr, ++vitr)
    *itr = clamped_convert<T>((*itr) + (*vitr));
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator-=(const RealVoxel<T2>& ivox)
{
  ASSERT(isSameSize(ivox));
  typename RealVoxel<T2>::const_iterator vitr = ivox.begin();
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr, ++vitr)
    *itr = clamped_convert<T>((*itr) + (*vitr));
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator*=(const RealVoxel<T2>& ivox)
{
  ASSERT(isSameSize(ivox));
  typename RealVoxel<T2>::const_iterator vitr = ivox.begin();
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr, ++vitr)
    *itr = clamped_convert<T>((*itr) + (*vitr));
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator/=(const RealVoxel<T2>& ivox)
{
  ASSERT(isSameSize(ivox));
  typename RealVoxel<T2>::const_iterator vitr = ivox.begin();
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr, ++vitr)
    *itr = clamped_convert<T>((*itr) + (*vitr));
  return *this;
}

// ######################################################################
template <class T> inline RealVoxel<T> &
RealVoxel<T>::operator=(const T& val)
{
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr)
    *itr = val;
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator+(const T2& val) const
{ return valueOp<RealVoxelAdd>(val); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator-(const T2& val) const
{ return valueOp<RealVoxelSub>(val); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator*(const T2& val) const
{ return valueOp<RealVoxelMult>(val); }

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::operator/(const T2& val) const
{ return valueOp<RealVoxelDiv>(val); }

// ######################################################################
template <class T> template <class T2, class MATH_OP>
inline RealVoxel<typename promote_trait<T,T2>::TP>
RealVoxel<T>::valueOp(const T2& val) const
{
  MATH_OP math_op;
  typedef typename promote_trait<T,T2>::TP TPRO;
  RealVoxel<TPRO> result(itsSizeX,itsSizeY,itsSizeZ,NO_INIT);
  typename RealVoxel<TPRO>::iterator ritr = result.beginw();
  for(typename RealVoxel<T>::const_iterator itr = begin(), stop = end();
      itr != stop; ++itr, ++ritr)
    *ritr = math_op.op((*itr),val);
  return result;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator+=(const T2& val)
{
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr)
    *itr = clamped_convert<T>((*itr) + val);
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator-=(const T2& val)
{
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr)
    *itr = clamped_convert<T>((*itr) - val);
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator*=(const T2& val)
{
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr)
    *itr = clamped_convert<T>((*itr) * val);
  return *this;
}

// ######################################################################
template <class T> template <class T2>
inline RealVoxel<T> &
RealVoxel<T>::operator/=(const T2& val)
{
  for (typename RealVoxel<T>::iterator itr = beginw(), stop = endw();
       itr != stop; ++itr)
    *itr = clamped_convert<T>((*itr) / val);
  return *this;
}

// ######################################################################
// ##### Image Slice Overloaded RealVoxel Operators
// ######################################################################

template <class T> template <class T2>
inline void
RealVoxel<T>::addImage(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelAdd>(axis,img); }

// ######################################################################
template <class T> template <class T2>
inline void
RealVoxel<T>::subtractImage(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelSub>(axis,img); }

// ######################################################################
template <class T> template <class T2>
inline void
RealVoxel<T>::subtractImageFrom(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelSubInv>(axis,img); }

// ######################################################################
template <class T> template <class T2>
inline void
RealVoxel<T>::multiplyImage(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelMult>(axis,img); }

// ######################################################################
template <class T> template <class T2>
inline void
RealVoxel<T>::divideImage(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelDiv>(axis,img); }

// ######################################################################
template <class T> template <class T2>
inline void
RealVoxel<T>::divideImageFrom(const imgvox::axis axis, const Image<T2>& img)
{ sliceOp<RealVoxelDivInv>(axis,img); }

// ######################################################################
// Operate on the VoxelImage using an Image Slice
template <class T> template <class T2, class MATH_OP>
inline void
RealVoxel<T>::sliceOp(const imgvox::axis axis, const Image<T2>& img)
{
  /* This method works very similar to getSlice and setSlice
  */

  // pre compute the image X by Y slice size for Z offset computation

  // This is the operator like a + b or a * b
  // it is used as a template so that all this stuff only has to be
  // put in once and we just define new ops for each math function
  MATH_OP math_op;

  const unsigned int XY = itsSizeX * itsSizeY;
  if(axis == imgvox::X_AXIS) // op slice by x-axis
  {
    ASSERT(img.getWidth()  == (int)itsSizeZ);
    ASSERT(img.getHeight() == (int)itsSizeY);

    for(unsigned int i = 0; i < itsSizeX; i++)
    {
      typename Image<T2>::const_iterator itr = img.begin();
      // There is no natural alignement. The data offset
      // into the array needs to be computed new each time
      for(unsigned int k = 0; k < itsSizeZ; k++)
      {
        const unsigned int xZ = i + k * XY;
        for(unsigned int j = 0; j < itsSizeY; j++)
        {
          const unsigned int pos = j * itsSizeX + xZ;
          beginw()[pos] = math_op.op(begin()[pos], *itr++);
        }
      }
    }
  }
  else if(axis == imgvox::Y_AXIS) // op slice by y-axis
  {
    ASSERT(img.getWidth()  == (int)itsSizeX);
    ASSERT(img.getHeight() == (int)itsSizeZ);

    for(unsigned int j = 0; j < itsSizeY; j++)
    {
      typename Image<T2>::const_iterator itr = img.begin();

      const unsigned int Y = j * itsSizeX;
      // We will need to compute the z coordinate in jumps
      for(unsigned int k = 0; k < itsSizeZ; k++)
      {
        const unsigned int yZ = Y + k * XY;
        const T *a            = &begin()[yZ];
        T *b                  = &beginw()[yZ];
        // The x axis alignment is natural and we can just
        // do a standard iterator copy
        for(unsigned int i = 0; i < itsSizeX; i++)
        {
          *b++ = math_op.op(*itr++,*a++);
        }
      }
    }
  }
  else if(axis == imgvox::Z_AXIS) // op slice by z-axis
  {
    ASSERT(img.getWidth()  == (int)itsSizeX);
    ASSERT(img.getHeight() == (int)itsSizeY);

    for(unsigned int k = 0; k < itsSizeZ; k++)
    {
      typename Image<T2>::const_iterator itr = img.begin();

      const unsigned int Z = k * XY;
      const T *a           = &begin()[Z];
      T *b                 = &beginw()[Z];
      // Each point is naturally aligned so we can do an easy copy
      // once we have computed the z alignment offset
      for(unsigned int p = 0; p < XY; p++)
      {
        *b++ = math_op.op(*itr++,*a++);
      }
    }
  }
  else
    LFATAL("Cannot get a slice without specifying axis for slice using -1");
}

// ######################################################################
// ##### Functions for testing/debugging only:
// ######################################################################

template <class T> inline
bool RealVoxel<T>::hasSameData(const RealVoxel<T>& b) const
{ return itsVox.hasSameData(b.itsVox); }

// ######################################################################
template <class T> inline
long RealVoxel<T>::refCount() const throw() { return itsVox.refCount(); }

// ######################################################################
template <class T> inline
bool RealVoxel<T>::isShared() const throw() { return itsVox.isShared(); }

// ######################################################################
// ##### Protected methods:
// ######################################################################

template <class T> inline
const ArrayData<T>& RealVoxel<T>::impl() const
{ return itsVox.get(); }

// ######################################################################
template <class T> inline
ArrayData<T>& RealVoxel<T>::uniq()
{ return itsVox.uniq(); }

// ######################################################################
// ######################################################################
// ##### Special Helper Classes
// ######################################################################
// ######################################################################
template <class T> inline RealVoxelAdd<T>::RealVoxelAdd() {}
template <class T> inline RealVoxelAdd<T>::~RealVoxelAdd() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelAdd<T>::op(const T& a, const T2& b) const
{ return a + b; }

// ######################################################################

template <class T> inline RealVoxelSub<T>::RealVoxelSub() {}
template <class T> inline RealVoxelSub<T>::~RealVoxelSub() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelSub<T>::op(const T& a, const T2& b) const
{ return a - b; }

// ######################################################################

template <class T> inline RealVoxelSubInv<T>::RealVoxelSubInv() {}
template <class T> inline RealVoxelSubInv<T>::~RealVoxelSubInv() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelSubInv<T>::op(const T& a, const T2& b) const
{ return b - a; }

// ######################################################################

template <class T> inline RealVoxelMult<T>::RealVoxelMult() {}
template <class T> inline RealVoxelMult<T>::~RealVoxelMult() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelMult<T>::op(const T& a, const T2& b) const
{ return a * b; }

// ######################################################################

template <class T> inline RealVoxelDiv<T>::RealVoxelDiv() {}
template <class T> inline RealVoxelDiv<T>::~RealVoxelDiv() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelDiv<T>::op(const T& a, const T2& b) const
{ return a / b; }

// ######################################################################

template <class T> inline RealVoxelDivInv<T>::RealVoxelDivInv() {}
template <class T> inline RealVoxelDivInv<T>::~RealVoxelDivInv() {}

template <class T> template <class T2>
inline typename promote_trait<T,T2>::TP
RealVoxelDivInv<T>::op(const T& a, const T2& b) const
{ return b / a; }

#endif

---