lib/maracasVisuLib/src/kernel/volume.cxx

   1 /*=========================================================================
   2
   3   Program:   wxMaracas
   4   Module:    $RCSfile: volume.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2008/10/31 16:32:56 $
   7   Version:   $Revision: 1.1 $
   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 #ifdef KGFO_USE_IDO
 689
 690     if( _creator == SELF || _creator == VTK ) {
 691
 692         size += sizeof( PRIVATE_VOLUME );
 693
 694         _privateIdo = ( PRIVATE_VOLUME* ) new uchar[ size ];
 695         _privateIdo->UsedNbZ = 0;  // Warning, I don't really know
 696         _privateIdo->UsedNbY = 0;  // the reason to use these three
 697         _privateIdo->UsedNbX = 0;  // fields. - lflorez
 698         _privateIdo->subObject = 1;
 699         _privateIdo->DimX = _dims[ CX ];
 700         _privateIdo->DimY = _dims[ CY ];
 701         _privateIdo->DimZ = _dims[ CZ ];
 702         _privateIdo->_message = 0;
 703         _privateIdo->_fichier = 0;
 704         _privateIdo->Type = IDOTypes[ _type ];
 705
 706         _images = ( void*** ) &( _privateIdo[ 1 ] );
 707
 708 #else
 709
 710         _images = ( void*** )new uchar[ size ];
 711
 712 #endif // KGFO_USE_IDO
 713
 714         _columns = ( void** )( _images + _dims[ CZ ] );
 715         void** plane = _columns;
 716         for( uint z = 0; z < _dims[ CZ ]; z++ ) {
 717
 718             _images[ z ] = plane;
 719             plane += _dims[ CY ];
 720
 721         } // rof
 722         void* line = _raw;
 723         for( uint y = 0; y < _dims[ CZ ] * _dims[ CY ]; y++ ) {
 724
 725             _columns[ y ] = line;
 726             line = ( void* )( ( uchar* ) line +
 727                               _dims[ CX ] * SIZETypes[ _type ] );
 728
 729         } // rof
 730
 731 #ifdef KGFO_USE_IDO
 732
 733     } else
 734         _images = ( void*** )( &_privateIdo[ 1 ] );
 735
 736 #endif // KGFO_USE_IDO
 737
 738 #ifdef KGFO_USE_VTK
 739
 740     vtkCharArray* carray;
 741     vtkDoubleArray* darray;
 742     vtkFloatArray* farray;
 743     vtkIntArray* iarray;
 744     vtkShortArray* sarray;
 745     vtkUnsignedCharArray* ucarray;
 746     vtkUnsignedIntArray* uiarray;
 747     vtkUnsignedShortArray* usarray;
 748
 749     size = _dims[ CX ] * _dims[ CY ] * _dims[ CZ ];
 750
 751     if( _creator == SELF || _creator == IDO ) {
 752
 753         _vtk = vtkImageData::New( );
 754         _vtk->SetDimensions( _dims[ CX ], _dims[ CY ], _dims[ CZ ] );
 755         _vtk->SetSpacing( _sizes[ CX ], _sizes[ CY ], _sizes[ CZ ] );
 756
 757         switch( _type ) {
 758
 759         case CHAR:
 760             carray = vtkCharArray::New( );
 761             carray->SetNumberOfComponents( 1 );
 762             carray->SetArray( ( char* )( _raw ), size, 1 );
 763             _vtk->SetScalarType( VTK_CHAR );
 764             _vtk->GetPointData( )->SetScalars( carray );
 765             carray->Delete( );
 766             break;
 767
 768         case UCHAR:
 769             ucarray = vtkUnsignedCharArray::New( );
 770             ucarray->SetNumberOfComponents( 1 );
 771             ucarray->SetArray( ( uchar* )( _raw ), size, 1 );
 772             _vtk->SetScalarType( VTK_UNSIGNED_CHAR );
 773             _vtk->GetPointData( )->SetScalars( ucarray );
 774             ucarray->Delete( );
 775             break;
 776
 777         case SHORT:
 778             sarray = vtkShortArray::New( );
 779             sarray->SetNumberOfComponents( 1 );
 780             sarray->SetArray( ( short* )( _raw ), size, 1 );
 781             _vtk->SetScalarType( VTK_SHORT );
 782             _vtk->GetPointData( )->SetScalars( sarray );
 783             sarray->Delete( );
 784             break;
 785
 786         case INT:
 787             iarray = vtkIntArray::New( );
 788             iarray->SetNumberOfComponents( 1 );
 789             iarray->SetArray( ( int* )( _raw ), size, 1 );
 790             _vtk->SetScalarType( VTK_INT );
 791             _vtk->GetPointData( )->SetScalars( iarray );
 792             iarray->Delete( );
 793             break;
 794
 795         case USHORT:
 796             usarray = vtkUnsignedShortArray::New( );
 797             usarray->SetNumberOfComponents( 1 );
 798             usarray->SetArray( ( ushort* )( _raw ), size, 1 );
 799             _vtk->SetScalarType( VTK_UNSIGNED_SHORT );
 800             _vtk->GetPointData( )->SetScalars( usarray );
 801             usarray->Delete( );
 802             break;
 803
 804         case UINT:
 805             uiarray = vtkUnsignedIntArray::New( );
 806             uiarray->SetNumberOfComponents( 1 );
 807             uiarray->SetArray( ( uint* )( _raw ), size, 1 );
 808             _vtk->SetScalarType( VTK_UNSIGNED_INT );
 809             _vtk->GetPointData( )->SetScalars( uiarray );
 810             uiarray->Delete( );
 811             break;
 812
 813         case FLOAT:
 814             farray = vtkFloatArray::New( );
 815             farray->SetNumberOfComponents( 1 );
 816             farray->SetArray( ( float* )( _raw ), size, 1 );
 817             _vtk->SetScalarType( VTK_FLOAT );
 818             _vtk->GetPointData( )->SetScalars( farray );
 819             farray->Delete( );
 820             break;
 821
 822         case DOUBLE:
 823             darray = vtkDoubleArray::New( );
 824             darray->SetNumberOfComponents( 1 );
 825             darray->SetArray( ( double* )( _raw ), size, 1 );
 826             _vtk->SetScalarType( VTK_DOUBLE );
 827             _vtk->GetPointData( )->SetScalars( darray );
 828             darray->Delete( );
 829             break;
 830
 831         } // fswitch
 832
 833     } // fi
 834
 835 #endif // KGFO_USE_VTK
 836 }
 837
 838 // -------------------------------------------------------------------------
 839 void kVolume::deallocate( )
 840 {
 841 #ifdef KGFO_USE_VTK
 842
 843     if( _vtk ) _vtk->Delete();
 844     _vtk = NULL;
 845
 846 #endif // KGFO_USE_VTK
 847
 848 #ifdef KGFO_USE_IDO
 849
 850     if( _creator == SELF || _creator == VTK ) {
 851
 852         delete[] _privateIdo;
 853         _privateIdo = NULL;
 854
 855     } // fi
 856
 857 #else
 858
 859     delete[] ( uchar* )_images;
 860
 861 #endif // KGFO_USE_IDO
 862
 863     if( _raw && _creator == SELF )
 864
 865 //EED purify 12/sept/2006
 866 //      delete[] ( uchar* )_raw;
 867
 868         free ( _raw );
 869
 870     _creator    = SELF;
 871     _raw                = NULL;
 872     _columns    = NULL;
 873     _images             = NULL;
 874 }
 875
 876 #ifdef KGFO_USE_VTK
 877
 878 // -------------------------------------------------------------------------
 879 kVolume::kVolume( vtkImageData* org )
 880     : _raw( 0 ), _columns( 0 ), _images( 0 ),
 881       _creator( VTK ),
 882 #ifdef KGFO_USE_IDO
 883       _privateIdo( NULL ),
 884 #endif // KGFO_USE_IDO
 885       _vtk( 0 )
 886 {
 887     //int i, j, k, y;
 888     int itmp[ 3 ];
 889     double ftmp[ 3 ];
 890     //double v;
 891
 892     switch( org->GetScalarType( ) ) {
 893
 894     case VTK_CHAR:           _type = CHAR;   break;
 895     case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
 896     case VTK_SHORT:          _type = SHORT;  break;
 897     case VTK_INT:            _type = INT;    break;
 898     case VTK_UNSIGNED_SHORT: _type = USHORT; break;
 899     case VTK_UNSIGNED_INT:   _type = UINT;   break;
 900     case VTK_FLOAT:          _type = FLOAT;  break;
 901     case VTK_DOUBLE:         _type = DOUBLE; break;
 902     default: break;
 903
 904     } // fswitch
 905
 906     org->GetDimensions( itmp );
 907     _dims[ CX ] = ( uint )itmp[ 0 ];
 908     _dims[ CY ] = ( uint )itmp[ 1 ];
 909     _dims[ CZ ] = ( uint )itmp[ 2 ];
 910     org->GetSpacing( ftmp );
 911     _sizes[ CX ] = ( double )ftmp[ 0 ];
 912     _sizes[ CY ] = ( double )ftmp[ 1 ];
 913     _sizes[ CZ ] = ( double )ftmp[ 2 ];
 914
 915     _raw = org->GetPointData( )->GetScalars( )->GetVoidPointer( 0 );
 916     //_vtk = org;
 917     _vtk = vtkImageData::New();
 918     _vtk->ShallowCopy( org );
 919     buildIndex( );
 920 }
 921
 922 // -------------------------------------------------------------------------
 923 kVolume& kVolume::operator=( vtkImageData* org )
 924 {
 925     copyFrom( org );
 926     return( *this );
 927 }
 928
 929 // -------------------------------------------------------------------------
 930 void kVolume::copyFrom( vtkImageData* org )
 931 {
 932     int i, j, k;//, y;
 933     int itmp[ 3 ];
 934     int ext[ 6 ];
 935     double ftmp[ 3 ];
 936     double v;
 937
 938     deallocate( );
 939
 940 #ifdef KGFO_USE_IDO
 941
 942     _privateIdo = NULL;
 943
 944 #endif // KGFO_USE_IDO
 945
 946     _vtk                = NULL;
 947     _raw                = NULL;
 948     _columns    = NULL;
 949     _images             = NULL;
 950     _creator    = SELF;
 951
 952     switch( org->GetScalarType( ) ) {
 953
 954     case VTK_CHAR:           _type = CHAR;   break;
 955     case VTK_UNSIGNED_CHAR:  _type = UCHAR;  break;
 956     case VTK_SHORT:          _type = SHORT;  break;
 957     case VTK_INT:            _type = INT;    break;
 958     case VTK_UNSIGNED_SHORT: _type = USHORT; break;
 959     case VTK_UNSIGNED_INT:   _type = UINT;   break;
 960     case VTK_FLOAT:          _type = FLOAT;  break;
 961     case VTK_DOUBLE:         _type = DOUBLE; break;
 962     default: break;
 963
 964     } // fswitch
 965
 966     org->GetDimensions( itmp );
 967     _dims[ CX ] = ( uint )itmp[ 0 ];
 968     _dims[ CY ] = ( uint )itmp[ 1 ];
 969     _dims[ CZ ] = ( uint )itmp[ 2 ];
 970     org->GetSpacing( ftmp );
 971     _sizes[ CX ] = ( double )ftmp[ 0 ];
 972     _sizes[ CY ] = ( double )ftmp[ 1 ];
 973     _sizes[ CZ ] = ( double )ftmp[ 2 ];
 974
 975     allocate( );
 976     buildIndex( );
 977
 978     // This avoids vtk extent crap conversion...
 979     org->GetExtent( ext );
 980     for( i = ext[ 0 ]; i <= ext[ 1 ]; i++ ) {
 981
 982         for( j = ext[ 2 ]; j <= ext[ 3 ]; j++ ) {
 983
 984             for( k = ext[ 4 ]; k <= ext[ 5 ]; k++ ) {
 985
 986                 v = org->GetScalarComponentAsDouble( i, j, k, 0 );
 987                 setPixel( v, i - ext[ 0 ], j - ext[ 2 ], k - ext[ 4 ] );
 988
 989             } // rof
 990
 991         } // rof
 992
 993     } // rof
 994 }
 995
 996 #endif // KGFO_USE_VTK
 997
 998 #ifdef KGFO_USE_IDO
 999
1000 // -------------------------------------------------------------------------
1001 kVolume::kVolume( PPPVOLUME org )
1002     : _raw( 0 ), _columns( 0 ), _images( 0 ),
1003       _creator( IDO ),
1004 #ifdef KGFO_USE_VTK
1005       _vtk( 0 ),
1006 #endif // KGFO_USE_VTK
1007       _privateIdo( 0 )
1008 {
1009     switch( IdVolType( org ) ) {
1010
1011     case VOL_UCHAR:  _type = UCHAR;  break;
1012     case VOL_CHAR:   _type = CHAR;   break;
1013     case VOL_FLOAT:  _type = FLOAT;  break;
1014     case VOL_DOUBLE: _type = DOUBLE; break;
1015     case VOL_SHORT:  _type = SHORT;  break;
1016     case VOL_USHORT: _type = USHORT; break;
1017     case VOL_LONG:   _type = INT;    break;
1018     case VOL_ULONG:  _type = UINT;   break;
1019     default: break;
1020
1021     } // fswitch
1022
1023     _dims[ CX ] = IdVolDimX( org );
1024     _dims[ CY ] = IdVolDimY( org );
1025     _dims[ CZ ] = IdVolDimZ( org );
1026     _sizes[ CX ] = 1;
1027     _sizes[ CY ] = 1;
1028     _sizes[ CZ ] = 1;
1029
1030     _privateIdo = _IdVolPrivate( org );
1031     _raw = ( ( void*** ) &( _privateIdo[ 1 ] ) )[ 0 ][ 0 ];
1032     buildIndex( );
1033 }
1034
1035 // -------------------------------------------------------------------------
1036 kVolume& kVolume::operator=( PPPVOLUME org )
1037 {
1038     copyFrom( org );
1039     return( *this );
1040 }
1041
1042 // -------------------------------------------------------------------------
1043 void kVolume::copyFrom( PPPVOLUME org )
1044 {
1045     void* buffer;
1046
1047     deallocate( );
1048
1049     _raw                = NULL;
1050     _columns    = NULL;
1051     _images             = NULL;
1052     _creator    = SELF;
1053 #ifdef KGFO_USE_VTK
1054     _vtk = NULL;
1055 #endif // KGFO_USE_VTK
1056     _privateIdo = NULL;
1057
1058     switch( IdVolType( org ) ) {
1059
1060     case VOL_UCHAR:  _type = UCHAR;  break;
1061     case VOL_CHAR:   _type = CHAR;   break;
1062     case VOL_FLOAT:  _type = FLOAT;  break;
1063     case VOL_DOUBLE: _type = DOUBLE; break;
1064     case VOL_SHORT:  _type = SHORT;  break;
1065     case VOL_USHORT: _type = USHORT; break;
1066     case VOL_LONG:   _type = INT;    break;
1067     case VOL_ULONG:  _type = UINT;   break;
1068     default: break;
1069
1070     } // fswitch
1071
1072     _dims[ CX ] = IdVolDimX( org );
1073     _dims[ CY ] = IdVolDimY( org );
1074     _dims[ CZ ] = IdVolDimZ( org );
1075     _sizes[ CX ] = 1;
1076     _sizes[ CY ] = 1;
1077     _sizes[ CZ ] = 1;
1078
1079     allocate( );
1080     buildIndex( );
1081     buffer = ( ( void*** ) &( org[ 1 ] ) )[ 0 ][ 0 ];
1082     memcpy( _raw, buffer, getRawSizeInBytes( ) );
1083 }
1084
1085 #endif // KGFO_USE_IDO
1086
1087 // eof - volume.cxx