shoot.c

Go to the documentation of this file.
00001 /*                         S H O O T . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 2000-2006 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * as published by the Free Software Foundation; either version 2 of
00010  * the License, or (at your option) any later version.
00011  *
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Library General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public
00018  * License along with this file; see the file named COPYING for more
00019  * information.
00020  */
00021 
00022 /** @addtogroup ray */
00023 /*@{*/
00024 
00025 /** @file shoot.c
00026  *      Ray Tracing program shot coordinator.
00027  *
00028  *  This is the heart of LIBRT's ray-tracing capability.
00029  *
00030  *  Given a ray, shoot it at all the relevant parts of the model,
00031  *  (building the finished_segs chain), and then call rt_boolregions()
00032  *  to build and evaluate the partition chain.
00033  *  If the ray actually hit anything, call the application's
00034  *  a_hit() routine with a pointer to the partition chain,
00035  *  otherwise, call the application's a_miss() routine.
00036  *
00037  *  It is important to note that rays extend infinitely only in the
00038  *  positive direction.  The ray is composed of all points P, where
00039  *
00040  *      P = r_pt + K * r_dir
00041  *
00042  *  for K ranging from 0 to +infinity.  There is no looking backwards.
00043  *
00044  *  Authors -
00045  *      Michael John Muuss
00046  *      Glenn Durfee
00047  *
00048  *  Source -
00049  *      The U. S. Army Research Laboratory
00050  *      Aberdeen Proving Ground, Maryland  21005-5068  USA
00051  */
00052 
00053 #ifndef lint
00054 static const char RCSshoot[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/shoot.c,v 14.16 2006/09/16 02:04:26 lbutler Exp $ (BRL)";
00055 #endif
00056 
00057 #include "common.h"
00058 
00059 #include <stdio.h>
00060 #include <math.h>
00061 #ifdef HAVE_STRING_H
00062 #  include <string.h>
00063 #else
00064 #  include <strings.h>
00065 #endif
00066 #include "machine.h"
00067 #include "vmath.h"
00068 #include "bu.h"
00069 #include "raytrace.h"
00070 #include "plot3.h"
00071 #include "./debug.h"
00072 
00073 struct resource rt_uniresource;         /* Resources for uniprocessor */
00074 
00075 extern void     rt_plot_cell(const union cutter *cutp, const struct rt_shootray_status *ssp, struct bu_list *waiting_segs_hd, struct rt_i *rtip);               /* at end of file */
00076 
00077 #define V3PT_DEPARTING_RPP(_step, _lo, _hi, _pt ) \
00078                 PT_DEPARTING_RPP(_step, _lo, _hi, (_pt)[X], (_pt)[Y], (_pt)[Z] )
00079 #define PT_DEPARTING_RPP(_step, _lo, _hi, _px, _py, _pz ) \
00080                 (   ((_step)[X] <= 0 && (_px) < (_lo)[X]) || \
00081                     ((_step)[X] >= 0 && (_px) > (_hi)[X]) || \
00082                     ((_step)[Y] <= 0 && (_py) < (_lo)[Y]) || \
00083                     ((_step)[Y] >= 0 && (_py) > (_hi)[Y]) || \
00084                     ((_step)[Z] <= 0 && (_pz) < (_lo)[Z]) || \
00085                     ((_step)[Z] >= 0 && (_pz) > (_hi)[Z])   )
00086 
00087 /*
00088  *                      R T _ R E S _ P I E C E S _ I N I T
00089  *
00090  *  Allocate the per-processor state variables needed
00091  *  to support rt_shootray()'s use of 'solid pieces'.
00092  */
00093 void
00094 rt_res_pieces_init(struct resource *resp, struct rt_i *rtip)
00095 {
00096         struct rt_piecestate    *psptab;
00097         struct rt_piecestate    *psp;
00098         struct soltab           *stp;
00099 
00100         RT_CK_RESOURCE(resp);
00101         RT_CK_RTI(rtip);
00102 
00103         psptab = bu_calloc( sizeof(struct rt_piecestate),
00104                 rtip->rti_nsolids_with_pieces, "re_pieces[]" );
00105         resp->re_pieces = psptab;
00106 
00107         RT_VISIT_ALL_SOLTABS_START( stp, rtip )  {
00108                 RT_CK_SOLTAB(stp);
00109                 if( stp->st_npieces <= 1 )  continue;
00110                 psp = &psptab[stp->st_piecestate_num];
00111                 psp->magic = RT_PIECESTATE_MAGIC;
00112                 psp->stp = stp;
00113                 psp->shot = bu_bitv_new(stp->st_npieces);
00114                 rt_htbl_init( &psp->htab, 8, "psp->htab" );
00115                 psp->cutp = CUTTER_NULL;
00116         } RT_VISIT_ALL_SOLTABS_END
00117 }
00118 
00119 /*
00120  *                      R T _ R E S _ P I E C E S _ C L E A N
00121  */
00122 void
00123 rt_res_pieces_clean(struct resource *resp, struct rt_i *rtip)
00124 {
00125         struct rt_piecestate    *psp;
00126         int                     i;
00127 
00128         RT_CK_RESOURCE(resp);
00129         RT_CK_RTI(rtip);
00130 
00131         if( !resp->re_pieces ) {
00132                 rtip->rti_nsolids_with_pieces = 0;
00133                 return;
00134         }
00135         for( i = rtip->rti_nsolids_with_pieces-1; i >= 0; i-- )  {
00136                 psp = &resp->re_pieces[i];
00137                 RT_CK_PIECESTATE(psp);
00138                 rt_htbl_free(&psp->htab);
00139                 bu_bitv_free(psp->shot);
00140                 psp->shot = NULL;       /* sanity */
00141                 psp->magic = 0;
00142         }
00143         bu_free( (char *)resp->re_pieces, "re_pieces[]" );
00144         resp->re_pieces = NULL;
00145 
00146         bu_ptbl_free( &resp->re_pieces_pending );
00147 
00148         rtip->rti_nsolids_with_pieces = 0;
00149 }
00150 
00151 /*
00152  *                      R T _ F I N D _ N U G R I D
00153  *
00154  *  Along the given axis, find which NUgrid cell this value lies in.
00155  *  Use method of binary subdivision.
00156  */
00157 int
00158 rt_find_nugrid(const struct nugridnode *nugnp, int axis, fastf_t val)
00159 {
00160         int     min;
00161         int     max;
00162         int     lim = nugnp->nu_cells_per_axis[axis]-1;
00163         int     cur;
00164 
00165         if( val < nugnp->nu_axis[axis][0].nu_spos ||
00166             val > nugnp->nu_axis[axis][lim].nu_epos )
00167                 return -1;
00168 
00169         min = 0;
00170         max = lim;
00171 again:
00172         cur = (min + max) / 2;
00173 
00174         if( cur <= 0 )  return 0;
00175         if( cur >= lim )  return lim;
00176 
00177         if( val < nugnp->nu_axis[axis][cur].nu_spos )  {
00178                 max = cur;
00179                 goto again;
00180         }
00181         if( val > nugnp->nu_axis[axis][cur+1].nu_epos )  {
00182                 min = cur+1;
00183                 goto again;
00184         }
00185         if( val < nugnp->nu_axis[axis][cur].nu_epos )
00186                 return cur;
00187         else
00188                 return cur+1;
00189 }
00190 
00191 
00192 /*
00193  *                      R T _ A D V A N C E _ T O _ N E X T _ C E L L
00194  */
00195 
00196 const union cutter *
00197 rt_advance_to_next_cell(register struct rt_shootray_status *ssp)
00198 {
00199         register const union cutter             *cutp, *curcut = ssp->curcut;
00200         register const struct application       *ap = ssp->ap;
00201         register fastf_t                        t0, px, py, pz;
00202         int                                     push_flag = 0;
00203         double                                  fraction;
00204         int                                     exponent;
00205 
00206         ssp->box_num++;
00207 
00208         if( curcut == &ssp->ap->a_rt_i->rti_inf_box )  {
00209                 /* Last pass did the infinite solids, there is nothing more */
00210                 ssp->curcut = CUTTER_NULL;
00211                 return CUTTER_NULL;
00212         }
00213 
00214 #if EXTRA_SAFETY
00215         if( curcut == CUTTER_NULL ) {
00216                 bu_log(
00217                    "rt_advance_to_next_cell: warning: ssp->curcut not set\n" );
00218                 ssp->curcut = curcut = &ap->a_rt_i->rti_CutHead;
00219         }
00220 #endif
00221 
00222         for( ;; ) {
00223                 /* Set cutp to CUTTER_NULL.  If it fails to become set in the
00224                    following switch statement, we know that we have exited the
00225                    subnode.  If this subnode is the highest-level node, then
00226                    we are done advancing the ray through the model. */
00227                 cutp = CUTTER_NULL;
00228                 push_flag = 0;
00229 
00230                 /*
00231                  *  The point corresponding to the box_start distance
00232                  *  may not be in the "right" place,
00233                  *  due to the effects of floating point fuzz:
00234                  *  1)  The point might lie just outside
00235                  *      the model RPP, resulting in the point not
00236                  *      falling within the RPP of the indicated cell,
00237                  *      or
00238                  *  2)  The poing might lie just a little bit on the
00239                  *      wrong side of the cell wall, resulting in
00240                  *      the ray getting "stuck", and needing rescuing
00241                  *      all the time by the error recovery code below.
00242                  *  Therefore, "nudge" the point just slightly into the
00243                  *  next cell by adding OFFSET_DIST.
00244                  *  XXX At present, a cell is never less than 1mm wide.
00245                  *  XXX The value of OFFSET_DIST should be some
00246                  *      percentage of the cell's smallest dimension,
00247                  *      rather than an absolute distance in mm.
00248                  *      This will prevent doing microscopic models.
00249                  */
00250                 t0 = ssp->box_start;
00251                 /* NB: can't compute px,py,pz here since t0 may advance
00252                    in the following statement! */
00253 
00254 top:            switch( curcut->cut_type ) {
00255                 case CUT_NUGRIDNODE: {
00256                         /*
00257                          ****************************************************************************************
00258                          *
00259                          *  This portion implements Gigante's non-uniform 3-D space grid/mesh discretization.
00260                          *
00261                          ****************************************************************************************
00262                          */
00263                         register int out_axis;
00264                         register const struct nu_axis  **nu_axis =
00265                              (const struct nu_axis **)&curcut->nugn.nu_axis[0];
00266                         register const int              *nu_stepsize =
00267                              &curcut->nugn.nu_stepsize[0];
00268                         register const int              *nu_cells_per_axis =
00269                              &curcut->nugn.nu_cells_per_axis[0];
00270                         register const union cutter     *nu_grid =
00271                              curcut->nugn.nu_grid;
00272 
00273                         if( ssp->lastcell == CUTTER_NULL ) {
00274                                 /* We have just started into this NUgrid.  We
00275                                    must find our location and set up the
00276                                    NUgrid traversal state variables. */
00277                                 register int x, y, z;
00278 
00279                                 px = ap->a_ray.r_pt[X] + t0*ap->a_ray.r_dir[X];
00280                                 py = ap->a_ray.r_pt[Y] + t0*ap->a_ray.r_dir[Y];
00281                                 pz = ap->a_ray.r_pt[Z] + t0*ap->a_ray.r_dir[Z];
00282 
00283                                 /* Must find cell that contains newray.r_pt.
00284                                    We do this by binary subdivision.
00285                                    If any are out of bounds, we have left the
00286                                    NUgrid and will pop a level off the stack
00287                                    in the outer loop (if applicable).  */
00288                                 x = rt_find_nugrid( &curcut->nugn, X, px );
00289                                 if( x<0 ) break;
00290                                 y = rt_find_nugrid( &curcut->nugn, Y, py );
00291                                 if( y<0 ) break;
00292                                 z = rt_find_nugrid( &curcut->nugn, Z, pz );
00293                                 if( z<0 ) break;
00294 
00295                                 cutp = &nu_grid[z*nu_stepsize[Z] +
00296                                                 y*nu_stepsize[Y] +
00297                                                 x*nu_stepsize[X]];
00298 
00299                                 ssp->igrid[X] = x;
00300                                 ssp->igrid[Y] = y;
00301                                 ssp->igrid[Z] = z;
00302 
00303                                 NUGRID_T_SETUP( X, px, x );
00304                                 NUGRID_T_SETUP( Y, py, y );
00305                                 NUGRID_T_SETUP( Z, pz, z );
00306                         } else {
00307                                 /* Advance from previous cell to next cell */
00308                                 /* Take next step, finding ray entry distance*/
00309                                 cutp = ssp->lastcell;
00310                                 out_axis = ssp->out_axis;
00311 
00312                                 /* We may be simply advancing to the next box
00313                                    in the *same* NUgrid cell (if, for instance,
00314                                    the NUgrid cell is a cutnode with
00315                                    boxnode leaves).  So if t0 hasn't advanced
00316                                    past the end of the box, advance
00317                                    a tiny bit (less than rt_ct_optim makes
00318                                    boxnodes) and be handled by tree-traversing
00319                                    code below. */
00320 
00321                                 if( cutp->cut_type == CUT_CUTNODE &&
00322                                     t0 + OFFSET_DIST < ssp->tv[out_axis] ) {
00323                                         t0 += OFFSET_DIST;
00324                                         break;
00325                                 }
00326 
00327                                 /* Advance to the next cell as appropriate,
00328                                    bailing out with cutp=CUTTER_NULL
00329                                    if we run past the end of the NUgrid
00330                                    array. */
00331 
00332 again:                          t0 = ssp->tv[out_axis];
00333                                 if( ssp->rstep[out_axis] > 0 ) {
00334                                         if( ++(ssp->igrid[out_axis]) >=
00335                                             nu_cells_per_axis[out_axis] ) {
00336                                                 cutp = CUTTER_NULL;
00337                                                 break;
00338                                         }
00339                                         cutp += nu_stepsize[out_axis];
00340                                 } else {
00341                                         if( --(ssp->igrid[out_axis]) < 0 ) {
00342                                                 cutp = CUTTER_NULL;
00343                                                 break;
00344                                         }
00345                                         cutp -= nu_stepsize[out_axis];
00346                                 }
00347 
00348                                 /* Update t advancement value for this axis */
00349                                 NUGRID_T_ADV( out_axis, ssp->igrid[out_axis] );
00350                         }
00351 
00352                         /* find minimum exit t value */
00353                         if( ssp->tv[X] < ssp->tv[Y] )  {
00354                                 if( ssp->tv[Z] < ssp->tv[X] )  {
00355                                         out_axis = Z;
00356                                 } else {
00357                                         out_axis = X;
00358                                 }
00359                         } else {
00360                                 if( ssp->tv[Z] < ssp->tv[Y] )  {
00361                                         out_axis = Z;
00362                                 } else {
00363                                         out_axis = Y;
00364                                 }
00365                         }
00366 
00367                         /* Zip through empty cells. */
00368                         if( cutp->cut_type == CUT_BOXNODE &&
00369                             cutp->bn.bn_len <= 0 &&
00370                             cutp->bn.bn_piecelen <= 0 ) {
00371                                 ++ssp->resp->re_nempty_cells;
00372                                 goto again;
00373                         }
00374 
00375                         ssp->out_axis = out_axis;
00376                         ssp->lastcell = cutp;
00377 
00378                         if( RT_G_DEBUG&DEBUG_ADVANCE )
00379                                 bu_log( "t0=%g found in cell (%d,%d,%d), out_axis=%c at %g; cell is %s\n",
00380                                         t0,
00381                                         ssp->igrid[X], ssp->igrid[Y],
00382                                         ssp->igrid[Z],
00383                                         "XYZ*"[ssp->out_axis],
00384                                         ssp->tv[out_axis],
00385                                         cutp->cut_type==CUT_CUTNODE?"cut":
00386                                         cutp->cut_type==CUT_BOXNODE?"box":
00387                                         cutp->cut_type==CUT_NUGRIDNODE?"nu":
00388                                         "?" );
00389                         break; }
00390                 case CUT_CUTNODE:
00391                         /* fall through */
00392                 case CUT_BOXNODE:
00393                         /*
00394                          ****************************************************************************************
00395                          *
00396                          *  This portion implements Muuss' non-uniform binary space partitioning tree.
00397                          *
00398                          ****************************************************************************************
00399                          */
00400                         t0 += OFFSET_DIST;
00401                         cutp = curcut;
00402                         break;
00403                 default:
00404                         rt_bomb(
00405                        "rt_advance_to_next_cell: unknown high-level cutnode" );
00406                 }
00407 
00408                 if( cutp==CUTTER_NULL ) {
00409 pop_space_stack:
00410                         /*
00411                          *  Pop the stack of nested space partitioning methods.
00412                          *  Move up out of the current node, or return if there
00413                          *  is nothing left to do.
00414                          */
00415                         {
00416                                 register struct rt_shootray_status *old = ssp->old_status;
00417 
00418                                 if( old == NULL ) goto escaped_from_model;
00419                                 *ssp = *old;            /* struct copy */
00420                                 bu_free( old, "old rt_shootray_status" );
00421                                 curcut = ssp->curcut;
00422                                 cutp = CUTTER_NULL;
00423                                 continue;
00424                         }
00425                 }
00426 
00427                 /* Compute position and bail if we're outside of the current level. */
00428                 px = ap->a_ray.r_pt[X] + t0*ap->a_ray.r_dir[X];
00429                 py = ap->a_ray.r_pt[Y] + t0*ap->a_ray.r_dir[Y];
00430                 pz = ap->a_ray.r_pt[Z] + t0*ap->a_ray.r_dir[Z];
00431 
00432                 /* Optimization: when it's a boxnode in a nugrid, just return. */
00433                 if( cutp->cut_type == CUT_BOXNODE &&
00434                     curcut->cut_type == CUT_NUGRIDNODE ) {
00435                         ssp->newray.r_pt[X] = px;
00436                         ssp->newray.r_pt[Y] = py;
00437                         ssp->newray.r_pt[Z] = pz;
00438                         ssp->newray.r_min = 0;
00439                         ssp->newray.r_max = ssp->tv[ssp->out_axis] - t0;
00440                         goto done_return_cutp;
00441                 }
00442 
00443                 /*  Given direction of travel, see if point is outside bound.
00444                  *  This will be the model RPP for NUBSP.
00445                  */
00446                 if( PT_DEPARTING_RPP( ssp->rstep, ssp->curmin, ssp->curmax, px, py, pz ) )
00447                         goto pop_space_stack;
00448 
00449                 if( RT_G_DEBUG&DEBUG_ADVANCE ) {
00450                         bu_log(
00451                    "rt_advance_to_next_cell() dist_corr=%g, pt=(%g, %g, %g)\n",
00452                                 t0 /*ssp->dist_corr*/, px, py, pz );
00453                 }
00454 
00455                 while( cutp->cut_type == CUT_CUTNODE ) {
00456                         switch( cutp->cn.cn_axis )  {
00457                         case X:
00458                                 if( !(px < cutp->cn.cn_point) )  {
00459                                         cutp=cutp->cn.cn_r;
00460                                 }  else  {
00461                                         cutp=cutp->cn.cn_l;
00462                                 }
00463                                 break;
00464                         case Y:
00465                                 if( !(py < cutp->cn.cn_point) )  {
00466                                         cutp=cutp->cn.cn_r;
00467                                 }  else  {
00468                                         cutp=cutp->cn.cn_l;
00469                                 }
00470                                 break;
00471                         case Z:
00472                                 if( !(pz < cutp->cn.cn_point) )  {
00473                                         cutp=cutp->cn.cn_r;
00474                                 }  else  {
00475                                         cutp=cutp->cn.cn_l;
00476                                 }
00477                                 break;
00478                         }
00479                 }
00480 
00481                 if( cutp == CUTTER_NULL )
00482                         rt_bomb( "rt_advance_to_next_cell: leaf is NULL!" );
00483 
00484                 switch( cutp->cut_type ) {
00485                 case CUT_BOXNODE:
00486                         if( RT_G_DEBUG&DEBUG_ADVANCE &&
00487                             PT_DEPARTING_RPP( ssp->rstep, ssp->curmin, ssp->curmax, px, py, pz )
00488                         ) {
00489                                 /* This cell is old news. */
00490                                 bu_log(
00491           "rt_advance_to_next_cell(): point not in cell, advancing\n   pt (%.20e,%.20e,%.20e)\n",
00492                                         px, py, pz );
00493                                 bu_log( "  min (%.20e,%.20e,%.20e)\n",
00494                                         V3ARGS(cutp->bn.bn_min) );
00495                                 bu_log( "  max (%.20e,%.20e,%.20e)\n",
00496                                         V3ARGS(cutp->bn.bn_max) );
00497                                 bu_log( "pt=(%g,%g,%g)\n",px, py, pz );
00498                                 rt_pr_cut( cutp, 0 );
00499 
00500                                 /*
00501                                  * Move newray point further into new box.
00502                                  * Try again.
00503                                  */
00504                                 t0 += OFFSET_DIST;
00505                                 goto top;
00506                         }
00507                         /* Don't get stuck within the same box for long */
00508                         if( cutp==ssp->lastcut ) {
00509                                 fastf_t delta;
00510 push_to_next_box:                               ;
00511                                 if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00512                                         bu_log(
00513                                                 "%d,%d box push odist_corr=%.20e n=%.20e model_end=%.20e\n",
00514                                                 ap->a_x, ap->a_y,
00515                                                 ssp->odist_corr,
00516                                                 t0 /*ssp->dist_corr*/,
00517                                                 ssp->model_end );
00518                                         bu_log(
00519                                                 "box_start o=%.20e n=%.20e\nbox_end   o=%.20e n=%.20e\n",
00520                                                 ssp->obox_start,
00521                                                 ssp->box_start,
00522                                                 ssp->obox_end,
00523                                                 ssp->box_end );
00524                                         bu_log( "Point=(%g,%g,%g)\n",
00525                                                 px, py, pz );
00526                                         VPRINT( "Dir",
00527                                                 ssp->newray.r_dir );
00528                                         rt_pr_cut( cutp, 0 );
00529                                 }
00530 
00531                                 /* Advance 1mm, or smallest value that hardware
00532                                  * floating point resolution will allow.
00533                                  */
00534                                 fraction = frexp( ssp->box_end,
00535                                                   &exponent );
00536 
00537                                 if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00538                                         bu_log(
00539                                                 "exp=%d, fraction=%.20e\n",
00540                                                 exponent, fraction );
00541                                 }
00542                                 if( sizeof(fastf_t) <= 4 )
00543                                         fraction += 1.0e-5;
00544                                 else
00545                                         fraction += 1.0e-14;
00546                                 delta = ldexp( fraction, exponent );
00547                                 if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00548                                         bu_log(
00549                                                 "ldexp: delta=%g, fract=%g, exp=%d\n",
00550                                                 delta,
00551                                                 fraction,
00552                                                 exponent );
00553                                 }
00554 
00555                                 /* Never advance less than 1mm */
00556                                 if( delta < 1 ) delta = 1.0;
00557                                 ssp->box_start = ssp->box_end + delta;
00558                                 ssp->box_end = ssp->box_start + delta;
00559 
00560                                 if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00561                                         bu_log(
00562                                                 "push%d: was=%.20e, now=%.20e\n\n",
00563                                                 push_flag,
00564                                                 ssp->box_end,
00565                                                 ssp->box_start );
00566                                 }
00567                                 push_flag++;
00568                                 if( push_flag > 3 ) {
00569                                         bu_log( "rt_advance_to_next_cell(): INTERNAL ERROR: infinite loop aborted, ray %d,%d truncated\n",
00570                                                 ap->a_x, ap->a_y );
00571                                         goto pop_space_stack;
00572                                 }
00573                                 /* See if point marched outside model RPP */
00574                                 if( ssp->box_start > ssp->model_end )
00575                                         goto pop_space_stack;
00576                                 t0 = ssp->box_start + OFFSET_DIST;
00577                                 goto top;
00578                         }
00579                         if( push_flag ) {
00580                                 push_flag = 0;
00581                                 if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00582                                         bu_log(
00583                                                 "%d,%d Escaped %d. dist_corr=%g, box_start=%g, box_end=%g\n",
00584                                                 ap->a_x, ap->a_y,
00585                                                 push_flag,
00586                                                 t0 /*ssp->dist_corr*/,
00587                                                 ssp->box_start,
00588                                                 ssp->box_end );
00589                                 }
00590                         }
00591                         if( RT_G_DEBUG & DEBUG_ADVANCE ) {
00592                                 bu_log(
00593                                         "rt_advance_to_next_cell()=x%x lastcut=x%x\n",
00594                                         cutp, ssp->lastcut);
00595                         }
00596 
00597                         ssp->newray.r_pt[X] = px;
00598                         ssp->newray.r_pt[Y] = py;
00599                         ssp->newray.r_pt[Z] = pz;
00600                         if( !rt_in_rpp( &ssp->newray, ssp->inv_dir,
00601                                         cutp->bn.bn_min,
00602                                         cutp->bn.bn_max) )  {
00603                                 bu_log("rt_advance_to_next_cell():  MISSED BOX\nrmin,rmax(%.20e,%.20e) box(%.20e,%.20e)\n",
00604                                        ssp->newray.r_min,
00605                                        ssp->newray.r_max,
00606                                        ssp->box_start, ssp->box_end );
00607                                 goto push_to_next_box;
00608                         }
00609 
00610 done_return_cutp:       ssp->lastcut = cutp;
00611 #if EXTRA_SAFETY
00612                         /* Diagnostic purposes only */
00613                         ssp->odist_corr = ssp->dist_corr;
00614                         ssp->obox_start = ssp->box_start;
00615                         ssp->obox_end = ssp->box_end;
00616 #endif
00617                         ssp->dist_corr = t0;
00618                         ssp->box_start = t0 + ssp->newray.r_min;
00619                         ssp->box_end = t0 + ssp->newray.r_max;
00620                         if( RT_G_DEBUG & DEBUG_ADVANCE )  {
00621                                 bu_log(
00622                                 "rt_advance_to_next_cell() box=(%g, %g)\n",
00623                                         ssp->box_start, ssp->box_end );
00624                         }
00625                         return cutp;
00626                 case CUT_NUGRIDNODE: {
00627                         struct rt_shootray_status *old;
00628 
00629                         BU_GETSTRUCT( old, rt_shootray_status );
00630                         *old = *ssp;    /* struct copy */
00631 
00632                         /* Descend into node */
00633                         ssp->old_status = old;
00634                         ssp->lastcut = ssp->lastcell = CUTTER_NULL;
00635                         curcut = ssp->curcut = cutp;
00636                         VSET( ssp->curmin, cutp->nugn.nu_axis[X][0].nu_spos,
00637                                            cutp->nugn.nu_axis[Y][0].nu_spos,
00638                                            cutp->nugn.nu_axis[Z][0].nu_spos );
00639                         VSET( ssp->curmax, ssp->curcut->nugn.nu_axis[X][ssp->curcut->nugn.nu_cells_per_axis[X]-1].nu_epos,
00640                                  ssp->curcut->nugn.nu_axis[Y][ssp->curcut->nugn.nu_cells_per_axis[Y]-1].nu_epos,
00641                                  ssp->curcut->nugn.nu_axis[Z][ssp->curcut->nugn.nu_cells_per_axis[Z]-1].nu_epos );
00642                         break; }
00643                 case CUT_CUTNODE:
00644                         rt_bomb( "rt_advance_to_next_cell: impossible: cutnote as leaf!" );
00645                         break;
00646                 default:
00647                         rt_bomb( "rt_advance_to_next_cell: unknown spt node" );
00648                         break;
00649                 }
00650 
00651                 /* Continue with the current space partitioning algorithm. */
00652         }
00653         /* NOTREACHED */
00654         /*      bu_bomb("rt_advance_to_next_cell: escaped for(;;) loop: impossible!"); */
00655 
00656         /*
00657          *  If ray has escaped from model RPP, and there are infinite solids
00658          *  in the model, there is one more (special) BOXNODE for the
00659          *  caller to process.
00660          */
00661 escaped_from_model:
00662         curcut = &ssp->ap->a_rt_i->rti_inf_box;
00663         if( curcut->bn.bn_len <= 0 && curcut->bn.bn_piecelen <= 0 )
00664                 curcut = CUTTER_NULL;
00665         ssp->curcut = curcut;
00666         return curcut;
00667 }
00668 
00669 /*              R T _ F I N D _ B A C K I N G _ D I S T
00670  *
00671  *      This routine traces a ray from its start point to model exit through the space partitioning tree.
00672  *      The objective is to find all primitives that use "pieces" and also have a bounding box that
00673  *      extends behind the ray start point. The minimum (most negative) intersection with such a bounding
00674  *      box is returned. The "backbits" bit vector (provided by the caller) gets a bit set for every
00675  *      primitive that meets the above criteria.
00676  *      No primitive intersections are performed.
00677  *
00678  *      XXXX This routine does not work for NUGRID XXXX
00679  */
00680 fastf_t
00681 rt_find_backing_dist( struct rt_shootray_status *ss, struct bu_bitv *backbits ) {
00682         fastf_t min_backing_dist=BACKING_DIST;
00683         fastf_t cur_dist=0.0;
00684         point_t cur_pt;
00685         union cutter *cutp;
00686         struct bu_bitv *solidbits;
00687         struct xray ray;
00688         struct resource *resp;
00689         struct rt_i *rtip;
00690         int i;
00691 
00692         resp = ss->ap->a_resource;
00693         rtip = ss->ap->a_rt_i;
00694 
00695         /* get a bit vector of our own to avoid duplicate bounding box intersection calculations */
00696         solidbits = get_solidbitv( rtip->nsolids, resp );
00697         bu_bitv_clear(solidbits);
00698 
00699         ray = ss->ap->a_ray;    /* struct copy, don't mess with the original */
00700 
00701         /* cur_dist keeps track of where we are along the ray */
00702         /* stop when cur_dist reaches far intersection of ray and model bounding box */
00703         while( cur_dist <= ss->ap->a_ray.r_max ) {
00704                 /* calculate the current point along the ray */
00705                 VJOIN1( cur_pt, ss->ap->a_ray.r_pt, cur_dist, ss->ap->a_ray.r_dir );
00706 
00707                 /* descend into the space partitioning tree based on this point */
00708                 cutp = &ss->ap->a_rt_i->rti_CutHead;
00709                 while( cutp->cut_type == CUT_CUTNODE ) {
00710                         if( cur_pt[cutp->cn.cn_axis] >= cutp->cn.cn_point )  {
00711                                 cutp=cutp->cn.cn_r;
00712                         }  else  {
00713                                 cutp=cutp->cn.cn_l;
00714                         }
00715                 }
00716 
00717                 /* we are now at the box node for the current point */
00718                 /* check if the ray intersects this box */
00719                 if( !rt_in_rpp( &ray, ss->inv_dir, cutp->bn.bn_min, cutp->bn.bn_max ) ) {
00720                         /* ray does not intersect this cell
00721                          * one of two situations must exist:
00722                          *      1. ray starts outside model bounding box (no need for these calculations)
00723                          *      2. we have proceeded beyond end of model bounding box (we are done)
00724                          * in either case, we are finished
00725                          */
00726                         goto done;
00727                 } else {
00728                         /* increment cur_dist into next cell for next execution of this loop */
00729                         cur_dist = ray.r_max + OFFSET_DIST;
00730                 }
00731 
00732                 /* process this box node (look at all the pieces) */
00733                 for( i=0 ; i<cutp->bn.bn_piecelen ; i++ ) {
00734                         struct rt_piecelist *plp=&cutp->bn.bn_piecelist[i];
00735 
00736                         if( BU_BITTEST( solidbits, plp->stp->st_bit ) == 0 ) {
00737                                 /* we haven't looked at this primitive before */
00738                                 if( rt_in_rpp( &ray, ss->inv_dir, plp->stp->st_min, plp->stp->st_max ) ) {
00739                                         /* ray intersects this primitive bounding box */
00740 
00741                                         if( ray.r_min < BACKING_DIST ) {
00742                                                 if( ray.r_min < min_backing_dist ) {
00743                                                         /* move our backing distance back to catch this one */
00744                                                         min_backing_dist = ray.r_min;
00745                                                 }
00746 
00747                                                 /* add this one to our list of primitives to check */
00748                                                 BU_BITSET( backbits, plp->stp->st_bit );
00749                                         }
00750                                 }
00751                                 /* set bit so we don't repeat this calculation */
00752                                 BU_BITSET( solidbits, plp->stp->st_bit );
00753                         }
00754                 }
00755         }
00756  done:
00757         /* put our bit vector on the resource list */
00758         BU_CK_BITV(solidbits);
00759         BU_LIST_APPEND( &resp->re_solid_bitv, &solidbits->l );
00760 
00761         /* return our minimum backing distance */
00762         return min_backing_dist;
00763 }
00764 
00765 /*
00766  *                      R T _ S H O O T R A Y
00767  *
00768  *  Note that the direction vector r_dir
00769  *  must have unit length;  this is mandatory, and is not ordinarily checked,
00770  *  in the name of efficiency.
00771  *
00772  *  Input:  Pointer to an application structure, with these mandatory fields:
00773  *      a_ray.r_pt      Starting point of ray to be fired
00774  *      a_ray.r_dir     UNIT VECTOR with direction to fire in (dir cosines)
00775  *      a_hit           Routine to call when something is hit
00776  *      a_miss          Routine to call when ray misses everything
00777  *
00778  *  Calls user's a_miss() or a_hit() routine as appropriate.
00779  *  Passes a_hit() routine list of partitions, with only hit_dist
00780  *  fields valid.  Normal computation deferred to user code,
00781  *  to avoid needless computation here.
00782  *
00783  *  Formal Return: whatever the application function returns (an int).
00784  *
00785  *  NOTE:  The appliction functions may call rt_shootray() recursively.
00786  *      Thus, none of the local variables may be static.
00787  *
00788  *  To prevent having to lock the statistics variables in a PARALLEL
00789  *  environment, all the statistics variables have been moved into
00790  *  the 'resource' structure, which is allocated per-CPU.
00791  */
00792 
00793 int
00794 rt_shootray(register struct application *ap)
00795 {
00796         struct rt_shootray_status       ss;
00797         struct seg              new_segs;       /* from solid intersections */
00798         struct seg              waiting_segs;   /* awaiting rt_boolweave() */
00799         struct seg              finished_segs;  /* processed by rt_boolweave() */
00800         fastf_t                 last_bool_start;
00801         struct bu_bitv          *solidbits;     /* bits for all solids shot so far */
00802         struct bu_bitv          *backbits=NULL; /* bits for all solids using pieces that need to be intersected behind
00803                                                    the the ray start point */
00804         struct bu_ptbl          *regionbits;    /* table of all involved regions */
00805         char                    *status;
00806         auto struct partition   InitialPart;    /* Head of Initial Partitions */
00807         auto struct partition   FinalPart;      /* Head of Final Partitions */
00808         struct soltab           **stpp;
00809         register const union cutter *cutp;
00810         struct resource         *resp;
00811         struct rt_i             *rtip;
00812         const int               debug_shoot = RT_G_DEBUG & DEBUG_SHOOT;
00813         fastf_t                 pending_hit = 0; /* dist of closest odd hit pending */
00814 
00815         RT_AP_CHECK(ap);
00816         if( ap->a_magic )  {
00817                 RT_CK_AP(ap);
00818         } else {
00819                 ap->a_magic = RT_AP_MAGIC;
00820         }
00821         if( ap->a_ray.magic )  {
00822                 RT_CK_RAY(&(ap->a_ray));
00823         } else {
00824                 ap->a_ray.magic = RT_RAY_MAGIC;
00825         }
00826         if( ap->a_resource == RESOURCE_NULL )  {
00827                 ap->a_resource = &rt_uniresource;
00828                 if(RT_G_DEBUG)bu_log("rt_shootray:  defaulting a_resource to &rt_uniresource\n");
00829                 if( rt_uniresource.re_magic == 0 )
00830                         rt_init_resource( &rt_uniresource, 0, ap->a_rt_i );
00831         }
00832         ss.ap = ap;
00833         rtip = ap->a_rt_i;
00834         RT_CK_RTI( rtip );
00835         resp = ap->a_resource;
00836         RT_CK_RESOURCE(resp);
00837         ss.resp = resp;
00838 
00839         if(RT_G_DEBUG) {
00840             /* only test extensively if something in run-time debug is enabled */
00841             if (RT_G_DEBUG & (DEBUG_ALLRAYS|DEBUG_SHOOT|DEBUG_PARTITION|DEBUG_ALLHITS)) {
00842                 bu_log_indent_delta(2);
00843                 bu_log("\n**********shootray cpu=%d  %d,%d lvl=%d a_onehit=%d (%s)\n",
00844                         resp->re_cpu,
00845                         ap->a_x, ap->a_y,
00846                         ap->a_level,
00847                         ap->a_onehit,
00848                         ap->a_purpose != (char *)0 ? ap->a_purpose : "?" );
00849                 bu_log("Pnt (%.20G, %.20G, %.20G)\nDir (%.20G, %.20G, %.20G)\n",
00850                         V3ARGS(ap->a_ray.r_pt),
00851                         V3ARGS(ap->a_ray.r_dir) );
00852             }
00853         }
00854 #ifndef NO_BADRAY_CHECKING
00855         if(RT_BADVEC(ap->a_ray.r_pt)||RT_BADVEC(ap->a_ray.r_dir))  {
00856                 bu_log("\n**********shootray cpu=%d  %d,%d lvl=%d (%s)\n",
00857                         resp->re_cpu,
00858                         ap->a_x, ap->a_y,
00859                         ap->a_level,
00860                         ap->a_purpose != (char *)0 ? ap->a_purpose : "?" );
00861                 VPRINT(" r_pt", ap->a_ray.r_pt);
00862                 VPRINT("r_dir", ap->a_ray.r_dir);
00863                 rt_bomb("rt_shootray() bad ray\n");
00864         }
00865 #endif
00866 
00867         if( rtip->needprep )
00868                 rt_prep_parallel(rtip, 1);      /* Stay on our CPU */
00869 
00870         InitialPart.pt_forw = InitialPart.pt_back = &InitialPart;
00871         InitialPart.pt_magic = PT_HD_MAGIC;
00872         FinalPart.pt_forw = FinalPart.pt_back = &FinalPart;
00873         FinalPart.pt_magic = PT_HD_MAGIC;
00874         ap->a_Final_Part_hdp = &FinalPart;
00875 
00876         BU_LIST_INIT( &new_segs.l );
00877         BU_LIST_INIT( &waiting_segs.l );
00878         BU_LIST_INIT( &finished_segs.l );
00879         ap->a_finished_segs_hdp = &finished_segs;
00880 
00881         if( BU_LIST_UNINITIALIZED( &resp->re_parthead ) )  {
00882                 /* XXX This shouldn't happen any more */
00883                 bu_log("rt_shootray() resp=x%x uninitialized, fixing it\n", resp);
00884                 /*
00885                  *  We've been handed a mostly un-initialized resource struct,
00886                  *  with only a magic number and a cpu number filled in.
00887                  *  Init it and add it to the table.
00888                  *  This is how application-provided resource structures
00889                  *  are remembered for later cleanup by the library.
00890                  */
00891                 rt_init_resource( resp, resp->re_cpu, rtip );
00892         }
00893         /* Ensure that this CPU's resource structure is registered */
00894         if( resp != &rt_uniresource )
00895                 BU_ASSERT_PTR( BU_PTBL_GET(&rtip->rti_resources, resp->re_cpu), !=, NULL );
00896 
00897         solidbits = get_solidbitv( rtip->nsolids, resp );
00898         bu_bitv_clear(solidbits);
00899 
00900         if( BU_LIST_IS_EMPTY( &resp->re_region_ptbl ) )  {
00901                 BU_GETSTRUCT( regionbits, bu_ptbl );
00902                 bu_ptbl_init( regionbits, 7, "rt_shootray() regionbits ptbl" );
00903         } else {
00904                 regionbits = BU_LIST_FIRST( bu_ptbl, &resp->re_region_ptbl );
00905                 BU_LIST_DEQUEUE( &regionbits->l );
00906                 BU_CK_PTBL(regionbits);
00907         }
00908 
00909         if( !resp->re_pieces && rtip->rti_nsolids_with_pieces > 0 )  {
00910                 /* Initialize this processors 'solid pieces' state */
00911                 rt_res_pieces_init( resp, rtip );
00912         }
00913         if( BU_LIST_MAGIC_WRONG( &resp->re_pieces_pending.l, BU_PTBL_MAGIC ) )
00914                 bu_ptbl_init( &resp->re_pieces_pending, 100, "re_pieces_pending" );
00915         bu_ptbl_reset( &resp->re_pieces_pending );
00916 
00917         /* Verify that direction vector has unit length */
00918         if(RT_G_DEBUG) {
00919                 FAST fastf_t f, diff;
00920                 /* Fancy version of BN_VEC_NON_UNIT_LEN() */
00921                 f = MAGSQ(ap->a_ray.r_dir);
00922                 if( NEAR_ZERO(f, 0.0001) )  {
00923                         rt_bomb("rt_shootray:  zero length dir vector\n");
00924                         return(0);
00925                 }
00926                 diff = f - 1;
00927                 if( !NEAR_ZERO( diff, 0.0001 ) )  {
00928                         bu_log("rt_shootray: non-unit dir vect (x%d y%d lvl%d)\n",
00929                                 ap->a_x, ap->a_y, ap->a_level );
00930                         f = 1/f;
00931                         VSCALE( ap->a_ray.r_dir, ap->a_ray.r_dir, f );
00932                 }
00933         }
00934 
00935         /*
00936          *  Record essential statistics in per-processor data structure.
00937          */
00938         resp->re_nshootray++;
00939 
00940         /* Compute the inverse of the direction cosines */
00941         if( ap->a_ray.r_dir[X] < -SQRT_SMALL_FASTF )  {
00942                 ss.abs_inv_dir[X] = -(ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
00943                 ss.rstep[X] = -1;
00944         } else if( ap->a_ray.r_dir[X] > SQRT_SMALL_FASTF )  {
00945                 ss.abs_inv_dir[X] =  (ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
00946                 ss.rstep[X] = 1;
00947         } else {
00948                 ap->a_ray.r_dir[X] = 0.0;
00949                 ss.abs_inv_dir[X] = ss.inv_dir[X] = INFINITY;
00950                 ss.rstep[X] = 0;
00951         }
00952         if( ap->a_ray.r_dir[Y] < -SQRT_SMALL_FASTF )  {
00953                 ss.abs_inv_dir[Y] = -(ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
00954                 ss.rstep[Y] = -1;
00955         } else if( ap->a_ray.r_dir[Y] > SQRT_SMALL_FASTF )  {
00956                 ss.abs_inv_dir[Y] =  (ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
00957                 ss.rstep[Y] = 1;
00958         } else {
00959                 ap->a_ray.r_dir[Y] = 0.0;
00960                 ss.abs_inv_dir[Y] = ss.inv_dir[Y] = INFINITY;
00961                 ss.rstep[Y] = 0;
00962         }
00963         if( ap->a_ray.r_dir[Z] < -SQRT_SMALL_FASTF )  {
00964                 ss.abs_inv_dir[Z] = -(ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
00965                 ss.rstep[Z] = -1;
00966         } else if( ap->a_ray.r_dir[Z] > SQRT_SMALL_FASTF )  {
00967                 ss.abs_inv_dir[Z] =  (ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
00968                 ss.rstep[Z] = 1;
00969         } else {
00970                 ap->a_ray.r_dir[Z] = 0.0;
00971                 ss.abs_inv_dir[Z] = ss.inv_dir[Z] = INFINITY;
00972                 ss.rstep[Z] = 0;
00973         }
00974         VMOVE( ap->a_inv_dir, ss.inv_dir );
00975 
00976         /*
00977          *  If ray does not enter the model RPP, skip on.
00978          *  If ray ends exactly at the model RPP, trace it.
00979          */
00980         if( !rt_in_rpp( &ap->a_ray, ss.inv_dir, rtip->mdl_min, rtip->mdl_max )  ||
00981             ap->a_ray.r_max < 0.0 )  {
00982                 cutp = &ap->a_rt_i->rti_inf_box;
00983                 if( cutp->bn.bn_len > 0 )  {
00984                         /* Model has infinite solids, need to fire at them. */
00985                         ss.box_start = BACKING_DIST;
00986                         ss.model_start = 0;
00987                         ss.box_end = ss.model_end = INFINITY;
00988                         ss.lastcut = CUTTER_NULL;
00989                         ss.old_status = (struct rt_shootray_status *)NULL;
00990                         ss.curcut = cutp;
00991                         ss.lastcell = ss.curcut;
00992                         VMOVE( ss.curmin, rtip->mdl_min );
00993                         VMOVE( ss.curmax, rtip->mdl_max );
00994                         last_bool_start = BACKING_DIST;
00995                         ss.newray = ap->a_ray;          /* struct copy */
00996                         ss.odist_corr = ss.obox_start = ss.obox_end = -99;
00997                         ss.dist_corr = 0.0;
00998                         goto start_cell;
00999                 }
01000                 resp->re_nmiss_model++;
01001                 ap->a_return = ap->a_miss( ap );
01002                 status = "MISS model";
01003                 goto out;
01004         }
01005 
01006         /*
01007          *  The interesting part of the ray starts at distance 0.
01008          *  If the ray enters the model at a negative distance,
01009          *  (ie, the ray starts within the model RPP),
01010          *  we only look at little bit behind (BACKING_DIST) to see if we are
01011          *  just coming out of something, but never further back than
01012          *  the intersection with the model RPP.
01013          *  If the ray enters the model at a positive distance,
01014          *  we always start there.
01015          *  It is vital that we never pick a start point outside the
01016          *  model RPP, or the space partitioning tree will pick the
01017          *  wrong box and the ray will miss it.
01018          *
01019          *  BACKING_DIST should probably be determined by floating point
01020          *  noise factor due to model RPP size -vs- number of bits of
01021          *  floating point mantissa significance, rather than a constant,
01022          *  but that is too hideous to think about here.
01023          *  Also note that applications that really depend on knowing
01024          *  what region they are leaving from should probably back their
01025          *  own start-point up, rather than depending on it here, but
01026          *  it isn't much trouble here.
01027          *
01028          *  Modification by JRA for pieces methodology:
01029          *      The original algorithm here assumed that if we encountered any primitive
01030          *      along the positive direction of the ray, ALL its intersections would be calculated.
01031          *      With pieces, we may see only an exit hit if the entrance piece is in a space partition cell
01032          *      that is more than "BACKING_DIST" behind the ray start point (leading to incorrect results).
01033          *      I have modified the setting of "ss.box_start" (when pieces are present and the ray start point is
01034          *      inside the model bounding box) as follows (see rt_find_backing_dist()):
01035          *              The ray is traced through the space partitioning tree
01036          *              The ray is intersected with the bounding box of each primitive using pieces in each cell
01037          *              The minimum of all these intersections is set as the initial "ss.box_start".
01038          *              The "backbits" bit vector has a bit set for each of the primitives using pieces that have
01039          *                      bounding boxes that extend behind the ray start point
01040          *      Further below (in the "pieces" loop), I have added code to ignore primitives that do not have a bit
01041          *      set in the backbits vector when we are behind the ray start point.
01042          */
01043 
01044         /* these two values set the point where the ray tracing actually begins and ends */
01045         ss.box_start = ss.model_start = ap->a_ray.r_min;
01046         ss.box_end = ss.model_end = ap->a_ray.r_max;
01047 
01048         if( ap->a_rt_i->rti_nsolids_with_pieces > 0 ) {
01049                 /* pieces are present */
01050                 if( ss.box_start < BACKING_DIST ) {
01051                         /* the first ray intersection with the model bounding box is more than BACKING_DIST
01052                          * behind the ray start point
01053                          */
01054 
01055                         /* get a bit vector to keep track of which primitives need to be intersected
01056                          * behind the ray start point (those having bounding boxes extending behind the
01057                          * ray start point and using pieces)
01058                          */
01059                         backbits = get_solidbitv( rtip->nsolids, resp );
01060                         bu_bitv_clear(backbits);
01061 
01062                         /* call "rt_find_backing_dist()" to calculate the required
01063                          * start point for calculation, and to fill in the "backbits" bit vector
01064                          */
01065                         ss.box_start = rt_find_backing_dist( &ss, backbits );
01066                 }
01067         } else {
01068                 /* no pieces present, use the old scheme */
01069                 if( ss.box_start < BACKING_DIST )
01070                         ss.box_start = BACKING_DIST; /* Only look a little bit behind */
01071         }
01072 
01073         ss.lastcut = CUTTER_NULL;
01074         ss.old_status = (struct rt_shootray_status *)NULL;
01075         ss.curcut = &ap->a_rt_i->rti_CutHead;
01076         if( ss.curcut->cut_type == CUT_NUGRIDNODE ) {
01077                 ss.lastcell = CUTTER_NULL;
01078                 VSET( ss.curmin, ss.curcut->nugn.nu_axis[X][0].nu_spos,
01079                                  ss.curcut->nugn.nu_axis[Y][0].nu_spos,
01080                                  ss.curcut->nugn.nu_axis[Z][0].nu_spos );
01081                 VSET( ss.curmax, ss.curcut->nugn.nu_axis[X][ss.curcut->nugn.nu_cells_per_axis[X]-1].nu_epos,
01082                                  ss.curcut->nugn.nu_axis[Y][ss.curcut->nugn.nu_cells_per_axis[Y]-1].nu_epos,
01083                                  ss.curcut->nugn.nu_axis[Z][ss.curcut->nugn.nu_cells_per_axis[Z]-1].nu_epos );
01084         } else if( ss.curcut->cut_type == CUT_CUTNODE ||
01085                    ss.curcut->cut_type == CUT_BOXNODE ) {
01086                 ss.lastcell = ss.curcut;
01087                 VMOVE( ss.curmin, rtip->mdl_min );
01088                 VMOVE( ss.curmax, rtip->mdl_max );
01089         }
01090 
01091         last_bool_start = BACKING_DIST;
01092         ss.newray = ap->a_ray;          /* struct copy */
01093         ss.odist_corr = ss.obox_start = ss.obox_end = -99;
01094         ss.dist_corr = 0.0;
01095         ss.box_num = 0;
01096 
01097         /*
01098          *  While the ray remains inside model space,
01099          *  push from box to box until ray emerges from
01100          *  model space again (or first hit is found, if user is impatient).
01101          *  It is vitally important to always stay within the model RPP, or
01102          *  the space partitoning tree will pick wrong boxes & miss them.
01103          */
01104         while( (cutp = rt_advance_to_next_cell( &ss )) != CUTTER_NULL )  {
01105 start_cell:
01106                 if(debug_shoot) {
01107                         bu_log("BOX #%d interval is %g..%g\n", ss.box_num, ss.box_start, ss.box_end);
01108                         rt_pr_cut( cutp, 0 );
01109                 }
01110 
01111                 if( cutp->bn.bn_len <= 0 && cutp->bn.bn_piecelen <= 0 )  {
01112                         /* Push ray onwards to next box */
01113                         ss.box_start = ss.box_end;
01114                         resp->re_nempty_cells++;
01115                         continue;
01116                 }
01117 
01118                 /* Consider all "pieces" of all solids within the box */
01119                 pending_hit = ss.box_end;
01120                 if( cutp->bn.bn_piecelen > 0 )  {
01121                         register struct rt_piecelist *plp;
01122 
01123                         plp = &(cutp->bn.bn_piecelist[cutp->bn.bn_piecelen-1]);
01124                         for( ; plp >= cutp->bn.bn_piecelist; plp-- )  {
01125                                 struct rt_piecestate *psp;
01126                                 struct soltab   *stp;
01127                                 int ret;
01128                                 int had_hits_before;
01129 
01130                                 RT_CK_PIECELIST(plp);
01131 
01132                                 /* Consider all pieces of this one solid in this cell */
01133                                 stp = plp->stp;
01134                                 RT_CK_SOLTAB(stp);
01135 
01136                                 if( backbits && ss.box_end < BACKING_DIST && BU_BITTEST( backbits, stp->st_bit ) == 0 ) {
01137                                         /* we are behind the ray start point and this primitive is not one
01138                                          * that we need to intersect back here
01139                                          */
01140                                         continue;
01141                                 }
01142 
01143                                 psp = &(resp->re_pieces[stp->st_piecestate_num]);
01144                                 RT_CK_PIECESTATE(psp);
01145                                 if( psp->ray_seqno != resp->re_nshootray )  {
01146                                         /* state is from an earlier ray, scrub */
01147                                         BU_BITV_ZEROALL(psp->shot);
01148                                         psp->ray_seqno = resp->re_nshootray;
01149                                         rt_htbl_reset( &psp->htab );
01150 
01151                                         /* Compute ray entry and exit to entire solid's bounding box */
01152                                         if( !rt_in_rpp( &ss.newray, ss.inv_dir,
01153                                             stp->st_min, stp->st_max ) )  {
01154                                                 if(debug_shoot)bu_log("rpp miss %s (all pieces)\n", stp->st_name);
01155                                                 resp->re_prune_solrpp++;
01156                                                 BU_BITSET( solidbits, stp->st_bit );
01157                                                 continue;       /* MISS */
01158                                         }
01159                                         psp->mindist = ss.newray.r_min + ss.dist_corr;
01160                                         psp->maxdist = ss.newray.r_max + ss.dist_corr;
01161                                         if(debug_shoot) bu_log("%s mindist=%g, maxdist=%g\n", stp->st_name, psp->mindist, psp->maxdist);
01162                                         had_hits_before = 0;
01163                                 } else {
01164                                         if( BU_BITTEST( solidbits, stp->st_bit ) )  {
01165                                                 /* we missed the solid RPP in an earlier cell */
01166                                                 resp->re_ndup++;
01167                                                 continue;       /* already shot */
01168                                         }
01169                                         had_hits_before = psp->htab.end;
01170                                 }
01171 
01172                                 /*
01173                                  *  Allow this solid to shoot at all of its
01174                                  *  'pieces' in this cell, all at once.
01175                                  *  'newray' has been transformed to be near
01176                                  *  to this cell, and
01177                                  *  'dist_corr' is the additive correction
01178                                  *  factor that ft_piece_shot() must apply
01179                                  *  to hits calculated using 'newray'.
01180                                  */
01181                                 resp->re_piece_shots++;
01182                                 psp->cutp = cutp;
01183                                 if( (ret = stp->st_meth->ft_piece_shot(
01184                                     psp, plp, ss.dist_corr, &ss.newray, ap, &waiting_segs )) <= 0 )  {
01185                                         /* No hits at all */
01186                                         resp->re_piece_shot_miss++;
01187                                 } else {
01188                                         resp->re_piece_shot_hit++;
01189                                 }
01190                                 if(debug_shoot)bu_log("shooting %s pieces, nhit=%d\n", stp->st_name, ret);
01191 
01192                                 /*  See if this solid has been fully processed yet.
01193                                  *  If ray has passed through bounding volume, we're done.
01194                                  *  ft_piece_hitsegs() will only be called once per ray.
01195                                  */
01196                                 if( ss.box_end > psp->maxdist && psp->htab.end > 0 ) {
01197                                         /* Convert hits into segs */
01198                                         if(debug_shoot)bu_log("shooting %s pieces complete, making segs\n", stp->st_name);
01199                                         /* Distance correction was handled in ft_piece_shot */
01200                                         stp->st_meth->ft_piece_hitsegs( psp, &waiting_segs, ap );
01201                                         rt_htbl_reset( &psp->htab );
01202                                         BU_BITSET( solidbits, stp->st_bit );
01203 
01204                                         if( had_hits_before )
01205                                                 bu_ptbl_rm( &resp->re_pieces_pending, (long *)psp );
01206                                 } else {
01207                                         if( !had_hits_before )
01208                                                 bu_ptbl_ins_unique( &resp->re_pieces_pending, (long *)psp );
01209                                 }
01210                         }
01211                 }
01212 
01213                 /* Consider all solids within the box */
01214                 if( cutp->bn.bn_len > 0 && ss.box_end >= BACKING_DIST )  {
01215                         stpp = &(cutp->bn.bn_list[cutp->bn.bn_len-1]);
01216                         for( ; stpp >= cutp->bn.bn_list; stpp-- )  {
01217                                 register struct soltab *stp = *stpp;
01218 
01219                                 if( BU_BITTEST( solidbits, stp->st_bit ) )  {
01220                                         resp->re_ndup++;
01221                                         continue;       /* already shot */
01222                                 }
01223 
01224                                 /* Shoot a ray */
01225                                 BU_BITSET( solidbits, stp->st_bit );
01226 
01227                                 /* Check against bounding RPP, if desired by solid */
01228                                 if( stp->st_meth->ft_use_rpp )  {
01229                                         if( !rt_in_rpp( &ss.newray, ss.inv_dir,
01230                                             stp->st_min, stp->st_max ) )  {
01231                                                 if(debug_shoot)bu_log("rpp miss %s\n", stp->st_name);
01232                                                 resp->re_prune_solrpp++;
01233                                                 continue;       /* MISS */
01234                                         }
01235                                         if( ss.dist_corr + ss.newray.r_max < BACKING_DIST )  {
01236                                                 if(debug_shoot)bu_log("rpp skip %s, dist_corr=%g, r_max=%g\n", stp->st_name, ss.dist_corr, ss.newray.r_max);
01237                                                 resp->re_prune_solrpp++;
01238                                                 continue;       /* MISS */
01239                                         }
01240                                 }
01241 
01242                                 if(debug_shoot)bu_log("shooting %s\n", stp->st_name);
01243                                 resp->re_shots++;
01244                                 BU_LIST_INIT( &(new_segs.l) );
01245                                 if( stp->st_meth->ft_shot(
01246                                     stp, &ss.newray, ap, &new_segs ) <= 0 )  {
01247                                         resp->re_shot_miss++;
01248                                         continue;       /* MISS */
01249                                 }
01250 
01251                                 /* Add seg chain to list awaiting rt_boolweave() */
01252                                 {
01253                                         register struct seg *s2;
01254                                         while(BU_LIST_WHILE(s2,seg,&(new_segs.l)))  {
01255                                                 BU_LIST_DEQUEUE( &(s2->l) );
01256                                                 /* Restore to original distance */
01257                                                 s2->seg_in.hit_dist += ss.dist_corr;
01258                                                 s2->seg_out.hit_dist += ss.dist_corr;
01259                                                 s2->seg_in.hit_rayp = s2->seg_out.hit_rayp = &ap->a_ray;
01260                                                 BU_LIST_INSERT( &(waiting_segs.l), &(s2->l) );
01261                                         }
01262                                 }
01263                                 resp->re_shot_hit++;
01264                         }
01265                 }
01266                 if( RT_G_DEBUG & DEBUG_ADVANCE )
01267                         rt_plot_cell( cutp, &ss, &(waiting_segs.l), rtip);
01268 
01269                 /*
01270                  *  If a_onehit == 0 and a_ray_length <= 0, then the ray
01271                  *  is traced to +infinity.
01272                  *
01273                  *  If a_onehit != 0, then it indicates how many hit points
01274                  *  (which are greater than the ray start point of 0.0)
01275                  *  the application requires, ie, partitions with inhit >= 0.
01276                  *  (If negative, indicates number of non-air hits needed).
01277                  *  If this box yielded additional segments,
01278                  *  immediately weave them into the partition list,
01279                  *  and perform final boolean evaluation.
01280                  *  If this results in the required number of final
01281                  *  partitions, then cease ray-tracing and hand the
01282                  *  partitions over to the application.
01283                  *  All partitions will have valid in and out distances.
01284                  *  a_ray_length is treated similarly to a_onehit.
01285                  */
01286                 if( BU_LIST_NON_EMPTY( &(waiting_segs.l) ) )  {
01287                         if(debug_shoot)  {
01288                                 struct seg *segp;
01289                                 bu_log("Waiting segs:\n");
01290                                 for( BU_LIST_FOR( segp, seg, &(waiting_segs.l) ) )  {
01291                                         rt_pr_seg(segp);
01292                                 }
01293                         }
01294                         if( ap->a_onehit != 0 )  {
01295                                 int     done;
01296 
01297                                 /* Weave these segments into partition list */
01298                                 rt_boolweave( &finished_segs, &waiting_segs, &InitialPart, ap );
01299 
01300                                 if( BU_PTBL_LEN( &resp->re_pieces_pending ) > 0 )  {
01301                                         /* Find the lowest pending mindist, that's as far as boolfinal can progress to */
01302                                         struct rt_piecestate **psp;
01303                                         for( BU_PTBL_FOR( psp, (struct rt_piecestate **), &resp->re_pieces_pending ) )  {
01304                                                 FAST fastf_t dist;
01305 
01306                                                 dist = (*psp)->mindist;
01307                                                 BU_ASSERT_DOUBLE( dist, <, INFINITY );
01308                                                 if( dist < pending_hit )  {
01309                                                         pending_hit = dist;
01310                                                         if(debug_shoot) bu_log("pending_hit lowered to %g by %s\n", pending_hit, (*psp)->stp->st_name);
01311                                                 }
01312                                         }
01313                                 }
01314 
01315                                 /* Evaluate regions upto end of good segs */
01316                                 if( ss.box_end < pending_hit )  pending_hit = ss.box_end;
01317                                 done = rt_boolfinal( &InitialPart, &FinalPart,
01318                                         last_bool_start, pending_hit, regionbits, ap, solidbits );
01319                                 last_bool_start = pending_hit;
01320 
01321                                 /* See if enough partitions have been acquired */
01322                                 if( done > 0 )  goto hitit;
01323                         }
01324                 }
01325 
01326                 if( ap->a_ray_length > 0.0 &&
01327                     ss.box_end >= ap->a_ray_length &&
01328                     ap->a_ray_length < pending_hit )
01329                         goto weave;
01330 
01331                 /* Push ray onwards to next box */
01332                 ss.box_start = ss.box_end;
01333         }
01334 
01335         /*
01336          *  Ray has finally left known space --
01337          *  Weave any remaining segments into the partition list.
01338          */
01339 weave:
01340         if( RT_G_DEBUG&DEBUG_ADVANCE )
01341                 bu_log( "rt_shootray: ray has left known space\n" );
01342 
01343         /* Process any pending hits into segs */
01344         if( BU_PTBL_LEN( &resp->re_pieces_pending ) > 0 )  {
01345                 struct rt_piecestate **psp;
01346                 for( BU_PTBL_FOR( psp, (struct rt_piecestate **), &resp->re_pieces_pending ) )  {
01347                         if( (*psp)->htab.end > 0 )  {
01348                                 /* Convert any pending hits into segs */
01349                                 /* Distance correction was handled in ft_piece_shot */
01350                                 (*psp)->stp->st_meth->ft_piece_hitsegs( *psp, &waiting_segs, ap );
01351                                 rt_htbl_reset( &(*psp)->htab );
01352                         }
01353                         *psp = NULL;
01354                 }
01355                 bu_ptbl_reset( &resp->re_pieces_pending );
01356         }
01357 
01358         if( BU_LIST_NON_EMPTY( &(waiting_segs.l) ) )  {
01359                 rt_boolweave( &finished_segs, &waiting_segs, &InitialPart, ap );
01360         }
01361 
01362         /* finished_segs chain now has all segments hit by this ray */
01363         if( BU_LIST_IS_EMPTY( &(finished_segs.l) ) )  {
01364                 ap->a_return = ap->a_miss( ap );
01365                 status = "MISS prims";
01366                 goto out;
01367         }
01368 
01369         /*
01370          *  All intersections of the ray with the model have
01371          *  been computed.  Evaluate the boolean trees over each partition.
01372          */
01373         (void)rt_boolfinal( &InitialPart, &FinalPart, BACKING_DIST,
01374                 INFINITY,
01375                 regionbits, ap, solidbits);
01376 
01377         if( FinalPart.pt_forw == &FinalPart )  {
01378                 ap->a_return = ap->a_miss( ap );
01379                 status = "MISS bool";
01380                 RT_FREE_PT_LIST( &InitialPart, resp );
01381                 RT_FREE_SEG_LIST( &finished_segs, resp );
01382                 goto out;
01383         }
01384 
01385         /*
01386          *  Ray/model intersections exist.  Pass the list to the
01387          *  user's a_hit() routine.  Note that only the hit_dist
01388          *  elements of pt_inhit and pt_outhit have been computed yet.
01389          *  To compute both hit_point and hit_normal, use the
01390          *
01391          *      RT_HIT_NORM( hitp, stp, rayp )
01392          *
01393          *  macro.  To compute just hit_point, use
01394          *
01395          *  VJOIN1( hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir );
01396          */
01397 hitit:
01398         if(debug_shoot)  rt_pr_partitions(rtip,&FinalPart,"a_hit()");
01399 
01400         /*
01401          *  Before recursing, release storage for unused Initial partitions.
01402          *  finished_segs can not be released yet, because FinalPart
01403          *  partitions will point to hits in those segments.
01404          */
01405         RT_FREE_PT_LIST( &InitialPart, resp );
01406 
01407         /*
01408          *  finished_segs is only used by special hit routines
01409          *  which don't follow the traditional solid modeling paradigm.
01410          */
01411         if(RT_G_DEBUG&DEBUG_ALLHITS) rt_pr_partitions(rtip,&FinalPart,"Partition list passed to a_hit() routine");
01412         ap->a_return = ap->a_hit( ap, &FinalPart, &finished_segs );
01413         status = "HIT";
01414 
01415         RT_FREE_SEG_LIST( &finished_segs, resp );
01416         RT_FREE_PT_LIST( &FinalPart, resp );
01417 
01418         /*
01419          * Processing of this ray is complete.
01420          */
01421 out:
01422         /*  Return dynamic resources to their freelists.  */
01423         BU_CK_BITV(solidbits);
01424         BU_LIST_APPEND( &resp->re_solid_bitv, &solidbits->l );
01425         if( backbits ) {
01426                 BU_CK_BITV(backbits);
01427                 BU_LIST_APPEND( &resp->re_solid_bitv, &backbits->l );
01428         }
01429         BU_CK_PTBL(regionbits);
01430         BU_LIST_APPEND( &resp->re_region_ptbl, &regionbits->l );
01431 
01432         /* Clean up any pending hits */
01433         if( BU_PTBL_LEN( &resp->re_pieces_pending ) > 0 )  {
01434                 struct rt_piecestate **psp;
01435                 for( BU_PTBL_FOR( psp, (struct rt_piecestate **), &resp->re_pieces_pending ) )  {
01436                         if( (*psp)->htab.end > 0 )
01437                                 rt_htbl_reset( &(*psp)->htab );
01438                 }
01439                 bu_ptbl_reset( &resp->re_pieces_pending );
01440         }
01441 
01442         /* Terminate any logging */
01443         if(RT_G_DEBUG&(DEBUG_ALLRAYS|DEBUG_SHOOT|DEBUG_PARTITION|DEBUG_ALLHITS))  {
01444                 bu_log_indent_delta(-2);
01445                 bu_log("----------shootray cpu=%d  %d,%d lvl=%d (%s) %s ret=%d\n",
01446                         resp->re_cpu,
01447                         ap->a_x, ap->a_y,
01448                         ap->a_level,
01449                         ap->a_purpose != (char *)0 ? ap->a_purpose : "?",
01450                         status, ap->a_return);
01451         }
01452         return( ap->a_return );
01453 }
01454 
01455 /*
01456  *                      R T _ C E L L _ N _ O N _ R A Y
01457  *
01458  *  Return pointer to cell 'n' along a given ray.
01459  *  Used for debugging of how space partitioning interacts with shootray.
01460  *  Intended to mirror the operation of rt_shootray().
01461  *  The first cell is 0.
01462  */
01463 const union cutter *
01464 rt_cell_n_on_ray(register struct application *ap, int n)
01465 
01466                         /* First cell is #0 */
01467 {
01468         struct rt_shootray_status       ss;
01469         register const union cutter *cutp;
01470         struct resource         *resp;
01471         struct rt_i             *rtip;
01472         const int               debug_shoot = RT_G_DEBUG & DEBUG_SHOOT;
01473 
01474         RT_AP_CHECK(ap);
01475         if( ap->a_magic )  {
01476                 RT_CK_AP(ap);
01477         } else {
01478                 ap->a_magic = RT_AP_MAGIC;
01479         }
01480         if( ap->a_ray.magic )  {
01481                 RT_CK_RAY(&(ap->a_ray));
01482         } else {
01483                 ap->a_ray.magic = RT_RAY_MAGIC;
01484         }
01485         if( ap->a_resource == RESOURCE_NULL )  {
01486                 ap->a_resource = &rt_uniresource;
01487                 if(RT_G_DEBUG)bu_log("rt_cell_n_on_ray:  defaulting a_resource to &rt_uniresource\n");
01488                 if( rt_uniresource.re_magic == 0 )
01489                         rt_init_resource( &rt_uniresource, 0, ap->a_rt_i );
01490         }
01491         ss.ap = ap;
01492         rtip = ap->a_rt_i;
01493         RT_CK_RTI( rtip );
01494         resp = ap->a_resource;
01495         RT_CK_RESOURCE(resp);
01496         ss.resp = resp;
01497 
01498         if(RT_G_DEBUG&(DEBUG_ALLRAYS|DEBUG_SHOOT|DEBUG_PARTITION|DEBUG_ALLHITS)) {
01499                 bu_log_indent_delta(2);
01500                 bu_log("\n**********cell_n_on_ray cpu=%d  %d,%d lvl=%d (%s), n=%d\n",
01501                         resp->re_cpu,
01502                         ap->a_x, ap->a_y,
01503                         ap->a_level,
01504                         ap->a_purpose != (char *)0 ? ap->a_purpose : "?", n );
01505                 bu_log("Pnt (%g, %g, %g) a_onehit=%d\n",
01506                         V3ARGS(ap->a_ray.r_pt),
01507                         ap->a_onehit );
01508                 VPRINT("Dir", ap->a_ray.r_dir);
01509         }
01510 #ifndef NO_BADRAY_CHECKING
01511         if(RT_BADVEC(ap->a_ray.r_pt)||RT_BADVEC(ap->a_ray.r_dir))  {
01512                 bu_log("\n**********cell_n_on_ray cpu=%d  %d,%d lvl=%d (%s)\n",
01513                         resp->re_cpu,
01514                         ap->a_x, ap->a_y,
01515                         ap->a_level,
01516                         ap->a_purpose != (char *)0 ? ap->a_purpose : "?" );
01517                 VPRINT(" r_pt", ap->a_ray.r_pt);
01518                 VPRINT("r_dir", ap->a_ray.r_dir);
01519                 rt_bomb("rt_cell_n_on_ray() bad ray\n");
01520         }
01521 #endif
01522 
01523         if( rtip->needprep )
01524                 rt_prep_parallel(rtip, 1);      /* Stay on our CPU */
01525 
01526         if( BU_LIST_UNINITIALIZED( &resp->re_parthead ) )  {
01527                 /* XXX This shouldn't happen any more */
01528                 bu_log("rt_cell_n_on_ray() resp=x%x uninitialized, fixing it\n", resp);
01529                 /*
01530                  *  We've been handed a mostly un-initialized resource struct,
01531                  *  with only a magic number and a cpu number filled in.
01532                  *  Init it and add it to the table.
01533                  *  This is how application-provided resource structures
01534                  *  are remembered for later cleanup by the library.
01535                  */
01536                 rt_init_resource( resp, resp->re_cpu, rtip );
01537         }
01538         /* Ensure that this CPU's resource structure is registered */
01539         BU_ASSERT_PTR( BU_PTBL_GET(&rtip->rti_resources, resp->re_cpu), !=, NULL );
01540 
01541         /* Verify that direction vector has unit length */
01542         if(RT_G_DEBUG) {
01543                 FAST fastf_t f, diff;
01544                 /* Fancy version of BN_VEC_NON_UNIT_LEN() */
01545                 f = MAGSQ(ap->a_ray.r_dir);
01546                 if( NEAR_ZERO(f, 0.0001) )  {
01547                         rt_bomb("rt_cell_n_on_ray:  zero length dir vector\n");
01548                         return CUTTER_NULL;
01549                 }
01550                 diff = f - 1;
01551                 if( !NEAR_ZERO( diff, 0.0001 ) )  {
01552                         bu_log("rt_cell_n_on_ray: non-unit dir vect (x%d y%d lvl%d)\n",
01553                                 ap->a_x, ap->a_y, ap->a_level );
01554                         f = 1/f;
01555                         VSCALE( ap->a_ray.r_dir, ap->a_ray.r_dir, f );
01556                 }
01557         }
01558 
01559         /* Compute the inverse of the direction cosines */
01560         if( ap->a_ray.r_dir[X] < -SQRT_SMALL_FASTF )  {
01561                 ss.abs_inv_dir[X] = -(ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
01562                 ss.rstep[X] = -1;
01563         } else if( ap->a_ray.r_dir[X] > SQRT_SMALL_FASTF )  {
01564                 ss.abs_inv_dir[X] =  (ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
01565                 ss.rstep[X] = 1;
01566         } else {
01567                 ap->a_ray.r_dir[X] = 0.0;
01568                 ss.abs_inv_dir[X] = ss.inv_dir[X] = INFINITY;
01569                 ss.rstep[X] = 0;
01570         }
01571         if( ap->a_ray.r_dir[Y] < -SQRT_SMALL_FASTF )  {
01572                 ss.abs_inv_dir[Y] = -(ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
01573                 ss.rstep[Y] = -1;
01574         } else if( ap->a_ray.r_dir[Y] > SQRT_SMALL_FASTF )  {
01575                 ss.abs_inv_dir[Y] =  (ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
01576                 ss.rstep[Y] = 1;
01577         } else {
01578                 ap->a_ray.r_dir[Y] = 0.0;
01579                 ss.abs_inv_dir[Y] = ss.inv_dir[Y] = INFINITY;
01580                 ss.rstep[Y] = 0;
01581         }
01582         if( ap->a_ray.r_dir[Z] < -SQRT_SMALL_FASTF )  {
01583                 ss.abs_inv_dir[Z] = -(ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
01584                 ss.rstep[Z] = -1;
01585         } else if( ap->a_ray.r_dir[Z] > SQRT_SMALL_FASTF )  {
01586                 ss.abs_inv_dir[Z] =  (ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
01587                 ss.rstep[Z] = 1;
01588         } else {
01589                 ap->a_ray.r_dir[Z] = 0.0;
01590                 ss.abs_inv_dir[Z] = ss.inv_dir[Z] = INFINITY;
01591                 ss.rstep[Z] = 0;
01592         }
01593 
01594         /*
01595          *  If ray does not enter the model RPP, skip on.
01596          *  If ray ends exactly at the model RPP, trace it.
01597          */
01598         if( !rt_in_rpp( &ap->a_ray, ss.inv_dir, rtip->mdl_min, rtip->mdl_max )  ||
01599             ap->a_ray.r_max < 0.0 )  {
01600                 cutp = &ap->a_rt_i->rti_inf_box;
01601                 if( cutp->bn.bn_len > 0 )  {
01602                         if( n == 0 )  return cutp;
01603                 }
01604                 return CUTTER_NULL;
01605         }
01606 
01607         /*
01608          *  The interesting part of the ray starts at distance 0.
01609          *  If the ray enters the model at a negative distance,
01610          *  (ie, the ray starts within the model RPP),
01611          *  we only look at little bit behind (BACKING_DIST) to see if we are
01612          *  just coming out of something, but never further back than
01613          *  the intersection with the model RPP.
01614          *  If the ray enters the model at a positive distance,
01615          *  we always start there.
01616          *  It is vital that we never pick a start point outside the
01617          *  model RPP, or the space partitioning tree will pick the
01618          *  wrong box and the ray will miss it.
01619          *
01620          *  BACKING_DIST should probably be determined by floating point
01621          *  noise factor due to model RPP size -vs- number of bits of
01622          *  floating point mantissa significance, rather than a constant,
01623          *  but that is too hideous to think about here.
01624          *  Also note that applications that really depend on knowing
01625          *  what region they are leaving from should probably back their
01626          *  own start-point up, rather than depending on it here, but
01627          *  it isn't much trouble here.
01628          */
01629         ss.box_start = ss.model_start = ap->a_ray.r_min;
01630         ss.box_end = ss.model_end = ap->a_ray.r_max;
01631 
01632         if( ss.box_start < BACKING_DIST )
01633                 ss.box_start = BACKING_DIST; /* Only look a little bit behind */
01634 
01635         ss.lastcut = CUTTER_NULL;
01636         ss.old_status = (struct rt_shootray_status *)NULL;
01637         ss.curcut = &ap->a_rt_i->rti_CutHead;
01638         if( ss.curcut->cut_type == CUT_NUGRIDNODE ) {
01639                 ss.lastcell = CUTTER_NULL;
01640                 VSET( ss.curmin, ss.curcut->nugn.nu_axis[X][0].nu_spos,
01641                                  ss.curcut->nugn.nu_axis[Y][0].nu_spos,
01642                                  ss.curcut->nugn.nu_axis[Z][0].nu_spos );
01643                 VSET( ss.curmax, ss.curcut->nugn.nu_axis[X][ss.curcut->nugn.nu_cells_per_axis[X]-1].nu_epos,
01644                                  ss.curcut->nugn.nu_axis[Y][ss.curcut->nugn.nu_cells_per_axis[Y]-1].nu_epos,
01645                                  ss.curcut->nugn.nu_axis[Z][ss.curcut->nugn.nu_cells_per_axis[Z]-1].nu_epos );
01646         } else if( ss.curcut->cut_type == CUT_CUTNODE ||
01647                    ss.curcut->cut_type == CUT_BOXNODE ) {
01648                 ss.lastcell = ss.curcut;
01649                 VMOVE( ss.curmin, rtip->mdl_min );
01650                 VMOVE( ss.curmax, rtip->mdl_max );
01651         }
01652 
01653         ss.newray = ap->a_ray;          /* struct copy */
01654         ss.odist_corr = ss.obox_start = ss.obox_end = -99;
01655         ss.dist_corr = 0.0;
01656 
01657         /*
01658          *  While the ray remains inside model space,
01659          *  push from box to box until ray emerges from
01660          *  model space again (or first hit is found, if user is impatient).
01661          *  It is vitally important to always stay within the model RPP, or
01662          *  the space partitoning tree will pick wrong boxes & miss them.
01663          */
01664         while( (cutp = rt_advance_to_next_cell( &ss )) != CUTTER_NULL )  {
01665                 if(debug_shoot) {
01666                         rt_pr_cut( cutp, 0 );
01667                 }
01668                 if( --n <= 0 )  return cutp;
01669 
01670                 /* Push ray onwards to next box */
01671                 ss.box_start = ss.box_end;
01672         }
01673         return CUTTER_NULL;
01674 }
01675 
01676 /*
01677  *                      R T _ I N _ R P P
01678  *
01679  *  Compute the intersections of a ray with a rectangular parallelpiped (RPP)
01680  *  that has faces parallel to the coordinate planes
01681  *
01682  *  The algorithm here was developed by Gary Kuehl for GIFT.
01683  *  A good description of the approach used can be found in
01684  *  "??" by XYZZY and Barsky,
01685  *  ACM Transactions on Graphics, Vol 3 No 1, January 1984.
01686  *
01687  * Note -
01688  *  The computation of entry and exit distance is mandatory, as the final
01689  *  test catches the majority of misses.
01690  *
01691  * Note -
01692  *  A hit is returned if the intersect is behind the start point.
01693  *
01694  *  Returns -
01695  *       0  if ray does not hit RPP,
01696  *      !0  if ray hits RPP.
01697  *
01698  *  Implicit return -
01699  *      rp->r_min = dist from start of ray to point at which ray ENTERS solid
01700  *      rp->r_max = dist from start of ray to point at which ray LEAVES solid
01701  */
01702 int
01703 rt_in_rpp(struct xray           *rp,
01704           register const fastf_t *invdir,       /* inverses of rp->r_dir[] */
01705           register const fastf_t *min,
01706           register const fastf_t *max)
01707 {
01708         register const fastf_t  *pt = &rp->r_pt[0];
01709         register fastf_t        sv;
01710 #define st sv                   /* reuse the register */
01711         register fastf_t        rmin = -INFINITY;
01712         register fastf_t        rmax =  INFINITY;
01713 
01714         /* Start with infinite ray, and trim it down */
01715 
01716         /* X axis */
01717         if( *invdir < 0.0 )  {
01718                 /* Heading towards smaller numbers */
01719                 /* if( *min > *pt )  miss */
01720                 if(rmax > (sv = (*min - *pt) * *invdir) )
01721                         rmax = sv;
01722                 if( rmin < (st = (*max - *pt) * *invdir) )
01723                         rmin = st;
01724         }  else if( *invdir > 0.0 )  {
01725                 /* Heading towards larger numbers */
01726                 /* if( *max < *pt )  miss */
01727                 if(rmax > (st = (*max - *pt) * *invdir) )
01728                         rmax = st;
01729                 if( rmin < ((sv = (*min - *pt) * *invdir)) )
01730                         rmin = sv;
01731         }  else  {
01732                 /*
01733                  *  Direction cosines along this axis is NEAR 0,
01734                  *  which implies that the ray is perpendicular to the axis,
01735                  *  so merely check position against the boundaries.
01736                  */
01737                 if( (*min > *pt) || (*max < *pt) )
01738                         return(0);      /* MISS */
01739         }
01740 
01741         /* Y axis */
01742         pt++; invdir++; max++; min++;
01743         if( *invdir < 0.0 )  {
01744                 if(rmax > (sv = (*min - *pt) * *invdir) )
01745                         rmax = sv;
01746                 if( rmin < (st = (*max - *pt) * *invdir) )
01747                         rmin = st;
01748         }  else if( *invdir > 0.0 )  {
01749                 if(rmax > (st = (*max - *pt) * *invdir) )
01750                         rmax = st;
01751                 if( rmin < ((sv = (*min - *pt) * *invdir)) )
01752                         rmin = sv;
01753         }  else  {
01754                 if( (*min > *pt) || (*max < *pt) )
01755                         return(0);      /* MISS */
01756         }
01757 
01758         /* Z axis */
01759         pt++; invdir++; max++; min++;
01760         if( *invdir < 0.0 )  {
01761                 if(rmax > (sv = (*min - *pt) * *invdir) )
01762                         rmax = sv;
01763                 if( rmin < (st = (*max - *pt) * *invdir) )
01764                         rmin = st;
01765         }  else if( *invdir > 0.0 )  {
01766                 if(rmax > (st = (*max - *pt) * *invdir) )
01767                         rmax = st;
01768                 if( rmin < ((sv = (*min - *pt) * *invdir)) )
01769                         rmin = sv;
01770         }  else  {
01771                 if( (*min > *pt) || (*max < *pt) )
01772                         return(0);      /* MISS */
01773         }
01774 
01775         /* If equal, RPP is actually a plane */
01776         if( rmin > rmax )
01777                 return(0);      /* MISS */
01778 
01779         /* HIT.  Only now do rp->r_min and rp->r_max have to be written */
01780         rp->r_min = rmin;
01781         rp->r_max = rmax;
01782         return(1);              /* HIT */
01783 }
01784 
01785 /* For debugging */
01786 int
01787 rt_DB_rpp(register struct xray *rp, register const fastf_t *invdir, register const fastf_t *min, register const fastf_t *max)
01788 
01789                                 /* inverses of rp->r_dir[] */
01790 
01791 
01792 {
01793         register const fastf_t *pt = &rp->r_pt[0];
01794         FAST fastf_t sv;
01795 
01796         /* Start with infinite ray, and trim it down */
01797         rp->r_min = -INFINITY;
01798         rp->r_max = INFINITY;
01799 
01800 VPRINT("r_pt ", pt);
01801 VPRINT("r_dir", rp->r_dir);
01802 VPRINT("min  ", min);
01803 VPRINT("max  ", max);
01804         /* X axis */
01805         if( rp->r_dir[X] < 0.0 )  {
01806                 /* Heading towards smaller numbers */
01807                 /* if( *min > *pt )  miss */
01808                 sv = (*min - *pt) * *invdir;
01809 bu_log("sv=%g, r_max=%g\n", sv, rp->r_max);
01810                 if(rp->r_max > sv)
01811                         rp->r_max = sv;
01812                 st = (*max - *pt) * *invdir;
01813 bu_log("st=%g, r_min=%g\n", st, rp->r_min);
01814                 if( rp->r_min < st )
01815                         rp->r_min = st;
01816 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01817         }  else if( rp->r_dir[X] > 0.0 )  {
01818                 /* Heading towards larger numbers */
01819                 /* if( *max < *pt )  miss */
01820                 st = (*max - *pt) * *invdir;
01821 bu_log("st=%g, r_max=%g\n", st, rp->r_max);
01822                 if(rp->r_max > st)
01823                         rp->r_max = st;
01824                 sv = (*min - *pt) * *invdir;
01825 bu_log("sv=%g, r_min=%g\n", sv, rp->r_min);
01826                 if( rp->r_min < sv )
01827                         rp->r_min = sv;
01828 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01829         }  else  {
01830                 /*
01831                  *  Direction cosines along this axis is NEAR 0,
01832                  *  which implies that the ray is perpendicular to the axis,
01833                  *  so merely check position against the boundaries.
01834                  */
01835                 if( (*min > *pt) || (*max < *pt) )
01836                         goto miss;
01837         }
01838 
01839         /* Y axis */
01840         pt++; invdir++; max++; min++;
01841         if( rp->r_dir[Y] < 0.0 )  {
01842                 /* Heading towards smaller numbers */
01843                 /* if( *min > *pt )  miss */
01844                 sv = (*min - *pt) * *invdir;
01845 bu_log("sv=%g, r_max=%g\n", sv, rp->r_max);
01846                 if(rp->r_max > sv)
01847                         rp->r_max = sv;
01848                 st = (*max - *pt) * *invdir;
01849 bu_log("st=%g, r_min=%g\n", st, rp->r_min);
01850                 if( rp->r_min < st )
01851                         rp->r_min = st;
01852 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01853         }  else if( rp->r_dir[Y] > 0.0 )  {
01854                 /* Heading towards larger numbers */
01855                 /* if( *max < *pt )  miss */
01856                 st = (*max - *pt) * *invdir;
01857 bu_log("st=%g, r_max=%g\n", st, rp->r_max);
01858                 if(rp->r_max > st)
01859                         rp->r_max = st;
01860                 sv = (*min - *pt) * *invdir;
01861 bu_log("sv=%g, r_min=%g\n", sv, rp->r_min);
01862                 if( rp->r_min < sv )
01863                         rp->r_min = sv;
01864 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01865         }  else  {
01866                 if( (*min > *pt) || (*max < *pt) )
01867                         goto miss;
01868         }
01869 
01870         /* Z axis */
01871         pt++; invdir++; max++; min++;
01872         if( rp->r_dir[Z] < 0.0 )  {
01873                 /* Heading towards smaller numbers */
01874                 /* if( *min > *pt )  miss */
01875                 sv = (*min - *pt) * *invdir;
01876 bu_log("sv=%g, r_max=%g\n", sv, rp->r_max);
01877                 if(rp->r_max > sv)
01878                         rp->r_max = sv;
01879                 st = (*max - *pt) * *invdir;
01880 bu_log("st=%g, r_min=%g\n", st, rp->r_min);
01881                 if( rp->r_min < st )
01882                         rp->r_min = st;
01883 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01884         }  else if( rp->r_dir[Z] > 0.0 )  {
01885                 /* Heading towards larger numbers */
01886                 /* if( *max < *pt )  miss */
01887                 st = (*max - *pt) * *invdir;
01888 bu_log("st=%g, r_max=%g\n", st, rp->r_max);
01889                 if(rp->r_max > st)
01890                         rp->r_max = st;
01891                 sv = (*min - *pt) * *invdir;
01892 bu_log("sv=%g, r_min=%g\n", sv, rp->r_min);
01893                 if( rp->r_min < sv )
01894                         rp->r_min = sv;
01895 bu_log("r_min=%g, r_max=%g\n", rp->r_min, rp->r_max);
01896         }  else  {
01897                 if( (*min > *pt) || (*max < *pt) )
01898                         goto miss;
01899         }
01900 
01901         /* If equal, RPP is actually a plane */
01902         if( rp->r_min > rp->r_max )
01903                 goto miss;
01904         bu_log("HIT:  %g..%g\n", rp->r_min, rp->r_max );
01905         return(1);              /* HIT */
01906 miss:
01907         bu_log("MISS\n");
01908         return(0);              /* MISS */
01909 }
01910 
01911 #undef st
01912 
01913 
01914 #define SEG_MISS(SEG)           (SEG).seg_stp=(struct soltab *) 0;
01915 
01916 /* Stub function which will "similate" a call to a vector shot routine */
01917 void
01918 rt_vstub(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
01919                                /* An array of solid pointers */
01920                                /* An array of ray pointers */
01921                                /* array of segs (results returned) */
01922                                /* Number of ray/object pairs */
01923                              /* pointer to an application */
01924 {
01925         register int    i;
01926         register struct seg *tmp_seg;
01927         struct seg      seghead;
01928 
01929         BU_LIST_INIT( &(seghead.l) );
01930 
01931         /* go through each ray/solid pair and call a scalar function */
01932         for (i = 0; i < n; i++) {
01933                 if (stp[i] != 0){ /* skip call if solid table pointer is NULL */
01934                         /* do scalar call, place results in segp array */
01935                         if( rt_functab[stp[i]->st_id].ft_shot(stp[i], rp[i], ap, &seghead) <= 0 )  {
01936                                 SEG_MISS(segp[i]);
01937                         } else {
01938                                 tmp_seg = BU_LIST_FIRST(seg, &(seghead.l) );
01939                                 BU_LIST_DEQUEUE( &(tmp_seg->l) );
01940                                 segp[i] = *tmp_seg; /* structure copy */
01941                                 RT_FREE_SEG(tmp_seg, ap->a_resource);
01942                         }
01943                 }
01944         }
01945 }
01946 
01947 /*
01948  *                      R T _ P R _ L I B R A R Y _ V E R S I O N
01949  *
01950  *  In case anyone actually cares, print out library's compilation version.
01951  */
01952 void
01953 rt_pr_library_version(void)
01954 {
01955         bu_log("%s", rt_version);
01956 }
01957 
01958 void
01959 rt_zero_res_stats( struct resource *resp )
01960 {
01961         RT_CK_RESOURCE( resp );
01962 
01963         resp->re_nshootray = 0;
01964         resp->re_nmiss_model = 0;
01965 
01966         resp->re_shots = 0;
01967         resp->re_shot_hit = 0;
01968         resp->re_shot_miss = 0;
01969 
01970         resp->re_prune_solrpp = 0;
01971 
01972         resp->re_ndup = 0;
01973         resp->re_nempty_cells = 0;
01974 
01975         resp->re_piece_shots = 0;
01976         resp->re_piece_shot_hit = 0;
01977         resp->re_piece_shot_miss = 0;
01978         resp->re_piece_ndup = 0;
01979 }
01980 
01981 /*
01982  *                      R T _ A D D _ R E S _ S T A T S
01983  *
01984  *  To be called only in non-parallel mode, to tally up the statistics
01985  *  from the resource structure(s) into the rt instance structure.
01986  *
01987  *  Non-parallel programs should call
01988  *      rt_add_res_stats( rtip, RESOURCE_NULL );
01989  *  to have the default resource results tallied in.
01990  */
01991 void
01992 rt_add_res_stats(register struct rt_i *rtip, register struct resource *resp)
01993 {
01994         RT_CK_RTI( rtip );
01995 
01996         if( resp == RESOURCE_NULL )  resp = &rt_uniresource;
01997         RT_CK_RESOURCE( resp );
01998 
01999         rtip->rti_nrays += resp->re_nshootray;
02000         rtip->nmiss_model += resp->re_nmiss_model;
02001 
02002         rtip->nshots += resp->re_shots + resp->re_piece_shots;
02003         rtip->nhits += resp->re_shot_hit + resp->re_piece_shot_hit;
02004         rtip->nmiss += resp->re_shot_miss + resp->re_piece_shot_miss;
02005 
02006         rtip->nmiss_solid += resp->re_prune_solrpp;
02007 
02008         rtip->ndup += resp->re_ndup + resp->re_piece_ndup;
02009         rtip->nempty_cells += resp->re_nempty_cells;
02010 
02011         /* Zero out resource totals, so repeated calls are not harmful */
02012         rt_zero_res_stats( resp );
02013 }
02014 
02015 /*
02016  *  Routines for plotting the progress of one ray through the model.  -Mike
02017  */
02018 void
02019 rt_3move_raydist(FILE *fp, struct xray *rayp, double dist)
02020 {
02021         point_t p;
02022 
02023         VJOIN1( p, rayp->r_pt, dist, rayp->r_dir );
02024         pdv_3move( fp, p );
02025 }
02026 
02027 void
02028 rt_3cont_raydist(FILE *fp, struct xray *rayp, double dist)
02029 {
02030         point_t p;
02031 
02032         VJOIN1( p, rayp->r_pt, dist, rayp->r_dir );
02033         pdv_3cont( fp, p );
02034 }
02035 
02036 /*
02037                 rt_plot_cell( cutp, &ss, &(waiting_segs.l), rtip);
02038  */
02039 void
02040 rt_plot_cell(const union cutter *cutp, const struct rt_shootray_status *ssp, struct bu_list *waiting_segs_hd, struct rt_i *rtip)
02041 {
02042         char            buf[128];
02043         static int      fnum = 0;
02044         FILE            *fp;
02045         struct soltab   **stpp;
02046         struct application      *ap;
02047 
02048         RT_CK_RTI(rtip);
02049         RT_AP_CHECK(ssp->ap);
02050         RT_CK_RTI(ssp->ap->a_rt_i);
02051         ap = ssp->ap;
02052 
02053         sprintf( buf, "cell%d.pl", fnum++ );
02054         if( (fp = fopen( buf, "w" )) == NULL )  {
02055                 perror(buf);
02056         }
02057 
02058         pl_color( fp, 0, 100, 0 );              /* green box for model RPP */
02059 
02060         /* Plot the model RPP, to provide some context */
02061         pdv_3space( fp, rtip->rti_pmin, rtip->rti_pmax );
02062         pdv_3box( fp, rtip->rti_pmin, rtip->rti_pmax );
02063 
02064         /* Plot the outline of this one cell */
02065         pl_color( fp, 80, 80, 250 );
02066         switch( cutp->cut_type )  {
02067         case CUT_BOXNODE:
02068                 pdv_3box( fp, cutp->bn.bn_min, cutp->bn.bn_max );
02069                 break;
02070         default:
02071                 bu_log("cut_type = %d\n", cutp->cut_type );
02072                 bu_bomb("Unknown cut_type\n");
02073         }
02074 
02075         if( cutp->bn.bn_len > 0 ) {
02076                 /* Plot every solid listed in this cell */
02077                 stpp = &(cutp->bn.bn_list[cutp->bn.bn_len-1]);
02078                 for( ; stpp >= cutp->bn.bn_list; stpp-- )  {
02079                         register struct soltab *stp = *stpp;
02080 
02081                         rt_plot_solid( fp, rtip, stp, ap->a_resource );
02082                 }
02083         }
02084 
02085         /* Plot interval of ray in box, in green */
02086         pl_color( fp, 100, 255, 200 );
02087         rt_3move_raydist( fp, &ap->a_ray, ssp->box_start );
02088         rt_3cont_raydist( fp, &ap->a_ray, ssp->box_end );
02089 
02090         if( bu_list_len( waiting_segs_hd ) <= 0 )  {
02091                 /* No segments, just plot the whole ray */
02092                 pl_color( fp, 255, 255, 0 );    /* yellow -- no segs */
02093                 rt_3move_raydist( fp, &ap->a_ray, ssp->model_start );
02094                 rt_3cont_raydist( fp, &ap->a_ray, ssp->box_start );
02095                 rt_3move_raydist( fp, &ap->a_ray, ssp->box_end );
02096                 rt_3cont_raydist( fp, &ap->a_ray, ssp->model_end );
02097         } else {
02098                 /* Plot the segments awaiting boolweave. */
02099                 struct seg      *segp;
02100 
02101                 for( BU_LIST_FOR( segp, seg, waiting_segs_hd ) )  {
02102                         RT_CK_SEG(segp);
02103                         pl_color( fp, 255, 0, 0 );      /* red */
02104                         rt_3move_raydist( fp, &ap->a_ray, segp->seg_in.hit_dist );
02105                         rt_3cont_raydist( fp, &ap->a_ray, segp->seg_out.hit_dist );
02106                 }
02107         }
02108 
02109         fclose(fp);
02110         bu_log("wrote %s\n", buf);
02111 }
02112 
02113 /*@}*/
02114 /*
02115  * Local Variables:
02116  * mode: C
02117  * tab-width: 8
02118  * c-basic-offset: 4
02119  * indent-tabs-mode: t
02120  * End:
02121  * ex: shiftwidth=4 tabstop=8
02122  */

Generated on Mon Sep 18 01:24:56 2006 for BRL-CAD by  doxygen 1.4.6