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

  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 )
{
   /**
    * unsigned short range : 0 -> 65535
    * unsigned int   range : 0 -> 2147483647 // 2^31 - 1
    */
    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