00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
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;
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);
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
00089
00090
00091
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
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;
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
00153
00154
00155
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
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
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
00224
00225
00226
00227 cutp = CUTTER_NULL;
00228 push_flag = 0;
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 t0 = ssp->box_start;
00251
00252
00253
00254 top: switch( curcut->cut_type ) {
00255 case CUT_NUGRIDNODE: {
00256
00257
00258
00259
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
00275
00276
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
00284
00285
00286
00287
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
00308
00309 cutp = ssp->lastcell;
00310 out_axis = ssp->out_axis;
00311
00312
00313
00314
00315
00316
00317
00318
00319
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
00328
00329
00330
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
00349 NUGRID_T_ADV( out_axis, ssp->igrid[out_axis] );
00350 }
00351
00352
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
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
00392 case CUT_BOXNODE:
00393
00394
00395
00396
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
00412
00413
00414
00415 {
00416 register struct rt_shootray_status *old = ssp->old_status;
00417
00418 if( old == NULL ) goto escaped_from_model;
00419 *ssp = *old;
00420 bu_free( old, "old rt_shootray_status" );
00421 curcut = ssp->curcut;
00422 cutp = CUTTER_NULL;
00423 continue;
00424 }
00425 }
00426
00427
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
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
00444
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 , 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
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
00502
00503
00504 t0 += OFFSET_DIST;
00505 goto top;
00506 }
00507
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 ,
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
00532
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
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
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 ,
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
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;
00631
00632
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
00652 }
00653
00654
00655
00656
00657
00658
00659
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
00670
00671
00672
00673
00674
00675
00676
00677
00678
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
00696 solidbits = get_solidbitv( rtip->nsolids, resp );
00697 bu_bitv_clear(solidbits);
00698
00699 ray = ss->ap->a_ray;
00700
00701
00702
00703 while( cur_dist <= ss->ap->a_ray.r_max ) {
00704
00705 VJOIN1( cur_pt, ss->ap->a_ray.r_pt, cur_dist, ss->ap->a_ray.r_dir );
00706
00707
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
00718
00719 if( !rt_in_rpp( &ray, ss->inv_dir, cutp->bn.bn_min, cutp->bn.bn_max ) ) {
00720
00721
00722
00723
00724
00725
00726 goto done;
00727 } else {
00728
00729 cur_dist = ray.r_max + OFFSET_DIST;
00730 }
00731
00732
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
00738 if( rt_in_rpp( &ray, ss->inv_dir, plp->stp->st_min, plp->stp->st_max ) ) {
00739
00740
00741 if( ray.r_min < BACKING_DIST ) {
00742 if( ray.r_min < min_backing_dist ) {
00743
00744 min_backing_dist = ray.r_min;
00745 }
00746
00747
00748 BU_BITSET( backbits, plp->stp->st_bit );
00749 }
00750 }
00751
00752 BU_BITSET( solidbits, plp->stp->st_bit );
00753 }
00754 }
00755 }
00756 done:
00757
00758 BU_CK_BITV(solidbits);
00759 BU_LIST_APPEND( &resp->re_solid_bitv, &solidbits->l );
00760
00761
00762 return min_backing_dist;
00763 }
00764
00765
00766
00767
00768
00769
00770
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793 int
00794 rt_shootray(register struct application *ap)
00795 {
00796 struct rt_shootray_status ss;
00797 struct seg new_segs;
00798 struct seg waiting_segs;
00799 struct seg finished_segs;
00800 fastf_t last_bool_start;
00801 struct bu_bitv *solidbits;
00802 struct bu_bitv *backbits=NULL;
00803
00804 struct bu_ptbl *regionbits;
00805 char *status;
00806 auto struct partition InitialPart;
00807 auto struct partition FinalPart;
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;
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
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);
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
00883 bu_log("rt_shootray() resp=x%x uninitialized, fixing it\n", resp);
00884
00885
00886
00887
00888
00889
00890
00891 rt_init_resource( resp, resp->re_cpu, rtip );
00892 }
00893
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( ®ionbits->l );
00906 BU_CK_PTBL(regionbits);
00907 }
00908
00909 if( !resp->re_pieces && rtip->rti_nsolids_with_pieces > 0 ) {
00910
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
00918 if(RT_G_DEBUG) {
00919 FAST fastf_t f, diff;
00920
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
00937
00938 resp->re_nshootray++;
00939
00940
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
00978
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
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;
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
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042
01043
01044
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
01050 if( ss.box_start < BACKING_DIST ) {
01051
01052
01053
01054
01055
01056
01057
01058
01059 backbits = get_solidbitv( rtip->nsolids, resp );
01060 bu_bitv_clear(backbits);
01061
01062
01063
01064
01065 ss.box_start = rt_find_backing_dist( &ss, backbits );
01066 }
01067 } else {
01068
01069 if( ss.box_start < BACKING_DIST )
01070 ss.box_start = BACKING_DIST;
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;
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
01099
01100
01101
01102
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
01113 ss.box_start = ss.box_end;
01114 resp->re_nempty_cells++;
01115 continue;
01116 }
01117
01118
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
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
01138
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
01147 BU_BITV_ZEROALL(psp->shot);
01148 psp->ray_seqno = resp->re_nshootray;
01149 rt_htbl_reset( &psp->htab );
01150
01151
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;
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
01166 resp->re_ndup++;
01167 continue;
01168 }
01169 had_hits_before = psp->htab.end;
01170 }
01171
01172
01173
01174
01175
01176
01177
01178
01179
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
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
01193
01194
01195
01196 if( ss.box_end > psp->maxdist && psp->htab.end > 0 ) {
01197
01198 if(debug_shoot)bu_log("shooting %s pieces complete, making segs\n", stp->st_name);
01199
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
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;
01222 }
01223
01224
01225 BU_BITSET( solidbits, stp->st_bit );
01226
01227
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;
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;
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;
01249 }
01250
01251
01252 {
01253 register struct seg *s2;
01254 while(BU_LIST_WHILE(s2,seg,&(new_segs.l))) {
01255 BU_LIST_DEQUEUE( &(s2->l) );
01256
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
01271
01272
01273
01274
01275
01276
01277
01278
01279
01280
01281
01282
01283
01284
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
01298 rt_boolweave( &finished_segs, &waiting_segs, &InitialPart, ap );
01299
01300 if( BU_PTBL_LEN( &resp->re_pieces_pending ) > 0 ) {
01301
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
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
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
01332 ss.box_start = ss.box_end;
01333 }
01334
01335
01336
01337
01338
01339 weave:
01340 if( RT_G_DEBUG&DEBUG_ADVANCE )
01341 bu_log( "rt_shootray: ray has left known space\n" );
01342
01343
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
01349
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
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
01371
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
01387
01388
01389
01390
01391
01392
01393
01394
01395
01396
01397 hitit:
01398 if(debug_shoot) rt_pr_partitions(rtip,&FinalPart,"a_hit()");
01399
01400
01401
01402
01403
01404
01405 RT_FREE_PT_LIST( &InitialPart, resp );
01406
01407
01408
01409
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
01420
01421 out:
01422
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, ®ionbits->l );
01431
01432
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
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
01457
01458
01459
01460
01461
01462
01463 const union cutter *
01464 rt_cell_n_on_ray(register struct application *ap, int n)
01465
01466
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);
01525
01526 if( BU_LIST_UNINITIALIZED( &resp->re_parthead ) ) {
01527
01528 bu_log("rt_cell_n_on_ray() resp=x%x uninitialized, fixing it\n", resp);
01529
01530
01531
01532
01533
01534
01535
01536 rt_init_resource( resp, resp->re_cpu, rtip );
01537 }
01538
01539 BU_ASSERT_PTR( BU_PTBL_GET(&rtip->rti_resources, resp->re_cpu), !=, NULL );
01540
01541
01542 if(RT_G_DEBUG) {
01543 FAST fastf_t f, diff;
01544
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
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
01596
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
01609
01610
01611
01612
01613
01614
01615
01616
01617
01618
01619
01620
01621
01622
01623
01624
01625
01626
01627
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;
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;
01654 ss.odist_corr = ss.obox_start = ss.obox_end = -99;
01655 ss.dist_corr = 0.0;
01656
01657
01658
01659
01660
01661
01662
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
01671 ss.box_start = ss.box_end;
01672 }
01673 return CUTTER_NULL;
01674 }
01675
01676
01677
01678
01679
01680
01681
01682
01683
01684
01685
01686
01687
01688
01689
01690
01691
01692
01693
01694
01695
01696
01697
01698
01699
01700
01701
01702 int
01703 rt_in_rpp(struct xray *rp,
01704 register const fastf_t *invdir,
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
01711 register fastf_t rmin = -INFINITY;
01712 register fastf_t rmax = INFINITY;
01713
01714
01715
01716
01717 if( *invdir < 0.0 ) {
01718
01719
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
01726
01727 if(rmax > (st = (*max - *pt) * *invdir) )
01728 rmax = st;
01729 if( rmin < ((sv = (*min - *pt) * *invdir)) )
01730 rmin = sv;
01731 } else {
01732
01733
01734
01735
01736
01737 if( (*min > *pt) || (*max < *pt) )
01738 return(0);
01739 }
01740
01741
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);
01756 }
01757
01758
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);
01773 }
01774
01775
01776 if( rmin > rmax )
01777 return(0);
01778
01779
01780 rp->r_min = rmin;
01781 rp->r_max = rmax;
01782 return(1);
01783 }
01784
01785
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
01790
01791
01792 {
01793 register const fastf_t *pt = &rp->r_pt[0];
01794 FAST fastf_t sv;
01795
01796
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
01805 if( rp->r_dir[X] < 0.0 ) {
01806
01807
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
01819
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
01832
01833
01834
01835 if( (*min > *pt) || (*max < *pt) )
01836 goto miss;
01837 }
01838
01839
01840 pt++; invdir++; max++; min++;
01841 if( rp->r_dir[Y] < 0.0 ) {
01842
01843
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
01855
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
01871 pt++; invdir++; max++; min++;
01872 if( rp->r_dir[Z] < 0.0 ) {
01873
01874
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
01886
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
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);
01906 miss:
01907 bu_log("MISS\n");
01908 return(0);
01909 }
01910
01911 #undef st
01912
01913
01914 #define SEG_MISS(SEG) (SEG).seg_stp=(struct soltab *) 0;
01915
01916
01917 void
01918 rt_vstub(struct soltab **stp, struct xray **rp, struct seg *segp, int n, struct application *ap)
01919
01920
01921
01922
01923
01924 {
01925 register int i;
01926 register struct seg *tmp_seg;
01927 struct seg seghead;
01928
01929 BU_LIST_INIT( &(seghead.l) );
01930
01931
01932 for (i = 0; i < n; i++) {
01933 if (stp[i] != 0){
01934
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;
01941 RT_FREE_SEG(tmp_seg, ap->a_resource);
01942 }
01943 }
01944 }
01945 }
01946
01947
01948
01949
01950
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
01983
01984
01985
01986
01987
01988
01989
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
02012 rt_zero_res_stats( resp );
02013 }
02014
02015
02016
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
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 );
02059
02060
02061 pdv_3space( fp, rtip->rti_pmin, rtip->rti_pmax );
02062 pdv_3box( fp, rtip->rti_pmin, rtip->rti_pmax );
02063
02064
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
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
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
02092 pl_color( fp, 255, 255, 0 );
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
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 );
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
02116
02117
02118
02119
02120
02121
02122