lib/maracasVisuLib/src/kernel/volume.cxx

   1 /*=========================================================================
   2
   3   Program:   wxMaracas
   4   Module:    $RCSfile: volume.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2008/11/21 18:13:38 $
   7   Version:   $Revision: 1.3 $
   8
   9   Copyright: (c) 2002, 2003
  10   License:
  11
  12      This software is distributed WITHOUT ANY WARRANTY; without even
  13      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14      PURPOSE.  See the above copyright notice for more information.
  15
  16 =========================================================================*/
  17
  18 // PS -> //#include <gsl/gsl_math.h>//PS
  19 #include <memory.h>
  20 #include "volume.hxx"
  21
  22 #ifdef KGFO_USE_VTK
  23
  24 #include <vtkCharArray.h>
  25 #include <vtkDataArray.h>
  26 #include <vtkDoubleArray.h>
  27 #include <vtkFloatArray.h>
  28 #include <vtkIntArray.h>
  29 #include <vtkPointData.h>
  30 #include <vtkShortArray.h>
  31 #include <vtkUnsignedCharArray.h>
  32 #include <vtkUnsignedIntArray.h>
  33 #include <vtkUnsignedShortArray.h>
  34
  35 // -------------------------------------------------------------------------
  36 const vtkIdType kVolume::VTKTypes[] = { VTK_CHAR, VTK_FLOAT, VTK_DOUBLE,
  37                                         VTK_INT, VTK_SHORT, VTK_UNSIGNED_CHAR,
  38                                         VTK_UNSIGNED_INT, VTK_UNSIGNED_SHORT };
  39
  40 #endif // KGFO_USE_VTK
  41
  42 //#ifdef KGFO_USE_IDO
  43
  44 // -------------------------------------------------------------------------
  45 //const int kVolume::IDOTypes[] = { VOL_CHAR, VOL_FLOAT, VOL_DOUBLE, VOL_LONG,
  46 //                                VOL_SHORT, VOL_UCHAR, VOL_ULONG,
  47 //                                VOL_USHORT };
  48
  49 //#endif // KGFO_USE_IDO
  50
  51 // -------------------------------------------------------------------------
  52 const void* kVolume::BLANK   = ( void* ) 0;
  53 const void* kVolume::NOALLOC = ( void* )-1;
  54 const int kVolume::SIZETypes[] = { sizeof( char ), sizeof( float ),
  55                                    sizeof( double ), sizeof( int ),
  56                                    sizeof( short ), sizeof( uchar ),
  57                                    sizeof( uint ), sizeof( ushort ) };
  58
  59 // ---------------------------------------------------------------------------
  60 template< class FROM, class TO >
  61     static void convertCastT( FROM* src, TO* dest, ulong size )
  62 {
  63     FROM* end = src + size;
  64     for( ; src < end; src++, dest++ ) *dest = ( TO )*src;
  65 }
  66
  67
  68 // ---------------------------------------------------------------------------
  69 template< class FROM, class TO >
  70     static void convertScaleT( FROM *src, TO *dest, ulong size,
  71                                double smin, double tmin, double slope )
  72 {
  73     FROM* end = src + size;
  74     for( ; src < end; src++, dest++ )
  75         *dest = ( TO )( ( double( *src ) - smin ) * slope + tmin );
  76 }
  77
  78 // ---------------------------------------------------------------------------
  79 template< class TYPE >
  80     static void getMinMaxT( TYPE *src, ulong size,
  81                             double& min, double& max )
  82 {
  83     TYPE* end = src + size;
  84     TYPE m, M;
  85     bool st;
  86
  87     m = ( TYPE )0;
  88     M = ( TYPE )0;
  89     for( st = true; src < end; src++, st = false ) {
  90
  91         if( *src < m || st ) m = *src;
  92         if( *src > M || st ) M = *src;
  93
  94     } // rof
  95
  96     min = ( double )m;
  97     max = ( double )M;
  98 }
  99
 100 // -------------------------------------------------------------------------
 101 kVolume::kVolume( )
 102     : _type( UCHAR ), _creator( SELF ),
 103 #ifdef KGFO_USE_VTK
 104       _vtk( NULL ),
 105 #endif // KGFO_USE_VTK
 106 //#ifdef KGFO_USE_IDO
 107 //      _privateIdo( NULL ),
 108 //#endif // KGFO_USE_IDO
 109       _raw( NULL ), _columns( NULL ),
 110       _images( NULL )
 111 {
 112     _dims[ CX ] = 1; _dims[ CY ] = 1; _dims[ CZ ] = 1;
 113     _sizes[ CX ] = 1; _sizes[ CY ] = 1; _sizes[ CZ ] = 1;
 114     allocate( );
 115     buildIndex( );
 116 }
 117
 118 // -------------------------------------------------------------------------
 119 kVolume::kVolume( Type type,
 120                   uint xdim, uint ydim, uint zdim,
 121                   double xsize, double ysize, double zsize,
 122                   void* data )
 123     : _type( type ), _creator( SELF ),
 124 #ifdef KGFO_USE_VTK
 125       _vtk( NULL ),
 126 #endif // KGFO_USE_VTK
 127 //#ifdef KGFO_USE_IDO
 128 //      _privateIdo( NULL ),
 129 //#endif // KGFO_USE_IDO
 130       _raw( NULL ), _columns( NULL ),
 131       _images( NULL )
 132 {
 133     _dims[ CX ] = xdim; _dims[ CY ] = ydim; _dims[ CZ ] = zdim;
 134     _sizes[ CX ] = xsize; _sizes[ CY ] = ysize; _sizes[ CZ ] = zsize;
 135     if( data != NOALLOC ) {
 136
 137         if( data != BLANK )
 138             _raw = data;
 139         else
 140             allocate( );
 141         buildIndex( );
 142
 143     } // fi
 144 }
 145
 146 // -------------------------------------------------------------------------
 147 kVolume::kVolume( Type type,
 148                   const uint *dims,
 149                   const double *sizes,
 150                   void* data )
 151     : _type( type ), _creator( SELF ),
 152 #ifdef KGFO_USE_VTK
 153       _vtk( NULL ),
 154 #endif // KGFO_USE_VTK
 155 //#ifdef KGFO_USE_IDO
 156 //      _privateIdo( NULL ),
 157 //#endif // KGFO_USE_IDO
 158       _raw( NULL ), _columns( NULL ),
 159       _images( NULL )
 160 {
 161     memcpy( _dims, dims, 3 * sizeof( uint ) );
 162     memcpy( _sizes, sizes, 3 * sizeof( double ) );
 163     if( data != NOALLOC ) {
 164
 165         if( data != BLANK )
 166             _raw = data;
 167         else
 168             allocate( );
 169         buildIndex( );
 170
 171     } // fi
 172 }
 173
 174 // -------------------------------------------------------------------------
 175 kVolume::kVolume( const kVolume& org )
 176     : _type( UCHAR ), _creator( SELF ),
 177 #ifdef KGFO_USE_VTK
 178       _vtk( NULL ),
 179 #endif // KGFO_USE_VTK
 180 //#ifdef KGFO_USE_IDO
 181 //      _privateIdo( NULL ),
 182 //#endif // KGFO_USE_IDO
 183       _raw( NULL ), _columns( NULL ),
 184       _images( NULL )
 185 {
 186     copyFrom( org );
 187 }
 188
 189 // -------------------------------------------------------------------------
 190 kVolume& kVolume::operator=( const kVolume& org )
 191 {
 192     copyFrom( org );
 193     return( *this );
 194 }
 195
 196 // -------------------------------------------------------------------------
 197 void kVolume::copyFrom( const kVolume& org )
 198 {
 199     _type = org._type;
 200     _creator = SELF;
 201     _raw = 0;
 202     _columns = 0;
 203     _images = 0;
 204
 205     memcpy( _dims, org._dims, 3 * sizeof( uint ) );
 206     memcpy( _sizes, org._sizes, 3 * sizeof( double ) );
 207     if( org._raw ) {
 208
 209         allocate( );
 210         buildIndex( );
 211         memcpy( _raw, org._raw, getRawSizeInBytes( ) );
 212
 213     } // fi
 214 }
 215
 216 // -------------------------------------------------------------------------
 217 bool kVolume::sameDimension( const kVolume& comp )
 218 {
 219     return( ( _dims[ CX ] == comp._dims[ CX ] &&
 220               _dims[ CY ] == comp._dims[ CY ] &&
 221               _dims[ CZ ] == comp._dims[ CZ ] ) );
 222 }
 223
 224 // -------------------------------------------------------------------------
 225 double kVolume::getPixel( uint x, uint y, uint z ) const
 226 {
 227     double p;
 228
 229     switch( _type ) {
 230
 231     case CHAR:   p = ( double )(   ( char*** )_images )[ z ][ y ][ x ]; break;
 232     case FLOAT:  p = ( double )(  ( float*** )_images )[ z ][ y ][ x ]; break;
 233     case DOUBLE: p = ( double )( ( double*** )_images )[ z ][ y ][ x ]; break;
 234     case INT:    p = ( double )(    ( int*** )_images )[ z ][ y ][ x ]; break;
 235     case SHORT:  p = ( double )(  ( short*** )_images )[ z ][ y ][ x ]; break;
 236     case UCHAR:  p = ( double )(  ( uchar*** )_images )[ z ][ y ][ x ]; break;
 237     case UINT:   p = ( double )(   ( uint*** )_images )[ z ][ y ][ x ]; break;
 238     case USHORT: p = ( double )( ( ushort*** )_images )[ z ][ y ][ x ]; break;
 239     default: p = 0.0; break;
 240
 241     } // fswitch
 242
 243     return( p );
 244 }
 245
 246 // -------------------------------------------------------------------------
 247 void kVolume::setPixel( double v, uint x, uint y, uint z )
 248 {
 249
 250     switch( _type ) {
 251
 252     case CHAR:     ( ( char*** )_images )[ z ][ y ][ x ] = ( char )v; break;
 253     case FLOAT:   ( ( float*** )_images )[ z ][ y ][ x ] = ( float )v; break;
 254     case DOUBLE: ( ( double*** )_images )[ z ][ y ][ x ] = ( double )v; break;
 255     case INT:       ( ( int*** )_images )[ z ][ y ][ x ] = ( int )v; break;
 256     case SHORT:   ( ( short*** )_images )[ z ][ y ][ x ] = ( short )v; break;
 257     case UCHAR:   ( ( uchar*** )_images )[ z ][ y ][ x ] = ( uchar )v; break;
 258     case UINT:     ( ( uint*** )_images )[ z ][ y ][ x ] = ( uint )v; break;
 259     case USHORT: ( ( ushort*** )_images )[ z ][ y ][ x ] = ( ushort )v; break;
 260     default: break;
 261
 262     } // fswitch
 263 }
 264
 265 // -------------------------------------------------------------------------
 266 void kVolume::convertCast( Type type )
 267 {
 268     if( _type != type ) {
 269
 270         void* buffer;
 271         ulong size = getRawSize( );
 272
 273         buffer = ( void* )new uchar[ size * SIZETypes[ type ] ];
 274
 275         switch( _type ) {
 276
 277         case CHAR:
 278             switch( type ) {
 279
 280             case SHORT:  convertCastT( ( char* )_raw, ( short* )buffer, size ); break;
 281             case INT:    convertCastT( ( char* )_raw, ( int* )buffer, size ); break;
 282             case USHORT: convertCastT( ( char* )_raw, ( ushort* )buffer, size ); break;
 283             case UINT:   convertCastT( ( char* )_raw, ( uint* )buffer, size ); break;
 284             case FLOAT:  convertCastT( ( char* )_raw, ( float* )buffer, size ); break;
 285             case DOUBLE: convertCastT( ( char* )_raw, ( double* )buffer, size ); break;
 286             case UCHAR:  convertCastT( ( char* )_raw, ( uchar* )buffer, size ); break;
 287
 288             } // fswitch
 289             break;
 290
 291         case SHORT:
 292             switch( type ) {
 293
 294             case CHAR:   convertCastT( ( short* )_raw, ( char* )buffer, size ); break;
 295             case INT:    convertCastT( ( short* )_raw, ( int* )buffer, size ); break;
 296             case USHORT: convertCastT( ( short* )_raw, ( ushort* )buffer, size ); break;
 297             case UINT:   convertCastT( ( short* )_raw, ( uint* )buffer, size ); break;
 298             case FLOAT:  convertCastT( ( short* )_raw, ( float* )buffer, size ); break;
 299             case DOUBLE: convertCastT( ( short* )_raw, ( double* )buffer, size ); break;
 300             case UCHAR:  convertCastT( ( short* )_raw, ( uchar* )buffer, size ); break;
 301
 302             } // fswitch
 303             break;
 304
 305         case INT:
 306             switch( type ) {
 307
 308             case CHAR:   convertCastT( ( int* )_raw, ( char* )buffer, size ); break;
 309             case SHORT:  convertCastT( ( int* )_raw, ( short* )buffer, size ); break;
 310             case USHORT: convertCastT( ( int* )_raw, ( ushort* )buffer, size ); break;
 311             case UINT:   convertCastT( ( int* )_raw, ( uint* )buffer, size ); break;
 312             case FLOAT:  convertCastT( ( int* )_raw, ( float* )buffer, size ); break;
 313             case DOUBLE: convertCastT( ( int* )_raw, ( double* )buffer, size ); break;
 314             case UCHAR:  convertCastT( ( int* )_raw, ( uchar* )buffer, size ); break;
 315
 316             } // fswitch
 317             break;
 318
 319         case USHORT:
 320             switch( type ) {
 321
 322             case CHAR:   convertCastT( ( ushort* )_raw, ( char* )buffer, size ); break;
 323             case SHORT:  convertCastT( ( ushort* )_raw, ( short* )buffer, size ); break;
 324             case INT:    convertCastT( ( ushort* )_raw, ( int* )buffer, size ); break;
 325             case UINT:   convertCastT( ( ushort* )_raw, ( uint* )buffer, size ); break;
 326             case FLOAT:  convertCastT( ( ushort* )_raw, ( float* )buffer, size ); break;
 327             case DOUBLE: convertCastT( ( ushort* )_raw, ( double* )buffer, size ); break;
 328             case UCHAR:  convertCastT( ( ushort* )_raw, ( uchar* )buffer, size ); break;
 329
 330             } // fswitch
 331             break;
 332
 333         case UINT:
 334             switch( type ) {
 335
 336             case CHAR:   convertCastT( ( uint* )_raw, ( char* )buffer, size ); break;
 337             case SHORT:  convertCastT( ( uint* )_raw, ( short* )buffer, size ); break;
 338             case INT:    convertCastT( ( uint* )_raw, ( int* )buffer, size ); break;
 339             case USHORT: convertCastT( ( uint* )_raw, ( ushort* )buffer, size ); break;
 340             case FLOAT:  convertCastT( ( uint* )_raw, ( float* )buffer, size ); break;
 341             case DOUBLE: convertCastT( ( uint* )_raw, ( double* )buffer, size ); break;
 342             case UCHAR:  convertCastT( ( uint* )_raw, ( uchar* )buffer, size ); break;
 343
 344             } // fswitch
 345             break;
 346
 347         case FLOAT:
 348             switch( type ) {
 349
 350             case CHAR:   convertCastT( ( float* )_raw, ( char* )buffer, size ); break;
 351             case SHORT:  convertCastT( ( float* )_raw, ( short* )buffer, size ); break;
 352             case INT:    convertCastT( ( float* )_raw, ( int* )buffer, size ); break;
 353             case USHORT: convertCastT( ( float* )_raw, ( ushort* )buffer, size ); break;
 354             case UINT:   convertCastT( ( float* )_raw, ( uint* )buffer, size ); break;
 355             case DOUBLE: convertCastT( ( float* )_raw, ( double* )buffer, size ); break;
 356             case UCHAR:  convertCastT( ( float* )_raw, ( uchar* )buffer, size ); break;
 357
 358             } // fswitch
 359             break;
 360
 361         case DOUBLE:
 362             switch( type ) {
 363
 364             case CHAR:   convertCastT( ( double* )_raw, ( char* )buffer, size ); break;
 365             case SHORT:  convertCastT( ( double* )_raw, ( short* )buffer, size ); break;
 366             case INT:    convertCastT( ( double* )_raw, ( int* )buffer, size ); break;
 367             case USHORT: convertCastT( ( double* )_raw, ( ushort* )buffer, size ); break;
 368             case UINT:   convertCastT( ( double* )_raw, ( uint* )buffer, size ); break;
 369             case FLOAT:  convertCastT( ( double* )_raw, ( float* )buffer, size ); break;
 370             case UCHAR:  convertCastT( ( double* )_raw, ( uchar* )buffer, size ); break;
 371
 372             } // fswitch
 373             break;
 374
 375         case UCHAR:
 376             switch( type ) {
 377
 378             case CHAR:   convertCastT( ( uchar* )_raw, ( char* )buffer, size ); break;
 379             case SHORT:  convertCastT( ( uchar* )_raw, ( short* )buffer, size ); break;
 380             case INT:    convertCastT( ( uchar* )_raw, ( int* )buffer, size ); break;
 381             case USHORT: convertCastT( ( uchar* )_raw, ( ushort* )buffer, size ); break;
 382             case UINT:   convertCastT( ( uchar* )_raw, ( uint* )buffer, size ); break;
 383             case FLOAT:  convertCastT( ( uchar* )_raw, ( float* )buffer, size ); break;
 384             case DOUBLE: convertCastT( ( uchar* )_raw, ( double* )buffer, size ); break;
 385
 386             } // fswitch
 387             break;
 388
 389         } // fswitch
 390
 391         _type = type;
 392         deallocate( );
 393         _creator = SELF;
 394         _raw = buffer;
 395         buildIndex( );
 396
 397     } // fi
 398 }
 399
 400 // -------------------------------------------------------------------------
 401 void kVolume::convertScale( Type type, double min, double max )
 402 {
 403     double tmin, tmax, smin, smax;
 404
 405 // PS ->     //tmin = GSL_MIN( min, max ); //PS
 406         tmin= ((min<max)?min:max);
 407 // PS ->     //tmax = GSL_MAX( min, max ); //PS
 408         tmax= ((min>max)?min:max);
 409
 410     getMinMax( smin, smax );
 411
 412     // compute scaling slope
 413     double a = ( tmax - tmin ) / ( smax - smin );
 414
 415     void* buffer;
 416     ulong size = getRawSize( );
 417
 418     buffer = ( void* )new uchar[ size * SIZETypes[ type ] ];
 419
 420     switch( _type ) {
 421
 422     case CHAR:
 423         switch( type ) {
 424
 425         case CHAR:   convertScaleT( ( char* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 426         case SHORT:  convertScaleT( ( char* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 427         case INT:    convertScaleT( ( char* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 428         case USHORT: convertScaleT( ( char* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 429         case UINT:   convertScaleT( ( char* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 430         case FLOAT:  convertScaleT( ( char* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 431         case DOUBLE: convertScaleT( ( char* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 432         case UCHAR:  convertScaleT( ( char* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 433
 434         } // fswitch
 435         break;
 436     case SHORT:
 437         switch( type ) {
 438
 439         case CHAR:   convertScaleT( ( short* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 440         case SHORT:  convertScaleT( ( short* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 441         case INT:    convertScaleT( ( short* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 442         case USHORT: convertScaleT( ( short* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 443         case UINT:   convertScaleT( ( short* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 444         case FLOAT:  convertScaleT( ( short* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 445         case DOUBLE: convertScaleT( ( short* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 446         case UCHAR:  convertScaleT( ( short* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 447
 448         } // fswitch
 449         break;
 450     case INT:
 451         switch( type ) {
 452
 453         case CHAR:   convertScaleT( ( int* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 454         case SHORT:  convertScaleT( ( int* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 455         case INT:    convertScaleT( ( int* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 456         case USHORT: convertScaleT( ( int* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 457         case UINT:   convertScaleT( ( int* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 458         case FLOAT:  convertScaleT( ( int* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 459         case DOUBLE: convertScaleT( ( int* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 460         case UCHAR:  convertScaleT( ( int* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 461
 462         } // fswitch
 463         break;
 464     case USHORT:
 465         switch( type ) {
 466
 467         case CHAR:   convertScaleT( ( ushort* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 468         case SHORT:  convertScaleT( ( ushort* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 469         case INT:    convertScaleT( ( ushort* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 470         case USHORT: convertScaleT( ( ushort* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 471         case UINT:   convertScaleT( ( ushort* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 472         case FLOAT:  convertScaleT( ( ushort* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 473         case DOUBLE: convertScaleT( ( ushort* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 474         case UCHAR:  convertScaleT( ( ushort* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 475
 476         } // fswitch
 477         break;
 478     case UINT:
 479         switch( type ) {
 480
 481         case CHAR:   convertScaleT( ( uint* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 482         case SHORT:  convertScaleT( ( uint* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 483         case INT:    convertScaleT( ( uint* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 484         case USHORT: convertScaleT( ( uint* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 485         case UINT:   convertScaleT( ( uint* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 486         case FLOAT:  convertScaleT( ( uint* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 487         case DOUBLE: convertScaleT( ( uint* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 488         case UCHAR:  convertScaleT( ( uint* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 489
 490         } // fswitch
 491         break;
 492     case UCHAR:
 493         switch( type ) {
 494
 495         case CHAR:   convertScaleT( ( uchar* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 496         case SHORT:  convertScaleT( ( uchar* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 497         case INT:    convertScaleT( ( uchar* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 498         case USHORT: convertScaleT( ( uchar* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 499         case UINT:   convertScaleT( ( uchar* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 500         case FLOAT:  convertScaleT( ( uchar* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 501         case DOUBLE: convertScaleT( ( uchar* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 502         case UCHAR:  convertScaleT( ( uchar* )_raw, ( uchar* )buffer, size, smin, tmin, a ); break;
 503
 504         } // fswitch
 505         break;
 506     case DOUBLE:
 507         switch( type ) {
 508
 509         case CHAR:   convertScaleT( ( double* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 510         case SHORT:  convertScaleT( ( double* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 511         case INT:    convertScaleT( ( double* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 512         case USHORT: convertScaleT( ( double* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 513         case UINT:   convertScaleT( ( double* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 514         case FLOAT:  convertScaleT( ( double* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 515         case DOUBLE: convertScaleT( ( double* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 516         case UCHAR:  convertScaleT( ( uchar* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 517
 518         } // fswitch
 519         break;
 520     case FLOAT:
 521         switch( type ) {
 522
 523         case CHAR:   convertScaleT( ( float* )_raw, ( char* )buffer, size, smin, tmin, a ); break;
 524         case SHORT:  convertScaleT( ( float* )_raw, ( short* )buffer, size, smin, tmin, a ); break;
 525         case INT:    convertScaleT( ( float* )_raw, ( int* )buffer, size, smin, tmin, a ); break;
 526         case USHORT: convertScaleT( ( float* )_raw, ( ushort* )buffer, size, smin, tmin, a ); break;
 527         case UINT:   convertScaleT( ( float* )_raw, ( uint* )buffer, size, smin, tmin, a ); break;
 528         case FLOAT:  convertScaleT( ( float* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
 529         case DOUBLE: convertScaleT( ( float* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 530         case UCHAR:  convertScaleT( ( float* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
 531
 532         } // fswitch
 533         break;
 534
 535     } // fswitch
 536
 537     _type = type;
 538     deallocate( );
 539     _creator = SELF;
 540     _raw = buffer;
 541     buildIndex( );
 542 }
 543
 544 // -------------------------------------------------------------------------
 545 void kVolume::getMinMax( double& min, double& max ) const
 546 {
 547     ulong size = getRawSize( );
 548
 549     switch( _type ) {
 550
 551     case CHAR:   getMinMaxT( ( char* )_raw, size, min, max ); break;
 552     case FLOAT:  getMinMaxT( ( float* )_raw, size, min, max ); break;
 553     case DOUBLE: getMinMaxT( ( double* )_raw, size, min, max ); break;
 554     case INT:    getMinMaxT( ( int* )_raw, size, min, max ); break;
 555     case SHORT:  getMinMaxT( ( short* )_raw, size, min, max ); break;
 556     case UCHAR:  getMinMaxT( ( uchar* )_raw, size, min, max ); break;
 557     case UINT:   getMinMaxT( ( uint* )_raw, size, min, max ); break;
 558     case USHORT: getMinMaxT( ( ushort* )_raw, size, min, max ); break;
 559
 560     } // fswitch
 561 }
 562
 563 // -------------------------------------------------------------------------
 564 double kVolume::getMin( ) const
 565 {
 566     double m, M;
 567
 568     getMinMax( m, M );
 569     return( m );
 570 }
 571
 572 // -------------------------------------------------------------------------
 573 double kVolume::getMax( ) const
 574 {
 575     double m, M;
 576
 577     getMinMax( m, M );
 578     return( M );
 579 }
 580
 581 // -------------------------------------------------------------------------
 582 /*double kVolume::GetMaxIntSphere( double* p, double r )
 583 {
 584     int minX, minY, minZ, maxX, maxY, maxZ;
 585     gslobj_vector vP( p, 3 ), v( 3 );
 586     double maxint, tmp;
 587     bool start;
 588
 589     minX = int( floor( p[ 0 ] - r ) );
 590     minY = int( floor( p[ 1 ] - r ) );
 591     minZ = int( floor( p[ 2 ] - r ) );
 592     maxX = int( ceil( p[ 0 ] + r ) );
 593     maxY = int( ceil( p[ 1 ] + r ) );
 594     maxZ = int( ceil( p[ 2 ] + r ) );
 595
 596     minX = GSL_MAX( minX, 0 );
 597     minY = GSL_MAX( minY, 0 );
 598     minZ = GSL_MAX( minZ, 0 );
 599
 600     maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );
 601     maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );
 602     maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );
 603
 604     maxint = 0.0;
 605     start = true;
 606     for( v( 0 ) = minX; v( 0 ) <= maxX; v( 0 )++ )
 607       for( v( 1 ) = minY; v( 1 ) <= maxY; v( 1 )++ )
 608         for( v( 2 ) = minZ; v( 2 ) <= maxZ; v( 2 )++ )
 609           if( ( v - vP ).norm2( ) <= r ) {
 610             tmp = this->getPixel( ( uint )v(0), ( uint )v(1), ( uint )v(2));
 611             maxint = ( tmp > maxint || start )? tmp: maxint;
 612             start = false;
 613           } // fi
 614
 615     return( maxint );
 616 }*/
 617 // -------------------------------------------------------------------------
 618 unsigned short kVolume::GetMaxIntSphere2( double* p, double r )
 619 {
 620    /**
 621     * unsigned short range : 0 -> 65535
 622     * unsigned int   range : 0 -> 2147483647 // 2^31 - 1
 623     */
 624     int minX, minY, minZ, maxX, maxY, maxZ;
 625     unsigned int v[3];
 626     unsigned short maxint = 0, tmp;
 627
 628     minX = int( floor( p[ 0 ] - r ) );
 629     minY = int( floor( p[ 1 ] - r ) );
 630     minZ = int( floor( p[ 2 ] - r ) );
 631     maxX = int( ceil( p[ 0 ] + r ) );
 632     maxY = int( ceil( p[ 1 ] + r ) );
 633     maxZ = int( ceil( p[ 2 ] + r ) );
 634
 635 // PS ->     //minX = GSL_MAX( minX, 0 );//PS
 636 // PS ->     //minY = GSL_MAX( minY, 0 );//PS
 637 // PS ->     //minZ = GSL_MAX( minZ, 0 );//PS
 638         minX=((minX>0)?minX:0);
 639         minY=((minY>0)?minY:0);
 640         minZ=((minZ>0)?minZ:0);
 641
 642 // PS ->     //maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );//PS
 643 // PS ->     //maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );//PS
 644 // PS ->     //maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );//PS
 645         int xDim=this->getXdim( ) - 1;
 646         maxX= (maxX<xDim?maxX:xDim);
 647
 648         int yDim=this->getYdim( ) - 1;
 649         maxY= (maxY<yDim?maxY:yDim);
 650
 651         int zDim=this->getZdim( ) - 1;
 652         maxZ= (maxZ<zDim?maxZ:zDim);
 653
 654     double r2 = r*r;  //need to do comparison in double ... don't know why...
 655     for( v[0] = minX; v[0] <= maxX; v[0]++ )
 656       for( v[1] = minY; v[1] <= maxY; v[1]++ )
 657         for( v[2] = minZ; v[2] <= maxZ; v[2]++ )
 658           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 )
 659           {
 660             tmp = (unsigned short)this->getPixel( v[0], v[1], v[2] );
 661             maxint = ( tmp > maxint ) ? tmp : maxint;
 662           }
 663
 664     return( maxint );
 665 }
 666
 667 // -------------------------------------------------------------------------
 668 void kVolume::allocate( )
 669 {
 670     ulong size = getRawSizeInBytes( );
 671
 672     if( _creator == SELF ) {
 673
 674         _raw = ( void* )new uchar[ size ];
 675         memset( _raw, 0, size );
 676
 677     } // fi
 678 }
 679
 680 // -------------------------------------------------------------------------
 681 void kVolume::buildIndex( )
 682 {
 683     ulong size;
 684
 685     size = ( _dims[ CZ ] * sizeof( uchar** ) ) +
 686         ( _dims[ CZ ] * _dims[ CY ] * sizeof( void* ) );
 687
 688         _images = ( void*** )new uchar[ size ];
 689
 690         _columns = ( void** )( _images + _dims[ CZ ] );
 691         void** plane = _columns;
 692         for( uint z = 0; z < _dims[ CZ ]; z++ ) {
 693
 694             _images[ z ] = plane;
 695             plane += _dims[ CY ];
 696
 697         } // rof
 698         void* line = _raw;
 699         for( uint y = 0; y < _dims[ CZ ] * _dims[ CY ]; y++ ) {
 700
 701             _columns[ y ] = line;
 702             line = ( void* )( ( uchar* ) line +
 703                               _dims[ CX ] * SIZETypes[ _type ] );
 704
 705         } // rof
 706
 707 #ifdef KGFO_USE_VTK
 708
 709     vtkCharArray* carray;
 710     vtkDoubleArray* darray;
 711     vtkFloatArray* farray;
 712     vtkIntArray* iarray;
 713     vtkShortArray* sarray;
 714     vtkUnsignedCharArray* ucarray;
 715     vtkUnsignedIntArray* uiarray;
 716     vtkUnsignedShortArray* usarray;
 717
 718     size = _dims[ CX ] * _dims[ CY ] * _dims[ CZ ];
 719
 720     if( _creator == SELF || _creator == IDO ) {
 721
 722         _vtk = vtkImageData::New( );
 723         _vtk->SetDimensions( _dims[ CX ], _dims[ CY ], _dims[ CZ ] );
 724         _vtk->SetSpacing( _sizes[ CX ], _sizes[ CY ], _sizes[ CZ ] );
 725
 726         switch( _type ) {
 727
 728         case CHAR:
 729             carray = vtkCharArray::New( );
 730             carray->SetNumberOfComponents( 1 );
 731             carray->SetArray( ( char* )( _raw ), size, 1 );
 732             _vtk->SetScalarType( VTK_CHAR );
 733             _vtk->GetPointData( )->SetScalars( carray );
 734             carray->Delete( );
 735             break;
 736
 737         case UCHAR:
 738             ucarray = vtkUnsignedCharArray::New( );
 739             ucarray->SetNumberOfComponents( 1 );
 740             ucarray->SetArray( ( uchar* )( _raw ), size, 1 );
 741             _vtk->SetScalarType( VTK_UNSIGNED_CHAR );
 742             _vtk->GetPointData( )->SetScalars( ucarray );
 743             ucarray->Delete( );
 744             break;
 745
 746         case SHORT:
 747             sarray = vtkShortArray::New( );
 748             sarray->SetNumberOfComponents( 1 );
 749             sarray->SetArray( ( short* )( _raw ), size, 1 );
 750             _vtk->SetScalarType( VTK_SHORT );
 751             _vtk->GetPointData( )->SetScalars( sarray );
 752             sarray->Delete( );
 753             break;
 754
 755         case INT:
 756             iarray = vtkIntArray::New( );
 757             iarray->SetNumberOfComponents( 1 );
 758             iarray->SetArray( ( int* )( _raw ), size, 1 );
 759             _vtk->SetScalarType( VTK_INT );
 760             _vtk->GetPointData( )->SetScalars( iarray );
 761             iarray->Delete( );
 762             break;
 763
 764         case USHORT:
 765             usarray = vtkUnsignedShortArray::New( );
 766             usarray->SetNumberOfComponents( 1 );
 767             usarray->SetArray( ( ushort* )( _raw ), size, 1 );
 768             _vtk->SetScalarType( VTK_UNSIGNED_SHORT );
 769             _vtk->GetPointData( )->SetScalars( usarray );
 770             usarray->Delete( );
 771             break;
 772
 773         case UINT:
 774             uiarray = vtkUnsignedIntArray::New( );
 775             uiarray->SetNumberOfComponents( 1 );
 776             uiarray->SetArray( ( uint* )( _raw ), size, 1 );
 777             _vtk->SetScalarType( VTK_UNSIGNED_INT );
 778             _vtk->GetPointData( )->SetScalars( uiarray );
 779             uiarray->Delete( );
 780             break;
 781
 782         case FLOAT:
 783             farray = vtkFloatArray::New( );
 784             farray->SetNumberOfComponents( 1 );
 785             farray->SetArray( ( float* )( _raw ), size, 1 );
 786             _vtk->SetScalarType( VTK_FLOAT );
 787             _vtk->GetPointData( )->SetScalars( farray );
 788             farray->Delete( );
 789             break;
 790
 791         case DOUBLE:
 792             darray = vtkDoubleArray::New( );
 793             darray->SetNumberOfComponents( 1 );
 794             darray->SetArray( ( double* )( _raw ), size, 1 );
 795             _vtk->SetScalarType( VTK_DOUBLE );
 796             _vtk->GetPointData( )->SetScalars( darray );
 797             darray->Delete( );
 798             break;
 799
 800         } // fswitch
 801
 802     } // fi
 803
 804 #endif // KGFO_USE_VTK
 805 }
 806
 807 // -------------------------------------------------------------------------
 808 void kVolume::deallocate( )
 809 {
 810 #ifdef KGFO_USE_VTK
 811
 812     if( _vtk ) _vtk->Delete();
 813     _vtk = NULL;
 814
 815 #endif // KGFO_USE_VTK
 816
 817 //#ifdef KGFO_USE_IDO
 818
 819 //    if( _creator == SELF || _creator == VTK ) {
 820
 821 //      delete[] _privateIdo;
 822 //      _privateIdo = NULL;
 823
 824 //    } // fi
 825
 826 //#else
 827
 828     delete[] ( uchar* )_images;
 829
 830 //#endif // KGFO_USE_IDO
 831
 832     if( _raw && _creator == SELF )
 833
 834 //EED purify 12/sept/2006
 835 //      delete[] ( uchar* )_raw;
 836
 837         free ( _raw );
 838
 839     _creator    = SELF;
 840     _raw                = NULL;
 841     _columns    = NULL;
 842     _images             = NULL;
 843 }
 844
 845 #ifdef KGFO_USE_VTK
 846
 847 // -------------------------------------------------------------------------
 848 kVolume::kVolume( vtkImageData* org )
 849     : _raw( 0 ), _columns( 0 ), _images( 0 ),
 850       _creator( VTK ),
 851 //#ifdef KGFO_USE_IDO
 852 //      _privateIdo( NULL ),
 853 //#endif // KGFO_USE_IDO
 854       _vtk( 0 )
 855 {
 856     //int i, j, k, y;
 857     int itmp[ 3 ];
 858     double ftmp[ 3 ];
 859     //double v;
 860
 861     switch( org->GetScalarType( ) ) {
 862
 863     case VTK_CHAR:           _type = CHAR;   break;
 864     case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
 865     case VTK_SHORT:          _type = SHORT;  break;
 866     case VTK_INT:            _type = INT;    break;
 867     case VTK_UNSIGNED_SHORT: _type = USHORT; break;
 868     case VTK_UNSIGNED_INT:   _type = UINT;   break;
 869     case VTK_FLOAT:          _type = FLOAT;  break;
 870     case VTK_DOUBLE:         _type = DOUBLE; break;
 871     default: break;
 872
 873     } // fswitch
 874
 875     org->GetDimensions( itmp );
 876     _dims[ CX ] = ( uint )itmp[ 0 ];
 877     _dims[ CY ] = ( uint )itmp[ 1 ];
 878     _dims[ CZ ] = ( uint )itmp[ 2 ];
 879     org->GetSpacing( ftmp );
 880     _sizes[ CX ] = ( double )ftmp[ 0 ];
 881     _sizes[ CY ] = ( double )ftmp[ 1 ];
 882     _sizes[ CZ ] = ( double )ftmp[ 2 ];
 883
 884     _raw = org->GetPointData( )->GetScalars( )->GetVoidPointer( 0 );
 885     //_vtk = org;
 886     _vtk = vtkImageData::New();
 887     _vtk->ShallowCopy( org );
 888     buildIndex( );
 889 }
 890
 891 // -------------------------------------------------------------------------
 892 kVolume& kVolume::operator=( vtkImageData* org )
 893 {
 894     copyFrom( org );
 895     return( *this );
 896 }
 897
 898 // -------------------------------------------------------------------------
 899 void kVolume::copyFrom( vtkImageData* org )
 900 {
 901     int i, j, k;//, y;
 902     int itmp[ 3 ];
 903     int ext[ 6 ];
 904     double ftmp[ 3 ];
 905     double v;
 906
 907     deallocate( );
 908
 909 //#ifdef KGFO_USE_IDO
 910
 911 //    _privateIdo = NULL;
 912
 913 //#endif // KGFO_USE_IDO
 914
 915     _vtk                = NULL;
 916     _raw                = NULL;
 917     _columns    = NULL;
 918     _images             = NULL;
 919     _creator    = SELF;
 920
 921     switch( org->GetScalarType( ) ) {
 922
 923     case VTK_CHAR:           _type = CHAR;   break;
 924     case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
 925     case VTK_SHORT:          _type = SHORT;  break;
 926     case VTK_INT:            _type = INT;    break;
 927     case VTK_UNSIGNED_SHORT: _type = USHORT; break;
 928     case VTK_UNSIGNED_INT:   _type = UINT;   break;
 929     case VTK_FLOAT:          _type = FLOAT;  break;
 930     case VTK_DOUBLE:         _type = DOUBLE; break;
 931     default: break;
 932
 933     } // fswitch
 934
 935     org->GetDimensions( itmp );
 936     _dims[ CX ] = ( uint )itmp[ 0 ];
 937     _dims[ CY ] = ( uint )itmp[ 1 ];
 938     _dims[ CZ ] = ( uint )itmp[ 2 ];
 939     org->GetSpacing( ftmp );
 940     _sizes[ CX ] = ( double )ftmp[ 0 ];
 941     _sizes[ CY ] = ( double )ftmp[ 1 ];
 942     _sizes[ CZ ] = ( double )ftmp[ 2 ];
 943
 944     allocate( );
 945     buildIndex( );
 946
 947     // This avoids vtk extent crap conversion...
 948     org->GetExtent( ext );
 949     for( i = ext[ 0 ]; i <= ext[ 1 ]; i++ ) {
 950         for( j = ext[ 2 ]; j <= ext[ 3 ]; j++ ) {
 951             for( k = ext[ 4 ]; k <= ext[ 5 ]; k++ ) {
 952                 v = org->GetScalarComponentAsDouble( i, j, k, 0 );
 953                 setPixel( v, i - ext[ 0 ], j - ext[ 2 ], k - ext[ 4 ] );
 954             } // rof
 955         } // rof
 956     } // rof
 957 }
 958
 959 #endif // KGFO_USE_VTK
 960
 961
 962 // eof - volume.cxx