BRL-CAD
shoot.c
Go to the documentation of this file.
1 /* S H O O T . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2000-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 
21 
22 #include "common.h"
23 
24 #include <math.h>
25 #include <string.h>
26 #include "bio.h"
27 
28 #include "vmath.h"
29 
30 #include "raytrace.h"
31 #include "plot3.h"
32 
33 
34 #define V3PT_DEPARTING_RPP(_step, _lo, _hi, _pt) \
35  PT_DEPARTING_RPP(_step, _lo, _hi, (_pt)[X], (_pt)[Y], (_pt)[Z])
36 #define PT_DEPARTING_RPP(_step, _lo, _hi, _px, _py, _pz) \
37  (((_step)[X] <= 0 && (_px) < (_lo)[X]) || \
38  ((_step)[X] >= 0 && (_px) > (_hi)[X]) || \
39  ((_step)[Y] <= 0 && (_py) < (_lo)[Y]) || \
40  ((_step)[Y] >= 0 && (_py) > (_hi)[Y]) || \
41  ((_step)[Z] <= 0 && (_pz) < (_lo)[Z]) || \
42  ((_step)[Z] >= 0 && (_pz) > (_hi)[Z]))
43 
44 
45 void
46 rt_res_pieces_init(struct resource *resp, struct rt_i *rtip)
47 {
48  struct rt_piecestate *psptab;
49  struct rt_piecestate *psp;
50  struct soltab *stp;
51 
52  RT_CK_RESOURCE(resp);
53  RT_CK_RTI(rtip);
54 
55  psptab = (struct rt_piecestate *)bu_calloc(rtip->rti_nsolids_with_pieces, sizeof(struct rt_piecestate), "re_pieces[]");
56  resp->re_pieces = psptab;
57 
58  RT_VISIT_ALL_SOLTABS_START(stp, rtip) {
59  RT_CK_SOLTAB(stp);
60  if (stp->st_npieces <= 1) continue;
61  psp = &psptab[stp->st_piecestate_num];
63  psp->stp = stp;
64  psp->shot = bu_bitv_new(stp->st_npieces);
65  rt_htbl_init(&psp->htab, 8, "psp->htab");
66  psp->cutp = CUTTER_NULL;
68  }
69 
70 
71 void
72 rt_res_pieces_clean(struct resource *resp, struct rt_i *rtip)
73 {
74  struct rt_piecestate *psp;
75  int i;
76 
77  RT_CK_RESOURCE(resp);
78  if (rtip) {
79  RT_CK_RTI(rtip);
80  }
81  if (BU_PTBL_TEST(&resp->re_pieces_pending)) {
83  }
84 
85  if (!resp->re_pieces) {
86  /* no pieces allocated, nothing to do */
87  if (rtip) {
88  rtip->rti_nsolids_with_pieces = 0;
89  }
90  return;
91  }
92 
93  /* pieces allocated, need to clean up */
94  if (rtip) {
95  for (i = rtip->rti_nsolids_with_pieces-1; i >= 0; i--) {
96  psp = &resp->re_pieces[i];
97 
98  /*
99  * Skip uninitialized structures.
100  *
101  * Doing this until we figure out why all "struct
102  * rt_piecestate" array members are NOT getting
103  * initialized in rt_res_pieces_init() above. Initial
104  * glance looks like tree.c/rt_gettrees_muves() can be
105  * called in such a way as to assign
106  * stp->st_piecestate_num multiple times, each time with a
107  * different value. This value is an index into the
108  * resp->re_pieces array. When this happens, it causes all
109  * but the last referenced "struct rt_piecestate" member
110  * to be initialized.
111  */
112  if (psp->magic == 0)
113  continue;
114 
115  RT_CK_PIECESTATE(psp);
116  rt_htbl_free(&psp->htab);
117  bu_bitv_free(psp->shot);
118  psp->shot = NULL; /* sanity */
119  psp->magic = 0;
120  }
121 
122  rtip->rti_nsolids_with_pieces = 0;
123  }
124  bu_free((char *)resp->re_pieces, "re_pieces[]");
125  resp->re_pieces = NULL;
126 }
127 
128 
129 /**
130  * Along the given axis, find which NUgrid cell this value lies in.
131  * Use method of binary subdivision.
132  */
133 int
134 rt_find_nugrid(const struct nugridnode *nugnp, int axis, fastf_t val)
135 {
136  int min;
137  int max;
138  int lim = nugnp->nu_cells_per_axis[axis]-1;
139  int cur;
140 
141  if (val < nugnp->nu_axis[axis][0].nu_spos ||
142  val > nugnp->nu_axis[axis][lim].nu_epos)
143  return -1;
144 
145  min = 0;
146  max = lim;
147 again:
148  cur = (min + max) / 2;
149 
150  if (cur <= 0) return 0;
151  if (cur >= lim) return lim;
152 
153  if (val < nugnp->nu_axis[axis][cur].nu_spos) {
154  max = cur;
155  goto again;
156  }
157  if (val > nugnp->nu_axis[axis][cur+1].nu_epos) {
158  min = cur+1;
159  goto again;
160  }
161 
162  if (val < nugnp->nu_axis[axis][cur].nu_epos)
163  return cur;
164 
165  return cur+1;
166 }
167 
168 
169 const union cutter *
171 {
172  register const union cutter *cutp, *curcut = ssp->curcut;
173  register const struct application *ap = ssp->ap;
174  register fastf_t t0, px, py, pz;
175  int push_flag = 0;
176  double fraction;
177  int exponent;
178 
179  ssp->box_num++;
180 
181  if (curcut == &ssp->ap->a_rt_i->rti_inf_box) {
182  /* Last pass did the infinite solids, there is nothing more */
183  ssp->curcut = CUTTER_NULL;
184  return CUTTER_NULL;
185  }
186 
187  for (;;) {
188  /* Set cutp to CUTTER_NULL. If it fails to become set in the
189  * following switch statement, we know that we have exited the
190  * subnode. If this subnode is the highest-level node, then
191  * we are done advancing the ray through the model.
192  */
193  cutp = CUTTER_NULL;
194  push_flag = 0;
195 
196  /* The point corresponding to the box_start distance may not
197  * be in the "right" place, due to the effects of floating
198  * point fuzz:
199  *
200  * 1) The point might lie just outside the model RPP,
201  * resulting in the point not falling within the RPP of the
202  * indicated cell, or
203  *
204  * 2) The poing might lie just a little bit on the wrong side
205  * of the cell wall, resulting in the ray getting "stuck", and
206  * needing rescuing all the time by the error recovery code
207  * below.
208  *
209  * Therefore, "nudge" the point just slightly into the next
210  * cell by adding our distance tolerance.
211  *
212  * XXX At present, a cell is never less than 1mm wide.
213  *
214  * XXX The "nudge" value was based on an absolute value defined
215  * by OFFSET_DIST but has been changed to use distance tolerance
216  * specified in mm and can now be overridden by a user.
217  */
218 
219  t0 = ssp->box_start;
220  /* NB: can't compute px, py, pz here since t0 may advance in
221  * the following statement!
222  */
223 
224  top:
225  switch (curcut->cut_type) {
226  case CUT_NUGRIDNODE: {
227  /*********************************************************
228  * NOTE: This portion implements Gigante's non-uniform 3-D
229  * space grid/mesh discretization.
230  *********************************************************/
231  register int out_axis;
232  register const struct nu_axis **nu_axis =
233  (const struct nu_axis **)&curcut->nugn.nu_axis[0];
234  register const int *nu_stepsize =
235  &curcut->nugn.nu_stepsize[0];
236  register const int *nu_cells_per_axis =
237  &curcut->nugn.nu_cells_per_axis[0];
238  register const union cutter *nu_grid =
239  curcut->nugn.nu_grid;
240 
241  if (ssp->lastcell == CUTTER_NULL || ssp->lastcell->cut_type != CUT_NUGRIDNODE) {
242  /* We have just started into this NUgrid. We must
243  * find our location and set up the NUgrid traversal
244  * state variables.
245  */
246  register int x, y, z;
247 
248  px = ap->a_ray.r_pt[X] + t0*ap->a_ray.r_dir[X];
249  py = ap->a_ray.r_pt[Y] + t0*ap->a_ray.r_dir[Y];
250  pz = ap->a_ray.r_pt[Z] + t0*ap->a_ray.r_dir[Z];
251 
252  /* Must find cell that contains newray.r_pt. We do
253  * this by binary subdivision. If any are out of
254  * bounds, we have left the NUgrid and will pop a
255  * level off the stack in the outer loop (if
256  * applicable).
257  */
258  x = rt_find_nugrid(&curcut->nugn, X, px);
259  if (x<0) break;
260  y = rt_find_nugrid(&curcut->nugn, Y, py);
261  if (y<0) break;
262  z = rt_find_nugrid(&curcut->nugn, Z, pz);
263  if (z<0) break;
264 
265  cutp = &nu_grid[z*nu_stepsize[Z] +
266  y*nu_stepsize[Y] +
267  x*nu_stepsize[X]];
268 
269  ssp->igrid[X] = x;
270  ssp->igrid[Y] = y;
271  ssp->igrid[Z] = z;
272 
273  NUGRID_T_SETUP(X, px, x);
274  NUGRID_T_SETUP(Y, py, y);
275  NUGRID_T_SETUP(Z, pz, z);
276  } else {
277  /* Advance from previous cell to next cell. Take next
278  * step, finding ray entry distance.
279  */
280  cutp = ssp->lastcell;
281  out_axis = ssp->out_axis;
282 
283  /* We may be simply advancing to the next box in the
284  * *same* NUgrid cell (if, for instance, the NUgrid
285  * cell is a cutnode with boxnode leaves). So if t0
286  * hasn't advanced past the end of the box, advance a
287  * tiny bit (less than rt_ct_optim makes boxnodes) and
288  * be handled by tree-traversing code below.
289  */
290 
291  if (cutp->cut_type == CUT_CUTNODE &&
292  t0 + ssp->ap->a_rt_i->rti_tol.dist < ssp->tv[out_axis]) {
293  t0 += ssp->ap->a_rt_i->rti_tol.dist;
294  break;
295  }
296 
297  /* Advance to the next cell as appropriate, bailing
298  * out with cutp=CUTTER_NULL if we run past the end of
299  * the NUgrid array.
300  */
301 
302  nugrid_again:
303  t0 = ssp->tv[out_axis];
304  if (ssp->rstep[out_axis] > 0) {
305  if (++(ssp->igrid[out_axis]) >=
306  nu_cells_per_axis[out_axis]) {
307  cutp = CUTTER_NULL;
308  break;
309  }
310 
311  if (cutp != CUTTER_NULL && cutp != ssp->curcut) {
312  cutp += nu_stepsize[out_axis];
313  }
314  } else {
315  if (--(ssp->igrid[out_axis]) < 0) {
316  cutp = CUTTER_NULL;
317  break;
318  }
319  if (cutp != CUTTER_NULL && cutp != ssp->curcut) {
320  cutp -= nu_stepsize[out_axis];
321  }
322  }
323 
324  /* Update t advancement value for this axis */
325  NUGRID_T_ADV(out_axis, ssp->igrid[out_axis]);
326  }
327 
328  /* find minimum exit t value */
329  if (ssp->tv[X] < ssp->tv[Y]) {
330  if (ssp->tv[Z] < ssp->tv[X]) {
331  out_axis = Z;
332  } else {
333  out_axis = X;
334  }
335  } else {
336  if (ssp->tv[Z] < ssp->tv[Y]) {
337  out_axis = Z;
338  } else {
339  out_axis = Y;
340  }
341  }
342 
343  /* Zip through empty cells. */
344  if (cutp->cut_type == CUT_BOXNODE &&
345  cutp->bn.bn_len <= 0 &&
346  cutp->bn.bn_piecelen <= 0) {
347  ++ssp->resp->re_nempty_cells;
348  goto nugrid_again;
349  }
350 
351  ssp->out_axis = out_axis;
352  ssp->lastcell = cutp;
353 
355  bu_log("t0=%g found in cell (%d, %d, %d), out_axis=%c at %g; cell is %s\n",
356  t0,
357  ssp->igrid[X], ssp->igrid[Y],
358  ssp->igrid[Z],
359  "XYZ*"[ssp->out_axis],
360  ssp->tv[out_axis],
361  cutp->cut_type==CUT_CUTNODE?"cut":
362  cutp->cut_type==CUT_BOXNODE?"box":
363  cutp->cut_type==CUT_NUGRIDNODE?"nu":
364  "?");
365  break;
366  } /* CUT_NUGRIDNODE */
367 
368  case CUT_CUTNODE:
369  /* fall through */
370 
371  case CUT_BOXNODE:
372  /*********************************************************
373  * NOTE: This portion implements Muuss' non-uniform binary
374  * space partitioning tree.
375  *********************************************************/
376  t0 += ssp->ap->a_rt_i->rti_tol.dist;
377  cutp = curcut;
378  break;
379  default:
380  bu_bomb("rt_advance_to_next_cell: unknown high-level cutnode");
381  }
382 
383  if (cutp==CUTTER_NULL) {
384  pop_space_stack:
385  /*
386  * Pop the stack of nested space partitioning methods.
387  * Move up out of the current node, or return if there is
388  * nothing left to do.
389  */
390  {
391  register struct rt_shootray_status *old = ssp->old_status;
392 
393  if (old == NULL) goto escaped_from_model;
394  *ssp = *old; /* struct copy */
395  bu_free(old, "old rt_shootray_status");
396  curcut = ssp->curcut;
397  cutp = CUTTER_NULL;
398  continue;
399  }
400  }
401 
402  /* Compute position and bail if we're outside of the current
403  * level.
404  */
405  px = ap->a_ray.r_pt[X] + t0*ap->a_ray.r_dir[X];
406  py = ap->a_ray.r_pt[Y] + t0*ap->a_ray.r_dir[Y];
407  pz = ap->a_ray.r_pt[Z] + t0*ap->a_ray.r_dir[Z];
408 
409  /* Optimization: when it's a boxnode in a nugrid, just
410  * return.
411  */
412  if (cutp->cut_type == CUT_BOXNODE &&
413  curcut->cut_type == CUT_NUGRIDNODE) {
414  ssp->newray.r_pt[X] = px;
415  ssp->newray.r_pt[Y] = py;
416  ssp->newray.r_pt[Z] = pz;
417  ssp->newray.r_min = 0;
418  ssp->newray.r_max = ssp->tv[ssp->out_axis] - t0;
419  goto done_return_cutp;
420  }
421 
422  /* Given direction of travel, see if point is outside bound.
423  * This will be the model RPP for NUBSP.
424  */
425  if (PT_DEPARTING_RPP(ssp->rstep, ssp->curmin, ssp->curmax, px, py, pz))
426  goto pop_space_stack;
427 
429  bu_log(
430  "rt_advance_to_next_cell() dist_corr=%g, pt=(%g, %g, %g)\n",
431  t0 /*ssp->dist_corr*/, px, py, pz);
432  }
433 
434  while (cutp->cut_type == CUT_CUTNODE) {
435  switch (cutp->cn.cn_axis) {
436  case X:
437  if (!(px < cutp->cn.cn_point)) {
438  cutp=cutp->cn.cn_r;
439  } else {
440  cutp=cutp->cn.cn_l;
441  }
442  break;
443  case Y:
444  if (!(py < cutp->cn.cn_point)) {
445  cutp=cutp->cn.cn_r;
446  } else {
447  cutp=cutp->cn.cn_l;
448  }
449  break;
450  case Z:
451  if (!(pz < cutp->cn.cn_point)) {
452  cutp=cutp->cn.cn_r;
453  } else {
454  cutp=cutp->cn.cn_l;
455  }
456  break;
457  }
458  }
459 
460  if (cutp == CUTTER_NULL)
461  bu_bomb("rt_advance_to_next_cell: leaf is NULL!");
462 
463  switch (cutp->cut_type) {
464  case CUT_BOXNODE:
465  if (RT_G_DEBUG&DEBUG_ADVANCE &&
466  PT_DEPARTING_RPP(ssp->rstep, ssp->curmin, ssp->curmax, px, py, pz)
467  ) {
468  /* This cell is old news. */
469  bu_log("rt_advance_to_next_cell(): point not in cell, advancing\n pt (%.20e, %.20e, %.20e)\n",
470  px, py, pz);
471  bu_log(" min (%.20e, %.20e, %.20e)\n",
472  V3ARGS(cutp->bn.bn_min));
473  bu_log(" max (%.20e, %.20e, %.20e)\n",
474  V3ARGS(cutp->bn.bn_max));
475  bu_log("pt=(%g, %g, %g)\n", px, py, pz);
476  rt_pr_cut(cutp, 0);
477 
478  /*
479  * Move newray point further into new box.
480  * Try again.
481  */
482  t0 += ssp->ap->a_rt_i->rti_tol.dist;
483  goto top;
484  }
485  /* Don't get stuck within the same box for long */
486  if (cutp==ssp->lastcut) {
487  fastf_t delta;
488  push_to_next_box: ;
489  if (RT_G_DEBUG & DEBUG_ADVANCE) {
490  bu_log(
491  "%d, %d box push odist_corr=%.20e n=%.20e model_end=%.20e\n",
492  ap->a_x, ap->a_y,
493  ssp->odist_corr,
494  t0 /*ssp->dist_corr*/,
495  ssp->model_end);
496  bu_log(
497  "box_start o=%.20e n=%.20e\nbox_end o=%.20e n=%.20e\n",
498  ssp->obox_start,
499  ssp->box_start,
500  ssp->obox_end,
501  ssp->box_end);
502  bu_log("Point=(%g, %g, %g)\n",
503  px, py, pz);
504  VPRINT("Dir",
505  ssp->newray.r_dir);
506  rt_pr_cut(cutp, 0);
507  }
508 
509  /* Advance 1mm, or (square of) smallest value that
510  * hardware floating point resolution will allow.
511  */
512  fraction = frexp(ssp->box_end,
513  &exponent);
514 
515  if (RT_G_DEBUG & DEBUG_ADVANCE) {
516  bu_log(
517  "exp=%d, fraction=%.20e\n",
518  exponent, fraction);
519  }
520  fraction += SQRT_SMALL_FASTF;
521  delta = ldexp(fraction, exponent);
522  if (RT_G_DEBUG & DEBUG_ADVANCE) {
523  bu_log(
524  "ldexp: delta=%g, fract=%g, exp=%d\n",
525  delta,
526  fraction,
527  exponent);
528  }
529 
530  /* Never advance less than 1mm */
531  if (delta < 1) delta = 1.0;
532  ssp->box_start = ssp->box_end + delta;
533  ssp->box_end = ssp->box_start + delta;
534 
535  if (RT_G_DEBUG & DEBUG_ADVANCE) {
536  bu_log(
537  "push%d: was=%.20e, now=%.20e\n\n",
538  push_flag,
539  ssp->box_end,
540  ssp->box_start);
541  }
542  push_flag++;
543  if (push_flag > 3) {
544  bu_log("rt_advance_to_next_cell(): INTERNAL ERROR: infinite loop aborted, ray %d, %d truncated\n",
545  ap->a_x, ap->a_y);
546  goto pop_space_stack;
547  }
548  /* See if point marched outside model RPP */
549  if (ssp->box_start > ssp->model_end)
550  goto pop_space_stack;
551  t0 = ssp->box_start + ssp->ap->a_rt_i->rti_tol.dist;
552  goto top;
553  }
554  if (push_flag) {
555  push_flag = 0;
556  if (RT_G_DEBUG & DEBUG_ADVANCE) {
557  bu_log(
558  "%d, %d Escaped %d. dist_corr=%g, box_start=%g, box_end=%g\n",
559  ap->a_x, ap->a_y,
560  push_flag,
561  t0 /*ssp->dist_corr*/,
562  ssp->box_start,
563  ssp->box_end);
564  }
565  }
566  if (RT_G_DEBUG & DEBUG_ADVANCE) {
567  bu_log(
568  "rt_advance_to_next_cell()=%p lastcut=%p\n",
569  (void *)cutp, (void *)ssp->lastcut);
570  }
571 
572  ssp->newray.r_pt[X] = px;
573  ssp->newray.r_pt[Y] = py;
574  ssp->newray.r_pt[Z] = pz;
575  if (!rt_in_rpp(&ssp->newray, ssp->inv_dir,
576  cutp->bn.bn_min,
577  cutp->bn.bn_max)) {
578  bu_log("rt_advance_to_next_cell(): MISSED BOX\nrmin, rmax(%.20e, %.20e) box(%.20e, %.20e)\n",
579  ssp->newray.r_min,
580  ssp->newray.r_max,
581  ssp->box_start, ssp->box_end);
582  goto push_to_next_box;
583  }
584 
585  done_return_cutp: ssp->lastcut = cutp;
586  ssp->dist_corr = t0;
587  ssp->box_start = t0 + ssp->newray.r_min;
588  ssp->box_end = t0 + ssp->newray.r_max;
589  if (RT_G_DEBUG & DEBUG_ADVANCE) {
590  bu_log(
591  "rt_advance_to_next_cell() box=(%g, %g)\n",
592  ssp->box_start, ssp->box_end);
593  }
594  return cutp;
595  case CUT_NUGRIDNODE: {
596  struct rt_shootray_status *old;
597 
598  BU_ALLOC(old, struct rt_shootray_status);
599  *old = *ssp; /* struct copy */
600 
601  /* Descend into node */
602  ssp->old_status = old;
603  ssp->lastcut = ssp->lastcell = CUTTER_NULL;
604  curcut = ssp->curcut = cutp;
605  VSET(ssp->curmin, cutp->nugn.nu_axis[X][0].nu_spos,
606  cutp->nugn.nu_axis[Y][0].nu_spos,
607  cutp->nugn.nu_axis[Z][0].nu_spos);
611  break; }
612  case CUT_CUTNODE:
613  bu_bomb("rt_advance_to_next_cell: impossible: cutnote as leaf!");
614  break;
615  default:
616  bu_bomb("rt_advance_to_next_cell: unknown spt node");
617  break;
618  }
619 
620  /* Continue with the current space partitioning algorithm. */
621  }
622  /* NOTREACHED */
623 
624  /*
625  * If ray has escaped from model RPP, and there are infinite
626  * solids in the model, there is one more (special) BOXNODE for
627  * the caller to process.
628  */
629 escaped_from_model:
630  curcut = &ssp->ap->a_rt_i->rti_inf_box;
631  if (curcut->bn.bn_len <= 0 && curcut->bn.bn_piecelen <= 0)
632  curcut = CUTTER_NULL;
633  ssp->curcut = curcut;
634  return curcut;
635 }
636 
637 
638 /**
639  * This routine traces a ray from its start point to model exit
640  * through the space partitioning tree. The objective is to find all
641  * primitives that use "pieces" and also have a bounding box that
642  * extends behind the ray start point. The minimum (most negative)
643  * intersection with such a bounding box is returned. The "backbits"
644  * bit vector (provided by the caller) gets a bit set for every
645  * primitive that meets the above criteria. No primitive
646  * intersections are performed.
647  *
648  * XXX This routine does not work for NUGRID XXX
649  */
650 fastf_t
651 rt_find_backing_dist(struct rt_shootray_status *ss, struct bu_bitv *backbits) {
652  fastf_t min_backing_dist=BACKING_DIST;
653  fastf_t cur_dist=0.0;
654  point_t cur_pt;
655  union cutter *cutp;
656  struct bu_bitv *solidbits;
657  struct xray ray;
658  struct resource *resp;
659  struct rt_i *rtip;
660  size_t i;
661 
662  resp = ss->ap->a_resource;
663  rtip = ss->ap->a_rt_i;
664 
665  /* get a bit vector of our own to avoid duplicate bounding box
666  * intersection calculations
667  */
668  solidbits = rt_get_solidbitv(rtip->nsolids, resp);
669 
670  ray = ss->ap->a_ray; /* struct copy, don't mess with the original */
671 
672  /* cur_dist keeps track of where we are along the ray. stop when
673  * cur_dist reaches far intersection of ray and model bounding box
674  */
675  while (cur_dist <= ss->ap->a_ray.r_max) {
676  /* calculate the current point along the ray */
677  VJOIN1(cur_pt, ss->ap->a_ray.r_pt, cur_dist, ss->ap->a_ray.r_dir);
678 
679  /* descend into the space partitioning tree based on this
680  * point.
681  */
682  cutp = &ss->ap->a_rt_i->rti_CutHead;
683  while (cutp->cut_type == CUT_CUTNODE) {
684  if (cur_pt[cutp->cn.cn_axis] >= cutp->cn.cn_point) {
685  cutp=cutp->cn.cn_r;
686  } else {
687  cutp=cutp->cn.cn_l;
688  }
689  }
690 
691  /* we are now at the box node for the current point */
692  /* check if the ray intersects this box */
693  if (!rt_in_rpp(&ray, ss->inv_dir, cutp->bn.bn_min, cutp->bn.bn_max)) {
694  /* ray does not intersect this cell
695  * one of two situations must exist:
696  * 1. ray starts outside model bounding box (no need for these calculations)
697  * 2. we have proceeded beyond end of model bounding box (we are done)
698  * in either case, we are finished
699  */
700  goto done;
701  } else {
702  /* increment cur_dist into next cell for next execution of this loop */
703  cur_dist = ray.r_max + ss->ap->a_rt_i->rti_tol.dist;
704  }
705 
706  /* process this box node (look at all the pieces) */
707  for (i=0; i<cutp->bn.bn_piecelen; i++) {
708  struct rt_piecelist *plp=&cutp->bn.bn_piecelist[i];
709 
710  if (BU_BITTEST(solidbits, plp->stp->st_bit) == 0) {
711  /* we haven't looked at this primitive before */
712  if (rt_in_rpp(&ray, ss->inv_dir, plp->stp->st_min, plp->stp->st_max)) {
713  /* ray intersects this primitive bounding box */
714 
715  if (ray.r_min < BACKING_DIST) {
716  if (ray.r_min < min_backing_dist) {
717  /* move our backing distance back to catch this one */
718  min_backing_dist = ray.r_min;
719  }
720 
721  /* add this one to our list of primitives to check */
722  BU_BITSET(backbits, plp->stp->st_bit);
723  }
724  }
725  /* set bit so we don't repeat this calculation */
726  BU_BITSET(solidbits, plp->stp->st_bit);
727  }
728  }
729  }
730 done:
731  /* put our bit vector on the resource list */
732  BU_CK_BITV(solidbits);
733  BU_LIST_APPEND(&resp->re_solid_bitv, &solidbits->l);
734 
735  /* return our minimum backing distance */
736  return min_backing_dist;
737 }
738 
739 
740 /**
741  * Routines for plotting the progress of one ray through the model.
742  */
743 void
744 rt_3move_raydist(FILE *fp, struct xray *rayp, double dist)
745 {
746  point_t p;
747 
748  VJOIN1(p, rayp->r_pt, dist, rayp->r_dir);
749  pdv_3move(fp, p);
750 }
751 
752 
753 void
754 rt_3cont_raydist(FILE *fp, struct xray *rayp, double dist)
755 {
756  point_t p;
757 
758  VJOIN1(p, rayp->r_pt, dist, rayp->r_dir);
759  pdv_3cont(fp, p);
760 }
761 
762 
763 void
764 rt_plot_cell(const union cutter *cutp, const struct rt_shootray_status *ssp, struct bu_list *waiting_segs_hd, struct rt_i *rtip)
765 {
766  char buf[128];
767  static int fnum = 0;
768  FILE *fp;
769  struct soltab **stpp;
770  struct application *ap;
771 
772  RT_CK_RTI(rtip);
773  RT_AP_CHECK(ssp->ap);
774  RT_CK_RTI(ssp->ap->a_rt_i);
775  ap = ssp->ap;
776 
777  sprintf(buf, "cell%d.plot3", fnum++);
778  fp = fopen(buf, "wb");
779  if (fp == NULL) {
780  perror(buf);
781  } else {
782  /* green box for model RPP */
783  pl_color(fp, 0, 100, 0);
784 
785  /* Plot the model RPP, to provide some context */
786  pdv_3space(fp, rtip->rti_pmin, rtip->rti_pmax);
787  pdv_3box(fp, rtip->rti_pmin, rtip->rti_pmax);
788 
789  /* Plot the outline of this one cell */
790  pl_color(fp, 80, 80, 250);
791  switch (cutp->cut_type) {
792  case CUT_BOXNODE:
793  pdv_3box(fp, cutp->bn.bn_min, cutp->bn.bn_max);
794  break;
795  default:
796  bu_log("cut_type = %d\n", cutp->cut_type);
797  bu_bomb("Unknown cut_type\n");
798  }
799 
800  if (cutp->bn.bn_len > 0) {
801  /* Plot every solid listed in this cell */
802  stpp = &(cutp->bn.bn_list[cutp->bn.bn_len-1]);
803  for (; stpp >= cutp->bn.bn_list; stpp--) {
804  register struct soltab *stp = *stpp;
805 
806  rt_plot_solid(fp, rtip, stp, ap->a_resource);
807  }
808  }
809 
810  /* Plot interval of ray in box, in green */
811  pl_color(fp, 100, 255, 200);
812  rt_3move_raydist(fp, &ap->a_ray, ssp->box_start);
813  rt_3cont_raydist(fp, &ap->a_ray, ssp->box_end);
814 
815  if (bu_list_len(waiting_segs_hd) <= 0) {
816  /* No segments, just plot the whole ray */
817  pl_color(fp, 255, 255, 0); /* yellow -- no segs */
818  rt_3move_raydist(fp, &ap->a_ray, ssp->model_start);
819  rt_3cont_raydist(fp, &ap->a_ray, ssp->box_start);
820  rt_3move_raydist(fp, &ap->a_ray, ssp->box_end);
821  rt_3cont_raydist(fp, &ap->a_ray, ssp->model_end);
822  } else {
823  /* Plot the segments awaiting boolweave. */
824  struct seg *segp;
825 
826  for (BU_LIST_FOR(segp, seg, waiting_segs_hd)) {
827  RT_CK_SEG(segp);
828  pl_color(fp, 255, 0, 0); /* red */
829  rt_3move_raydist(fp, &ap->a_ray, segp->seg_in.hit_dist);
830  rt_3cont_raydist(fp, &ap->a_ray, segp->seg_out.hit_dist);
831  }
832  }
833 
834  fclose(fp);
835  bu_log("wrote %s\n", buf);
836  }
837 }
838 
839 
840 int
841 rt_shootray(register struct application *ap)
842 {
843  struct rt_shootray_status ss;
844  struct seg new_segs; /* from solid intersections */
845  struct seg waiting_segs; /* awaiting rt_boolweave() */
846  struct seg finished_segs; /* processed by rt_boolweave() */
847  fastf_t last_bool_start;
848  struct bu_bitv *solidbits; /* bits for all solids shot so far */
849  struct bu_bitv *backbits=NULL; /* bits for all solids using pieces that need to be intersected behind
850  the ray start point */
851  struct bu_ptbl *regionbits; /* table of all involved regions */
852  char *status;
853  struct partition InitialPart; /* Head of Initial Partitions */
854  struct partition FinalPart; /* Head of Final Partitions */
855  struct soltab **stpp;
856  register const union cutter *cutp;
857  struct resource *resp;
858  struct rt_i *rtip;
859  const int debug_shoot = RT_G_DEBUG & DEBUG_SHOOT;
860  fastf_t pending_hit = 0; /* dist of closest odd hit pending */
861 
862  RT_AP_CHECK(ap);
863  if (ap->a_magic) {
864  RT_CK_AP(ap);
865  } else {
866  ap->a_magic = RT_AP_MAGIC;
867  }
868  if (ap->a_ray.magic) {
869  RT_CK_RAY(&(ap->a_ray));
870  } else {
871  ap->a_ray.magic = RT_RAY_MAGIC;
872  }
873 #if 0
874  /* FIXME: resolve the rt_uniresource initialization dilemma, wrap
875  * in a function so we can call rt_init_resource().
876  */
877  if (!ap->a_resource) {
878  ap->a_resource = &rt_uniresource;
879  }
881 #endif
882  if (ap->a_resource == RESOURCE_NULL) {
883  ap->a_resource = &rt_uniresource;
884  if (RT_G_DEBUG)
885  bu_log("rt_shootray: defaulting a_resource to &rt_uniresource\n");
886  if (rt_uniresource.re_magic == 0)
888  }
889  ss.ap = ap;
890  rtip = ap->a_rt_i;
891  RT_CK_RTI(rtip);
892  resp = ap->a_resource;
893  RT_CK_RESOURCE(resp);
894  ss.resp = resp;
895 
896  if (RT_G_DEBUG) {
897  /* only test extensively if something in run-time debug is enabled */
900  bu_log("\n**********shootray cpu=%d %d, %d lvl=%d a_onehit=%d (%s)\n",
901  resp->re_cpu,
902  ap->a_x, ap->a_y,
903  ap->a_level,
904  ap->a_onehit,
905  ap->a_purpose != (char *)0 ? ap->a_purpose : "?");
906  bu_log("Pnt (%.20G, %.20G, %.20G)\nDir (%.20G, %.20G, %.20G)\n",
907  V3ARGS(ap->a_ray.r_pt),
908  V3ARGS(ap->a_ray.r_dir));
909  }
910  }
911 
912  if (rtip->needprep)
913  rt_prep_parallel(rtip, 1); /* Stay on our CPU */
914 
915  InitialPart.pt_forw = InitialPart.pt_back = &InitialPart;
916  InitialPart.pt_magic = PT_HD_MAGIC;
917  FinalPart.pt_forw = FinalPart.pt_back = &FinalPart;
918  FinalPart.pt_magic = PT_HD_MAGIC;
919  ap->a_Final_Part_hdp = &FinalPart;
920 
921  BU_LIST_INIT(&new_segs.l);
922  BU_LIST_INIT(&waiting_segs.l);
923  BU_LIST_INIT(&finished_segs.l);
924  ap->a_finished_segs_hdp = &finished_segs;
925 
926  if (!BU_LIST_IS_INITIALIZED(&resp->re_parthead)) {
927  /* XXX This shouldn't happen any more */
928  bu_log("rt_shootray() resp=%p uninitialized, fixing it\n", (void *)resp);
929  /*
930  * We've been handed a mostly un-initialized resource struct,
931  * with only a magic number and a cpu number filled in. Init
932  * it and add it to the table. This is how
933  * application-provided resource structures are remembered for
934  * later cleanup by the library.
935  */
936  rt_init_resource(resp, resp->re_cpu, rtip);
937  }
938  /* Ensure that this CPU's resource structure is registered */
939  if (resp != &rt_uniresource)
940  BU_ASSERT_PTR(BU_PTBL_GET(&rtip->rti_resources, resp->re_cpu), !=, NULL);
941 
942  solidbits = rt_get_solidbitv(rtip->nsolids, resp);
943 
944  if (BU_LIST_IS_EMPTY(&resp->re_region_ptbl)) {
945  BU_ALLOC(regionbits, struct bu_ptbl);
946  bu_ptbl_init(regionbits, 7, "rt_shootray() regionbits ptbl");
947  } else {
948  regionbits = BU_LIST_FIRST(bu_ptbl, &resp->re_region_ptbl);
949  BU_LIST_DEQUEUE(&regionbits->l);
950  BU_CK_PTBL(regionbits);
951  }
952 
953  if (!resp->re_pieces && rtip->rti_nsolids_with_pieces > 0) {
954  /* Initialize this processors 'solid pieces' state */
955  rt_res_pieces_init(resp, rtip);
956  }
958  /* only happens first time through */
959  bu_ptbl_init(&resp->re_pieces_pending, 100, "re_pieces_pending");
960  }
962 
963  /* Verify that direction vector has unit length */
964  if (RT_G_DEBUG) {
965  fastf_t f, diff;
966 
967  f = MAGSQ(ap->a_ray.r_dir);
968  if (NEAR_ZERO(f, ap->a_rt_i->rti_tol.dist)) {
969  bu_bomb("rt_shootray: zero length dir vector\n");
970  }
971  diff = f - 1;
972  if (!NEAR_ZERO(diff, ap->a_rt_i->rti_tol.dist)) {
973  bu_log("rt_shootray: non-unit dir vect (x%d y%d lvl%d)\n",
974  ap->a_x, ap->a_y, ap->a_level);
975  f = 1/f;
976  VSCALE(ap->a_ray.r_dir, ap->a_ray.r_dir, f);
977  }
978  }
979 
980  /*
981  * Record essential statistics in per-processor data structure.
982  */
983  resp->re_nshootray++;
984 
985  /* Compute the inverse of the direction cosines */
986  if (ap->a_ray.r_dir[X] < -SQRT_SMALL_FASTF) {
987  ss.abs_inv_dir[X] = -(ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
988  ss.rstep[X] = -1;
989  } else if (ap->a_ray.r_dir[X] > SQRT_SMALL_FASTF) {
990  ss.abs_inv_dir[X] = (ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
991  ss.rstep[X] = 1;
992  } else {
993  ap->a_ray.r_dir[X] = 0.0;
994  ss.abs_inv_dir[X] = ss.inv_dir[X] = INFINITY;
995  ss.rstep[X] = 0;
996  }
997  if (ap->a_ray.r_dir[Y] < -SQRT_SMALL_FASTF) {
998  ss.abs_inv_dir[Y] = -(ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
999  ss.rstep[Y] = -1;
1000  } else if (ap->a_ray.r_dir[Y] > SQRT_SMALL_FASTF) {
1001  ss.abs_inv_dir[Y] = (ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
1002  ss.rstep[Y] = 1;
1003  } else {
1004  ap->a_ray.r_dir[Y] = 0.0;
1005  ss.abs_inv_dir[Y] = ss.inv_dir[Y] = INFINITY;
1006  ss.rstep[Y] = 0;
1007  }
1008  if (ap->a_ray.r_dir[Z] < -SQRT_SMALL_FASTF) {
1009  ss.abs_inv_dir[Z] = -(ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
1010  ss.rstep[Z] = -1;
1011  } else if (ap->a_ray.r_dir[Z] > SQRT_SMALL_FASTF) {
1012  ss.abs_inv_dir[Z] = (ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
1013  ss.rstep[Z] = 1;
1014  } else {
1015  ap->a_ray.r_dir[Z] = 0.0;
1016  ss.abs_inv_dir[Z] = ss.inv_dir[Z] = INFINITY;
1017  ss.rstep[Z] = 0;
1018  }
1019  VMOVE(ap->a_inv_dir, ss.inv_dir);
1020 
1021  /*
1022  * If ray does not enter the model RPP, skip on. If ray ends
1023  * exactly at the model RPP, trace it.
1024  */
1025  if (!rt_in_rpp(&ap->a_ray, ss.inv_dir, rtip->mdl_min, rtip->mdl_max) ||
1026  ap->a_ray.r_max < 0.0) {
1027  cutp = &ap->a_rt_i->rti_inf_box;
1028  if (cutp->bn.bn_len > 0) {
1029  /* Model has infinite solids, need to fire at them. */
1030  ss.box_start = BACKING_DIST;
1031  ss.model_start = 0;
1032  ss.box_end = ss.model_end = INFINITY;
1033  ss.lastcut = CUTTER_NULL;
1034  ss.old_status = (struct rt_shootray_status *)NULL;
1035  ss.curcut = cutp;
1036  ss.lastcell = ss.curcut;
1037  VMOVE(ss.curmin, rtip->mdl_min);
1038  VMOVE(ss.curmax, rtip->mdl_max);
1039  last_bool_start = BACKING_DIST;
1040  ss.newray = ap->a_ray; /* struct copy */
1041  ss.odist_corr = ss.obox_start = ss.obox_end = -99;
1042  ss.dist_corr = 0.0;
1043  goto start_cell;
1044  }
1045  resp->re_nmiss_model++;
1046  if (ap->a_miss)
1047  ap->a_return = ap->a_miss(ap);
1048  else
1049  ap->a_return = 0;
1050  status = "MISS model";
1051  goto out;
1052  }
1053 
1054  /*
1055  * The interesting part of the ray starts at distance 0. If the
1056  * ray enters the model at a negative distance, (i.e., the ray
1057  * starts within the model RPP), we only look at little bit behind
1058  * (BACKING_DIST) to see if we are just coming out of something,
1059  * but never further back than the intersection with the model
1060  * RPP. If the ray enters the model at a positive distance, we
1061  * always start there. It is vital that we never pick a start
1062  * point outside the model RPP, or the space partitioning tree
1063  * will pick the wrong box and the ray will miss it.
1064  *
1065  * BACKING_DIST should probably be determined by floating point
1066  * noise factor due to model RPP size -vs- number of bits of
1067  * floating point mantissa significance, rather than a constant,
1068  * but that is too hideous to think about here. Also note that
1069  * applications that really depend on knowing what region they are
1070  * leaving from should probably back their own start-point up,
1071  * rather than depending on it here, but it isn't much trouble
1072  * here.
1073  *
1074  * Modification by JRA for pieces methodology:
1075  *
1076  * The original algorithm here assumed that if we encountered any
1077  * primitive along the positive direction of the ray, ALL its
1078  * intersections would be calculated. With pieces, we may see
1079  * only an exit hit if the entrance piece is in a space partition
1080  * cell that is more than "BACKING_DIST" behind the ray start
1081  * point (leading to incorrect results). I have modified the
1082  * setting of "ss.box_start" (when pieces are present and the ray
1083  * start point is inside the model bounding box) as follows (see
1084  * rt_find_backing_dist()):
1085  *
1086  * - The ray is traced through the space partitioning tree.
1087  *
1088  * - The ray is intersected with the bounding box of each
1089  * primitive using pieces in each cell
1090  *
1091  * - The minimum of all these intersections is set as the initial
1092  * "ss.box_start".
1093  *
1094  * - The "backbits" bit vector has a bit set for each of the
1095  * primitives using pieces that have bounding boxes that extend
1096  * behind the ray start point
1097  *
1098  * Further below (in the "pieces" loop), I have added code to
1099  * ignore primitives that do not have a bit set in the backbits
1100  * vector when we are behind the ray start point.
1101  */
1102 
1103  /* these two values set the point where the ray tracing actually
1104  * begins and ends
1105  */
1106  ss.box_start = ss.model_start = ap->a_ray.r_min;
1107  ss.box_end = ss.model_end = ap->a_ray.r_max;
1108 
1109  if (ap->a_rt_i->rti_nsolids_with_pieces > 0) {
1110  /* pieces are present */
1111  if (ss.box_start < BACKING_DIST) {
1112  /* the first ray intersection with the model bounding box
1113  * is more than BACKING_DIST behind the ray start point
1114  */
1115 
1116  /* get a bit vector to keep track of which primitives need
1117  * to be intersected behind the ray start point (those
1118  * having bounding boxes extending behind the ray start
1119  * point and using pieces)
1120  */
1121  backbits = rt_get_solidbitv(rtip->nsolids, resp);
1122 
1123  /* call "rt_find_backing_dist()" to calculate the required
1124  * start point for calculation, and to fill in the
1125  * "backbits" bit vector
1126  */
1127  ss.box_start = rt_find_backing_dist(&ss, backbits);
1128  }
1129  } else {
1130  /* no pieces present, use the old scheme */
1131  if (ss.box_start < BACKING_DIST)
1132  ss.box_start = BACKING_DIST; /* Only look a little bit behind */
1133  }
1134 
1135  ss.lastcut = CUTTER_NULL;
1136  ss.old_status = (struct rt_shootray_status *)NULL;
1137  ss.curcut = &ap->a_rt_i->rti_CutHead;
1138  if (ss.curcut->cut_type == CUT_NUGRIDNODE) {
1139  ss.lastcell = CUTTER_NULL;
1140  VSET(ss.curmin, ss.curcut->nugn.nu_axis[X][0].nu_spos,
1141  ss.curcut->nugn.nu_axis[Y][0].nu_spos,
1142  ss.curcut->nugn.nu_axis[Z][0].nu_spos);
1146  } else if (ss.curcut->cut_type == CUT_CUTNODE ||
1147  ss.curcut->cut_type == CUT_BOXNODE) {
1148  ss.lastcell = ss.curcut;
1149  VMOVE(ss.curmin, rtip->mdl_min);
1150  VMOVE(ss.curmax, rtip->mdl_max);
1151  }
1152 
1153  last_bool_start = BACKING_DIST;
1154  ss.newray = ap->a_ray; /* struct copy */
1155  ss.odist_corr = ss.obox_start = ss.obox_end = -99;
1156  ss.dist_corr = 0.0;
1157  ss.box_num = 0;
1158 
1159  /*
1160  * While the ray remains inside model space, push from box to box
1161  * until ray emerges from model space again (or first hit is
1162  * found, if user is impatient). It is vitally important to
1163  * always stay within the model RPP, or the space partitioning tree
1164  * will pick wrong boxes & miss them.
1165  */
1166  while ((cutp = rt_advance_to_next_cell(&ss)) != CUTTER_NULL) {
1167  start_cell:
1168  if (debug_shoot) {
1169  bu_log("BOX #%d interval is %g..%g\n", ss.box_num, ss.box_start, ss.box_end);
1170  rt_pr_cut(cutp, 0);
1171  }
1172 
1173  if (cutp->bn.bn_len <= 0 && cutp->bn.bn_piecelen <= 0) {
1174  /* Push ray onwards to next box */
1175  ss.box_start = ss.box_end;
1176  resp->re_nempty_cells++;
1177  continue;
1178  }
1179 
1180  /* Consider all "pieces" of all solids within the box */
1181  pending_hit = ss.box_end;
1182  if (cutp->bn.bn_piecelen > 0) {
1183  register struct rt_piecelist *plp;
1184 
1185  plp = &(cutp->bn.bn_piecelist[cutp->bn.bn_piecelen-1]);
1186  for (; plp >= cutp->bn.bn_piecelist; plp--) {
1187  struct rt_piecestate *psp;
1188  struct soltab *stp;
1189  int ret;
1190  int had_hits_before;
1191 
1192  RT_CK_PIECELIST(plp);
1193 
1194  /* Consider all pieces of this one solid in this
1195  * cell.
1196  */
1197  stp = plp->stp;
1198  RT_CK_SOLTAB(stp);
1199 
1200  if (backbits && ss.box_end < BACKING_DIST && BU_BITTEST(backbits, stp->st_bit) == 0) {
1201  /* we are behind the ray start point and this
1202  * primitive is not one that we need to intersect
1203  * back here.
1204  */
1205  continue;
1206  }
1207 
1208  psp = &(resp->re_pieces[stp->st_piecestate_num]);
1209  RT_CK_PIECESTATE(psp);
1210  if (psp->ray_seqno != resp->re_nshootray) {
1211  /* state is from an earlier ray, scrub */
1212  BU_BITV_ZEROALL(psp->shot);
1213  psp->ray_seqno = resp->re_nshootray;
1214  rt_htbl_reset(&psp->htab);
1215 
1216  /* Compute ray entry and exit to entire solid's
1217  * bounding box.
1218  */
1219  if (!rt_in_rpp(&ss.newray, ss.inv_dir,
1220  stp->st_min, stp->st_max)) {
1221  if (debug_shoot)bu_log("rpp miss %s (all pieces)\n", stp->st_name);
1222  resp->re_prune_solrpp++;
1223  BU_BITSET(solidbits, stp->st_bit);
1224  continue; /* MISS */
1225  }
1226  psp->mindist = ss.newray.r_min + ss.dist_corr;
1227  psp->maxdist = ss.newray.r_max + ss.dist_corr;
1228  if (debug_shoot) bu_log("%s mindist=%g, maxdist=%g\n", stp->st_name, psp->mindist, psp->maxdist);
1229  had_hits_before = 0;
1230  } else {
1231  if (BU_BITTEST(solidbits, stp->st_bit)) {
1232  /* we missed the solid RPP in an earlier cell */
1233  resp->re_ndup++;
1234  continue; /* already shot */
1235  }
1236  had_hits_before = psp->htab.end;
1237  }
1238 
1239  /*
1240  * Allow this solid to shoot at all of its 'pieces' in
1241  * this cell, all at once. 'newray' has been
1242  * transformed to be near to this cell, and
1243  * 'dist_corr' is the additive correction factor that
1244  * ft_piece_shot() must apply to hits calculated using
1245  * 'newray'.
1246  */
1247  resp->re_piece_shots++;
1248  psp->cutp = cutp;
1249 
1250  ret = -1;
1251  if (stp->st_meth->ft_piece_shot) {
1252  ret = stp->st_meth->ft_piece_shot(psp, plp, ss.dist_corr, &ss.newray, ap, &waiting_segs);
1253  }
1254  if (ret <= 0) {
1255  /* No hits at all */
1256  resp->re_piece_shot_miss++;
1257  } else {
1258  resp->re_piece_shot_hit++;
1259  }
1260  if (debug_shoot)bu_log("shooting %s pieces, nhit=%d\n", stp->st_name, ret);
1261 
1262  /* See if this solid has been fully processed yet. If
1263  * ray has passed through bounding volume, we're done.
1264  * ft_piece_hitsegs() will only be called once per
1265  * ray.
1266  */
1267  if (ss.box_end > psp->maxdist && psp->htab.end > 0) {
1268  /* Convert hits into segs */
1269  if (debug_shoot)bu_log("shooting %s pieces complete, making segs\n", stp->st_name);
1270  /* Distance correction was handled in ft_piece_shot */
1271  if (stp->st_meth->ft_piece_hitsegs)
1272  stp->st_meth->ft_piece_hitsegs(psp, &waiting_segs, ap);
1273  rt_htbl_reset(&psp->htab);
1274  BU_BITSET(solidbits, stp->st_bit);
1275 
1276  if (had_hits_before)
1277  bu_ptbl_rm(&resp->re_pieces_pending, (long *)psp);
1278  } else {
1279  if (!had_hits_before)
1280  bu_ptbl_ins_unique(&resp->re_pieces_pending, (long *)psp);
1281  }
1282  }
1283  }
1284 
1285  /* Consider all solids within the box */
1286  if (cutp->bn.bn_len > 0 && ss.box_end >= BACKING_DIST) {
1287  stpp = &(cutp->bn.bn_list[cutp->bn.bn_len-1]);
1288  for (; stpp >= cutp->bn.bn_list; stpp--) {
1289  register struct soltab *stp = *stpp;
1290  int ret;
1291 
1292  if (BU_BITTEST(solidbits, stp->st_bit)) {
1293  resp->re_ndup++;
1294  continue; /* already shot */
1295  }
1296 
1297  /* Shoot a ray */
1298  BU_BITSET(solidbits, stp->st_bit);
1299 
1300  /* Check against bounding RPP, if desired by solid */
1301  if (stp->st_meth->ft_use_rpp) {
1302  if (!rt_in_rpp(&ss.newray, ss.inv_dir,
1303  stp->st_min, stp->st_max)) {
1304  if (debug_shoot)bu_log("rpp miss %s\n", stp->st_name);
1305  resp->re_prune_solrpp++;
1306  continue; /* MISS */
1307  }
1308  if (ss.dist_corr + ss.newray.r_max < BACKING_DIST) {
1309  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);
1310  resp->re_prune_solrpp++;
1311  continue; /* MISS */
1312  }
1313  }
1314 
1315  if (debug_shoot)bu_log("shooting %s\n", stp->st_name);
1316  resp->re_shots++;
1317  BU_LIST_INIT(&(new_segs.l));
1318 
1319  ret = -1;
1320  if (stp->st_meth->ft_shot) {
1321  ret = stp->st_meth->ft_shot(stp, &ss.newray, ap, &new_segs);
1322  }
1323  if (ret <= 0) {
1324  resp->re_shot_miss++;
1325  continue; /* MISS */
1326  }
1327 
1328  /* Add seg chain to list awaiting rt_boolweave() */
1329  {
1330  register struct seg *s2;
1331  while (BU_LIST_WHILE(s2, seg, &(new_segs.l))) {
1332  BU_LIST_DEQUEUE(&(s2->l));
1333  /* Restore to original distance */
1334  s2->seg_in.hit_dist += ss.dist_corr;
1335  s2->seg_out.hit_dist += ss.dist_corr;
1336  s2->seg_in.hit_rayp = s2->seg_out.hit_rayp = &ap->a_ray;
1337  BU_LIST_INSERT(&(waiting_segs.l), &(s2->l));
1338  }
1339  }
1340  resp->re_shot_hit++;
1341  }
1342  }
1343  if (RT_G_DEBUG & DEBUG_ADVANCE)
1344  rt_plot_cell(cutp, &ss, &(waiting_segs.l), rtip);
1345 
1346  /*
1347  * If a_onehit == 0 and a_ray_length <= 0, then the ray is
1348  * traced to +infinity.
1349  *
1350  * If a_onehit != 0, then it indicates how many hit points
1351  * (which are greater than the ray start point of 0.0) the
1352  * application requires, i.e., partitions with inhit >= 0. (If
1353  * negative, indicates number of non-air hits needed). If
1354  * this box yielded additional segments, immediately weave
1355  * them into the partition list, and perform final boolean
1356  * evaluation. If this results in the required number of
1357  * final partitions, then cease ray-tracing and hand the
1358  * partitions over to the application. All partitions will
1359  * have valid in and out distances. a_ray_length is treated
1360  * similarly to a_onehit.
1361  */
1362  if (BU_LIST_NON_EMPTY(&(waiting_segs.l))) {
1363  if (debug_shoot) {
1364  struct seg *segp;
1365  bu_log("Waiting segs:\n");
1366  for (BU_LIST_FOR(segp, seg, &(waiting_segs.l))) {
1367  rt_pr_seg(segp);
1368  }
1369  }
1370  if (ap->a_onehit != 0) {
1371  int done;
1372 
1373  /* Weave these segments into partition list */
1374  rt_boolweave(&finished_segs, &waiting_segs, &InitialPart, ap);
1375 
1376  if (BU_PTBL_LEN(&resp->re_pieces_pending) > 0) {
1377 
1378  /* Find the lowest pending mindist, that's as far
1379  * as boolfinal can progress to.
1380  */
1381  struct rt_piecestate **psp;
1382  for (BU_PTBL_FOR(psp, (struct rt_piecestate **), &resp->re_pieces_pending)) {
1383  register fastf_t dist;
1384 
1385  dist = (*psp)->mindist;
1386  BU_ASSERT_DOUBLE(dist, <, INFINITY);
1387  if (dist < pending_hit) {
1388  pending_hit = dist;
1389  if (debug_shoot) bu_log("pending_hit lowered to %g by %s\n", pending_hit, (*psp)->stp->st_name);
1390  }
1391  }
1392  }
1393 
1394  /* Evaluate regions up to end of good segs */
1395  if (ss.box_end < pending_hit) pending_hit = ss.box_end;
1396  done = rt_boolfinal(&InitialPart, &FinalPart,
1397  last_bool_start, pending_hit, regionbits, ap, solidbits);
1398  last_bool_start = pending_hit;
1399 
1400  /* See if enough partitions have been acquired */
1401  if (done > 0) goto hitit;
1402  }
1403  }
1404 
1405  if (ap->a_ray_length > 0.0 &&
1406  ss.box_end >= ap->a_ray_length &&
1407  ap->a_ray_length < pending_hit)
1408  goto weave;
1409 
1410  /* Push ray onwards to next box */
1411  ss.box_start = ss.box_end;
1412  }
1413 
1414  /*
1415  * Ray has finally left known space -- Weave any remaining
1416  * segments into the partition list.
1417  */
1418 weave:
1420  bu_log("rt_shootray: ray has left known space\n");
1421 
1422  /* Process any pending hits into segs */
1423  if (BU_PTBL_LEN(&resp->re_pieces_pending) > 0) {
1424  struct rt_piecestate **psp;
1425  for (BU_PTBL_FOR(psp, (struct rt_piecestate **), &resp->re_pieces_pending)) {
1426  if ((*psp)->htab.end > 0) {
1427  /* Convert any pending hits into segs */
1428  /* Distance correction was handled in ft_piece_shot */
1429  (*psp)->stp->st_meth->ft_piece_hitsegs(*psp, &waiting_segs, ap);
1430  rt_htbl_reset(&(*psp)->htab);
1431  }
1432  *psp = NULL;
1433  }
1435  }
1436 
1437  if (BU_LIST_NON_EMPTY(&(waiting_segs.l))) {
1438  rt_boolweave(&finished_segs, &waiting_segs, &InitialPart, ap);
1439  }
1440 
1441  /* finished_segs chain now has all segments hit by this ray */
1442  if (BU_LIST_IS_EMPTY(&(finished_segs.l))) {
1443  if (ap->a_miss)
1444  ap->a_return = ap->a_miss(ap);
1445  else
1446  ap->a_return = 0;
1447  status = "MISS primitives";
1448  goto out;
1449  }
1450 
1451  /*
1452  * All intersections of the ray with the model have been computed.
1453  * Evaluate the boolean trees over each partition.
1454  */
1455  (void)rt_boolfinal(&InitialPart, &FinalPart, BACKING_DIST,
1456  INFINITY,
1457  regionbits, ap, solidbits);
1458 
1459  if (FinalPart.pt_forw == &FinalPart) {
1460  if (ap->a_miss)
1461  ap->a_return = ap->a_miss(ap);
1462  else
1463  ap->a_return = 0;
1464  status = "MISS bool";
1465  RT_FREE_PT_LIST(&InitialPart, resp);
1466  RT_FREE_SEG_LIST(&finished_segs, resp);
1467  goto out;
1468  }
1469 
1470  /*
1471  * Ray/model intersections exist. Pass the list to the user's
1472  * a_hit() routine. Note that only the hit_dist elements of
1473  * pt_inhit and pt_outhit have been computed yet. To compute both
1474  * hit_point and hit_normal, use the
1475  *
1476  * RT_HIT_NORMAL(NULL, hitp, stp, rayp, 0);
1477  *
1478  * macro. To compute just hit_point, use
1479  *
1480  * VJOIN1(hitp->hit_point, rp->r_pt, hitp->hit_dist, rp->r_dir);
1481  */
1482 hitit:
1483  if (debug_shoot) rt_pr_partitions(rtip, &FinalPart, "a_hit()");
1484 
1485  /*
1486  * Before recursing, release storage for unused Initial
1487  * partitions. finished_segs can not be released yet, because
1488  * FinalPart partitions will point to hits in those segments.
1489  */
1490  RT_FREE_PT_LIST(&InitialPart, resp);
1491 
1492  /*
1493  * finished_segs is only used by special hit routines which don't
1494  * follow the traditional solid modeling paradigm.
1495  */
1496  if (RT_G_DEBUG&DEBUG_ALLHITS) rt_pr_partitions(rtip, &FinalPart, "Partition list passed to a_hit() routine");
1497  if (ap->a_hit) {
1498  ap->a_return = ap->a_hit(ap, &FinalPart, &finished_segs);
1499  status = "HIT";
1500  } else {
1501  ap->a_return = 0;
1502  status = "MISS (unexpected)";
1503  }
1504 
1505  RT_FREE_SEG_LIST(&finished_segs, resp);
1506  RT_FREE_PT_LIST(&FinalPart, resp);
1507 
1508  /*
1509  * Processing of this ray is complete.
1510  */
1511 out:
1512  /* Return dynamic resources to their freelists. */
1513  BU_CK_BITV(solidbits);
1514  BU_LIST_APPEND(&resp->re_solid_bitv, &solidbits->l);
1515  if (backbits) {
1516  BU_CK_BITV(backbits);
1517  BU_LIST_APPEND(&resp->re_solid_bitv, &backbits->l);
1518  }
1519  BU_CK_PTBL(regionbits);
1520  BU_LIST_APPEND(&resp->re_region_ptbl, &regionbits->l);
1521 
1522  /* Clean up any pending hits */
1523  if (BU_PTBL_LEN(&resp->re_pieces_pending) > 0) {
1524  struct rt_piecestate **psp;
1525  for (BU_PTBL_FOR(psp, (struct rt_piecestate **), &resp->re_pieces_pending)) {
1526  if ((*psp)->htab.end > 0)
1527  rt_htbl_reset(&(*psp)->htab);
1528  }
1530  }
1531 
1532  /* Terminate any logging */
1533  if (RT_G_DEBUG&(DEBUG_ALLRAYS|DEBUG_SHOOT|DEBUG_PARTITION|DEBUG_ALLHITS)) {
1534  bu_log_indent_delta(-2);
1535  bu_log("----------shootray cpu=%d %d, %d lvl=%d (%s) %s ret=%d\n",
1536  resp->re_cpu,
1537  ap->a_x, ap->a_y,
1538  ap->a_level,
1539  ap->a_purpose != (char *)0 ? ap->a_purpose : "?",
1540  status, ap->a_return);
1541  }
1542  return ap->a_return;
1543 }
1544 
1545 
1546 const union cutter *
1547 rt_cell_n_on_ray(register struct application *ap, int n)
1548 
1549 /* First cell is #0 */
1550 {
1551  struct rt_shootray_status ss;
1552  register const union cutter *cutp;
1553  struct resource *resp;
1554  struct rt_i *rtip;
1555  const int debug_shoot = RT_G_DEBUG & DEBUG_SHOOT;
1556 
1557  memset(&ss, 0, sizeof(struct rt_shootray_status));
1558 
1559  RT_AP_CHECK(ap);
1560  if (ap->a_magic) {
1561  RT_CK_AP(ap);
1562  } else {
1563  ap->a_magic = RT_AP_MAGIC;
1564  }
1565  if (ap->a_ray.magic) {
1566  RT_CK_RAY(&(ap->a_ray));
1567  } else {
1568  ap->a_ray.magic = RT_RAY_MAGIC;
1569  }
1570  if (ap->a_resource == RESOURCE_NULL) {
1571  ap->a_resource = &rt_uniresource;
1572  if (RT_G_DEBUG)bu_log("rt_cell_n_on_ray: defaulting a_resource to &rt_uniresource\n");
1573  if (rt_uniresource.re_magic == 0)
1575  }
1576  ss.ap = ap;
1577  rtip = ap->a_rt_i;
1578  RT_CK_RTI(rtip);
1579  resp = ap->a_resource;
1580  RT_CK_RESOURCE(resp);
1581  ss.resp = resp;
1582 
1585  bu_log("\n**********cell_n_on_ray cpu=%d %d, %d lvl=%d (%s), n=%d\n",
1586  resp->re_cpu,
1587  ap->a_x, ap->a_y,
1588  ap->a_level,
1589  ap->a_purpose != (char *)0 ? ap->a_purpose : "?", n);
1590  bu_log("Pnt (%g, %g, %g) a_onehit=%d\n",
1591  V3ARGS(ap->a_ray.r_pt),
1592  ap->a_onehit);
1593  VPRINT("Dir", ap->a_ray.r_dir);
1594  }
1595 
1596  if (rtip->needprep)
1597  rt_prep_parallel(rtip, 1); /* Stay on our CPU */
1598 
1599  if (!BU_LIST_IS_INITIALIZED(&resp->re_parthead)) {
1600  /* XXX This shouldn't happen any more */
1601  bu_log("rt_cell_n_on_ray() resp=%p uninitialized, fixing it\n", (void *)resp);
1602  /*
1603  * We've been handed a mostly un-initialized resource struct,
1604  * with only a magic number and a cpu number filled in.
1605  * Init it and add it to the table.
1606  * This is how application-provided resource structures
1607  * are remembered for later cleanup by the library.
1608  */
1609  rt_init_resource(resp, resp->re_cpu, rtip);
1610  }
1611  /* Ensure that this CPU's resource structure is registered */
1612  BU_ASSERT_PTR(BU_PTBL_GET(&rtip->rti_resources, resp->re_cpu), !=, NULL);
1613 
1614  /* Verify that direction vector has unit length */
1615  if (RT_G_DEBUG) {
1616  fastf_t f, diff;
1617 
1618  f = MAGSQ(ap->a_ray.r_dir);
1619  if (NEAR_ZERO(f, ap->a_rt_i->rti_tol.dist)) {
1620  bu_bomb("rt_cell_n_on_ray: zero length dir vector\n");
1621  }
1622  diff = f - 1;
1623  if (!NEAR_ZERO(diff, ap->a_rt_i->rti_tol.dist)) {
1624  bu_log("rt_cell_n_on_ray: non-unit dir vect (x%d y%d lvl%d)\n",
1625  ap->a_x, ap->a_y, ap->a_level);
1626  f = 1/f;
1627  VSCALE(ap->a_ray.r_dir, ap->a_ray.r_dir, f);
1628  }
1629  }
1630 
1631  /* Compute the inverse of the direction cosines */
1632  if (ap->a_ray.r_dir[X] < -SQRT_SMALL_FASTF) {
1633  ss.abs_inv_dir[X] = -(ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
1634  ss.rstep[X] = -1;
1635  } else if (ap->a_ray.r_dir[X] > SQRT_SMALL_FASTF) {
1636  ss.abs_inv_dir[X] = (ss.inv_dir[X]=1.0/ap->a_ray.r_dir[X]);
1637  ss.rstep[X] = 1;
1638  } else {
1639  ap->a_ray.r_dir[X] = 0.0;
1640  ss.abs_inv_dir[X] = ss.inv_dir[X] = INFINITY;
1641  ss.rstep[X] = 0;
1642  }
1643  if (ap->a_ray.r_dir[Y] < -SQRT_SMALL_FASTF) {
1644  ss.abs_inv_dir[Y] = -(ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
1645  ss.rstep[Y] = -1;
1646  } else if (ap->a_ray.r_dir[Y] > SQRT_SMALL_FASTF) {
1647  ss.abs_inv_dir[Y] = (ss.inv_dir[Y]=1.0/ap->a_ray.r_dir[Y]);
1648  ss.rstep[Y] = 1;
1649  } else {
1650  ap->a_ray.r_dir[Y] = 0.0;
1651  ss.abs_inv_dir[Y] = ss.inv_dir[Y] = INFINITY;
1652  ss.rstep[Y] = 0;
1653  }
1654  if (ap->a_ray.r_dir[Z] < -SQRT_SMALL_FASTF) {
1655  ss.abs_inv_dir[Z] = -(ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
1656  ss.rstep[Z] = -1;
1657  } else if (ap->a_ray.r_dir[Z] > SQRT_SMALL_FASTF) {
1658  ss.abs_inv_dir[Z] = (ss.inv_dir[Z]=1.0/ap->a_ray.r_dir[Z]);
1659  ss.rstep[Z] = 1;
1660  } else {
1661  ap->a_ray.r_dir[Z] = 0.0;
1662  ss.abs_inv_dir[Z] = ss.inv_dir[Z] = INFINITY;
1663  ss.rstep[Z] = 0;
1664  }
1665 
1666  /*
1667  * If ray does not enter the model RPP, skip on.
1668  * If ray ends exactly at the model RPP, trace it.
1669  */
1670  if (!rt_in_rpp(&ap->a_ray, ss.inv_dir, rtip->mdl_min, rtip->mdl_max) ||
1671  ap->a_ray.r_max < 0.0) {
1672  cutp = &ap->a_rt_i->rti_inf_box;
1673  if (cutp->bn.bn_len > 0) {
1674  if (n == 0) return cutp;
1675  }
1676  return CUTTER_NULL;
1677  }
1678 
1679  /*
1680  * The interesting part of the ray starts at distance 0. If the
1681  * ray enters the model at a negative distance, (i.e., the ray
1682  * starts within the model RPP), we only look at little bit behind
1683  * (BACKING_DIST) to see if we are just coming out of something,
1684  * but never further back than the intersection with the model
1685  * RPP. If the ray enters the model at a positive distance, we
1686  * always start there. It is vital that we never pick a start
1687  * point outside the model RPP, or the space partitioning tree
1688  * will pick the wrong box and the ray will miss it.
1689  *
1690  * BACKING_DIST should probably be determined by floating point
1691  * noise factor due to model RPP size -vs- number of bits of
1692  * floating point mantissa significance, rather than a constant,
1693  * but that is too hideous to think about here. Also note that
1694  * applications that really depend on knowing what region they are
1695  * leaving from should probably back their own start-point up,
1696  * rather than depending on it here, but it isn't much trouble
1697  * here.
1698  */
1699  ss.box_start = ss.model_start = ap->a_ray.r_min;
1700  ss.box_end = ss.model_end = ap->a_ray.r_max;
1701 
1702  if (ss.box_start < BACKING_DIST)
1703  ss.box_start = BACKING_DIST; /* Only look a little bit behind */
1704 
1705  ss.lastcut = CUTTER_NULL;
1706  ss.old_status = (struct rt_shootray_status *)NULL;
1707  ss.curcut = &ap->a_rt_i->rti_CutHead;
1708  if (ss.curcut->cut_type == CUT_NUGRIDNODE) {
1709  ss.lastcell = CUTTER_NULL;
1710  VSET(ss.curmin, ss.curcut->nugn.nu_axis[X][0].nu_spos,
1711  ss.curcut->nugn.nu_axis[Y][0].nu_spos,
1712  ss.curcut->nugn.nu_axis[Z][0].nu_spos);
1716  } else if (ss.curcut->cut_type == CUT_CUTNODE ||
1717  ss.curcut->cut_type == CUT_BOXNODE) {
1718  ss.lastcell = ss.curcut;
1719  VMOVE(ss.curmin, rtip->mdl_min);
1720  VMOVE(ss.curmax, rtip->mdl_max);
1721  }
1722 
1723  ss.newray = ap->a_ray; /* struct copy */
1724  ss.odist_corr = ss.obox_start = ss.obox_end = -99;
1725  ss.dist_corr = 0.0;
1726 
1727  /*
1728  * While the ray remains inside model space, push from box to box
1729  * until ray emerges from model space again (or first hit is
1730  * found, if user is impatient). It is vitally important to
1731  * always stay within the model RPP, or the space partitioning tree
1732  * will pick wrong boxes & miss them.
1733  */
1734  while ((cutp = rt_advance_to_next_cell(&ss)) != CUTTER_NULL) {
1735  if (debug_shoot) {
1736  rt_pr_cut(cutp, 0);
1737  }
1738  if (--n <= 0) return cutp;
1739 
1740  /* Push ray onwards to next box */
1741  ss.box_start = ss.box_end;
1742  }
1743  return CUTTER_NULL;
1744 }
1745 
1746 
1747 void
1749 {
1750  RT_CK_RESOURCE(resp);
1751 
1752  resp->re_nshootray = 0;
1753  resp->re_nmiss_model = 0;
1754 
1755  resp->re_shots = 0;
1756  resp->re_shot_hit = 0;
1757  resp->re_shot_miss = 0;
1758 
1759  resp->re_prune_solrpp = 0;
1760 
1761  resp->re_ndup = 0;
1762  resp->re_nempty_cells = 0;
1763 
1764  resp->re_piece_shots = 0;
1765  resp->re_piece_shot_hit = 0;
1766  resp->re_piece_shot_miss = 0;
1767  resp->re_piece_ndup = 0;
1768 }
1769 
1770 
1771 void
1772 rt_add_res_stats(register struct rt_i *rtip, register struct resource *resp)
1773 {
1774  RT_CK_RTI(rtip);
1775 
1776  if (!resp) {
1777  resp = &rt_uniresource;
1778  }
1779  RT_CK_RESOURCE(resp);
1780 
1781  rtip->rti_nrays += resp->re_nshootray;
1782  rtip->nmiss_model += resp->re_nmiss_model;
1783 
1784  rtip->nshots += resp->re_shots + resp->re_piece_shots;
1785  rtip->nhits += resp->re_shot_hit + resp->re_piece_shot_hit;
1786  rtip->nmiss += resp->re_shot_miss + resp->re_piece_shot_miss;
1787 
1788  rtip->nmiss_solid += resp->re_prune_solrpp;
1789 
1790  rtip->ndup += resp->re_ndup + resp->re_piece_ndup;
1791  rtip->nempty_cells += resp->re_nempty_cells;
1792 
1793  /* Zero out resource totals, so repeated calls are not harmful */
1794  rt_zero_res_stats(resp);
1795 }
1796 
1797 static int
1798 rt_shootray_simple_hit(struct application *a, struct partition *PartHeadp, struct seg *UNUSED(s))
1799 {
1800  struct partition *p = NULL, *c = NULL, *pp;
1801 
1802  for (pp = PartHeadp->pt_forw; pp != PartHeadp; pp = pp->pt_forw) {
1803  if (p) {
1804  BU_ALLOC(c->pt_forw, struct partition);
1805  c->pt_forw->pt_back = c;
1806  c = c->pt_forw;
1807  } else {
1808  BU_ALLOC(p, struct partition);
1809  c = p;
1810  c->pt_forw = NULL;
1811  c->pt_back = NULL;
1812  }
1813  /* partial deep copy of the partition */
1814  c->pt_magic = pp->pt_magic;
1815  c->pt_inseg = NULL;
1816 
1817  BU_ALLOC(c->pt_inhit, struct hit);
1818  c->pt_inhit->hit_magic = pp->pt_inhit->hit_magic;
1819  c->pt_inhit->hit_dist = pp->pt_inhit->hit_dist;
1820  c->pt_inhit->hit_surfno = pp->pt_inhit->hit_surfno;
1821  c->pt_outseg = NULL;
1822 
1823  BU_ALLOC(c->pt_outhit, struct hit);
1824  c->pt_outhit->hit_magic = pp->pt_outhit->hit_magic;
1825  c->pt_outhit->hit_dist = pp->pt_outhit->hit_dist;
1826  c->pt_outhit->hit_surfno = pp->pt_outhit->hit_surfno;
1827  c->pt_regionp = pp->pt_regionp;
1828  c->pt_inflip = pp->pt_inflip;
1829  c->pt_outflip = pp->pt_outflip;
1830  c->pt_overlap_reg = NULL;
1831  }
1832  a->a_uptr = p;
1833  return 0;
1834 }
1835 
1836 static int
1837 rt_shootray_simple_miss(struct application *a)
1838 {
1839  a->a_uptr = NULL;
1840  return 0;
1841 }
1842 
1843 
1844 struct partition *
1845 rt_shootray_simple(struct application *a, point_t origin, vect_t direction)
1846 {
1847  int (*hit)(struct application *, struct partition *, struct seg *);
1848  int (*miss)(struct application *);
1849  void (*logoverlap)(struct application *, const struct partition *, const struct bu_ptbl *, const struct partition *);
1850 
1851  hit = a->a_hit;
1852  miss = a->a_miss;
1853  logoverlap = a->a_logoverlap;
1854 
1856  a->a_hit = rt_shootray_simple_hit;
1857  a->a_miss = rt_shootray_simple_miss;
1858  VMOVE(a->a_ray.r_pt, origin);
1859  VMOVE(a->a_ray.r_dir, direction);
1860  rt_shootray(a);
1861 
1862  a->a_hit = hit;
1863  a->a_miss = miss;
1864  a->a_logoverlap = logoverlap;
1865 
1866  return (struct partition *)a->a_uptr;
1867 }
1868 
1869 
1870 /*
1871  * Local Variables:
1872  * mode: C
1873  * tab-width: 8
1874  * indent-tabs-mode: t
1875  * c-file-style: "stroustrup"
1876  * End:
1877  * ex: shiftwidth=4 tabstop=8
1878  */
Definition: db_flip.c:35
struct xray a_ray
Actual ray to be shot.
Definition: raytrace.h:1583
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
const union cutter * rt_advance_to_next_cell(register struct rt_shootray_status *ssp)
Definition: shoot.c:170
struct bu_bitv * shot
Definition: raytrace.h:1381
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
void rt_zero_res_stats(struct resource *resp)
Definition: shoot.c:1748
struct seg * a_finished_segs_hdp
Definition: raytrace.h:1612
void rt_res_pieces_clean(struct resource *resp, struct rt_i *rtip)
Definition: shoot.c:72
fastf_t bn_min[3]
Definition: raytrace.h:694
int igrid[3]
integer cell coordinates
Definition: raytrace.h:2266
struct hit seg_in
IN information.
Definition: raytrace.h:370
struct rt_piecestate * re_pieces
array [rti_nsolids_with_pieces]
Definition: raytrace.h:1470
long re_piece_shot_miss
ft_piece_shot() returned a hit
Definition: raytrace.h:1474
#define DEBUG_ALLHITS
2 Print partitions passed to a_hit()
Definition: raytrace.h:83
struct bu_ptbl re_pieces_pending
pieces with an odd hit pending
Definition: raytrace.h:1475
uint32_t pt_magic
sanity check
Definition: raytrace.h:573
Definition: list.h:118
struct bu_bitv * bu_bitv_new(size_t nbits)
Definition: bitv.c:91
point_t mdl_min
min corner of model bounding RPP
Definition: raytrace.h:1769
void(* ft_piece_hitsegs)(struct rt_piecestate *psp, struct seg *seghead, struct application *ap)
Definition: raytrace.h:2074
struct soltab ** bn_list
bn_list[bn_len]
Definition: raytrace.h:696
size_t nhits
solid ft_shot() returned a hit
Definition: raytrace.h:1788
int rt_plot_solid(FILE *fp, struct rt_i *rtip, const struct soltab *stp, struct resource *resp)
#define RT_RAY_MAGIC
Definition: magic.h:165
#define RT_CK_PIECELIST(_p)
Definition: raytrace.h:1410
uint32_t magic
Definition: raytrace.h:1378
#define RT_CK_RTI(_p)
Definition: raytrace.h:1833
const union cutter * cutp
current bounding volume
Definition: raytrace.h:1385
void pdv_3move(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:618
void rt_silent_logoverlap(struct application *ap, const struct partition *pp, const struct bu_ptbl *regiontable, const struct partition *InputHdp)
Definition: bool.c:1122
double dist
>= 0
Definition: tol.h:73
#define BU_PTBL_FOR(ip, cast, ptbl)
Definition: ptbl.h:125
if lu s
Definition: nmg_mod.c:3860
#define PT_DEPARTING_RPP(_step, _lo, _hi, _px, _py, _pz)
Definition: shoot.c:36
void bu_ptbl_init(struct bu_ptbl *b, size_t len, const char *str)
Definition: ptbl.c:32
#define CUT_CUTNODE
Definition: raytrace.h:718
#define VSET(a, b, c, d)
Definition: color.c:53
long re_nempty_cells
number of empty NUgrid cells passed through
Definition: raytrace.h:1468
Definition: raytrace.h:215
#define RT_CK_SEG(_p)
Definition: raytrace.h:377
point_t rti_pmin
for plotting, min RPP
Definition: raytrace.h:1771
#define BU_LIST_IS_EMPTY(hp)
Definition: list.h:295
Definition: raytrace.h:368
#define BU_LIST_IS_INITIALIZED(_hp)
Definition: list.h:175
fastf_t bn_max[3]
Definition: raytrace.h:695
Definition: raytrace.h:248
long re_nmiss_model
Rays pruned by model RPP.
Definition: raytrace.h:1460
void logoverlap(struct application *ap, const struct partition *pp, const struct bu_ptbl *regiontable, const struct partition *InputHdp)
Definition: gqa.c:915
int ft_use_rpp
Definition: raytrace.h:2046
long re_nshootray
Calls to rt_shootray()
Definition: raytrace.h:1459
int bu_ptbl_rm(struct bu_ptbl *b, const long *p)
#define DEBUG_SHOOT
3 Info about rt_shootray() processing
Definition: raytrace.h:84
long ray_seqno
res_nshootray
Definition: raytrace.h:1379
void rt_htbl_init(struct rt_htbl *b, size_t len, const char *str)
Definition: htbl.c:42
fastf_t rt_find_backing_dist(struct rt_shootray_status *ss, struct bu_bitv *backbits)
Definition: shoot.c:651
struct bu_list re_parthead
Head of freelist.
Definition: raytrace.h:1448
Header file for the BRL-CAD common definitions.
fastf_t odist_corr
Definition: raytrace.h:2250
struct xray newray
closer ray start
Definition: raytrace.h:2257
int bu_ptbl_ins_unique(struct bu_ptbl *b, long *p)
#define NUGRID_T_SETUP(_ax, _cval, _cno)
Definition: raytrace.h:2274
#define BU_CK_PTBL(_p)
Definition: ptbl.h:74
struct nu_axis * nu_axis[3]
Definition: raytrace.h:712
#define RT_CK_RAY(_p)
Definition: raytrace.h:224
#define BU_LIST_APPEND(old, new)
Definition: list.h:197
vect_t abs_inv_dir
absolute values of inv_dir
Definition: raytrace.h:2261
struct resource * a_resource
dynamic memory resources
Definition: raytrace.h:1591
struct soltab * stp
Definition: raytrace.h:1380
fastf_t model_start
Definition: raytrace.h:2255
int rt_in_rpp(struct xray *rp, const fastf_t *invdir, const fastf_t *min, const fastf_t *max)
#define BU_LIST_NON_EMPTY(hp)
Definition: list.h:296
void bu_log_indent_delta(int delta)
Definition: log.c:54
int a_onehit
flag to stop on first hit
Definition: raytrace.h:1586
union cutter rti_inf_box
List of infinite solids.
Definition: raytrace.h:1794
fastf_t maxdist
dist ray leaves solids bounding volume
Definition: raytrace.h:1383
void bu_ptbl_reset(struct bu_ptbl *b)
Definition: ptbl.c:49
#define RESOURCE_NULL
Definition: raytrace.h:1489
int(* a_hit)(struct application *, struct partition *, struct seg *)
called when shot hits model
Definition: raytrace.h:1584
fastf_t nu_spos
cell start position
Definition: raytrace.h:705
#define BU_BITSET(_bv, bit)
Definition: bitv.h:183
struct bu_list l
Definition: raytrace.h:369
size_t nempty_cells
number of empty NUgrid cells
Definition: raytrace.h:1792
Definition: ptbl.h:62
int rt_find_nugrid(const struct nugridnode *nugnp, int axis, fastf_t val)
Definition: shoot.c:134
void rt_init_resource(struct resource *resp, int cpu_num, struct rt_i *rtip)
Definition: prep.c:585
long st_bit
solids bit vector index (const)
Definition: raytrace.h:439
#define RT_PIECESTATE_MAGIC
Definition: magic.h:164
struct partition * a_Final_Part_hdp
Definition: raytrace.h:1613
if(share_geom)
Definition: nmg_mod.c:3829
#define RT_AP_CHECK(_ap)
Definition: raytrace.h:1685
long re_prune_solrpp
shot missed solid RPP, ft_shot skipped
Definition: raytrace.h:1466
uint32_t re_magic
Magic number.
Definition: raytrace.h:1441
long re_piece_shot_hit
ft_piece_shot() returned a miss
Definition: raytrace.h:1473
const union cutter * curcut
Definition: raytrace.h:2264
Definition: color.c:49
uint32_t a_magic
Definition: raytrace.h:1581
struct rt_i * a_rt_i
this librt instance
Definition: raytrace.h:1588
#define RT_AP_MAGIC
Definition: magic.h:152
void * memset(void *s, int c, size_t n)
#define DEBUG_PARTITION
14 Info about bool_weave()
Definition: raytrace.h:99
struct resource rt_uniresource
default. Defined in librt/globals.c
Definition: globals.c:41
long re_piece_ndup
ft_piece_shot() calls skipped for already-ft_shot() solids
Definition: raytrace.h:1471
#define RT_G_DEBUG
Definition: raytrace.h:1718
vect_t inv_dir
inverses of ap->a_ray.r_dir
Definition: raytrace.h:2260
#define BU_BITTEST(_bv, bit)
Definition: bitv.h:168
vect_t a_inv_dir
filled in by rt_shootray(), inverse of ray direction cosines
Definition: raytrace.h:1614
size_t bn_len
of solids in list
Definition: raytrace.h:697
struct bu_list l
Definition: bitv.h:106
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
#define CUT_NUGRIDNODE
Definition: raytrace.h:720
size_t rti_nsolids_with_pieces
solids using pieces
Definition: raytrace.h:1819
union cutter rti_CutHead
Head of cut tree.
Definition: raytrace.h:1793
int nu_stepsize[3]
number of cells to jump for one step in each axis
Definition: raytrace.h:714
#define BU_PTBL_GET(ptbl, i)
Definition: ptbl.h:108
struct rt_htbl htab
accumulating hits here
Definition: raytrace.h:1384
int box_num
which cell along ray
Definition: raytrace.h:2270
void pdv_3space(register FILE *plotfp, const fastf_t *min, const fastf_t *max)
Definition: plot3.c:560
int(* ft_shot)(struct soltab *stp, struct xray *rp, struct application *ap, struct seg *seghead)
Definition: raytrace.h:2052
#define RT_CK_AP(_p)
Definition: raytrace.h:1674
#define V3ARGS(a)
Definition: color.c:56
int a_x
Screen X of ray, if applicable.
Definition: raytrace.h:1596
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
static void top()
point_t st_max
max X, Y, Z of bounding RPP
Definition: raytrace.h:438
struct rt_piecelist * bn_piecelist
[] solids with pieces
Definition: raytrace.h:699
#define SQRT_SMALL_FASTF
Definition: defines.h:346
const union cutter * lastcut
Definition: raytrace.h:2263
struct hit seg_out
OUT information.
Definition: raytrace.h:371
#define DEBUG_ADVANCE
24 Cell-to-cell space partitioning
Definition: raytrace.h:109
fastf_t r_max
exit dist from bounding sphere
Definition: raytrace.h:221
struct xray * hit_rayp
pointer to defining ray
Definition: raytrace.h:256
fastf_t dist_corr
correction distance
Definition: raytrace.h:2249
void pdv_3cont(register FILE *plotfp, const fastf_t *pt)
Definition: plot3.c:630
int a_level
recursion level (for printing)
Definition: raytrace.h:1595
struct bu_list l
Definition: ptbl.h:63
const char * a_purpose
Debug string: purpose of ray.
Definition: raytrace.h:1598
#define UNUSED(parameter)
Definition: common.h:239
size_t nmiss
solid ft_shot() returned a miss
Definition: raytrace.h:1787
long re_shots
calls to ft_shot()
Definition: raytrace.h:1462
void rt_pr_partitions(const struct rt_i *rtip, const struct partition *phead, const char *title)
vect_t tv
next t intercept values
Definition: raytrace.h:2267
#define PT_HD_MAGIC
Definition: magic.h:208
#define BU_PTBL_MAGIC
Definition: magic.h:59
goto out
Definition: nmg_mod.c:3846
void rt_pr_seg(const struct seg *segp)
struct partition * pt_back
backwards link
Definition: raytrace.h:575
fastf_t a_ray_length
distance from ray start to end intersections
Definition: raytrace.h:1587
fastf_t nu_epos
cell end position
Definition: raytrace.h:706
struct bu_list re_region_ptbl
head of freelist
Definition: raytrace.h:1453
#define BU_LIST_WHILE(p, structure, hp)
Definition: list.h:410
#define BU_BITV_ZEROALL(_bv)
Definition: bitv.h:191
long st_piecestate_num
re_pieces[] subscript
Definition: raytrace.h:445
void rt_plot_cell(const union cutter *cutp, const struct rt_shootray_status *ssp, struct bu_list *waiting_segs_hd, struct rt_i *rtip)
Definition: shoot.c:764
#define NUGRID_T_ADV(_ax, _cno)
Definition: raytrace.h:2284
fastf_t obox_start
Definition: raytrace.h:2252
size_t bn_piecelen
of piecelists used
Definition: raytrace.h:700
void pl_color(register FILE *plotfp, int r, int g, int b)
Definition: plot3.c:325
vect_t r_dir
Direction of ray (UNIT Length)
Definition: raytrace.h:219
#define BU_PTBL_LEN(ptbl)
Definition: ptbl.h:107
void bu_ptbl_free(struct bu_ptbl *b)
Definition: ptbl.c:226
struct bn_tol rti_tol
Math tolerances for this model.
Definition: raytrace.h:1765
int a_return
Return of a_hit()/a_miss()
Definition: raytrace.h:1601
fastf_t r_min
entry dist to bounding sphere
Definition: raytrace.h:220
void rt_add_res_stats(register struct rt_i *rtip, register struct resource *resp)
Definition: shoot.c:1772
struct cutnode cn
Definition: raytrace.h:725
struct bu_bitv * rt_get_solidbitv(size_t nbits, struct resource *resp)
Definition: prep.c:828
long re_ndup
ft_shot() calls skipped for already-ft_shot() solids
Definition: raytrace.h:1467
size_t end
index of first available location
Definition: raytrace.h:1363
#define BU_LIST_INIT(_hp)
Definition: list.h:148
#define RT_VISIT_ALL_SOLTABS_END
Definition: raytrace.h:1850
void(* a_logoverlap)(struct application *, const struct partition *, const struct bu_ptbl *, const struct partition *)
called to log overlaps
Definition: raytrace.h:1594
long re_shot_hit
ft_shot() returned a miss
Definition: raytrace.h:1463
point_t r_pt
Point at which ray starts.
Definition: raytrace.h:218
struct rt_shootray_status * old_status
Definition: raytrace.h:2269
int rstep[3]
-/0/+ dir of ray in axis
Definition: raytrace.h:2262
point_t st_min
min X, Y, Z of bounding RPP
Definition: raytrace.h:437
struct bu_list re_solid_bitv
head of freelist
Definition: raytrace.h:1452
void rt_boolweave(struct seg *out_hd, struct seg *in_hd, struct partition *PartHeadp, struct application *ap)
Definition: bool.c:145
union cutter * cn_r
val >= point
Definition: raytrace.h:689
#define RT_CK_RESOURCE(_p)
Definition: raytrace.h:1490
#define BU_ASSERT_PTR(_lhs, _relation, _rhs)
Definition: defines.h:227
int hit(struct application *ap, struct partition *PartHeadp, struct seg *segs)
Definition: gqa.c:963
#define DEBUG_ALLRAYS
1 Print calls to rt_shootray()
Definition: raytrace.h:82
#define RT_VISIT_ALL_SOLTABS_START(_s, _rti)
Definition: raytrace.h:1845
axis
Definition: color.c:48
size_t nshots
of calls to ft_shot()
Definition: raytrace.h:1786
const union cutter * rt_cell_n_on_ray(register struct application *ap, int n)
Definition: shoot.c:1547
size_t nsolids
total # of solids participating
Definition: raytrace.h:1783
point_t mdl_max
max corner of model bounding RPP
Definition: raytrace.h:1770
int a_y
Screen Y of ray, if applicable.
Definition: raytrace.h:1597
void rt_htbl_free(struct rt_htbl *b)
Definition: htbl.c:69
#define RT_CK_SOLTAB(_p)
Definition: raytrace.h:453
int rt_shootray(register struct application *ap)
Definition: shoot.c:841
void rt_htbl_reset(struct rt_htbl *b)
Definition: htbl.c:59
int cut_type
Definition: raytrace.h:723
size_t nmiss_solid
shots missed solid RPP
Definition: raytrace.h:1790
long re_shot_miss
ft_shot() returned a hit
Definition: raytrace.h:1464
void bu_bitv_free(struct bu_bitv *bv)
Definition: bitv.c:113
void rt_3cont_raydist(FILE *fp, struct xray *rayp, double dist)
Definition: shoot.c:754
union cutter * cn_l
val < point
Definition: raytrace.h:688
int bu_list_len(const struct bu_list *hd)
int out_axis
axis ray will leave through
Definition: raytrace.h:2268
#define BU_ASSERT_DOUBLE(_lhs, _relation, _rhs)
Definition: defines.h:279
int miss(struct application *ap)
Definition: gqa.c:1278
int(* ft_piece_shot)(struct rt_piecestate *psp, struct rt_piecelist *plp, double dist, struct xray *ray, struct application *ap, struct seg *seghead)
Definition: raytrace.h:2066
Definition: color.c:51
void rt_prep_parallel(struct rt_i *rtip, int ncpu)
void * a_uptr
application-specific pointer
Definition: raytrace.h:1618
struct partition * rt_shootray_simple(struct application *a, point_t origin, vect_t direction)
Definition: shoot.c:1845
struct resource * resp
Definition: raytrace.h:2259
#define RT_CK_PIECESTATE(_p)
Definition: raytrace.h:1387
struct boxnode bn
Definition: raytrace.h:726
size_t ndup
duplicate shots at a given solid
Definition: raytrace.h:1791
int re_cpu
processor number, for ID
Definition: raytrace.h:1442
struct nugridnode nugn
Definition: raytrace.h:727
const union cutter * lastcell
Definition: raytrace.h:2263
#define CUTTER_NULL
Definition: raytrace.h:731
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
long re_piece_shots
calls to ft_piece_shot()
Definition: raytrace.h:1472
int(* a_miss)(struct application *)
called when shot misses
Definition: raytrace.h:1585
int cn_axis
0, 1, 2 = cut along X, Y, Z
Definition: raytrace.h:686
fastf_t mindist
dist ray enters solids bounding volume
Definition: raytrace.h:1382
struct bu_ptbl rti_resources
list of 'struct resource's encountered
Definition: raytrace.h:1811
void rt_pr_cut(const union cutter *cutp, int lvl)
Definition: cut.c:1703
int rt_boolfinal(struct partition *InputHdp, struct partition *FinalHdp, fastf_t startdist, fastf_t enddist, struct bu_ptbl *regionbits, struct application *ap, const struct bu_bitv *solidbits)
Definition: bool.c:1490
const struct rt_functab * st_meth
pointer to per-solid methods
Definition: raytrace.h:428
struct soltab * stp
ref back to solid
Definition: raytrace.h:1408
int nu_cells_per_axis[3]
number of slabs
Definition: raytrace.h:713
void pdv_3box(register FILE *plotfp, const fastf_t *a, const fastf_t *b)
Definition: plot3.c:688
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
#define BU_LIST_MAGIC_EQUAL(_l, _magic)
Definition: list.h:135
#define CUT_BOXNODE
Definition: raytrace.h:719
Definition: bitv.h:105
struct partition * pt_forw
forwards link
Definition: raytrace.h:574
#define BU_PTBL_TEST(ptbl)
Definition: ptbl.h:110
fastf_t hit_dist
dist from r_pt to hit_point
Definition: raytrace.h:250
#define RT_FREE_PT_LIST(_headp, _res)
Definition: raytrace.h:630
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300
#define BU_CK_BITV(_bp)
Definition: bitv.h:116
#define VPRINT(a, b)
Definition: raytrace.h:1881
union cutter * nu_grid
3-D array of boxnodes
Definition: raytrace.h:715
struct application * ap
Definition: raytrace.h:2258
point_t rti_pmax
for plotting, max RPP
Definition: raytrace.h:1772
int needprep
needs rt_prep
Definition: raytrace.h:1776
Definition: color.c:50
size_t rti_nrays
calls to rt_shootray()
Definition: raytrace.h:1784
uint32_t magic
Definition: raytrace.h:216
size_t nmiss_model
rays missed model RPP
Definition: raytrace.h:1785
#define BU_LIST_FIRST(structure, hp)
Definition: list.h:312
void rt_3move_raydist(FILE *fp, struct xray *rayp, double dist)
Definition: shoot.c:744
void rt_res_pieces_init(struct resource *resp, struct rt_i *rtip)
Definition: shoot.c:46
fastf_t cn_point
cut through axis==point
Definition: raytrace.h:687
#define RT_FREE_SEG_LIST(_segheadp, _res)
Definition: raytrace.h:401
long st_npieces
pieces used by this solid
Definition: raytrace.h:444
#define BACKING_DIST
mm to look behind start point
Definition: raytrace.h:2290
#define UNLIKELY(expression)
Definition: common.h:282