| /*************************************************************************/ |
| /* */ |
| /* Copyright (c) 1994 Stanford University */ |
| /* */ |
| /* All rights reserved. */ |
| /* */ |
| /* Permission is given to use, copy, and modify this software for any */ |
| /* non-commercial purpose as long as this copyright notice is not */ |
| /* removed. All other uses, including redistribution in whole or in */ |
| /* part, are forbidden without prior written permission. */ |
| /* */ |
| /* This software is provided with absolutely no warranty and no */ |
| /* support. */ |
| /* */ |
| /*************************************************************************/ |
| |
| /****************************************************************************** |
| * * |
| * view.c: Compute viewing direction. * |
| * * |
| ******************************************************************************/ |
| |
| #include "incl.h" |
| |
| #define XAXIS 1 /* Positive X-axis points rightwards */ |
| #define YAXIS 2 /* Positive Y-axis points upwards */ |
| #define ZAXIS 3 /* Positive Z-axis points into screen */ |
| |
| float transformation_matrix[4][4];/* current transformation matrix of floats */ |
| float out_invvertex[2][2][2][NM]; /* Image and object space centers */ |
| /* of outer voxels in output map */ |
| float uout_invvertex[2][2][2][NM];/* Image and object space vertices */ |
| /* of output map unit voxel */ |
| |
| long frust_len; /* Size of clipping frustum */ |
| /* (mins will be 0 in this program, */ |
| /* {x,y}len will be <= IK{X,Y}SIZE) */ |
| float out_diag_len[NM]; /* Per-axis lengths and combined length of */ |
| float out_diag_length; /* diagonal of input map in image space */ |
| float depth_cueing[MAX_OUTLEN]; /* Pre-computed table of depth cueing */ |
| long image_zlen; /* number of samples along viewing ray */ |
| float in_max[NM]; /* Pre-computed clipping aids */ |
| long opc_xlen,opc_xylen; |
| long norm_xlen,norm_xylen; |
| |
| float invmatrix[4][4]; /* Inverse of viewing matrix: */ |
| /* 4 x 4 viewing transformation matrix */ |
| /* (orthographic projection used, so */ |
| /* matrix[][3] = 0, matrix[3][3] = 1) */ |
| |
| EXTERN_ENV |
| |
| void Compute_Pre_View() |
| { |
| long i, outz; |
| |
| for (i=0; i<NM; i++) |
| out_diag_len[i] = opc_len[i]-1; |
| out_diag_length = out_diag_len[X]*out_diag_len[X] + |
| out_diag_len[Y]*out_diag_len[Y] + |
| out_diag_len[Z]*out_diag_len[Z]; |
| out_diag_length = sqrt(out_diag_length); |
| frust_len = NINT(out_diag_length)+1; |
| image_zlen = frust_len - 1; |
| |
| /* Pre-compute table of depth cueing attenuation fractions as */ |
| /* exponential falloff from intensity of depth_hither*full at */ |
| /* hither (outz=0) to depth_yon*full at yon (frust_len-1). */ |
| /* Exponent > 1.0 falls off slower than linear, < 1.0 falls faster. */ |
| /* Clip fractions to valid range to allow agressive user */ |
| /* to set hither or yon outside the range 0.0..1.0. */ |
| if (image_zlen > MAX_OUTLEN) |
| Error ("MAX_OUTLEN exceeded in Ray_Trace.\n"); |
| for (outz=0; outz<=image_zlen; outz++) { |
| depth_cueing[outz] = depth_hither - |
| pow((float)(outz)/(float)image_zlen,depth_exponent) * |
| (depth_hither - depth_yon); |
| depth_cueing[outz] = MIN(MAX(depth_cueing[outz],0.0),1.0); |
| } |
| |
| /* Pre-compute clipping aids */ |
| in_max[X] = (float)(opc_len[X]-1)-SMALL-1.0/(float)MAX_PIXEL; |
| in_max[Y] = (float)(opc_len[Y]-1)-SMALL-1.0/(float)MAX_PIXEL; |
| in_max[Z] = (float)(opc_len[Z]-1)-SMALL-1.0/(float)MAX_PIXEL; |
| |
| /* Pre-compute subscripting aids */ |
| opc_xlen = opc_len[X] + 1; |
| opc_xylen = opc_len[X] * opc_len[Y] + 1; |
| |
| norm_xlen = norm_len[X] + 1; |
| norm_xylen = norm_len[X] * norm_len[Y] + 1; |
| } |
| |
| |
| void Select_View(delta_angle, axis) |
| float delta_angle; |
| long axis; |
| |
| { |
| Load_Identity_Matrix(invmatrix); |
| |
| /* Moves input map from all-positive octant to center of */ |
| /* coordinate system so subsequent rotations spin object */ |
| /* in place. Must be first after initialization to work. */ |
| Inverse_Concatenate_Translation(invmatrix, |
| (float)(out_diag_len[X])/2.0, |
| (float)(out_diag_len[Y])/2.0, |
| (float)(out_diag_len[Z])/2.0); |
| |
| /* scale dataset size in Z-direction for 256x256x113 dataset */ |
| Inverse_Concatenate_Scaling(invmatrix,1.0,1.0,1.0/(float)ZSCALE); |
| |
| /* rotation about axis by angle */ |
| if (frame != 0) |
| angle[axis] = angle[axis] + delta_angle; |
| Inverse_Concatenate_Rotation(invmatrix,XAXIS,-angle[X]); |
| Inverse_Concatenate_Rotation(invmatrix,YAXIS,-angle[Y]); |
| /* Inverse_Concatenate_Rotation(invmatrix,ZAXIS,-angle[Z]);*/ |
| Inverse_Concatenate_Rotation(invmatrix,XAXIS,30.0); |
| |
| /* Moves input map from center of coordinate system back */ |
| /* to within all-positive octant such that any rotation */ |
| /* (e.g. by using pre-matrix, rotations, and anpost-matrix) */ |
| /* exactly falls in the octant, preventing low-side clipping.*/ |
| /* Fails if non-isotropic image space scaling is applied */ |
| /* (i.e. different scaling in X,Y, or Z after rotations). */ |
| Inverse_Concatenate_Translation(invmatrix, |
| -out_diag_length/2.0, |
| -out_diag_length/2.0, |
| -out_diag_length/2.0); |
| |
| /* Insures that frustum entirely encloses any rotation of */ |
| /* map assuming they fall within the all-positive octant */ |
| /* (e.g. by using pre-matrix, rotations, and anpost-matrix), */ |
| Load_Transformation_Matrix(invmatrix); |
| Compute_Input_Dimensions(); |
| Compute_Input_Unit_Vector(); |
| } |
| |
| |
| void Compute_Input_Dimensions() |
| { |
| long x,y,z; |
| float in_invvertex[2][2][2][NM]; /* Image and object space centers */ |
| |
| in_invvertex[0][0][0][X] = 0; |
| in_invvertex[0][0][0][Y] = 0; |
| in_invvertex[0][0][0][Z] = 0; |
| |
| in_invvertex[0][0][1][X] = frust_len-1; |
| in_invvertex[0][0][1][Y] = 0; |
| in_invvertex[0][0][1][Z] = 0; |
| |
| in_invvertex[0][1][0][X] = 0; |
| in_invvertex[0][1][0][Y] = frust_len-1; |
| in_invvertex[0][1][0][Z] = 0; |
| |
| in_invvertex[0][1][1][X] = frust_len-1; |
| in_invvertex[0][1][1][Y] = frust_len-1; |
| in_invvertex[0][1][1][Z] = 0; |
| |
| in_invvertex[1][0][0][X] = 0; |
| in_invvertex[1][0][0][Y] = 0; |
| in_invvertex[1][0][0][Z] = frust_len-1; |
| |
| in_invvertex[1][0][1][X] = frust_len-1; |
| in_invvertex[1][0][1][Y] = 0; |
| in_invvertex[1][0][1][Z] = frust_len-1; |
| |
| in_invvertex[1][1][0][X] = 0; |
| in_invvertex[1][1][0][Y] = frust_len-1; |
| in_invvertex[1][1][0][Z] = frust_len-1; |
| |
| in_invvertex[1][1][1][X] = frust_len-1; |
| in_invvertex[1][1][1][Y] = frust_len-1; |
| in_invvertex[1][1][1][Z] = frust_len-1; |
| |
| for (z=0; z<2; z++) { |
| for (y=0; y<2; y++) { |
| for (x=0; x<2; x++) { |
| Transform_Point(in_invvertex[z][y][x][X], |
| in_invvertex[z][y][x][Y], |
| in_invvertex[z][y][x][Z], |
| &out_invvertex[z][y][x][X], |
| &out_invvertex[z][y][x][Y], |
| &out_invvertex[z][y][x][Z]); |
| } |
| } |
| } |
| } |
| |
| |
| void Compute_Input_Unit_Vector() |
| { |
| long x,y,z; |
| float uin_invvertex[2][2][2][NM]; |
| |
| uin_invvertex[0][0][0][X] = 0; |
| uin_invvertex[0][0][0][Y] = 0; |
| uin_invvertex[0][0][0][Z] = 0; |
| |
| uin_invvertex[0][0][1][X] = 1.0; |
| uin_invvertex[0][0][1][Y] = 0; |
| uin_invvertex[0][0][1][Z] = 0; |
| |
| uin_invvertex[0][1][0][X] = 0; |
| uin_invvertex[0][1][0][Y] = 1.0; |
| uin_invvertex[0][1][0][Z] = 0; |
| |
| uin_invvertex[0][1][1][X] = 1.0; |
| uin_invvertex[0][1][1][Y] = 1.0; |
| uin_invvertex[0][1][1][Z] = 0; |
| |
| uin_invvertex[1][0][0][X] = 0; |
| uin_invvertex[1][0][0][Y] = 0; |
| uin_invvertex[1][0][0][Z] = 1.0; |
| |
| uin_invvertex[1][0][1][X] = 1.0; |
| uin_invvertex[1][0][1][Y] = 0; |
| uin_invvertex[1][0][1][Z] = 1.0; |
| |
| uin_invvertex[1][1][0][X] = 0; |
| uin_invvertex[1][1][0][Y] = 1.0; |
| uin_invvertex[1][1][0][Z] = 1.0; |
| |
| uin_invvertex[1][1][1][X] = 1.0; |
| uin_invvertex[1][1][1][Y] = 1.0; |
| uin_invvertex[1][1][1][Z] = 1.0; |
| |
| for (z=0; z<2; z++) { |
| for (y=0; y<2; y++) { |
| for (x=0; x<2; x++) { |
| Transform_Point(uin_invvertex[z][y][x][X], |
| uin_invvertex[z][y][x][Y], |
| uin_invvertex[z][y][x][Z], |
| &uout_invvertex[z][y][x][X], |
| &uout_invvertex[z][y][x][Y], |
| &uout_invvertex[z][y][x][Z]); |
| } |
| } |
| } |
| } |
| |
| |
| void Load_Transformation_Matrix(matrix) |
| float matrix[4][4]; |
| { |
| Copy_Matrix(matrix,transformation_matrix); |
| } |
| |
| |
| void Transform_Point(xold,yold,zold,xnew,ynew,znew) |
| float xold,yold,zold; |
| float *xnew,*ynew,*znew; |
| { |
| *xnew = |
| xold * transformation_matrix[0][0] + |
| yold * transformation_matrix[1][0] + |
| zold * transformation_matrix[2][0] + |
| transformation_matrix[3][0]; |
| |
| *ynew = |
| xold * transformation_matrix[0][1] + |
| yold * transformation_matrix[1][1] + |
| zold * transformation_matrix[2][1] + |
| transformation_matrix[3][1]; |
| |
| *znew = |
| xold * transformation_matrix[0][2] + |
| yold * transformation_matrix[1][2] + |
| zold * transformation_matrix[2][2] + |
| transformation_matrix[3][2]; |
| } |
| |
| |
| void Inverse_Concatenate_Translation(matrix,xoffset,yoffset,zoffset) |
| float matrix[4][4],xoffset,yoffset,zoffset; |
| { |
| float translation_matrix[4][4]; |
| Load_Translation_Matrix(translation_matrix,xoffset,yoffset,zoffset); |
| Inverse_Concatenate_Transform(matrix,translation_matrix); |
| } |
| |
| |
| void Inverse_Concatenate_Scaling(matrix,xscale,yscale,zscale) |
| float matrix[4][4],xscale,yscale,zscale; |
| { |
| float scaling_matrix[4][4]; |
| Load_Scaling_Matrix(scaling_matrix,xscale,yscale,zscale); |
| Inverse_Concatenate_Transform(matrix,scaling_matrix); |
| } |
| |
| |
| void Inverse_Concatenate_Rotation(matrix,axis,angle) |
| float matrix[4][4],angle; |
| long axis; |
| { |
| float rotation_matrix[4][4]; |
| Load_Rotation_Matrix(rotation_matrix,axis,angle); |
| Inverse_Concatenate_Transform(matrix,rotation_matrix); |
| } |
| |
| |
| void Load_Identity_Matrix(matrix) |
| float matrix[4][4]; |
| { |
| long i,j; |
| for (i=0; i<4; i++) { |
| for (j=0; j<4; j++) { |
| matrix[i][j] = 0; |
| } |
| matrix[i][i] = 1.0; |
| } |
| } |
| |
| void Load_Translation_Matrix(matrix,xoffset,yoffset,zoffset) |
| float matrix[4][4],xoffset,yoffset,zoffset; |
| { |
| Load_Identity_Matrix(matrix); |
| matrix[3][0] = xoffset; |
| matrix[3][1] = yoffset; |
| matrix[3][2] = zoffset; |
| } |
| |
| |
| void Load_Scaling_Matrix(matrix,xscale,yscale,zscale) |
| float matrix[4][4],xscale,yscale,zscale; |
| { |
| Load_Identity_Matrix(matrix); |
| matrix[0][0] = xscale; |
| matrix[1][1] = yscale; |
| matrix[2][2] = zscale; |
| } |
| |
| |
| void Load_Rotation_Matrix(matrix,axis,angle) |
| float matrix[4][4],angle; |
| long axis; |
| { |
| Load_Identity_Matrix(matrix); |
| if (axis == XAXIS) { |
| matrix[1][1] = cos(angle/180.0*PI); |
| matrix[1][2] = sin(angle/180.0*PI); |
| matrix[2][1] = -sin(angle/180.0*PI); |
| matrix[2][2] = cos(angle/180.0*PI); |
| } |
| else if (axis == YAXIS) { |
| matrix[0][0] = cos(angle/180.0*PI); |
| matrix[0][2] = -sin(angle/180.0*PI); |
| matrix[2][0] = sin(angle/180.0*PI); |
| matrix[2][2] = cos(angle/180.0*PI); |
| } |
| else if (axis == ZAXIS) { |
| matrix[0][0] = cos(angle/180.0*PI); |
| matrix[0][1] = sin(angle/180.0*PI); |
| matrix[1][0] = -sin(angle/180.0*PI); |
| matrix[1][1] = cos(angle/180.0*PI); |
| } |
| } |
| |
| |
| void Concatenate_Transform(composite_matrix,transformation_matrix) |
| float composite_matrix[][4],transformation_matrix[][4]; |
| { |
| float temp_matrix[4][4]; |
| Multiply_Matrices(composite_matrix,transformation_matrix,temp_matrix); |
| Copy_Matrix(temp_matrix,composite_matrix); |
| } |
| |
| void Inverse_Concatenate_Transform(composite_matrix,transformation_matrix) |
| float composite_matrix[][4],transformation_matrix[][4]; |
| { |
| float temp_matrix[4][4]; |
| Multiply_Matrices(transformation_matrix,composite_matrix,temp_matrix); |
| Copy_Matrix(temp_matrix,composite_matrix); |
| } |
| |
| |
| void Multiply_Matrices(input_matrix1,input_matrix2,output_matrix) |
| float input_matrix1[][4],input_matrix2[][4],output_matrix[][4]; |
| { |
| long i,j; |
| for (i=0; i<4; i++) { |
| for (j=0; j<4; j++) { |
| output_matrix[i][j] = |
| input_matrix1[i][0] * |
| input_matrix2[0][j] + |
| input_matrix1[i][1] * |
| input_matrix2[1][j] + |
| input_matrix1[i][2] * |
| input_matrix2[2][j] + |
| input_matrix1[i][3] * |
| input_matrix2[3][j]; |
| } |
| } |
| } |
| |
| |
| void Copy_Matrix(input_matrix,output_matrix) |
| float input_matrix[][4],output_matrix[][4]; |
| { |
| long i,j; |
| for (i=0; i<4; i++) { |
| for (j=0; j<4; j++) { |
| output_matrix[i][j] = input_matrix[i][j]; |
| } |
| } |
| } |
