1 /*=========================================================================
4 Module: $RCSfile: volume.cxx,v $
6 Date: $Date: 2010/04/29 16:05:37 $
7 Version: $Revision: 1.9 $
9 Copyright: (c) 2002, 2003
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.
16 =========================================================================*/
18 // PS -> //#include <gsl/gsl_math.h>//PS
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>
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 };
40 #endif // KGFO_USE_VTK
42 // -------------------------------------------------------------------------
43 const void* kVolume::BLANK = ( void* ) 0;
44 const void* kVolume::NOALLOC = ( void* )-1;
45 const int kVolume::SIZETypes[] = { sizeof( char ), sizeof( float ),
46 sizeof( double ), sizeof( int ),
47 sizeof( short ), sizeof( uint8_t ),
48 sizeof( uint32_t ), sizeof( uint16_t ) };
50 // ---------------------------------------------------------------------------
51 template< class FROM, class TO >
52 static void convertCastT( FROM* src, TO* dest, ulong size )
54 FROM* end = src + size;
55 for( ; src < end; src++, dest++ ) *dest = ( TO )*src;
59 // ---------------------------------------------------------------------------
60 template< class FROM, class TO >
61 static void convertScaleT( FROM *src, TO *dest, ulong size,
62 double smin, double tmin, double slope )
64 FROM* end = src + size;
65 for( ; src < end; src++, dest++ )
66 *dest = ( TO )( ( double( *src ) - smin ) * slope + tmin );
69 // ---------------------------------------------------------------------------
70 template< class TYPE >
71 static void getMinMaxT( TYPE *src, ulong size,
72 double& min, double& max )
74 TYPE* end = src + size;
80 for( st = true; src < end; src++, st = false ) {
82 if( *src < m || st ) m = *src;
83 if( *src > M || st ) M = *src;
91 // -------------------------------------------------------------------------
101 #endif // KGFO_USE_VTK
103 _dims[ CX ] = 1; _dims[ CY ] = 1; _dims[ CZ ] = 1;
104 _sizes[ CX ] = 1; _sizes[ CY ] = 1; _sizes[ CZ ] = 1;
109 // -------------------------------------------------------------------------
110 kVolume::kVolume( Type type,
111 uint32_t xdim, uint32_t ydim, uint32_t zdim,
112 double xsize, double ysize, double zsize,
122 #endif // KGFO_USE_VTK
124 _dims[ CX ] = xdim; _dims[ CY ] = ydim; _dims[ CZ ] = zdim;
125 _sizes[ CX ] = xsize; _sizes[ CY ] = ysize; _sizes[ CZ ] = zsize;
126 if( data != NOALLOC ) {
137 // -------------------------------------------------------------------------
138 kVolume::kVolume( Type type,
139 const uint32_t *dims,
144 _raw( NULL ), _columns( NULL ),
149 #endif // KGFO_USE_VTK
152 memcpy( _dims, dims, 3 * sizeof( uint32_t ) );
153 memcpy( _sizes, sizes, 3 * sizeof( double ) );
154 if( data != NOALLOC ) {
165 // -------------------------------------------------------------------------
166 kVolume::kVolume( const kVolume& org )
175 #endif // KGFO_USE_VTK
181 // -------------------------------------------------------------------------
182 kVolume& kVolume::operator=( const kVolume& org )
188 // -------------------------------------------------------------------------
189 void kVolume::copyFrom( const kVolume& org )
197 memcpy( _dims, org._dims, 3 * sizeof( uint32_t ) );
198 memcpy( _sizes, org._sizes, 3 * sizeof( double ) );
203 memcpy( _raw, org._raw, getRawSizeInBytes( ) );
208 // -------------------------------------------------------------------------
209 bool kVolume::sameDimension( const kVolume& comp )
211 return( ( _dims[ CX ] == comp._dims[ CX ] &&
212 _dims[ CY ] == comp._dims[ CY ] &&
213 _dims[ CZ ] == comp._dims[ CZ ] ) );
216 // -------------------------------------------------------------------------
217 double kVolume::getPixel( uint32_t x, uint32_t y, uint32_t z ) const
223 case CHAR: p = ( double )( ( int8_t*** )_images )[ z ][ y ][ x ]; break;
224 case FLOAT: p = ( double )( ( float*** )_images )[ z ][ y ][ x ]; break;
225 case DOUBLE: p = ( double )( ( double*** )_images )[ z ][ y ][ x ]; break;
226 case INT: p = ( double )( ( int32_t*** )_images )[ z ][ y ][ x ]; break;
227 case SHORT: p = ( double )( ( int16_t*** )_images )[ z ][ y ][ x ]; break;
228 case UCHAR: p = ( double )( ( uint8_t*** )_images )[ z ][ y ][ x ]; break;
229 case UINT: p = ( double )( ( uint32_t*** )_images )[ z ][ y ][ x ]; break;
230 case USHORT: p = ( double )( ( uint16_t*** )_images )[ z ][ y ][ x ]; break;
231 default: p = 0.0; break;
238 // -------------------------------------------------------------------------
239 void kVolume::setPixel( double v, uint32_t x, uint32_t y, uint32_t z )
244 case CHAR: ( ( int8_t*** )_images )[ z ][ y ][ x ] = ( int8_t )v; break;
245 case FLOAT: ( ( float*** )_images )[ z ][ y ][ x ] = ( float )v; break;
246 case DOUBLE: ( ( double*** )_images )[ z ][ y ][ x ] = ( double )v; break;
247 case INT: ( ( int32_t*** )_images )[ z ][ y ][ x ] = ( int32_t )v; break;
248 case SHORT: ( ( int16_t*** )_images )[ z ][ y ][ x ] = ( int16_t )v; break;
249 case UCHAR: ( ( uint8_t*** )_images )[ z ][ y ][ x ] = ( uint8_t )v; break;
250 case UINT: ( ( uint32_t*** )_images )[ z ][ y ][ x ] = ( uint32_t )v; break;
251 case USHORT: ( ( uint16_t*** )_images )[ z ][ y ][ x ] = ( uint16_t )v; break;
257 // -------------------------------------------------------------------------
258 void kVolume::convertCast( Type type )
260 if( _type != type ) {
263 ulong size = getRawSize( );
265 buffer = ( void* )new uint8_t[ size * SIZETypes[ type ] ];
272 case SHORT: convertCastT( ( char* )_raw, ( int16_t* )buffer, size ); break;
273 case INT: convertCastT( ( char* )_raw, ( int32_t* )buffer, size ); break;
274 case USHORT: convertCastT( ( char* )_raw, ( uint16_t* )buffer, size ); break;
275 case UINT: convertCastT( ( char* )_raw, ( uint32_t* )buffer, size ); break;
276 case FLOAT: convertCastT( ( char* )_raw, ( float* )buffer, size ); break;
277 case DOUBLE: convertCastT( ( char* )_raw, ( double* )buffer, size ); break;
278 case UCHAR: convertCastT( ( char* )_raw, ( uint8_t* )buffer, size ); break;
287 case CHAR: convertCastT( ( int16_t* )_raw, ( uint8_t* )buffer, size ); break;
288 case INT: convertCastT( ( int16_t* )_raw, ( int32_t* )buffer, size ); break;
289 case USHORT: convertCastT( ( int16_t* )_raw, ( uint16_t* )buffer, size ); break;
290 case UINT: convertCastT( ( int16_t* )_raw, ( uint32_t* )buffer, size ); break;
291 case FLOAT: convertCastT( ( int16_t* )_raw, ( float* )buffer, size ); break;
292 case DOUBLE: convertCastT( ( int16_t* )_raw, ( double* )buffer, size ); break;
293 case UCHAR: convertCastT( ( int16_t* )_raw, ( uint8_t* )buffer, size ); break;
301 case CHAR: convertCastT( ( int32_t* )_raw, ( int8_t* )buffer, size ); break;
302 case SHORT: convertCastT( ( int32_t* )_raw, ( int16_t* )buffer, size ); break;
303 case USHORT: convertCastT( ( int32_t* )_raw, ( uint16_t* )buffer, size ); break;
304 case UINT: convertCastT( ( int32_t* )_raw, ( uint32_t* )buffer, size ); break;
305 case FLOAT: convertCastT( ( int32_t* )_raw, ( float* )buffer, size ); break;
306 case DOUBLE: convertCastT( ( int32_t* )_raw, ( double* )buffer, size ); break;
307 case UCHAR: convertCastT( ( int32_t* )_raw, ( uint8_t* )buffer, size ); break;
315 case CHAR: convertCastT( ( uint16_t* )_raw, ( int8_t* )buffer, size ); break;
316 case SHORT: convertCastT( ( uint16_t* )_raw, ( int16_t* )buffer, size ); break;
317 case INT: convertCastT( ( uint16_t* )_raw, ( int32_t* )buffer, size ); break;
318 case UINT: convertCastT( ( uint16_t* )_raw, ( uint32_t* )buffer, size ); break;
319 case FLOAT: convertCastT( ( uint16_t* )_raw, ( float* )buffer, size ); break;
320 case DOUBLE: convertCastT( ( uint16_t* )_raw, ( double* )buffer, size ); break;
321 case UCHAR: convertCastT( ( uint16_t* )_raw, ( uint8_t* )buffer, size ); break;
329 case CHAR: convertCastT( ( uint32_t* )_raw, ( int8_t* )buffer, size ); break;
330 case SHORT: convertCastT( ( uint32_t* )_raw, ( int16_t* )buffer, size ); break;
331 case INT: convertCastT( ( uint32_t* )_raw, ( int32_t* )buffer, size ); break;
332 case USHORT: convertCastT( ( uint32_t* )_raw, ( uint16_t* )buffer, size ); break;
333 case FLOAT: convertCastT( ( uint32_t* )_raw, ( float* )buffer, size ); break;
334 case DOUBLE: convertCastT( ( uint32_t* )_raw, ( double* )buffer, size ); break;
335 case UCHAR: convertCastT( ( uint32_t* )_raw, ( uint8_t* )buffer, size ); break;
343 case CHAR: convertCastT( ( float* )_raw, ( int8_t* )buffer, size ); break;
344 case SHORT: convertCastT( ( float* )_raw, ( int16_t* )buffer, size ); break;
345 case INT: convertCastT( ( float* )_raw, ( int32_t* )buffer, size ); break;
346 case USHORT: convertCastT( ( float* )_raw, ( uint16_t* )buffer, size ); break;
347 case UINT: convertCastT( ( float* )_raw, ( uint32_t* )buffer, size ); break;
348 case DOUBLE: convertCastT( ( float* )_raw, ( double* )buffer, size ); break;
349 case UCHAR: convertCastT( ( float* )_raw, ( uint8_t* )buffer, size ); break;
357 case CHAR: convertCastT( ( double* )_raw, ( int8_t* )buffer, size ); break;
358 case SHORT: convertCastT( ( double* )_raw, ( int16_t* )buffer, size ); break;
359 case INT: convertCastT( ( double* )_raw, ( int32_t* )buffer, size ); break;
360 case USHORT: convertCastT( ( double* )_raw, ( uint16_t* )buffer, size ); break;
361 case UINT: convertCastT( ( double* )_raw, ( uint32_t* )buffer, size ); break;
362 case FLOAT: convertCastT( ( double* )_raw, ( float* )buffer, size ); break;
363 case UCHAR: convertCastT( ( double* )_raw, ( uint8_t* )buffer, size ); break;
371 case CHAR: convertCastT( ( uint8_t* )_raw, ( int8_t* )buffer, size ); break;
372 case SHORT: convertCastT( ( uint8_t* )_raw, ( int16_t* )buffer, size ); break;
373 case INT: convertCastT( ( uint8_t* )_raw, ( int32_t* )buffer, size ); break;
374 case USHORT: convertCastT( ( uint8_t* )_raw, ( uint16_t* )buffer, size ); break;
375 case UINT: convertCastT( ( uint8_t* )_raw, ( uint32_t* )buffer, size ); break;
376 case FLOAT: convertCastT( ( uint8_t* )_raw, ( float* )buffer, size ); break;
377 case DOUBLE: convertCastT( ( uint8_t* )_raw, ( double* )buffer, size ); break;
393 // -------------------------------------------------------------------------
394 void kVolume::convertScale( Type type, double min, double max )
396 double tmin, tmax, smin, smax;
398 // PS -> //tmin = GSL_MIN( min, max ); //PS
399 tmin= ((min<max)?min:max);
400 // PS -> //tmax = GSL_MAX( min, max ); //PS
401 tmax= ((min>max)?min:max);
403 getMinMax( smin, smax );
405 // compute scaling slope
406 double a = ( tmax - tmin ) / ( smax - smin );
409 ulong size = getRawSize( );
411 buffer = ( void* )new uint8_t[ size * SIZETypes[ type ] ];
418 case CHAR: convertScaleT( ( int8_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
419 case SHORT: convertScaleT( ( int8_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
420 case INT: convertScaleT( ( int8_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
421 case USHORT: convertScaleT( ( int8_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
422 case UINT: convertScaleT( ( int8_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
423 case FLOAT: convertScaleT( ( int8_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
424 case DOUBLE: convertScaleT( ( int8_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
425 case UCHAR: convertScaleT( ( int8_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
432 case CHAR: convertScaleT( ( int16_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
433 case SHORT: convertScaleT( ( int16_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
434 case INT: convertScaleT( ( int16_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
435 case USHORT: convertScaleT( ( int16_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
436 case UINT: convertScaleT( ( int16_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
437 case FLOAT: convertScaleT( ( int16_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
438 case DOUBLE: convertScaleT( ( int16_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
439 case UCHAR: convertScaleT( ( int16_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
446 case CHAR: convertScaleT( ( int32_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
447 case SHORT: convertScaleT( ( int32_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
448 case INT: convertScaleT( ( int32_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
449 case USHORT: convertScaleT( ( int32_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
450 case UINT: convertScaleT( ( int32_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
451 case FLOAT: convertScaleT( ( int32_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
452 case DOUBLE: convertScaleT( ( int32_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
453 case UCHAR: convertScaleT( ( int32_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
460 case CHAR: convertScaleT( ( uint16_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
461 case SHORT: convertScaleT( ( uint16_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
462 case INT: convertScaleT( ( uint16_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
463 case USHORT: convertScaleT( ( uint16_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
464 case UINT: convertScaleT( ( uint16_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
465 case FLOAT: convertScaleT( ( uint16_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
466 case DOUBLE: convertScaleT( ( uint16_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
467 case UCHAR: convertScaleT( ( uint16_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
474 case CHAR: convertScaleT( ( uint32_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
475 case SHORT: convertScaleT( ( uint32_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
476 case INT: convertScaleT( ( uint32_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
477 case USHORT: convertScaleT( ( uint32_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
478 case UINT: convertScaleT( ( uint32_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
479 case FLOAT: convertScaleT( ( uint32_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
480 case DOUBLE: convertScaleT( ( uint32_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
481 case UCHAR: convertScaleT( ( uint32_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
488 case CHAR: convertScaleT( ( uint8_t* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
489 case SHORT: convertScaleT( ( uint8_t* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
490 case INT: convertScaleT( ( uint8_t* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
491 case USHORT: convertScaleT( ( uint8_t* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
492 case UINT: convertScaleT( ( uint8_t* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
493 case FLOAT: convertScaleT( ( uint8_t* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
494 case DOUBLE: convertScaleT( ( uint8_t* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
495 case UCHAR: convertScaleT( ( uint8_t* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
502 case CHAR: convertScaleT( ( double* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
503 case SHORT: convertScaleT( ( double* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
504 case INT: convertScaleT( ( double* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
505 case USHORT: convertScaleT( ( double* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
506 case UINT: convertScaleT( ( double* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
507 case FLOAT: convertScaleT( ( double* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
508 case DOUBLE: convertScaleT( ( double* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
509 case UCHAR: convertScaleT( ( double* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
516 case CHAR: convertScaleT( ( float* )_raw, ( int8_t* )buffer, size, smin, tmin, a ); break;
517 case SHORT: convertScaleT( ( float* )_raw, ( int16_t* )buffer, size, smin, tmin, a ); break;
518 case INT: convertScaleT( ( float* )_raw, ( int32_t* )buffer, size, smin, tmin, a ); break;
519 case USHORT: convertScaleT( ( float* )_raw, ( uint16_t* )buffer, size, smin, tmin, a ); break;
520 case UINT: convertScaleT( ( float* )_raw, ( uint32_t* )buffer, size, smin, tmin, a ); break;
521 case FLOAT: convertScaleT( ( float* )_raw, ( float* )buffer, size, smin, tmin, a ); break;
522 case DOUBLE: convertScaleT( ( float* )_raw, ( double* )buffer, size, smin, tmin, a ); break;
523 case UCHAR: convertScaleT( ( float* )_raw, ( uint8_t* )buffer, size, smin, tmin, a ); break;
537 // -------------------------------------------------------------------------
538 void kVolume::getMinMax( double& min, double& max ) const
540 ulong size = getRawSize( );
544 case CHAR: getMinMaxT( ( int8_t* )_raw, size, min, max ); break;
545 case FLOAT: getMinMaxT( ( float* )_raw, size, min, max ); break;
546 case DOUBLE: getMinMaxT( ( double* )_raw, size, min, max ); break;
547 case INT: getMinMaxT( ( int32_t* )_raw, size, min, max ); break;
548 case SHORT: getMinMaxT( ( int16_t* )_raw, size, min, max ); break;
549 case UCHAR: getMinMaxT( ( uint8_t* )_raw, size, min, max ); break;
550 case UINT: getMinMaxT( ( uint32_t* )_raw, size, min, max ); break;
551 case USHORT: getMinMaxT( ( uint16_t* )_raw, size, min, max ); break;
556 // -------------------------------------------------------------------------
557 double kVolume::getMin( ) const
565 // -------------------------------------------------------------------------
566 double kVolume::getMax( ) const
574 // -------------------------------------------------------------------------
575 /*double kVolume::GetMaxIntSphere( double* p, double r )
577 int minX, minY, minZ, maxX, maxY, maxZ;
578 gslobj_vector vP( p, 3 ), v( 3 );
582 minX = int( floor( p[ 0 ] - r ) );
583 minY = int( floor( p[ 1 ] - r ) );
584 minZ = int( floor( p[ 2 ] - r ) );
585 maxX = int( ceil( p[ 0 ] + r ) );
586 maxY = int( ceil( p[ 1 ] + r ) );
587 maxZ = int( ceil( p[ 2 ] + r ) );
589 minX = GSL_MAX( minX, 0 );
590 minY = GSL_MAX( minY, 0 );
591 minZ = GSL_MAX( minZ, 0 );
593 maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );
594 maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );
595 maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );
599 for( v( 0 ) = minX; v( 0 ) <= maxX; v( 0 )++ )
600 for( v( 1 ) = minY; v( 1 ) <= maxY; v( 1 )++ )
601 for( v( 2 ) = minZ; v( 2 ) <= maxZ; v( 2 )++ )
602 if( ( v - vP ).norm2( ) <= r ) {
603 tmp = this->getPixel( ( uint32_t )v(0), ( uint32_t )v(1), ( uint32_t )v(2));
604 maxint = ( tmp > maxint || start )? tmp: maxint;
610 // -------------------------------------------------------------------------
611 unsigned short kVolume::GetMaxIntSphere2( double* p, double r )
614 * unsigned short range : 0 -> 65535
615 * unsigned int range : 0 -> 2147483647 // 2^31 - 1
617 int minX, minY, minZ, maxX, maxY, maxZ;
619 unsigned short maxint = 0, tmp;
621 minX = int( floor( p[ 0 ] - r ) );
622 minY = int( floor( p[ 1 ] - r ) );
623 minZ = int( floor( p[ 2 ] - r ) );
624 maxX = int( ceil( p[ 0 ] + r ) );
625 maxY = int( ceil( p[ 1 ] + r ) );
626 maxZ = int( ceil( p[ 2 ] + r ) );
628 // PS -> //minX = GSL_MAX( minX, 0 );//PS
629 // PS -> //minY = GSL_MAX( minY, 0 );//PS
630 // PS -> //minZ = GSL_MAX( minZ, 0 );//PS
631 minX=((minX>0)?minX:0);
632 minY=((minY>0)?minY:0);
633 minZ=((minZ>0)?minZ:0);
635 // PS -> //maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );//PS
636 // PS -> //maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );//PS
637 // PS -> //maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );//PS
638 int xDim=this->getXdim( ) - 1;
639 maxX= (maxX<xDim?maxX:xDim);
641 int yDim=this->getYdim( ) - 1;
642 maxY= (maxY<yDim?maxY:yDim);
644 int zDim=this->getZdim( ) - 1;
645 maxZ= (maxZ<zDim?maxZ:zDim);
647 double r2 = r*r; //need to do comparison in double ... don't know why...
648 for( v[0] = minX; v[0] <= maxX; v[0]++ )
649 for( v[1] = minY; v[1] <= maxY; v[1]++ )
650 for( v[2] = minZ; v[2] <= maxZ; v[2]++ )
651 if( ((v[0]-p[0])*(v[0]-p[0])+(v[1]-p[1])*(v[1]-p[1])+(v[2]-p[2])*(v[2]-p[2])) <= r2 )
653 tmp = (unsigned short)this->getPixel( v[0], v[1], v[2] );
654 maxint = ( tmp > maxint ) ? tmp : maxint;
660 // -------------------------------------------------------------------------
661 void kVolume::allocate( )
663 ulong size = getRawSizeInBytes( );
665 if( _creator == SELF ) {
667 _raw = ( void* )new uint8_t[ size ];
668 memset( _raw, 0, size );
673 // -------------------------------------------------------------------------
674 void kVolume::buildIndex( )
678 size = ( _dims[ CZ ] * sizeof( uint8_t** ) ) +
679 ( _dims[ CZ ] * _dims[ CY ] * sizeof( void* ) );
681 _images = ( void*** )new uint8_t[ size ];
683 _columns = ( void** )( _images + _dims[ CZ ] );
684 void** plane = _columns;
685 for( uint32_t z = 0; z < _dims[ CZ ]; z++ ) {
687 _images[ z ] = plane;
688 plane += _dims[ CY ];
692 for( uint32_t y = 0; y < _dims[ CZ ] * _dims[ CY ]; y++ ) {
694 _columns[ y ] = line;
695 line = ( void* )( ( uint8_t* ) line +
696 _dims[ CX ] * SIZETypes[ _type ] );
702 vtkCharArray* carray;
703 vtkDoubleArray* darray;
704 vtkFloatArray* farray;
706 vtkShortArray* sarray;
707 vtkUnsignedCharArray* ucarray;
708 vtkUnsignedIntArray* uiarray;
709 vtkUnsignedShortArray* usarray;
711 size = _dims[ CX ] * _dims[ CY ] * _dims[ CZ ];
713 if( _creator == SELF || _creator == IDO ) {
715 _vtk = vtkImageData::New( );
716 _vtk->SetDimensions( _dims[ CX ], _dims[ CY ], _dims[ CZ ] );
717 _vtk->SetSpacing( _sizes[ CX ], _sizes[ CY ], _sizes[ CZ ] );
722 carray = vtkCharArray::New( );
723 carray->SetNumberOfComponents( 1 );
724 carray->SetArray( ( char* )( _raw ), size, 1 );
725 _vtk->SetScalarType( VTK_CHAR );
726 _vtk->GetPointData( )->SetScalars( carray );
731 ucarray = vtkUnsignedCharArray::New( );
732 ucarray->SetNumberOfComponents( 1 );
733 ucarray->SetArray( ( uint8_t* )( _raw ), size, 1 );
734 _vtk->SetScalarType( VTK_UNSIGNED_CHAR );
735 _vtk->GetPointData( )->SetScalars( ucarray );
740 sarray = vtkShortArray::New( );
741 sarray->SetNumberOfComponents( 1 );
742 sarray->SetArray( ( int16_t* )( _raw ), size, 1 );
743 _vtk->SetScalarType( VTK_SHORT );
744 _vtk->GetPointData( )->SetScalars( sarray );
749 iarray = vtkIntArray::New( );
750 iarray->SetNumberOfComponents( 1 );
751 iarray->SetArray( ( int32_t* )( _raw ), size, 1 );
752 _vtk->SetScalarType( VTK_INT );
753 _vtk->GetPointData( )->SetScalars( iarray );
758 usarray = vtkUnsignedShortArray::New( );
759 usarray->SetNumberOfComponents( 1 );
760 usarray->SetArray( ( uint16_t* )( _raw ), size, 1 );
761 _vtk->SetScalarType( VTK_UNSIGNED_SHORT );
762 _vtk->GetPointData( )->SetScalars( usarray );
767 uiarray = vtkUnsignedIntArray::New( );
768 uiarray->SetNumberOfComponents( 1 );
769 uiarray->SetArray( ( uint32_t* )( _raw ), size, 1 );
770 _vtk->SetScalarType( VTK_UNSIGNED_INT );
771 _vtk->GetPointData( )->SetScalars( uiarray );
776 farray = vtkFloatArray::New( );
777 farray->SetNumberOfComponents( 1 );
778 farray->SetArray( ( float* )( _raw ), size, 1 );
779 _vtk->SetScalarType( VTK_FLOAT );
780 _vtk->GetPointData( )->SetScalars( farray );
785 darray = vtkDoubleArray::New( );
786 darray->SetNumberOfComponents( 1 );
787 darray->SetArray( ( double* )( _raw ), size, 1 );
788 _vtk->SetScalarType( VTK_DOUBLE );
789 _vtk->GetPointData( )->SetScalars( darray );
797 #endif // KGFO_USE_VTK
800 // -------------------------------------------------------------------------
801 void kVolume::deallocate( )
804 if( _vtk ) _vtk->Delete();
806 #endif // KGFO_USE_VTK
807 delete[] ( uint8_t* )_images;
808 if( _raw && _creator == SELF )
810 //EED purify 12/sept/2006
811 // delete[] ( uint8_t* )_raw;
823 // -------------------------------------------------------------------------
824 kVolume::kVolume( vtkImageData* org )
837 switch( org->GetScalarType( ) ) {
839 case VTK_CHAR: _type = CHAR; break;
840 case VTK_UNSIGNED_CHAR: _type = UCHAR; break;
841 case VTK_SHORT: _type = SHORT; break;
842 case VTK_INT: _type = INT; break;
843 case VTK_UNSIGNED_SHORT: _type = USHORT; break;
844 case VTK_UNSIGNED_INT: _type = UINT; break;
845 case VTK_FLOAT: _type = FLOAT; break;
846 case VTK_DOUBLE: _type = DOUBLE; break;
851 org->GetDimensions( itmp );
852 _dims[ CX ] = ( uint32_t )itmp[ 0 ];
853 _dims[ CY ] = ( uint32_t )itmp[ 1 ];
854 _dims[ CZ ] = ( uint32_t )itmp[ 2 ];
855 org->GetSpacing( ftmp );
856 _sizes[ CX ] = ( double )ftmp[ 0 ];
857 _sizes[ CY ] = ( double )ftmp[ 1 ];
858 _sizes[ CZ ] = ( double )ftmp[ 2 ];
860 _raw = org->GetPointData( )->GetScalars( )->GetVoidPointer( 0 );
862 _vtk = vtkImageData::New();
863 _vtk->ShallowCopy( org );
867 // -------------------------------------------------------------------------
868 kVolume& kVolume::operator=( vtkImageData* org )
874 // -------------------------------------------------------------------------
875 void kVolume::copyFrom( vtkImageData* org )
885 //#ifdef KGFO_USE_IDO
887 // _privateIdo = NULL;
889 //#endif // KGFO_USE_IDO
897 switch( org->GetScalarType( ) ) {
899 case VTK_CHAR: _type = CHAR; break;
900 case VTK_UNSIGNED_CHAR: _type = UCHAR; break;
901 case VTK_SHORT: _type = SHORT; break;
902 case VTK_INT: _type = INT; break;
903 case VTK_UNSIGNED_SHORT: _type = USHORT; break;
904 case VTK_UNSIGNED_INT: _type = UINT; break;
905 case VTK_FLOAT: _type = FLOAT; break;
906 case VTK_DOUBLE: _type = DOUBLE; break;
911 org->GetDimensions( itmp );
912 _dims[ CX ] = ( uint32_t )itmp[ 0 ];
913 _dims[ CY ] = ( uint32_t )itmp[ 1 ];
914 _dims[ CZ ] = ( uint32_t )itmp[ 2 ];
915 org->GetSpacing( ftmp );
916 _sizes[ CX ] = ( double )ftmp[ 0 ];
917 _sizes[ CY ] = ( double )ftmp[ 1 ];
918 _sizes[ CZ ] = ( double )ftmp[ 2 ];
923 // This avoids vtk extent crap conversion...
924 org->GetExtent( ext );
925 for( i = ext[ 0 ]; i <= ext[ 1 ]; i++ ) {
926 for( j = ext[ 2 ]; j <= ext[ 3 ]; j++ ) {
927 for( k = ext[ 4 ]; k <= ext[ 5 ]; k++ ) {
928 v = org->GetScalarComponentAsDouble( i, j, k, 0 );
929 setPixel( v, i - ext[ 0 ], j - ext[ 2 ], k - ext[ 4 ] );
935 #endif // KGFO_USE_VTK