lib/maracasVisuLib/src/kernel/volume.cxx

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