lib/maracasVisuLib/src/kernel/volume.cxx

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