Mercurial > hg > minc-tools
changeset 62:40f6b15f786e
Initial revision
author | david <david> |
---|---|
date | Mon, 11 Jan 1993 11:30:47 +0000 |
parents | 4bdb79d75e0a |
children | 637853a7276a |
files | volume_io/Geometry/points.c volume_io/Geometry/transforms.c |
diffstat | 2 files changed, 752 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/volume_io/Geometry/points.c @@ -0,0 +1,149 @@ + +#include <def_mni.h> + +public Boolean null_Point( + Point *point ) +{ + return( Point_x(*point) == 0.0 && + Point_y(*point) == 0.0 && + Point_z(*point) == 0.0 ); +} + +public Boolean null_Vector( + Vector *vector ) +{ + return( Vector_x(*vector) == 0.0 && + Vector_y(*vector) == 0.0 && + Vector_z(*vector) == 0.0 ); +} + +public Real distance_between_points( + Point *p1, + Point *p2 ) +{ + Vector diff; + + SUB_POINTS( diff, *p1, *p2 ); + + return( MAGNITUDE( diff ) ); +} + +public Boolean points_within_distance( + Point *p1, + Point *p2, + Real distance ) +{ + Real dx; + + dx = Point_x(*p1) - Point_x(*p2); + if( dx < 0.0 ) dx = -dx; + + if( dx > distance ) + return( FALSE ); + + return( distance_between_points( p1, p2 ) <= distance ); +} + +public void create_noncolinear_vector( + Vector *v, + Vector *not_v ) +{ + Real x, y, z; + + x = Vector_x(*v); + y = Vector_y(*v); + z = Vector_z(*v); + + if( x != 0.0 || y != 0.0 ) + { + fill_Vector( *not_v, x, y, z + 1.0 ); + } + else + { + fill_Vector( *not_v, 1.0, y, z ); + } +} + +public void apply_point_to_min_and_max( + Point *point, + Point *min_point, + Point *max_point ) +{ + int c; + + for_less( c, 0, N_DIMENSIONS ) + { + if( Point_coord(*point,c) < Point_coord(*min_point,c) ) + { + Point_coord(*min_point,c) = Point_coord(*point,c); + } + if( Point_coord(*point,c) > Point_coord(*max_point,c) ) + { + Point_coord(*max_point,c) = Point_coord(*point,c); + } + } +} + +public void expand_min_and_max_points( + Point *min_point, + Point *max_point, + Point *min_to_check, + Point *max_to_check ) +{ + int c; + + for_less( c, 0, N_DIMENSIONS ) + { + if( Point_coord(*min_to_check,c) < Point_coord(*min_point,c) ) + { + Point_coord(*min_point,c) = Point_coord(*min_to_check,c); + } + if( Point_coord(*max_to_check,c) > Point_coord(*max_point,c) ) + { + Point_coord(*max_point,c) = Point_coord(*max_to_check,c); + } + } +} + +public void get_range_points( + int n_points, + Point points[], + Point *min_corner, + Point *max_corner ) +{ + int i; + + if( n_points > 0 ) + { + *min_corner = points[0]; + *max_corner = points[0]; + + for_less( i, 1, n_points ) + { + apply_point_to_min_and_max( &points[i], min_corner, max_corner ); + } + } +} + +public void get_points_centroid( + int n_points, + Point points[], + Point *centroid ) +{ + int i; + + Point_x(*centroid) = 0.0; + Point_y(*centroid) = 0.0; + Point_z(*centroid) = 0.0; + + for_less( i, 0, n_points ) + { + Point_x(*centroid) += Point_x(points[i]); + Point_y(*centroid) += Point_y(points[i]); + Point_z(*centroid) += Point_z(points[i]); + } + + Point_x(*centroid) /= (Real) n_points; + Point_y(*centroid) /= (Real) n_points; + Point_z(*centroid) /= (Real) n_points; +}
new file mode 100644 --- /dev/null +++ b/volume_io/Geometry/transforms.c @@ -0,0 +1,603 @@ + +#include <def_mni.h> + +public void make_identity_transform( Transform *transform ) +{ + Transform_elem( *transform, 0, 0 ) = 1.0; + Transform_elem( *transform, 0, 1 ) = 0.0; + Transform_elem( *transform, 0, 2 ) = 0.0; + Transform_elem( *transform, 0, 3 ) = 0.0; + Transform_elem( *transform, 1, 0 ) = 0.0; + Transform_elem( *transform, 1, 1 ) = 1.0; + Transform_elem( *transform, 1, 2 ) = 0.0; + Transform_elem( *transform, 1, 3 ) = 0.0; + Transform_elem( *transform, 2, 0 ) = 0.0; + Transform_elem( *transform, 2, 1 ) = 0.0; + Transform_elem( *transform, 2, 2 ) = 1.0; + Transform_elem( *transform, 2, 3 ) = 0.0; + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void make_scale_transform( + Real sx, + Real sy, + Real sz, + Transform *transform ) +{ + make_identity_transform( transform ); + + Transform_elem( *transform, 0, 0 ) = sx; + Transform_elem( *transform, 1, 1 ) = sy; + Transform_elem( *transform, 2, 2 ) = sz; +} + +public void get_transform_origin( + Transform *transform, + Point *origin ) +{ + fill_Point( *origin, + Transform_elem(*transform,0,3), + Transform_elem(*transform,1,3), + Transform_elem(*transform,2,3) ); +} + +public void set_transform_origin( + Transform *transform, + Point *origin ) +{ + Transform_elem(*transform,0,3) = Point_x(*origin); + Transform_elem(*transform,1,3) = Point_y(*origin); + Transform_elem(*transform,2,3) = Point_z(*origin); +} + +public void get_transform_x_axis( + Transform *transform, + Vector *x_axis ) +{ + fill_Vector( *x_axis, + Transform_elem(*transform,0,0), + Transform_elem(*transform,1,0), + Transform_elem(*transform,2,0) ); +} + +public void set_transform_x_axis( + Transform *transform, + Vector *x_axis ) +{ + Transform_elem(*transform,0,0) = Vector_x(*x_axis); + Transform_elem(*transform,1,0) = Vector_y(*x_axis); + Transform_elem(*transform,2,0) = Vector_z(*x_axis); +} + +public void get_transform_y_axis( + Transform *transform, + Vector *y_axis ) +{ + fill_Vector( *y_axis, + Transform_elem(*transform,0,1), + Transform_elem(*transform,1,1), + Transform_elem(*transform,2,1) ); +} + +public void set_transform_y_axis( + Transform *transform, + Vector *y_axis ) +{ + Transform_elem(*transform,0,1) = Vector_x(*y_axis); + Transform_elem(*transform,1,1) = Vector_y(*y_axis); + Transform_elem(*transform,2,1) = Vector_z(*y_axis); +} + +public void get_transform_z_axis( + Transform *transform, + Vector *z_axis ) +{ + fill_Vector( *z_axis, + Transform_elem(*transform,0,2), + Transform_elem(*transform,1,2), + Transform_elem(*transform,2,2) ); +} + +public void set_transform_z_axis( + Transform *transform, + Vector *z_axis ) +{ + Transform_elem(*transform,0,2) = Vector_x(*z_axis); + Transform_elem(*transform,1,2) = Vector_y(*z_axis); + Transform_elem(*transform,2,2) = Vector_z(*z_axis); +} + +public void set_transform_x_and_z_axes( + Transform *transform, + Vector *x_axis, + Vector *z_axis ) +{ + Vector n_z, n_y, n_x; + + NORMALIZE_VECTOR( n_z, *z_axis ); + CROSS_VECTORS( n_y, n_z, *x_axis ); + NORMALIZE_VECTOR( n_y, n_y ); + CROSS_VECTORS( n_x, n_z, n_y ); + NORMALIZE_VECTOR( n_x, n_x ); + + set_transform_x_axis( transform, &n_x ); + set_transform_y_axis( transform, &n_y ); + set_transform_z_axis( transform, &n_z ); +} + +public void make_translation_transform( + Vector *trans, + Transform *transform ) +{ + Transform_elem( *transform, 0, 0 ) = 1.0; + Transform_elem( *transform, 0, 1 ) = 0.0; + Transform_elem( *transform, 0, 2 ) = 0.0; + Transform_elem( *transform, 0, 3 ) = Vector_x(*trans); + Transform_elem( *transform, 1, 0 ) = 0.0; + Transform_elem( *transform, 1, 1 ) = 1.0; + Transform_elem( *transform, 1, 2 ) = 0.0; + Transform_elem( *transform, 1, 3 ) = Vector_y(*trans); + Transform_elem( *transform, 2, 0 ) = 0.0; + Transform_elem( *transform, 2, 1 ) = 0.0; + Transform_elem( *transform, 2, 2 ) = 1.0; + Transform_elem( *transform, 2, 3 ) = Vector_z(*trans); + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void make_origin_transform( + Point *origin, + Transform *transform ) +{ + Transform_elem( *transform, 0, 0 ) = 1.0; + Transform_elem( *transform, 0, 1 ) = 0.0; + Transform_elem( *transform, 0, 2 ) = 0.0; + Transform_elem( *transform, 0, 3 ) = Point_x(*origin); + Transform_elem( *transform, 1, 0 ) = 0.0; + Transform_elem( *transform, 1, 1 ) = 1.0; + Transform_elem( *transform, 1, 2 ) = 0.0; + Transform_elem( *transform, 1, 3 ) = Point_y(*origin); + Transform_elem( *transform, 2, 0 ) = 0.0; + Transform_elem( *transform, 2, 1 ) = 0.0; + Transform_elem( *transform, 2, 2 ) = 1.0; + Transform_elem( *transform, 2, 3 ) = Point_z(*origin); + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void make_change_to_bases_transform( + Point *origin, + Vector *x_axis, + Vector *y_axis, + Vector *z_axis, + Transform *transform ) +{ + Transform_elem( *transform, 0, 0 ) = Vector_x( *x_axis ); + Transform_elem( *transform, 0, 1 ) = Vector_x( *y_axis ); + Transform_elem( *transform, 0, 2 ) = Vector_x( *z_axis ); + Transform_elem( *transform, 0, 3 ) = Point_x( *origin ); + Transform_elem( *transform, 1, 0 ) = Vector_y( *x_axis ); + Transform_elem( *transform, 1, 1 ) = Vector_y( *y_axis ); + Transform_elem( *transform, 1, 2 ) = Vector_y( *z_axis ); + Transform_elem( *transform, 1, 3 ) = Point_y( *origin ); + Transform_elem( *transform, 2, 0 ) = Vector_z( *x_axis ); + Transform_elem( *transform, 2, 1 ) = Vector_z( *y_axis ); + Transform_elem( *transform, 2, 2 ) = Vector_z( *z_axis ); + Transform_elem( *transform, 2, 3 ) = Point_z( *origin ); + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void make_change_from_bases_transform( + Point *origin, + Vector *x_axis, + Vector *y_axis, + Vector *z_axis, + Transform *transform ) +{ + Transform_elem( *transform, 0, 0 ) = Vector_x( *x_axis ); + Transform_elem( *transform, 0, 1 ) = Vector_y( *x_axis ); + Transform_elem( *transform, 0, 2 ) = Vector_z( *x_axis ); + Transform_elem( *transform, 0, 3 ) = - DOT_POINT_VECTOR( *origin, *x_axis ); + Transform_elem( *transform, 1, 0 ) = Vector_x( *y_axis ); + Transform_elem( *transform, 1, 1 ) = Vector_y( *y_axis ); + Transform_elem( *transform, 1, 2 ) = Vector_z( *y_axis ); + Transform_elem( *transform, 1, 3 ) = - DOT_POINT_VECTOR( *origin, *y_axis ); + Transform_elem( *transform, 2, 0 ) = Vector_x( *z_axis ); + Transform_elem( *transform, 2, 1 ) = Vector_y( *z_axis ); + Transform_elem( *transform, 2, 2 ) = Vector_z( *z_axis ); + Transform_elem( *transform, 2, 3 ) = - DOT_POINT_VECTOR( *origin, *z_axis ); + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void concat_transforms( + Transform *result, + Transform *t1, + Transform *t2 ) +{ + int i, j, k; + Real sum; + Boolean result_is_also_an_arg; + Transform tmp, *t; + + if( result == t1 || result == t2 ) + { + result_is_also_an_arg = TRUE; + t = &tmp; + } + else + { + result_is_also_an_arg = FALSE; + t = result; + } + + for_less( i, 0, 4 ) + { + for_less( j, 0, 4 ) + { + sum = 0.0; + for_less( k, 0, 4 ) + { + sum += Transform_elem( *t2, i, k ) * + Transform_elem( *t1, k, j ); + } + Transform_elem( *t, i, j ) = sum; + } + } + + if( result_is_also_an_arg ) + { + *result = tmp; + } +} + +public void make_rotation_transform( + Real radians, + int axis, + Transform *transform ) +{ + int a1, a2; + Real c, s; + + a1 = (axis + 1) % N_DIMENSIONS; + a2 = (axis + 2) % N_DIMENSIONS; + + make_identity_transform( transform ); + + c = cos( (double) radians ); + s = sin( (double) radians ); + + Transform_elem( *transform, a1, a1 ) = c; + Transform_elem( *transform, a1, a2 ) = s; + Transform_elem( *transform, a2, a1 ) = -s; + Transform_elem( *transform, a2, a2 ) = c; +} + +public void make_transform_relative_to_point( + Point *point, + Transform *transform, + Transform *rel_transform ) +{ + Vector vect; + Transform to_origin, to_point; + + fill_Vector( vect, Point_x(*point), Point_y(*point), Point_z(*point) ); + + make_translation_transform( &vect, &to_point ); + + SCALE_VECTOR( vect, vect, -1.0 ); + + make_translation_transform( &vect, &to_origin ); + + concat_transforms( rel_transform, &to_origin, transform ); + concat_transforms( rel_transform, rel_transform, &to_point ); +} + +public void make_transform_in_coordinate_system( + Point *origin, + Vector *x_axis, + Vector *y_axis, + Vector *z_axis, + Transform *transform, + Transform *rel_transform ) +{ + Transform to_bases, from_bases; + + make_change_to_bases_transform( origin, x_axis, y_axis, z_axis, &to_bases ); + make_change_from_bases_transform( origin, x_axis, y_axis, z_axis, + &from_bases ); + + concat_transforms( rel_transform, &from_bases, transform ); + concat_transforms( rel_transform, rel_transform, &to_bases ); +} + +public void make_rotation_about_axis( + Vector *axis, + Real angle, + Transform *transform ) +{ + Real c, s, t; + Real txy, txz, tyz, sx, sy, sz; + Real x, y, z; + + c = cos( (double) -angle ); + s = sin( (double) -angle ); + t = 1.0 - c; + + x = Point_x( *axis ); + y = Point_y( *axis ); + z = Point_z( *axis ); + + txy = t * x * y; + txz = t * x * z; + tyz = t * y * z; + + sx = s * x; + sy = s * y; + sz = s * z; + + Transform_elem( *transform, 0, 0 ) = t * x * x + c; + Transform_elem( *transform, 0, 1 ) = txy + sz; + Transform_elem( *transform, 0, 2 ) = txz - sy; + Transform_elem( *transform, 0, 3 ) = 0.0; + + Transform_elem( *transform, 1, 0 ) = txy - sz; + Transform_elem( *transform, 1, 1 ) = t * y * y + c; + Transform_elem( *transform, 1, 2 ) = tyz + sx; + Transform_elem( *transform, 1, 3 ) = 0.0; + + Transform_elem( *transform, 2, 0 ) = txz + sy; + Transform_elem( *transform, 2, 1 ) = tyz - sx; + Transform_elem( *transform, 2, 2 ) = t * z * z + c; + Transform_elem( *transform, 2, 3 ) = 0.0; + + Transform_elem( *transform, 3, 0 ) = 0.0; + Transform_elem( *transform, 3, 1 ) = 0.0; + Transform_elem( *transform, 3, 2 ) = 0.0; + Transform_elem( *transform, 3, 3 ) = 1.0; +} + +public void transform_point( + Transform *transform, + Point *p, + Point *transformed_point ) +{ + Boolean args_same; + Real w; + Point *t, tmp; + + if( p == transformed_point ) + { + args_same = TRUE; + t = &tmp; + } + else + { + args_same = FALSE; + t = transformed_point; + } + + Point_x(*t) = Transform_elem(*transform,0,0) * Point_x(*p) + + Transform_elem(*transform,0,1) * Point_y(*p) + + Transform_elem(*transform,0,2) * Point_z(*p) + + Transform_elem(*transform,0,3); + + Point_y(*t) = Transform_elem(*transform,1,0) * Point_x(*p) + + Transform_elem(*transform,1,1) * Point_y(*p) + + Transform_elem(*transform,1,2) * Point_z(*p) + + Transform_elem(*transform,1,3); + + Point_z(*t) = Transform_elem(*transform,2,0) * Point_x(*p) + + Transform_elem(*transform,2,1) * Point_y(*p) + + Transform_elem(*transform,2,2) * Point_z(*p) + + Transform_elem(*transform,2,3); + + w = Transform_elem(*transform,3,0) * Point_x(*p) + + Transform_elem(*transform,3,1) * Point_y(*p) + + Transform_elem(*transform,3,2) * Point_z(*p) + + Transform_elem(*transform,3,3); + + if( w != 0.0 && w != 1.0 ) + { + Point_x(*t) /= w; + Point_y(*t) /= w; + Point_z(*t) /= w; + } + + if( args_same ) + { + *transformed_point = tmp; + } +} + +public void transform_vector( + Transform *transform, + Vector *v, + Vector *transformed_vector ) +{ + Boolean args_same; + Vector *t, tmp; + + if( v == transformed_vector ) + { + args_same = TRUE; + t = &tmp; + } + else + { + args_same = FALSE; + t = transformed_vector; + } + + Vector_x(*t) = Transform_elem(*transform,0,0) * Vector_x(*v) + + Transform_elem(*transform,0,1) * Vector_y(*v) + + Transform_elem(*transform,0,2) * Vector_z(*v); + + Vector_y(*t) = Transform_elem(*transform,1,0) * Vector_x(*v) + + Transform_elem(*transform,1,1) * Vector_y(*v) + + Transform_elem(*transform,1,2) * Vector_z(*v); + + Vector_z(*t) = Transform_elem(*transform,2,0) * Vector_x(*v) + + Transform_elem(*transform,2,1) * Vector_y(*v) + + Transform_elem(*transform,2,2) * Vector_z(*v); + + if( args_same ) + { + *transformed_vector = tmp; + } +} + +public void inverse_transform_point( + Transform *transform, + Point *p, + Point *transformed_point ) +{ + Real x, y, z; + + x = Point_x(*p) - Transform_elem(*transform,0,3); + y = Point_y(*p) - Transform_elem(*transform,1,3); + z = Point_z(*p) - Transform_elem(*transform,2,3); + + Point_x(*transformed_point) = Transform_elem(*transform,0,0) * x + + Transform_elem(*transform,1,0) * y + + Transform_elem(*transform,2,0) * z; + + Point_y(*transformed_point) = Transform_elem(*transform,0,1) * x + + Transform_elem(*transform,1,1) * y + + Transform_elem(*transform,2,1) * z; + + Point_z(*transformed_point) = Transform_elem(*transform,0,2) * x + + Transform_elem(*transform,1,2) * y + + Transform_elem(*transform,2,2) * z; +} + +public void inverse_transform_vector( + Transform *transform, + Vector *v, + Vector *transformed_vector ) +{ + Real x, y, z; + + x = Vector_x(*v); + y = Vector_y(*v); + z = Vector_z(*v); + + Vector_x(*transformed_vector) = Transform_elem(*transform,0,0) * x + + Transform_elem(*transform,1,0) * y + + Transform_elem(*transform,2,0) * z; + + Vector_y(*transformed_vector) = Transform_elem(*transform,0,1) * x + + Transform_elem(*transform,1,1) * y + + Transform_elem(*transform,2,1) * z; + + Vector_z(*transformed_vector) = Transform_elem(*transform,0,2) * x + + Transform_elem(*transform,1,2) * y + + Transform_elem(*transform,2,2) * z; +} + +public void convert_2d_transform_to_rotation_translation( + Transform *transform, + Real *degrees_clockwise, + Real *x_trans, + Real *y_trans ) +{ + Real x, y; + + x = Transform_elem(*transform, X, X ); + y = Transform_elem(*transform, Y, X ); + + *degrees_clockwise = RAD_TO_DEG * compute_clockwise_rotation( x, y ); + *x_trans = Transform_elem( *transform, X, 3 ); + *y_trans = Transform_elem( *transform, Y, 3 ); +} + +public Real compute_clockwise_rotation( Real x, Real y ) +{ + if( x == 0.0 ) + { + if( y < 0.0 ) + return( PI / 2.0 ); + else if( y > 0.0 ) + return( 3.0 * PI / 2.0 ); + else + return( 0.0 ); + } + else if( y == 0.0 ) + { + if( x > 0.0 ) + return( 0.0 ); + else + return( PI ); + } + else + return( (Real) atan2( (double) y, (double) x ) ); +} + +public void make_identity_transform_2d( Transform_2d *transform ) +{ + Transform_2d_elem( *transform, 0, 0 ) = 1.0; + Transform_2d_elem( *transform, 0, 1 ) = 0.0; + Transform_2d_elem( *transform, 0, 2 ) = 0.0; + Transform_2d_elem( *transform, 1, 0 ) = 0.0; + Transform_2d_elem( *transform, 1, 1 ) = 1.0; + Transform_2d_elem( *transform, 1, 2 ) = 0.0; +} + +public void get_inverse_transform_2d( + Transform_2d *transform, + Transform_2d *inverse ) +{ + Real determinant; + + determinant = Transform_2d_elem(*transform,0,0) * + Transform_2d_elem(*transform,1,1) - + Transform_2d_elem(*transform,1,0) * + Transform_2d_elem(*transform,0,1); + + Transform_2d_elem(*inverse,0,0) = Transform_2d_elem(*transform,1,1) / + determinant; + Transform_2d_elem(*inverse,0,1) = -Transform_2d_elem(*transform,0,1) / + determinant; + Transform_2d_elem(*inverse,0,2) = (Transform_2d_elem(*transform,1,2) * + Transform_2d_elem(*transform,0,1) - + Transform_2d_elem(*transform,1,1) * + Transform_2d_elem(*transform,0,2)) / + determinant; + + Transform_2d_elem(*inverse,1,0) = -Transform_2d_elem(*transform,1,0) / + determinant; + Transform_2d_elem(*inverse,1,1) = Transform_2d_elem(*transform,0,0) / + determinant; + Transform_2d_elem(*inverse,1,2) = (Transform_2d_elem(*transform,1,0) * + Transform_2d_elem(*transform,0,2) - + Transform_2d_elem(*transform,0,0) * + Transform_2d_elem(*transform,1,2)) / + determinant; +} + +public void transform_point_2d( + Transform_2d *transform, + Real x, + Real y, + Real *x_trans, + Real *y_trans ) +{ + *x_trans = Transform_2d_elem(*transform,0,0) * x + + Transform_2d_elem(*transform,0,1) * y + + Transform_2d_elem(*transform,0,2); + *y_trans = Transform_2d_elem(*transform,1,0) * x + + Transform_2d_elem(*transform,1,1) * y + + Transform_2d_elem(*transform,1,2); +}