UT/UT_VoxelArray.h

Go to the documentation of this file.

/*
 * PROPRIETARY INFORMATION.  This software is proprietary to
 * Side Effects Software Inc., and is not to be reproduced,
 * transmitted, or disclosed in any way without written permission.
 *
 * Produced by:
 *      Jeff Lait
 *      Side Effects Software Inc
 *      477 Richmond Street West
 *      Toronto, Ontario
 *      Canada   M5V 3E7
 *      416-504-9876
 *
 * NAME:        UT_VoxelArray.h ( UT Library, C++)
 *
 * COMMENTS:
 *      This provides support for transparently tiled voxel arrays of data.
 *      The given type, T, should support normal arithmatic operations.
 *
 *      The created array has elements indexed from 0, ie: [0..xdiv-1].
 */

#ifndef __UT_VoxelArray__
#define __UT_VoxelArray__

#include "UT_API.h"
#include <SYS/SYS_Types.h>
#include "UT_Vector2.h"
#include "UT_Vector3.h"
#include "UT_Vector4.h"
#include "UT_PtrArray.h"
#include "UT_RefArray.h"
#include "UT_FilterType.h"
#include "UT_COW.h"
#include "UT_ThreadedAlgorithm.h"
#include <VM/VM_SIMD.h>

#include <boost/shared_ptr.hpp>

class UT_BoundingBox;
class UT_Filter;


static const int        TILEBITS = 4;
static const int        TILESIZE = 1 << TILEBITS;
static const int        TILEMASK = TILESIZE-1;

///
/// Behaviour of out of bound reads.
///
enum UT_VoxelBorderType
{
    UT_VOXELBORDER_CONSTANT,
    UT_VOXELBORDER_REPEAT,
    UT_VOXELBORDER_STREAK
};

template <typename T> class UT_VoxelTile;
template <typename T, bool DoRead, bool DoWrite, bool TestForWrite> class UT_VoxelProbe;
template <typename T> class UT_VoxelProbeCube;
template <typename T> class UT_VoxelProbeFace;

///
/// UT_VoxelTileCompress
///
/// A compression engine for UT_VoxelTiles of a specific type.  This
/// is a verb class which is invoked from the voxeltile class.
///
template <typename T>
class UT_VoxelTileCompress
{
public:
                 UT_VoxelTileCompress() {}
    virtual     ~UT_VoxelTileCompress() {}
    
    /// Attempts to write data directly to the compressed tile.
    /// Returns false if not possible.
    virtual bool writeThrough(UT_VoxelTile<T> &tile, 
                              int x, int y, int z, T t) const = 0;

    /// Reads directly from the compressed data.
    /// Cannot alter the tile in any way because it must be threadsafe.
    virtual T    getValue(const UT_VoxelTile<T> &tile,
                          int x, int y, int z) const = 0;

    /// Attempts to compress the data according to the given tolerance.
    /// If succesful, returns true.
    virtual bool tryCompress(UT_VoxelTile<T> &tile,
                          fpreal tol, T min, T max) const = 0;

    /// Returns the length in bytes of the data in the tile.
    /// It must be at least one byte long.
    virtual int  getDataLength(const UT_VoxelTile<T> &tile) const = 0;

    /// Determines the min & max values of the tile.  A default
    /// implementation uses getValue() on all voxels.
    virtual void findMinMax(const UT_VoxelTile<T> &tile, T &min, T &max) const;

    /// Does this engine support saving and loading?
    virtual bool canSave() const { return false; }
    virtual void save(ostream &os, const UT_VoxelTile<T> &tile) const {}
    virtual void load(UT_IStream &is, UT_VoxelTile<T> &tile) const {}

    /// Returns the unique name of this compression engine so
    /// we can look up engines by name (the index of the compression
    /// engine is assigned at load time so isn't constant)
    virtual const char *getName() = 0;
};

UT_API UT_PtrArray<UT_VoxelTileCompress<fpreal16> *> &UTvoxelTileGetCompressionEngines(fpreal16 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<fpreal32> *> &UTvoxelTileGetCompressionEngines(fpreal32 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<fpreal64> *> &UTvoxelTileGetCompressionEngines(fpreal64 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<uint8> *> &UTvoxelTileGetCompressionEngines(uint8 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<int8> *> &UTvoxelTileGetCompressionEngines(int8 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<int16> *> &UTvoxelTileGetCompressionEngines(int16 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<int32> *> &UTvoxelTileGetCompressionEngines(int32 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<int64> *> &UTvoxelTileGetCompressionEngines(int64 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<UT_Vector2> *> &UTvoxelTileGetCompressionEngines(UT_Vector2 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<UT_Vector3> *> &UTvoxelTileGetCompressionEngines(UT_Vector3 *dummy);
UT_API UT_PtrArray<UT_VoxelTileCompress<UT_Vector4> *> &UTvoxelTileGetCompressionEngines(UT_Vector4 *dummy);

#define DEFINE_STD_FUNC(TYPE)                                   \
inline void                                                     \
UTvoxelTileExpandMinMax(TYPE v, TYPE &min, TYPE &max)           \
{                                                               \
    if (v < min)                                                \
        min = v;                                                \
    else if (v > max)                                           \
        max = v;                                                \
}                                                               \
                                                                \
inline fpreal                                                   \
UTvoxelTileDist(TYPE a, TYPE b)                                 \
{                                                               \
    return (fpreal) SYSabs(a - b);                              \
}

DEFINE_STD_FUNC(fpreal16)
DEFINE_STD_FUNC(fpreal32)
DEFINE_STD_FUNC(fpreal64)
DEFINE_STD_FUNC(uint8)
DEFINE_STD_FUNC(int8)
DEFINE_STD_FUNC(int16)
DEFINE_STD_FUNC(int32)
DEFINE_STD_FUNC(int64)

#undef DEFINE_STD_FUNC

inline void 
UTvoxelTileExpandMinMax(UT_Vector2 v, UT_Vector2 &min, UT_Vector2 &max)
{
    min.x() = SYSmin(v.x(), min.x());
    max.x() = SYSmax(v.x(), max.x());

    min.y() = SYSmin(v.y(), min.y());
    max.y() = SYSmax(v.y(), max.y());
}

inline void 
UTvoxelTileExpandMinMax(UT_Vector3 v, UT_Vector3 &min, UT_Vector3 &max)
{
    min.x() = SYSmin(v.x(), min.x());
    max.x() = SYSmax(v.x(), max.x());

    min.y() = SYSmin(v.y(), min.y());
    max.y() = SYSmax(v.y(), max.y());

    min.z() = SYSmin(v.z(), min.z());
    max.z() = SYSmax(v.z(), max.z());
}

inline void 
UTvoxelTileExpandMinMax(UT_Vector4 v, UT_Vector4 &min, UT_Vector4 &max)
{
    min.x() = SYSmin(v.x(), min.x());
    max.x() = SYSmax(v.x(), max.x());

    min.y() = SYSmin(v.y(), min.y());
    max.y() = SYSmax(v.y(), max.y());

    min.z() = SYSmin(v.z(), min.z());
    max.z() = SYSmax(v.z(), max.z());

    min.w() = SYSmin(v.w(), min.w());
    max.w() = SYSmax(v.w(), max.w());
}

inline fpreal
UTvoxelTileDist(const UT_Vector2 &a, const UT_Vector2 &b)
{
    return SYSabs(a.x() - b.x()) + SYSabs(a.y() - b.y());
}

inline fpreal
UTvoxelTileDist(const UT_Vector3 &a, const UT_Vector3 &b)
{
    return SYSabs(a.x() - b.x()) + SYSabs(a.y() - b.y())
         + SYSabs(a.z() - b.z());
}

inline fpreal
UTvoxelTileDist(const UT_Vector4 &a, const UT_Vector4 &b)
{
    return SYSabs(a.x() - b.x()) + SYSabs(a.y() - b.y())
         + SYSabs(a.z() - b.z()) + SYSabs(a.w() - b.w());
}

///
/// UT_VoxelTile
///
/// A UT_VoxelArray is composed of a number of these tiles.  This is
/// done for two reasons:
/// 1) Increased memory locality when processing neighbouring points.
/// 2) Ability to compress or page out unneeded tiles.
/// Currently, the only special ability is the ability to create constant
/// tiles.
///
/// To the end user of the UT_VoxelArray, the UT_VoxelTile should be
/// usually transparent.  The only exception may be if they want to do
/// a FOR_ALL_TILES in order to ensure an optimal traversal order.
///
template <typename T>
class UT_VoxelTile
{
public:
             UT_VoxelTile(int xres, int yres, int zres);
    virtual ~UT_VoxelTile();

    // Copy constructor:
             UT_VoxelTile(const UT_VoxelTile<T> &src);

    // Assignment operator:
    const UT_VoxelTile<T> &operator=(const UT_VoxelTile<T> &src);

    enum CompressionType
    {
        COMPRESS_RAW,
        COMPRESS_RAWFULL,
        COMPRESS_CONSTANT,
        COMPRESS_ENGINE
    };

    /// Fetch a given local value.  (x,y,z) should be local to
    /// this tile.
    T            operator()(int x, int y, int z) const
                 {
                     UT_ASSERT_P(x >= 0 && y >= 0 && z >= 0);
                     UT_ASSERT_P(x < myRes[0] && y < myRes[1] && z < myRes[2]);

                     switch (myCompressionType)
                     {
                         case COMPRESS_RAW:
                             return ((T *)myData)[
                                 ((z * myRes[1]) + y) * myRes[0] + x ];

                         case COMPRESS_CONSTANT:
                             return *(T *)myData;

                         case COMPRESS_RAWFULL:
                             return ((T *)myData)[
                                 ((z * TILESIZE) + y) * TILESIZE + x ];
                     }

                     // By default use the compression engine.
                     UT_VoxelTileCompress<T>    *engine;

                     engine = getCompressionEngine(myCompressionType);
                     return engine->getValue(*this, x, y, z);
                 }

    /// Lerps two numbers, templated to work with T.
    inline static T      lerpValues(T v1, T v2, fpreal32 bias)
    {
        return v1 + (v2 - v1) * bias;
    }

    /// Does a trilinear interpolation.  x,y,z should be local to this
    /// as should x+1, y+1, and z+1.  fx-fz should be 0..1.
    T            lerp(int x, int y, int z, float fx, float fy, float fz) const;
#if 0
    /// MSVC can't handle aligned parameters after the third so
    /// frac must come first.
    T            lerp(v4uf frac, int x, int y, int z) const;
#endif

    /// Returns a cached line to our internal data, at local address x,y,z.
    /// cacheline is a caller allocated structure to fill out if we have
    /// to decompress.  If forcecopy isn't set and we can, the result may
    /// be an internal pointer.  stride is set to the update for moving one
    /// x position in the cache.
    /// strideofone should be set to true if you want to prevent 0 stride
    /// results for constant tiles.
    T           *fillCacheLine(T *cacheline, int &stride, int x, int y, int z, bool forcecopy, bool strideofone) const;

    /// Fills a cache line from an external buffer into our own data.
    void         writeCacheLine(T *cacheline, int y, int z);

    /// The setData is intentionally seperate so we can avoid
    /// expanding constant data when we write the same value to it.
    void         setValue(int x, int y, int z, T t);

    /// Finds the minimum and maximum T values
    void         findMinMax(T &min, T &max) const;

    /// Returns if this tile is constant.
    bool         isConstant() const 
                    { return myCompressionType == COMPRESS_CONSTANT; }

    /// Returns if this tile is in raw format.
    bool         isRaw() const
                    { return myCompressionType == COMPRESS_RAW; }

    /// Returns if this tile is in raw full format.
    bool         isRawFull() const
                    { return myCompressionType == COMPRESS_RAWFULL; }

    /// Attempts to compress this tile.  Returns true if any
    /// compression performed.
    bool         tryCompress(fpreal tolerance);

    /// Turns this tile into a constant tile of the given value.
    void         makeConstant(T t);

    /// Turns a compressed tile into a raw tile.
    void         uncompress();

    /// Turns a tile into a raw full tile.
    void         uncompressFull();

    /// Returns the raw full data of the tile.
    T           *rawFullData()
                 {
                     uncompressFull();
                     return (T *)myData;
                 }

    /// Read the current resolution.
    int          xres() const { return myRes[0]; }
    int          yres() const { return myRes[1]; }
    int          zres() const { return myRes[2]; }

    int          getRes(int dim) const { return myRes[dim]; }

    /// Returns the amount of memory used by this tile.
    int64        getMemoryUsage() const;

    /// Returns the amount of data used by the tile myData pointer.
    int64        getDataLength() const;

    /// A routine used by filtered evaluation to accumulated a partial
    /// filtered sum in this tile.
    ///         pstart, pend - voxel bounds (in UT_VoxelArray coordinates)
    ///         weights - weight array
    ///         start - UT_VoxelArray coordinates at [0] in the weight array
    void         weightedSum(int pstart[3], int pend[3],
                             float *weights[3], int start[3],
                             T &result);

    /// Designed to be specialized according to T
    
    /// Update min & max to encompass T itself.
    static void  expandMinMax(T v, T &min, T &max)
    {
        UTvoxelTileExpandMinMax(v, min, max);
    }
    
    /// Return the "distance" of a & b.  This is used for
    /// tolerance checks on equality comparisons.
    static fpreal        dist(T a, T b)
    {
        return UTvoxelTileDist(a, b);
    }
    
    static void  registerCompressionEngine(UT_VoxelTileCompress<T> *engine);

    // Returns the index of the bound compression engine.
    static int   lookupCompressionEngine(const char *name);
    // Given an index, gets the compression engine.
    static UT_VoxelTileCompress<T> *getCompressionEngine(int index);

    /// Saves this tile's data, in compressed form.  May uncompress itself
    /// if the engine doesn't support saving.
    void                save(ostream &os); 

    /// Loads tile data.  Uses the compression index to map the saved
    /// compression types into the correct loading compression types.
    void                load(UT_IStream &is, const UT_IntArray &compression);

    /// Stores a list of compresson engines to os.
    static void          saveCompressionTypes(ostream &os);

    /// Builds a translation table from the given stream's compression types
    /// into our own valid compression types.
    static void          loadCompressionTypes(UT_IStream &is, UT_IntArray &compressions);
    
protected:
    // Attempts to set the value to the native compressed format
    // Some compression types allow some values to be written
    // without decompression. Eg, you can write to a constant tile
    // the tile's own value without decompression.
    // If this returns true, t has been written.
    bool         writeThrough(int x, int y, int z, T t);

public:
    // This is only public so the compression engines can get to it.
    // It is blind data, do not alter!
    void        *myData;
private:

    /// Resolutions.
    int8         myRes[3];

    /// Am I a constant tile?
    int8         myCompressionType;
    
    static UT_PtrArray<UT_VoxelTileCompress<T> *> &getCompressionEngines()
    {
        return UTvoxelTileGetCompressionEngines((T *) 0);
    }

    friend class UT_VoxelTileCompress<T>;
    template <typename S, bool DoWrite, bool DoRead, bool TestForWrites>
    friend class UT_VoxelProbe;
};

///
/// UT_VoxelArray
///
/// This provides data structure to hold a three dimmensional array
/// of data.  The data should be some simple arithmatic type, such
/// as uint8, fpreal16, or UT_Vector3.
///
/// Some operations, such as gradiants, may make less sense with uint8.
/// 
template <typename T>
class UT_VoxelArray
{
public:
             UT_VoxelArray();
    virtual ~UT_VoxelArray();

    /// Copy constructor:
             UT_VoxelArray(const UT_VoxelArray<T> &src);

    /// Assignment operator:
    const UT_VoxelArray<T> &operator=(const UT_VoxelArray<T> &src);

    /// This sets the voxelarray to have the given resolution, resetting
    /// all elements to 0.
    void         size(int xres, int yres, int zres);

    /// This will ensure this voxel array matches the given voxel array
    /// in terms of dimensions & border conditions.  It may invoke
    /// a size() and hence reset the field to 0.
    void         match(const UT_VoxelArray<T> &src);

    int          getXRes() const { return myRes[0]; }
    int          getYRes() const { return myRes[1]; }
    int          getZRes() const { return myRes[2]; }
    int          getRes(int axis) const { return myRes[axis]; }

    /// Return the amount of memory used by this array.
    int64        getMemoryUsage() const;

    /// Sets this voxel array to the given constant value.  All tiles
    /// are turned into constant tiles.
    void         constant(T t);

    /// If this voxel array is all constant tiles, returns true.
    /// The optional pointer is initialized to the constant value iff
    /// the array is constant.  (Note by constant we mean made of constant
    /// tiles of the same value - if some tiles are uncompressed but
    /// constant, it will still return false)
    bool         isConstant(T *cval = 0) const;

    /// This convience function lets you sample the voxel array.
    /// pos is in the range [0..1]^3.
    /// T value trilinearly interpolated.   Edges are determined by the border
    /// mode.
    /// The cells are sampled at the center of the voxels.
    T            operator()(UT_Vector3 pos) const;
#if 0
    T            operator()(v4uf pos) const;
#endif

    /// Filtered evaluation of the voxel array.  This operation should
    /// exhibit the same behavior as IMG3D_Channel::evaluate.  Borders
    /// evaluate to 0 and the border type is ignored.  radius is 
    /// measured in voxels.
    T            evaluate(const UT_Vector3 &pos, const UT_Filter &filter,
                          fpreal radius) const;

    /// Fills this by resampling the given voxel array.  Uses evaluate
    /// internally, so the same caveats apply regarding border behaviour.
    void         resample(const UT_VoxelArray<T> &src,
                          UT_FilterType filtertype = UT_FILTER_POINT);

    /// Converts a 3d position in range [0..1]^3 into the closest
    /// index value.
    /// Returns false if the resulting index was out of range.  The index
    /// will still be set.
    bool         posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const;
    /// Convertes a 3d position in [0..1]^3 into the equivalent in
    /// the integer cell space.  Does not clamp to the closest value.
    bool         posToIndex(UT_Vector3 pos, UT_Vector3 &ipos) const;
    /// Converts an index into a position.
    /// Returns false if the source index was out of range, in which case
    /// pos will be outside [0..1]^3
    bool         indexToPos(int x, int y, int z, UT_Vector3 &pos) const; 

    /// Clamps the given x, y, and z values to lie inside the valid index
    /// range.
    void         clampIndex(int &x, int &y, int &z) const
                 {
                     x = SYSclamp(x, 0, myRes[0]-1);
                     y = SYSclamp(y, 0, myRes[1]-1);
                     z = SYSclamp(z, 0, myRes[2]-1);
                 }

    /// Returns true if the given x, y, z values lie inside the valid index.
    bool         isValidIndex(int x, int y, int z) const
                 {
                     return x >= 0 && x < myRes[0] &&
                            y >= 0 && y < myRes[1] &&
                            z >= 0 && z < myRes[2];
                 }
    
    /// This allows you to read & write the raw data.
    /// Out of bound reads are illegal.
    T            operator()(int x, int y, int z) const
                 {
                     UT_ASSERT_P(isValidIndex(x, y, z));
                     return (*getTile(x >> TILEBITS,
                                      y >> TILEBITS,
                                      z >> TILEBITS))
                         (x & TILEMASK, y & TILEMASK, z & TILEMASK);
                 }
    void         setValue(int x, int y, int z, T t)
                 {
                     UT_ASSERT_P(isValidIndex(x, y, z));
                     getTile(x >> TILEBITS,
                             y >> TILEBITS,
                             z >> TILEBITS)->setValue(
                                 x & TILEMASK, y & TILEMASK, z & TILEMASK, t);
                 }

    /// This will clamp the bounds to fit within the voxel array,
    /// using the border type to resolve out of range values.
    T            getValue(int x, int y, int z) const
                 {
                     // Verify our voxel array is non-empty.
                     if (!myRes[0] && !myRes[1] && !myRes[2])
                     {
                         UT_ASSERT(!"getValue on empty voxel array!");
                         return myBorderValue;
                     }

                     // Verify we are within range.  Behaviour on out of
                     // range depends on border type.
                     switch (myBorderType)
                     {
                     case UT_VOXELBORDER_CONSTANT:
                         if (!isValidIndex(x, y, z))
                             return myBorderValue;
                         break;

                     case UT_VOXELBORDER_REPEAT:
                         if (x < 0 || x >= myRes[0])
                         {
                             x %= myRes[0];
                             if (x < 0)
                                 x += myRes[0];
                         }
                         if (y < 0 || y >= myRes[1])
                         {
                             y %= myRes[1];
                             if (y < 0)
                                 y += myRes[1];
                         }
                         if (z < 0 || z >= myRes[2])
                         {
                             z %= myRes[2];
                             if (z < 0)
                                 z += myRes[2];
                         }
                         break;

                     case UT_VOXELBORDER_STREAK:
                         clampIndex(x, y, z);
                         break;
                     }

                     // It is now within bounds, do normal fetch.
                     return (*this)(x, y, z);
                 }

    void         setBorder(UT_VoxelBorderType type, T t);
    UT_VoxelBorderType getBorder() const { return myBorderType; }
    T            getBorderValue() const { return myBorderValue; }

    /// This tries to compress or collapse each tile.  This can
    /// be expensive (ie, converting a tile to constant), so
    /// should be saved until modifications are complete.
    void         collapseAllTiles();

    /// Uncompresses all tiles into non-constant tiles.  Useful
    /// if you have a multithreaded algorithm that may need to
    /// both read and write, if you write to a collapsed tile
    /// while someone else reads from it, bad stuff happens.
    /// Instead, you can expandAllTiles.  This may have serious
    /// consequences in memory use, however.
    void         expandAllTiles();
    
    /// The direct tile access methods are to make TBF writing a bit
    /// more efficient.
    UT_VoxelTile<T>     *getTile(int tx, int ty, int tz) const
                         { return myTiles[xyzTileToLinear(tx, ty, tz)]; }
    UT_VoxelTile<T>     *getLinearTile(int idx) const
                         { return myTiles[idx]; }
    void                 linearTileToXYZ(int idx, int &x, int &y, int &z) const
                         {
                             x = idx % myTileRes[0];
                             idx -= x;
                             idx /= myTileRes[0];
                             y = idx % myTileRes[1];
                             idx -= y;
                             idx /= myTileRes[1];
                             z = idx;
                         }
    int                  xyzTileToLinear(int x, int y, int z) const
                         { return (z * myTileRes[1] + y) * myTileRes[0] + x; }

    /// Number of tiles along that axis.  Not to be confused with
    /// the resolution of the individual tiles.
    int                  getTileRes(int dim) const { return myTileRes[dim]; }
    int                  numTiles() const 
                 { return myTileRes[0] * myTileRes[1] * myTileRes[2]; }
    int                  numVoxels() const 
                 { return myRes[0] * myRes[1] * myRes[2]; }

    fpreal               getCompressionTolerance() const
                 { return myCompressionTolerance; }
    void                 setCompressionTolerance(fpreal tol)
                 { myCompressionTolerance = tol; }

    /// Saves only the data of this array to the given stream.
    /// To reload it you will have to have a matching array in tiles
    /// dimensions and size.
    /// Not consta as unsupported compression engines will be uncompressed.
    void                 saveData(ostream &os);

    /// Load an array, requires you have already size()d this array.
    void                 loadData(UT_IStream &is);
    
private:
    THREADED_METHOD3(UT_VoxelArray<T>, numTiles() > 1,
                     resamplethread,
                     const UT_VoxelArray<T> &, src,
                     const UT_Filter *, filter,
                     float, radius)
    void         resamplethreadPartial(const UT_VoxelArray<T> &src,
                                       const UT_Filter *filter,
                                       float radius,
                                       const UT_JobInfo &info);

    void         deleteVoxels();
    
    /// Number of elements in each dimension.
    int                  myRes[3];

    /// Number of tiles in each dimension.
    int                  myTileRes[3];

    /// Compression tolerance for lossy compression.
    fpreal               myCompressionTolerance;

    /// Double dereferenced so we can theoritically resize easily.
    UT_VoxelTile<T>     **myTiles;
    
    /// Outside values get this if constant borders are used
    T                    myBorderValue;
    UT_VoxelBorderType   myBorderType;
};


///
/// UT_VoxelMipMap
///
/// This provides a mip-map type structure for a voxel array.
/// It manages the different levels of voxels arrays that are needed.
/// You can create different types of mip maps: average, maximum, etc,
/// which can allow different tricks.
/// Each level is one half the previous level, rounded up.
/// Out of bound voxels are ignored from the lower levels.
///
template <typename T>
class UT_VoxelMipMap
{
public:
    /// The different types of functions that can be used for
    /// constructing a mip map.
    enum mipmaptype { MIPMAP_MAXIMUM=0, MIPMAP_AVERAGE=1, MIPMAP_MINIMUM=2 };
    
             UT_VoxelMipMap();
    virtual ~UT_VoxelMipMap();

    /// Copy constructor.
             UT_VoxelMipMap(const UT_VoxelMipMap<T> &src);

    /// Assignment operator:
    const UT_VoxelMipMap<T> &operator=(const UT_VoxelMipMap<T> &src);

    /// Builds from a given voxel array.  The ownership flag determines
    /// if we gain ownership of the voxel array and should delete it.
    /// In any case, the new levels are owned by us.
    void                 build(UT_VoxelArray<T> *baselevel, 
                               mipmaptype function);

    /// Same as above but construct mipmaps simultaneously for more than
    /// one function.  The order of the functions will correspond to the
    /// order of the data values passed to the traversal callback.
    void                 build(UT_VoxelArray<T> *baselevel, 
                               const UT_RefArray<mipmaptype> &functions);

    /// This does a top down traversal of the implicit octtree defined
    /// by the voxel array.  Returning false will abort that
    /// branch of the octree.
    /// The bounding box given is in cell space and is an exclusive
    /// box of the included cells (ie: (0..1)^3 means just cell 0,0,0)
    /// Note that each bounding box will not be square, unless you
    /// have the good fortune of starting with a power of 2 cube.
    /// The boolean goes true when the the callback is invoked on a
    /// base level.
    void                 traverseTopDown(bool (*function)(
                                                const UT_RefArray<T> &,
                                                const UT_BoundingBox &,
                                                bool,
                                                void *),
                                         void *data) const;

    /// Return the amount of memory used by this mipmap.
    int64                getMemoryUsage() const;

private:
    void                 doTraverse(int x, int y, int z, int level,
                                    bool (*function)(const UT_RefArray<T> &,
                                                     const UT_BoundingBox &,
                                                     bool,
                                                     void *),
                                    void *data) const;
    
    void                 initializePrivate();
    void                 destroyPrivate();

protected:
    T                    mixValues(T t1, T t2, mipmaptype function) const;
    
    /// This stores the base most level that was provided
    /// externally.
    UT_VoxelArray<T>    *myBaseLevel;
    /// If true, we will delete the base level when we are done.
    bool                 myOwnBase;

    /// Tracks the number of levels which we used to represent
    /// this hierarchy.
    int                  myNumLevels;
    /// The array of VoxelArrays, one per level.
    /// myLevels[0] is a 1x1x1 array. Each successive layer is twice
    /// as big in each each dimension.  However, every layer is clamped
    /// against the resolution of the base layer.
    /// We own all these layers.
    UT_PtrArray<UT_VoxelArray<T> **>     myLevels;
};


/// Iterator for Voxel Arrays
///
/// This class eliminates the need for having
/// for (z = 0; z < zres; z++)
///   ...
///     for (x = 0; x < xres; x++)
/// loops everywhere.
///
/// Note that the order of iteration is undefined!  (The actual order is
/// to complete each tile in turn, thereby hopefully improving cache
/// coherency)
///
/// It is safe to write to the voxel array while this iterator is active.
/// It is *not* safe to resize the voxel array (or destroy it)
///
/// The iterator is similar in principal to an STL iterator, but somewhat
/// simpler. The classic STL loop
///   for ( it = begin(); it != end(); ++it )
/// is done using
///   for ( it.rewind(); !it.atEnd(); it.advance() )
/// 
template <typename T>
class UT_VoxelArrayIterator
{
public:
             UT_VoxelArrayIterator();
             UT_VoxelArrayIterator(UT_VoxelArray<T> *vox);
             UT_VoxelArrayIterator(UT_COWReadHandle<UT_VoxelArray<T> > handle);
    virtual ~UT_VoxelArrayIterator();       

    void        setArray(UT_VoxelArray<T> *vox)
                {
                    myCurTile = -1;
                    myHandle.resetHandle();
                    myArray = vox;
                    // Reset the range
                    setPartialRange(0, 1);
                }

    /// Iterates over the array pointed to by the handle.  Only
    /// supports read access during the iteration as it does
    /// a read lock.
    void        setHandle(UT_COWReadHandle<UT_VoxelArray<T> > handle)
                {
                    myHandle = handle;
                    // Ideally we'd have a separate const iterator
                    // from our non-const iterator so this would
                    // only be exposed in the const version.
                    myArray = const_cast<UT_VoxelArray<T> *>(&*myHandle);

                    // Reset our range.
                    myCurTile = -1;
                    setPartialRange(0, 1);
                }


    /// Restricts this iterator to only run over a subset
    /// of the tiles.  The tiles will be divided into approximately
    /// numrange equal groups, this will be the idx'th.
    /// The resulting iterator may have zero tiles.
    void        setPartialRange(int idx, int numranges);

    /// Restricts this iterator to the tiles that intersect
    /// the given bounding box of voxel coordinates.
    /// Note that this will not be a precise restriction as
    /// each tile is either included or not.
    /// You should setPartialRange() after setting the bbox range
    /// The bounding box is on the [0..1]^3 range.
    void        restrictToBBox(const UT_BoundingBox &bbox);
    /// The [xmin, xmax] are inclusive and measured in voxels.
    void        restrictToBBox(int xmin, int xmax,
                               int ymin, int ymax,
                               int zmin, int zmax);

    /// Resets the iterator to point to the first voxel.
    void        rewind();

    /// Returns true if we have iterated over all of the voxels.
    bool        atEnd() const
                { return myCurTile < 0; }

    /// Advances the iterator to point to the next voxel.
    void        advance()
                {
                    // We try to advance each axis, rolling over to the next.
                    // If we exhaust this tile, we call advanceTile.
                    myPos[0]++;
                    myTileLocalPos[0]++;
                    if (myTileLocalPos[0] >= myTileSize[0])
                    {
                        // Wrapped in X.
                        myPos[0] -= myTileLocalPos[0];
                        myTileLocalPos[0] = 0;

                        myPos[1]++;
                        myTileLocalPos[1]++;
                        if (myTileLocalPos[1] >= myTileSize[1])
                        {
                            // Wrapped in Y.
                            myPos[1] -= myTileLocalPos[1];
                            myTileLocalPos[1] = 0;

                            myPos[2]++;
                            myTileLocalPos[2]++;
                            if (myTileLocalPos[2] >= myTileSize[2])
                            {
                                // Wrapped in Z!  Finished this tile!
                                advanceTile();
                            }
                        }
                    }
                }

    /// Retrieve the current location of the iterator.
    int         x() const { return myPos[0]; }
    int         y() const { return myPos[1]; }
    int         z() const { return myPos[2]; }
    int         idx(int idx) const { return myPos[idx]; }
    
    /// Retrieves the value that we are currently pointing at.
    /// This is faster than an operator(x,y,z) as we already know
    /// our current tile and that bounds checking isn't needed.
    T            getValue() const
                 {
                     UT_ASSERT_P(myCurTile >= 0);

                     UT_VoxelTile<T>    *tile;

                     tile = myArray->getLinearTile(myCurTile);
                     return (*tile)(myTileLocalPos[0],
                                    myTileLocalPos[1],
                                    myTileLocalPos[2]);
                 }
    
    /// Sets the voxel we are currently pointing to the given value.
    void         setValue(T t) const
                 {
                     UT_ASSERT_P(myCurTile >= 0);

                     UT_VoxelTile<T>    *tile;

                     tile = myArray->getLinearTile(myCurTile);

                     tile->setValue(myTileLocalPos[0],
                                    myTileLocalPos[1],
                                    myTileLocalPos[2], t);
                 }

    /// Returns true if the tile we are currently in is a constant tile.
    bool         isTileConstant() const
                 {
                     UT_ASSERT_P(myCurTile >= 0);

                     UT_VoxelTile<T>    *tile;

                     tile = myArray->getLinearTile(myCurTile);
                     return tile->isConstant();
                 }

    /// Returns true if we are at the start of a new tile.
    bool         isStartOfTile() const
                 { return !(myTileLocalPos[0] ||
                            myTileLocalPos[1] ||
                            myTileLocalPos[2]); }

    /// Returns the VoxelTile we are currently processing
    UT_VoxelTile<T>     *getTile() const
                         {
                             UT_ASSERT_P(myCurTile >= 0);
                             return myArray->getLinearTile(myCurTile);
                         }

    /// Advances the iterator to point to the next tile.  Useful if the
    /// constant test showed that you didn't need to deal with this one.
    void         advanceTile();

    /// Sets a flag which causes the iterator to tryCompress()
    /// tiles when it is done with them.  
    bool         getCompressOnExit() const { return myShouldCompressOnExit; }
    void         setCompressOnExit(bool shouldcompress)
                 { myShouldCompressOnExit = shouldcompress; }

protected:
    /// The array we belong to.
    UT_VoxelArray<T>    *myArray;
    /// The handle that we have locked to get our array.  It is null
    /// by default which makes the lock/unlock nops.
    UT_COWReadHandle<UT_VoxelArray<T> > myHandle;
    
    /// Absolute index into voxel array.
    int          myPos[3];

    /// Flag determining if we should compress tiles whenever we
    /// advance out of them.
    bool         myShouldCompressOnExit;

    bool         myUseTileList;
    UT_IntArray  myTileList;

public:
    /// Our current linear tile idx.  A value of -1 implies at end.
    int          myCurTile;

    /// Our current index into the tile list
    int          myCurTileListIdx;

    /// Our start & end tiles for processing a subrange.
    /// The tile range is half open [start, end)
    int          myTileStart, myTileEnd;

    /// Which tile we are as per tx,ty,tz rather than linear index.
    int          myTilePos[3];

    /// Our position within the current tile.
    int          myTileLocalPos[3];

    /// The size of the current tile
    int          myTileSize[3];
};

/// Probe for Voxel Arrays
///
/// This class is designed to allow for efficient evaluation
/// of aligned indices of a voxel array, provided the voxels are iterated
/// in a tile-by-tile, x-inner most, manner.
///
/// This class will create a local copy of the voxel data where needed,
/// uncompressing the information once for every 16 queries.  It will
/// also create an aligned buffer so you can safely use v4uf on fpreal32
/// data.
///
/// For queries where you need surrounding values, the prex and postx can
/// specify padding on the probe.  prex should be -1 to allow reading
/// -1 offset, postx 1 to allow reading a 1 offset.
///

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
class UT_VoxelProbe
{
public:
                 UT_VoxelProbe();
                 UT_VoxelProbe(UT_VoxelArray<T> *vox, int prex = 0, int postx = 0);
    virtual     ~UT_VoxelProbe();

    void         setArray(UT_VoxelArray<T> *vox, int prex = 0, int postx = 0);

    inline T     getValue() const 
    { 
        return *myCurLine; 
    }
    inline T     getValue(int offset) const 
    { 
        return myCurLine[myStride*offset]; 
    }

    inline void  setValue(T value) 
    {
        UT_ASSERT_P(DoWrite);
        *myCurLine = value;
        if (TestForWrites)
            myDirty = true;
    }
                

    /// Resets where we currently point to.
    /// Returns true if we had to reset our cache line.  If we didn't,
    /// and you have multiple probes acting in-step, you can just
    /// advanceX() the other probes
    template <typename S> 
    bool         setIndex(UT_VoxelArrayIterator<S> &vit)
                 { return setIndex(vit.x(), vit.y(), vit.z()); }

    bool         setIndex(int x, int y, int z);

    /// Blindly advances our current pointer.
    inline void  advanceX() 
    { 
        myCurLine += myStride; 
        myX++; 
        UT_ASSERT_P(myX < myMaxValidX); 
    }

    /// Adjusts our current pointer to the given absolute location,
    /// assumes the new value is inside our valid range.
    inline void  resetX(int x) 
    { 
        myCurLine += myStride * (x - myX); 
        myX = x; 
        UT_ASSERT_P(myX < myMaxValidX && myX >= myMinValidX); 
    }

protected:
    void         reloadCache(int x, int y, int z);

    void         writeCacheLine();

    void         buildConstantCache(T value);

    T           *myCurLine;
    /// myCacheLine[0] is the start of the cache line, so -1 would be
    /// the first pre-rolled value
    T           *myCacheLine;
    /// Where we actually allocated our cache line, aligned to 4x multiple
    /// to ensure SSE compatible.
    T           *myAllocCacheLine;

    int          myX, myY, myZ;
    int          myPreX, myPostX;
    int          myStride;
    bool         myForceCopy;
    /// Half inclusive [,) range of valid x queries for current cache.
    int          myMinValidX, myMaxValidX;

    /// Determines if we have anything to write back, only
    /// valid if TestForWrites is enabled.
    bool         myDirty;

    UT_VoxelArray<T>    *myArray;

    friend class UT_VoxelProbeCube<T>;
    friend class UT_VoxelProbeFace<T>;
};

///
/// The vector probe is three normal probes into separate voxel arrays
/// making it easier to read and write to aligned vector fields.
/// If the vector field is face-centered, see the UT_VoxelProbeFace.
///
template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
class UT_VoxelVectorProbe
{
public:
                 UT_VoxelVectorProbe()
                 { }
                 UT_VoxelVectorProbe(UT_VoxelArray<T> *vx, UT_VoxelArray<T> *vy, UT_VoxelArray<T> *vz)
                 { setArray(vx, vy, vz); }
    virtual     ~UT_VoxelVectorProbe()
                 {}

    void         setArray(UT_VoxelArray<T> *vx, UT_VoxelArray<T> *vy, UT_VoxelArray<T> *vz)
    {
        myLines[0].setArray(vx);
        myLines[1].setArray(vy);
        myLines[2].setArray(vz);
    }

    inline UT_Vector3    getValue() const 
    { 
        return UT_Vector3(myLines[0].getValue(), myLines[1].getValue(), myLines[2].getValue()); 
    }
    inline T     getValue(int axis) const 
    { 
        return myLines[axis].getValue(); 
    }

    inline void  setValue(const UT_Vector3 &v) 
    {
        myLines[0].setValue(v.x());
        myLines[1].setValue(v.y());
        myLines[2].setValue(v.z());
    }
                
    inline void  setComponent(int axis, T val) 
    {
        myLines[axis].setValue(val);
    }

    /// Resets where we currently point to.
    /// Returns true if we had to reset our cache line.  If we didn't,
    /// and you have multiple probes acting in-step, you can just
    /// advanceX() the other probes
    template <typename S> 
    bool         setIndex(UT_VoxelArrayIterator<S> &vit)
                 { return setIndex(vit.x(), vit.y(), vit.z()); }

    bool         setIndex(int x, int y, int z)
                 { 
                    if (myLines[0].setIndex(x, y, z))
                    {
                       myLines[1].setIndex(x, y, z);
                       myLines[2].setIndex(x, y, z);
                       return true;
                    }
                    myLines[1].advanceX();
                    myLines[2].advanceX();
                    return false;
                 }

    void         advanceX()
                 { myLines[0].advanceX(); myLines[1].advanceX(); myLines[2].advanceX(); }

protected:
    UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>    myLines[3];
};

template <typename T>
class
UT_VoxelProbeCube
{
public:
                 UT_VoxelProbeCube();
    virtual     ~UT_VoxelProbeCube();

    void         setCubeArray(UT_VoxelArray<T> *vox);
    void         setPlusArray(UT_VoxelArray<T> *vox);

    /// Allows you to query +/-1 in each direction.  In cube update,
    /// all are valid.  In plus update, only one of x y and z may be
    /// non zero.
    inline T     getValue(int x, int y, int z) const
                 { return myLines[y+1][z+1].getValue(x); }

    template <typename S> 
    bool         setIndexCube(UT_VoxelArrayIterator<S> &vit)
                 { return setIndexCube(vit.x(), vit.y(), vit.z()); }
    bool         setIndexCube(int x, int y, int z);

    template <typename S> 
    bool         setIndexPlus(UT_VoxelArrayIterator<S> &vit)
                 { return setIndexPlus(vit.x(), vit.y(), vit.z()); }
    bool         setIndexPlus(int x, int y, int z);

    /// Computes central difference gradient, does not scale
    /// by the step size (which is twice voxelsize)
    UT_Vector3   gradient() const
                 { return UT_Vector3(getValue(1,0,0) - getValue(-1,0,0),
                                     getValue(0,1,0) - getValue(0,-1,0),
                                     getValue(0,0,1) - getValue(0,0,-1)); }

protected:
    /// Does an rotation of our cache lines, ym becomes y0 and y0 becomes yp,
    /// so further queries with y+1 will be cache hits for 2 out of 3.
    static void  rotateLines(UT_VoxelProbe<T, true, false, false> &ym, UT_VoxelProbe<T, true, false, false> &y0, UT_VoxelProbe<T, true, false, false> &yp);

    UT_VoxelProbe<T, true, false, false>        myLines[3][3];
    /// Cached look up position.  myValid stores if they are
    /// valid values or not
    bool                        myValid;
    int                         myX, myY, myZ;
    /// Half inclusive [,) range of valid x queries for current cache.
    int                         myMinValidX, myMaxValidX;
};

///
/// UT_VoxelProbeAverage
/// 
/// When working with MAC grids one often has slightly misalgined
/// fields.  Ie, one field is at the half-grid spacing of another field.
/// The step values are 0 if the dimension is algined, -1 for half a step
/// back (ie, (val(-1)+val(0))/2) and 1 for half a step forward
/// (ie, (val(0)+val(1))/2)
///
template <typename T, int XStep, int YStep, int ZStep>
class
UT_VoxelProbeAverage
{
public:
                 UT_VoxelProbeAverage() {}
    virtual     ~UT_VoxelProbeAverage() {}

    void         setArray(UT_VoxelArray<T> *vox);

    template <typename S> 
    bool         setIndex(UT_VoxelArrayIterator<S> &vit)
                 { return setIndex(vit.x(), vit.y(), vit.z()); }
    bool         setIndex(int x, int y, int z);

    /// Returns the velocity centered at this index, thus an average
    /// of the values in each of our internal probes.
    inline T     getValue() const
    {
        if (ZStep)
            return (valueZ(1) + valueZ(0)) * 0.5;
        return valueZ(0);
    }

protected:
    inline T     valueZ(int z) const
    {
        if (YStep)
            return (valueYZ(1, z) + valueYZ(0, z)) * 0.5;
        return valueYZ(0, z);
    }

    inline T     valueYZ(int y, int z) const
    {
        if (XStep > 0)
            return (myLines[y][z].getValue(1) + myLines[y][z].getValue(0)) * 0.5;
        if (XStep < 0)
            return (myLines[y][z].getValue(-1) + myLines[y][z].getValue(0)) * 0.5;
        return myLines[y][z].getValue();
    }

    // Stores [Y][Z] lines.
    UT_VoxelProbe<T, true, false, false>        myLines[2][2];
};


///
/// UT_VoxelProbeFace is designed to walk over three velocity
/// fields that store face-centered values.  The indices refer
/// to the centers of the voxels.
///
template <typename T>
class
UT_VoxelProbeFace
{
public:
                 UT_VoxelProbeFace();
    virtual     ~UT_VoxelProbeFace();

    void         setArray(UT_VoxelArray<T> *vx, UT_VoxelArray<T> *vy, UT_VoxelArray<T> *vz);
    void         setVoxelSize(const UT_Vector3 &voxelsize);

    template <typename S> 
    bool         setIndex(UT_VoxelArrayIterator<S> &vit)
                 { return setIndex(vit.x(), vit.y(), vit.z()); }
    bool         setIndex(int x, int y, int z);

    /// Get the face values on each face component.
    /// Parameters are axis then side.
    /// 0 is the lower face, 1 the higher face.
    inline T    face(int axis, int side) const
    {
        if (axis == 0)
            return myLines[0][0].getValue(side);
        else
            return myLines[axis][side].getValue();
    }

    /// Returns the velocity centered at this index, thus an average
    /// of the values in each of our internal probes.
    inline UT_Vector3    value() const
    {
        return UT_Vector3(0.5f * (face(0, 0) + face(0, 1)),
                          0.5f * (face(1, 0) + face(1, 1)),
                          0.5f * (face(2, 0) + face(2, 1)));
    }

    /// Returns the divergence of this cell.
    inline T            divergence() const
    {
        return    (face(0,1)-face(0,0)) * myVoxelSize.x()
                + (face(1,1)-face(1,0)) * myVoxelSize.y()
                + (face(2,1)-face(2,0)) * myVoxelSize.z();

    }

protected:
    
    static void swapLines(UT_VoxelProbe<T, true, false, false> &ym, 
                            UT_VoxelProbe<T, true, false, false> &yp);
    
    
    UT_VoxelProbe<T, true, false, false>        myLines[3][2];

    /// Cached look up position.  myValid stores if they are
    /// valid values or not
    bool                        myValid;
    int                         myX, myY, myZ;
    /// Half inclusive [,) range of valid x queries for current cache.
    int                         myMinValidX, myMaxValidX;

    UT_Vector3                  myVoxelSize, myInvVoxelSize;
};


#if defined( WIN32 ) || defined( LINUX ) || defined( MBSD ) || defined(GAMEOS)
    #include "UT_VoxelArray.C"
#endif


// Typedefs for common voxel array types
typedef UT_VoxelArray<fpreal32>                 UT_VoxelArrayF;
typedef UT_VoxelMipMap<fpreal32>                UT_VoxelMipMapF;
typedef UT_VoxelArrayIterator<fpreal32>         UT_VoxelArrayIteratorF;
// Read only probe
typedef UT_VoxelProbe<fpreal32, true, false, false>     UT_VoxelProbeF;
typedef UT_VoxelVectorProbe<fpreal32, true, false, false> UT_VoxelVectorProbeF;
// Write only
typedef UT_VoxelProbe<fpreal32, false, true, false>     UT_VoxelWOProbeF;
typedef UT_VoxelVectorProbe<fpreal32, false, true, false> UT_VoxelVectorWOProbeF;
// Read/Write always writeback.
typedef UT_VoxelProbe<fpreal32, true, true, false>      UT_VoxelRWProbeF;
typedef UT_VoxelVectorProbe<fpreal32, true, true, false> UT_VoxelVectorRWProbeF;
// Read/Write with testing
typedef UT_VoxelProbe<fpreal32, true, true, true>       UT_VoxelRWTProbeF;
typedef UT_VoxelVectorProbe<fpreal32, true, true, true> UT_VoxelVectorRWTProbeF;


typedef UT_VoxelProbeCube<fpreal32>             UT_VoxelProbeCubeF;

#if 0
typedef UT_VoxelArrayHandle<fpreal32>           UT_VoxelArrayHandleF;
typedef UT_VoxelArrayHandleAutoReadLock<fpreal32>       UT_VoxelArrayHandleAutoReadLockF;
typedef UT_VoxelArrayHandleAutoWriteLock<fpreal32>      UT_VoxelArrayHandleAutoWriteLockF;
#endif

typedef UT_COWHandle<UT_VoxelArray<fpreal32> >  UT_VoxelArrayHandleF;
typedef UT_COWReadHandle<UT_VoxelArray<fpreal32> >      UT_VoxelArrayReadHandleF;
typedef UT_COWWriteHandle<UT_VoxelArray<fpreal32> >     UT_VoxelArrayWriteHandleF;

#endif

---