lib/maracasVisuLib/src/kernel/volume.cxx

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