C$Procedure XDDA ( list voxels intersected by a ray ) SUBROUTINE XDDA ( OBS, RAY, NVOX, NVXLST, NVX, VXLST2 ) C$ Abstract C C Given a ray and a voxel grid, return a list of unique C voxels the ray intersects. C C$ Disclaimer C C THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE C CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S. C GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE C ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE C PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS" C TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY C WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A C PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC C SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE C SOFTWARE AND RELATED MATERIALS, HOWEVER USED. C C IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA C BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT C LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND, C INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS, C REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE C REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY. C C RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF C THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY C CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE C ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE. C C$ Required_Reading C C None. C C$ Keywords C C PLATE VOXEL C C$ Declarations IMPLICIT NONE DOUBLE PRECISION OBS ( 3 ) DOUBLE PRECISION RAY ( 3 ) INTEGER NVOX ( 3 ) INTEGER NVXLST INTEGER NVX INTEGER VXLST2 ( * ) C$ Brief_I/O C C VARIABLE I/O DESCRIPTION C -------- --- -------------------------------------------------- C OBS I Voxel grid coordinates of ray's starting point. C RAY I Direction vector of ray. C NVOX I Dimensions of grid in voxel units. C NVXLST I Order if VXLST2. C NVX O Number of voxels in the VXLST2 list. C VXLST2 O List of voxels intersected by ray. C C$ Detailed_Input C C OBS Voxel grid coordinates of ray's starting point. C C RAY Direction vector of ray from OBS, note ray direction C is invarient with regards to the model to voxel C space mapping. C C NVOX The integer 3-vector containing the dimensions of C the voxel grid in voxel units, [Nx, Ny, Nz]. C C NVXLST The order of VXLST2 as defined in the calling routine. C The counter NVX will not exceed this value. C C$ Detailed_Output C C NVX Number of voxel IDs contained by the VXLST2 list. C C VXLST2 List of voxel IDs intersected by ray. Voxel ID C ordering corresponds to order in which ray C intersected the voxels. C C VXLST2(1) - first intersected voxel C VXLST2(2) - second intersected voxel C C ...etc... C C$ Parameters C C None. C C$ Files C C None. C C$ Exceptions C C 1) The SPICE error SPICE(RAYISZEROVECTOR) signals when C the input RAY has all zero components. C C 2) The SPICE error SPICE(NOTINGRID) signals when C the input OBS describes a position not in or on the C voxel grid. C C 3) The SPICE error SPICE(ARRAYTOOSMALL) signals when C the value of the NVX counter (number of intersected C voxels) exceeds the order of the VXLST2 input vector. C C$ Particulars C C Later. C C$ Examples C C C C C Give a position within the voxel grid, COORD, and a C C 3-vector, BSIGHT, centered on COORD, return the list of C C voxel IDS the BSIGHT vector intersets. C C CALL XDDA ( COORD, BSIGHT, NVOX, NVXLST, NVXOUT, VOXLST ) C C$ Restrictions C C None. C C$ Author_and_Institution C C J.A. Bytof (JPL) C N.J. Bachman (JPL) C C$ Literature_References C C Woo, A. "Fast Ray-Box Intersection", Graphic Gems I, 395-396. C C$ Version C C- SPICELIB Version 1.2.0, 15-JUL-2006 (EDW) C C Algorithm now implements Nat Bachman's C modification - the driver direction corresponds C to ray component with largest magnitude. C C RMDUPI call removed as it is no longer needed. C C- SPICELIB Version 1.1.0, 19-OCT-2004 (EDW) C C Added logic to remove duplicate voxel IDs from C the return list. Extended programing comments. C C- SPICELIB Version 1.0.1, 26-AUG-2002 (BVS) C C Replaced WRITE with normal error reporting calls. C C- SPICELIB Version 1.0.0, 03-FEB-1999 (JAB) C C-& C$ Index_Entries C C list voxels intersected by a ray C C-& C C SPICELIB functions C DOUBLE PRECISION DPMAX LOGICAL RETURN LOGICAL VZERO C C Local variables C DOUBLE PRECISION AX2ERR DOUBLE PRECISION AX3ERR DOUBLE PRECISION MAXCMP DOUBLE PRECISION OBSLOC ( 3 ) DOUBLE PRECISION S12 DOUBLE PRECISION S13 INTEGER I INTEGER IAXIS ( 3 ) INTEGER ICOORD ( 3 ) INTEGER NEXT ( 3 ) INTEGER STEP ( 3 ) INTEGER VOX2ID INTEGER VXLST ( 3 ) LOGICAL INGRD C C Saved variables C SAVE NEXT C C Initial values C DATA NEXT / 2, 3, 1 / C C Standard SPICELIB error handling. C IF ( RETURN () ) THEN RETURN END IF CALL CHKIN ( 'XDDA' ) C C Check if ray direction vector is a zero vector. C IF ( VZERO( RAY ) ) THEN CALL SETMSG ( 'Ray is the zero vector.' ) CALL SIGERR ( 'SPICE(RAYISZEROVECTOR)' ) CALL CHKOUT ( 'XDDA' ) RETURN END IF C C Find the largest component of the direction vector. C IAXIS(1) = 1 MAXCMP = ABS( RAY(1) ) DO I = 2, 3 IF ( ABS( RAY(I) ) .GT. MAXCMP ) THEN IAXIS(1) = I MAXCMP = ABS( RAY(I) ) END IF END DO C C Set the indices of the orthogonal components of the direction C vector. We maintain a right-handed relationship between the axes C labeled by IAXIS(1), IAXIS(2), and IAXIS(3): the third axis is C the cross product of the first and second. C IAXIS(2) = NEXT ( IAXIS(1) ) IAXIS(3) = NEXT ( IAXIS(2) ) C C Which voxel contains the observer? Truncate the observer C coordinates to integers. Add 1 to each coord to compensate C for 1 based counting. C C WRITE (*,*) 'IAXIS = ', IAXIS DO I = 1, 3 ICOORD(I) = 1 + INT ( OBS(IAXIS(I)) ) VXLST ( IAXIS(I) ) = ICOORD(I) END DO C C Check the observer coordinates exist as part C of the voxel grid. C IF ( .NOT. INGRD(NVOX, VXLST ) ) THEN CALL SETMSG ( 'Ray origin not in voxel grid.' ) CALL SIGERR ( 'SPICE(NOTINGRID)' ) CALL CHKOUT ( 'XDDA' ) RETURN END IF VXLST2(1) = VOX2ID ( VXLST, NVOX ) C C Initialize the counter for number of voxels the ray C intercepts. The test above confirmed the observer C position lies in a voxel belonging to the grid. C NVX = 1 C C Calculate the relative location of starting point C within the voxel. The coordinates of a voxel's corners C are integer values with each voxel side length 1 C (in voxel coords). C C The DMOD calls equate to the function: C C OBS(I) - DBLE( INT( OBS(IAXIS(I)) ) ) C C So the X0, Y0, Z0 variables have the fractional C component values of the corresponding OBS(I). C DO I = 1, 3 OBSLOC(I) = DMOD( OBS(IAXIS(I)), 1.D0 ) END DO C C Determine the error term initial values and increments. C C Different quadrant directions result in slightly C different initial conditions; the only effect of C quadrant in the main loop is in the step direction. C AX2ERR = DPMAX() AX3ERR = AX2ERR IF ( RAY( IAXIS(2) ) .NE. 0.D0 ) THEN S12 = ABS ( RAY( IAXIS(1) )/ RAY( IAXIS(2) ) ) IF ( RAY(IAXIS(1)) .GT. 0.D0 ) THEN IF( RAY(IAXIS(2)) .GT. 0.D0 ) THEN AX2ERR = S12*(1.D0 - OBSLOC(2)) + OBSLOC(1) - 1.D0 ELSE AX2ERR = S12*OBSLOC(2) + OBSLOC(1) - 1.D0 END IF ELSE IF( RAY(IAXIS(2)) .GT. 0.D0 ) THEN AX2ERR = S12*(1.D0 - OBSLOC(2)) - OBSLOC(1) ELSE AX2ERR = S12*OBSLOC(2) - OBSLOC(1) END IF END IF ELSE S12 = 0.D0 END IF IF ( RAY( IAXIS(3) ) .NE. 0.D0 ) THEN S13 = ABS( RAY( IAXIS(1) ) / RAY( IAXIS(3) ) ) IF ( RAY( IAXIS(1) ) .GT. 0.D0 ) THEN IF( RAY( IAXIS(3) ) .GT. 0.D0 ) THEN AX3ERR = S13*(1.D0 - OBSLOC(3) ) + OBSLOC(1) - 1.D0 ELSE AX3ERR = S13*OBSLOC(3) + OBSLOC(1) - 1.D0 END IF ELSE IF( RAY( IAXIS(3) ) .GT. 0.D0 ) THEN AX3ERR = S13*(1.D0 - OBSLOC(3)) - OBSLOC(1) ELSE AX3ERR = S13*OBSLOC(3) - OBSLOC(1) END IF END IF ELSE S13 = 0.D0 END IF C WRITE (*,*) 'S12, S13 = ', S12, S13 C C A voxel, in voxel space, has length one for each size, C a step has magnitude one or zero, zero only when the C magnitude of a component of RAY has value zero. C DO I = 1, 3 IF ( RAY( IAXIS(I) ) .GT. 0.D0 ) THEN C C Positive component direction, positive step. C STEP(I) = 1 ELSE IF ( RAY( IAXIS(I) ) .LT. 0.D0 ) THEN C C Negative component direction, negative step. C STEP(I) = -1 ELSE IF ( RAY( IAXIS(I) ) .EQ. 0.D0 ) THEN C C No component in this direction, no step. C STEP(I) = 0 END IF END DO C C Follow the ray until it exits the voxel grid. C DO WHILE ( INGRD( NVOX, VXLST ) ) IF ( (AX2ERR .LT. 0.D0) .OR. (AX3ERR .LT. 0.D0) ) THEN C C Ray has crossed over into the next voxel in IAXIS(2) or C IAXIS(3) C IF ( AX2ERR .LT. AX3ERR ) THEN ICOORD(2) = ICOORD(2) + STEP(2) AX2ERR = AX2ERR + S12 NVX = NVX + 1 ELSE ICOORD(3) = ICOORD(3) + STEP(3) AX3ERR = AX3ERR + S13 NVX = NVX + 1 END IF ELSE C C No change in IAXIS(2) or IAXIS(3), step in IAXIS(1). C ICOORD(1) = ICOORD(1) + STEP(1) NVX = NVX + 1 AX2ERR = AX2ERR - 1.D0 AX3ERR = AX3ERR - 1.D0 C WRITE (*,*) 'AX2ERR, AX3ERR = ', AX2ERR, AX3ERR END IF C C Check we have room in VXLST2. C IF( NVX .GT. NVXLST ) THEN CALL SETMSG ( 'Index larger than array. '// . 'Index = #1. Array size = #2.' ) CALL ERRINT ( '#1', NVX ) CALL ERRINT ( '#2', NVXLST ) CALL SIGERR ( 'SPICE(ARRAYTOOSMALL)' ) CALL CHKOUT ( 'XDDA' ) RETURN END IF C C Pack the voxel indices into VXLST using C the values calculated in this loop pass. C DO I = 1, 3 VXLST( IAXIS(I) ) = ICOORD(I) END DO C C Store the voxel ID corresponding to the C VXLST voxel indices. C VXLST2(NVX) = VOX2ID( VXLST, NVOX ) END DO C C Subtract one off the voxel count since the final voxel C exists outside the grid. C NVX = NVX - 1 C C Standard SPICE error handling. C CALL CHKOUT ( 'XDDA' ) RETURN END