1 /*# ---------------------------------------------------------------------
3 # Copyright (c) CREATIS (Centre de Recherche en Acquisition et Traitement de l'Image
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
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.
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
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 # ------------------------------------------------------------------------ */
26 /*=========================================================================
29 Module: $RCSfile: volume.cxx,v $
31 Date: $Date: 2012/11/15 14:16:13 $
32 Version: $Revision: 1.10 $
34 Copyright: (c) 2002, 2003
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.
41 =========================================================================*/
43 // PS -> //#include <gsl/gsl_math.h>//PS
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>
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 };
65 #endif // KGFO_USE_VTK
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 ) };
75 // ---------------------------------------------------------------------------
76 template< class FROM, class TO >
77 static void convertCastT( FROM* src, TO* dest, ulong size )
79 FROM* end = src + size;
80 for( ; src < end; src++, dest++ ) *dest = ( TO )*src;
84 // ---------------------------------------------------------------------------
85 template< class FROM, class TO >
86 static void convertScaleT( FROM *src, TO *dest, ulong size,
87 double smin, double tmin, double slope )
89 FROM* end = src + size;
90 for( ; src < end; src++, dest++ )
91 *dest = ( TO )( ( double( *src ) - smin ) * slope + tmin );
94 // ---------------------------------------------------------------------------
95 template< class TYPE >
96 static void getMinMaxT( TYPE *src, ulong size,
97 double& min, double& max )
99 TYPE* end = src + size;
105 for( st = true; src < end; src++, st = false ) {
107 if( *src < m || st ) m = *src;
108 if( *src > M || st ) M = *src;
116 // -------------------------------------------------------------------------
126 #endif // KGFO_USE_VTK
128 _dims[ CX ] = 1; _dims[ CY ] = 1; _dims[ CZ ] = 1;
129 _sizes[ CX ] = 1; _sizes[ CY ] = 1; _sizes[ CZ ] = 1;
134 // -------------------------------------------------------------------------
135 kVolume::kVolume( Type type,
136 uint32_t xdim, uint32_t ydim, uint32_t zdim,
137 double xsize, double ysize, double zsize,
147 #endif // KGFO_USE_VTK
149 _dims[ CX ] = xdim; _dims[ CY ] = ydim; _dims[ CZ ] = zdim;
150 _sizes[ CX ] = xsize; _sizes[ CY ] = ysize; _sizes[ CZ ] = zsize;
151 if( data != NOALLOC ) {
162 // -------------------------------------------------------------------------
163 kVolume::kVolume( Type type,
164 const uint32_t *dims,
169 _raw( NULL ), _columns( NULL ),
174 #endif // KGFO_USE_VTK
177 memcpy( _dims, dims, 3 * sizeof( uint32_t ) );
178 memcpy( _sizes, sizes, 3 * sizeof( double ) );
179 if( data != NOALLOC ) {
190 // -------------------------------------------------------------------------
191 kVolume::kVolume( const kVolume& org )
200 #endif // KGFO_USE_VTK
206 // -------------------------------------------------------------------------
207 kVolume& kVolume::operator=( const kVolume& org )
213 // -------------------------------------------------------------------------
214 void kVolume::copyFrom( const kVolume& org )
222 memcpy( _dims, org._dims, 3 * sizeof( uint32_t ) );
223 memcpy( _sizes, org._sizes, 3 * sizeof( double ) );
228 memcpy( _raw, org._raw, getRawSizeInBytes( ) );
233 // -------------------------------------------------------------------------
234 bool kVolume::sameDimension( const kVolume& comp )
236 return( ( _dims[ CX ] == comp._dims[ CX ] &&
237 _dims[ CY ] == comp._dims[ CY ] &&
238 _dims[ CZ ] == comp._dims[ CZ ] ) );
241 // -------------------------------------------------------------------------
242 double kVolume::getPixel( uint32_t x, uint32_t y, uint32_t z ) const
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;
263 // -------------------------------------------------------------------------
264 void kVolume::setPixel( double v, uint32_t x, uint32_t y, uint32_t z )
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;
282 // -------------------------------------------------------------------------
283 void kVolume::convertCast( Type type )
285 if( _type != type ) {
288 ulong size = getRawSize( );
290 buffer = ( void* )new uint8_t[ size * SIZETypes[ type ] ];
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;
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;
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;
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;
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;
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;
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;
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;
418 // -------------------------------------------------------------------------
419 void kVolume::convertScale( Type type, double min, double max )
421 double tmin, tmax, smin, smax;
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);
428 getMinMax( smin, smax );
430 // compute scaling slope
431 double a = ( tmax - tmin ) / ( smax - smin );
434 ulong size = getRawSize( );
436 buffer = ( void* )new uint8_t[ size * SIZETypes[ type ] ];
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;
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;
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;
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;
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;
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;
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;
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;
562 // -------------------------------------------------------------------------
563 void kVolume::getMinMax( double& min, double& max ) const
565 ulong size = getRawSize( );
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;
581 // -------------------------------------------------------------------------
582 double kVolume::getMin( ) const
590 // -------------------------------------------------------------------------
591 double kVolume::getMax( ) const
599 // -------------------------------------------------------------------------
600 /*double kVolume::GetMaxIntSphere( double* p, double r )
602 int minX, minY, minZ, maxX, maxY, maxZ;
603 gslobj_vector vP( p, 3 ), v( 3 );
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 ) );
614 minX = GSL_MAX( minX, 0 );
615 minY = GSL_MAX( minY, 0 );
616 minZ = GSL_MAX( minZ, 0 );
618 maxX = GSL_MIN( maxX, this->getXdim( ) - 1 );
619 maxY = GSL_MIN( maxY, this->getYdim( ) - 1 );
620 maxZ = GSL_MIN( maxZ, this->getZdim( ) - 1 );
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;
635 // -------------------------------------------------------------------------
636 unsigned short kVolume::GetMaxIntSphere2( double* p, double r )
639 * unsigned short range : 0 -> 65535
640 * unsigned int range : 0 -> 2147483647 // 2^31 - 1
642 int minX, minY, minZ, maxX, maxY, maxZ;
644 unsigned short maxint = 0, tmp;
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 ) );
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);
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);
666 int yDim=this->getYdim( ) - 1;
667 maxY= (maxY<yDim?maxY:yDim);
669 int zDim=this->getZdim( ) - 1;
670 maxZ= (maxZ<zDim?maxZ:zDim);
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 )
678 tmp = (unsigned short)this->getPixel( v[0], v[1], v[2] );
679 maxint = ( tmp > maxint ) ? tmp : maxint;
685 // -------------------------------------------------------------------------
686 void kVolume::allocate( )
688 ulong size = getRawSizeInBytes( );
690 if( _creator == SELF ) {
692 _raw = ( void* )new uint8_t[ size ];
693 memset( _raw, 0, size );
698 // -------------------------------------------------------------------------
699 void kVolume::buildIndex( )
703 size = ( _dims[ CZ ] * sizeof( uint8_t** ) ) +
704 ( _dims[ CZ ] * _dims[ CY ] * sizeof( void* ) );
706 _images = ( void*** )new uint8_t[ size ];
708 _columns = ( void** )( _images + _dims[ CZ ] );
709 void** plane = _columns;
710 for( uint32_t z = 0; z < _dims[ CZ ]; z++ ) {
712 _images[ z ] = plane;
713 plane += _dims[ CY ];
717 for( uint32_t y = 0; y < _dims[ CZ ] * _dims[ CY ]; y++ ) {
719 _columns[ y ] = line;
720 line = ( void* )( ( uint8_t* ) line +
721 _dims[ CX ] * SIZETypes[ _type ] );
727 vtkCharArray* carray;
728 vtkDoubleArray* darray;
729 vtkFloatArray* farray;
731 vtkShortArray* sarray;
732 vtkUnsignedCharArray* ucarray;
733 vtkUnsignedIntArray* uiarray;
734 vtkUnsignedShortArray* usarray;
736 size = _dims[ CX ] * _dims[ CY ] * _dims[ CZ ];
738 if( _creator == SELF || _creator == IDO ) {
740 _vtk = vtkImageData::New( );
741 _vtk->SetDimensions( _dims[ CX ], _dims[ CY ], _dims[ CZ ] );
742 _vtk->SetSpacing( _sizes[ CX ], _sizes[ CY ], _sizes[ CZ ] );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
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 );
822 #endif // KGFO_USE_VTK
825 // -------------------------------------------------------------------------
826 void kVolume::deallocate( )
829 if( _vtk ) _vtk->Delete();
831 #endif // KGFO_USE_VTK
832 delete[] ( uint8_t* )_images;
833 if( _raw && _creator == SELF )
835 //EED purify 12/sept/2006
836 // delete[] ( uint8_t* )_raw;
848 // -------------------------------------------------------------------------
849 kVolume::kVolume( vtkImageData* org )
862 switch( org->GetScalarType( ) ) {
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;
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 ];
885 _raw = org->GetPointData( )->GetScalars( )->GetVoidPointer( 0 );
887 _vtk = vtkImageData::New();
888 _vtk->ShallowCopy( org );
892 // -------------------------------------------------------------------------
893 kVolume& kVolume::operator=( vtkImageData* org )
899 // -------------------------------------------------------------------------
900 void kVolume::copyFrom( vtkImageData* org )
910 //#ifdef KGFO_USE_IDO
912 // _privateIdo = NULL;
914 //#endif // KGFO_USE_IDO
922 switch( org->GetScalarType( ) ) {
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;
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 ];
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 ] );
960 #endif // KGFO_USE_VTK