UT/UT_VoxelArray.C

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.C ( UT Library, C++)
 *
 * COMMENTS:
 *      Tiled Voxel Array Implementation.
 */

#include "UT_VoxelArray.h"
#include <SYS/SYS_Align.h>
#include "UT_BoundingBox.h"
#include "UT_Filter.h"
#include "UT_StackAlloc.h"

///
/// fpreal16 conversion functions
///
inline fpreal16
UTvoxelConvertFP16(fpreal16 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(fpreal32 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(fpreal64 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(int8 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(int16 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(int32 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(int64 a) { return a; }
inline fpreal16
UTvoxelConvertFP16(UT_Vector2 a) { return a.x(); }
inline fpreal16
UTvoxelConvertFP16(UT_Vector3 a) { return a.x(); }
inline fpreal16
UTvoxelConvertFP16(UT_Vector4 a) { return a.x(); }

///
/// VoxelTileCompress definitions
///
template <typename T>
void
UT_VoxelTileCompress<T>::findMinMax(const UT_VoxelTile<T> &tile,
                    T &min, T &max) const
{
    int         x, y, z;

    min = max = getValue(tile, 0, 0, 0);

    // Generic approach is to just use the getValue() ability.
    for (z = 0; z < tile.zres(); z++)
    {
        for (y = 0; y < tile.yres(); y++)
        {
            for (x = 0; x < tile.xres(); x++)
            {
                tile.expandMinMax(getValue(tile, x, y, z), min, max);
            }
        }
    }
    return;
}

//
// VoxelTile definitions.
//

template <typename T>
UT_VoxelTile<T>::UT_VoxelTile(int xres, int yres, int zres)
{
    myRes[0] = xres;
    myRes[1] = yres;
    myRes[2] = zres;

    myCompressionType = COMPRESS_CONSTANT;
    myData = SYSamalloc(sizeof(T));
    
    // It is not accidental that this is not a static_cast!
    // There isn't a UT_Vector3(fpreal) constructor (as we have a foolish
    // UT_Vector3(fpreal *) constructor that would cause massive problems)
    // but there is a UT_Vector3 operator=(fpreal)!
    ((T *)myData)[0] = 0;
}

template <typename T>
UT_VoxelTile<T>::~UT_VoxelTile()
{
    SYSafree(myData);
}

template <typename T>
UT_VoxelTile<T>::UT_VoxelTile(const UT_VoxelTile<T> &src)
{
    myData = 0;

    // Use assignment operator.
    *this = src;
}

template <typename T>
const UT_VoxelTile<T> &
UT_VoxelTile<T>::operator=(const UT_VoxelTile<T> &src)
{
    if (&src == this)
        return *this;
    
    if (myData)
        SYSafree(myData);
    myData = 0;
    
    myRes[0] = src.myRes[0];
    myRes[1] = src.myRes[1];
    myRes[2] = src.myRes[2];

    myCompressionType = src.myCompressionType;
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            myData = SYSamalloc(
                    sizeof(T) * myRes[0] * myRes[1] * myRes[2], 128);
            memcpy(myData, src.myData, sizeof(T) * myRes[0] * myRes[1] * myRes[2]);
            break;
        case COMPRESS_CONSTANT:
            myData = SYSamalloc(sizeof(T));
            memcpy(myData, src.myData, sizeof(T));
            break;
        case COMPRESS_RAWFULL:
            myData = SYSamalloc(
                    sizeof(T) * TILESIZE * TILESIZE * TILESIZE, 128);
            memcpy(myData, src.myData,
                    sizeof(T) * TILESIZE * TILESIZE * TILESIZE);
            break;
        case COMPRESS_FPREAL16:
            myData = SYSamalloc(
                    sizeof(fpreal16) * myRes[0] * myRes[1] * myRes[2], 128);
            memcpy(myData, src.myData,
                    sizeof(fpreal16) * myRes[0] * myRes[1] * myRes[2]);
            break;
        default:
        {
            UT_VoxelTileCompress<T>     *engine;

            engine = getCompressionEngine(myCompressionType);
            myData = SYSamalloc(engine->getDataLength(src));
            memcpy(myData, src.myData, engine->getDataLength(src));
            break;
        }
    }

    return *this;
}

template <typename T>
bool
UT_VoxelTile<T>::writeThrough(int x, int y, int z, T t)
{
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            // Do the assignment.
            ((T *)myData)[ ((z * myRes[1]) + y) * myRes[0] + x ] = t;
            return true;

        case COMPRESS_CONSTANT:
            if (*(T *)myData == t)
            {
                // Trivially true.
                return true;
            }
            return false;

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

        case COMPRESS_FPREAL16:
            return false;
    }

    // Use the compression engine.
    UT_VoxelTileCompress<T>     *engine;

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

template <typename T>
T
UT_VoxelTile<T>::lerp(int x, int y, int z, fpreal32 fx, fpreal32 fy, fpreal32 fz) const
{
    T           vx, vx1, vy, vy1, vz;

    switch (myCompressionType)
    {
        case COMPRESS_RAW:
        case COMPRESS_RAWFULL:
        {
            T           *data = (T *) myData;
            int          offset = (z * myRes[1] + y) * myRes[0] + x;
            int          yinc = myRes[0];
            int          zinc = myRes[0] * myRes[1];

            // Lerp x:x+1, y, z
            vx = lerpValues(data[offset], data[offset+1], fx);
            // Lerp x:x+1, y+1, z
            vx1 = lerpValues(data[offset+yinc], data[offset+yinc+1], fx);
            
            // Lerp x:x+1, y:y+1, z
            vy = lerpValues(vx, vx1, fy);
            
            // Lerp x:x+1, y, z+1
            vx = lerpValues(data[offset+zinc], data[offset+zinc+1], fx);
            // Lerp x:x+1, y+1, z+1
            vx1 = lerpValues(data[offset+zinc+yinc], data[offset+zinc+yinc+1], fx);

            // Lerp x:x+1, y:y+1, z+1
            vy1 = lerpValues(vx, vx1, fy);

            // Lerp x:x+1, y:y+1, z:z+1
            vz = lerpValues(vy, vy1, fz);
            break;
        }

        case COMPRESS_FPREAL16:
        {
            fpreal16    *data = (fpreal16 *) myData;
            int          offset = (z * myRes[1] + y) * myRes[0] + x;
            int          yinc = myRes[0];
            int          zinc = myRes[0] * myRes[1];
            fpreal16     vx, vx1, vy, vy1, vz;

            // Lerp x:x+1, y, z
            vx = SYSlerp(data[offset], data[offset+1], fx);
            // Lerp x:x+1, y+1, z
            vx1 = SYSlerp(data[offset+yinc], data[offset+yinc+1], fx);
            
            // Lerp x:x+1, y:y+1, z
            vy = SYSlerp(vx, vx1, fy);
            
            // Lerp x:x+1, y, z+1
            vx = SYSlerp(data[offset+zinc], data[offset+zinc+1], fx);
            // Lerp x:x+1, y+1, z+1
            vx1 = SYSlerp(data[offset+zinc+yinc], data[offset+zinc+yinc+1], fx);

            // Lerp x:x+1, y:y+1, z+1
            vy1 = SYSlerp(vx, vx1, fy);

            // Lerp x:x+1, y:y+1, z:z+1
            vz = lerpValues(vy, vy1, fz);

            return T(vz);
        }
        case COMPRESS_CONSTANT:
        {
            // This is quite trivial to do a trilerp on.
            vz = *(T *) myData;
            break;
        }
            
        default:
        {
            UT_VoxelTileCompress<T>     *engine;

            engine = getCompressionEngine(myCompressionType);
            // Lerp x:x+1, y, z
            vx  = lerpValues(engine->getValue(*this, x,   y, z),
                             engine->getValue(*this, x+1, y, z),
                             fx);
            // Lerp x:x+1, y+1, z
            vx1 = lerpValues(engine->getValue(*this, x,   y+1, z),
                             engine->getValue(*this, x+1, y+1, z),
                             fx);
            
            // Lerp x:x+1, y:y+1, z
            vy = lerpValues(vx, vx1, fy);
            
            // Lerp x:x+1, y, z+1
            vx  = lerpValues(engine->getValue(*this, x,   y, z+1),
                             engine->getValue(*this, x+1, y, z+1),
                             fx);
            // Lerp x:x+1, y+1, z+1
            vx1 = lerpValues(engine->getValue(*this, x,   y+1, z+1),
                             engine->getValue(*this, x+1, y+1, z+1),
                             fx);

            // Lerp x:x+1, y:y+1, z+1
            vy1 = lerpValues(vx, vx1, fy);

            // Lerp x:x+1, y:y+1, z:z+1
            vz = lerpValues(vy, vy1, fz);
            break;
        }
    }

    return vz;
}

#if 0
template <typename T>
T
UT_VoxelTile<T>::lerp(v4uf frac, int x, int y, int z) const
{
    v4uf                a, b;

    switch (myCompressionType)
    {
        case COMPRESS_RAW:
        case COMPRESS_RAWFULL:
        {
            T           *data = (T *) myData;
            int          offset = (z * myRes[1] + y) * myRes[0] + x;
            int          yinc = myRes[0];
            int          zinc = myRes[0] * myRes[1];

            a = v4uf( data[offset],
                      data[offset+zinc],
                      data[offset+yinc],
                      data[offset+yinc+zinc] );
            b = v4uf( data[offset+1],
                      data[offset+zinc+1],
                      data[offset+yinc+1],
                      data[offset+yinc+zinc+1] );
            break;
        }

        case COMPRESS_CONSTANT:
        {
            // This is quite trivial to do a trilerp on.
            return *(T *) myData;
        }
            
        default:
        {
            UT_VoxelTileCompress<T>     *engine;

            engine = getCompressionEngine(myCompressionType);
            // Lerp x:x+1, y, z
            a = v4uf( engine->getValue(*this, x, y, z),
                      engine->getValue(*this, x, y, z+1),
                      engine->getValue(*this, x, y+1, z),
                      engine->getValue(*this, x, y+1, z+1) );
            b = v4uf( engine->getValue(*this, x+1, y, z),
                      engine->getValue(*this, x+1, y, z+1),
                      engine->getValue(*this, x+1, y+1, z),
                      engine->getValue(*this, x+1, y+1, z+1) );
            break;
        }
    }

    v4uf                fx, fy, fz;

    fx = frac.swizzle<0, 0, 0, 0>();
    fy = frac.swizzle<1, 1, 1, 1>();
    fz = frac.swizzle<2, 2, 2, 2>();

    b -= a;
    a = madd(b, fx, a);

    b = a.swizzle<2, 3, 0, 1>();
    b -= a;
    a = madd(b, fy, a);

    b = a.swizzle<1, 2, 3, 0>();
    b -= a;
    a = madd(b, fz, a);

    return a[0];
}
#endif

template <typename T>
T *
UT_VoxelTile<T>::fillCacheLine(T *cacheline, int &stride, int x, int y, int z, bool forcecopy, bool strideofone) const
{
    UT_ASSERT_P(x >= 0 && y >= 0 && z >= 0);
    UT_ASSERT_P(x < myRes[0] && y < myRes[1] && z < myRes[2]);

    T           *src;
    int          i, xres = myRes[0];

    // All the directly handled types exit directly from this switch.
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            stride = 1;
            src = (T *)myData + (z * myRes[1] + y) * xres;
            if (!forcecopy)
                return &src[x];

            for (i = 0; i < xres; i++)
                cacheline[i] = src[i];

            return &cacheline[x];

        case COMPRESS_FPREAL16:
        {
            fpreal16            *src;

            stride = 1;
            src = (fpreal16 *)myData + (z * myRes[1] + y) * xres;

            for (i = 0; i < xres; i++)
                cacheline[i] = src[i];

            return &cacheline[x];
        }


        case COMPRESS_CONSTANT:
            src = (T *) myData;
            if (!forcecopy && !strideofone)
            {
                stride = 0;
                return src;
            }
            stride = 1;
            // Fill out a constant value
            for (i = 0; i < xres; i++)
                cacheline[i] = *src;

            return &cacheline[x];


        case COMPRESS_RAWFULL:
            stride = 1;
            src = (T *)myData + (z * TILESIZE + y) * TILESIZE;
            if (!forcecopy)
                return &src[x];

            for (i = 0; i < xres; i++)
                cacheline[i] = src[i];

            return &cacheline[x];
    }

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

    engine = getCompressionEngine(myCompressionType);

    // We could add support for a direct cacheline fill to accelerate
    // this case as well.
    stride = 1;
    for (i = 0; i < xres; i++)
        cacheline[i] = engine->getValue(*this, i, y, z);

    return &cacheline[x];
}

template <typename T>
void
UT_VoxelTile<T>::writeCacheLine(T *cacheline, int y, int z) 
{
    UT_ASSERT_P(y >= 0 && z >= 0);
    UT_ASSERT_P(y < myRes[1] && z < myRes[2]);

    T           *dst, value;
    int          i, xres = myRes[0];

    // All the directly handled types exit directly from this switch.
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            dst = (T *)myData + (z * myRes[1] + y) * xres;
            for (i = 0; i < xres; i++)
                *dst++ = *cacheline++;
            return;

        case COMPRESS_CONSTANT:
            value = *(T *) myData;
            for (i = 0; i < xres; i++)
                if (cacheline[i] != value)
                    break;
            // If everything was equal, our write is trivial.
            if (i == xres)
                return;
            
            break;

        case COMPRESS_RAWFULL:
            dst = (T *)myData + (z * TILESIZE + y) * TILESIZE;
            for (i = 0; i < TILESIZE; i++)
                *dst++ = *cacheline++;

            return;
    }

    // Switch back to invoking writeThrough.  Ideally we'd have
    // a version that can handle a whole cache line at once
    for (i = 0; i < xres; i++)
        if (!writeThrough(i, y, z, cacheline[i]))
            break;

    // Determine if we failed to write everything through
    if (i != xres)
    {
        // Uncompress and reinvoke ourselves.
        uncompress();
        writeCacheLine(cacheline, y, z);
    }
}

template <typename T>
void
UT_VoxelTile<T>::setValue(int x, int y, int z, T t)
{
    UT_ASSERT_P(x >= 0 && y >= 0 && z >= 0);
    UT_ASSERT_P(x < myRes[0] && y < myRes[1] && z < myRes[2]);

    // Determine if assignment is compatible with current
    // compression technique.
    if (writeThrough(x, y, z, t))
    {
        return;
    }
    // Can't write to our current type of tile with the
    // given value, so abandon.
    uncompress();

    // Attempt to write through again.  This must succeed
    // as we should now be uncompressed.
    bool        success;
    
    success = writeThrough(x, y, z, t);
    UT_ASSERT_P(success);
}

template <typename T>
void
UT_VoxelTile<T>::uncompress()
{
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            // Trivial.
            return;

        case COMPRESS_CONSTANT:
        {
            // Must expand the tile!
            T           cval;
            int         i, n;

            myCompressionType = COMPRESS_RAW;

            if (myRes[0] == TILESIZE &&
                myRes[1] == TILESIZE &&
                myRes[2] == TILESIZE)
            {
                // Flag that we can do fast lookups on this tile.
                myCompressionType = COMPRESS_RAWFULL;
            }

            cval = *(T *)myData;
            SYSafree(myData);
            n = myRes[0] * myRes[1] * myRes[2];
            myData = SYSamalloc(sizeof(T) * n, 128);

            for (i = 0; i < n; i++)
            {
                ((T *)myData)[i] = cval;
            }
            return;
        }
        case COMPRESS_RAWFULL:
        {
            T           *raw;
            int          x, y, z, i, n;

            if (myRes[0] == TILESIZE &&
                myRes[1] == TILESIZE &&
                myRes[2] == TILESIZE)
            {
                // Trivial
                return;
            }

            // Need to contract to the actual compact size.
            myCompressionType = COMPRESS_RAW;

            n = myRes[0] * myRes[1] * myRes[2];
            raw = (T *)SYSamalloc(sizeof(T) * n, 128);
            i = 0;
            for (z = 0; z < myRes[2]; z++)
            {
                for (y = 0; y < myRes[1]; y++)
                {
                    for (x = 0; x < myRes[0]; x++)
                    {
                        raw[i++] = ((T *)myData)[x+(y+z*TILESIZE)*TILESIZE];
                    }
                }
            }
            SYSafree(myData);
            myData = raw;

            return;
        }

        case COMPRESS_FPREAL16:
        {
            T                   *raw;
            int                 x, y, z, i, n;
            fpreal16            *src = (fpreal16 *) myData;

            myCompressionType = COMPRESS_RAW;

            n = myRes[0] * myRes[1] * myRes[2];
            raw = (T *)SYSamalloc(sizeof(T) * n, 128);
            i = 0;
            for (z = 0; z < myRes[2]; z++)
            {
                for (y = 0; y < myRes[1]; y++)
                {
                    for (x = 0; x < myRes[0]; x++)
                    {
                        raw[i] = src[i];
                        i++;
                    }
                }
            }
            SYSafree(myData);
            myData = raw;
            return;
        }
    }
    
    // Use the compression engine.
    UT_VoxelTileCompress<T>     *engine;

    engine = getCompressionEngine(myCompressionType);

    // We use the read ability to set our values.
    int                  x, y, z, i;
    T                   *raw;

    raw = (T *) SYSamalloc(sizeof(T) * myRes[0] * myRes[1] * myRes[2], 128);
    i = 0;
    for (z = 0; z < myRes[2]; z++)
    {
        for (y = 0; y < myRes[1]; y++)
        {
            for (x = 0; x < myRes[0]; x++)
            {
                raw[i++] = engine->getValue(*this, x, y, z);
            }
        }
    }

    // Now copy in the result
    myCompressionType = COMPRESS_RAW;
    if (myRes[0] == TILESIZE &&
        myRes[1] == TILESIZE &&
        myRes[2] == TILESIZE)
    {
        // Flag that we can do fast lookups on this tile.
        myCompressionType = COMPRESS_RAWFULL;
    }

    SYSafree(myData);
    myData = raw;
}

template <typename T>
void
UT_VoxelTile<T>::uncompressFull()
{
    T           *raw;
    int          x, y, z, i;

    if (myCompressionType == COMPRESS_RAWFULL)
        return;

    uncompress();

    UT_ASSERT(myCompressionType == COMPRESS_RAW);

    myCompressionType = COMPRESS_RAWFULL;

    // The RAWFULL format only differs from RAW when the tile dimensions
    // are smaller than the maximum tile size.
    if (myRes[0] < TILESIZE || myRes[1] < TILESIZE || myRes[2] < TILESIZE)
    {
        raw = (T *)SYSamalloc(sizeof(T) * TILESIZE * TILESIZE * TILESIZE, 128);
        i = 0;
        for (z = 0; z < myRes[2]; z++)
        {
            for (y = 0; y < myRes[1]; y++)
            {
                for (x = 0; x < myRes[0]; x++)
                {
                    raw[x+(y+z*TILESIZE)*TILESIZE] = ((T *)myData)[i++];
                }
            }
        }
        SYSafree(myData);
        myData = raw;
    }
}

template <typename T>
void
UT_VoxelTile<T>::findAverage(T &avg) const
{
    float               irx, iry, irz;

    irx = 1.0 / myRes[0];
    iry = 1.0 / myRes[1];
    irz = 1.0 / myRes[2];
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
        {
            int         i;
            const T     *data = (const T*) myData;

            i = 0;
            T           zavg;
            zavg = 0;
            for (int z = 0; z < myRes[2]; z++)
            {
                T               yavg;
                yavg = 0;
                for (int y = 0; y < myRes[1]; y++)
                {
                    T           xavg;
                    xavg = 0;
                    for (int x = 0; x < myRes[0]; x++)
                    {
                        xavg += data[i++];
                    }
                    xavg *= irx;
                    yavg += xavg;
                }
                yavg *= iry;
                zavg += yavg;
            }
            zavg *= irz;

            avg = zavg;
            return;
        }
            
        case COMPRESS_FPREAL16:
        {
            int         i;
            const fpreal16 *data = (fpreal16 *) myData;

            i = 0;
            fpreal16            zavg = 0;
            for (int z = 0; z < myRes[2]; z++)
            {
                fpreal16                yavg = 0;
                for (int y = 0; y < myRes[1]; y++)
                {
                    fpreal16            xavg = 0;
                    for (int x = 0; x < myRes[0]; x++)
                    {
                        xavg += data[i++];
                    }
                    xavg *= irx;
                    yavg += xavg;
                }
                yavg *= iry;
                zavg += yavg;
            }
            zavg *= irz;

            avg = zavg;
            return;
        }
            
        case COMPRESS_CONSTANT:
            avg = *(T*)myData;
            return;

        case COMPRESS_RAWFULL:
        {
            int         offset;
            const T     *data = (const T*) myData;

            offset = 0;
            T           zavg;
            zavg = 0;
            for (int z = 0; z < myRes[2]; z++)
            {
                T               yavg;
                yavg = 0;
                for (int y = 0; y < myRes[1]; y++)
                {
                    T           xavg;
                    xavg = 0;
                    for (int x = 0; x < myRes[0]; x++)
                    {
                        xavg += data[x+offset];
                    }
                    xavg *= irx;
                    yavg += xavg;
                    offset += TILESIZE;
                }
                yavg *= iry;
                zavg += yavg;
            }
            zavg *= irz;

            avg = zavg;
            return;
        }
            
        default:
        {
            T           zavg;
            zavg = 0;
            for (int z = 0; z < myRes[2]; z++)
            {
                T               yavg;
                yavg = 0;
                for (int y = 0; y < myRes[1]; y++)
                {
                    T           xavg;
                    xavg = 0;
                    for (int x = 0; x < myRes[0]; x++)
                    {
                        xavg += (*this)(x, y, z);
                    }
                    xavg *= irx;
                    yavg += xavg;
                }
                yavg *= iry;
                zavg += yavg;
            }
            zavg *= irz;

            avg = zavg;
            return;
        }
    }
}

template <typename T>
void
UT_VoxelTile<T>::findMinMax(T &min, T &max) const
{
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
        {
            int         n = myRes[0] * myRes[1] * myRes[2];
            int         i;

            min = max = *(T*)myData;
            for (i = 1; i < n; i++)
            {
                expandMinMax( ((T*)myData)[i], min, max );
            }
            return;
        }
            
        case COMPRESS_FPREAL16:
        {
            int         n = myRes[0] * myRes[1] * myRes[2];
            int         i;
            fpreal16    *src = (fpreal16 *)myData;

            min = max = *src;
            for (i = 1; i < n; i++)
            {
                T               val;
                val = src[i];
                expandMinMax( val, min, max );
            }
            return;
        }
            
        case COMPRESS_CONSTANT:
            min = max = *(T*)myData;
            return;

        case COMPRESS_RAWFULL:
        {
            int         x, y, z, offset;

            min = max = *(T*)myData;
            offset = 0;
            for (z = 0; z < myRes[2]; z++)
            {
                for (y = 0; y < myRes[1]; y++)
                {
                    for (x = 0; x < myRes[0]; x++)
                    {
                        expandMinMax(
                                ((T*)myData)[x+offset],
                                min, max );
                    }
                    offset += TILESIZE;
                }
            }
            return;
        }
            
        default:
        {
            // Use the compression engine.
            UT_VoxelTileCompress<T>     *engine;

            engine = getCompressionEngine(myCompressionType);

            engine->findMinMax(*this, min, max);
            return;
        }
    }
}

template <typename T>
bool
UT_VoxelTile<T>::hasNan() const
{
    switch (myCompressionType)
    {
        case COMPRESS_RAW:
        case COMPRESS_RAWFULL:
        {
            int         n = myRes[0] * myRes[1] * myRes[2];
            int         i;

            for (i = 0; i < n; i++)
            {
                if (SYSisNan( ((T*)myData)[i] ))
                    return true;
            }
            return false;
        }

        case COMPRESS_FPREAL16:
        {
            return false;
        }
            
        case COMPRESS_CONSTANT:
            if (SYSisNan(*(T*)myData))
                return true;
            return false;

        default:
        {
            // Use the compression engine.
            UT_VoxelTileCompress<T>     *engine;
            int         x, y, z;

            engine = getCompressionEngine(myCompressionType);

            for (z = 0; z < myRes[2]; z++)
                for (y = 0; y < myRes[1]; y++)
                    for (x = 0; x < myRes[0]; x++)
                        if (SYSisNan(engine->getValue(*this, x, y, z)))
                        {
                            return true;
                        }

            return false;
        }
    }
}

template <typename T>
bool
UT_VoxelTile<T>::tryCompress(const UT_VoxelCompressOptions &options)
{
    T           min, max;
    bool        losslessonly = (options.myQuantizeTol == 0.0);
    int         i;
    UT_VoxelTileCompress<T>     *engine;

    // This is as easy as it gets.
    if (myCompressionType == COMPRESS_CONSTANT)
        return false;

    findMinMax(min, max);

    // See if we can be made into a constant tile.
    if (dist(min, max) <= options.myConstantTol)
    {
        // Ignore if we are already constant.
        if (myCompressionType == COMPRESS_CONSTANT)
            return false;
        
        // We are fully constant!
        if (min != max)
        {
            T           avg;
            findAverage(avg);
            makeConstant(avg);
        }
        else
        {
            makeConstant(min);
        }
        return true;
    }

    bool                result = false;

    if (options.myAllowFP16)
    {
        if (myCompressionType == COMPRESS_RAW ||
            myCompressionType == COMPRESS_RAWFULL)
        {
            makeFpreal16();
            result = true;
        }
    }

    for (i = 0; i < getCompressionEngines().entries(); i++)
    {
        engine = getCompressionEngines()(i);

        // Ignore possibly expensive lossy compressions.
        if (losslessonly && !engine->isLossless())
            continue;

        if (engine->tryCompress(*this, options, min, max))
        {
            myCompressionType = i + COMPRESS_ENGINE;
            result = true;
            // We keep testing in case another compression engine
            // can get a better number...
        }
    }

    // No suitable compression found.
    return result;
}

template <typename T>
void
UT_VoxelTile<T>::makeConstant(T t)
{
    if (!isConstant())
    {
        SYSafree(myData);
        myData = SYSamalloc(sizeof(T));
    }

    *(T*)myData = t;
    myCompressionType = COMPRESS_CONSTANT;
}

template <typename T>
void
UT_VoxelTile<T>::makeFpreal16()
{
    if (isConstant())
    {
        return;
    }

    if (myCompressionType == COMPRESS_FPREAL16)
        return;

    // Get our new data.
    int                  len = myRes[2] * myRes[1] * myRes[0];
    fpreal16            *data = (fpreal16 *)SYSamalloc(sizeof(fpreal16) * len, 128);

    if (myCompressionType == COMPRESS_RAW || 
        myCompressionType == COMPRESS_RAWFULL)
    {
        for (int i = 0; i < len; i++)
        {
            data[i] = UTvoxelConvertFP16( ((T*)myData)[i] );
        }
    }
    else
    {
        // Apply any converters.
        int             i = 0;

        for (int z = 0; z < myRes[2]; z++)
        {
            for (int y = 0; y < myRes[1]; y++)
            {
                for (int x = 0; x < myRes[0]; x++)
                {
                    data[i++] = UTvoxelConvertFP16( (*this)(x, y, z));
                }
            }
        }
    }

    SYSafree(myData);
    myData = data;
    myCompressionType = COMPRESS_FPREAL16;
}

template <typename T>
int64
UT_VoxelTile<T>::getMemoryUsage() const
{
    int64               usage;

    usage = sizeof(*this);
    usage += getDataLength();
    return usage;
}

template <typename T>
int64
UT_VoxelTile<T>::getDataLength() const
{
    int64               usage;

    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            usage = sizeof(T) * xres() * yres() * zres();
            break;
        case COMPRESS_FPREAL16:
            usage = sizeof(fpreal16) * xres() * yres() * zres();
            break;
        case COMPRESS_CONSTANT:
            usage = sizeof(T);
            break;
        case COMPRESS_RAWFULL:
            usage = sizeof(T) * TILESIZE * TILESIZE * TILESIZE;
            break;
        default:
        {
            // Use the compression engine.
            UT_VoxelTileCompress<T>     *engine;
            engine = getCompressionEngine(myCompressionType);
            usage = engine->getDataLength(*this);
            break;
        }
    }
    return usage;
}

template <typename T>
void
UT_VoxelTile<T>::weightedSum(int pstart[3], int pend[3],
                             float *weights[3], int start[3],
                             T &result)
{
    int         ix, iy, iz, i;
    int         px, py, pz;
    int         ixstart, ixend;
    int         tstart[3];
    fpreal      w, pw;
    T           psumx, psumy;

    switch (myCompressionType)
    {
        case COMPRESS_CONSTANT:
        {
            w = 1;
            for (i = 0; i < 3; i++)
            {
                pw = 0;
                for (ix = 0; ix < pend[i]-pstart[i]; ix++)
                    pw += weights[i][ix+pstart[i]-start[i]];
                w *= pw;
            }
            result += w * ((T*)myData)[0];
            break;
        }

        case COMPRESS_RAW:
        {
            tstart[0] = pstart[0] & TILEMASK;
            tstart[1] = pstart[1] & TILEMASK;
            tstart[2] = pstart[2] & TILEMASK;
            ixstart     = pstart[0]-start[0];
            ixend       = pend[0]-start[0];
            pz = tstart[2];
            UT_ASSERT(pz < myRes[2]);
            UT_ASSERT(ixend - ixstart <= myRes[0]);
            for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
            {
                psumy = 0;
                py = tstart[1];
                UT_ASSERT(py < myRes[1]);
                for (iy = pstart[1]; iy < pend[1]; iy++, py++)
                {
                    psumx = 0;
                    px = ((pz * myRes[1]) + py) * myRes[0] + tstart[0];
                    for (ix = ixstart; ix < ixend; ix++, px++)
                    {
                        psumx += weights[0][ix]* ((T*)myData)[px];
                    }
                    psumy += weights[1][iy-start[1]] * psumx;
                }
                result += weights[2][iz-start[2]] * psumy;
            }
            break;
        }

        case COMPRESS_FPREAL16:
        {
            fpreal16            *src = (fpreal16 *) myData;

            tstart[0] = pstart[0] & TILEMASK;
            tstart[1] = pstart[1] & TILEMASK;
            tstart[2] = pstart[2] & TILEMASK;
            ixstart     = pstart[0]-start[0];
            ixend       = pend[0]-start[0];
            pz = tstart[2];
            UT_ASSERT(pz < myRes[2]);
            UT_ASSERT(ixend - ixstart <= myRes[0]);
            for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
            {
                psumy = 0;
                py = tstart[1];
                UT_ASSERT(py < myRes[1]);
                for (iy = pstart[1]; iy < pend[1]; iy++, py++)
                {
                    psumx = 0;
                    px = ((pz * myRes[1]) + py) * myRes[0] + tstart[0];
                    for (ix = ixstart; ix < ixend; ix++, px++)
                    {
                        T               val;
                        val = src[px];
                        psumx += weights[0][ix]* val;
                    }
                    psumy += weights[1][iy-start[1]] * psumx;
                }
                result += weights[2][iz-start[2]] * psumy;
            }
            break;
        }
        case COMPRESS_RAWFULL:
        {
            tstart[0] = pstart[0] & TILEMASK;
            tstart[1] = pstart[1] & TILEMASK;
            tstart[2] = pstart[2] & TILEMASK;
            ixstart     = pstart[0]-start[0];
            ixend       = pend[0]-start[0];
            pz = tstart[2];
            for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
            {
                psumy = 0;
                py = tstart[1];
                for (iy = pstart[1]; iy < pend[1]; iy++, py++)
                {
                    psumx = 0;
                    px = ((pz * TILESIZE) + py) * TILESIZE + tstart[0];
                    for (ix = ixstart; ix < ixend; ix++, px++)
                    {
                        psumx += weights[0][ix]* ((T*)myData)[px];
                    }
                    psumy += weights[1][iy-start[1]] * psumx;
                }
                result += weights[2][iz-start[2]] * psumy;
            }
            break;
        }

        default:
        {
            // For all other compression types, we use our
            // getValue() accessor rather than trying to
            // do anything fancy.
            tstart[0] = pstart[0] & TILEMASK;
            tstart[1] = pstart[1] & TILEMASK;
            tstart[2] = pstart[2] & TILEMASK;
            ixstart     = pstart[0]-start[0];
            ixend       = pend[0]-start[0];
            pz = tstart[2];
            for (iz = pstart[2]; iz < pend[2]; iz++, pz++)
            {
                psumy = 0;
                py = tstart[1];
                for (iy = pstart[1]; iy < pend[1]; iy++, py++)
                {
                    psumx = 0;
                    px = tstart[0];
                    for (ix = ixstart; ix < ixend; ix++, px++)
                    {
                        psumx += weights[0][ix] * 
                                  (*this)(px, py, pz);
                    }
                    psumy += weights[1][iy-start[1]] * psumx;
                }
                result += weights[2][iz-start[2]] * psumy;
            }
            break;
        }
    }
}

template <typename T>
void
UT_VoxelTile<T>::registerCompressionEngine(UT_VoxelTileCompress<T> *engine)
{
    int                          i;

    // Repalce old copy, if any.
    for (i = 0; i < getCompressionEngines().entries(); i++)
    {
        if (!strcmp(engine->getName(), getCompressionEngines()(i)->getName()))
        {
            getCompressionEngines()(i) = engine;
            return;
        }
    }
            
    getCompressionEngines().append(engine);
}

template <typename T>
int
UT_VoxelTile<T>::lookupCompressionEngine(const char *name)
{
    int                          i;

    if (!name)
        return -1;

    for (i = 0; i < getCompressionEngines().entries(); i++)
    {
        if (!strcmp(name, getCompressionEngines()(i)->getName()))
        {
            return i + COMPRESS_ENGINE;
        }
    }

    return -1;
}

template <typename T>
UT_VoxelTileCompress<T> *
UT_VoxelTile<T>::getCompressionEngine(int index)
{
    index -= COMPRESS_ENGINE;

    return getCompressionEngines()(index);
}

template <typename T>
void
UT_VoxelTile<T>::save(ostream &os)
{
    // Always save in our native format...
    if (myCompressionType >= COMPRESS_ENGINE)
    {
        UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);

        if (engine->canSave())
        {
            char        type = myCompressionType;

            UTwrite(os, &type, 1);
            engine->save(os, *this);
            return;
        }

        // Can't save, must uncompress.
        uncompress();
    }

    int         len;
    char        type = myCompressionType;

    UT_ASSERT(type >= 0 && type < COMPRESS_ENGINE);

    UTwrite(os, &type, 1);

    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            len = myRes[2] * myRes[1] * myRes[0];
            UTwrite(os, (T *) myData, len);
            break;

        case COMPRESS_FPREAL16:
            len = myRes[2] * myRes[1] * myRes[0];
            UTwrite(os, (int16 *) myData, len);
            break;

        case COMPRESS_RAWFULL:
            len = TILESIZE * TILESIZE * TILESIZE;
            UTwrite(os, (T *) myData, len);
            break;

        case COMPRESS_CONSTANT:
            UTwrite(os, (T *) myData, 1);
            break;
    }
}

template <typename T>
void
UT_VoxelTile<T>::load(UT_IStream &is, const UT_IntArray &compress)
{
    char                type, otype;
    int                 len;

    is.readChar(type);

    // Perform the indirection to find out our native type.
    if (type >= 0 && type < compress.entries())
    {
        otype = type;
        type = compress(type);
        if (type == -1)
        {
            cerr << "Missing compression engine " << (int) otype << endl;
        }
    }

    if (type >= COMPRESS_ENGINE)
    {
        myCompressionType = type;

        UT_VoxelTileCompress<T> *engine = getCompressionEngine(myCompressionType);

        engine->load(is, *this);

        return;
    }

    UT_ASSERT(type >= 0);
    if (type < 0)
    {
        return;
    }

    if (myData)
        SYSafree(myData);

    myCompressionType = type;

    switch (myCompressionType)
    {
        case COMPRESS_RAW:
            len = myRes[2] * myRes[1] * myRes[0];
            myData = SYSamalloc(sizeof(T) * len, 128);
            is.read((T *) myData, len);
            break;

        case COMPRESS_FPREAL16:
            len = myRes[2] * myRes[1] * myRes[0];
            myData = SYSamalloc(sizeof(fpreal16) * len, 128);
            is.read((int16 *) myData, len);
            break;

        case COMPRESS_RAWFULL:
            len = TILESIZE * TILESIZE * TILESIZE;
            myData = SYSamalloc(sizeof(T) * len, 128);
            is.read((T *) myData, len);
            break;

        case COMPRESS_CONSTANT:
            myData = SYSamalloc(sizeof(T));
            is.read((T *) myData, 1);
            break;
    }
}

template <typename T>
void
UT_VoxelTile<T>::saveCompressionTypes(ostream &os)
{
    int16       ntype = getCompressionEngines().entries();
    int         i;

    ntype += COMPRESS_ENGINE;

    UTwrite(os, &ntype);

    UT_ASSERT(COMPRESS_ENGINE == 4);    // Change lower lines!
    UTsaveStringBinary(os, "raw", UT_STRING_8BIT_IO);
    UTsaveStringBinary(os, "rawfull", UT_STRING_8BIT_IO);
    UTsaveStringBinary(os, "constant", UT_STRING_8BIT_IO);
    UTsaveStringBinary(os, "fpreal16", UT_STRING_8BIT_IO);

    ntype -= COMPRESS_ENGINE;
    for (i = 0; i < ntype; i++)
    {
        UTsaveStringBinary(os, getCompressionEngines()(i)->getName(), UT_STRING_8BIT_IO);
    }
}

template <typename T>
void
UT_VoxelTile<T>::loadCompressionTypes(UT_IStream &is, UT_IntArray &compress)
{
    int16       ntype;
    int         i, idx;

    compress.entries(0);

    is.read(&ntype);

    for (i = 0; i < ntype; i++)
    {
        UT_String               name;

        is.readBinaryString(name, UT_ISTREAM_8BIT_IO);
        if (name == "raw")
            compress.append(COMPRESS_RAW);
        else if (name == "rawfull")
            compress.append(COMPRESS_RAWFULL);
        else if (name == "constant")
            compress.append(COMPRESS_CONSTANT);
        else if (name == "fpreal16")
            compress.append(COMPRESS_FPREAL16);
        else
        {
            idx = lookupCompressionEngine(name);

            // -1 means a failure to find it in our set..
            // this is only a bad thing if a tile actually uses this engine.

            compress.append(idx);
        }
    }
}



//
// VoxelArray definitions.
//

template <typename T>
UT_VoxelArray<T>::UT_VoxelArray()
{
    myRes[0] = 0;
    myRes[1] = 0;
    myRes[2] = 0;
    myTileRes[0] = 0;
    myTileRes[1] = 0;
    myTileRes[2] = 0;

    myTiles = 0;

    myBorderValue = 0;
    myBorderScale[0] = 0;
    myBorderScale[1] = 0;
    myBorderScale[2] = 0;
    myBorderType = UT_VOXELBORDER_STREAK;
}

template <typename T>
UT_VoxelArray<T>::~UT_VoxelArray()
{
    deleteVoxels();
}

template <typename T>
UT_VoxelArray<T>::UT_VoxelArray(const UT_VoxelArray<T> &src)
{
    myRes[0] = 0;
    myRes[1] = 0;
    myRes[2] = 0;
    myTileRes[0] = 0;
    myTileRes[1] = 0;
    myTileRes[2] = 0;

    myTiles = 0;

    *this = src;
}

template <typename T>
const UT_VoxelArray<T> &
UT_VoxelArray<T>::operator=(const UT_VoxelArray<T> &src)
{
    // Paranoid code:
    if (&src == this)
        return *this;

    // Clear out my own data (this is done again by size, but not
    // maybe with resize?)
    deleteVoxels();

    myBorderScale[0] = src.myBorderScale[0];
    myBorderScale[1] = src.myBorderScale[1];
    myBorderScale[2] = src.myBorderScale[2];
    myBorderValue = src.myBorderValue;
    myBorderType = src.myBorderType;

    myCompressionOptions = src.myCompressionOptions;

    // Allocate proper size.
    size(src.myRes[0], src.myRes[1], src.myRes[2]);

    int                 i, ntiles;

    ntiles = numTiles();
    
    for (i = 0; i < ntiles; i++)
    {
        *myTiles[i] = *src.myTiles[i];
    }

    return *this;
}

template <typename T>
void
UT_VoxelArray<T>::deleteVoxels()
{
    int         i, n;

    n = numTiles();
    for (i = 0; i < n; i++)
        delete myTiles[i];
    delete [] myTiles;
    myTiles = 0;

    myTileRes[0] = myTileRes[1] = myTileRes[2] = 0;
    myRes[0] = myRes[1] = myRes[2] = 0;
}

template <typename T>
void
UT_VoxelArray<T>::size(int xres, int yres, int zres)
{
    deleteVoxels();
    
    myRes[0] = xres;
    myRes[1] = yres;
    myRes[2] = zres;

    // Initialize our tiles.
    myTileRes[0] = (xres + TILEMASK) >> TILEBITS;
    myTileRes[1] = (yres + TILEMASK) >> TILEBITS;
    myTileRes[2] = (zres + TILEMASK) >> TILEBITS;

    int                 ntiles;

    ntiles = myTileRes[0] * myTileRes[1] * myTileRes[2];
    
    if (ntiles)
    {
        myTiles = new UT_VoxelTile<T> *[ntiles];

        int             tx, ty, tz, i;
        int             xr, yr, zr;     

        i = 0;
        for (tz = 0; tz < myTileRes[2]; tz++)
        {
            if (tz < myTileRes[2]-1)
                zr = TILESIZE;
            else
                zr = zres - tz * TILESIZE;

            for (ty = 0; ty < myTileRes[1]; ty++)
            {
                if (ty < myTileRes[1]-1)
                    yr = TILESIZE;
                else
                    yr = yres - ty * TILESIZE;
                
                xr = TILESIZE;
                for (tx = 0; tx < myTileRes[0]-1; tx++)
                {
                    myTiles[i] = new UT_VoxelTile<T>(xr, yr, zr);
                    
                    i++;
                }
                xr = xres - tx * TILESIZE;
                myTiles[i] = new UT_VoxelTile<T>(xr, yr, zr);
                i++;
            }
        }
    }
    else
        myTiles = 0;
}

template <typename T>
void
UT_VoxelArray<T>::match(const UT_VoxelArray<T> &src)
{
    // Only resize if our res is different.
    if (src.getXRes() != getXRes() ||
        src.getYRes() != getYRes() ||
        src.getZRes() != getZRes())
    {
        size(src.getXRes(), src.getYRes(), src.getZRes());
    }

    // Update border conditions.
    myBorderType = src.myBorderType;
    myBorderScale[0] = src.myBorderScale[0];
    myBorderScale[1] = src.myBorderScale[1];
    myBorderScale[2] = src.myBorderScale[2];
    myBorderValue = src.myBorderValue;

    // Match our compression tolerances
    myCompressionOptions = src.myCompressionOptions;
}

template <typename T>
int64
UT_VoxelArray<T>::getMemoryUsage() const
{
    int64        totaltilesize = 0;
    int          i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
        totaltilesize += myTiles[i]->getMemoryUsage();

    return sizeof(*this) + totaltilesize;
}

template <typename T>
void
UT_VoxelArray<T>::constant(T t)
{
    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        myTiles[i]->makeConstant(t);
    }
}

template <typename T>
bool
UT_VoxelArray<T>::isConstant(T *t) const
{
    int         i, ntiles;
    T           cval;
    const fpreal        tol = UT_FTOLERANCE;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        if (!myTiles[i]->isConstant())
        {
            return false;
        }

        if (!i)
        {
            // First tile, get the constant value.
            cval = (*myTiles[i])(0, 0, 0);
        }
        else
        {
            // See if we have deviated too much.
            if (UT_VoxelTile<T>::dist(cval, (*myTiles[i])(0, 0, 0)) > tol)
            {
                return false;
            }
        }
    }

    // All tiles are both constant and within tolerance of each
    // other.  Write out our constant value and return true.
    if (t)
        *t = cval;

    return true;
}

template <typename T>
bool
UT_VoxelArray<T>::hasNan() const
{
    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        if (myTiles[i]->hasNan())
        {
            return true;
        }
    }

    return false;
}

template <typename T>
T
UT_VoxelArray<T>::operator()(UT_Vector3F pos) const
{
#if 0
    v4uf                pos4(pos.x(), pos.y(), pos.z(), 0);


    return (*this)(pos4);
#else
    int                 x, y, z;
    // Yes, these have to be 32 becaues split float requires 32!
    fpreal32            fx, fy, fz;

    // We go from the position in the unit cube into the index.
    // The center of cells must map to the exact integer indices.
    pos.x() *= myRes[0];
    pos.y() *= myRes[1];
    pos.z() *= myRes[2];
    pos.x() -= 0.5;
    pos.y() -= 0.5;
    pos.z() -= 0.5;

    // Determine integer & fractional components.
    fx = pos.x();
    SYSfastSplitFloat(fx, x);
    fy = pos.y();
    SYSfastSplitFloat(fy, y);
    fz = pos.z();
    SYSfastSplitFloat(fz, z);

    // Do trilinear interpolation.
    T                   vx, vx1, vy, vy1, vz;

    // NOTE:
    // If you get a crash reading one of these getValues, note
    // that if you are reading from the same voxel array that you
    // are writing to, and doing so in a multi-threaded, tile-by-tile
    // approach, just because your positions match doesn't mean you
    // won't access neighbouring tiles.  To avoid this, perform
    // the obvious optimization of checking for aligned tiles
    // and using direct indices.
    // (I have avoided using if (!fx && !fy && !fz) as the test
    // as that is both potentially inaccurate (round off in conversion
    // may mess us up) and said optimzation is a useful one anyways)

    // Optimization for common cases (values are within the voxel range and
    // are all within the same tile)
    if ( (x > 0) && (y > 0) && (z > 0) &&
        // Do not use x+1 < foo as if x is MAX_INT that will falsely
        // report in bounds!
         (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
    {
        int                     xm, ym, zm;

        xm = x & TILEMASK;
        ym = y & TILEMASK;
        zm = z & TILEMASK;

        if ((xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
        {
            const UT_VoxelTile<T>       *tile =
                getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);

            vz = tile->lerp(xm, ym, zm, fx, fy, fz);
        }
        else
        {
            // We cross tile boundaries but we remain within
            // the voxel grid.  We can thus avoid any bounds checking
            // and use operator() rather than getValue.

            // Lerp x:x+1, y, z
            vx = UT_VoxelTile<T>::lerpValues((*this)(x,   y, z),
                                             (*this)(x+1, y, z),
                                             fx);
            // Lerp x:x+1, y+1, z
            vx1= UT_VoxelTile<T>::lerpValues((*this)(x,   y+1, z),
                                             (*this)(x+1, y+1, z),
                                             fx);
            // Lerp x:x+1, y:y+1, z
            vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
            
            // Lerp x:x+1, y, z+1
            vx = UT_VoxelTile<T>::lerpValues((*this)(x,   y, z+1),
                                             (*this)(x+1, y, z+1),
                                             fx);
            // Lerp x:x+1, y+1, z+1
            vx1= UT_VoxelTile<T>::lerpValues((*this)(x,   y+1, z+1),
                                             (*this)(x+1, y+1, z+1),
                                             fx);

            // Lerp x:x+1, y:y+1, z+1
            vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);

            // Lerp x:x+1, y:y+1, z:z+1
            vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
        }
    }
    else
    {
        // Lerp x:x+1, y, z
        vx = UT_VoxelTile<T>::lerpValues(getValue(x,   y, z),
                                         getValue(x+1, y, z),
                                         fx);
        // Lerp x:x+1, y+1, z
        vx1= UT_VoxelTile<T>::lerpValues(getValue(x,   y+1, z),
                                         getValue(x+1, y+1, z),
                                         fx);
        
        // Lerp x:x+1, y:y+1, z
        vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
        
        // Lerp x:x+1, y, z+1
        vx = UT_VoxelTile<T>::lerpValues(getValue(x,   y, z+1),
                                         getValue(x+1, y, z+1),
                                         fx);
        // Lerp x:x+1, y+1, z+1
        vx1= UT_VoxelTile<T>::lerpValues(getValue(x,   y+1, z+1),
                                         getValue(x+1, y+1, z+1),
                                         fx);

        // Lerp x:x+1, y:y+1, z+1
        vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);

        // Lerp x:x+1, y:y+1, z:z+1
        vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
    }

    return vz;
#endif
}

template <typename T>
T
UT_VoxelArray<T>::operator()(UT_Vector3D pos) const
{
    int                 x, y, z;
    // Yes, these have to be 32 becaues split float requires 32!
    fpreal32            fx, fy, fz;

    // We go from the position in the unit cube into the index.
    // The center of cells must map to the exact integer indices.
    pos.x() *= myRes[0];
    pos.y() *= myRes[1];
    pos.z() *= myRes[2];
    pos.x() -= 0.5;
    pos.y() -= 0.5;
    pos.z() -= 0.5;

    // Determine integer & fractional components.
    fx = pos.x();
    SYSfastSplitFloat(fx, x);
    fy = pos.y();
    SYSfastSplitFloat(fy, y);
    fz = pos.z();
    SYSfastSplitFloat(fz, z);

    // Do trilinear interpolation.
    T                   vx, vx1, vy, vy1, vz;

    // NOTE:
    // If you get a crash reading one of these getValues, note
    // that if you are reading from the same voxel array that you
    // are writing to, and doing so in a multi-threaded, tile-by-tile
    // approach, just because your positions match doesn't mean you
    // won't access neighbouring tiles.  To avoid this, perform
    // the obvious optimization of checking for aligned tiles
    // and using direct indices.
    // (I have avoided using if (!fx && !fy && !fz) as the test
    // as that is both potentially inaccurate (round off in conversion
    // may mess us up) and said optimzation is a useful one anyways)

    // Optimization for common cases (values are within the voxel range and
    // are all within the same tile)
    if ( (x > 0) && (y > 0) && (z > 0) &&
        // Do not use x+1 < foo as if x is MAX_INT that will falsely
        // report in bounds!
         (x < myRes[0]-1) && (y < myRes[1]-1) && (z < myRes[2]-1) )
    {
        int                     xm, ym, zm;

        xm = x & TILEMASK;
        ym = y & TILEMASK;
        zm = z & TILEMASK;

        if ((xm != TILEMASK) && (ym != TILEMASK) && (zm != TILEMASK))
        {
            const UT_VoxelTile<T>       *tile =
                getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);

            vz = tile->lerp(xm, ym, zm, fx, fy, fz);
        }
        else
        {
            // We cross tile boundaries but we remain within
            // the voxel grid.  We can thus avoid any bounds checking
            // and use operator() rather than getValue.

            // Lerp x:x+1, y, z
            vx = UT_VoxelTile<T>::lerpValues((*this)(x,   y, z),
                                             (*this)(x+1, y, z),
                                             fx);
            // Lerp x:x+1, y+1, z
            vx1= UT_VoxelTile<T>::lerpValues((*this)(x,   y+1, z),
                                             (*this)(x+1, y+1, z),
                                             fx);
            // Lerp x:x+1, y:y+1, z
            vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
            
            // Lerp x:x+1, y, z+1
            vx = UT_VoxelTile<T>::lerpValues((*this)(x,   y, z+1),
                                             (*this)(x+1, y, z+1),
                                             fx);
            // Lerp x:x+1, y+1, z+1
            vx1= UT_VoxelTile<T>::lerpValues((*this)(x,   y+1, z+1),
                                             (*this)(x+1, y+1, z+1),
                                             fx);

            // Lerp x:x+1, y:y+1, z+1
            vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);

            // Lerp x:x+1, y:y+1, z:z+1
            vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
        }
    }
    else
    {
        // Lerp x:x+1, y, z
        vx = UT_VoxelTile<T>::lerpValues(getValue(x,   y, z),
                                         getValue(x+1, y, z),
                                         fx);
        // Lerp x:x+1, y+1, z
        vx1= UT_VoxelTile<T>::lerpValues(getValue(x,   y+1, z),
                                         getValue(x+1, y+1, z),
                                         fx);
        
        // Lerp x:x+1, y:y+1, z
        vy = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);
        
        // Lerp x:x+1, y, z+1
        vx = UT_VoxelTile<T>::lerpValues(getValue(x,   y, z+1),
                                         getValue(x+1, y, z+1),
                                         fx);
        // Lerp x:x+1, y+1, z+1
        vx1= UT_VoxelTile<T>::lerpValues(getValue(x,   y+1, z+1),
                                         getValue(x+1, y+1, z+1),
                                         fx);

        // Lerp x:x+1, y:y+1, z+1
        vy1 = UT_VoxelTile<T>::lerpValues(vx, vx1, fy);

        // Lerp x:x+1, y:y+1, z:z+1
        vz = UT_VoxelTile<T>::lerpValues(vy, vy1, fz);
    }

    return vz;
}

#if 0
template <typename T>
T
UT_VoxelArray<T>::operator()(v4uf pos) const
{
    int                 x, y, z;

    v4ui                idx;

    // We go from the position in the unit cube into the index.
    // The center of cells must map to the exact integer indices.
    pos *= v4uf((float) myRes[0], (float) myRes[1], (float) myRes[2], 0.0f);

    // Because we use truncation we'll get inaccurate results
    // at the zero boundary, we thus bump up by two
    pos += 1.5;

    idx = pos.splitFloat();

    // And counteract that bump
    idx -= 2;

    x = idx[0];
    y = idx[1];
    z = idx[2];

    // Do trilinear interpolation.
    //      A      |        B
    // 0 +z +y +yz | +x +xz +xy +xyz

    v4uf                a, b, fx, fy, fz;

    int                 xm, ym, zm;

    xm = x & TILEMASK;
    ym = y & TILEMASK;
    zm = z & TILEMASK;

    // Optimization for common case (values are within the voxel range and
    // are all within the same tile)
    if ( (x > 0) & (y > 0) & (z > 0) &
         (x < myRes[0]-1) & (y < myRes[1]-1) & (z < myRes[2]-1) &
         (xm != TILEMASK) & (ym != TILEMASK) & (zm != TILEMASK) )
/*
    if (isValidIndex(x, y, z) && isValidIndex(x+1, y+1, z+1) &&
            (x & TILEMASK) != TILEMASK &&
            (y & TILEMASK) != TILEMASK &&
            (z & TILEMASK) != TILEMASK)
*/
    {
        const UT_VoxelTile<T>   *tile =
            getTile(x >> TILEBITS, y >> TILEBITS, z >> TILEBITS);

        return tile->lerp(pos, xm, ym, zm);
    }
    else
    {
        a = v4uf( getValue(x, y, z),
                  getValue(x, y, z+1),
                  getValue(x, y+1, z),
                  getValue(x, y+1, z+1) );
        b = v4uf( getValue(x+1, y, z),
                  getValue(x+1, y, z+1),
                  getValue(x+1, y+1, z),
                  getValue(x+1, y+1, z+1) );
    }

    fx = pos.swizzle<0, 0, 0, 0>();
    fy = pos.swizzle<1, 1, 1, 1>();
    fz = pos.swizzle<2, 2, 2, 2>();

    b -= a;
    a = madd(b, fx, a);

    b = a.swizzle<2, 3, 0, 1>();
    b -= a;
    a = madd(b, fy, a);

    b = a.swizzle<1, 2, 3, 0>();
    b -= a;
    a = madd(b, fz, a);

    return a[0];
}
#endif

static inline int
firstTile(int &start, int &end, int res)
{
    if (start < 0)
    {
        start += res;
        end += res;
    }
    return start;
}

static inline int
nextTile(int &tile, int &pstart, int &start, int &end, int res)
{
    int         pend;

    if (pstart >= res)
    {
        pstart -= res;
        start -= res;
        end -= res;
    }
    tile = pstart >> TILEBITS;
    pend = SYSmin((tile+1) * TILESIZE, end, res);

    UT_ASSERT(pstart >= 0 && pstart < res);
    UT_ASSERT(pend > 0 && pstart <= res);
    UT_ASSERT(pend-pstart > 0);
    UT_ASSERT(pend-pstart <= TILESIZE);

    return pend;
}

template <typename T>
T
UT_VoxelArray<T>::evaluate(const UT_Vector3 &pos, const UT_Filter &filter,
                           fpreal radius) const
{
    UT_Vector3           tpos;
    UT_FilterWindow      win[3];
    UT_FilterWrap        wrap;
    UT_VoxelTile<T>     *tile;
    float               *weights[3];
    float               *wdata;
    fpreal               iradius, visible;
    int                  start[3], end[3], size[3];
    int                  pstart[3], pend[3];
    int                  tx, ty, tz, i;
    T                    result;

    if (!myTiles)
    {
        // If the array is empty, just use the border value.
        return myBorderValue;
    }

    switch (myBorderType)
    {
        case UT_VOXELBORDER_CONSTANT:   wrap = UT_WRAP_BORDER; break;
        case UT_VOXELBORDER_REPEAT:     wrap = UT_WRAP_REPEAT; break;
        case UT_VOXELBORDER_STREAK:     wrap = UT_WRAP_CLAMP; break;
        
        // NOTE: This is incorrect!
        case UT_VOXELBORDER_EXTRAP:     wrap = UT_WRAP_CLAMP; break;
        default:                        wrap = UT_WRAP_CLAMP; break;
    }

    
    result = 0;

    radius *= filter.getSupport();
    iradius = radius == 0 ? 1000 : 0.5F / radius;

    // Make a local copy of the position so that we can modify it.
    tpos = pos;
    for (i = 0; i < 3; i++)
    {
        tpos[i] = tpos[i]*myRes[i];
        win[i].initWeights(tpos[i], radius, myRes[i], wrap);
        if (win[i].getSize() <= 0)
            return result;
    }

    // Create stack storage for weight arrays
    wdata = (float *)UTstackAlloc(
            (win[0].getSize()+win[1].getSize()+win[2].getSize())*
            sizeof(float));
    weights[0] = wdata;
    weights[1] = weights[0] + win[0].getSize();
    weights[2] = weights[1] + win[1].getSize();

    visible = 1;
    for (i = 0; i < 3; i++)
    {
        win[i].fillWeights(weights[i], filter, tpos[i], iradius,
                myRes[i], wrap);
        start[i] = win[i].getStart() % myRes[i];
        size[i] = win[i].getSize();

        UT_ASSERT(start[i] >= 0);
        UT_ASSERT(size[i] <= myRes[i]);

        end[i] = start[i] + size[i];
        visible *= win[i].getVisible();
    }

    // Accumulate filtered results
    if (weights[0] && weights[1] && weights[2])
    {
        pstart[2] = firstTile(start[2], end[2], myRes[2]);
        while (pstart[2] < end[2])
        {
            pend[2] = nextTile(tz, pstart[2], start[2], end[2], myRes[2]);
            pstart[1] = firstTile(start[1], end[1], myRes[1]);
            while (pstart[1] < end[1])
            {
                pend[1] = nextTile(ty, pstart[1], start[1], end[1], myRes[1]);
                pstart[0] = firstTile(start[0], end[0], myRes[0]);
                while (pstart[0] < end[0])
                {
                    pend[0] = nextTile(tx, pstart[0], start[0], end[0], myRes[0]);
                    tile = getTile(tx, ty, tz);
                    UT_ASSERT(tile);
                    tile->weightedSum(pstart, pend, weights, start, result);
                    pstart[0] = pend[0];
                }
                pstart[1] = pend[1];
            }
            pstart[2] = pend[2];
        }
    }
    UTstackFree(wdata);

    if (visible < 1)
    {
        result += (1-visible)*myBorderValue;
    }

    return result;
}

template <typename T>
void
UT_VoxelArray<T>::resample(const UT_VoxelArray<T> &src,
                          UT_FilterType filtertype,
                          float filterwidthscale)
{
    fpreal               radius;
    UT_Filter           *filter;
    UT_Vector3           pos;

    filter = UT_Filter::getFilter(filtertype);

    radius = SYSmax( getXRes() / (fpreal)src.getXRes(),
                     getYRes() / (fpreal)src.getYRes(),
                     getZRes() / (fpreal)src.getZRes(),
                     1.0f );
    radius *= 0.5 * filterwidthscale;

    resamplethread(src, filter, radius);

    UT_Filter::releaseFilter(filter);
}

template <typename T>
void
UT_VoxelArray<T>::flattenPartial(T *flatarray, int ystride32, int zstride32,
                                const UT_JobInfo &info) const
{
    UT_VoxelArrayIterator<T>    vit;
    int64                       ystride = ystride32;
    int64                       zstride = zstride32;

    // Const cast.
    vit.setArray((UT_VoxelArray<T> *)this);
    vit.splitByTile(info);

    for (vit.rewind(); !vit.atEnd(); vit.advance())
    {
        flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride] = vit.getValue();
    }
}

template <typename T>
void
UT_VoxelArray<T>::extractFromFlattenedPartial(const T *flatarray, 
                                int ystride32, int zstride32,
                                const UT_JobInfo &info)
{
    UT_VoxelArrayIterator<T>    vit;
    int64                       ystride = ystride32;
    int64                       zstride = zstride32;

    vit.setArray(this);
    vit.splitByTile(info);
    vit.setCompressOnExit(true);

    for (vit.rewind(); !vit.atEnd(); vit.advance())
    {
        vit.setValue(flatarray[vit.x() + vit.y()*ystride + vit.z()*zstride]);
    }
}

template <typename T>
void
UT_VoxelArray<T>::resamplethreadPartial(const UT_VoxelArray<T> &src,
                        const UT_Filter *filter, float radius,
                        const UT_JobInfo &info)
{
    UT_VoxelArrayIterator<T>    vit;
    UT_Vector3                  pos;
    UT_Vector3                  ratio;

    vit.setArray(this);
    vit.splitByTile(info);
    vit.setCompressOnExit(true);

    ratio.x() = 1.0f / getXRes();
    ratio.y() = 1.0f / getYRes();
    ratio.z() = 1.0f / getZRes();

    for (vit.rewind(); !vit.atEnd(); vit.advance())
    {
        pos.x() = vit.x()+0.5f;
        pos.y() = vit.y()+0.5f;
        pos.z() = vit.z()+0.5f;
        pos *= ratio;

        vit.setValue(src.evaluate(pos, *filter, radius));
    }
}

template <typename T>
bool
UT_VoxelArray<T>::posToIndex(UT_Vector3 pos, int &x, int &y, int &z) const
{
    // We go from the position in the unit cube into the index.
    // The center of cells must map to the exact integer indices.
    pos.x() *= myRes[0];
    pos.y() *= myRes[1];
    pos.z() *= myRes[2];

    // The centers of cells are now mapped .5 too high.  Ie, the true
    // center of cell index (0,0,0) would be (0.5,0.5,0.5)  This, however,
    // is exactly what we want for rounding.
    x = (int) SYSfloor(pos.x());
    y = (int) SYSfloor(pos.y());
    z = (int) SYSfloor(pos.z());

    // Determine if out of bounds.
    return isValidIndex(x, y, z);
}

template <typename T>
bool
UT_VoxelArray<T>::posToIndex(UT_Vector3 pos, UT_Vector3 &ipos) const
{
    // We go from the position in the unit cube into the index.
    // The center of cells must map to the exact integer indices.
    pos.x() *= myRes[0];
    pos.y() *= myRes[1];
    pos.z() *= myRes[2];

    // The centers of cells are now mapped .5 too high.  Ie, the true
    // center of cell index (0,0,0) would be (0.5,0.5,0.5)
    pos.x() -= 0.5;
    pos.y() -= 0.5;
    pos.z() -= 0.5;

    ipos = pos;

    // Determine if out of bounds.
    if (pos.x() < 0 || pos.x() >= myRes[0] ||
        pos.y() < 0 || pos.y() >= myRes[1] ||
        pos.z() < 0 || pos.z() >= myRes[2])
        return false;

    return true;
}

template <typename T>
bool
UT_VoxelArray<T>::indexToPos(int x, int y, int z, UT_Vector3F &pos) const
{
    pos.x() = x;
    pos.y() = y;
    pos.z() = z;

    // Move the indices to the centers of the cells.
    pos.x() += 0.5;
    pos.y() += 0.5;
    pos.z() += 0.5;

    // And scale into the unit cube.
    if (myRes[0])
        pos.x() /= (fpreal) myRes[0];
    if (myRes[1])
        pos.y() /= (fpreal) myRes[1];
    if (myRes[2])
        pos.z() /= (fpreal) myRes[2];

    // Return true if the original coordinates were in range.
    return isValidIndex(x, y, z);
}

template <typename T>
bool
UT_VoxelArray<T>::indexToPos(int x, int y, int z, UT_Vector3D &pos) const
{
    pos.x() = x;
    pos.y() = y;
    pos.z() = z;

    // Move the indices to the centers of the cells.
    pos.x() += 0.5;
    pos.y() += 0.5;
    pos.z() += 0.5;

    // And scale into the unit cube.
    if (myRes[0])
        pos.x() /= (fpreal) myRes[0];
    if (myRes[1])
        pos.y() /= (fpreal) myRes[1];
    if (myRes[2])
        pos.z() /= (fpreal) myRes[2];

    // Return true if the original coordinates were in range.
    return isValidIndex(x, y, z);
}

template <typename T>
void
UT_VoxelArray<T>::setBorder(UT_VoxelBorderType type, T t)
{
    myBorderType = type;
    myBorderValue = t;
}

template <typename T>
void
UT_VoxelArray<T>::setBorderScale(T sx, T sy, T sz)
{
    myBorderScale[0] = sx;
    myBorderScale[1] = sy;
    myBorderScale[2] = sz;
}

template <typename T>
void
UT_VoxelArray<T>::collapseAllTiles()
{
    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        myTiles[i]->tryCompress(getCompressionOptions());
    }
}

template <typename T>
void
UT_VoxelArray<T>::expandAllTiles()
{
    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        myTiles[i]->uncompress();
    }
}

template <typename T>
void
UT_VoxelArray<T>::expandAllNonConstTiles()
{
    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        if (!myTiles[i]->isConstant())
            myTiles[i]->uncompress();
    }
}

template <typename T>
void
UT_VoxelArray<T>::saveData(ostream &os)
{
    T           cval;
    char        version;

    // First determine if we are fully constant.
    if (isConstant(&cval))
    {
        // Save a constant array.
        version = 0;
        UTwrite(os, &version, 1);
        UTwrite(os, &cval);
        return;
    }

    // Compressed tiles.
    version = 1;
    UTwrite(os, &version, 1);

    // Store list of compression types
    UT_VoxelTile<T>::saveCompressionTypes(os);

    int         i, ntiles;

    ntiles = numTiles();
    for (i = 0; i < ntiles; i++)
    {
        myTiles[i]->save(os);
    }
}


template <typename T>
void
UT_VoxelArray<T>::loadData(UT_IStream &is)
{
    T           cval;
    char        version;

    is.readChar(version);

    // First determine if we are fully constant.
    if (version == 0)
    {
        // Save a constant array.
        is.read(&cval);

        constant(cval);
        return;
    }

    if (version == 1)
    {
        UT_IntArray             compressions;

        // Store list of compression types
        UT_VoxelTile<T>::loadCompressionTypes(is, compressions);

        int             i, ntiles;

        ntiles = numTiles();
        for (i = 0; i < ntiles; i++)
        {
            myTiles[i]->load(is, compressions);
        }
    }
}



//
// UT_VoxelMipMap functions
//
template <typename T>
UT_VoxelMipMap<T>::UT_VoxelMipMap()
{
    initializePrivate();
}

template <typename T>
UT_VoxelMipMap<T>::~UT_VoxelMipMap()
{
    destroyPrivate();
}

template <typename T>
UT_VoxelMipMap<T>::UT_VoxelMipMap(const UT_VoxelMipMap<T> &src)
{
    initializePrivate();

    *this = src;
}

template <typename T>
const UT_VoxelMipMap<T> &
UT_VoxelMipMap<T>::operator=(const UT_VoxelMipMap<T> &src)
{
    UT_VoxelArray<T>    *levels;
    int                  level, i;
    
    if (&src == this)
        return;

    destroyPrivate();

    // We do a deep copy.  We have learned from UT_String that
    // shallow copies will destroy us in the end.
    myOwnBase = true;
    myBaseLevel = new UT_VoxelArray<T>;
    *myBaseLevel = *src.myBaseLevel;

    myNumLevels = src.myNumLevels;

    for (i = 0; i < src.myLevels.entries(); i++)
    {
        levels = new UT_VoxelArray<T> *[myNumLevels];
        myLevels.append(levels);
        
        for (level = 0; level < myNumLevels; level++)
        {
            levels[level] = new UT_VoxelArray<T>;
            *levels[level] = *src.myLevels[i][level];
        }
    }
}

template <typename T>
void
UT_VoxelMipMap<T>::build(UT_VoxelArray<T> *baselevel,
                         mipmaptype function)
{
    UT_RefArray<mipmaptype>     functions;

    functions.append(function);
    build(baselevel, functions);
}

template <typename T>
void
UT_VoxelMipMap<T>::build(UT_VoxelArray<T> *baselevel,
                         const UT_RefArray<mipmaptype> &functions)
{
    destroyPrivate();

    // Setup our baselevel.
    myBaseLevel = baselevel;
    myOwnBase = false;

    // Now build from bottom up.  First, calculate the number
    // of levels we'll need.
    myNumLevels = 0;
    int                 maxres, level;

    maxres = SYSmax(myBaseLevel->getXRes(), myBaseLevel->getYRes());
    maxres = SYSmax(maxres, myBaseLevel->getZRes());

    //cerr << "Max res " << maxres << endl;

    while (maxres > 1)
    {
        myNumLevels++;
        maxres++;
        maxres >>= 1;
    }

    //cerr << "Levels: " << myNumLevels << endl;

    int                  x, y, z, ix, iy, iz;
    int                  xres, yres, zres;
    UT_VoxelArray<T>    *lastlevel;
    UT_VoxelArray<T>    **levels;
    T                    x1, x2, y1, y2, z1, z2;
    mipmaptype           function;

    // Create our desired voxel array list.
    for (int i = 0; i < functions.entries(); i++)
    {
        function = functions(i);

        levels = new UT_VoxelArray<T> *[myNumLevels];
        myLevels.append(levels);

        lastlevel = myBaseLevel;
        for (level = myNumLevels-1; level >= 0; level--)
        {
            // Find our resolutions.
            xres = (lastlevel->getXRes() + 1) >> 1;
            yres = (lastlevel->getYRes() + 1) >> 1;
            zres = (lastlevel->getZRes() + 1) >> 1;

            levels[level] = new UT_VoxelArray<T>;
            levels[level]->size(xres, yres, zres);

            iz = 1;
            for (z = 0; z < zres; z++)
            {
                if (z == zres-1)
                {
                    if (2*z + 1 >= lastlevel->getZRes())
                        iz = 0;
                }
                iy = 1;
                for (y = 0; y < yres; y++)
                {
                    if (y == yres-1)
                    {
                        if (2*y + 1 >= lastlevel->getYRes())
                            iy = 0;
                    }
                    ix = 1;
                    for (x = 0; x < xres; x++)
                    {
                        if (x == xres-1)
                        {
                            if (2*x + 1 >= lastlevel->getXRes())
                                ix = 0;
                        }
                        x1 = (*lastlevel)(2*x,    2*y,    2*z);
                        x2 = (*lastlevel)(2*x+ix, 2*y,    2*z);
                        y1 = mixValues(x1, x2, function);
                        
                        x1 = (*lastlevel)(2*x,    2*y+iy, 2*z);
                        x2 = (*lastlevel)(2*x+ix, 2*y+iy, 2*z);
                        y2 = mixValues(x1, x2, function);

                        z1 = mixValues(y1, y2, function);

                        x1 = (*lastlevel)(2*x,    2*y,    2*z+iz);
                        x2 = (*lastlevel)(2*x+ix, 2*y,    2*z+iz);
                        y1 = mixValues(x1, x2, function);
                        
                        x1 = (*lastlevel)(2*x,    2*y+iy, 2*z+iz);
                        x2 = (*lastlevel)(2*x+ix, 2*y+iy, 2*z+iz);
                        y2 = mixValues(x1, x2, function);

                        z2 = mixValues(y1, y2, function);

                        levels[level]->setValue(x, y, z, 
                                                  mixValues(z1, z2, function));
                    }
                }
            }

            lastlevel = levels[level];    
        }
    }
}

template <typename T>
T
UT_VoxelMipMap<T>::mixValues(T t1, T t2, mipmaptype function) const
{
    switch (function)
    {
        case MIPMAP_MAXIMUM:
            return SYSmax(t1, t2);

        case MIPMAP_AVERAGE:
            return (t1 + t2) / 2;

        case MIPMAP_MINIMUM:
            return SYSmin(t1, t2);
    }

    return t1;
}

template <typename T>
int64
UT_VoxelMipMap<T>::getMemoryUsage() const
{
    int64        totalarraysize = 0;

    for (int j = 0; j < myLevels.entries(); j++)
    {
        for( int i = 0; i < myNumLevels; i++ )
            totalarraysize += myLevels[j][i]->getMemoryUsage();
    }

    return sizeof(*this) + totalarraysize;
}

template <typename T>
void
UT_VoxelMipMap<T>::traverseTopDown(bool (*function)(const UT_RefArray<T> &,
                                                    const UT_BoundingBox &,
                                                    bool, void *),
                                    void *data) const
{
    doTraverse(0, 0, 0, 0, function, data);
}

template <typename T>
void
UT_VoxelMipMap<T>::doTraverse(int x, int y, int z, int level,
                                bool (*function)(const UT_RefArray<T> &,
                                                 const UT_BoundingBox &,
                                                 bool, void *),
                                void *data) const
{
    UT_RefArray<T>       tval;
    bool                 isfinal;
    UT_VoxelArray<T>    *vox;
    int                  shift;
    int                  i;

    if (level == myNumLevels)
    {
        isfinal = true;
        vox = myBaseLevel;
        tval.append((*vox)(x, y, z));
        for (i = 1; i < myLevels.entries(); i++)
            tval.append(tval(0));
    }
    else
    {
        isfinal = false;
        for (i = 0; i < myLevels.entries(); i++)
        {
            vox = myLevels[i][level];
            tval.append((*vox)(x, y, z));
        }
    }

    shift = myNumLevels - level;

    UT_BoundingBox              box;

    box.initBounds(SYSmin(x << shift, myBaseLevel->getXRes()),
                   SYSmin(y << shift, myBaseLevel->getYRes()),
                   SYSmin(z << shift, myBaseLevel->getZRes()));
    box.enlargeBounds(SYSmin((x+1) << shift, myBaseLevel->getXRes()),
                      SYSmin((y+1) << shift, myBaseLevel->getYRes()),
                      SYSmin((z+1) << shift, myBaseLevel->getZRes()));

    if (!function(tval, box, isfinal, data))
    {
        // Function asked for an early exit.
        // Give it.
        return;
    }

    // We now want to proceed to the next level.
    if (isfinal)
        return;

    level++;
    shift--;
    x <<= 1;
    y <<= 1;
    z <<= 1;

    bool                xinc, yinc, zinc;

    // Determine which of our children are valid.
    if ( ((x+1) << shift) < myBaseLevel->getXRes() )
        xinc = true;
    else
        xinc = false;
    if ( ((y+1) << shift) < myBaseLevel->getYRes() )
        yinc = true;
    else
        yinc = false;
    if ( ((z+1) << shift) < myBaseLevel->getZRes() )
        zinc = true;
    else
        zinc = false;
    
    // Traverse all valid children.
    doTraverse(x, y, z, level, function, data);

    if (xinc)
    {
        doTraverse(x+1, y, z, level, function, data);
        if (yinc)
        {
            doTraverse(x+1, y+1, z, level, function, data);
            if (zinc)
                doTraverse(x+1, y+1, z+1, level, function, data);
        }
        if (zinc)
        {
            doTraverse(x+1, y, z+1, level, function, data);
        }
    }
    if (yinc)
    {
        doTraverse(x, y+1, z, level, function, data);
        if (zinc)
            doTraverse(x, y+1, z+1, level, function, data);
    }
    if (zinc)
        doTraverse(x, y, z+1, level, function, data);
}

template <typename T>
void
UT_VoxelMipMap<T>::initializePrivate()
{
    myNumLevels = 0;
    myBaseLevel = 0;
    myOwnBase = false;
}

template <typename T>
void
UT_VoxelMipMap<T>::destroyPrivate()
{
    int         level, i;

    for (i = 0; i < myLevels.entries(); i++)
    {
        for (level = 0; level < myNumLevels; level++)
            delete myLevels[i][level];
        delete [] myLevels[i];
    }
    myLevels.entries(0);

    if (myOwnBase)
        delete myBaseLevel;

    initializePrivate();
}

template <typename T>
UT_VoxelArrayIterator<T>::UT_VoxelArrayIterator()
{
    myArray = 0;
    myHandle.resetHandle();
    myCurTile = -1;
    myShouldCompressOnExit = false;
    myUseTileList = false;
    myJobInfo = 0;
}

template <typename T>
UT_VoxelArrayIterator<T>::UT_VoxelArrayIterator(UT_VoxelArray<T> *vox)
{
    myShouldCompressOnExit = false;
    myUseTileList = false;
    myJobInfo = 0;
    setArray(vox);
}

template <typename T>
UT_VoxelArrayIterator<T>::UT_VoxelArrayIterator(UT_COWReadHandle<UT_VoxelArray<T> > handle)
{
    myShouldCompressOnExit = false;
    myUseTileList = false;
    myJobInfo = 0;
    setHandle(handle);
}

template <typename T>
UT_VoxelArrayIterator<T>::~UT_VoxelArrayIterator()
{
}

template <typename T>
void
UT_VoxelArrayIterator<T>::setPartialRange(int idx, int numranges)
{
    int         numtiles;
    fpreal      tileperrange;
    
    // Must already have been set.
    UT_ASSERT(myArray);
    if (!myArray)
        return;

    // Divide our tiles in to num ranges groups.
    numtiles = myArray->numTiles();

    // If we are using a tile list, use it instead
    if (myUseTileList)
        numtiles = myTileList.entries();

    // Trivial case.
    if (!numtiles)
    {
        myTileStart = 0;
        myTileEnd = 0;
        return;
    }

    // Sanity check idx.
    if (idx < 0 || idx >= numranges)
    {
        UT_ASSERT(!"Idx out of bounds");
        myTileStart = -1;
        myTileEnd = -1;
        return;
    }
    
    // Sanity test numranges
    if (numranges < 1)
    {
        UT_ASSERT(!"Invalid range count!");
        numranges = 1;
    }

    // Compute tile per range.
    // We use floating point so, if you have 15 tiles and
    // 16 processors, you don't decide to do all 15
    // tiles on one processor.
    tileperrange = (fpreal)numtiles / (fpreal)numranges;
    
    myTileStart = (int) SYSfloor(idx * tileperrange);
    myTileEnd = (int) SYSfloor((idx+1) * tileperrange); 

    // Special case to ensure we always get the last tile,
    // despite the evils of floating point.
    if (idx == numranges-1)
        myTileEnd = numtiles;

    UT_ASSERT(myTileStart >= 0);
    UT_ASSERT(myTileEnd <= numtiles);
}

template <typename T>
void
UT_VoxelArrayIterator<T>::splitByTile(const UT_JobInfo &info)
{
    // Must already have been set.
    UT_ASSERT(myArray);
    if (!myArray)
        return;

    // If there is one thread, don't bother
    if (info.numJobs() == 1)
        myJobInfo = 0;
    else
        myJobInfo = &info;
}

template <typename T>
void
UT_VoxelArrayIterator<T>::restrictToBBox(const UT_BoundingBox &bbox)
{
    UT_Vector3          pmin, pmax, vmin, vmax;

    pmin = bbox.minvec();
    pmax = bbox.maxvec();

    pmin.x() = SYSmax(pmin.x(), 0.0f);
    pmin.y() = SYSmax(pmin.y(), 0.0f);
    pmin.z() = SYSmax(pmin.z(), 0.0f);
    pmax.x() = SYSmin(pmax.x(), 1.0f);
    pmax.y() = SYSmin(pmax.y(), 1.0f);
    pmax.z() = SYSmin(pmax.z(), 1.0f);

    myArray->posToIndex(pmin, vmin);
    myArray->posToIndex(pmax, vmax);

    restrictToBBox(SYSfloor(vmin.x()), SYSceil(vmax.x()),
                    SYSfloor(vmin.y()), SYSceil(vmax.y()),
                    SYSfloor(vmin.z()), SYSceil(vmax.z()));
}

template <typename T>
void
UT_VoxelArrayIterator<T>::restrictToBBox(int xmin, int xmax,
                                        int ymin, int ymax,
                                        int zmin, int zmax)
{
    int                 xres, yres, zres, x, y, z;

    xres = myArray->getXRes();
    yres = myArray->getYRes();
    zres = myArray->getZRes();

    myTileList.entries(0);
    myUseTileList = true;
    // Make sure we have any tiles that are valid.
    if (xmin < xres && xmax >= 0 &&
        ymin < yres && ymax >= 0 &&
        zmin < zres && zmax >= 0)
    {
        // Clamp into valid range.
        myArray->clampIndex(xmin, ymin, zmin);
        myArray->clampIndex(xmax, ymax, zmax);

        // Convert to tiles...
        xmin >>= TILEBITS;
        ymin >>= TILEBITS;
        zmin >>= TILEBITS;
        xmax >>= TILEBITS;
        ymax >>= TILEBITS;
        zmax >>= TILEBITS;

        // Check if we are all accounted for.
        if (myArray->numTiles() == (xmax-xmin+1)*(ymax-ymin+1)*(zmax-zmin+1))
        {
            UT_ASSERT(xmin == 0 && ymin == 0 && zmin == 0);
            // No need for a tile list, just run everything!
            myUseTileList = false;
        }
        else
        {
            // Iterate over all active tiles, adding to our list.
            for (z = zmin; z <= zmax; z++)
            {
                for (y = ymin; y <= ymax; y++)
                {
                    for (x = xmin; x <= xmax; x++)
                    {
                        myTileList.append(myArray->xyzTileToLinear(x, y, z));
                    }
                }
            }
        }
    }

    // Recompute our ranges.
    myCurTile = -1;
    setPartialRange(0, 1);
}
    

template <typename T>
void
UT_VoxelArrayIterator<T>::rewind()
{
    // Ensure we have at least one tile in each direction
    if (!myArray ||
        !myArray->getRes(0) || !myArray->getRes(1) || !myArray->getRes(2))
    {
        myCurTile = -1;
        return;
    }

    if (myUseTileList)
    {
        if (myJobInfo)
            myCurTileListIdx = myJobInfo->nextTask();
        else
            myCurTileListIdx = myTileStart;
        if (myCurTileListIdx < 0 || myCurTileListIdx >= myTileEnd)
        {
            myCurTile = -1;
            return;
        }
        myCurTile = myTileList(myCurTileListIdx);
    }
    else
    {
        if (myJobInfo)
            myCurTile = myJobInfo->nextTask();
        else
            myCurTile = myTileStart;
        // If this is an empty range, we quit right away.
        if (myCurTile < 0 || myCurTile >= myTileEnd)
        {
            myCurTile = -1;
            return;
        }
    }

    UT_VoxelTile<T>     *tile;
    
    tile = myArray->getLinearTile(myCurTile);

    // Get the tile position
    if (myCurTile)
    {
        myArray->linearTileToXYZ(myCurTile, 
                                 myTilePos[0], myTilePos[1], myTilePos[2]);
        myPos[0] = TILESIZE * myTilePos[0];
        myPos[1] = TILESIZE * myTilePos[1];
        myPos[2] = TILESIZE * myTilePos[2];
    }
    else
    {
        myTilePos[0] = 0;
        myTilePos[1] = 0;
        myTilePos[2] = 0;
        myPos[0] = 0;
        myPos[1] = 0;
        myPos[2] = 0;
    }

    myTileLocalPos[0] = 0;
    myTileLocalPos[1] = 0;
    myTileLocalPos[2] = 0;

    myTileSize[0] = tile->xres();
    myTileSize[1] = tile->yres();
    myTileSize[2] = tile->zres();
}

template <typename T>
void
UT_VoxelArrayIterator<T>::advanceTile()
{
    if (myUseTileList)
    {
        if (getCompressOnExit())
        {
            // Verify our last tile was a legitimate one.
            if (myCurTile >= 0 && myCurTileListIdx < myTileEnd)
            {
                myArray->getLinearTile(myCurTile)->tryCompress(myArray->getCompressionOptions());
            }
        }

        // Advance our myCurTileListIdx and rebuild from 
        if (myJobInfo)
            myCurTileListIdx = myJobInfo->nextTask();
        else
            myCurTileListIdx++;
        if (myCurTileListIdx >= myTileEnd)
        {
            myCurTile = -1;
            return;
        }

        myCurTile = myTileList(myCurTileListIdx);

        myArray->linearTileToXYZ(myCurTile, 
                                 myTilePos[0], myTilePos[1], myTilePos[2]);

        UT_VoxelTile<T>         *tile;

        tile = myArray->getLinearTile(myCurTile);
        myTileLocalPos[0] = 0;
        myTileLocalPos[1] = 0;
        myTileLocalPos[2] = 0;
        myTileSize[0] = tile->xres();
        myTileSize[1] = tile->yres();
        myTileSize[2] = tile->zres();

        myPos[0] = TILESIZE * myTilePos[0];
        myPos[1] = TILESIZE * myTilePos[1];
        myPos[2] = TILESIZE * myTilePos[2];
        return;
    }

    // If requested, compress our last tile.
    if (getCompressOnExit())
    {
        // Verify our last tile was a legitimate one.
        if (myCurTile >= 0 && myCurTile < myTileEnd)
        {
            myArray->getLinearTile(myCurTile)->tryCompress(myArray->getCompressionOptions());
        }
    }
    
    if (myJobInfo)
    {
        myCurTile = myJobInfo->nextTask();
        if (myCurTile >= myTileEnd)
        {
            myCurTile = -1;
            return;
        }
        myArray->linearTileToXYZ(myCurTile, 
                                 myTilePos[0], myTilePos[1], myTilePos[2]);
    }
    else
    {
        // Advance our myTilePos, then rebuild everything from there.
        myTilePos[0]++;
        if (myTilePos[0] >= myArray->getTileRes(0))
        {
            myTilePos[0] = 0;
            myTilePos[1]++;
            if (myTilePos[1] >= myArray->getTileRes(1))
            {
                myTilePos[1] = 0;
                myTilePos[2]++;
                if (myTilePos[2] >= myArray->getTileRes(2))
                {
                    // All done!
                    myCurTile = -1;
                    return;
                }
            }
        }
        myCurTile = myArray->xyzTileToLinear(myTilePos[0], myTilePos[1], myTilePos[2]);
    }

    UT_VoxelTile<T>     *tile;

    // Rebuild our settings from this tile.

    // See if we hit our end.
    if (myCurTile >= myTileEnd)
    {
        myCurTile = -1;
        return;
    }

    tile = myArray->getLinearTile(myCurTile);
    myTileLocalPos[0] = 0;
    myTileLocalPos[1] = 0;
    myTileLocalPos[2] = 0;
    myTileSize[0] = tile->xres();
    myTileSize[1] = tile->yres();
    myTileSize[2] = tile->zres();

    myPos[0] = TILESIZE * myTilePos[0];
    myPos[1] = TILESIZE * myTilePos[1];
    myPos[2] = TILESIZE * myTilePos[2];
}


template <typename T>
void
UT_VoxelArrayIterator<T>::skipToEndOfTile()
{
    myTileLocalPos[0] = myTileSize[0];
    myTileLocalPos[1] = myTileSize[1];
    myTileLocalPos[2] = myTileSize[2];
}

template <typename T>
template <typename OP>
void
UT_VoxelArrayIterator<T>::applyOperation(OP op,
                            const UT_VoxelArray<T> &a)
{
    UT_ASSERT(myArray->isMatching(a));

    rewind();

    while (!atEnd())
    {
        UT_VoxelTile<T>         *atile, *tile;

        tile = myArray->getLinearTile(myCurTile);
        atile = a.getLinearTile(myCurTile);

        if (!atile->isSimpleCompression())
        {
            // Have to use getValue...
            for (int z = 0; z < tile->zres(); z++)
                for (int y = 0; y < tile->yres(); y++)
                    for (int x = 0; x < tile->xres(); x++)
                    {
                        T               aval, val;

                        val = tile->operator()(x, y, z);
                        aval = atile->operator()(x, y, z);

                        val = op(val, aval);

                        tile->setValue(x, y, z, val);
                    }
        }
        else
        {
            if (!tile->isSimpleCompression())
                tile->uncompress();

            int         ainc = atile->isConstant() ? 0 : 1;

            if (tile->isConstant())
            {
                if (!ainc)
                {
                    // Woohoo!  Too simple!
                    T           aval, val;

                    val = tile->rawData()[0];
                    aval = atile->rawData()[0];

                    val = op(val, aval);

                    tile->rawData()[0] = val;
                }
                else
                {
                    tile->uncompress();
                }
            }

            // In case we uncompressed ourself above...
            // we have to test again
            if (!tile->isConstant())
            {
                const T         *aval;
                T               *val;

                val = tile->rawData();
                aval = atile->rawData();

                int             n = tile->numVoxels();
                for (int i = 0; i < n; i++)
                {
                    val[i] = op(val[i], *aval);
                    aval += ainc;
                }
            }
        }

        advanceTile();
    }
}

template <typename T>
template <typename OP>
void
UT_VoxelArrayIterator<T>::applyOperation(OP op, T a)
{
    rewind();

    while (!atEnd())
    {
        UT_VoxelTile<T>         *tile;

        tile = myArray->getLinearTile(myCurTile);

        if (!tile->isSimpleCompression())
            tile->uncompress();

        if (tile->isConstant())
        {
            T           val;

            val = tile->rawData()[0];

            val = op(val, a);

            tile->rawData()[0] = val;
        }
        else
        {
            T           *val;

            val = tile->rawData();

            int         n = tile->numVoxels();
            for (int i = 0; i < n; i++)
            {
                val[i] = op(val[i], a);
            }
        }

        advanceTile();
    }
}

template <typename T>
template <typename OP>
void
UT_VoxelArrayIterator<T>::applyOperation(OP op,
                            const UT_VoxelArray<T> &a,
                            const UT_VoxelArray<T> &b)
{
    UT_ASSERT(myArray->isMatching(a));
    UT_ASSERT(myArray->isMatching(b));

    rewind();

    while (!atEnd())
    {
        UT_VoxelTile<T>         *atile, *btile, *tile;

        tile = myArray->getLinearTile(myCurTile);
        atile = a.getLinearTile(myCurTile);
        btile = b.getLinearTile(myCurTile);

        if (!atile->isSimpleCompression() || !btile->isSimpleCompression())
        {
            // Have to use getValue...
            for (int z = 0; z < tile->zres(); z++)
                for (int y = 0; y < tile->yres(); y++)
                    for (int x = 0; x < tile->xres(); x++)
                    {
                        T               aval, bval, val;

                        val = tile->operator()(x, y, z);
                        aval = atile->operator()(x, y, z);
                        bval = btile->operator()(x, y, z);

                        val = op(val, aval, bval);

                        tile->setValue(x, y, z, val);
                    }
        }
        else
        {
            if (!tile->isSimpleCompression())
                tile->uncompress();

            int         ainc = atile->isConstant() ? 0 : 1;
            int         binc = btile->isConstant() ? 0 : 1;

            if (tile->isConstant())
            {
                if (!ainc && !binc)
                {
                    // Woohoo!  Too simple!
                    T           aval, bval, val;

                    val = tile->rawData()[0];
                    aval = atile->rawData()[0];
                    bval = btile->rawData()[0];

                    val = op(val, aval, bval);

                    tile->rawData()[0] = val;
                }
                else
                {
                    tile->uncompress();
                }
            }

            // In case we uncompressed ourself above...
            // we have to test again
            if (!tile->isConstant())
            {
                const T         *aval, *bval;
                T               *val;

                val = tile->rawData();
                aval = atile->rawData();
                bval = btile->rawData();

                int             n = tile->numVoxels();
                for (int i = 0; i < n; i++)
                {
                    val[i] = op(val[i], *aval, *bval);
                    aval += ainc;
                    bval += binc;
                }
            }
        }

        advanceTile();
    }
}

///
/// UT_VoxelTileIterator implementation
///
template <typename T>
UT_VoxelTileIterator<T>::UT_VoxelTileIterator()
{
    myCurTile = 0;
    myArray = 0;
    myAtEnd = true;
    myShouldCompressOnExit = false;
}

template <typename T>
UT_VoxelTileIterator<T>::UT_VoxelTileIterator(const UT_VoxelArrayIterator<T> &vit)
{
    myCurTile = 0;
    myArray = 0;
    myAtEnd = true;
    myShouldCompressOnExit = false;
    setTile(vit);
}

template <typename T>
template <typename S>
UT_VoxelTileIterator<T>::UT_VoxelTileIterator(const UT_VoxelArrayIterator<S> &vit, UT_VoxelArray<T> *array)
{
    myCurTile = 0;
    myArray = 0;
    myAtEnd = true;
    myShouldCompressOnExit = false;
    setTile(vit, array);
}

template <typename T>
UT_VoxelTileIterator<T>::~UT_VoxelTileIterator()
{
}

template <typename T>
void
UT_VoxelTileIterator<T>::rewind()
{
    // Ensure we have at least one voxel in each direction
    if (!myCurTile ||
        !myCurTile->xres() || !myCurTile->yres() || !myCurTile->zres())
    {
        myCurTile = 0;
        return;
    }

    myPos[0] = myTileStart[0];
    myPos[1] = myTileStart[1];
    myPos[2] = myTileStart[2];

    myTileLocalPos[0] = 0;
    myTileLocalPos[1] = 0;
    myTileLocalPos[2] = 0;

    myTileSize[0] = myCurTile->xres();
    myTileSize[1] = myCurTile->yres();
    myTileSize[2] = myCurTile->zres();

    myAtEnd = false;
}

template <typename T>
void
UT_VoxelTileIterator<T>::advanceTile()
{
    if (getCompressOnExit())
    {
        // Verify our last tile was a legitimate one.
        if (myCurTile)
        {
            myCurTile->tryCompress(myArray->getCompressionOptions());
        }
    }
    myAtEnd = true;
}

template <typename T>
template <typename OP>
bool
UT_VoxelTileIterator<T>::reduceOperation(OP &op)
{
    rewind();

    if (!myCurTile->isSimpleCompression())
    {
        // Have to use getValue...
        for (int z = 0; z < myTileSize[2]; z++)
            for (int y = 0; y < myTileSize[1]; y++)
                for (int x = 0; x < myTileSize[0]; x++)
                {
                    T           val;

                    val = myCurTile->operator()(x, y, z);

                    if (!op.reduce(val))
                        return false;
                }
    }
    else if (myCurTile->isConstant())
    {
        // Woohoo!  Too simple!
        T               val;

        val = myCurTile->rawData()[0];
        int             n = myCurTile->numVoxels();

        if (!op.reduceMany(val, n))
            return false;
    }
    else
    {
        T               *val;

        val = myCurTile->rawData();

        int             n = myCurTile->numVoxels();
        for (int i = 0; i < n; i++)
        {
            if (!op.reduce(val[i]))
                return false;
        }
    }
    return true;
}

///
/// UT_VoxelProbe methods
///

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::UT_VoxelProbe()
{
    myCurLine = 0;
    myAllocCacheLine = 0;
    myDirty = false;

    myArray = 0;
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::UT_VoxelProbe(UT_VoxelArray<T> *vox, int prex, int postx)
{
    // A sure signal it hasn't been reset.
    myCurLine = 0;
    myAllocCacheLine = 0;
    myDirty = false;

    setArray(vox, prex, postx);
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::~UT_VoxelProbe()
{
    if (DoWrite)
    {
        if (!TestForWrites || myDirty)
        {
            // Final write...
            writeCacheLine();
        }
    }
    delete [] myAllocCacheLine;
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
void
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::setArray(UT_VoxelArray<T> *vox, int prex, int postx)
{
    // Round up our pre and post
    int                 prepad, postpad;

    myCurLine = 0;

    prepad = (prex - 3) / 4;
    postpad = (postx + 3) / 4;

    myAllocCacheLine = new T [TILESIZE - prepad*4 + postpad*4];
    myCacheLine = &myAllocCacheLine[-prepad];

    myPreX = prex;
    myPostX = postx;

    myForceCopy = false;
    if (myPreX || myPostX)
        myForceCopy = true;

    myDirty = false;

    myArray = vox;
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
bool
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::setIndex(int x, int y, int z)
{
    // Check if we have to reload our cache.
    if (myCurLine && y == myY && z == myZ)
    {
        if (x < myMaxValidX)
        {
            if (x == myX+1)
            {
                // A simple advanceX sufficies
                advanceX();
                return false;
            }

            if (x >= myMinValidX)
            {
                // We can just recenter our search location.
                resetX(x);
                // Other functions can't just do advanceX as that
                // just does ++ in X.
                return true;
            }
        }
    }

    // Store existing cache...
    if (DoWrite)
    {
        if (!TestForWrites || myDirty)
            writeCacheLine();
    }

    // Everything failed, return to reloading the cache.
    reloadCache(x, y, z);

    if (TestForWrites)
        myDirty = false;

    return true;
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
void
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::reloadCache(int x, int y, int z)
{
    UT_VoxelTile<T>     *tile;
    bool                 xout = false, yout = false, zout = false;
    bool                 manualbuild = false;

    myX = x;
    myY = y;
    myZ = z;
    myMinValidX = x & ~TILEMASK;
    myMaxValidX = myMinValidX + TILESIZE;
    if (myMaxValidX > myArray->getXRes())
        myMaxValidX = myArray->getXRes(); 

    if (x < 0 || x >= myArray->getXRes())
        xout = true;
    if (y < 0 || y >= myArray->getYRes())
        yout = true;
    if (z < 0 || z >= myArray->getZRes())
        zout = true;

    // If y or z are invalid, they will be invalid for every voxel
    if (yout || zout)
    {
        // We can often handle this by clamping...
        switch (myArray->getBorder())
        {
            case UT_VOXELBORDER_CONSTANT:
                buildConstantCache(myArray->getBorderValue());

                // ALL DONE
                return;

            case UT_VOXELBORDER_REPEAT:
                // Simply modulate our lookup.
                if (yout)
                {
                    y %= myArray->getYRes();
                    if (y < 0)
                        y += myArray->getYRes();
                }
                if (zout)
                {
                    z %= myArray->getZRes();
                    if (z < 0)
                        z += myArray->getZRes();
                }
                break;

            case UT_VOXELBORDER_STREAK:
            {
                // Clamp
                int             tx = 0;
                myArray->clampIndex(tx, y, z);
                break;
            }

            case UT_VOXELBORDER_EXTRAP:
            {
                // Force a manual build.
                manualbuild = true;
                break;
            }
        }
    }

    // Note y and z may no longer equal myY and myZ, this is not
    // a problem however as we have set up our cached versions
    // to the unclamped versions.

    if (xout || manualbuild)
    {
        // If there is no pre & post, this will always be constant.
        if (!myPreX && !myPostX && !manualbuild)
        {
            buildConstantCache(myArray->getValue(x, y, z));

            // ALL DONE
            return;
        }
        else
        {
            // If we are STREAK or CONSTANT, we have a constant
            // value in our own cache.
            // If we are REPEAT, we want to modulo our x value
            // and run the normal code path.
            if (myArray->getBorder() == UT_VOXELBORDER_REPEAT)
            {
                x %= myArray->getXRes();
                if (x < 0)
                    x += myArray->getXRes();
            }
            else
            {
                int             i;

                for (i = myPreX; i < 0; i++)
                {
                    myCacheLine[i] = myArray->getValue(myMinValidX+i, y, z);
                }

                T               value = myArray->getValue(x, y, z);
                for (; i < TILESIZE; i++)
                {
                    myCacheLine[i] = value;
                }

                for (; i < TILESIZE + myPostX; i++)
                {
                    myCacheLine[i] = myArray->getValue(myMaxValidX+i, y, z);
                }

                myCurLine = &myCacheLine[x & TILEMASK];
                myStride = 1;

                // ALL DONE
                return;
            }
        }
    }

    int                          xtile, ytile, ztile, tileidx;
    int                          lx, ly, lz;
    int                          i;

    xtile = x >> TILEBITS;
    ytile = y >> TILEBITS;
    ztile = z >> TILEBITS;

    // Get our local indices
    lx = x & TILEMASK;
    ly = y & TILEMASK;
    lz = z & TILEMASK;

    tileidx = (ztile * myArray->getTileRes(1) + ytile) * myArray->getTileRes(0);

    if (myPreX)
    {
        if (xtile)
        {
            // Simple to fetch...
            tile = myArray->getLinearTile(tileidx+xtile-1);
            for (i = myPreX; i < 0; i++)
            {
                // Safe to to & TILEMASK as we know earlier tiles
                // are always full...
                myCacheLine[i] = (*tile)(i & TILEMASK, ly, lz);
            }
        }
        else
        {
            if (myArray->getBorder() == UT_VOXELBORDER_REPEAT)
            {
                int             resx = myArray->getXRes();
                int             xpos;

                xpos = myPreX;
                xpos %= resx;
                // We add resx to guarantee we are in range.
                xpos += resx;

                // Manually invoke getValue()
                for (i = myPreX; i < 0; i++)
                {
                    myCacheLine[i] = (*myArray)(xpos, ly, lz);
                    xpos++;
                    if (xpos > resx)
                        xpos -= resx;
                }
            }
            else
            {
                T               value;

                if  (myArray->getBorder() == UT_VOXELBORDER_STREAK)
                {
                    tile = myArray->getLinearTile(tileidx+xtile);
                    value = (*tile)(0, ly, lz);
                }
                else
                    value = myArray->getBorderValue();

                // Fill in value.
                for (i = myPreX; i < 0; i++)
                    myCacheLine[i] = value;
            }
        }
    }

    if (myPostX)
    {
        int              cachelen;
        int              resx = myArray->getXRes();

        // This is where we start writing to.
        cachelen = myMaxValidX - myMinValidX;

        // Append our end part in.
        // First, determine if we read past the end...
        if (myMaxValidX + myPostX > myArray->getXRes())
        {
            // This can be very messy.  For example, we may have
            // a 1 wide tile after this tile and be looking two voxels
            // ahead, which means we can't guarantee our lookup
            // is entirely within one tile.
            // However, we can break it into two loops.
            int         xpos = myMaxValidX;

            // Portion that still fits in the next tile...
            i = 0;
            if (xpos < resx)
            {
                tile = myArray->getLinearTile(tileidx+xtile+1);
                for (; i < myPostX && xpos < resx; i++)
                {
                    myCacheLine[i + cachelen] = (*tile)(i, ly, lz);
                    xpos++;
                }
            }
            // Portion that reads past the end.
            if (i < myPostX)
            {
                if (myArray->getBorder() == UT_VOXELBORDER_REPEAT)
                {
                    xpos = xpos % resx;

                    // Revert to the array operator.
                    for (; i < myPostX; i++)
                    {
                        myCacheLine[i + cachelen] = (*myArray)(xpos, y, z);
                        xpos++;
                        if (xpos > resx)
                            xpos -= resx;
                    }
                }
                else
                {
                    T           value;

                    if (myArray->getBorder() == UT_VOXELBORDER_STREAK)
                    {
                        tile = myArray->getLinearTile(tileidx+xtile);
                        value = (*tile)(tile->xres()-1, ly, lz);
                    }
                    else
                        value = myArray->getBorderValue();

                    for (; i < myPostX; i++)
                        myCacheLine[i + cachelen] = value;
                }
            }
        }
        else
        {
            // All groovy, we fit in so thus must fit in the next tile
            tile = myArray->getLinearTile(tileidx+xtile+1);
            for (i = 0; i < myPostX; i++)
            {
                // Safe to to & TILEMASK as we know earlier tiles
                // are always full...
                myCacheLine[i + cachelen] = (*tile)(i, ly, lz);
            }
        }

    }

    tile = myArray->getLinearTile(tileidx+xtile);
    myCurLine = tile->fillCacheLine(myCacheLine, myStride, lx, ly, lz, myForceCopy, DoWrite);
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
void
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::buildConstantCache(T value)
{
    if (DoWrite)
    {
        // Force a full copy.
        myStride = 1;

        int             i;

        for (i = myPreX; i < TILESIZE+myPostX; i++)
            myCacheLine[i] = value;
        myCurLine = myCacheLine;
    }
    else
    {
        myCacheLine[0] = value;
        // These are to ensure our SSE is compatible.
        myCacheLine[1] = value;
        myCacheLine[2] = value;
        myCacheLine[3] = value;

        myCurLine = myCacheLine;
        myStride = 0;
    }
}

template <typename T, bool DoRead, bool DoWrite, bool TestForWrites>
void
UT_VoxelProbe<T, DoRead, DoWrite, TestForWrites>::writeCacheLine()
{
    if (!DoWrite)
    {
        UT_ASSERT(0);
        return;
    }
    // Ensure we have a valid loaded line, otherwise no point
    // doing a write back...
    if (!myCurLine)
        return;

    // Reset our current line...
    myCurLine -= myX - myMinValidX;

    // Determine if we actually have to write back,
    // if we had a pointer inside the tile we don't have to.
    if (myCurLine != myCacheLine)
        return;

    // Look up our voxel
    int                 xtile, ytile, ztile, y, z;
    UT_VoxelTile<T>     *tile;

    xtile = myMinValidX >> TILEBITS;
    ytile = myY >> TILEBITS;
    ztile = myZ >> TILEBITS;
    y = myY & TILEMASK;
    z = myZ & TILEMASK;

    tile = myArray->getTile(xtile, ytile, ztile);

    // Write back our results
    tile->writeCacheLine(myCurLine, y, z);
}

///
/// VoxelProbeCube functions
///
template <typename T>
UT_VoxelProbeCube<T>::UT_VoxelProbeCube()
{
    myValid = false;
}

template <typename T>
UT_VoxelProbeCube<T>::~UT_VoxelProbeCube()
{
}

template <typename T>
void
UT_VoxelProbeCube<T>::setCubeArray(UT_VoxelArray<T> *vox)
{
    myLines[0][0].setArray(vox, -1, 1);
    myLines[0][1].setArray(vox, -1, 1);
    myLines[0][2].setArray(vox, -1, 1);

    myLines[1][0].setArray(vox, -1, 1);
    myLines[1][1].setArray(vox, -1, 1);
    myLines[1][2].setArray(vox, -1, 1);

    myLines[2][0].setArray(vox, -1, 1);
    myLines[2][1].setArray(vox, -1, 1);
    myLines[2][2].setArray(vox, -1, 1);

    myValid = false;
}

template <typename T>
void
UT_VoxelProbeCube<T>::setPlusArray(UT_VoxelArray<T> *vox)
{
    /// This coudl be 0,0, but by keeping it the full range
    /// we ensure it is legal to rotate when we do a +1
    myLines[0][1].setArray(vox, -1, 1);

    myLines[1][0].setArray(vox, 0, 0);
    myLines[1][1].setArray(vox, -1, 1);
    myLines[1][2].setArray(vox, 0, 0);

    myLines[2][1].setArray(vox, -1, 1);

    myValid = false;
}

template <typename T>
bool
UT_VoxelProbeCube<T>::setIndexCube(int x, int y, int z)
{
    if (myValid && myZ == z)
    {
        if (myY == y)
        {
            // Potential for a simple advance...
            if (x < myMaxValidX && x == myX+1)
            {
                // AdvanceX.
                myLines[0][0].advanceX();
                myLines[0][1].advanceX();
                myLines[0][2].advanceX();

                myLines[1][0].advanceX();
                myLines[1][1].advanceX();
                myLines[1][2].advanceX();

                myLines[2][0].advanceX();
                myLines[2][1].advanceX();
                myLines[2][2].advanceX();

                // Update our cache.
                myX = x;

                return false;
            }
        }
#if 1
        else if (y == myY+1 && x < myMaxValidX && x >= myMinValidX)
        {
            // We have finished our x pass and just incremented y by one
            // Rather than resetting all our lines we can just swap
            // our y+1 lines into our current lines and then run the
            // normal reset.  
            rotateLines(myLines[0][0], myLines[1][0], myLines[2][0]);
            rotateLines(myLines[0][1], myLines[1][1], myLines[2][1]);
            rotateLines(myLines[0][2], myLines[1][2], myLines[2][2]);

            // The first 6 lines can just reset their X values
            // directly
            myLines[0][0].resetX(x);
            myLines[0][1].resetX(x);
            myLines[0][2].resetX(x);

            myLines[1][0].resetX(x);
            myLines[1][1].resetX(x);
            myLines[1][2].resetX(x);

            // Only the new lines need a reload.
            myLines[2][0].setIndex(x, y+1, z-1);
            myLines[2][1].setIndex(x, y+1, z);
            myLines[2][2].setIndex(x, y+1, z+1);

            // Update the cache values that have changed.
            myX = x;
            myY = y;

            return true;
        }
#endif
    }

    // Now just invoke setIndex on all children
    myLines[0][0].setIndex(x, y-1, z-1);
    myLines[0][1].setIndex(x, y-1, z);
    myLines[0][2].setIndex(x, y-1, z+1);

    myLines[1][0].setIndex(x, y, z-1);
    myLines[1][1].setIndex(x, y, z);
    myLines[1][2].setIndex(x, y, z+1);

    myLines[2][0].setIndex(x, y+1, z-1);
    myLines[2][1].setIndex(x, y+1, z);
    myLines[2][2].setIndex(x, y+1, z+1);

    // update our cache values
    myX = x;
    myY = y;
    myZ = z;
    myValid = true;
    myMinValidX = myLines[1][1].myMinValidX;
    myMaxValidX = myLines[1][1].myMaxValidX;

    return true;
}

template <typename T>
bool
UT_VoxelProbeCube<T>::setIndexPlus(int x, int y, int z)
{
    if (myValid && myZ == z)
    {
        if (myY == y)
        {
            // Potential for a simple advance...
            if (x < myMaxValidX && x == myX+1)
            {
                // AdvanceX.
                myLines[0][1].advanceX();

                myLines[1][0].advanceX();
                myLines[1][1].advanceX();
                myLines[1][2].advanceX();

                myLines[2][1].advanceX();

                // Update our cache.
                myX = x;

                return false;
            }
        }
        else if (y == myY+1 && x < myMaxValidX && x >= myMinValidX)
        {
            // We have finished our x pass and just incremented y by one
            // We can thus rotate the meaning of our central
            // cache lines and just reset their x pointers, leaving
            // only three real resets to be done.
            rotateLines(myLines[0][1], myLines[1][1], myLines[2][1]);

            myLines[0][1].resetX(x);
            myLines[1][1].resetX(x);

            myLines[1][0].setIndex(x, y, z-1);
            myLines[1][2].setIndex(x, y, z+1);

            myLines[2][1].setIndex(x, y+1, z);

            myX = x;
            myY = y;
            return true;
        }
    }

    // Now just invoke setIndex on all children
    myLines[0][1].setIndex(x, y-1, z);

    myLines[1][0].setIndex(x, y, z-1);
    myLines[1][1].setIndex(x, y, z);
    myLines[1][2].setIndex(x, y, z+1);

    myLines[2][1].setIndex(x, y+1, z);

    // update our cache values
    myX = x;
    myY = y;
    myZ = z;
    myValid = true;
    myMinValidX = myLines[1][1].myMinValidX;
    myMaxValidX = myLines[1][1].myMaxValidX;

    return true;
}

template <typename T>
fpreal64
UT_VoxelProbeCube<T>::curvature(const UT_Vector3 &invvoxelsize) const
{
    // These are our derivatives of Phi.
    fpreal64            Px, Py, Pz;
    fpreal64            Pxx, Pyy, Pzz;
    fpreal64            Pxy, Pxz, Pyz;
    fpreal64            gradlen;
    fpreal64            k;

    // Compute first derivatives.
    // dPhi = (Phi+1 - Phi-1) / 2 * dx
    
    Px = getValue(1, 0, 0) - getValue(-1, 0, 0);
    Px *= 0.5 * invvoxelsize.x();
    
    Py = getValue(0, 1, 0) - getValue(0, -1, 0);
    Py *= 0.5 * invvoxelsize.y();

    Pz = getValue(0, 0, 1) - getValue(0, 0, -1);
    Pz *= 0.5 * invvoxelsize.z();

    // Compute second derivatives.  (Note Pxy == Pyx)

    // d^2Phi = (Phi+1 - 2 Phi + Phi-1) / (dx*dx)
    Pxx = getValue(1, 0, 0)
          - 2 * getValue(0, 0, 0) 
          + getValue(-1, 0, 0);
    Pxx *= invvoxelsize.x() * invvoxelsize.x();

    Pyy = getValue(0, 1, 0)
          - 2 * getValue(0, 0, 0) 
          + getValue(0, -1, 0);
    Pyy *= invvoxelsize.y() * invvoxelsize.y();

    Pzz = getValue(0, 0, 1)
          - 2 * getValue(0, 0, 0) 
          + getValue(0, 0, -1);
    Pzz *= invvoxelsize.z() * invvoxelsize.z();

    // A bit more complicated :>
    Pxy  = getValue(1, 1,0) - getValue(-1, 1,0);
    Pxy -= getValue(1,-1,0) - getValue(-1,-1,0);
    Pxy *= 0.25 * invvoxelsize.x() * invvoxelsize.y();

    Pxz  = getValue(1,0, 1) - getValue(-1,0, 1);
    Pxz -= getValue(1,0,-1) - getValue(-1,0,-1);
    Pxz *= 0.25 * invvoxelsize.x() * invvoxelsize.z();

    Pyz  = getValue(0,1, 1) - getValue(0,-1, 1);
    Pyz -= getValue(0,1,-1) - getValue(0,-1,-1);
    Pyz *= 0.25 * invvoxelsize.y() * invvoxelsize.z();

    // Calculate the |grad(phi)| term;
    gradlen = SYSsqrt(Px * Px + Py * Py + Pz * Pz);

    // Finally, our curvature!
    // This is equation 1.8 from the holy book.
    // The problem is that it implies that 0 gradient means 0 curvature.
    // This is not true!
    // Even if Px,Py,Pz == 0, if Pxx != 0, we have a curved surface
    // consider a point at the maxima of a circle.
    k = Px*Px * (Pyy + Pzz) + Py*Py * (Pxx + Pzz) + Pz*Pz * (Pxx + Pyy);
    k -= 2 * (Pxy*Px*Py + Pyz*Py*Pz + Pxz*Px*Pz);

    // Avoid #IND in places with exactly zero curvature.
    if (!gradlen)
        k = 0;
    else
        k /= gradlen * gradlen * gradlen;

    // Clamp our curvature...
    fpreal64            maxk;

    maxk = invvoxelsize.maxComponent();
    if (k < -maxk)
        k = -maxk;
    if (k > maxk)
        k = maxk;

    return k;
}

template <typename T>
fpreal64
UT_VoxelProbeCube<T>::laplacian(const UT_Vector3 &invvoxelsize) const
{
    fpreal64            Pxx, Pyy, Pzz;
    fpreal64            centralval;

    centralval = getValue(0, 0, 0);

    // d^2Phi = (Phi+1 - 2 Phi + Phi-1) / (dx*dx)
    Pxx = getValue(1, 0, 0)
          - 2 * centralval
          + getValue(-1, 0, 0);
    Pxx *= invvoxelsize.x() * invvoxelsize.x();

    Pyy = getValue(0, 1, 0)
          - 2 * centralval
          + getValue(0, -1, 0);
    Pyy *= invvoxelsize.y() * invvoxelsize.y();

    Pzz = getValue(0, 0, +1)
          - 2 * centralval
          + getValue(0, 0, -1);
    Pzz *= invvoxelsize.z() * invvoxelsize.z();

    return Pxx + Pyy + Pzz;
}

template <typename T>
void
UT_VoxelProbeCube<T>::rotateLines(UT_VoxelProbe<T, true, false, false> &ym, UT_VoxelProbe<T, true, false, false> &y0, UT_VoxelProbe<T, true, false, false> &yp)
{
    T                           *tmpcache, *tmpalloc;
    const T                     *tmpcur;

    // We take advantage of the fact we know only a small portion
    // of the cache lines needs to be copied.
    tmpcache = ym.myCacheLine;
    tmpalloc = ym.myAllocCacheLine;
    tmpcur = ym.myCurLine;

    ym.myCacheLine = y0.myCacheLine;
    ym.myAllocCacheLine = y0.myAllocCacheLine;
    ym.myCurLine = y0.myCurLine;
    ym.myStride = y0.myStride;
    ym.myY++;

    y0.myCacheLine = yp.myCacheLine;
    y0.myAllocCacheLine = yp.myAllocCacheLine;
    y0.myCurLine = yp.myCurLine;
    y0.myStride = yp.myStride;
    y0.myY++;

    yp.myCacheLine = tmpcache;
    yp.myAllocCacheLine = tmpalloc;
    // Setting to zero will force a rebuild.
    yp.myCurLine = 0;
}

///
/// UT_VoxelProbeFace methods
///
template <typename T>
UT_VoxelProbeFace<T>::UT_VoxelProbeFace()
{
    myValid = false;
}

template <typename T>
UT_VoxelProbeFace<T>::~UT_VoxelProbeFace()
{
}


template <typename T>
void
UT_VoxelProbeFace<T>::setArray(UT_VoxelArray<T> *vx, UT_VoxelArray<T> *vy, UT_VoxelArray<T> *vz)
{
    // We need one extra to the right on the X probe
    myLines[0][0].setArray(vx, 0, 1);
    
    // The rest can be direct reads
    myLines[1][0].setArray(vy, 0, 0);
    myLines[1][1].setArray(vy, 0, 0);

    myLines[2][0].setArray(vz, 0, 0);
    myLines[2][1].setArray(vz, 0, 0);
    
    myValid = false;
}

template <typename T>
void
UT_VoxelProbeFace<T>::setVoxelSize(const UT_Vector3 &size)
{
    myVoxelSize = size;
    myInvVoxelSize = 1;
    myInvVoxelSize /= myVoxelSize;
}

template <typename T>
bool
UT_VoxelProbeFace<T>::setIndex(int x, int y, int z)
{
    if (myValid && myZ == z)
    {
        if (myY == y)
        {
            // Potential for a simple advance...
            if (x < myMaxValidX && x == myX+1)
            {
                // AdvanceX.
                myLines[0][0].advanceX();

                myLines[1][0].advanceX();
                myLines[1][1].advanceX();

                myLines[2][0].advanceX();
                myLines[2][1].advanceX();

                // Update our cache.
                myX = x;

                return false;
            }
        }
        else if (y == myY+1 && x < myMaxValidX && x >= myMinValidX)
        {
            // We have finished our x pass and just incremented y by one
            // We can swap our y lines to get to the next read for
            // those lines.
            swapLines(myLines[1][0], myLines[1][1]);

            myLines[1][0].resetX(x);

            // All the other lines need to be reloaded.
            myLines[0][0].setIndex(x, y, z);
            myLines[1][1].setIndex(x, y+1, z);

            myLines[2][0].setIndex(x, y, z);
            myLines[2][1].setIndex(x, y, z+1);

            myX = x;
            myY = y;
            return true;
        }
    }

    // Now just invoke setIndex on all children
    myLines[0][0].setIndex(x, y, z);

    myLines[1][0].setIndex(x, y, z);
    myLines[1][1].setIndex(x, y+1, z);

    myLines[2][0].setIndex(x, y, z);
    myLines[2][1].setIndex(x, y, z+1);

    // update our cache values
    myX = x;
    myY = y;
    myZ = z;
    myValid = true;
    myMinValidX = myLines[0][0].myMinValidX;
    myMaxValidX = myLines[0][0].myMaxValidX;

    return true;
}

template <typename T>
void
UT_VoxelProbeFace<T>::swapLines(UT_VoxelProbe<T, true, false, false> &ym, 
                                UT_VoxelProbe<T, true, false, false> &yp)
{
    T                           *tmpcache, *tmpalloc;
    const T                     *tmpcur;

    // We take advantage of the fact we know only a small portion
    // of the cache lines needs to be copied.
    tmpcache = ym.myCacheLine;
    tmpalloc = ym.myAllocCacheLine;
    tmpcur = ym.myCurLine;

    ym.myCacheLine = yp.myCacheLine;
    ym.myAllocCacheLine = yp.myAllocCacheLine;
    ym.myCurLine = yp.myCurLine;
    ym.myStride = yp.myStride;
    ym.myY++;

    yp.myCacheLine = tmpcache;
    yp.myAllocCacheLine = tmpalloc;
    // Setting to zero will force a rebuild.
    yp.myCurLine = 0;
}

///
/// VoxelProbeAverage methods
///
template <typename T, int XStep, int YStep, int ZStep>
void
UT_VoxelProbeAverage<T, XStep, YStep, ZStep>::setArray(UT_VoxelArray<T> *vox)
{
    int         prex = (XStep < 0) ? XStep : 0;
    int         postx = (XStep > 0) ? XStep : 0;

    myLines[0][0].setArray(vox, prex, postx);
    if (YStep)
    {
        myLines[1][0].setArray(vox, prex, postx);
        if (ZStep)
        {
            myLines[1][1].setArray(vox, prex, postx);
        }
    }
    if (ZStep)
        myLines[0][1].setArray(vox, prex, postx);
}

template <typename T, int XStep, int YStep, int ZStep>
bool
UT_VoxelProbeAverage<T, XStep, YStep, ZStep>::setIndex(int x, int y, int z)
{
    bool                result = false;

    // Adjust x, y, and z according to our half step.
    // y and z negative steps require us decrementing.  x steps
    // do not require a change as we use the pre/post to affect this,
    // if we adjusted the actual x we would get twice the cache misses.
    if (YStep < 0)
        y--;
    if (ZStep < 0)
        z--;

    result |= myLines[0][0].setIndex(x, y, z);
    if (YStep)
    {
        result |= myLines[1][0].setIndex(x, y+1, z);
        if (ZStep)
            result |= myLines[1][1].setIndex(x, y+1, z+1);
    }
    if (ZStep)
        result |= myLines[0][1].setIndex(x, y, z+1);

    return result;
}

---