lib/maracasVisuLib/src/kernel/volume.cxx

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