volume.cxx

Go to the documentation of this file.

/*=========================================================================

  Program:   wxMaracas
  Module:    $RCSfile: volume.cxx,v $
  Language:  C++
  Date:      $Date: 2009/05/14 13:54:43 $
  Version:   $Revision: 1.5 $

  Copyright: (c) 2002, 2003
  License:
  
     This software is distributed WITHOUT ANY WARRANTY; without even 
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

// PS -> //#include <gsl/gsl_math.h>//PS
#include <memory.h>
#include "volume.hxx"

#ifdef KGFO_USE_VTK

#include <vtkCharArray.h>
#include <vtkDataArray.h>
#include <vtkDoubleArray.h>
#include <vtkFloatArray.h>
#include <vtkIntArray.h>
#include <vtkPointData.h>
#include <vtkShortArray.h>
#include <vtkUnsignedCharArray.h>
#include <vtkUnsignedIntArray.h>
#include <vtkUnsignedShortArray.h>

// -------------------------------------------------------------------------
const vtkIdType kVolume::VTKTypes[] = { VTK_CHAR, VTK_FLOAT, VTK_DOUBLE,
                                        VTK_INT, VTK_SHORT, VTK_UNSIGNED_CHAR,
                                        VTK_UNSIGNED_INT, VTK_UNSIGNED_SHORT };

#endif // KGFO_USE_VTK

//#ifdef KGFO_USE_IDO

// -------------------------------------------------------------------------
//const int kVolume::IDOTypes[] = { VOL_CHAR, VOL_FLOAT, VOL_DOUBLE, VOL_LONG,
//                                VOL_SHORT, VOL_UCHAR, VOL_ULONG,
//                                VOL_USHORT };

//#endif // KGFO_USE_IDO

// -------------------------------------------------------------------------
const void* kVolume::BLANK   = ( void* ) 0;
const void* kVolume::NOALLOC = ( void* )-1;
const int kVolume::SIZETypes[] = { sizeof( char ), sizeof( float ), 
                                   sizeof( double ), sizeof( int ),
                                   sizeof( short ), sizeof( uchar ),
                                   sizeof( uint ), sizeof( ushort ) };

// ---------------------------------------------------------------------------
template< class FROM, class TO >
    static void convertCastT( FROM* src, TO* dest, ulong size )
{
    FROM* end = src + size;
    for( ; src < end; src++, dest++ ) *dest = ( TO )*src;
}


// ---------------------------------------------------------------------------
template< class FROM, class TO >
    static void convertScaleT( FROM *src, TO *dest, ulong size,
                               double smin, double tmin, double slope )
{
    FROM* end = src + size;
    for( ; src < end; src++, dest++ )
        *dest = ( TO )( ( double( *src ) - smin ) * slope + tmin );
}

// ---------------------------------------------------------------------------
template< class TYPE >
    static void getMinMaxT( TYPE *src, ulong size,
                            double& min, double& max )
{
    TYPE* end = src + size;
    TYPE m, M;
    bool st;
        
    m = ( TYPE )0;
    M = ( TYPE )0;
    for( st = true; src < end; src++, st = false ) {
                
        if( *src < m || st ) m = *src;
        if( *src > M || st ) M = *src;
                
    } // rof
        
    min = ( double )m;
    max = ( double )M;
}

// -------------------------------------------------------------------------
kVolume::kVolume( )
    : _type( UCHAR ), _creator( SELF ),
#ifdef KGFO_USE_VTK
      _vtk( NULL ),
#endif // KGFO_USE_VTK
//#ifdef KGFO_USE_IDO
//      _privateIdo( NULL ),
//#endif // KGFO_USE_IDO
      _raw( NULL ), _columns( NULL ),
      _images( NULL )
{
    _dims[ CX ] = 1; _dims[ CY ] = 1; _dims[ CZ ] = 1;
    _sizes[ CX ] = 1; _sizes[ CY ] = 1; _sizes[ CZ ] = 1;
    allocate( );
    buildIndex( );
}

// -------------------------------------------------------------------------
kVolume::kVolume( Type type,
                  uint xdim, uint ydim, uint zdim,
                  double xsize, double ysize, double zsize,
                  void* data )
    : _type( type ), _creator( SELF ),
#ifdef KGFO_USE_VTK
      _vtk( NULL ),
#endif // KGFO_USE_VTK
//#ifdef KGFO_USE_IDO
//      _privateIdo( NULL ),
//#endif // KGFO_USE_IDO
      _raw( NULL ), _columns( NULL ),
      _images( NULL )
{
    _dims[ CX ] = xdim; _dims[ CY ] = ydim; _dims[ CZ ] = zdim;
    _sizes[ CX ] = xsize; _sizes[ CY ] = ysize; _sizes[ CZ ] = zsize;
    if( data != NOALLOC ) {
        
        if( data != BLANK )
            _raw = data;
        else
            allocate( );
        buildIndex( );
                        
    } // fi
}

// -------------------------------------------------------------------------
kVolume::kVolume( Type type,
                  const uint *dims,
                  const double *sizes,
                  void* data )
    : _type( type ), _creator( SELF ),
#ifdef KGFO_USE_VTK
      _vtk( NULL ),
#endif // KGFO_USE_VTK
//#ifdef KGFO_USE_IDO
//      _privateIdo( NULL ),
//#endif // KGFO_USE_IDO
      _raw( NULL ), _columns( NULL ),
      _images( NULL )
{
    memcpy( _dims, dims, 3 * sizeof( uint ) );
    memcpy( _sizes, sizes, 3 * sizeof( double ) );
    if( data != NOALLOC ) {
        
        if( data != BLANK )
            _raw = data;
        else
            allocate( );
        buildIndex( );
                        
    } // fi
}

// -------------------------------------------------------------------------
kVolume::kVolume( const kVolume& org )
    : _type( UCHAR ), _creator( SELF ),
#ifdef KGFO_USE_VTK
      _vtk( NULL ),
#endif // KGFO_USE_VTK
//#ifdef KGFO_USE_IDO
//      _privateIdo( NULL ),
//#endif // KGFO_USE_IDO
      _raw( NULL ), _columns( NULL ),
      _images( NULL )
{
    copyFrom( org );
}

// -------------------------------------------------------------------------
kVolume& kVolume::operator=( const kVolume& org )
{
    copyFrom( org );
    return( *this );
}

// -------------------------------------------------------------------------
void kVolume::copyFrom( const kVolume& org )
{
    _type = org._type;
    _creator = SELF;
    _raw = 0;
    _columns = 0;
    _images = 0;

    memcpy( _dims, org._dims, 3 * sizeof( uint ) );
    memcpy( _sizes, org._sizes, 3 * sizeof( double ) );
    if( org._raw ) {

        allocate( );
        buildIndex( );
        memcpy( _raw, org._raw, getRawSizeInBytes( ) );

    } // fi
}

// -------------------------------------------------------------------------
bool kVolume::sameDimension( const kVolume& comp )
{
    return( ( _dims[ CX ] == comp._dims[ CX ] &&
              _dims[ CY ] == comp._dims[ CY ] &&
              _dims[ CZ ] == comp._dims[ CZ ] ) );
}

// -------------------------------------------------------------------------
double kVolume::getPixel( uint x, uint y, uint z ) const
{
    double p;

    switch( _type ) {
                
    case CHAR:   p = ( double )(   ( char*** )_images )[ z ][ y ][ x ]; break;
    case FLOAT:  p = ( double )(  ( float*** )_images )[ z ][ y ][ x ]; break;
    case DOUBLE: p = ( double )( ( double*** )_images )[ z ][ y ][ x ]; break;
    case INT:    p = ( double )(    ( int*** )_images )[ z ][ y ][ x ]; break;
    case SHORT:  p = ( double )(  ( short*** )_images )[ z ][ y ][ x ]; break;
    case UCHAR:  p = ( double )(  ( uchar*** )_images )[ z ][ y ][ x ]; break;
    case UINT:   p = ( double )(   ( uint*** )_images )[ z ][ y ][ x ]; break;
    case USHORT: p = ( double )( ( ushort*** )_images )[ z ][ y ][ x ]; break;
    default: p = 0.0; break;

    } // fswitch

    return( p );
}

// -------------------------------------------------------------------------
void kVolume::setPixel( double v, uint x, uint y, uint z )
{

    switch( _type ) {
                
    case CHAR:     ( ( char*** )_images )[ z ][ y ][ x ] = ( char )v; break;
    case FLOAT:   ( ( float*** )_images )[ z ][ y ][ x ] = ( float )v; break;
    case DOUBLE: ( ( double*** )_images )[ z ][ y ][ x ] = ( double )v; break;
    case INT:       ( ( int*** )_images )[ z ][ y ][ x ] = ( int )v; break;
    case SHORT:   ( ( short*** )_images )[ z ][ y ][ x ] = ( short )v; break;
    case UCHAR:   ( ( uchar*** )_images )[ z ][ y ][ x ] = ( uchar )v; break;
    case UINT:     ( ( uint*** )_images )[ z ][ y ][ x ] = ( uint )v; break;
    case USHORT: ( ( ushort*** )_images )[ z ][ y ][ x ] = ( ushort )v; break;
    default: break;

    } // fswitch
}

// -------------------------------------------------------------------------
void kVolume::convertCast( Type type )
{
    if( _type != type ) {

        void* buffer;
        ulong size = getRawSize( );

        buffer = ( void* )new uchar[ size * SIZETypes[ type ] ];

        switch( _type ) {

        case CHAR:
            switch( type ) {

            case SHORT:  convertCastT( ( char* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( char* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( char* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( char* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( char* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( char* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( char* )_raw, ( uchar* )buffer, size ); break;
            default : break;

            } // fswitch
            break;

        case SHORT:
            switch( type ) {

            case CHAR:   convertCastT( ( short* )_raw, ( char* )buffer, size ); break;
            case INT:    convertCastT( ( short* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( short* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( short* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( short* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( short* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( short* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case INT:
            switch( type ) {

            case CHAR:   convertCastT( ( int* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( int* )_raw, ( short* )buffer, size ); break;
            case USHORT: convertCastT( ( int* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( int* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( int* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( int* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( int* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case USHORT:
            switch( type ) {

            case CHAR:   convertCastT( ( ushort* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( ushort* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( ushort* )_raw, ( int* )buffer, size ); break;
            case UINT:   convertCastT( ( ushort* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( ushort* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( ushort* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( ushort* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case UINT:
            switch( type ) {

            case CHAR:   convertCastT( ( uint* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( uint* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( uint* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( uint* )_raw, ( ushort* )buffer, size ); break;
            case FLOAT:  convertCastT( ( uint* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( uint* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( uint* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case FLOAT:
            switch( type ) {

            case CHAR:   convertCastT( ( float* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( float* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( float* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( float* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( float* )_raw, ( uint* )buffer, size ); break;
            case DOUBLE: convertCastT( ( float* )_raw, ( double* )buffer, size ); break;
            case UCHAR:  convertCastT( ( float* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case DOUBLE:
            switch( type ) {

            case CHAR:   convertCastT( ( double* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( double* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( double* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( double* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( double* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( double* )_raw, ( float* )buffer, size ); break;
            case UCHAR:  convertCastT( ( double* )_raw, ( uchar* )buffer, size ); break;
            default : break;
            } // fswitch
            break;

        case UCHAR:
            switch( type ) {

            case CHAR:   convertCastT( ( uchar* )_raw, ( char* )buffer, size ); break;
            case SHORT:  convertCastT( ( uchar* )_raw, ( short* )buffer, size ); break;
            case INT:    convertCastT( ( uchar* )_raw, ( int* )buffer, size ); break;
            case USHORT: convertCastT( ( uchar* )_raw, ( ushort* )buffer, size ); break;
            case UINT:   convertCastT( ( uchar* )_raw, ( uint* )buffer, size ); break;
            case FLOAT:  convertCastT( ( uchar* )_raw, ( float* )buffer, size ); break;
            case DOUBLE: convertCastT( ( uchar* )_raw, ( double* )buffer, size ); break;
            default : break;
            } // fswitch
            break;
     
        } // fswitch

        _type = type;
        deallocate( );
        _creator = SELF;
        _raw = buffer;
        buildIndex( );

    } // fi
}

// -------------------------------------------------------------------------
void kVolume::convertScale( Type type, double min, double max )
{
    double tmin, tmax, smin, smax;

// PS ->     //tmin = GSL_MIN( min, max ); //PS
        tmin= ((min<max)?min:max);
// PS ->     //tmax = GSL_MAX( min, max ); //PS
        tmax= ((min>max)?min:max);

    getMinMax( smin, smax );
        
    // compute scaling slope
    double a = ( tmax - tmin ) / ( smax - smin );

    void* buffer;
    ulong size = getRawSize( );

    buffer = ( void* )new uchar[ size * SIZETypes[ type ] ];

    switch( _type ) {

    case CHAR:
        switch( type ) {

        case CHAR:   convertScaleT( ( char* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( char* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( char* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( char* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( char* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( char* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( char* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( char* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case SHORT:
        switch( type ) {

        case CHAR:   convertScaleT( ( short* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( short* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( short* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( short* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( short* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( short* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( short* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( short* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case INT:
        switch( type ) {

        case CHAR:   convertScaleT( ( int* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( int* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( int* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( int* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( int* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( int* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( int* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( int* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case USHORT:
        switch( type ) {

        case CHAR:   convertScaleT( ( ushort* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( ushort* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( ushort* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( ushort* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( ushort* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( ushort* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( ushort* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( ushort* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case UINT:
        switch( type ) {

        case CHAR:   convertScaleT( ( uint* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( uint* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( uint* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( uint* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( uint* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( uint* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( uint* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( uint* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case UCHAR:
        switch( type ) {

        case CHAR:   convertScaleT( ( uchar* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( uchar* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( uchar* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( uchar* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( uchar* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( uchar* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( uchar* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( uchar* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case DOUBLE:
        switch( type ) {

        case CHAR:   convertScaleT( ( double* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( double* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( double* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( double* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( double* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( double* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( double* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( uchar* )_raw, ( double* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;
    case FLOAT:
        switch( type ) {

        case CHAR:   convertScaleT( ( float* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
        case SHORT:  convertScaleT( ( float* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
        case INT:    convertScaleT( ( float* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
        case USHORT: convertScaleT( ( float* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
        case UINT:   convertScaleT( ( float* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
        case FLOAT:  convertScaleT( ( float* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
        case DOUBLE: convertScaleT( ( float* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
        case UCHAR:  convertScaleT( ( float* )_raw, ( double* )buffer, size, smin, tmin, a ); break;

        } // fswitch
        break;

    } // fswitch

    _type = type;
    deallocate( );
    _creator = SELF;
    _raw = buffer;
    buildIndex( );
}

// -------------------------------------------------------------------------
void kVolume::getMinMax( double& min, double& max ) const
{
    ulong size = getRawSize( );

    switch( _type ) {
                
    case CHAR:   getMinMaxT( ( char* )_raw, size, min, max ); break;
    case FLOAT:  getMinMaxT( ( float* )_raw, size, min, max ); break;
    case DOUBLE: getMinMaxT( ( double* )_raw, size, min, max ); break;
    case INT:    getMinMaxT( ( int* )_raw, size, min, max ); break;
    case SHORT:  getMinMaxT( ( short* )_raw, size, min, max ); break;
    case UCHAR:  getMinMaxT( ( uchar* )_raw, size, min, max ); break;
    case UINT:   getMinMaxT( ( uint* )_raw, size, min, max ); break;
    case USHORT: getMinMaxT( ( ushort* )_raw, size, min, max ); break;

    } // fswitch
}

// -------------------------------------------------------------------------
double kVolume::getMin( ) const
{
    double m, M;

    getMinMax( m, M );
    return( m );
}

// -------------------------------------------------------------------------
double kVolume::getMax( ) const
{
    double m, M;

    getMinMax( m, M );
    return( M );
}

// -------------------------------------------------------------------------
/*double kVolume::GetMaxIntSphere( double* p, double r )
{
    int minX, minY, minZ, maxX, maxY, maxZ;
    gslobj_vector vP( p, 3 ), v( 3 );
    double maxint, tmp;
    bool start;
        
    minX = int( floor( p[ 0 ] - r ) );
    minY = int( floor( p[ 1 ] - r ) );
    minZ = int( floor( p[ 2 ] - r ) );
    maxX = int( ceil( p[ 0 ] + r ) );
    maxY = int( ceil( p[ 1 ] + r ) );
    maxZ = int( ceil( p[ 2 ] + r ) );
        
    minX = GSL_MAX( minX, 0 );
    minY = GSL_MAX( minY, 0 );
    minZ = GSL_MAX( minZ, 0 );
        
    maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );
    maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );
    maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );

    maxint = 0.0;
    start = true;
    for( v( 0 ) = minX; v( 0 ) <= maxX; v( 0 )++ )
      for( v( 1 ) = minY; v( 1 ) <= maxY; v( 1 )++ )
        for( v( 2 ) = minZ; v( 2 ) <= maxZ; v( 2 )++ )
          if( ( v - vP ).norm2( ) <= r ) {
            tmp = this->getPixel( ( uint )v(0), ( uint )v(1), ( uint )v(2));
            maxint = ( tmp > maxint || start )? tmp: maxint;
            start = false;
          } // fi

    return( maxint );
}*/
// -------------------------------------------------------------------------
unsigned short kVolume::GetMaxIntSphere2( double* p, double r )
{
    int minX, minY, minZ, maxX, maxY, maxZ;
    unsigned int v[3];
    unsigned short maxint = 0, tmp;

    minX = int( floor( p[ 0 ] - r ) );
    minY = int( floor( p[ 1 ] - r ) );
    minZ = int( floor( p[ 2 ] - r ) );
    maxX = int( ceil( p[ 0 ] + r ) );
    maxY = int( ceil( p[ 1 ] + r ) );
    maxZ = int( ceil( p[ 2 ] + r ) );

// PS ->     //minX = GSL_MAX( minX, 0 );//PS
// PS ->     //minY = GSL_MAX( minY, 0 );//PS
// PS ->     //minZ = GSL_MAX( minZ, 0 );//PS
        minX=((minX>0)?minX:0);
        minY=((minY>0)?minY:0);
        minZ=((minZ>0)?minZ:0);
        
// PS ->     //maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );//PS
// PS ->     //maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );//PS
// PS ->     //maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );//PS
        int xDim=this->getXdim( ) - 1;
        maxX= (maxX<xDim?maxX:xDim);

        int yDim=this->getYdim( ) - 1;
        maxY= (maxY<yDim?maxY:yDim);

        int zDim=this->getZdim( ) - 1;
        maxZ= (maxZ<zDim?maxZ:zDim);

    double r2 = r*r;  //need to do comparison in double ... don't know why...
    for( v[0] = minX; v[0] <= maxX; v[0]++ )
      for( v[1] = minY; v[1] <= maxY; v[1]++ )
        for( v[2] = minZ; v[2] <= maxZ; v[2]++ )
          if( ((v[0]-p[0])*(v[0]-p[0])+(v[1]-p[1])*(v[1]-p[1])+(v[2]-p[2])*(v[2]-p[2])) <= r2 )
          {
            tmp = (unsigned short)this->getPixel( v[0], v[1], v[2] );
            maxint = ( tmp > maxint ) ? tmp : maxint;
          }

    return( maxint );
}

// -------------------------------------------------------------------------
void kVolume::allocate( )
{
    ulong size = getRawSizeInBytes( );

    if( _creator == SELF ) {

        _raw = ( void* )new uchar[ size ];
        memset( _raw, 0, size );

    } // fi
}

// -------------------------------------------------------------------------
void kVolume::buildIndex( )
{
    ulong size;

    size = ( _dims[ CZ ] * sizeof( uchar** ) ) +
        ( _dims[ CZ ] * _dims[ CY ] * sizeof( void* ) );
                
        _images = ( void*** )new uchar[ size ];
                
        _columns = ( void** )( _images + _dims[ CZ ] );
        void** plane = _columns;
        for( uint z = 0; z < _dims[ CZ ]; z++ ) {
                        
            _images[ z ] = plane;
            plane += _dims[ CY ];
                        
        } // rof
        void* line = _raw;
        for( uint y = 0; y < _dims[ CZ ] * _dims[ CY ]; y++ ) {
                        
            _columns[ y ] = line;
            line = ( void* )( ( uchar* ) line +
                              _dims[ CX ] * SIZETypes[ _type ] );
                        
        } // rof
        
#ifdef KGFO_USE_VTK
        
    vtkCharArray* carray;
    vtkDoubleArray* darray;
    vtkFloatArray* farray;
    vtkIntArray* iarray;
    vtkShortArray* sarray;
    vtkUnsignedCharArray* ucarray;
    vtkUnsignedIntArray* uiarray;
    vtkUnsignedShortArray* usarray;
        
    size = _dims[ CX ] * _dims[ CY ] * _dims[ CZ ];

    if( _creator == SELF || _creator == IDO ) {
          
        _vtk = vtkImageData::New( );
        _vtk->SetDimensions( _dims[ CX ], _dims[ CY ], _dims[ CZ ] );
        _vtk->SetSpacing( _sizes[ CX ], _sizes[ CY ], _sizes[ CZ ] );

        switch( _type ) {

        case CHAR:
            carray = vtkCharArray::New( );
            carray->SetNumberOfComponents( 1 );
            carray->SetArray( ( char* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_CHAR );
            _vtk->GetPointData( )->SetScalars( carray );
            carray->Delete( );
            break;

        case UCHAR:
            ucarray = vtkUnsignedCharArray::New( );
            ucarray->SetNumberOfComponents( 1 );
            ucarray->SetArray( ( uchar* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_UNSIGNED_CHAR );
            _vtk->GetPointData( )->SetScalars( ucarray );
            ucarray->Delete( );
            break;

        case SHORT:
            sarray = vtkShortArray::New( );
            sarray->SetNumberOfComponents( 1 );
            sarray->SetArray( ( short* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_SHORT );
            _vtk->GetPointData( )->SetScalars( sarray );
            sarray->Delete( );
            break;

        case INT:
            iarray = vtkIntArray::New( );
            iarray->SetNumberOfComponents( 1 );
            iarray->SetArray( ( int* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_INT );
            _vtk->GetPointData( )->SetScalars( iarray );
            iarray->Delete( );
            break;

        case USHORT:
            usarray = vtkUnsignedShortArray::New( );
            usarray->SetNumberOfComponents( 1 );
            usarray->SetArray( ( ushort* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_UNSIGNED_SHORT );
            _vtk->GetPointData( )->SetScalars( usarray );
            usarray->Delete( );
            break;

        case UINT:
            uiarray = vtkUnsignedIntArray::New( );
            uiarray->SetNumberOfComponents( 1 );
            uiarray->SetArray( ( uint* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_UNSIGNED_INT );
            _vtk->GetPointData( )->SetScalars( uiarray );
            uiarray->Delete( );
            break;

        case FLOAT:
            farray = vtkFloatArray::New( );
            farray->SetNumberOfComponents( 1 );
            farray->SetArray( ( float* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_FLOAT );
            _vtk->GetPointData( )->SetScalars( farray );
            farray->Delete( );
            break;
                  
        case DOUBLE:
            darray = vtkDoubleArray::New( );
            darray->SetNumberOfComponents( 1 );
            darray->SetArray( ( double* )( _raw ), size, 1 );
            _vtk->SetScalarType( VTK_DOUBLE );
            _vtk->GetPointData( )->SetScalars( darray );
            darray->Delete( );
            break;
                  
        } // fswitch
          
    } // fi

#endif // KGFO_USE_VTK
}

// -------------------------------------------------------------------------
void kVolume::deallocate( )
{
#ifdef KGFO_USE_VTK

    if( _vtk ) _vtk->Delete();
    _vtk = NULL;

#endif // KGFO_USE_VTK

//#ifdef KGFO_USE_IDO

//    if( _creator == SELF || _creator == VTK ) {

//      delete[] _privateIdo;
//      _privateIdo = NULL;
        
//    } // fi

//#else

    delete[] ( uchar* )_images;

//#endif // KGFO_USE_IDO

    if( _raw && _creator == SELF )

//EED purify 12/sept/2006
//      delete[] ( uchar* )_raw;

        free ( _raw );

    _creator    = SELF;
    _raw                = NULL;
    _columns    = NULL;
    _images             = NULL;
}

#ifdef KGFO_USE_VTK

// -------------------------------------------------------------------------
kVolume::kVolume( vtkImageData* org )
    : _raw( 0 ), _columns( 0 ), _images( 0 ),
      _creator( VTK ),
//#ifdef KGFO_USE_IDO
//      _privateIdo( NULL ),
//#endif // KGFO_USE_IDO
      _vtk( 0 )
{
    //int i, j, k, y;
    int itmp[ 3 ];
    double ftmp[ 3 ];
    //double v;

    switch( org->GetScalarType( ) ) {

    case VTK_CHAR:           _type = CHAR;   break;
    case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
    case VTK_SHORT:          _type = SHORT;  break;
    case VTK_INT:            _type = INT;    break;
    case VTK_UNSIGNED_SHORT: _type = USHORT; break;
    case VTK_UNSIGNED_INT:   _type = UINT;   break;
    case VTK_FLOAT:          _type = FLOAT;  break;
    case VTK_DOUBLE:         _type = DOUBLE; break;
    default: break;
                
    } // fswitch
        
    org->GetDimensions( itmp );
    _dims[ CX ] = ( uint )itmp[ 0 ];
    _dims[ CY ] = ( uint )itmp[ 1 ];
    _dims[ CZ ] = ( uint )itmp[ 2 ];
    org->GetSpacing( ftmp );
    _sizes[ CX ] = ( double )ftmp[ 0 ];
    _sizes[ CY ] = ( double )ftmp[ 1 ];
    _sizes[ CZ ] = ( double )ftmp[ 2 ];

    _raw = org->GetPointData( )->GetScalars( )->GetVoidPointer( 0 );
    //_vtk = org;
    _vtk = vtkImageData::New();
    _vtk->ShallowCopy( org );
    buildIndex( );
}

// -------------------------------------------------------------------------
kVolume& kVolume::operator=( vtkImageData* org )
{
    copyFrom( org );
    return( *this );
}

// -------------------------------------------------------------------------
void kVolume::copyFrom( vtkImageData* org )
{
    int i, j, k;//, y;
    int itmp[ 3 ];
    int ext[ 6 ];
    double ftmp[ 3 ];
    double v;

    deallocate( );
        
//#ifdef KGFO_USE_IDO

//    _privateIdo = NULL;

//#endif // KGFO_USE_IDO

    _vtk                = NULL;
    _raw                = NULL;
    _columns    = NULL;
    _images             = NULL;
    _creator    = SELF;

    switch( org->GetScalarType( ) ) {

    case VTK_CHAR:           _type = CHAR;   break;
    case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
    case VTK_SHORT:          _type = SHORT;  break;
    case VTK_INT:            _type = INT;    break;
    case VTK_UNSIGNED_SHORT: _type = USHORT; break;
    case VTK_UNSIGNED_INT:   _type = UINT;   break;
    case VTK_FLOAT:          _type = FLOAT;  break;
    case VTK_DOUBLE:         _type = DOUBLE; break;
    default: break;
                
    } // fswitch
        
    org->GetDimensions( itmp );
    _dims[ CX ] = ( uint )itmp[ 0 ];
    _dims[ CY ] = ( uint )itmp[ 1 ];
    _dims[ CZ ] = ( uint )itmp[ 2 ];
    org->GetSpacing( ftmp );
    _sizes[ CX ] = ( double )ftmp[ 0 ];
    _sizes[ CY ] = ( double )ftmp[ 1 ];
    _sizes[ CZ ] = ( double )ftmp[ 2 ];

    allocate( );
    buildIndex( );

    // This avoids vtk extent crap conversion...
    org->GetExtent( ext );
    for( i = ext[ 0 ]; i <= ext[ 1 ]; i++ ) {
        for( j = ext[ 2 ]; j <= ext[ 3 ]; j++ ) {
            for( k = ext[ 4 ]; k <= ext[ 5 ]; k++ ) {           
                v = org->GetScalarComponentAsDouble( i, j, k, 0 );
                setPixel( v, i - ext[ 0 ], j - ext[ 2 ], k - ext[ 4 ] );                
            } // rof
        } // rof
    } // rof
}

#endif // KGFO_USE_VTK


// eof - volume.cxx

---