Mercurial > hg > minc-tools
changeset 1089:b37a3ef7098a
check_in_all
| author | david <david> |
|---|---|
| date | Wed, 13 Aug 1997 13:21:28 +0000 |
| parents | b609c92cd668 |
| children | 4dd5e5ebf714 |
| files | volume_io/Geometry/transforms.c volume_io/Include/volume_io/alloc.h volume_io/Include/volume_io/basic.h volume_io/Include/volume_io/system_dependent.h volume_io/Include/volume_io/volume.h volume_io/Prog_utils/files.c volume_io/Volumes/get_hyperslab.c volume_io/Volumes/input_free.c volume_io/Volumes/input_mnc.c volume_io/Volumes/output_mnc.c volume_io/Volumes/output_volume.c volume_io/Volumes/volumes.c |
| diffstat | 12 files changed, 978 insertions(+), 393 deletions(-) [+] |
line wrap: on
line diff
--- a/volume_io/Geometry/transforms.c +++ b/volume_io/Geometry/transforms.c @@ -15,7 +15,7 @@ #include <internal_volume_io.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Geometry/transforms.c,v 1.18 1996-05-17 19:36:08 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Geometry/transforms.c,v 1.19 1997-08-13 13:21:28 david Exp $"; #endif /* ----------------------------- MNI Header ----------------------------------- @@ -142,6 +142,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : get_transform_origin_real +@INPUT : transform +@OUTPUT : origin +@RETURNS : +@DESCRIPTION: Passes back the origin of the transform, i.e., where the + point (0,0,0) would be transformed to. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void get_transform_origin_real( + Transform *transform, + Real origin[] ) +{ + origin[X] = Transform_elem(*transform,0,3); + origin[Y] = Transform_elem(*transform,1,3); + origin[Z] = Transform_elem(*transform,2,3); +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : get_transform_x_axis @INPUT : transform @OUTPUT : x_axis @@ -166,6 +189,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : get_transform_x_axis_real +@INPUT : transform +@OUTPUT : x_axis +@RETURNS : +@DESCRIPTION: Passes back the x axis of the transform, i.e., the vector + to which the vector (1,0,0) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void get_transform_x_axis_real( + Transform *transform, + Real x_axis[] ) +{ + x_axis[X] = Transform_elem(*transform,0,0); + x_axis[Y] = Transform_elem(*transform,1,0); + x_axis[Z] = Transform_elem(*transform,2,0); +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : set_transform_x_axis @INPUT : x_axis @OUTPUT : transform @@ -189,6 +235,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : set_transform_x_axis_real +@INPUT : x_axis +@OUTPUT : transform +@RETURNS : +@DESCRIPTION: Sets the x axis of the transform, i.e., the vector + to which the vector (1,0,0) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void set_transform_x_axis_real( + Transform *transform, + Real x_axis[] ) +{ + Transform_elem(*transform,0,0) = (Transform_elem_type) x_axis[X]; + Transform_elem(*transform,1,0) = (Transform_elem_type) x_axis[Y]; + Transform_elem(*transform,2,0) = (Transform_elem_type) x_axis[Z]; +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : get_transform_y_axis @INPUT : transform @OUTPUT : y_axis @@ -213,6 +282,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : get_transform_y_axis_real +@INPUT : transform +@OUTPUT : y_axis +@RETURNS : +@DESCRIPTION: Passes back the y axis of the transform, i.e., the vector + to which the vector (0,1,0) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void get_transform_y_axis_real( + Transform *transform, + Real y_axis[] ) +{ + y_axis[X] = Transform_elem(*transform,0,1); + y_axis[Y] = Transform_elem(*transform,1,1); + y_axis[Z] = Transform_elem(*transform,2,1); +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : set_transform_y_axis @INPUT : y_axis @OUTPUT : transform @@ -236,6 +328,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : set_transform_y_axis_real +@INPUT : y_axis +@OUTPUT : transform +@RETURNS : +@DESCRIPTION: Sets the y axis of the transform, i.e., the vector + to which the vector (0,1,0) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void set_transform_y_axis_real( + Transform *transform, + Real y_axis[] ) +{ + Transform_elem(*transform,0,1) = (Transform_elem_type) y_axis[X]; + Transform_elem(*transform,1,1) = (Transform_elem_type) y_axis[Y]; + Transform_elem(*transform,2,1) = (Transform_elem_type) y_axis[Z]; +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : get_transform_z_axis @INPUT : transform @OUTPUT : z_axis @@ -260,6 +375,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : get_transform_z_axis_real +@INPUT : transform +@OUTPUT : z_axis +@RETURNS : +@DESCRIPTION: Passes back the z axis of the transform, i.e., the vector + to which the vector (0,0,1) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void get_transform_z_axis_real( + Transform *transform, + Real z_axis[] ) +{ + z_axis[X] = Transform_elem(*transform,0,2); + z_axis[Y] = Transform_elem(*transform,1,2); + z_axis[Z] = Transform_elem(*transform,2,2); +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : set_transform_z_axis @INPUT : z_axis @OUTPUT : transform @@ -283,6 +421,29 @@ } /* ----------------------------- MNI Header ----------------------------------- +@NAME : set_transform_z_axis_real +@INPUT : z_axis +@OUTPUT : transform +@RETURNS : +@DESCRIPTION: Sets the z axis of the transform, i.e., the vector + to which the vector (0,0,1) would be transformed. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void set_transform_z_axis_real( + Transform *transform, + Real z_axis[] ) +{ + Transform_elem(*transform,0,2) = (Transform_elem_type) z_axis[X]; + Transform_elem(*transform,1,2) = (Transform_elem_type) z_axis[Y]; + Transform_elem(*transform,2,2) = (Transform_elem_type) z_axis[Z]; +} + +/* ----------------------------- MNI Header ----------------------------------- @NAME : make_change_to_bases_transform @INPUT : origin x_axis
--- a/volume_io/Include/volume_io/alloc.h +++ b/volume_io/Include/volume_io/alloc.h @@ -29,7 +29,7 @@ ---------------------------------------------------------------------------- */ #ifndef lint -static char alloc_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/alloc.h,v 1.14 1996-11-15 16:09:43 david Exp $"; +static char alloc_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/alloc.h,v 1.15 1997-08-13 13:21:30 david Exp $"; #endif #include <basic.h> @@ -275,7 +275,7 @@ ---------------------------------------------------------------------------- */ #define ALLOC5D( ptr, n1, n2, n3, n4, n5 ) \ - ASSIGN(ptr) = alloc_memory_5d( (size_t) (n1), (size_t) (n2), \ + ASSIGN_PTR(ptr) = alloc_memory_5d( (size_t) (n1), (size_t) (n2), \ (size_t) (n3), (size_t) (n4), (size_t) (n5), \ sizeof(*****(ptr)) _ALLOC_SOURCE_LINE )
--- a/volume_io/Include/volume_io/basic.h +++ b/volume_io/Include/volume_io/basic.h @@ -29,7 +29,7 @@ ---------------------------------------------------------------------------- */ #ifndef lint -static char basic_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/basic.h,v 1.29 1997-03-23 21:11:29 david Exp $"; +static char basic_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/basic.h,v 1.30 1997-08-13 13:21:30 david Exp $"; #endif #include <stdlib.h> @@ -89,14 +89,14 @@ #endif #define MAX( x, y ) ( ((x) >= (y)) ? (x) : (y) ) -#define MAX3( x, y, z ) MAX( x, MAX(y,z) ) +#define MAX3( x, y, z ) ( ((x) >= (y)) ? MAX( x, z ) : MAX( y, z ) ) #ifdef MIN #undef MIN #endif #define MIN( x, y ) ( ((x) <= (y)) ? (x) : (y) ) -#define MIN3( x, y, z ) MIN( x, MIN(y,z) ) +#define MIN3( x, y, z ) ( ((x) <= (y)) ? MIN( x, z ) : MIN( y, z ) ) /* --------- gets the address of a 2-d array element in a 1-d array ----- */
--- a/volume_io/Include/volume_io/system_dependent.h +++ b/volume_io/Include/volume_io/system_dependent.h @@ -14,6 +14,13 @@ express or implied warranty. ---------------------------------------------------------------------------- */ +#ifdef sun +#ifndef __GNUC__ +#define NO_DBL_MAX /*--- sun os does not define DBL_MAX */ +#endif +#define NO_STRERROR /*--- sun os does not define strerror */ +#endif + #ifdef NO_DBL_MAX #ifndef DBL_MAX #include <values.h>
--- a/volume_io/Include/volume_io/volume.h +++ b/volume_io/Include/volume_io/volume.h @@ -16,7 +16,7 @@ ---------------------------------------------------------------------------- */ #ifndef lint -static char volume_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/volume.h,v 1.48 1996-12-09 20:20:23 david Exp $"; +static char volume_rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Include/volume_io/volume.h,v 1.49 1997-08-13 13:21:29 david Exp $"; #endif /* ----------------------------- MNI Header ----------------------------------- @@ -42,6 +42,8 @@ { Real global_image_range[2]; STRING dimension_names[MAX_DIMENSIONS]; + BOOLEAN use_starts_set; + BOOLEAN use_volume_starts_and_steps; } minc_output_options; #include <volume_cache.h> @@ -75,11 +77,10 @@ Real real_value_translation; Real separations[MAX_DIMENSIONS]; - Real translation_voxel[MAX_DIMENSIONS]; + Real starts[MAX_DIMENSIONS]; Real direction_cosines[MAX_DIMENSIONS][N_DIMENSIONS]; - Real world_space_for_translation_voxel[N_DIMENSIONS]; - + BOOLEAN voxel_to_world_transform_uptodate; General_transform voxel_to_world_transform; STRING coordinate_system_name;
--- a/volume_io/Prog_utils/files.c +++ b/volume_io/Prog_utils/files.c @@ -19,7 +19,7 @@ #include <errno.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Prog_utils/files.c,v 1.35 1997-06-29 19:00:21 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Prog_utils/files.c,v 1.36 1997-08-13 13:21:31 david Exp $"; #endif private BOOLEAN has_no_extension( STRING );
--- a/volume_io/Volumes/get_hyperslab.c +++ b/volume_io/Volumes/get_hyperslab.c @@ -15,7 +15,7 @@ #include <internal_volume_io.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/get_hyperslab.c,v 1.5 1996-11-15 16:09:48 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/get_hyperslab.c,v 1.6 1997-08-13 13:21:35 david Exp $"; #endif public void convert_voxels_to_values( @@ -198,42 +198,32 @@ } } -private Real *unsigned_byte_to_real = NULL; -private Real *signed_byte_to_real = NULL; -private Real *unsigned_short_to_real = NULL; -private Real *signed_short_to_real = NULL; +private Real *int_to_real_conversion = NULL; -private void check_real_conversion_lookup( - Data_types data_type ) +private void check_real_conversion_lookup( void ) { - Real min_value, max_value, **ptr; + Real min_value1, max_value1, min_value2, max_value2; long i, long_min, long_max; - switch( data_type ) - { - case UNSIGNED_BYTE: ptr = &unsigned_byte_to_real; break; - case SIGNED_BYTE: ptr = &signed_byte_to_real; break; - case UNSIGNED_SHORT: ptr = &unsigned_short_to_real; break; - case SIGNED_SHORT: ptr = &signed_short_to_real; break; - default: return; - } - - if( *ptr != NULL ) + if( int_to_real_conversion != NULL ) return; - get_type_range( data_type, &min_value, &max_value ); - long_min = (long) min_value; - long_max = (long) max_value; + get_type_range( UNSIGNED_SHORT, &min_value1, &max_value1 ); + get_type_range( SIGNED_SHORT, &min_value2, &max_value2 ); - ALLOC( *ptr, long_max - long_min + 1 ); + long_min = (long) MIN( min_value1, min_value2 ); + long_max = (long) MAX( max_value1, max_value2 ); + + ALLOC( int_to_real_conversion, long_max - long_min + 1 ); #ifndef NO_DEBUG_ALLOC - (void) unrecord_ptr_alloc_check( *ptr, __FILE__, __LINE__ ); + (void) unrecord_ptr_alloc_check( int_to_real_conversion, + __FILE__, __LINE__ ); #endif - (*ptr) -= long_min; + int_to_real_conversion -= long_min; for_inclusive( i, long_min, long_max ) - (*ptr)[i] = (Real) i; + int_to_real_conversion[i] = (Real) i; } public void get_voxel_values_5d( @@ -270,7 +260,7 @@ step2 -= n3 * step3; step3 -= n4 * step4; - check_real_conversion_lookup( data_type ); + check_real_conversion_lookup(); switch( data_type ) { @@ -286,7 +276,7 @@ { for_less( i4, 0, n4 ) { - *values = unsigned_byte_to_real[ + *values = int_to_real_conversion[ (long) *unsigned_byte_ptr]; ++values; unsigned_byte_ptr += step4; @@ -313,7 +303,7 @@ { for_less( i4, 0, n4 ) { - *values = signed_byte_to_real[ + *values = int_to_real_conversion[ (long) *signed_byte_ptr]; ++values; signed_byte_ptr += step4; @@ -340,7 +330,7 @@ { for_less( i4, 0, n4 ) { - *values = unsigned_short_to_real[ + *values = int_to_real_conversion[ (long) *unsigned_short_ptr]; ++values; unsigned_short_ptr += step4; @@ -367,7 +357,7 @@ { for_less( i4, 0, n4 ) { - *values = signed_short_to_real[ + *values = int_to_real_conversion[ (long) *signed_short_ptr]; ++values; signed_short_ptr += step4; @@ -518,7 +508,7 @@ step0 -= n1 * step1; step1 -= n2 * step2; step2 -= n3 * step3; - check_real_conversion_lookup( data_type ); + check_real_conversion_lookup(); switch( data_type ) { @@ -532,7 +522,7 @@ { for_less( i3, 0, n3 ) { - *values = unsigned_byte_to_real[ + *values = int_to_real_conversion[ (long) *unsigned_byte_ptr]; ++values; unsigned_byte_ptr += step3; @@ -555,7 +545,7 @@ { for_less( i3, 0, n3 ) { - *values = signed_byte_to_real[(long) *signed_byte_ptr]; + *values = int_to_real_conversion[(long) *signed_byte_ptr]; ++values; signed_byte_ptr += step3; } @@ -577,7 +567,7 @@ { for_less( i3, 0, n3 ) { - *values = unsigned_short_to_real[ + *values = int_to_real_conversion[ (long) *unsigned_short_ptr]; ++values; unsigned_short_ptr += step3; @@ -600,7 +590,7 @@ { for_less( i3, 0, n3 ) { - *values = signed_short_to_real[ + *values = int_to_real_conversion[ (long) *signed_short_ptr]; ++values; signed_short_ptr += step3; @@ -730,7 +720,7 @@ step2 = steps[2]; step0 -= n1 * step1; step1 -= n2 * step2; - check_real_conversion_lookup( data_type ); + check_real_conversion_lookup(); switch( data_type ) { @@ -742,7 +732,7 @@ { for_less( i2, 0, n2 ) { - *values = unsigned_byte_to_real[(long) *unsigned_byte_ptr]; + *values = int_to_real_conversion[(long) *unsigned_byte_ptr]; ++values; unsigned_byte_ptr += step2; } @@ -760,7 +750,7 @@ { for_less( i2, 0, n2 ) { - *values = signed_byte_to_real[(long) *signed_byte_ptr]; + *values = int_to_real_conversion[(long) *signed_byte_ptr]; ++values; signed_byte_ptr += step2; } @@ -778,7 +768,7 @@ { for_less( i2, 0, n2 ) { - *values = unsigned_short_to_real[(long) *unsigned_short_ptr]; + *values = int_to_real_conversion[(long) *unsigned_short_ptr]; ++values; unsigned_short_ptr += step2; } @@ -796,7 +786,7 @@ { for_less( i2, 0, n2 ) { - *values = signed_short_to_real[(long) *signed_short_ptr]; + *values = int_to_real_conversion[(long) *signed_short_ptr]; ++values; signed_short_ptr += step2; } @@ -904,7 +894,7 @@ step0 = steps[0]; step1 = steps[1]; step0 -= n1 * step1; - check_real_conversion_lookup( data_type ); + check_real_conversion_lookup(); switch( data_type ) { @@ -914,7 +904,7 @@ { for_less( i1, 0, n1 ) { - *values = unsigned_byte_to_real[(long) *unsigned_byte_ptr]; + *values = int_to_real_conversion[(long) *unsigned_byte_ptr]; ++values; unsigned_byte_ptr += step1; } @@ -928,7 +918,7 @@ { for_less( i1, 0, n1 ) { - *values = signed_byte_to_real[(long) *signed_byte_ptr]; + *values = int_to_real_conversion[(long) *signed_byte_ptr]; ++values; signed_byte_ptr += step1; } @@ -942,7 +932,7 @@ { for_less( i1, 0, n1 ) { - *values = unsigned_short_to_real[(long) *unsigned_short_ptr]; + *values = int_to_real_conversion[(long) *unsigned_short_ptr]; ++values; unsigned_short_ptr += step1; } @@ -956,7 +946,7 @@ { for_less( i1, 0, n1 ) { - *values = signed_short_to_real[(long) *signed_short_ptr]; + *values = int_to_real_conversion[(long) *signed_short_ptr]; ++values; signed_short_ptr += step1; } @@ -1039,7 +1029,7 @@ float *float_ptr; double *double_ptr; - check_real_conversion_lookup( data_type ); + check_real_conversion_lookup(); switch( data_type ) { @@ -1047,7 +1037,7 @@ ASSIGN_PTR(unsigned_byte_ptr) = void_ptr; for_less( i0, 0, n0 ) { - *values = unsigned_byte_to_real[(long) *unsigned_byte_ptr]; + *values = int_to_real_conversion[(long) *unsigned_byte_ptr]; ++values; unsigned_byte_ptr += step0; } @@ -1057,7 +1047,7 @@ ASSIGN_PTR(signed_byte_ptr) = void_ptr; for_less( i0, 0, n0 ) { - *values = signed_byte_to_real[(long) *signed_byte_ptr]; + *values = int_to_real_conversion[(long) *signed_byte_ptr]; ++values; signed_byte_ptr += step0; } @@ -1067,7 +1057,7 @@ ASSIGN_PTR(unsigned_short_ptr) = void_ptr; for_less( i0, 0, n0 ) { - *values = unsigned_short_to_real[(long) *unsigned_short_ptr]; + *values = int_to_real_conversion[(long) *unsigned_short_ptr]; ++values; unsigned_short_ptr += step0; } @@ -1077,7 +1067,7 @@ ASSIGN_PTR(signed_short_ptr) = void_ptr; for_less( i0, 0, n0 ) { - *values = signed_short_to_real[(long) *signed_short_ptr]; + *values = int_to_real_conversion[(long) *signed_short_ptr]; ++values; signed_short_ptr += step0; } @@ -1442,3 +1432,40 @@ get_voxel_values( volume, void_ptr, 1 - dim, &steps[dim], &counts[dim], values ); } + +public void get_volume_voxel_hyperslab( + Volume volume, + int v0, + int v1, + int v2, + int v3, + int v4, + int n0, + int n1, + int n2, + int n3, + int n4, + Real voxels[] ) +{ + switch( get_volume_n_dimensions(volume) ) + { + case 1: + get_volume_voxel_hyperslab_1d( volume, v0, n0, voxels ); + break; + case 2: + get_volume_voxel_hyperslab_2d( volume, v0, v1, n0, n1, voxels ); + break; + case 3: + get_volume_voxel_hyperslab_3d( volume, v0, v1, v2, n0, n1, n2, voxels ); + break; + case 4: + get_volume_voxel_hyperslab_4d( volume, v0, v1, v2, v3, + n0, n1, n2, n3, voxels ); + break; + case 5: + get_volume_voxel_hyperslab_5d( volume, v0, v1, v2, v3, v4, + n0, n1, n2, n3, n4, voxels ); + break; + } +} +
--- a/volume_io/Volumes/input_free.c +++ b/volume_io/Volumes/input_free.c @@ -16,7 +16,7 @@ #include <minc.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/input_free.c,v 1.27 1996-05-17 19:36:20 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/input_free.c,v 1.28 1997-08-13 13:21:32 david Exp $"; #endif #define DEFAULT_SUFFIX "fre" @@ -37,7 +37,7 @@ in volume_input, then the file is converted to bytes on input. This function assumes that volume->filename has been assigned. @CREATED : David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - now uses volume starts ---------------------------------------------------------------------------- */ public Status initialize_free_format_input( @@ -57,7 +57,6 @@ char ch; Real file_separations[MAX_DIMENSIONS]; Real volume_separations[MAX_DIMENSIONS]; - Real origin_voxel[MAX_DIMENSIONS]; Real trans[N_DIMENSIONS]; FILE *file; BOOLEAN axis_valid; @@ -261,11 +260,8 @@ } } - for_less( axis, 0, N_DIMENSIONS ) - origin_voxel[axis] = 0.0; - set_volume_separations( volume, volume_separations ); - set_volume_translation( volume, origin_voxel, trans ); + set_volume_starts( volume, trans ); } set_volume_type( volume, desired_data_type, FALSE, 0.0, 0.0 );
--- a/volume_io/Volumes/input_mnc.c +++ b/volume_io/Volumes/input_mnc.c @@ -16,7 +16,7 @@ #include <minc.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/input_mnc.c,v 1.57 1997-04-17 17:25:47 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/input_mnc.c,v 1.58 1997-08-13 13:21:33 david Exp $"; #endif #define INVALID_AXIS -1 @@ -92,6 +92,9 @@ @CALLS : @CREATED : 1993 David MacDonald @MODIFIED : Nov. 15, 1996 D. MacDonald - added handling of space type +@MODIFIED : May 20, 1997 D. MacDonald - removed float arithmetic in + transformations and made volumes + store starts/steps/dircos ---------------------------------------------------------------------------- */ public Minc_file initialize_minc_input_from_minc_id( @@ -117,16 +120,13 @@ int sizes[MAX_VAR_DIMS]; double file_separations[MAX_VAR_DIMS]; Real volume_separations[MI_NUM_SPACE_DIMS]; + Real volume_starts[MI_NUM_SPACE_DIMS]; Real default_voxel_min, default_voxel_max; - Real world_space[N_DIMENSIONS], voxel_zero; + Real voxel_zero; double start_position[MAX_VAR_DIMS]; double dir_cosines[MAX_VAR_DIMS][MI_NUM_SPACE_DIMS]; double tmp_cosines[MI_NUM_SPACE_DIMS]; BOOLEAN spatial_dim_flags[MAX_VAR_DIMS]; - Vector offset; - Point origin; - Real zero_voxel[MAX_DIMENSIONS]; - Vector spatial_axis; double real_min, real_max; int d, dimvar, which_valid_axis, axis; int spatial_axis_indices[MAX_VAR_DIMS]; @@ -356,67 +356,30 @@ sizes[file->to_volume_index[d]] = (int) file->sizes_in_file[d]; volume_separations[file->to_volume_index[d]] = file_separations[d]; + + if( file->to_volume_index[d] != INVALID_AXIS ) + { + volume_starts[file->to_volume_index[d]] = start_position[d]; + set_volume_direction_unit_cosine( volume, + file->to_volume_index[d], dir_cosines[d] ); + } } } + /* --- create the world transform stored in the volume */ + + set_volume_separations( volume, volume_separations ); + set_volume_starts( volume, volume_starts ); + if( space_type_consensus ) set_volume_space_type( volume, prev_space_type ); /* --- create the file world transform */ - fill_Point( origin, 0.0, 0.0, 0.0 ); - - for_less( d, 0, MAX_DIMENSIONS ) - zero_voxel[d] = 0.0; - - for_less( d, 0, N_DIMENSIONS ) - { - axis = file->spatial_axes[d]; - if( axis != INVALID_AXIS ) - { - fill_Vector( spatial_axis, - dir_cosines[axis][0], - dir_cosines[axis][1], - dir_cosines[axis][2] ); - NORMALIZE_VECTOR( spatial_axis, spatial_axis ); - - SCALE_VECTOR( offset, spatial_axis, start_position[axis] ); - ADD_POINT_VECTOR( origin, origin, offset ); - } - } - - world_space[X] = (Real) Point_x(origin); - world_space[Y] = (Real) Point_y(origin); - world_space[Z] = (Real) Point_z(origin); - compute_world_transform( file->spatial_axes, file_separations, - zero_voxel, world_space, dir_cosines, + dir_cosines, start_position, &file->voxel_to_world_transform ); - /* --- create the world transform stored in the volume */ - - fill_Point( origin, 0.0, 0.0, 0.0 ); - - for_less( d, 0, file->n_file_dimensions ) - { - if( file->to_volume_index[d] != INVALID_AXIS ) - { - set_volume_direction_cosine( volume, - file->to_volume_index[d], - dir_cosines[d] ); - } - } - - general_transform_point( &file->voxel_to_world_transform, - 0.0, 0.0, 0.0, - &world_space[X], &world_space[Y], &world_space[Z]); - - for_less( d, 0, N_DIMENSIONS ) - zero_voxel[d] = 0.0; - - set_volume_translation( volume, zero_voxel, world_space ); - set_volume_separations( volume, volume_separations ); - /* --- decide on type conversion */ if( file->converting_to_colour ) @@ -1114,8 +1077,11 @@ public BOOLEAN advance_input_volume( Minc_file file ) { - int ind, c, axis; - Real voxel[MAX_DIMENSIONS], world_space[N_DIMENSIONS]; + int ind, c, axis; + Real voxel[MAX_DIMENSIONS], world_space[N_DIMENSIONS]; + Real vol_world_space[N_DIMENSIONS]; + Transform offset; + General_transform offset_transform, new_transform; ind = file->n_file_dimensions-1; @@ -1139,6 +1105,8 @@ for_less( ind, 0, get_volume_n_dimensions( file->volume ) ) file->indices[file->valid_file_axes[ind]] = 0; + /*--- update the volume's voxel-to-world transform */ + for_less( c, 0, N_DIMENSIONS ) { axis = file->spatial_axes[c]; @@ -1155,8 +1123,21 @@ for_less( c, 0, get_volume_n_dimensions(file->volume) ) voxel[c] = 0.0; + + convert_voxel_to_world( file->volume, voxel, + &vol_world_space[X], &vol_world_space[Y], + &vol_world_space[Z]); - set_volume_translation( file->volume, voxel, world_space ); + make_identity_transform( &offset ); + for_less( c, 0, N_DIMENSIONS ) + Transform_elem(offset,c,3) = world_space[c] - vol_world_space[c]; + create_linear_transform( &offset_transform, &offset ); + concat_general_transforms( get_voxel_to_world_transform(file->volume), + &offset_transform, &new_transform ); + set_voxel_to_world_transform( file->volume, &new_transform ); + delete_general_transform( &offset_transform ); + + /*--- update the volume if it is cached */ if( file->volume->is_cached_volume ) set_cache_volume_file_offset( &file->volume->cache, file->volume,
--- a/volume_io/Volumes/output_mnc.c +++ b/volume_io/Volumes/output_mnc.c @@ -16,7 +16,7 @@ #include <minc.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_mnc.c,v 1.49 1997-04-17 17:26:02 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_mnc.c,v 1.50 1997-08-13 13:21:34 david Exp $"; #endif #define INVALID_AXIS -1 @@ -84,103 +84,124 @@ @CALLS : @CREATED : Oct. 24, 1995 David MacDonald @MODIFIED : Nov. 15, 1996 D. MacDonald - added handling of space type +@MODIFIED : May. 20, 1997 D. MacDonald - removed float arithmetic +@MODIFIED : May 22, 1997 D. MacDonald - added use_volume_starts_and_steps ---------------------------------------------------------------------------- */ private Status output_world_transform( Minc_file file, STRING space_type, - General_transform *voxel_to_world_transform ) + General_transform *voxel_to_world_transform, + BOOLEAN use_volume_starts_and_steps_flag ) { - double separation[MAX_VAR_DIMS]; + double step[MAX_VAR_DIMS]; Real start[MAX_VAR_DIMS]; - double dir_cosines[MAX_VAR_DIMS][MI_NUM_SPACE_DIMS]; - int i, j, d, axis; - Point origin; - Vector axes[MAX_DIMENSIONS]; - Transform transform, t, inverse; + double dir_cosines[MAX_DIMENSIONS][N_DIMENSIONS]; + int dim, axis, spatial_axes[N_DIMENSIONS]; - if( voxel_to_world_transform == (General_transform *) NULL ) - { - make_identity_transform( &transform ); - } - else if( get_transform_type( voxel_to_world_transform ) == LINEAR ) + /*--- set all starts/steps/dir_cosines to default */ + + for_less( dim, 0, file->n_file_dimensions ) { - transform = *(get_linear_transform_ptr( voxel_to_world_transform )); - } - else - { - print_error( "Cannot output non-linear transforms. Using identity.\n"); - make_identity_transform( &transform ); + start[dim] = 0.0; + step[dim] = 1.0; + dir_cosines[dim][X] = 0.0; + dir_cosines[dim][Y] = 0.0; + dir_cosines[dim][Z] = 0.0; } - get_transform_origin( &transform, &origin ); - get_transform_x_axis( &transform, &axes[X] ); - get_transform_y_axis( &transform, &axes[Y] ); - get_transform_z_axis( &transform, &axes[Z] ); + /*--- if must use the volume's starts and steps */ - separation[X] = MAGNITUDE( axes[X] ); - separation[Y] = MAGNITUDE( axes[Y] ); - separation[Z] = MAGNITUDE( axes[Z] ); - SCALE_VECTOR( axes[X], axes[X], 1.0 / separation[X] ); - SCALE_VECTOR( axes[Y], axes[Y], 1.0 / separation[Y] ); - SCALE_VECTOR( axes[Z], axes[Z], 1.0 / separation[Z] ); + if( use_volume_starts_and_steps_flag ) + { + for_less( dim, 0, file->n_file_dimensions ) + { + if( convert_dim_name_to_spatial_axis( file->dim_names[dim], &axis )) + { + if( file->to_volume_index[dim] == INVALID_AXIS ) + dir_cosines[dim][axis] = 1.0; /*--- default */ + else + { + start[dim] = + file->volume->starts[file->to_volume_index[dim]]; + step[dim] = + file->volume->separations[file->to_volume_index[dim]]; + dir_cosines[dim][X] = file->volume->direction_cosines + [file->to_volume_index[dim]][X]; + dir_cosines[dim][Y] = file->volume->direction_cosines + [file->to_volume_index[dim]][Y]; + dir_cosines[dim][Z] = file->volume->direction_cosines + [file->to_volume_index[dim]][Z]; + } + } + } + } + else /*--- convert the linear transform to starts/steps/dir cosines */ + { + if( voxel_to_world_transform == NULL || + get_transform_type( voxel_to_world_transform ) != LINEAR ) + { + print_error( + "Cannot output null or non-linear transforms. Using identity.\n"); - for_less( i, 0, 3 ) - { - if( Vector_coord(axes[i],i) < 0.0f ) + for_less( dim, 0, file->n_file_dimensions ) + { + if( convert_dim_name_to_spatial_axis( file->dim_names[dim], + &axis )) + dir_cosines[dim][axis] = 1.0; + } + } + else { - SCALE_VECTOR( axes[i], axes[i], -1.0 ); - separation[i] *= -1.0; + spatial_axes[0] = INVALID_AXIS; + spatial_axes[1] = INVALID_AXIS; + spatial_axes[2] = INVALID_AXIS; + + for_less( dim, 0, file->n_file_dimensions ) + { + if( convert_dim_name_to_spatial_axis( file->dim_names[dim], + &axis )) + spatial_axes[axis] = dim; + } + + convert_transform_to_starts_and_steps( voxel_to_world_transform, + file->n_file_dimensions, + NULL, spatial_axes, + start, step, + dir_cosines ); } } - for_less( i, 0, 3 ) + for_less( dim, 0, file->n_file_dimensions ) { - for_less( j, 0, 3 ) - dir_cosines[i][j] = (Real) Vector_coord(axes[i],j); - } - - make_identity_transform( &t ); - set_transform_x_axis( &t, &axes[X] ); - set_transform_y_axis( &t, &axes[Y] ); - set_transform_z_axis( &t, &axes[Z] ); - - (void) compute_transform_inverse( &t, &inverse ); + if( convert_dim_name_to_spatial_axis( file->dim_names[dim], &axis ) ) + { + file->dim_ids[dim] = micreate_std_variable( file->cdfid, + file->dim_names[dim], NC_DOUBLE, 0, NULL); - transform_point( &inverse, - (Real) Point_x(origin), (Real) Point_y(origin), - (Real) Point_z(origin), - &start[X], &start[Y], &start[Z] ); - - for_less( d, 0, file->n_file_dimensions ) - { - if( convert_dim_name_to_spatial_axis( file->dim_names[d], &axis ) ) - { - file->dim_ids[d] = micreate_std_variable( file->cdfid, - file->dim_names[d], NC_DOUBLE, 0, NULL); - - if( file->dim_ids[d] < 0 ) + if( file->dim_ids[dim] < 0 ) return( ERROR ); - (void) miattputdbl( file->cdfid, file->dim_ids[d], MIstep, - separation[axis]); - (void) miattputdbl( file->cdfid, file->dim_ids[d], MIstart, - start[axis]); - if( !is_default_direction_cosine( axis, dir_cosines[axis] ) ) + (void) miattputdbl( file->cdfid, file->dim_ids[dim], MIstep, + step[dim]); + (void) miattputdbl( file->cdfid, file->dim_ids[dim], MIstart, + start[dim]); + if( !is_default_direction_cosine( axis, dir_cosines[dim] ) ) { - (void) ncattput( file->cdfid, file->dim_ids[d], + (void) ncattput( file->cdfid, file->dim_ids[dim], MIdirection_cosines, - NC_DOUBLE, N_DIMENSIONS, dir_cosines[axis]); + NC_DOUBLE, N_DIMENSIONS, dir_cosines[dim]); } - (void) miattputstr( file->cdfid, file->dim_ids[d], MIunits, UNITS ); + (void) miattputstr( file->cdfid, file->dim_ids[dim], MIunits, + UNITS ); if( !equal_strings( space_type, MI_UNKNOWN_SPACE ) ) { - (void) miattputstr( file->cdfid, file->dim_ids[d], MIspacetype, - space_type ); + (void) miattputstr( file->cdfid, file->dim_ids[dim], + MIspacetype, space_type ); } } else - file->dim_ids[d] = -1; + file->dim_ids[dim] = -1; } return( OK ); @@ -375,7 +396,8 @@ } if( output_world_transform( file, volume_to_attach->coordinate_system_name, - voxel_to_world_transform ) != OK ) + voxel_to_world_transform, + options->use_volume_starts_and_steps ) != OK ) { FREE( file ); return( NULL ); @@ -1556,7 +1578,7 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - added use_volume_starts_and_steps ---------------------------------------------------------------------------- */ public void set_default_minc_output_options( @@ -1569,6 +1591,9 @@ options->global_image_range[0] = 0.0; options->global_image_range[1] = -1.0; + + options->use_volume_starts_and_steps = FALSE; + options->use_starts_set = FALSE; } /* ----------------------------- MNI Header ----------------------------------- @@ -1680,3 +1705,28 @@ options->global_image_range[0] = real_min; options->global_image_range[1] = real_max; } + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : set_minc_output_use_volume_starts_and_steps_flag +@INPUT : options + flag +@OUTPUT : +@RETURNS : +@DESCRIPTION: Tells MINC output to use the exact starts and steps stored in + the volume, not the voxel-to-world-transform. This avoids + round-off errors in converting to transform on input, then + from transform on output. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May. 22, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void set_minc_output_use_volume_starts_and_steps_flag( + minc_output_options *options, + BOOLEAN flag ) +{ + options->use_volume_starts_and_steps = flag; + options->use_starts_set = TRUE; +}
--- a/volume_io/Volumes/output_volume.c +++ b/volume_io/Volumes/output_volume.c @@ -15,7 +15,7 @@ #include <internal_volume_io.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_volume.c,v 1.19 1995-11-17 20:25:43 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_volume.c,v 1.20 1997-08-13 13:21:34 david Exp $"; #endif /* ----------------------------- MNI Header ----------------------------------- @@ -292,7 +292,8 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - now sets + use_volume_starts_and_steps flag ---------------------------------------------------------------------------- */ public Status output_modified_volume( @@ -321,7 +322,7 @@ n_dims = get_volume_n_dimensions(volume); - if( options == (minc_output_options *) NULL ) + if( options == NULL ) set_default_minc_output_options( &used_options ); else used_options = *options; @@ -333,6 +334,17 @@ set_minc_output_real_range( &used_options, real_min, real_max ); } + /*--- if the user has not explicitly set the use_volume_starts_and_steps + flag, let's set it if the transform is linear, to output the + same starts as was input, and avoid round-off error */ + + if( !used_options.use_starts_set && + !used_options.use_volume_starts_and_steps && + get_transform_type(get_voxel_to_world_transform(volume)) == LINEAR ) + { + set_minc_output_use_volume_starts_and_steps_flag( &used_options, TRUE ); + } + minc_file = initialize_minc_output( filename, n_dims, dim_names, sizes, file_nc_data_type, file_signed_flag, @@ -340,7 +352,7 @@ get_voxel_to_world_transform(volume), volume, &used_options ); - if( minc_file == (Minc_file) NULL ) + if( minc_file == NULL ) return( ERROR ); status = copy_volume_auxiliary_and_history( minc_file, filename,
--- a/volume_io/Volumes/volumes.c +++ b/volume_io/Volumes/volumes.c @@ -17,7 +17,7 @@ #include <float.h> #ifndef lint -static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/volumes.c,v 1.65 1996-12-09 20:20:27 david Exp $"; +static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/volumes.c,v 1.66 1997-08-13 13:21:33 david Exp $"; #endif STRING XYZ_dimension_names[] = { MIxspace, MIyspace, MIzspace }; @@ -136,6 +136,7 @@ @CALLS : @CREATED : June, 1993 David MacDonald @MODIFIED : Nov. 15, 1996 D. MacDonald - handles space type +@MODIFIED : May 22, 1996 D. MacDonald - now stores starts/steps ---------------------------------------------------------------------------- */ public Volume create_volume( @@ -177,14 +178,11 @@ volume->real_value_translation = 0.0; for_less( i, 0, N_DIMENSIONS ) - { - volume->world_space_for_translation_voxel[i] = 0.0; volume->spatial_axes[i] = -1; - } for_less( i, 0, n_dimensions ) { - volume->translation_voxel[i] = 0.0; + volume->starts[i] = 0.0; volume->separations[i] = 1.0; volume->direction_cosines[i][X] = 0.0; volume->direction_cosines[i][Y] = 0.0; @@ -213,6 +211,7 @@ make_identity_transform( &identity ); create_linear_transform( &volume->voxel_to_world_transform, &identity ); + volume->voxel_to_world_transform_uptodate = TRUE; volume->coordinate_system_name = create_string( MI_UNKNOWN_SPACE ); @@ -571,76 +570,122 @@ get_volume_sizes( volume, sizes ); - for_less( i, 0, get_multidim_n_dimensions( &volume->array ) ) + for_less( i, 0, get_volume_n_dimensions( volume ) ) n *= (unsigned int) sizes[i]; return( n ); } /* ----------------------------- MNI Header ----------------------------------- -@NAME : set_voxel_to_world_transform +@NAME : assign_voxel_to_world_transform @INPUT : volume transform @OUTPUT : @RETURNS : -@DESCRIPTION: Sets the volume's transformation from voxel to world coords. +@DESCRIPTION: Updates the volume's transformation from voxel to world coords. @METHOD : @GLOBALS : @CALLS : -@CREATED : 1993 David MacDonald +@CREATED : May 20, 1997 D. MacDonald - created from + set_voxel_to_world_transform @MODIFIED : ---------------------------------------------------------------------------- */ -public void set_voxel_to_world_transform( +private void assign_voxel_to_world_transform( Volume volume, General_transform *transform ) { - Real sign; - int c, axis; - Vector axes[N_DIMENSIONS]; - Transform *linear_transform; - delete_general_transform( &volume->voxel_to_world_transform ); volume->voxel_to_world_transform = *transform; - - if( get_transform_type( transform ) == LINEAR ) - { - linear_transform = get_linear_transform_ptr( transform ); - get_transform_x_axis( linear_transform, &axes[X] ); - get_transform_y_axis( linear_transform, &axes[Y] ); - get_transform_z_axis( linear_transform, &axes[Z] ); +} - for_less( c, 0, N_DIMENSIONS ) - { - axis = volume->spatial_axes[c]; - if( volume->separations[axis] < 0.0 ) - sign = -1.0; - else - sign = 1.0; - if( axis >= 0 ) - volume->separations[axis] = sign * MAGNITUDE( axes[c] ); - } - } +/* ----------------------------- MNI Header ----------------------------------- +@NAME : dot_vectors +@INPUT : n + v1 + v2 +@OUTPUT : +@RETURNS : Dot product +@DESCRIPTION: Computes the dot product of 2 n-dimensional vectors. + This function should be moved to some vector routines. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : 1993 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +private Real dot_vectors( + int n, + Real v1[], + Real v2[] ) +{ + int i; + Real d; + + d = 0.0; + for_less( i, 0, n ) + d += v1[i] * v2[i]; + + return( d ); } /* ----------------------------- MNI Header ----------------------------------- -@NAME : get_voxel_to_world_transform -@INPUT : -@OUTPUT : -@RETURNS : transform -@DESCRIPTION: Returns a pointer to the voxel to world transform of the volume. -@METHOD : -@GLOBALS : -@CALLS : -@CREATED : 1993 David MacDonald -@MODIFIED : +@NAME : cross_3D_vector +@INPUT : v1 + v2 +@OUTPUT : cross +@RETURNS : +@DESCRIPTION: Computes the cross product of 2 n-dimensional vectors. + This function should be moved to some vector routines. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : 1993 David MacDonald +@MODIFIED : ---------------------------------------------------------------------------- */ -public General_transform *get_voxel_to_world_transform( - Volume volume ) +private void cross_3D_vector( + Real v1[], + Real v2[], + Real cross[] ) { - return( &volume->voxel_to_world_transform ); + cross[0] = v1[1] * v2[2] - v1[2] * v2[1]; + cross[1] = v1[2] * v2[0] - v1[0] * v2[2]; + cross[2] = v1[0] * v2[1] - v1[1] * v2[0]; +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : normalize_vector +@INPUT : n + v1 +@OUTPUT : v1_normalized +@RETURNS : +@DESCRIPTION: Normalizes the length of v1 to 1, placing result in + v1_normalized +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 22, 1997 D. MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +private void normalize_vector( + Real v1[], + Real v1_normalized[] ) +{ + int d; + Real mag; + + mag = dot_vectors( N_DIMENSIONS, v1, v1 ); + if( mag <= 0.0 ) + mag = 1.0; + + mag = sqrt( mag ); + + for_less( d, 0, N_DIMENSIONS ) + v1_normalized[d] = v1[d] / mag; } /* ----------------------------- MNI Header ----------------------------------- @@ -660,24 +705,26 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - now uses starts/steps ---------------------------------------------------------------------------- */ public void compute_world_transform( int spatial_axes[N_DIMENSIONS], Real separations[], - Real translation_voxel[], - Real world_space_for_translation_voxel[N_DIMENSIONS], Real direction_cosines[][N_DIMENSIONS], + Real starts[], General_transform *world_transform ) { Transform transform; Real separations_3D[N_DIMENSIONS]; - Vector directions[N_DIMENSIONS]; - Point point; - Real x_trans, y_trans, z_trans; - Real voxel[N_DIMENSIONS]; - int c, axis, n_axes, axis_list[N_DIMENSIONS]; + Real directions[N_DIMENSIONS][N_DIMENSIONS]; + Real starts_3D[N_DIMENSIONS]; + Real normal[N_DIMENSIONS]; + int dim, c, a1, a2, axis, n_axes; + int axis_list[N_DIMENSIONS]; + + /*--- find how many direction cosines are specified, and set the + 3d separations and starts */ n_axes = 0; @@ -687,19 +734,17 @@ if( axis >= 0 ) { separations_3D[c] = separations[axis]; - voxel[c] = translation_voxel[axis]; - fill_Vector( directions[c], - direction_cosines[axis][X], - direction_cosines[axis][Y], - direction_cosines[axis][Z]); - + starts_3D[c] = starts[axis]; + directions[c][X] = direction_cosines[axis][X]; + directions[c][Y] = direction_cosines[axis][Y]; + directions[c][Z] = direction_cosines[axis][Z]; axis_list[n_axes] = c; ++n_axes; } else { separations_3D[c] = 1.0; - voxel[c] = 0.0; + starts_3D[c] = 0.0; } } @@ -709,55 +754,81 @@ return; } + /*--- convert 1 or 2 axes to 3 axes */ + if( n_axes == 1 ) { - create_two_orthogonal_vectors( &directions[axis_list[0]], - &directions[(axis_list[0]+1) % N_DIMENSIONS], - &directions[(axis_list[0]+2) % N_DIMENSIONS] ); + a1 = (axis_list[0] + 1) % N_DIMENSIONS; + a2 = (axis_list[0] + 2) % N_DIMENSIONS; + + /*--- create an orthogonal vector */ + + directions[a1][X] = directions[axis_list[0]][Y] + + directions[axis_list[0]][Z]; + directions[a1][Y] = -directions[axis_list[0]][X] - + directions[axis_list[0]][Z]; + directions[a1][Z] = directions[axis_list[0]][Y] - + directions[axis_list[0]][X]; + + cross_3D_vector( directions[axis_list[0]], directions[a1], + directions[a2] ); + normalize_vector( directions[a1], directions[a1] ); + normalize_vector( directions[a2], directions[a2] ); } else if( n_axes == 2 ) { - axis = N_DIMENSIONS - axis_list[0] - axis_list[1]; + a2 = N_DIMENSIONS - axis_list[0] - axis_list[1]; - CROSS_VECTORS( directions[axis], - directions[(axis+1) % N_DIMENSIONS], - directions[(axis+2) % N_DIMENSIONS] ); + cross_3D_vector( directions[axis_list[0]], directions[axis_list[1]], + directions[a2] ); + + normalize_vector( directions[a2], directions[a2] ); } - if( EQUAL_VECTORS( directions[0], directions[1] ) || - EQUAL_VECTORS( directions[1], directions[2] ) || - EQUAL_VECTORS( directions[2], directions[0] ) ) + /*--- check to make sure that 3 axes are not a singular system */ + + for_less( dim, 0, N_DIMENSIONS ) + { + cross_3D_vector( directions[dim], directions[(dim+1)%N_DIMENSIONS], + normal ); + if( normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0 ) + break; + } + + if( dim < N_DIMENSIONS ) { - fill_Vector( directions[0], 1.0, 0.0, 0.0 ); - fill_Vector( directions[1], 0.0, 1.0, 0.0 ); - fill_Vector( directions[2], 0.0, 0.0, 1.0 ); + directions[0][0] = 1.0; + directions[0][1] = 0.0; + directions[0][2] = 0.0; + directions[1][0] = 0.0; + directions[1][1] = 1.0; + directions[1][2] = 0.0; + directions[2][0] = 0.0; + directions[2][1] = 0.0; + directions[2][2] = 1.0; } + /*--- make the linear transformation */ + + make_identity_transform( &transform ); + for_less( c, 0, N_DIMENSIONS ) { - NORMALIZE_VECTOR( directions[c], directions[c] ); - SCALE_VECTOR( directions[c], directions[c], separations_3D[c] ); - } - - fill_Point( point, 0.0, 0.0, 0.0 ); - - make_change_to_bases_transform( &point, - &directions[X], &directions[Y], &directions[Z], &transform ); + for_less( dim, 0, N_DIMENSIONS ) + { + Transform_elem(transform,dim,c) = directions[c][dim] * + separations_3D[c]; - transform_point( &transform, voxel[X], voxel[Y], voxel[Z], - &x_trans, &y_trans, &z_trans ); - - fill_Point( point, world_space_for_translation_voxel[X] - x_trans, - world_space_for_translation_voxel[Y] - y_trans, - world_space_for_translation_voxel[Z] - z_trans ); - - set_transform_origin( &transform, &point ); + Transform_elem(transform,dim,3) += directions[c][dim] * + starts_3D[c]; + } + } create_linear_transform( world_transform, &transform ); } /* ----------------------------- MNI Header ----------------------------------- -@NAME : recompute_world_transform +@NAME : check_recompute_world_transform @INPUT : volume @OUTPUT : @RETURNS : @@ -767,24 +838,277 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 20, 1997 D. MacDonald - now checks update flag ---------------------------------------------------------------------------- */ -private void recompute_world_transform( +private void check_recompute_world_transform( Volume volume ) { General_transform world_transform; - Real separations[MAX_DIMENSIONS]; + + if( !volume->voxel_to_world_transform_uptodate ) + { + volume->voxel_to_world_transform_uptodate = TRUE; + + compute_world_transform( volume->spatial_axes, + volume->separations, + volume->direction_cosines, + volume->starts, + &world_transform ); + + assign_voxel_to_world_transform( volume, &world_transform ); + } +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : convert_transform_origin_to_starts +@INPUT : origin + n_volume_dimensions + spatial_axes + dir_cosines +@OUTPUT : starts +@RETURNS : +@DESCRIPTION: Converts a transform origin into starts (multiples of the + dir_cosines). dir_cosines need not be mutually orthogonal +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May. 22, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +private void convert_transform_origin_to_starts( + Real origin[], + int n_volume_dimensions, + int spatial_axes[], + Real dir_cosines[][N_DIMENSIONS], + Real starts[] ) +{ + int axis, dim, which[N_DIMENSIONS], n_axes, i, j; + Real o_dot_c, c_dot_c; + Real x_dot_x, x_dot_y, x_dot_v, y_dot_y, y_dot_v, bottom; + Real **matrix, solution[N_DIMENSIONS]; + + for_less( dim, 0, n_volume_dimensions ) + starts[dim] = 0.0; + + /*--- get the list of valid axes (which) */ - get_volume_separations( volume, separations ); + n_axes = 0; + for_less( dim, 0, N_DIMENSIONS ) + { + axis = spatial_axes[dim]; + if( axis >= 0 ) + { + which[n_axes] = axis; + ++n_axes; + } + } + + /*--- get the starts: computed differently for 1, 2, or 3 axes */ + + if( n_axes == 1 ) + { + o_dot_c = dot_vectors( N_DIMENSIONS, origin, dir_cosines[which[0]] ); + c_dot_c = dot_vectors( N_DIMENSIONS, dir_cosines[which[0]], + dir_cosines[which[0]] ); + + if( c_dot_c != 0.0 ) + starts[which[0]] = o_dot_c / c_dot_c; + } + else if( n_axes == 2 ) + { + x_dot_x = dot_vectors( N_DIMENSIONS, dir_cosines[which[0]], + dir_cosines[which[0]] ); + x_dot_v = dot_vectors( N_DIMENSIONS, dir_cosines[which[0]], origin ); + x_dot_y = dot_vectors( N_DIMENSIONS, dir_cosines[which[0]], + dir_cosines[which[1]] ); + y_dot_y = dot_vectors( N_DIMENSIONS, dir_cosines[which[1]], + dir_cosines[which[1]] ); + y_dot_v = dot_vectors( N_DIMENSIONS, dir_cosines[which[1]], origin ); + + bottom = x_dot_x * y_dot_y - x_dot_y * x_dot_y; + + if( bottom != 0.0 ) + { + starts[which[0]] = (x_dot_v * y_dot_y - x_dot_y * y_dot_v) / bottom; + starts[which[1]] = (y_dot_v * x_dot_x - x_dot_y * x_dot_v) / bottom; + } + } + else if( n_axes == 3 ) + { + /*--- this is the usual case, solve the equations to find what + starts give the desired origin */ + + ALLOC2D( matrix, N_DIMENSIONS, N_DIMENSIONS ); + + for_less( i, 0, N_DIMENSIONS ) + for_less( j, 0, N_DIMENSIONS ) + { + matrix[i][j] = dir_cosines[which[j]][i]; + } + + if( solve_linear_system( N_DIMENSIONS, matrix, origin, solution ) ) + { + starts[which[0]] = solution[0]; + starts[which[1]] = solution[1]; + starts[which[2]] = solution[2]; + } + + FREE2D( matrix ); + } + else + { + print_error( + "Invalid number of axes in convert_transform_origin_to_starts\n"); + } +} - compute_world_transform( volume->spatial_axes, separations, - volume->translation_voxel, - volume->world_space_for_translation_voxel, - volume->direction_cosines, - &world_transform ); +/* ----------------------------- MNI Header ----------------------------------- +@NAME : convert_transform_to_starts_and_steps +@INPUT : transform + separation_signs +@OUTPUT : starts + steps + dir_cosines + spatial_axes +@RETURNS : +@DESCRIPTION: Converts a linear transform to a set of 3 starts, 3 steps, + and 3 direction cosines. The separation signs determine + the desired signs of each of the separations. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May. 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ + +public void convert_transform_to_starts_and_steps( + General_transform *transform, + int n_volume_dimensions, + Real step_signs[], + int spatial_axes[], + Real starts[], + Real steps[], + Real dir_cosines[][N_DIMENSIONS] ) +{ + Real sign, mag; + int axis, dim; + Real axes[N_DIMENSIONS][N_DIMENSIONS]; + Real origin[N_DIMENSIONS]; + Transform *linear_transform; + + if( get_transform_type( transform ) != LINEAR ) + { + print_error( "convert_transform_to_starts_and_steps(): non-linear transform found.\n" ); + return; + } + + linear_transform = get_linear_transform_ptr( transform ); + + get_transform_origin_real( linear_transform, origin ); + get_transform_x_axis_real( linear_transform, &axes[X][0] ); + get_transform_y_axis_real( linear_transform, &axes[Y][0] ); + get_transform_z_axis_real( linear_transform, &axes[Z][0] ); + + /*--- assign default steps */ + + for_less( dim, 0, n_volume_dimensions ) + steps[dim] = 1.0; + + /*--- assign the steps and dir_cosines for the spatial axes */ + + for_less( dim, 0, N_DIMENSIONS ) + { + axis = spatial_axes[dim]; + if( axis >= 0 ) + { + mag = dot_vectors( N_DIMENSIONS, axes[dim], axes[dim] ); + + if( mag <= 0.0 ) + mag = 1.0; + mag = sqrt( mag ); - set_voxel_to_world_transform( volume, &world_transform ); + if( step_signs == NULL ) + { + if( axes[dim][dim] < 0.0 ) + sign = -1.0; + else + sign = 1.0; + } + else /*--- make the sign of steps match the step_signs passed in */ + { + if( step_signs[axis] < 0.0 ) + sign = -1.0; + else + sign = 1.0; + } + + steps[axis] = sign * mag; + dir_cosines[axis][X] = axes[dim][X] / steps[axis]; + dir_cosines[axis][Y] = axes[dim][Y] / steps[axis]; + dir_cosines[axis][Z] = axes[dim][Z] / steps[axis]; + } + } + + /*--- finally, get the starts */ + + convert_transform_origin_to_starts( origin, n_volume_dimensions, + spatial_axes, dir_cosines, starts ); +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : set_voxel_to_world_transform +@INPUT : volume + transform +@OUTPUT : +@RETURNS : +@DESCRIPTION: Sets the volume's transformation from voxel to world coords. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : 1993 David MacDonald +@MODIFIED : May 22, 1997 D. MacDonald - recomputes the starts/steps +---------------------------------------------------------------------------- */ + +public void set_voxel_to_world_transform( + Volume volume, + General_transform *transform ) +{ + assign_voxel_to_world_transform( volume, transform ); + volume->voxel_to_world_transform_uptodate = TRUE; + + if( get_transform_type( transform ) == LINEAR ) + { + convert_transform_to_starts_and_steps( transform, + get_volume_n_dimensions(volume), + volume->separations, + volume->spatial_axes, + volume->starts, + volume->separations, + volume->direction_cosines ); + } +} + +/* ----------------------------- MNI Header ----------------------------------- +@NAME : get_voxel_to_world_transform +@INPUT : +@OUTPUT : +@RETURNS : transform +@DESCRIPTION: Returns a pointer to the voxel to world transform of the volume. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : 1993 David MacDonald +@MODIFIED : May 22, 1997 D. MacDonald - now delays recomputing transform +---------------------------------------------------------------------------- */ + +public General_transform *get_voxel_to_world_transform( + Volume volume ) +{ + check_recompute_world_transform( volume ); + + return( &volume->voxel_to_world_transform ); } /* ----------------------------- MNI Header ----------------------------------- @@ -917,7 +1241,7 @@ { int i; - for_less( i, 0, get_multidim_n_dimensions( &volume->array ) ) + for_less( i, 0, get_volume_n_dimensions( volume ) ) separations[i] = volume->separations[i]; } @@ -932,7 +1256,7 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - now delays recomputing transform ---------------------------------------------------------------------------- */ public void set_volume_separations( @@ -941,74 +1265,108 @@ { int i; - for_less( i, 0, get_multidim_n_dimensions( &volume->array ) ) + for_less( i, 0, get_volume_n_dimensions( volume ) ) volume->separations[i] = separations[i]; - recompute_world_transform( volume ); + volume->voxel_to_world_transform_uptodate = FALSE; } /* ----------------------------- MNI Header ----------------------------------- -@NAME : set_volume_translation +@NAME : set_volume_starts @INPUT : volume - voxel - world_space_voxel_maps_to + starts[] @OUTPUT : @RETURNS : @DESCRIPTION: Sets the translation portion of the voxel to world transform, - by specifying a point in voxel space (voxel), and a point - in world space (world_space_voxel_maps_to). + by specifying the start vector, as specified by the MINC format. @METHOD : @GLOBALS : @CALLS : -@CREATED : 1993 David MacDonald -@MODIFIED : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : May 22, 1997 D. MacDonald - now delays recomputing transform ---------------------------------------------------------------------------- */ -public void set_volume_translation( +public void set_volume_starts( Volume volume, - Real voxel[], - Real world_space_voxel_maps_to[] ) + Real starts[] ) { int c; - for_less( c, 0, get_multidim_n_dimensions( &volume->array ) ) - volume->translation_voxel[c] = voxel[c]; + for_less( c, 0, get_volume_n_dimensions( volume ) ) + volume->starts[c] = starts[c]; + + volume->voxel_to_world_transform_uptodate = FALSE; +} - for_less( c, 0, N_DIMENSIONS ) - volume->world_space_for_translation_voxel[c] = - world_space_voxel_maps_to[c]; +/* ----------------------------- MNI Header ----------------------------------- +@NAME : get_volume_starts +@INPUT : volume +@OUTPUT : starts +@RETURNS : +@DESCRIPTION: Passes back the start vector of the volume. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : May 20, 1997 David MacDonald +@MODIFIED : +---------------------------------------------------------------------------- */ - recompute_world_transform( volume ); +public void get_volume_starts( + Volume volume, + Real starts[] ) +{ + int c; + + for_less( c, 0, get_volume_n_dimensions( volume ) ) + starts[c] = volume->starts[c]; } /* ----------------------------- MNI Header ----------------------------------- -@NAME : get_volume_translation +@NAME : set_volume_direction_unit_cosine @INPUT : volume -@OUTPUT : voxel - world_space_voxel_maps_to + axis + dir +@OUTPUT : @RETURNS : -@DESCRIPTION: Passes back the translation portion of the voxel to world - transform of the volume. +@DESCRIPTION: Sets the direction cosine for one axis, assumed to be unit + length. @METHOD : @GLOBALS : @CALLS : -@CREATED : 1993 David MacDonald -@MODIFIED : +@CREATED : May 20, 1997 David MacDonald ---------------------------------------------------------------------------- */ -public void get_volume_translation( - Volume volume, - Real voxel[], - Real world_space_voxel_maps_to[] ) +public void set_volume_direction_unit_cosine( + Volume volume, + int axis, + Real dir[] ) { - int c; + int dim; + + if( axis < 0 || axis >= get_volume_n_dimensions(volume) ) + { + print_error( + "set_volume_direction_cosine: cannot set dir cosine for axis %d\n", + axis ); + return; + } - for_less( c, 0, get_multidim_n_dimensions( &volume->array ) ) - voxel[c] = volume->translation_voxel[c]; + /*--- check if this is a spatial axis */ + + for_less( dim, 0, N_DIMENSIONS ) + { + if( volume->spatial_axes[dim] == axis ) + break; + } - for_less( c, 0, N_DIMENSIONS ) - world_space_voxel_maps_to[c] = - volume->world_space_for_translation_voxel[c]; + if( dim == N_DIMENSIONS ) /* this is not a spatial axis, ignore the dir */ + return; + + volume->direction_cosines[axis][X] = dir[X]; + volume->direction_cosines[axis][Y] = dir[Y]; + volume->direction_cosines[axis][Z] = dir[Z]; + + volume->voxel_to_world_transform_uptodate = FALSE; } /* ----------------------------- MNI Header ----------------------------------- @@ -1023,7 +1381,8 @@ @GLOBALS : @CALLS : @CREATED : 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 20, 1997 D. MacDonald - split into + set_volume_direction_unit_cosine ---------------------------------------------------------------------------- */ public void set_volume_direction_cosine( @@ -1031,29 +1390,7 @@ int axis, Real dir[] ) { - int d; - Real len; - - if( axis < 0 || axis >= get_volume_n_dimensions(volume) ) - { - print_error( - "set_volume_direction_cosine: cannot set dir cosine for axis %d\n", - axis ); - return; - } - - for_less( d, 0, N_DIMENSIONS ) - { - if( volume->spatial_axes[d] == axis ) - break; - } - - if( d == N_DIMENSIONS ) /* this is not a spatial axis, ignore the dir */ - return; - - volume->direction_cosines[axis][X] = dir[X]; - volume->direction_cosines[axis][Y] = dir[Y]; - volume->direction_cosines[axis][Z] = dir[Z]; + Real len, unit_vector[N_DIMENSIONS]; len = dir[X] * dir[X] + dir[Y] * dir[Y] + dir[Z] * dir[Z]; @@ -1063,15 +1400,16 @@ return; } - if( len != 1.0 ) - { - len = sqrt( len ); - volume->direction_cosines[axis][X] /= len; - volume->direction_cosines[axis][Y] /= len; - volume->direction_cosines[axis][Z] /= len; - } + if( len <= 0.0 ) + len = 1.0; + + len = sqrt( len ); - recompute_world_transform( volume ); + unit_vector[X] = dir[X] / len; + unit_vector[Y] = dir[Y] / len; + unit_vector[Z] = dir[Z] / len; + + set_volume_direction_unit_cosine( volume, axis, unit_vector ); } /* ----------------------------- MNI Header ----------------------------------- @@ -1111,11 +1449,17 @@ } if( d == N_DIMENSIONS ) /* this is not a spatial axis, ignore the dir */ - return; - - dir[X] = volume->direction_cosines[axis][X]; - dir[Y] = volume->direction_cosines[axis][Y]; - dir[Z] = volume->direction_cosines[axis][Z]; + { + dir[X] = 0.0; + dir[Y] = 0.0; + dir[Z] = 0.0; + } + else + { + dir[X] = volume->direction_cosines[axis][X]; + dir[Y] = volume->direction_cosines[axis][Y]; + dir[Z] = volume->direction_cosines[axis][Z]; + } } /* ----------------------------- MNI Header ----------------------------------- @@ -1198,7 +1542,7 @@ Note that centre of first voxel corresponds to (0.0,0.0,0.0) in voxel coordinates. @CREATED : Mar 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - checks to recompute transform ---------------------------------------------------------------------------- */ public void convert_voxel_to_world( @@ -1210,6 +1554,8 @@ { Real xyz[N_DIMENSIONS]; + check_recompute_world_transform( volume ); + reorder_voxel_to_xyz( volume, voxel, xyz ); /* apply linear transform */ @@ -1275,7 +1621,7 @@ vector is a normal vector (ie. a derivative), so transforms by transpose of inverse transform. @CREATED : Mar 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - checks to recompute transform ---------------------------------------------------------------------------- */ public void convert_voxel_normal_vector_to_world( @@ -1288,7 +1634,7 @@ Real xyz[N_DIMENSIONS]; Transform *inverse; - reorder_voxel_to_xyz( volume, voxel_vector, xyz ); + check_recompute_world_transform( volume ); if( get_transform_type( &volume->voxel_to_world_transform ) != LINEAR ) handle_internal_error( "Cannot get normal vector of nonlinear xforms."); @@ -1298,6 +1644,8 @@ /* transform vector by transpose of inverse transformation */ + reorder_voxel_to_xyz( volume, voxel_vector, xyz ); + *x_world = Transform_elem(*inverse,0,0) * xyz[X] + Transform_elem(*inverse,1,0) * xyz[Y] + Transform_elem(*inverse,2,0) * xyz[Z]; @@ -1386,7 +1734,7 @@ @RETURNS : @DESCRIPTION: Converts from world coordinates to voxel coordinates. @CREATED : Mar 1993 David MacDonald -@MODIFIED : +@MODIFIED : May 22, 1997 D. MacDonald - checks to recompute transform ---------------------------------------------------------------------------- */ public void convert_world_to_voxel( @@ -1398,7 +1746,7 @@ { Real xyz[N_DIMENSIONS]; - /* apply linear transform */ + check_recompute_world_transform( volume ); general_inverse_transform_point( &volume->voxel_to_world_transform, x_world, y_world, z_world, @@ -1738,11 +2086,9 @@ { int c, sizes[MAX_DIMENSIONS]; Real separations[MAX_DIMENSIONS]; + Real starts[MAX_DIMENSIONS]; + Real dir_cosine[N_DIMENSIONS]; Volume copy; - General_transform transform; - - get_volume_sizes( volume, sizes ); - get_volume_separations( volume, separations ); if( nc_data_type == NC_UNSPECIFIED ) { @@ -1754,7 +2100,6 @@ copy = create_volume( get_volume_n_dimensions(volume), volume->dimension_names, nc_data_type, signed_flag, voxel_min, voxel_max ); - set_volume_sizes( copy, sizes ); for_less( c, 0, N_DIMENSIONS ) copy->spatial_axes[c] = volume->spatial_axes[c]; @@ -1763,15 +2108,20 @@ get_volume_real_min(volume), get_volume_real_max(volume) ); + get_volume_sizes( volume, sizes ); + set_volume_sizes( copy, sizes ); + + get_volume_separations( volume, separations ); set_volume_separations( copy, separations ); + get_volume_starts( volume, starts ); + set_volume_starts( copy, starts ); + for_less( c, 0, get_volume_n_dimensions(volume) ) - set_volume_direction_cosine( copy, c, volume->direction_cosines[c] ); - - copy_general_transform( get_voxel_to_world_transform(volume), - &transform ); - - set_voxel_to_world_transform( copy, &transform ); + { + get_volume_direction_cosine( volume, c, dir_cosine ); + set_volume_direction_unit_cosine( copy, c, dir_cosine ); + } set_volume_space_type( copy, volume->coordinate_system_name );