CVolvoxels.cpp

#include "CVolVoxels.h"
#include "lifeVariables.h"


void CVolVoxels::Init( unsigned short * pVoxelData, unsigned short dimX, unsigned short dimY, unsigned short dimZ )
{
    if ( _pVoxelData || !pVoxelData )
        throw;

    _pVoxelData = pVoxelData;
    _dimX = dimX;
    _dimY = dimY;
    _dimZ = dimZ;
    _size = _dimX * _dimY * _dimZ;

    _boundingDim = (int)((float)sqrt( _dimX * _dimX + _dimY * _dimY + _dimZ * _dimZ ));

    if ( _grad )
        throw;

    _grad = new sGrad[_size];
    if ( !_grad)
        throw;

    for ( int i = 0; i < _size; i++ )
    {
        _grad[i].x = 0;
        _grad[i].y = 0;
        _grad[i].z = 0;
    }

    ComputeGradient();

}

unsigned short CVolVoxels::GetMaxValue()
{
    unsigned uNewVal = 0;

    if ( _MaxValue != 0 )
        return _MaxValue;

    for ( int i = 0; i < _size; i++ )
    {
        if ( uNewVal < _pVoxelData[i] )
        {
             uNewVal = _pVoxelData[i];
        }
    }

    _MaxValue = uNewVal;

    if ( ! _MaxValue )
        throw;

    return _MaxValue;
}
  
/*
 * Computes normalized gradient
 *
 *
 */
void CVolVoxels::ComputeGradient()
{
    float gradX, gradY, gradZ;

    for ( int k = 1; k < _dimZ-1; k++ )
        for ( int j = 1; j < _dimY-1; j++ )
            for ( int i = 1; i < _dimX-1; i++ )
            {

                gradX = (float) ( GVolVoxels.GetVoxelValue(i+1,j,k) - GVolVoxels.GetVoxelValue(i-1,j,k) ) / 2;
                gradY = (float) ( GVolVoxels.GetVoxelValue(i,j+1,k) - GVolVoxels.GetVoxelValue(i,j-1,k) ) / 2;
                gradZ = (float) ( GVolVoxels.GetVoxelValue(i,j,k+1) - GVolVoxels.GetVoxelValue(i,j,k-1) ) / 2;

                CVector3f  tmpGrad(gradX, gradY,gradZ);

                SetGradient(i,j,k, tmpGrad.GetNormalized());

                sGrad   outGrad;

                GetGradientTriLinear(i,j,k, &outGrad);
//                GetGradient(i,j,k,&outGrad);
            }
}

/*
 * GetDensityTriLiner returns the density at a given location - based on a tri-linear interpolation.
 * Note: Although the operations done on each voxel density are performed as a floating operation, each 
 *       density value is an unsigned short - and the resulting density is then returned as an unsigned short.
 *
 */
float CVolVoxels::GetDensityTriLinear(sPoint3f pos)
{
    //
    // define begin and end interpolation points
    //
    unsigned x0,y0,z0, x1,y1,z1;
    float    dx, dy, dz;

    x0 = pos.x;     y0 = pos.y;     z0 = pos.z;
    x1 = x0+1;  y1 = y0+1;  z1 = z0+1;

    //
    // dont do interpolation if we're at the end of the data set.
    //
    if ( x1 > _dimX || y1 > _dimY || z1 > _dimZ )
        return GVolVoxels.GetVoxelValue(x0,y0,z0);

    //
    // compute ratios
    //
    dx = pos.x - x0; dy = pos.y - y0; dz = pos.z - z0;

    //
    // interpolate along x-axis, then along y-axis, then along z-axis
    // 
    float dx0, dx1, dy0, dy1, dz0, dz1, dxyz;

    dx0 = GVolVoxels.GetVoxelValue(x0,y0,z0) * (1-dx) + GVolVoxels.GetVoxelValue(x1,y0,z0) * dx;
    dx1 = GVolVoxels.GetVoxelValue(x0,y1,z0) * (1-dx) + GVolVoxels.GetVoxelValue(x1,y1,z0) * dx;
    dy0 = GVolVoxels.GetVoxelValue(x0,y0,z1) * (1-dx) + GVolVoxels.GetVoxelValue(x1,y0,z1) * dx;
    dy1 = GVolVoxels.GetVoxelValue(x0,y1,z1) * (1-dx) + GVolVoxels.GetVoxelValue(x1,y1,z1) * dx;

    dz0 = dx0 * (1-dy) + dx1 * dy;
    dz1 = dy0 * (1-dy) + dy1 * dy;

    dxyz = dz0 * (1-dz) + dz1 * dz;

    return dxyz; // return interpolated value
}


void CVolVoxels::GetGradientTriLinear(float i, float j, float k, sGrad * pGrad)
{
    //
    // define begin and end interpolation points
    //
    unsigned x0,y0,z0, x1,y1,z1;
    float    dx, dy, dz;

    x0 = i;     y0 = j;     z0 = k;
    x1 = x0+1;  y1 = y0+1;  z1 = z0+1;

    //
    // dont do interpolation if we're at the end of the data set.
    //
    if ( x1 > _dimX || y1 > _dimY || z1 > _dimZ )
    {
        GVolVoxels.GetGradient(x0,y0,z0,pGrad);
        return;
    }

    //
    // compute ratios
    //
    dx = i - x0; dy = j - y0; dz = k - z0;

    //
    // interpolate along x-axis, then along y-axis, then along z-axis
    // 
    sGrad  grad_x0, grad_x1, grad_y0, grad_y1, grad_z0, grad_z1, grad_xyz;

    grad_x0 = GVolVoxels.GetGradient(x0,y0,z0) * (1-dx) + GVolVoxels.GetGradient(x1,y0,z0) * dx;
    grad_x1 = GVolVoxels.GetGradient(x0,y1,z0) * (1-dx) + GVolVoxels.GetGradient(x1,y1,z0) * dx;
    grad_y0 = GVolVoxels.GetGradient(x0,y0,z1) * (1-dx) + GVolVoxels.GetGradient(x1,y0,z1) * dx;
    grad_y1 = GVolVoxels.GetGradient(x0,y1,z1) * (1-dx) + GVolVoxels.GetGradient(x1,y1,z1) * dx;

    grad_z0 = grad_x0 * (1-dy) + grad_x1 * dy;
    grad_z1 = grad_y0 * (1-dy) + grad_y1 * dy;

    grad_xyz = grad_z0 * (1-dz) + grad_z1 * dz;

    pGrad->x = grad_xyz.x;
    pGrad->y = grad_xyz.y;
    pGrad->z = grad_xyz.z;

    return; // return interpolated value
}

---