nmg_rt_isect.c

Go to the documentation of this file.
00001 /*                  N M G _ R T _ I S E C T . C
00002  * BRL-CAD
00003  *
00004  * Copyright (c) 1994-2006 United States Government as represented by
00005  * the U.S. Army Research Laboratory.
00006  *
00007  * This library is free software; you can redistribute it and/or
00008  * modify it under the terms of the GNU Lesser General Public License
00009  * as published by the Free Software Foundation; either version 2 of
00010  * the License, or (at your option) any later version.
00011  *
00012  * This library is distributed in the hope that it will be useful, but
00013  * WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00015  * Library General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU Lesser General Public
00018  * License along with this file; see the file named COPYING for more
00019  * information.
00020  */
00021 
00022 /** @addtogroup nmg */
00023 
00024 /*@{*/
00025 /** @file nmg_rt_isect.c
00026  *      Support routines for raytracing an NMG.
00027  *
00028  *  Author -
00029  *      Lee A. Butler
00030  *
00031  *  Source -
00032  *      The U. S. Army Research Laboratory
00033  *      Aberdeen Proving Ground, Maryland  21005-5068  USA
00034  */
00035 /*@}*/
00036 
00037 #ifndef lint
00038 static const char RCSid[] = "@(#)$Header: /cvsroot/brlcad/brlcad/src/librt/nmg_rt_isect.c,v 14.15 2006/09/16 02:04:25 lbutler Exp $ (ARL)";
00039 #endif
00040 
00041 #include "common.h"
00042 
00043 #include <stdlib.h>
00044 #include <stdio.h>
00045 #ifdef HAVE_STRING_H
00046 #  include <string.h>
00047 #endif
00048 #include <math.h>
00049 
00050 #include "machine.h"
00051 #include "vmath.h"
00052 #include "nmg.h"
00053 #include "raytrace.h"
00054 #include "nurb.h"
00055 #include "plot3.h"
00056 
00057 
00058 static void     vertex_neighborhood BU_ARGS((struct ray_data *rd, struct vertexuse *vu_p, struct hitmiss *myhit));
00059 
00060 const char *
00061 nmg_rt_inout_str(int code)
00062 {
00063         switch(code) {
00064         case HMG_HIT_IN_IN:
00065                 return("IN_IN");
00066         case HMG_HIT_IN_OUT:
00067                 return("IN_OUT");
00068         case HMG_HIT_OUT_IN:
00069                 return("OUT_IN");
00070         case HMG_HIT_OUT_OUT:
00071                 return("OUT_OUT");
00072         case HMG_HIT_IN_ON:
00073                 return("IN_ON");
00074         case HMG_HIT_ON_IN:
00075                 return("ON_IN");
00076         case HMG_HIT_OUT_ON:
00077                 return("OUT_ON");
00078         case HMG_HIT_ON_OUT:
00079                 return("ON_OUT");
00080         case HMG_HIT_ANY_ANY:
00081                 return("ANY_ANY");
00082         }
00083         return("?_?\n");
00084 }
00085 
00086 const char *
00087 nmg_rt_state_str(int code)
00088 {
00089         switch(code) {
00090         case NMG_RAY_STATE_INSIDE:
00091                 return "RS_INSIDE";
00092         case NMG_RAY_STATE_ON:
00093                 return "RS_ON";
00094         case NMG_RAY_STATE_OUTSIDE:
00095                 return "RS_OUTSIDE";
00096         case NMG_RAY_STATE_ANY:
00097                 return "RS_ANY";
00098         }
00099         return "RS_UNKNOWN";
00100 }
00101 
00102 /**
00103  *                              N M G _ C K _ H I T M I S S _ L I S T
00104  *
00105  *  Ensure that the ray makes a valid set of state transitions.
00106  */
00107 void
00108 nmg_ck_hitmiss_list(const struct bu_list *hd)
00109 {
00110         struct hitmiss  *hmp;
00111         int             state = NMG_RAY_STATE_OUTSIDE;
00112         int             istate;
00113         int             ostate;
00114 
00115         for( BU_LIST_FOR( hmp, hitmiss, hd ) )  {
00116 #ifndef FAST_NMG
00117                 NMG_CK_HITMISS(hmp);
00118 #endif
00119                 /* Skip hits on non-3-manifolds */
00120                 if( hmp->in_out == HMG_HIT_ANY_ANY )  continue;
00121 
00122                 istate = HMG_INBOUND_STATE(hmp);
00123                 ostate = HMG_OUTBOUND_STATE(hmp);
00124                 if( state != istate )  {
00125                         bu_log("ray state was=%s, transition=%s (istate=%s)\n",
00126                                 nmg_rt_state_str(state),
00127                                 nmg_rt_inout_str(hmp->in_out),
00128                                 nmg_rt_state_str(istate) );
00129                         rt_bomb("nmg_ck_hitmiss_list() NMG ray-tracer bad in/out state transition\n");
00130                 }
00131                 state = ostate;
00132         }
00133         if( state != NMG_RAY_STATE_OUTSIDE )  {
00134                 bu_log("ray ending state was %s, should have been RS_OUT\n", nmg_rt_state_str(state));
00135                 rt_bomb("nmg_ck_hitmiss_list() NMG ray-tracer bad ending state\n");
00136         }
00137 }
00138 
00139 /* Plot a faceuse and a line between pt and plane_pt */
00140 static int plot_file_number=0;
00141 
00142 static void
00143 nmg_rt_isect_plfu(struct faceuse *fu, fastf_t *pt, fastf_t *plane_pt)
00144 {
00145         FILE *fd;
00146         char name[25];
00147         long *b;
00148 
00149         NMG_CK_FACEUSE(fu);
00150 
00151         sprintf(name, "ray%02d.pl", plot_file_number++);
00152         if ((fd=fopen(name, "w")) == (FILE *)NULL) {
00153                 perror(name);
00154                 bu_bomb("unable to open file for writing");
00155         }
00156 
00157         bu_log("overlay %s\n", name);
00158         b = (long *)bu_calloc( fu->s_p->r_p->m_p->maxindex,
00159                 sizeof(long), "bit vec"),
00160 
00161         pl_erase(fd);
00162         pd_3space(fd,
00163                 fu->f_p->min_pt[0]-1.0,
00164                 fu->f_p->min_pt[1]-1.0,
00165                 fu->f_p->min_pt[2]-1.0,
00166                 fu->f_p->max_pt[0]+1.0,
00167                 fu->f_p->max_pt[1]+1.0,
00168                 fu->f_p->max_pt[2]+1.0);
00169 
00170         nmg_pl_fu(fd, fu, b, 255, 255, 255);
00171 
00172         pl_color(fd, 255, 50, 50);
00173         pdv_3line(fd, pt, plane_pt);
00174 
00175         bu_free((char *)b, "bit vec");
00176         fclose(fd);
00177 }
00178 
00179 static void
00180 pleu(struct edgeuse *eu, fastf_t *pt, fastf_t *plane_pt)
00181 {
00182         FILE *fd;
00183         char name[25];
00184         long *b;
00185         point_t min_pt, max_pt;
00186         int i;
00187         struct model *m;
00188 
00189         sprintf(name, "ray%02d.pl", plot_file_number++);
00190         if ((fd=fopen(name, "w")) == (FILE *)NULL) {
00191                 perror(name);
00192                 bu_bomb("unable to open file for writing");
00193         }
00194 
00195         bu_log("overlay %s\n", name);
00196         m = nmg_find_model( eu->up.magic_p );
00197         b = (long *)bu_calloc( m->maxindex, sizeof(long), "bit vec");
00198 
00199         pl_erase(fd);
00200 
00201         VMOVE(min_pt, eu->vu_p->v_p->vg_p->coord);
00202 
00203         for (i=0 ; i < 3 ; i++) {
00204                 if (eu->eumate_p->vu_p->v_p->vg_p->coord[i] < min_pt[i]) {
00205                         max_pt[i] = min_pt[i];
00206                         min_pt[i] = eu->eumate_p->vu_p->v_p->vg_p->coord[i];
00207                 } else {
00208                         max_pt[i] = eu->eumate_p->vu_p->v_p->vg_p->coord[i];
00209                 }
00210         }
00211         pd_3space(fd,
00212                 min_pt[0]-1.0, min_pt[1]-1.0, min_pt[2]-1.0,
00213                 max_pt[0]+1.0, max_pt[1]+1.0, max_pt[2]+1.0);
00214 
00215         nmg_pl_eu(fd, eu, b, 255, 255, 255);
00216         pl_color(fd, 255, 50, 50);
00217         pdv_3line(fd, pt, plane_pt);
00218         bu_free((char *)b, "bit vec");
00219         fclose(fd);
00220 }
00221 static void
00222 plvu(struct vertexuse *vu)
00223 {
00224 }
00225 
00226 void
00227 nmg_rt_print_hitmiss(struct hitmiss *a_hit)
00228 {
00229         NMG_CK_HITMISS(a_hit);
00230         bu_log("   dist:%12g pt=(%f %f %f) %s=x%x\n",
00231                 a_hit->hit.hit_dist,
00232                 a_hit->hit.hit_point[0],
00233                 a_hit->hit.hit_point[1],
00234                 a_hit->hit.hit_point[2],
00235                 bu_identify_magic(*(long *)a_hit->hit.hit_private),
00236                 a_hit->hit.hit_private
00237         );
00238         bu_log("\tstate:%s", nmg_rt_inout_str(a_hit->in_out));
00239 
00240         switch (a_hit->start_stop) {
00241         case NMG_HITMISS_SEG_IN:        bu_log(" SEG_START"); break;
00242         case NMG_HITMISS_SEG_OUT:       bu_log(" SEG_STOP"); break;
00243         }
00244 
00245         VPRINT("\n\tin_normal", a_hit->inbound_norm);
00246         VPRINT("\tout_normal", a_hit->outbound_norm);
00247 }
00248 void
00249 nmg_rt_print_hitlist(struct hitmiss *hl)
00250 {
00251         struct hitmiss *a_hit;
00252 
00253         bu_log("nmg/ray hit list:\n");
00254 
00255         for (BU_LIST_FOR(a_hit, hitmiss, &(hl->l))) {
00256                 nmg_rt_print_hitmiss(a_hit);
00257         }
00258 }
00259 
00260 
00261 
00262 
00263 /**
00264  *      H I T _ I N S
00265  *
00266  *      insert the new hit point in the correct place in the list of hits
00267  *      so that the list is always in sorted order
00268  */
00269 static void
00270 hit_ins(struct ray_data *rd, struct hitmiss *newhit)
00271 {
00272         struct hitmiss *a_hit;
00273 
00274         BU_LIST_MAGIC_SET(&newhit->l, NMG_RT_HIT_MAGIC);
00275         newhit->hit.hit_magic = RT_HIT_MAGIC;
00276 
00277         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00278                 bu_log("hit_ins()\n  inserting:");
00279                 nmg_rt_print_hitmiss(newhit);
00280         }
00281 
00282         for (BU_LIST_FOR(a_hit, hitmiss, &rd->rd_hit)) {
00283                 if (newhit->hit.hit_dist < a_hit->hit.hit_dist)
00284                         break;
00285         }
00286 
00287         /* a_hit now points to the item before which we should insert this
00288          * hit in the hit list.
00289          */
00290 
00291         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00292                 if (BU_LIST_NOT_HEAD(a_hit, &rd->rd_hit)) {
00293                         bu_log("   prior to:");
00294                         nmg_rt_print_hitmiss(a_hit);
00295                 } else {
00296                         bu_log("\tat end of list\n");
00297                 }
00298         }
00299 
00300         BU_LIST_INSERT(&a_hit->l, &newhit->l);
00301 
00302         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00303                 nmg_rt_print_hitlist((struct hitmiss *)&rd->rd_hit);
00304 
00305 }
00306 
00307 
00308 /*
00309  *  The ray missed this vertex.  Build the appropriate miss structure.
00310  */
00311 static struct hitmiss *
00312 ray_miss_vertex(struct ray_data *rd, struct vertexuse *vu_p)
00313 {
00314         struct hitmiss *myhit;
00315 
00316         NMG_CK_VERTEXUSE(vu_p);
00317         NMG_CK_VERTEX(vu_p->v_p);
00318         NMG_CK_VERTEX_G(vu_p->v_p->vg_p);
00319 
00320         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00321                 bu_log("ray_miss_vertex(%g %g %g)\n",
00322                         vu_p->v_p->vg_p->coord[0],
00323                         vu_p->v_p->vg_p->coord[1],
00324                         vu_p->v_p->vg_p->coord[2]);
00325         }
00326 
00327         if ( (myhit=NMG_INDEX_GET(rd->hitmiss, vu_p->v_p)) ) {
00328                 /* vertex previously processed */
00329                 if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_MAGIC)) {
00330                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00331                                 bu_log("ray_miss_vertex( vertex previously HIT!!!! )\n");
00332                 } else if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_SUB_MAGIC)) {
00333                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00334                                 bu_log("ray_miss_vertex( vertex previously HIT_SUB!?!? )\n");
00335                 }
00336                 return(myhit);
00337         }
00338         if ( (myhit=NMG_INDEX_GET(rd->hitmiss, vu_p->v_p)) ) {
00339                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00340                         bu_log("ray_miss_vertex( vertex previously missed )\n");
00341                 return(myhit);
00342         }
00343 
00344         GET_HITMISS(myhit, rd->ap);
00345         NMG_INDEX_ASSIGN(rd->hitmiss, vu_p->v_p, myhit);
00346         myhit->outbound_use = (long *)vu_p;
00347         myhit->inbound_use = (long *)vu_p;
00348         myhit->hit.hit_private = (genptr_t)vu_p->v_p;
00349 
00350         /* get build_vertex_miss() to compute this */
00351         myhit->dist_in_plane = -1.0;
00352 
00353         /* add myhit to the list of misses */
00354         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
00355         BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
00356 #ifndef FAST_NMG
00357         NMG_CK_HITMISS(myhit);
00358         NMG_CK_HITMISS_LISTS(myhit, rd);
00359 #endif
00360         return myhit;
00361 }
00362 
00363 /*
00364  * Support routine for vertex_neighborhood()
00365  *
00366  */
00367 static void
00368 get_pole_dist_to_face(struct ray_data *rd, struct vertexuse *vu, fastf_t *Pole, fastf_t *Pole_prj_pt, double *Pole_dist, fastf_t *Pole_pca, fastf_t *pointA, fastf_t *leftA, fastf_t *pointB, fastf_t *leftB, fastf_t *polar_height_vect, char *Pole_name)
00369 {
00370         vect_t pca_to_pole_vect;
00371         vect_t VtoPole_prj;
00372         point_t pcaA, pcaB;
00373         double distA, distB;
00374         int code, status;
00375 
00376 /* find the points of closest approach
00377  *  There are six distinct return values from bn_dist_pt3_lseg3():
00378  *
00379  *    Value     Condition
00380  *      -----------------------------------------------------------------
00381  *      0       P is within tolerance of lseg AB.  *dist isn't 0: (SPECIAL!!!)
00382  *                *dist = parametric dist = |PCA-A| / |B-A|.  pca=computed.
00383  *      1       P is within tolerance of point A.  *dist = 0, pca=A.
00384  *      2       P is within tolerance of point B.  *dist = 0, pca=B.
00385  *
00386  *      3       P is to the "left" of point A.  *dist=|P-A|, pca=A.
00387  *      4       P is to the "right" of point B.  *dist=|P-B|, pca=B.
00388  *      5       P is "above/below" lseg AB.  *dist=|PCA-P|, pca=computed.
00389  */
00390         code = bn_dist_pt3_lseg3( &distA, pcaA, vu->v_p->vg_p->coord, pointA,
00391                         Pole_prj_pt, rd->tol );
00392         if (code < 3) {
00393                 /* Point is on line */
00394                 *Pole_dist = MAGNITUDE(polar_height_vect);
00395                 VMOVE(Pole_pca, Pole_prj_pt);
00396                 return;
00397         }
00398 
00399         status = code << 4;
00400         code = bn_dist_pt3_lseg3( &distB, pcaB, vu->v_p->vg_p->coord, pointB,
00401                         Pole_prj_pt, rd->tol );
00402         if (code < 3) {
00403                 /* Point is on line */
00404                 *Pole_dist = MAGNITUDE(polar_height_vect);
00405                 VMOVE(Pole_pca, Pole_prj_pt);
00406                 return;
00407         }
00408 
00409         status |= code;
00410 
00411         /*        Status
00412          *      codeA CodeB
00413          *       3      3       Do the Tanenbaum patch thing
00414          *       3      4       This should not happen.
00415          *       3      5       compute dist from pole to pcaB
00416          *       4      3       This should not happen.
00417          *       4      4       This should not happen.
00418          *       4      5       This should not happen.
00419          *       5      3       compute dist from pole to pcaA
00420          *       5      4       This should not happen.
00421          *       5      5       pick the edge with smallest "dist"
00422          *                      and compute pole-pca distance.
00423          */
00424         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00425                 bu_log("get_pole_dist_to_face(%s) status from dist_pt_lseg == 0x%02x\n", Pole_name, status);
00426 
00427         switch (status) {
00428         case 0x35: /* compute dist from pole to pcaB */
00429 
00430                 /* if plane point is "inside" edge B, plane point is PCA */
00431                 VSUB2(VtoPole_prj, Pole_prj_pt, vu->v_p->vg_p->coord);
00432                 if (VDOT(leftB, VtoPole_prj) >= 0.0) {
00433                         /* plane point is "inside" edge B */
00434                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00435                                 bu_log("\tplane point inside face\n");
00436                         VMOVE(Pole_pca, Pole_prj_pt);
00437                         VSUB2(pca_to_pole_vect, Pole_prj_pt, Pole);
00438                         *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00439                         return;
00440                 }
00441 
00442                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00443                                 bu_log("\tplane point outside face\n");
00444                 VSUB2(pca_to_pole_vect, pcaB, Pole);
00445                 VMOVE(Pole_pca, pcaB);
00446                 break;
00447         case 0x53: /* compute dist from pole to pcaA */
00448 
00449                 /* if plane point is "inside" edge A, plane point is PCA */
00450                 VSUB2(VtoPole_prj, Pole_prj_pt, vu->v_p->vg_p->coord);
00451                 if (VDOT(leftA, VtoPole_prj) >= 0.0) {
00452                         /* plane point is "inside" edge A */
00453                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00454                                 bu_log("\tplane point inside face\n");
00455                         VMOVE(Pole_pca, Pole_prj_pt);
00456                         VSUB2(pca_to_pole_vect, Pole_prj_pt, Pole);
00457                         *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00458                         return;
00459                 }
00460 
00461                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00462                                 bu_log("\tplane point outside face\n");
00463                 VSUB2(pca_to_pole_vect, pcaA, Pole);
00464                 VMOVE(Pole_pca, pcaA);
00465                 break;
00466         case 0x55:/* pick the edge with smallest "dist"
00467                    * and compute pole-pca distance.
00468                    */
00469                 VSUB2(VtoPole_prj, Pole_prj_pt, vu->v_p->vg_p->coord);
00470 
00471                 if (distA < distB) {
00472                         VUNITIZE(VtoPole_prj);
00473                         VUNITIZE(leftA);
00474                         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00475                                 VPRINT("LeftA", leftA);
00476                                 VPRINT("VtoPole_prj", VtoPole_prj);
00477                         }
00478 
00479                         if (VDOT(leftA, VtoPole_prj) >= 0.0) {
00480                                 /* plane point is "inside" edge A */
00481                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00482                                         bu_log("\tplane point inside face\n");
00483                                 VSUB2(pca_to_pole_vect, Pole_prj_pt, Pole);
00484                                 *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00485                                 VMOVE(Pole_pca, Pole_prj_pt);
00486                                 return;
00487                         }
00488                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00489                                         bu_log("\tplane point outside face\n");
00490                         VSUB2(pca_to_pole_vect, pcaA, Pole);
00491                         VMOVE(Pole_pca, pcaA);
00492                 } else {
00493                         VUNITIZE(VtoPole_prj);
00494                         VUNITIZE(leftB);
00495                         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00496                                 VPRINT("LeftB", leftB);
00497                                 VPRINT("VtoPole_prj", VtoPole_prj);
00498                         }
00499                         if (VDOT(leftB, VtoPole_prj) >= 0.0) {
00500                                 /* plane point is "inside" edge B */
00501                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00502                                         bu_log("\tplane point inside face\n");
00503                                 VMOVE(Pole_pca, Pole_prj_pt);
00504                                 VSUB2(pca_to_pole_vect, Pole_prj_pt, Pole);
00505                                 *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00506                                 return;
00507                         }
00508                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00509                                         bu_log("\tplane point outside face\n");
00510                         VSUB2(pca_to_pole_vect, pcaB, Pole);
00511                         VMOVE(Pole_pca, pcaB);
00512                 }
00513                 break;
00514         case 0x33: {
00515                 /* The point is over the vertex shared by the points,
00516                  * and not over one edge or the other.  We now need to
00517                  * figure out which edge to classify this point against.
00518                  *
00519                  * Time to do the Tanenbaum algorithm.
00520                  */
00521                 double dotA;
00522                 double dotB;
00523 
00524                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00525                         bu_log("\tplane point beyond both edges.... Doing the Tanenbaum algorithm.\n");
00526 
00527                 VSUB2(VtoPole_prj, Pole_prj_pt, vu->v_p->vg_p->coord);
00528 
00529                 dotA = VDOT(leftA, VtoPole_prj);
00530                 dotB = VDOT(leftB, VtoPole_prj);
00531 
00532                 if (dotA < dotB) {
00533                         if (dotA >= 0.0) {
00534                                 /* Point is "inside" face,
00535                                  * PCA is plane projection point.
00536                                  */
00537                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00538                                         bu_log("\tpoint inside face\n");
00539                                 VMOVE(Pole_pca, Pole_prj_pt);
00540                                 VSUB2(pca_to_pole_vect, Pole, Pole_prj_pt);
00541                         } else {
00542                                 /* Point is "outside" face, PCA is at vertex. */
00543                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00544                                         bu_log("\tpoint outside face, PCA is vertex\n");
00545 
00546                                 VSUB2(pca_to_pole_vect, pcaA, Pole);
00547                                 VMOVE(Pole_pca, pcaA);
00548                         }
00549                 } else {
00550                         if (dotB >= 0.0) {
00551                                 /* Point is "inside" face,
00552                                  * PCA is plane projection point
00553                                  */
00554                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00555                                         bu_log("\tpoint is inside face\n");
00556                                 VMOVE(Pole_pca, Pole_prj_pt);
00557                                 VSUB2(pca_to_pole_vect, Pole, Pole_prj_pt);
00558                         } else {
00559 
00560                                 /* Point is "outside" face, PCA is at vertex. */
00561                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00562                                         bu_log("\tpoint is outside face, PCA is vertex\n");
00563                                 VSUB2(pca_to_pole_vect, pcaB, Pole);
00564                                 VMOVE(Pole_pca, pcaB);
00565                         }
00566                 }
00567                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00568                         VPRINT("\tpca_to_pole_vect", pca_to_pole_vect);
00569                 *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00570                 return;
00571         }
00572         case 0x34: /* fallthrough */
00573         case 0x54: /* fallthrough */
00574         case 0x43: /* fallthrough */
00575         case 0x44: /* fallthrough */
00576         case 0x45: /* fallthrough */
00577         default:
00578                 bu_log("%s %d: So-called 'Impossible' status codes\n",
00579                         __FILE__, __LINE__);
00580                 rt_bomb("get_pole_dist_to_face() Pretending NOT to bomb\n");
00581                 break;
00582         }
00583 
00584         *Pole_dist = MAGNITUDE(pca_to_pole_vect);
00585 
00586         return;
00587 }
00588 
00589 static void
00590 plot_neighborhood(fastf_t *North_Pole, fastf_t *North_pl_pt, fastf_t *North_pca, fastf_t *South_Pole, fastf_t *South_pl_pt, fastf_t *South_pca, fastf_t *pointA, fastf_t *pointB, fastf_t *norm, fastf_t *pt, fastf_t *leftA, fastf_t *leftB)
00591 {
00592         static int plotnum=0;
00593         FILE *pfd;
00594         char name[64];
00595         point_t my_pt;
00596         vect_t ray;
00597 
00598         sprintf(name, "vert%03d.pl", plotnum++);
00599         if ((pfd=fopen(name, "w")) == (FILE *)NULL) {
00600                 bu_log("Error opening %s\n", name);
00601                 return;
00602         } else
00603                 bu_log("overlay %s\n", name);
00604 
00605 
00606         /* draw the ray */
00607         pl_color(pfd, 255, 55, 55);
00608         pdv_3line(pfd, North_Pole, South_Pole);
00609 
00610         /* draw the area of the face */
00611         pl_color(pfd, 55, 255, 55);
00612         pdv_3move(pfd, pt);
00613         pdv_3cont(pfd, pointA);
00614         pdv_3cont(pfd, pointB);
00615         pdv_3cont(pfd, pt);
00616 
00617         /* draw the projections of the pole points */
00618         pl_color(pfd, 255, 255, 55);
00619         pdv_3line(pfd, North_Pole, North_pl_pt);
00620         if ( ! VEQUAL(North_pca, North_pl_pt) ) {
00621                 pdv_3line(pfd, North_Pole, North_pca);
00622                 pdv_3line(pfd, North_pl_pt, North_pca);
00623         }
00624         VSUB2(ray, South_Pole, North_Pole);
00625         VSCALE(ray, ray, -0.125);
00626         VADD2(my_pt, North_Pole, ray);
00627         pdv_3move(pfd, my_pt);
00628         pl_label(pfd, "N");
00629 
00630         pl_color(pfd, 55, 255, 255);
00631         pdv_3line(pfd, South_Pole, South_pl_pt);
00632         if ( ! VEQUAL(South_pca, South_pl_pt) ) {
00633                 pdv_3line(pfd, South_Pole, South_pca);
00634                 pdv_3line(pfd, South_pl_pt, South_pca);
00635         }
00636         VREVERSE(ray, ray);
00637         VADD2(my_pt, South_Pole, ray);
00638         pdv_3move(pfd, my_pt);
00639         pl_label(pfd, "S");
00640 
00641         pl_color(pfd, 128, 128, 128);
00642         VADD2(my_pt, pt, norm);
00643         pdv_3line(pfd, pt, my_pt);
00644 
00645         pl_color(pfd, 192, 192, 192);
00646         VADD2(my_pt, pointA, leftA);
00647         pdv_3line(pfd, pointA, my_pt);
00648 
00649         VADD2(my_pt, pointB, leftB);
00650         pdv_3line(pfd, pointB, my_pt);
00651 
00652 
00653         fclose(pfd);
00654 }
00655 
00656 
00657 
00658 
00659 
00660 
00661 
00662 
00663 /**     V E R T E X _ N E I G H B O R H O O D
00664  *
00665  *      Examine the vertex neighborhood to classify the ray as IN/ON/OUT of
00666  *      the NMG both before and after hitting the vertex.
00667  *
00668  *
00669  *      There is a conceptual sphere about the vertex.  For reasons associated
00670  *      with tolerancing, the sphere has a radius equal to the magnitude of
00671  *      the maximum dimension of the NMG bounding RPP.
00672  *
00673  *      The point where the ray enters this sphere is called the "North Pole"
00674  *      and the point where it exits the sphere is called the "South Pole"
00675  *
00676  *      For each face which uses this vertex we compute 2 "distance" metrics:
00677  *
00678  *      project the "north pole" and "south pole" down onto the plane of the
00679  *      face:
00680  *
00681  *
00682  *
00683  *                                          /
00684  *                                         /
00685  *                                        /North Pole
00686  *                                       / |
00687  *                      Face Normal  ^  /  |
00688  *                                   | /   | Projection of North Pole
00689  *      Plane of face                |/    V to plane
00690  *  ---------------------------------*-------------------------------
00691  *                Projection of ^   / Vertex
00692  *                   South Pole |  /
00693  *                     to plane | /
00694  *                              |/
00695  *                              /South Pole
00696  *                             /
00697  *                            /
00698  *                           /
00699  *                          /
00700  *                         /
00701  *                        /
00702  *                      |/_
00703  *                      Ray
00704  *
00705  *      If the projected polar point is inside the two edgeuses at this
00706  *              vertexuse, then
00707  *                      the distance to the face is the projection distance.
00708  *
00709  *              else
00710  *                      Each of the two edgeuses at this vertexuse
00711  *                      implies an infinite ray originating at the vertex.
00712  *                      Find the point of closest approach (PCA) of each ray
00713  *                      to the projection point.  For whichever ray passes
00714  *                      closer to the projection point, find the distance
00715  *                      from the original pole point to the PCA.  This is
00716  *                      the "distance to the face" metric for this face from
00717  *                      the given pole point.
00718  *
00719  *      We compute this metric for the North and South poles for all faces at
00720  *      the vertex.  The face with the smallest distance to the north (south)
00721  *      pole is used to determine the state of the ray before (after) it hits
00722  *      the vertex.
00723  *
00724  *      If the north (south) pole is "outside" the plane of the closest face,
00725  *      then the ray state is "outside" before (after) the ray hits the
00726  *      vertex.
00727  */
00728 static void
00729 vertex_neighborhood(struct ray_data *rd, struct vertexuse *vu_p, struct hitmiss *myhit)
00730 {
00731         struct vertexuse *vu;
00732         struct faceuse *fu;
00733         point_t South_Pole, North_Pole;
00734         struct faceuse *North_fu = (struct faceuse *)NULL;
00735         struct faceuse *South_fu = (struct faceuse *)NULL;
00736         struct vertexuse *North_vu = (struct vertexuse *)NULL;
00737         struct vertexuse *South_vu = (struct vertexuse *)NULL;
00738         point_t North_pl_pt, South_pl_pt;
00739         point_t North_pca, South_pca;
00740         double North_dist, South_dist;
00741         double North_min, South_min;
00742         double cos_angle;
00743         vect_t VtoPole;
00744         double scalar_dist;
00745         double dimen, t;
00746         point_t min_pt, max_pt;
00747         int i;
00748         vect_t norm, anti_norm;
00749         vect_t polar_height_vect;
00750         vect_t leftA, leftB;
00751         point_t pointA, pointB;
00752         struct edgeuse *eu;
00753         vect_t edge_vect;
00754         int found_faces;
00755 
00756         NMG_CK_VERTEXUSE(vu_p);
00757         NMG_CK_VERTEX(vu_p->v_p);
00758         NMG_CK_VERTEX_G(vu_p->v_p->vg_p);
00759 
00760         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00761                 bu_log("vertex_neighborhood\n");
00762 
00763         nmg_model_bb( min_pt, max_pt, nmg_find_model( vu_p->up.magic_p ));
00764         for (dimen= -MAX_FASTF,i=3 ; i-- ; ) {
00765                 t = max_pt[i]-min_pt[i];
00766                 if (t > dimen) dimen = t;
00767         }
00768 
00769         VJOIN1(North_Pole, vu_p->v_p->vg_p->coord, -dimen, rd->rp->r_dir);
00770         VJOIN1(South_Pole, vu_p->v_p->vg_p->coord, dimen, rd->rp->r_dir);
00771 
00772         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00773                 VPRINT("\tNorth Pole", North_Pole);
00774                 VPRINT("\tSouth Pole", South_Pole);
00775         }
00776 
00777         /* There is a conceptual sphere around the vertex
00778          * The max dimension of the bounding box for the NMG defines the size.
00779          * The point where the ray enters this sphere is
00780          *  called the North Pole, and the point where it
00781          *  exits is called the South Pole.
00782          */
00783 
00784         South_min = North_min = MAX_FASTF;
00785         found_faces = 0;
00786         myhit->in_out = 0;
00787 
00788         /* for every use of this vertex */
00789         for ( BU_LIST_FOR(vu, vertexuse, &vu_p->v_p->vu_hd) ) {
00790                 /* if the parent use is an (edge/loop)use of an
00791                  * OT_SAME faceuse that we haven't already processed...
00792                  */
00793                 NMG_CK_VERTEXUSE(vu);
00794                 NMG_CK_VERTEX(vu->v_p);
00795                 NMG_CK_VERTEX_G(vu->v_p->vg_p);
00796 
00797                 fu = nmg_find_fu_of_vu( vu );
00798                 if (fu) {
00799                   found_faces = 1;
00800                   if(*vu->up.magic_p == NMG_EDGEUSE_MAGIC &&
00801                     fu->orientation == OT_SAME ) {
00802 
00803                         /* The distance from each "Pole Point" to the faceuse
00804                          *  is the commodity in which we are interested.
00805                          *
00806                          * A pole point is projected
00807                          * This is either the distance to the plane (if the
00808                          * projection of the point into the plane lies within
00809                          * the angle of the two edgeuses at this vertex) or
00810                          * we take the distance from the pole to the PCA of
00811                          * the closest edge.
00812                          */
00813                         NMG_GET_FU_NORMAL( norm, fu );
00814                         VREVERSE(anti_norm, norm);
00815 
00816                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00817                                 VPRINT("\tchecking face", norm);
00818 
00819                         /* project the north pole onto the plane */
00820                         VSUB2(polar_height_vect, vu->v_p->vg_p->coord, North_Pole);
00821                         scalar_dist = VDOT(norm, polar_height_vect);
00822 
00823                         /* project the poles down onto the plane */
00824                         VJOIN1(North_pl_pt, North_Pole, scalar_dist, norm);
00825                         VJOIN1(South_pl_pt, South_Pole, scalar_dist, anti_norm);
00826                         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00827                                 VPRINT("\tNorth Plane Point", North_pl_pt);
00828                                 VPRINT("\tSouth Plane Point", South_pl_pt);
00829                         }
00830                         /* Find points on sphere in direction of edges
00831                          * (away from vertex along edge)
00832                          */
00833                         do
00834                                 eu = BU_LIST_PNEXT_CIRC(edgeuse, vu->up.eu_p);
00835                         while (eu->vu_p->v_p == vu->v_p);
00836                         nmg_find_eu_leftvec(leftA, vu->up.eu_p);
00837                         VSUB2(edge_vect, eu->vu_p->v_p->vg_p->coord,
00838                                 vu->v_p->vg_p->coord);
00839                         VUNITIZE(edge_vect);
00840                         VJOIN1(pointA, vu->v_p->vg_p->coord, dimen, edge_vect)
00841 
00842                         do
00843                                 eu = BU_LIST_PPREV_CIRC(edgeuse, vu->up.eu_p);
00844                         while (eu->vu_p->v_p == vu->v_p);
00845                         nmg_find_eu_leftvec(leftB, eu);
00846                         VSUB2(edge_vect, eu->vu_p->v_p->vg_p->coord,
00847                                 vu->v_p->vg_p->coord);
00848                         VUNITIZE(edge_vect);
00849                         VJOIN1(pointB, vu->v_p->vg_p->coord, dimen, edge_vect)
00850 
00851 
00852                         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00853                                 VPRINT("\tLeftA", leftA);
00854                                 VPRINT("\tLeftB", leftB);
00855                         }
00856                         /* find distance of face to North Pole */
00857                         get_pole_dist_to_face(rd, vu,
00858                                 North_Pole, North_pl_pt, &North_dist, North_pca,
00859                                 pointA, leftA, pointB, leftB,
00860                                 polar_height_vect, "North");
00861 
00862                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00863                                 bu_log("\tDist north pole<=>face %g\n", North_dist);
00864 
00865                         if (North_min > North_dist) {
00866                                 North_min = North_dist;
00867                                 North_fu = fu;
00868                                 North_vu = vu;
00869                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00870                                         bu_log("New North Pole Min: %g\n", North_min);
00871                         }
00872 
00873 
00874                         /* find distance of face to South Pole */
00875                         get_pole_dist_to_face(rd, vu,
00876                                 South_Pole, South_pl_pt, &South_dist, South_pca,
00877                                 pointA, leftA, pointB, leftB,
00878                                 polar_height_vect, "South");
00879 
00880                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00881                                 bu_log("\tDist south pole<=>face %g\n", South_dist);
00882 
00883                         if (South_min > South_dist) {
00884                                 South_min = South_dist;
00885                                 South_fu = fu;
00886                                 South_vu = vu;
00887                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00888                                         bu_log("New South Pole Min: %g\n", South_min);
00889                         }
00890 
00891 
00892                         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
00893                                 plot_neighborhood( North_Pole, North_pl_pt, North_pca,
00894                                                    South_Pole, South_pl_pt, South_pca,
00895                                                    pointA, pointB, norm,
00896                                                    vu_p->v_p->vg_p->coord,
00897                                                    leftA, leftB);
00898                         }
00899                     }
00900                 } else {
00901                         if (!found_faces)
00902                                 South_vu = North_vu = vu;
00903                 }
00904         }
00905 
00906         if (!found_faces) {
00907                 /* we've found a vertex floating in space */
00908                 myhit->outbound_use = myhit->inbound_use = (long *)North_vu;
00909                 myhit->in_out = HMG_HIT_ANY_ANY;
00910                 VREVERSE(myhit->hit.hit_normal, rd->rp->r_dir);
00911                 return;
00912         }
00913 
00914 
00915         NMG_GET_FU_NORMAL( norm, North_fu );
00916         VMOVE(myhit->inbound_norm, norm);
00917         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00918                 bu_log("North Pole Min: %g to %g %g %g\n", North_min,
00919                         norm[0], norm[1],norm[2]);
00920 
00921         /* compute status of ray as it is in-bound on the vertex */
00922         VSUB2(VtoPole, North_Pole, vu_p->v_p->vg_p->coord);
00923         cos_angle = VDOT(norm, VtoPole);
00924         if (BN_VECT_ARE_PERP(cos_angle, rd->tol))
00925                 myhit->in_out |= NMG_RAY_STATE_ON << 4;
00926         else if (cos_angle > 0.0)
00927                 myhit->in_out |= NMG_RAY_STATE_OUTSIDE << 4;
00928         else
00929                 myhit->in_out |= NMG_RAY_STATE_INSIDE << 4;
00930 
00931 
00932         /* compute status of ray as it is out-bound from the vertex */
00933         NMG_GET_FU_NORMAL( norm, South_fu );
00934         VMOVE(myhit->outbound_norm, norm);
00935         VSUB2(VtoPole, South_Pole, vu_p->v_p->vg_p->coord);
00936         cos_angle = VDOT(norm, VtoPole);
00937         if (BN_VECT_ARE_PERP(cos_angle, rd->tol))
00938                 myhit->in_out |= NMG_RAY_STATE_ON;
00939         else if (cos_angle > 0.0)
00940                 myhit->in_out |= NMG_RAY_STATE_OUTSIDE;
00941         else
00942                 myhit->in_out |= NMG_RAY_STATE_INSIDE;
00943 
00944 
00945         myhit->inbound_use = (long *)North_vu;
00946         myhit->outbound_use = (long *)South_vu;
00947 
00948         switch (myhit->in_out) {
00949 #if 1
00950         case HMG_HIT_ON_ON:     /* fallthrough???  -MJM??? */
00951 #endif
00952         case HMG_HIT_IN_IN:     /* fallthrough */
00953         case HMG_HIT_OUT_OUT:   /* fallthrough */
00954         case HMG_HIT_IN_ON:     /* fallthrough */
00955         case HMG_HIT_ON_IN:     /* two hits */
00956                 myhit->hit.hit_private = (genptr_t)North_vu;
00957                 break;
00958         case HMG_HIT_IN_OUT:    /* one hit - outbound */
00959         case HMG_HIT_ON_OUT:    /* one hit - outbound */
00960                 myhit->hit.hit_private = (genptr_t)South_vu;
00961                 break;
00962         case HMG_HIT_OUT_IN:    /* one hit - inbound */
00963         case HMG_HIT_OUT_ON:    /* one hit - inbound */
00964                 myhit->hit.hit_private = (genptr_t)North_vu;
00965                 break;
00966         default:
00967                 bu_log("%s %d: vertex_neighborhood() Bad vertex in_out state = x%x\n",
00968                         __FILE__, __LINE__, myhit->in_out);
00969                 rt_bomb("vertex_neighborhood() bad vertex in_out state\n");
00970                 break;
00971 
00972         }
00973 }
00974 
00975 
00976 
00977 
00978 
00979 
00980 /**
00981  *  Once it has been decided that the ray hits the vertex,
00982  *  this routine takes care of recording that fact.
00983  */
00984 static void
00985 ray_hit_vertex(struct ray_data *rd, struct vertexuse *vu_p, int status)
00986 {
00987         struct hitmiss *myhit;
00988         vect_t v;
00989 
00990         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
00991                 bu_log("ray_hit_vertex x%x (%g %g %g)\n",
00992                         vu_p->v_p, V3ARGS( vu_p->v_p->vg_p->coord ));
00993 
00994         if ( (myhit = NMG_INDEX_GET(rd->hitmiss, vu_p->v_p)) ) {
00995                 if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_MAGIC))
00996                         return;
00997                 /* oops, we have to change a MISS into a HIT */
00998                 BU_LIST_DEQUEUE(&myhit->l);
00999         } else {
01000                 GET_HITMISS(myhit, rd->ap);
01001                 NMG_INDEX_ASSIGN(rd->hitmiss, vu_p->v_p, myhit);
01002                 myhit->outbound_use = (long *)vu_p;
01003                 myhit->inbound_use = (long *)vu_p;
01004         }
01005 
01006         /* v = vector from ray point to hit vertex */
01007         VSUB2( v, vu_p->v_p->vg_p->coord, rd->rp->r_pt);
01008         myhit->hit.hit_dist = VDOT( v, rd->rp->r_dir ); /* distance along ray */
01009         VMOVE(myhit->hit.hit_point, vu_p->v_p->vg_p->coord);
01010         myhit->hit.hit_private = (genptr_t) vu_p->v_p;
01011 
01012         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01013                 bu_log( "\tray = ( %g %g %g ), dir=(%g %g %g ), dist=%g\n",
01014                         V3ARGS( rd->rp->r_pt ), V3ARGS( rd->rp->r_dir ), myhit->hit.hit_dist );
01015 
01016         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_MAGIC);
01017 
01018         /* XXX need to compute neighborhood of vertex so that ray
01019          * can be classified as in/on/out before and after the vertex.
01020          *
01021          */
01022         vertex_neighborhood(rd, vu_p, myhit);
01023 
01024         /* XXX we re really should temper the results of vertex_neighborhood()
01025          * with the knowledge of "status"
01026          */
01027 
01028         hit_ins(rd, myhit);
01029         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01030                 plvu(vu_p);
01031 #ifndef FAST_NMG
01032         NMG_CK_HITMISS(myhit);
01033 #endif
01034 }
01035 
01036 
01037 /**     I S E C T _ R A Y _ V E R T E X U S E
01038  *
01039  *      Called in one of the following situations:
01040  *              1)      Zero length edge
01041  *              2)      Vertexuse child of Loopuse
01042  *              3)      Vertexuse child of Shell
01043  *
01044  *      return:
01045  *              1 vertex was hit
01046  *              0 vertex was missed
01047  */
01048 static int
01049 isect_ray_vertexuse(struct ray_data *rd, struct vertexuse *vu_p)
01050 {
01051         struct hitmiss *myhit;
01052         double ray_vu_dist;
01053 
01054         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01055                 bu_log("isect_ray_vertexuse\n");
01056 
01057         NMG_CK_VERTEXUSE(vu_p);
01058         NMG_CK_VERTEX(vu_p->v_p);
01059         NMG_CK_VERTEX_G(vu_p->v_p->vg_p);
01060 
01061         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01062                 bu_log("nmg_isect_ray_vertexuse %g %g %g",
01063                         vu_p->v_p->vg_p->coord[0],
01064                         vu_p->v_p->vg_p->coord[1],
01065                         vu_p->v_p->vg_p->coord[2]);
01066 
01067         if ( (myhit = NMG_INDEX_GET(rd->hitmiss, vu_p->v_p)) ) {
01068                 if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_MAGIC)) {
01069                         /* we've previously hit this vertex */
01070                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01071                                 bu_log(" previously hit\n");
01072                         return(1);
01073                 } else {
01074                         /* we've previously missed this vertex */
01075                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01076                                 bu_log(" previously missed\n");
01077                         return 0;
01078                 }
01079         }
01080 
01081         /* intersect ray with vertex */
01082         ray_vu_dist = bn_dist_line3_pt3(rd->rp->r_pt, rd->rp->r_dir,
01083                                          vu_p->v_p->vg_p->coord);
01084 
01085         if (ray_vu_dist > rd->tol->dist) {
01086                 /* ray misses vertex */
01087                 ray_miss_vertex(rd, vu_p);
01088                 return 0;
01089         }
01090 
01091         /* ray hits vertex */
01092         ray_hit_vertex(rd, vu_p, NMG_VERT_UNKNOWN);
01093 
01094         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01095                 bu_log(" Ray hits vertex, dist %g (priv=x%x, v magic=x%x)\n",
01096                         myhit->hit.hit_dist,
01097                         myhit->hit.hit_private,
01098                         vu_p->v_p->magic);
01099 
01100                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01101                         nmg_rt_print_hitlist(rd->hitmiss[NMG_HIT_LIST]);
01102         }
01103 
01104         return 1;
01105 }
01106 
01107 
01108 /**
01109  * As the name implies, this routine is called when the ray and an edge are
01110  * colinear.  It handles marking the verticies as hit, remembering that this
01111  * is a seg_in/seg_out pair, and builds the hit on the edge.
01112  *
01113  */
01114 static void
01115 colinear_edge_ray(struct ray_data *rd, struct edgeuse *eu_p)
01116 {
01117         struct hitmiss *vhit1, *vhit2, *myhit;
01118 
01119 
01120         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01121                 bu_log("\t - colinear_edge_ray\n");
01122 
01123         vhit1 = NMG_INDEX_GET(rd->hitmiss, eu_p->vu_p->v_p);
01124         vhit2 = NMG_INDEX_GET(rd->hitmiss, eu_p->eumate_p->vu_p->v_p);
01125 
01126         /* record the hit on each vertex, and remember that these two hits
01127          * should be kept together.
01128          */
01129         if (vhit1->hit.hit_dist > vhit2->hit.hit_dist) {
01130                 ray_hit_vertex(rd, eu_p->vu_p, NMG_VERT_ENTER);
01131                 vhit1->start_stop = NMG_HITMISS_SEG_OUT;
01132 
01133                 ray_hit_vertex(rd, eu_p->eumate_p->vu_p, NMG_VERT_LEAVE);
01134                 vhit2->start_stop = NMG_HITMISS_SEG_IN;
01135         } else {
01136                 ray_hit_vertex(rd, eu_p->vu_p, NMG_VERT_LEAVE);
01137                 vhit1->start_stop = NMG_HITMISS_SEG_IN;
01138 
01139                 ray_hit_vertex(rd, eu_p->eumate_p->vu_p, NMG_VERT_ENTER);
01140                 vhit2->start_stop = NMG_HITMISS_SEG_OUT;
01141         }
01142         vhit1->other = vhit2;
01143         vhit2->other = vhit1;
01144 
01145         GET_HITMISS(myhit, rd->ap);
01146         NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit);
01147         myhit->hit.hit_private = (genptr_t)eu_p->e_p;
01148         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC);
01149         BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
01150 #ifndef FAST_NMG
01151         NMG_CK_HITMISS(myhit);
01152         NMG_CK_HITMISS_LISTS(myhit, rd);
01153 #endif
01154         return;
01155 }
01156 
01157 /**
01158  *  When a vertex at an end of an edge gets hit by the ray, this macro
01159  *  is used to build the hit structures for the vertex and the edge.
01160  */
01161 #ifndef FAST_NMG
01162 #define HIT_EDGE_VERTEX(rd, eu_p, vu_p) {\
01163         if (rt_g.NMG_debug & DEBUG_RT_ISECT) bu_log("hit_edge_vertex\n"); \
01164         if (*eu_p->up.magic_p == NMG_SHELL_MAGIC || \
01165             (*eu_p->up.magic_p == NMG_LOOPUSE_MAGIC && \
01166              *eu_p->up.lu_p->up.magic_p == NMG_SHELL_MAGIC)) \
01167                 ray_hit_vertex(rd, vu_p, NMG_VERT_ENTER_LEAVE); \
01168         else \
01169                 ray_hit_vertex(rd, vu_p, NMG_VERT_UNKNOWN); \
01170         GET_HITMISS(myhit); \
01171         NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit); \
01172         myhit->hit.hit_private = (genptr_t)eu_p->e_p; \
01173         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC); \
01174         BU_LIST_INSERT(&rd->rd_miss, &myhit->l); \
01175         NMG_CK_HITMISS(myhit); \
01176         NMG_CK_HITMISS_LISTS(myhit, rd); }
01177 #else
01178 #define HIT_EDGE_VERTEX(rd, eu_p, vu_p) {\
01179         if (rt_g.NMG_debug & DEBUG_RT_ISECT) bu_log("hit_edge_vertex\n"); \
01180         if (*eu_p->up.magic_p == NMG_SHELL_MAGIC || \
01181             (*eu_p->up.magic_p == NMG_LOOPUSE_MAGIC && \
01182              *eu_p->up.lu_p->up.magic_p == NMG_SHELL_MAGIC)) \
01183                 ray_hit_vertex(rd, vu_p, NMG_VERT_ENTER_LEAVE); \
01184         else \
01185                 ray_hit_vertex(rd, vu_p, NMG_VERT_UNKNOWN); \
01186         GET_HITMISS(myhit, rd->ap); \
01187         NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit); \
01188         myhit->hit.hit_private = (genptr_t)eu_p->e_p; \
01189         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC); \
01190         BU_LIST_INSERT(&rd->rd_miss, &myhit->l); }
01191 #endif
01192 
01193 
01194 /**
01195  *      Compute the "ray state" before and after the ray encounters a
01196  *      hit-point on an edge.
01197  */
01198 static void
01199 edge_hit_ray_state(struct ray_data *rd, struct edgeuse *eu, struct hitmiss *myhit)
01200 {
01201         double cos_angle;
01202         double inb_cos_angle = 2.0;
01203         double outb_cos_angle = -1.0;
01204         struct shell *s;
01205         struct faceuse *fu;
01206         struct faceuse *inb_fu = (struct faceuse *)NULL;
01207         struct faceuse *outb_fu = (struct faceuse *)NULL;
01208         struct edgeuse *inb_eu = (struct edgeuse *)NULL;
01209         struct edgeuse *outb_eu = (struct edgeuse *)NULL;
01210         struct edgeuse *eu_p;
01211         struct edgeuse *fu_eu;
01212         vect_t edge_left;
01213         vect_t eu_vec;
01214         vect_t norm;
01215         int faces_found;
01216 
01217 #ifndef FAST_NMG
01218         NMG_CK_HITMISS_LISTS(a_hit, rd);
01219 #endif
01220         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01221                 eu_p = BU_LIST_PNEXT_CIRC(edgeuse, eu);
01222                 bu_log("edge_hit_ray_state(%g %g %g -> %g %g %g _vs_ %g %g %g)\n",
01223                         eu->vu_p->v_p->vg_p->coord[0],
01224                         eu->vu_p->v_p->vg_p->coord[1],
01225                         eu->vu_p->v_p->vg_p->coord[2],
01226                         eu_p->vu_p->v_p->vg_p->coord[0],
01227                         eu_p->vu_p->v_p->vg_p->coord[1],
01228                         eu_p->vu_p->v_p->vg_p->coord[2],
01229                         rd->rp->r_dir[0],
01230                         rd->rp->r_dir[1],
01231                         rd->rp->r_dir[2]);
01232         }
01233 
01234         myhit->in_out = 0;
01235 
01236         s = nmg_find_s_of_eu( eu );
01237 
01238         faces_found = 0;
01239         eu_p = eu->e_p->eu_p;
01240         do {
01241                 if ( (fu=nmg_find_fu_of_eu(eu_p)) ) {
01242                         fu_eu = eu_p;
01243                         faces_found = 1;
01244                         if (fu->orientation == OT_OPPOSITE &&
01245                             fu->fumate_p->orientation == OT_SAME) {
01246                                 fu = fu->fumate_p;
01247                                 fu_eu = eu_p->eumate_p;
01248                             }
01249                         if (fu->orientation != OT_SAME) {
01250                                 bu_log("%s[%d]: I can't seem to find an OT_SAME faceuse\nThis must be a `dangling' face  I'll skip it\n", __FILE__, __LINE__);
01251                                 goto next_edgeuse;
01252                         }
01253 
01254                         if( fu->s_p != s )
01255                                 goto next_edgeuse;
01256 
01257                         if (nmg_find_eu_leftvec(edge_left, eu_p)) {
01258                                 bu_log("edgeuse not part of faceuse");
01259                                 goto next_edgeuse;
01260                         }
01261 
01262                         if (! (NMG_3MANIFOLD &
01263                                NMG_MANIFOLDS(rd->manifolds, fu->f_p)) ){
01264                                 bu_log("This is not a 3-Manifold face.  I'll skip it\n");
01265                                 goto next_edgeuse;
01266                         }
01267 
01268                         cos_angle = VDOT(edge_left, rd->rp->r_dir);
01269 
01270                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01271                         {
01272                                 bu_log("left_vect:(%g %g %g)  cos_angle:%g\n",
01273                                         edge_left[0], edge_left[1],
01274                                         edge_left[2], cos_angle);
01275                         }
01276 
01277                         if (cos_angle < inb_cos_angle) {
01278                                 inb_cos_angle = cos_angle;
01279                                 inb_fu = fu;
01280                                 inb_eu = fu_eu;
01281                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01282                                         bu_log("New inb cos_angle %g\n", inb_cos_angle);
01283                         }
01284                         if (cos_angle > outb_cos_angle) {
01285                                 outb_cos_angle = cos_angle;
01286                                 outb_fu = fu;
01287                                 outb_eu = fu_eu;
01288                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01289                                         bu_log("New outb cos_angle %g\n", outb_cos_angle);
01290                         }
01291                 }
01292 next_edgeuse:   eu_p = eu_p->eumate_p->radial_p;
01293         } while (eu_p != eu->e_p->eu_p);
01294 
01295 #ifndef FAST_NMG
01296         NMG_CK_HITMISS_LISTS(a_hit, rd);
01297 #endif
01298 
01299         if (!faces_found) {
01300                 /* we hit a wire edge */
01301                 myhit->in_out = HMG_HIT_ANY_ANY;
01302                 myhit->hit.hit_private = (genptr_t)eu;
01303                 myhit->inbound_use = myhit->outbound_use = (long *)eu;
01304 
01305                 eu_p = BU_LIST_PNEXT_CIRC(edgeuse, eu);
01306                 VSUB2(eu_vec, eu->vu_p->v_p->vg_p->coord,
01307                         eu_p->vu_p->v_p->vg_p->coord);
01308                 VCROSS(edge_left, eu_vec, rd->rp->r_dir);
01309                 VCROSS(myhit->inbound_norm, eu_vec, edge_left);
01310                 if (VDOT(myhit->inbound_norm, rd->rp->r_dir) > 0.0) {
01311                         VREVERSE(myhit->inbound_norm,myhit->inbound_norm);
01312                 }
01313                 VMOVE(myhit->outbound_norm, myhit->inbound_norm);
01314 
01315                 return;
01316         }
01317 
01318         /* inb_fu is closest to ray on outbound side
01319          * outb_fu is closest to ray on inbound side
01320          */
01321 
01322         /* Compute the ray state on the inbound side */
01323         NMG_GET_FU_NORMAL(norm, inb_fu);
01324         VMOVE(myhit->inbound_norm, norm);
01325         if (MAGSQ(norm) < VDIVIDE_TOL)
01326                 rt_bomb("edge_hit_ray_state() null normal!\n");
01327 
01328         cos_angle = VDOT(norm, rd->rp->r_dir);
01329 
01330         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01331                 VPRINT("\ninb face normal", norm);
01332                 bu_log("cos_angle wrt ray direction: %g\n", cos_angle);
01333         }
01334 
01335         if (BN_VECT_ARE_PERP(cos_angle, rd->tol))
01336                 myhit->in_out |= NMG_RAY_STATE_ON << 4;
01337         else if (cos_angle < 0.0)
01338                 myhit->in_out |= NMG_RAY_STATE_OUTSIDE << 4;
01339         else /* (cos_angle > 0.0) */
01340                 myhit->in_out |= NMG_RAY_STATE_INSIDE << 4;
01341 
01342 
01343         /* Compute the ray state on the outbound side */
01344         NMG_GET_FU_NORMAL(norm, outb_fu);
01345         VMOVE(myhit->outbound_norm, norm);
01346         cos_angle = VDOT(norm, rd->rp->r_dir);
01347 
01348         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01349                 VPRINT("\noutb face normal", norm);
01350                 bu_log("cos_angle wrt ray direction: %g\n", cos_angle);
01351         }
01352 
01353         if (BN_VECT_ARE_PERP(cos_angle, rd->tol))
01354                 myhit->in_out |= NMG_RAY_STATE_ON;
01355         else if (cos_angle > 0.0)
01356                 myhit->in_out |= NMG_RAY_STATE_OUTSIDE;
01357         else /* (cos_angle < 0.0) */
01358                 myhit->in_out |= NMG_RAY_STATE_INSIDE;
01359 
01360 
01361         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01362                 bu_log("myhit->in_out: 0x%02x/", myhit->in_out);
01363                 switch(myhit->in_out) {
01364                 case HMG_HIT_IN_IN:
01365                         bu_log("IN_IN\n"); break;
01366                 case HMG_HIT_IN_OUT:
01367                         bu_log("IN_OUT\n"); break;
01368                 case HMG_HIT_OUT_IN:
01369                         bu_log("OUT_IN\n"); break;
01370                 case HMG_HIT_OUT_OUT:
01371                         bu_log("OUT_OUT\n"); break;
01372                 case HMG_HIT_IN_ON:
01373                         bu_log("IN_ON\n"); break;
01374                 case HMG_HIT_ON_IN:
01375                         bu_log("ON_IN\n"); break;
01376                 case HMG_HIT_OUT_ON:
01377                         bu_log("OUT_ON\n"); break;
01378                 case HMG_HIT_ON_OUT:
01379                         bu_log("ON_OUT\n"); break;
01380                 case HMG_HIT_ON_ON:
01381                         bu_log("ON_ON\n"); break;
01382                 default:
01383                         bu_log("?_?\n"); break;
01384                 }
01385         }
01386 
01387 
01388         switch(myhit->in_out) {
01389         case HMG_HIT_ON_ON:     /* Another fall through?? JRA */
01390         case HMG_HIT_IN_IN:     /* fallthrough */
01391         case HMG_HIT_OUT_OUT:   /* fallthrough */
01392         case HMG_HIT_IN_ON:     /* fallthrough */
01393         case HMG_HIT_ON_IN:     /* two hits */
01394                 myhit->inbound_use = (long *)inb_eu;
01395                 myhit->outbound_use = (long *)outb_eu;
01396                 myhit->hit.hit_private = (genptr_t)inb_eu;
01397                 break;
01398         case HMG_HIT_IN_OUT:    /* one hit - outbound */
01399         case HMG_HIT_ON_OUT:    /* one hit - outbound */
01400                 myhit->inbound_use = (long *)outb_eu;
01401                 myhit->outbound_use = (long *)outb_eu;
01402                 myhit->hit.hit_private = (genptr_t)outb_eu;
01403                 break;
01404         case HMG_HIT_OUT_IN:    /* one hit - inbound */
01405         case HMG_HIT_OUT_ON:    /* one hit - inbound */
01406                 myhit->inbound_use = (long *)inb_eu;
01407                 myhit->outbound_use = (long *)inb_eu;
01408                 myhit->hit.hit_private = (genptr_t)inb_eu;
01409                 break;
01410         default:
01411                 bu_log("%s %d: Bad edge in/out state = x%x\n", __FILE__, __LINE__, myhit->in_out);
01412                 rt_bomb("edge_hit_ray_state() bad edge in_out state\n");
01413                 break;
01414         }
01415 #ifndef FAST_NMG
01416         NMG_CK_HITMISS_LISTS(a_hit, rd);
01417 #endif
01418 }
01419 
01420 /**
01421  *
01422  * record a hit on an edge.
01423  *
01424  */
01425 static void
01426 ray_hit_edge(struct ray_data *rd, struct edgeuse *eu_p, double dist_along_ray, fastf_t *pt)
01427 {
01428         struct hitmiss *myhit;
01429         ray_miss_vertex(rd, eu_p->vu_p);
01430         ray_miss_vertex(rd, eu_p->eumate_p->vu_p);
01431 #if DO_NITMISS_CHECKS
01432         NMG_CK_HITMISS_LISTS(a_hit, rd);
01433 #endif
01434 
01435         if (rt_g.NMG_debug & DEBUG_RT_ISECT) bu_log("\t - HIT edge 0x%08x (edgeuse=x%x)\n", eu_p->e_p, eu_p);
01436 
01437         if ( (myhit = NMG_INDEX_GET(rd->hitmiss, eu_p->e_p)) ) {
01438                 switch (((struct bu_list *)myhit)->magic) {
01439                 case NMG_RT_MISS_MAGIC:
01440                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01441                                 bu_log("\tedge previously missed, changing to hit\n");
01442                         BU_LIST_DEQUEUE(&myhit->l);
01443 #ifndef FAST_NMG
01444                         NMG_CK_HITMISS_LISTS(a_hit, rd);
01445 #endif
01446                         break;
01447                 case NMG_RT_HIT_SUB_MAGIC:
01448                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01449                                 bu_log("\tedge vertex previously hit\n");
01450                         return;
01451                 case NMG_RT_HIT_MAGIC:
01452                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01453                                 bu_log("\tedge previously hit\n");
01454                         return;
01455                 default:
01456                         break;
01457                 }
01458         } else {
01459                 GET_HITMISS(myhit, rd->ap);
01460         }
01461 
01462         /* create hit structure for this edge */
01463         NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit);
01464         myhit->outbound_use = (long *)eu_p;
01465         myhit->inbound_use = (long *)eu_p;
01466 
01467         /* build the hit structure */
01468         myhit->hit.hit_dist = dist_along_ray;
01469         VMOVE(myhit->hit.hit_point, pt)
01470         myhit->hit.hit_private = (genptr_t) eu_p->e_p;
01471 
01472         edge_hit_ray_state(rd, eu_p, myhit);
01473 #ifndef FAST_NMG
01474         NMG_CK_HITMISS_LISTS(a_hit, rd);
01475 #endif
01476         hit_ins(rd, myhit);
01477 #ifndef FAST_NMG
01478         NMG_CK_HITMISS_LISTS(a_hit, rd);
01479 #endif
01480 
01481         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01482                 register struct faceuse *fu;
01483                 if ((fu=nmg_find_fu_of_eu( eu_p )))
01484                         nmg_rt_isect_plfu(fu, rd->rp->r_pt, myhit->hit.hit_point);
01485                 else
01486                         pleu(eu_p, rd->rp->r_pt, myhit->hit.hit_point);
01487         }
01488 #ifndef FAST_NMG
01489         NMG_CK_HITMISS(myhit);
01490         NMG_CK_HITMISS_LISTS(a_hit, rd);
01491 #endif
01492 }
01493 
01494 void
01495 isect_ray_cnurb(struct ray_data *rd, struct edgeuse *eu_p)
01496 {
01497 }
01498 
01499 void
01500 isect_ray_lseg(struct ray_data *rd, struct edgeuse *eu_p)
01501 {
01502         int status;
01503         struct hitmiss *myhit;
01504         int vhit1, vhit2;
01505         double dist_along_ray;
01506 
01507         status = bn_isect_line_lseg( &dist_along_ray,
01508                         rd->rp->r_pt, rd->rp->r_dir,
01509                         eu_p->vu_p->v_p->vg_p->coord,
01510                         eu_p->eumate_p->vu_p->v_p->vg_p->coord,
01511                         rd->tol);
01512 
01513         switch (status) {
01514         case -4 :
01515                 /* Zero length edge.  The routine bn_isect_line_lseg() can't
01516                  * help us.  Intersect the ray with each vertex.  If either
01517                  * vertex is hit, then record that the edge has sub-elements
01518                  * which where hit.  Otherwise, record the edge as being
01519                  * missed.
01520                  */
01521 
01522                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01523                         bu_log("\t - Zero length edge\n");
01524 
01525                 vhit1 = isect_ray_vertexuse(rd, eu_p->vu_p);
01526                 vhit2 = isect_ray_vertexuse(rd, eu_p->eumate_p->vu_p);
01527 
01528                 GET_HITMISS(myhit, rd->ap);
01529                 NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit);
01530 
01531                 myhit->hit.hit_private = (genptr_t)eu_p->e_p;
01532 
01533                 if (vhit1 || vhit2) {
01534                         /* we hit the vertex */
01535                         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC);
01536                 } else {
01537                         /* both vertecies were missed, so edge is missed */
01538                         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
01539                 }
01540                 BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
01541 #ifndef FAST_NMG
01542                 NMG_CK_HITMISS(myhit);
01543                 NMG_CK_HITMISS_LISTS(myhit, rd);
01544 #endif
01545                 break;
01546         case -3 :       /* fallthrough */
01547         case -2 :
01548                 /* The ray misses the edge segment, but hits the infinite
01549                  * line of the edge to one end of the edge segment.  This
01550                  * is an exercise in tabulating the nature of the miss.
01551                  */
01552 
01553                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01554                         bu_log("\t - Miss edge, \"hit\" line\n");
01555                 /* record the fact that we missed each vertex */
01556                 (void)ray_miss_vertex(rd, eu_p->vu_p);
01557                 (void)ray_miss_vertex(rd, eu_p->eumate_p->vu_p);
01558 
01559                 /* record the fact that we missed the edge */
01560                 GET_HITMISS(myhit, rd->ap);
01561                 NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit);
01562                 myhit->hit.hit_private = (genptr_t)eu_p->e_p;
01563 
01564                 BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
01565                 BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
01566 #ifndef FAST_NMG
01567                 NMG_CK_HITMISS(myhit);
01568                 NMG_CK_HITMISS_LISTS(myhit, rd);
01569 #endif
01570                 break;
01571         case -1 : /* just plain missed the edge/line */
01572                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01573                         bu_log("\t - Miss edge/line\n");
01574 
01575                 ray_miss_vertex(rd, eu_p->vu_p);
01576                 ray_miss_vertex(rd, eu_p->eumate_p->vu_p);
01577 
01578                 GET_HITMISS(myhit, rd->ap);
01579                 NMG_INDEX_ASSIGN(rd->hitmiss, eu_p->e_p, myhit);
01580                 myhit->hit.hit_private = (genptr_t)eu_p->e_p;
01581 
01582                 BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
01583                 BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
01584 #ifndef FAST_NMG
01585                 NMG_CK_HITMISS(myhit);
01586                 NMG_CK_HITMISS_LISTS(myhit, rd);
01587 #endif
01588                 break;
01589         case 0 :  /* oh joy.  Lines are co-linear */
01590                 HIT_EDGE_VERTEX(rd, eu_p, eu_p->vu_p);
01591                 HIT_EDGE_VERTEX(rd, eu_p, eu_p->eumate_p->vu_p);
01592                 colinear_edge_ray(rd, eu_p);
01593                 break;
01594         case 1 :
01595                 HIT_EDGE_VERTEX(rd, eu_p, eu_p->vu_p);
01596                 break;
01597         case 2 :
01598                 HIT_EDGE_VERTEX(rd, eu_p, eu_p->eumate_p->vu_p);
01599                 break;
01600         case 3 : /* a hit on an edge */
01601                 {
01602                         point_t pt;
01603 
01604                         VJOIN1(pt, rd->rp->r_pt, dist_along_ray, rd->rp->r_dir);
01605 
01606                         ray_hit_edge(rd, eu_p, dist_along_ray, pt);
01607 
01608                         break;
01609                 }
01610         }
01611 }
01612 /**     I S E C T _ R A Y _ E D G E U S E
01613  *
01614  *      Intersect ray with edgeuse.  If they pass within tolerance, a hit
01615  *      is generated.
01616  *
01617  */
01618 static void
01619 isect_ray_edgeuse(struct ray_data *rd, struct edgeuse *eu_p)
01620 {
01621         struct hitmiss *myhit;
01622 
01623         NMG_CK_EDGEUSE(eu_p);
01624         NMG_CK_EDGEUSE(eu_p->eumate_p);
01625         NMG_CK_EDGE(eu_p->e_p);
01626         NMG_CK_VERTEXUSE(eu_p->vu_p);
01627         NMG_CK_VERTEXUSE(eu_p->eumate_p->vu_p);
01628         NMG_CK_VERTEX(eu_p->vu_p->v_p);
01629         NMG_CK_VERTEX(eu_p->eumate_p->vu_p->v_p);
01630 
01631         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01632                 bu_log("isect_ray_edgeuse (%g %g %g) -> (%g %g %g)",
01633                         eu_p->vu_p->v_p->vg_p->coord[0],
01634                         eu_p->vu_p->v_p->vg_p->coord[1],
01635                         eu_p->vu_p->v_p->vg_p->coord[2],
01636                         eu_p->eumate_p->vu_p->v_p->vg_p->coord[0],
01637                         eu_p->eumate_p->vu_p->v_p->vg_p->coord[1],
01638                         eu_p->eumate_p->vu_p->v_p->vg_p->coord[2]);
01639         }
01640 
01641         if (eu_p->e_p != eu_p->eumate_p->e_p)
01642                 rt_bomb("isect_ray_edgeuse() edgeuse mate has step-father\n");
01643 
01644         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01645                 bu_log("\n\tLooking for previous hit on edge 0x%08x ...\n",
01646                         eu_p->e_p);
01647 
01648         if ( (myhit = NMG_INDEX_GET(rd->hitmiss, eu_p->e_p)) ) {
01649                 if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_MAGIC)) {
01650                         /* previously hit vertex or edge */
01651                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01652                                 bu_log("\tedge previously hit\n");
01653                         return;
01654                 } else if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_HIT_SUB_MAGIC)) {
01655                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01656                                 bu_log("\tvertex of edge previously hit\n");
01657                         return;
01658                 } else if (BU_LIST_MAGIC_OK((struct bu_list *)myhit, NMG_RT_MISS_MAGIC)) {
01659                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01660                                 bu_log("\tedge previously missed\n");
01661                         return;
01662                 } else {
01663                         nmg_rt_bomb(rd, "what happened?\n");
01664                 }
01665         }
01666 
01667         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01668                 bu_log("\t No previous hit\n");
01669 
01670         if ( !eu_p->g.magic_p )
01671                 isect_ray_lseg(rd, eu_p);
01672         else
01673         {
01674                 switch (*eu_p->g.magic_p) {
01675                 case NMG_EDGE_G_LSEG_MAGIC:
01676                         isect_ray_lseg(rd, eu_p);
01677                         break;
01678                 case NMG_EDGE_G_CNURB_MAGIC:
01679                         isect_ray_cnurb(rd, eu_p);
01680                         break;
01681                 }
01682         }
01683 }
01684 
01685 
01686 
01687 
01688 
01689 
01690 /**     I S E C T _ R A Y _ L O O P U S E
01691  *
01692  */
01693 static void
01694 isect_ray_loopuse(struct ray_data *rd, struct loopuse *lu_p)
01695 {
01696         struct edgeuse *eu_p;
01697 
01698         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01699                 bu_log("isect_ray_loopuse( 0x%08x, loop:0x%08x)\n", rd, lu_p->l_p);
01700 
01701         NMG_CK_LOOPUSE(lu_p);
01702         NMG_CK_LOOP(lu_p->l_p);
01703         NMG_CK_LOOP_G(lu_p->l_p->lg_p);
01704 
01705         if (BU_LIST_FIRST_MAGIC(&lu_p->down_hd) == NMG_EDGEUSE_MAGIC) {
01706                 for (BU_LIST_FOR(eu_p, edgeuse, &lu_p->down_hd)) {
01707                         isect_ray_edgeuse(rd, eu_p);
01708                 }
01709                 return;
01710 
01711         } else if (BU_LIST_FIRST_MAGIC(&lu_p->down_hd)!=NMG_VERTEXUSE_MAGIC) {
01712                 bu_log("in %s at %d", __FILE__, __LINE__);
01713                 nmg_rt_bomb(rd, " bad loopuse child magic");
01714         }
01715 
01716         /* loopuse child is vertexuse */
01717 
01718         (void) isect_ray_vertexuse(rd,
01719                 BU_LIST_FIRST(vertexuse, &lu_p->down_hd));
01720 }
01721 
01722 
01723 
01724 #ifndef FAST_NMG
01725 #define FACE_MISS(rd, f_p) {\
01726         struct hitmiss *myhit; \
01727         GET_HITMISS(myhit, rd->ap); \
01728         NMG_INDEX_ASSIGN(rd->hitmiss, f_p, myhit); \
01729         myhit->hit.hit_private = (genptr_t)f_p; \
01730         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC); \
01731         BU_LIST_INSERT(&rd->rd_miss, &myhit->l); \
01732         NMG_CK_HITMISS(myhit); \
01733         NMG_CK_HITMISS_LISTS(myhit, rd); }
01734 #else
01735 #define FACE_MISS(rd, f_p) {\
01736         struct hitmiss *myhit; \
01737         GET_HITMISS(myhit, rd->ap); \
01738         NMG_INDEX_ASSIGN(rd->hitmiss, f_p, myhit); \
01739         myhit->hit.hit_private = (genptr_t)f_p; \
01740         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC); \
01741         BU_LIST_INSERT(&rd->rd_miss, &myhit->l); }
01742 #endif
01743 
01744 
01745 static void
01746 eu_touch_func(struct edgeuse *eu, fastf_t *pt, char *priv)
01747 {
01748         struct edgeuse *eu_next;
01749         struct ray_data *rd;
01750 
01751         NMG_CK_EDGEUSE(eu);
01752         NMG_CK_EDGEUSE(eu->eumate_p);
01753         NMG_CK_EDGE(eu->e_p);
01754         NMG_CK_VERTEXUSE(eu->vu_p);
01755         NMG_CK_VERTEXUSE(eu->eumate_p->vu_p);
01756         NMG_CK_VERTEX(eu->vu_p->v_p);
01757         NMG_CK_VERTEX(eu->eumate_p->vu_p->v_p);
01758 
01759         eu_next = BU_LIST_PNEXT_CIRC(edgeuse, eu);
01760 
01761         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01762                 bu_log("eu_touch(%g %g %g -> %g %g %g\n",
01763                         eu->vu_p->v_p->vg_p->coord[0],
01764                         eu->vu_p->v_p->vg_p->coord[1],
01765                         eu->vu_p->v_p->vg_p->coord[2],
01766                         eu_next->vu_p->v_p->vg_p->coord[0],
01767                         eu_next->vu_p->v_p->vg_p->coord[1],
01768                         eu_next->vu_p->v_p->vg_p->coord[2]);
01769 
01770         rd = (struct ray_data *)priv;
01771         rd->face_subhit = 1;
01772 
01773 #ifndef FAST_NMG
01774         NMG_CK_HITMISS_LISTS(myhit, rd);
01775 #endif
01776         ray_hit_edge(rd, eu, rd->ray_dist_to_plane, pt);
01777 #ifndef FAST_NMG
01778         NMG_CK_HITMISS_LISTS(myhit, rd);
01779 #endif
01780 
01781 }
01782 
01783 
01784 static void
01785 vu_touch_func(struct vertexuse *vu, fastf_t *pt, char *priv)
01786 {
01787         struct ray_data *rd;
01788 
01789         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01790                 bu_log("vu_touch(%g %g %g)\n",
01791                         vu->v_p->vg_p->coord[0],
01792                         vu->v_p->vg_p->coord[1],
01793                         vu->v_p->vg_p->coord[2]);
01794 
01795         rd = (struct ray_data *)priv;
01796 
01797         rd->face_subhit = 1;
01798         ray_hit_vertex(rd, vu, NMG_VERT_UNKNOWN);
01799 }
01800 
01801 static void
01802 record_face_hit(struct ray_data *rd, struct hitmiss *myhit, fastf_t *plane_pt, double dist, struct faceuse *fu_p, fastf_t *norm, int a_hit)
01803 {
01804         double cos_angle;
01805 
01806         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_MAGIC);
01807         myhit->outbound_use = (long *)fu_p;
01808         myhit->inbound_use = (long *)fu_p;
01809 
01810 
01811         VMOVE(myhit->hit.hit_point, plane_pt);
01812 
01813         /* also rd->ray_dist_to_plane */
01814         myhit->hit.hit_dist = dist;
01815         myhit->hit.hit_private = (genptr_t)fu_p->f_p;
01816 
01817 
01818         /* compute what the ray-state is before and after this
01819          * encountering this hit point.
01820          */
01821         cos_angle = VDOT(norm, rd->rp->r_dir);
01822         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
01823                 VPRINT("face Normal", norm);
01824                 bu_log("cos_angle wrt ray direction: %g\n", cos_angle);
01825                 bu_log( "fu x%x manifoldness = %d\n", fu_p, NMG_MANIFOLDS(rd->manifolds, fu_p) );
01826         }
01827 
01828 
01829         if ( ! (NMG_MANIFOLDS(rd->manifolds, fu_p) & NMG_3MANIFOLD) ) {
01830                 myhit->in_out = HMG_HIT_ANY_ANY;
01831 
01832                 if (cos_angle < rd->tol->perp) {
01833                         VMOVE(myhit->inbound_norm, norm);
01834                         VREVERSE(myhit->outbound_norm, norm);
01835                         myhit->inbound_use = (long *)fu_p;
01836                         myhit->outbound_use = (long *)fu_p->fumate_p;
01837                 } else {
01838                         VREVERSE(myhit->inbound_norm, norm);
01839                         VMOVE(myhit->outbound_norm, norm);
01840                         myhit->inbound_use = (long *)fu_p->fumate_p;
01841                         myhit->outbound_use = (long *)fu_p;
01842                 }
01843 
01844                 return;
01845         }
01846 
01847 
01848         switch (fu_p->orientation) {
01849         case OT_SAME:
01850                 if (BN_VECT_ARE_PERP(cos_angle, rd->tol)) {
01851                         /* perpendicular? */
01852                         bu_log("%s[%d]: Ray is in plane of face?\n",
01853                                 __FILE__, __LINE__);
01854                                 rt_bomb("record_face_hit() I quit\n");
01855                 } else if (cos_angle > 0.0) {
01856                         myhit->in_out = HMG_HIT_IN_OUT;
01857                         VREVERSE(myhit->outbound_norm, norm);
01858                         myhit->outbound_use = (long *)fu_p;
01859                         myhit->inbound_use = (long *)fu_p;
01860                 } else {
01861                         myhit->in_out = HMG_HIT_OUT_IN;
01862                         VMOVE(myhit->inbound_norm, norm);
01863                         myhit->inbound_use = (long *)fu_p;
01864                         myhit->outbound_use = (long *)fu_p;
01865                 }
01866                 break;
01867         case OT_OPPOSITE:
01868                 if (BN_VECT_ARE_PERP(cos_angle, rd->tol)) {
01869                         /* perpendicular? */
01870                         bu_log("%s[%d]: Ray is in plane of face?\n",
01871                                 __FILE__, __LINE__);
01872                                 rt_bomb("record_face_hit() I quit\n");
01873                 } else if (cos_angle > 0.0) {
01874                         myhit->in_out = HMG_HIT_OUT_IN;
01875                         VREVERSE(myhit->inbound_norm, norm);
01876                         myhit->inbound_use = (long *)fu_p;
01877                         myhit->outbound_use = (long *)fu_p;
01878                 } else {
01879                         myhit->in_out = HMG_HIT_IN_OUT;
01880                         VMOVE(myhit->outbound_norm, norm);
01881                         myhit->inbound_use = (long *)fu_p;
01882                         myhit->outbound_use = (long *)fu_p;
01883                 }
01884                 break;
01885         default:
01886                 bu_log("%s %d:face orientation not SAME/OPPOSITE\n",
01887                         __FILE__, __LINE__);
01888                 rt_bomb("record_face_hit() Crash and burn\n");
01889         }
01890 
01891         hit_ins(rd, myhit);
01892         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01893                 nmg_rt_isect_plfu(fu_p, rd->rp->r_pt, myhit->hit.hit_point);
01894 
01895 #ifndef FAST_NMG
01896         NMG_CK_HITMISS(myhit);
01897 #endif
01898 
01899 }
01900 
01901 /** this is the UV-space tolerance for the NURB intersector, a larger
01902 number will
01903  * improve speed, but too large will produce errors. Perhaps a routine to calculate
01904  * an appropriate UV_TOL could be constructed.
01905  */
01906 #define UV_TOL  1.0e-6
01907 
01908 void
01909 isect_ray_snurb_face(struct ray_data *rd, struct faceuse *fu, struct face_g_snurb *fg)
01910 {
01911         plane_t pl, pl1, pl2;
01912         struct rt_nurb_uv_hit *hp;
01913         struct bu_list bezier;
01914         struct bu_list hit_list;
01915         struct face_g_snurb *srf;
01916         struct face *f;
01917 
01918         NMG_CK_RD( rd );
01919         NMG_CK_FACE_G_SNURB( fg );
01920         NMG_CK_FACEUSE( fu );
01921 
01922         f = fu->f_p;
01923 
01924         /* calculate two orthogonal planes that intersect along ray */
01925         bn_vec_ortho( pl2, rd->rp->r_dir );
01926         VCROSS( pl1, pl2, rd->rp->r_dir );
01927         pl1[3] = VDOT( rd->rp->r_pt, pl1 );
01928         pl2[3] = VDOT( rd->rp->r_pt, pl2 );
01929 
01930         BU_LIST_INIT( &bezier );
01931         BU_LIST_INIT( &hit_list );
01932 
01933         rt_nurb_bezier( &bezier, fg, rd->ap->a_resource );
01934 
01935         while( BU_LIST_NON_EMPTY( &bezier ) )
01936         {
01937                 point_t srf_min, srf_max;
01938                 int planar;
01939                 point_t hit;
01940                 fastf_t dist;
01941 
01942                 srf = BU_LIST_FIRST( face_g_snurb,  &bezier );
01943                 BU_LIST_DEQUEUE( &srf->l );
01944 
01945                 /* calculate intersection points on NURB surface (in uv space) */
01946                 /* check if NURB surface is a simple planar surface */
01947                 if( srf->order[0] == 2 && srf->order[1] ==2 && srf->s_size[0] ==2 && srf->s_size[1] == 2 )
01948                         planar = 1;
01949                 else
01950                         planar = 0;
01951 
01952                 if( planar )
01953                 {
01954                         vect_t u_dir, v_dir;
01955                         point_t ctl_pt[4];
01956                         vect_t hit_dir;
01957                         int i,j;
01958                         int rational;
01959                         int coords;
01960                         fastf_t *pt;
01961                         fastf_t u, v;
01962 
01963                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01964                                 bu_log( "isect_ray_snurb_face: face is planar\n" );
01965 
01966                         rational = RT_NURB_IS_PT_RATIONAL( srf->pt_type );
01967                         coords = RT_NURB_EXTRACT_COORDS( srf->pt_type );
01968 
01969                         pt = srf->ctl_points;
01970                         for( i=0 ; i<4 ; i++ )
01971                         {
01972                                 for( j=0 ; j<coords ; j++ )
01973                                 {
01974                                         ctl_pt[i][j] = *pt;
01975                                         pt++;
01976                                 }
01977                                 if( rational )
01978                                 {
01979                                         for( j=0 ; j<coords-1 ; j++ )
01980                                                 ctl_pt[i][j] = ctl_pt[i][j]/ctl_pt[i][coords-1];
01981                                 }
01982                         }
01983                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01984                         {
01985                                 for( i=0 ; i<4 ; i++ )
01986                                         bu_log( "\tctl_point[%d] = (%g %g %g)\n", i, V3ARGS( ctl_pt[i] ) );
01987                                 bu_log( "uv range (%g %g) <-> (%g %g)\n", srf->u.knots[0], srf->v.knots[0], srf->u.knots[srf->u.k_size-1], srf->v.knots[srf->v.k_size-1] );
01988                         }
01989 
01990                         VSUB2( u_dir, ctl_pt[1], ctl_pt[0] );
01991                         VSUB2( v_dir, ctl_pt[2], ctl_pt[0] );
01992                         VCROSS( pl, u_dir, v_dir );
01993                         VUNITIZE( pl );
01994                         pl[3] = VDOT( pl, ctl_pt[0] );
01995                         hp = (struct rt_nurb_uv_hit *)NULL;
01996                         if( bn_isect_line3_plane( &dist,  rd->rp->r_pt,  rd->rp->r_dir, pl, rd->tol ) <= 0 )
01997                         {
01998                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
01999                                         bu_log( "\tNo intersection\n" );
02000 
02001                                 rt_nurb_free_snurb( srf, rd->ap->a_resource );
02002                                 continue;
02003                         }
02004 
02005                         VJOIN1( hit, rd->rp->r_pt, dist, rd->rp->r_dir )
02006                         VSUB2( hit_dir, hit, ctl_pt[0] )
02007                         u = srf->u.knots[0] + (srf->u.knots[srf->u.k_size-1] - srf->u.knots[0]) * VDOT( hit_dir, u_dir ) / MAGSQ( u_dir );
02008                         v = srf->v.knots[0] + (srf->v.knots[srf->v.k_size-1] - srf->v.knots[0]) * VDOT( hit_dir, v_dir ) / MAGSQ( v_dir );
02009 
02010                         if( u >= srf->u.knots[0] && u <= srf->u.knots[srf->u.k_size-1] &&
02011                             v >= srf->v.knots[0] && v <= srf->v.knots[srf->v.k_size-1] )
02012                         {
02013                                 hp = (struct rt_nurb_uv_hit *)bu_calloc( 1, sizeof( struct rt_nurb_uv_hit ), "struct rt_nurb_uv_hit" );
02014                                 hp->next = (struct rt_nurb_uv_hit *)NULL;
02015                                 hp->sub = 0;
02016                                 hp->u = u;
02017                                 hp->v = v;
02018 
02019                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02020                                         bu_log( "\thit at uv=(%g %g), xyz=(%g %g %g)\n", hp->u, hp->v, V3ARGS( hit ) );
02021                         }
02022                 }
02023                 else
02024                 {
02025                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02026                                 bu_log( "isect_ray_snurb_face: using planes (%g %g %g %g) (%g %g %g %g)\n",
02027                                         V4ARGS( pl1 ), V4ARGS( pl2 ) );
02028                         (void)rt_nurb_s_bound( srf, srf_min, srf_max );
02029                         if( !rt_in_rpp( rd->rp, rd->rd_invdir, srf_min, srf_max ) )
02030                         {
02031                                 rt_nurb_free_snurb( srf, rd->ap->a_resource );
02032                                 continue;
02033                         }
02034                         hp = rt_nurb_intersect( srf, pl1, pl2, UV_TOL, rd->ap->a_resource );
02035                 }
02036 
02037                 /* process each hit point */
02038                 while( hp != (struct rt_nurb_uv_hit *)NULL )
02039                 {
02040                         struct rt_nurb_uv_hit *next;
02041                         struct hitmiss *myhit;
02042                         vect_t to_hit;
02043                         fastf_t dot;
02044                         fastf_t homo_hit[4];
02045                         int ot_sames, ot_opps;
02046                         struct loopuse *lu;
02047 
02048                         next = hp->next;
02049 
02050                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02051                                 bu_log( "\tintersect snurb surface at uv=(%g %g)\n", hp->u, hp->v );
02052 
02053                         /* check if point is in face (trimming curves) */
02054                         ot_sames = 0;
02055                         ot_opps = 0;
02056 
02057                         for( BU_LIST_FOR( lu, loopuse, &fu->lu_hd ) )
02058                         {
02059                                 if( BU_LIST_FIRST_MAGIC( &lu->down_hd ) != NMG_EDGEUSE_MAGIC )
02060                                         continue;
02061 
02062                                 if( lu->orientation == OT_SAME )
02063                                         ot_sames += nmg_uv_in_lu( hp->u, hp->v, lu );
02064                                 else if( lu->orientation == OT_OPPOSITE )
02065                                         ot_opps += nmg_uv_in_lu( hp->u, hp->v, lu );
02066                                 else
02067                                 {
02068                                         bu_log( "isect_ray_snurb_face: lu orientation = %s!!\n",
02069                                                 nmg_orientation( lu->orientation ) );
02070                                         bu_bomb( "isect_ray_snurb_face: bad lu orientation\n" );
02071                                 }
02072                         }
02073 
02074                         if( ot_sames == 0 || ot_opps == ot_sames )
02075                         {
02076                                 /* not a hit */
02077 
02078                                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02079                                         bu_log( "\tNot a hit\n" );
02080 
02081                                 bu_free( (char *)hp, "hit" );
02082                                 hp = next;
02083                                 continue;
02084                         }
02085 
02086                         GET_HITMISS(myhit, rd->ap);
02087                         NMG_INDEX_ASSIGN(rd->hitmiss, fu->f_p, myhit);
02088                         myhit->hit.hit_private = (genptr_t)fu->f_p;
02089                         myhit->inbound_use = myhit->outbound_use = &fu->l.magic;
02090 
02091                         /* calculate actual hit point (x y z) */
02092                         if( planar )
02093                         {
02094                                 VMOVE( myhit->hit.hit_point, hit )
02095                                 myhit->hit.hit_dist = dist;
02096                                 VMOVE( myhit->hit.hit_normal, pl )
02097                         }
02098                         else
02099                         {
02100                                 rt_nurb_s_eval( srf, hp->u, hp->v, homo_hit );
02101                                 if( RT_NURB_IS_PT_RATIONAL( srf->pt_type ) )
02102                                 {
02103                                         fastf_t inv_homo;
02104 
02105                                         inv_homo = 1.0/homo_hit[3];
02106                                         VSCALE( myhit->hit.hit_point, homo_hit, inv_homo )
02107                                 }
02108                                 else
02109                                         VMOVE( myhit->hit.hit_point, homo_hit )
02110 
02111                                 VSUB2( to_hit, myhit->hit.hit_point, rd->rp->r_pt );
02112                                 myhit->hit.hit_dist = VDOT( to_hit, rd->rp->r_dir );
02113 
02114                                 /* get surface normal */
02115                                 rt_nurb_s_norm( srf, hp->u, hp->v, myhit->hit.hit_normal );
02116                         }
02117 
02118                         /* may need to reverse it */
02119                         if( f->flip )
02120                                 VREVERSE( myhit->hit.hit_normal, myhit->hit.hit_normal )
02121 
02122                         dot = VDOT( myhit->hit.hit_normal, rd->rp->r_dir );
02123 
02124                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02125                         {
02126                                 bu_log( "\thit at dist = %g (%g %g %g), norm = (%g %g %g)\n",
02127                                         myhit->hit.hit_dist,
02128                                         V3ARGS( myhit->hit.hit_point ),
02129                                         V3ARGS( myhit->hit.hit_normal ) );
02130                                 if( dot > 0.0 )
02131                                         bu_log( "\t\tdot = %g (exit point)\n", dot );
02132                                 else
02133                                         bu_log( "\t\tdot = %g (entrance point)\n", dot );
02134                         }
02135 
02136                         switch (fu->orientation) {
02137                         case OT_SAME:
02138                                 if (BN_VECT_ARE_PERP(dot, rd->tol)) {
02139                                         /* perpendicular? */
02140                                         bu_log("%s[%d]: Ray is in plane of face?\n",
02141                                                 __FILE__, __LINE__);
02142                                                 rt_bomb("record_face_hit() I quit\n");
02143                                 } else if (dot > 0.0) {
02144                                         myhit->in_out = HMG_HIT_IN_OUT;
02145                                         VMOVE(myhit->outbound_norm, myhit->hit.hit_normal);
02146                                         myhit->outbound_use = (long *)fu;
02147                                         myhit->inbound_use = (long *)fu;
02148                                 } else {
02149                                         myhit->in_out = HMG_HIT_OUT_IN;
02150                                         VMOVE(myhit->inbound_norm, myhit->hit.hit_normal);
02151                                         myhit->inbound_use = (long *)fu;
02152                                         myhit->outbound_use = (long *)fu;
02153                                 }
02154                                 break;
02155                         case OT_OPPOSITE:
02156                                 if (BN_VECT_ARE_PERP(dot, rd->tol)) {
02157                                         /* perpendicular? */
02158                                         bu_log("%s[%d]: Ray is in plane of face?\n",
02159                                                 __FILE__, __LINE__);
02160                                                 rt_bomb("record_face_hit() I quit\n");
02161                                 } else if (dot > 0.0) {
02162                                         myhit->in_out = HMG_HIT_OUT_IN;
02163                                         VREVERSE(myhit->inbound_norm, myhit->hit.hit_normal);
02164                                         myhit->inbound_use = (long *)fu;
02165                                         myhit->outbound_use = (long *)fu;
02166                                 } else {
02167                                         myhit->in_out = HMG_HIT_IN_OUT;
02168                                         VREVERSE(myhit->outbound_norm, myhit->hit.hit_normal);
02169                                         myhit->inbound_use = (long *)fu;
02170                                         myhit->outbound_use = (long *)fu;
02171                                 }
02172                                 break;
02173                         default:
02174                                 bu_log("%s %d:face orientation not SAME/OPPOSITE\n",
02175                                         __FILE__, __LINE__);
02176                                 rt_bomb("record_face_hit() Crash and burn\n");
02177                         }
02178 
02179                         hit_ins( rd, myhit );
02180 
02181                         bu_free( (char *)hp, "hit" );
02182                         hp = next;
02183                 }
02184                 rt_nurb_free_snurb( srf, rd->ap->a_resource );
02185         }
02186 }
02187 
02188 void
02189 isect_ray_planar_face(struct ray_data *rd, struct faceuse *fu_p, struct face_g_plane *fg_p)
02190 {
02191         plane_t                 norm;
02192         fastf_t                 dist;
02193         struct hitmiss          *myhit;
02194         point_t                 plane_pt;
02195         struct loopuse          *lu_p;
02196         int                     pt_class;
02197 
02198         /* the geometric intersection of the ray with the plane
02199          * of the face has already been done by isect_ray_faceuse().
02200          */
02201         NMG_GET_FU_PLANE( norm, fu_p );
02202 
02203         /* ray hits plane:
02204          *
02205          * Get the ray/plane intersection point.
02206          * Then compute whether this point lies within the area of the face.
02207          */
02208 
02209         VMOVE(plane_pt, rd->plane_pt )
02210         dist = rd->ray_dist_to_plane;
02211 
02212         if (DIST_PT_PLANE(plane_pt, norm) > rd->tol->dist) {
02213                 bu_log("%s:%d plane_pt (%g %g %g) @ dist (%g)out of tolerance\n",
02214                         __FILE__, __LINE__, V3ARGS(plane_pt), dist);
02215                 rt_bomb("isect_ray_planar_face() dist out of tol\n");
02216         }
02217 
02218         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
02219                 double new_dist;
02220                 bu_log("\tray (%16.10e %16.10e %16.10e) (-> %16.10e %16.10e %16.10e)\n",
02221                         rd->rp->r_pt[0],
02222                         rd->rp->r_pt[1],
02223                         rd->rp->r_pt[2],
02224                         rd->rp->r_dir[0],
02225                         rd->rp->r_dir[1],
02226                         rd->rp->r_dir[2]);
02227                 bu_log("\tplane/ray intersection point (%16.10e %16.10e %16.10e)\n",
02228                         V3ARGS(plane_pt));
02229                 bu_log("\tdistance along ray to intersection point %16.10e\n", dist);
02230 
02231                 new_dist=DIST_PT_PLANE(plane_pt, norm);
02232 
02233                 bu_log("\tDIST_PT_PLANE(%16.10e) 0x%08lx 0x%08lx\n", new_dist,
02234                         new_dist);
02235 
02236                 bn_isect_line3_plane(&new_dist, plane_pt, rd->rp->r_dir,
02237                         norm, rd->tol);
02238 
02239                 bu_log("Normal %16.10e %16.10e %16.10e %16.10e)\n",
02240                         V4ARGS(norm));
02241                 bu_log("recalculated plane/pt dist as %16.10e 0x%08lx 0x%08lx\n",
02242                         new_dist, new_dist);
02243                 bu_log("distance tol = %16.10e\n", rd->tol->dist);
02244         }
02245 
02246 
02247         /* determine if the plane point is in or out of the face, and
02248          * if it is within tolerance of any of the elements of the faceuse.
02249          *
02250          * The value of "rd->face_subhit" will be set non-zero if either
02251          * eu_touch_func or vu_touch_func is called.  They will be called
02252          * when an edge/vertex of the face is within tolerance of plane_pt.
02253          */
02254         rd->face_subhit = 0;
02255         rd->ray_dist_to_plane = dist;
02256         if( rd->classifying_ray )
02257                 pt_class = nmg_class_pt_fu_except(plane_pt, fu_p, (struct loopuse *)NULL,
02258                         0, 0, (char *)rd, NMG_FPI_PERGEOM, 1,
02259                         rd->tol);
02260         else
02261                 pt_class = nmg_class_pt_fu_except(plane_pt, fu_p, (struct loopuse *)NULL,
02262                         eu_touch_func, vu_touch_func, (char *)rd, NMG_FPI_PERGEOM, 0,
02263                         rd->tol);
02264 
02265 
02266         GET_HITMISS(myhit, rd->ap);
02267         NMG_INDEX_ASSIGN(rd->hitmiss, fu_p->f_p, myhit);
02268         myhit->hit.hit_private = (genptr_t)fu_p->f_p;
02269         myhit->inbound_use = myhit->outbound_use = &fu_p->l.magic;
02270 
02271 
02272 
02273         switch (pt_class) {
02274         case NMG_CLASS_Unknown  :
02275                 bu_log("%s[line:%d] ray/plane intercept point cannot be classified wrt face\n",
02276                         __FILE__, __LINE__);
02277                 rt_bomb("isect_ray_planar_face() class unknown\n");
02278                 break;
02279         case NMG_CLASS_AinB     :
02280         case NMG_CLASS_AonBshared :
02281                 /* if a sub-element of the face was within tolerance of the
02282                  * plane intercept, then a hit has already been recorded for
02283                  * that element, and we do NOT need to generate one for the
02284                  * face.
02285                  */
02286                 if (rd->face_subhit) {
02287                         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_HIT_SUB_MAGIC);
02288                         VMOVE(myhit->hit.hit_point, plane_pt);
02289                         /* also rd->ray_dist_to_plane */
02290                         myhit->hit.hit_dist = dist;
02291 
02292                         myhit->hit.hit_private = (genptr_t)fu_p->f_p;
02293                         BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
02294 #ifndef FAST_NMG
02295                         NMG_CK_HITMISS(myhit);
02296 #endif
02297                 } else {
02298                         /* The plane_pt was NOT within tolerance of a
02299                          * sub-element, but it WAS within the area of the
02300                          * face.  We need to record a hit on the face
02301                          */
02302                         record_face_hit(rd, myhit, plane_pt, dist, fu_p, norm, 0);
02303                 }
02304                 break;
02305         case NMG_CLASS_AoutB    :
02306                 BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
02307                 BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
02308 #ifndef FAST_NMG
02309                 NMG_CK_HITMISS(myhit);
02310 #endif
02311 
02312                 break;
02313         default :
02314                 bu_log("%s[line:%d] BIZZARE ray/plane intercept point classification\n",
02315                         __FILE__, __LINE__);
02316                 rt_bomb("isect_ray_planar_face() Bizz\n");
02317         }
02318 
02319         /* intersect the ray with the edges/verticies of the face */
02320         for ( BU_LIST_FOR(lu_p, loopuse, &fu_p->lu_hd) )
02321                 isect_ray_loopuse(rd, lu_p);
02322 
02323 }
02324 
02325 
02326 
02327 
02328 
02329 
02330 /**     I S E C T _ R A Y _ F A C E U S E
02331  *
02332  *      check to see if ray hits face.
02333  */
02334 static void
02335 isect_ray_faceuse(struct ray_data *rd, struct faceuse *fu_p)
02336 {
02337 
02338         struct hitmiss          *myhit;
02339         struct face             *fp;
02340         struct face_g_plane     *fgp;
02341 
02342         if (fu_p->orientation == OT_OPPOSITE) return;
02343 
02344         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02345                 bu_log("isect_ray_faceuse(0x%08x, faceuse:0x%08x/face:0x%08x)\n",
02346                         rd, fu_p, fu_p->f_p);
02347 
02348         NMG_CK_FACEUSE(fu_p);
02349         NMG_CK_FACEUSE(fu_p->fumate_p);
02350         fp = fu_p->f_p;
02351         NMG_CK_FACE(fp);
02352 
02353 
02354         /* if this face already processed, we are done. */
02355         if ( (myhit = NMG_INDEX_GET(rd->hitmiss, fp)) ) {
02356                 if (BU_LIST_MAGIC_OK((struct bu_list *)myhit,
02357                     NMG_RT_HIT_MAGIC)) {
02358                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02359                                 bu_log(" previously hit\n");
02360                 } else if (BU_LIST_MAGIC_OK((struct bu_list *)myhit,
02361                     NMG_RT_HIT_SUB_MAGIC)) {
02362                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02363                                 bu_log(" previously hit sub-element\n");
02364                 } else if (BU_LIST_MAGIC_OK((struct bu_list *)myhit,
02365                     NMG_RT_MISS_MAGIC)) {
02366                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02367                                 bu_log(" previously missed\n");
02368                 } else {
02369                         bu_log("%s %d:\n\tBad magic %ld (0x%08x) for hitmiss struct for faceuse 0x%08x\n",
02370                                 __FILE__, __LINE__,
02371                                 myhit->l.magic, myhit->l.magic, fu_p);
02372                         nmg_rt_bomb(rd, "Was I hit or not?\n");
02373                 }
02374                 return;
02375         }
02376 
02377 
02378         /* bounding box intersection */
02379         if( *fp->g.magic_p == NMG_FACE_G_PLANE_MAGIC )
02380         {
02381                 int code;
02382                 fastf_t dist;
02383                 point_t hit_pt;
02384 
02385                 fgp = fu_p->f_p->g.plane_p;
02386                 NMG_CK_FACE_G_PLANE(fgp);
02387 
02388                 code = bn_isect_line3_plane( &dist, rd->rp->r_pt, rd->rp->r_dir, fgp->N, rd->tol );
02389                 if( code < 1 )
02390                 {
02391                         GET_HITMISS(myhit, rd->ap);
02392                         NMG_INDEX_ASSIGN(rd->hitmiss, fu_p->f_p, myhit);
02393                         myhit->hit.hit_private = (genptr_t)fu_p->f_p;
02394                         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
02395                         BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
02396 #ifndef FAST_NMG
02397                         NMG_CK_HITMISS(myhit);
02398 #endif
02399 
02400                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02401                                 bu_log("missed bounding box\n");
02402                         return;
02403                 }
02404                 VJOIN1( hit_pt, rd->rp->r_pt, dist, rd->rp->r_dir )
02405                 if( !V3PT_IN_RPP( hit_pt, fp->min_pt, fp->max_pt ) )
02406                 {
02407                         GET_HITMISS(myhit, rd->ap);
02408                         NMG_INDEX_ASSIGN(rd->hitmiss, fu_p->f_p, myhit);
02409                         myhit->hit.hit_private = (genptr_t)fu_p->f_p;
02410                         BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
02411                         BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
02412 #ifndef FAST_NMG
02413                         NMG_CK_HITMISS(myhit);
02414 #endif
02415 
02416                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02417                                 bu_log("hit point not within face bounding box\n");
02418                         return;
02419                 }
02420 
02421                 VMOVE( rd->plane_pt, hit_pt )
02422                 rd->ray_dist_to_plane = dist;
02423         }
02424         else if (!rt_in_rpp(rd->rp, rd->rd_invdir,
02425             fu_p->f_p->min_pt, fu_p->f_p->max_pt) ) {
02426                 GET_HITMISS(myhit, rd->ap);
02427                 NMG_INDEX_ASSIGN(rd->hitmiss, fu_p->f_p, myhit);
02428                 myhit->hit.hit_private = (genptr_t)fu_p->f_p;
02429                 BU_LIST_MAGIC_SET(&myhit->l, NMG_RT_MISS_MAGIC);
02430                 BU_LIST_INSERT(&rd->rd_miss, &myhit->l);
02431 #ifndef FAST_NMG
02432                 NMG_CK_HITMISS(myhit);
02433 #endif
02434 
02435                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02436                         bu_log("missed bounding box\n");
02437                 return;
02438         }
02439 
02440         if (rt_g.NMG_debug & DEBUG_RT_ISECT) bu_log(" hit bounding box \n");
02441 
02442 
02443         switch (*fu_p->f_p->g.magic_p) {
02444         case NMG_FACE_G_PLANE_MAGIC:
02445                 isect_ray_planar_face(rd, fu_p, fu_p->f_p->g.plane_p);
02446                 break;
02447         case NMG_FACE_G_SNURB_MAGIC:
02448                 isect_ray_snurb_face(rd, fu_p, fu_p->f_p->g.snurb_p);
02449                 break;
02450         }
02451 }
02452 
02453 
02454 /**     I S E C T _ R A Y _ S H E L L
02455  *
02456  *      Implicit return:
02457  *              adds hit points to the hit-list "hl"
02458  */
02459 static void
02460 nmg_isect_ray_shell(struct ray_data *rd, const struct shell *s_p)
02461 {
02462         struct faceuse *fu_p;
02463         struct loopuse *lu_p;
02464         struct edgeuse *eu_p;
02465 
02466         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02467                 bu_log("nmg_isect_ray_shell( 0x%08x, 0x%08x)\n", rd, s_p);
02468 
02469         NMG_CK_SHELL(s_p);
02470         NMG_CK_SHELL_A(s_p->sa_p);
02471 
02472         /* does ray isect shell rpp ?
02473          * if not, we can just return.  there is no need to record the
02474          * miss for the shell, as there is only one "use" of a shell.
02475          */
02476         if (!rt_in_rpp(rd->rp, rd->rd_invdir,
02477             s_p->sa_p->min_pt, s_p->sa_p->max_pt) ) {
02478                 if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02479                         bu_log("nmg_isect_ray_shell( no RPP overlap)\n",
02480                                  rd, s_p);
02481                 return;
02482         }
02483 
02484         /* ray intersects shell, check sub-objects */
02485 
02486         for (BU_LIST_FOR(fu_p, faceuse, &(s_p->fu_hd)) )
02487                 isect_ray_faceuse(rd, fu_p);
02488 
02489         for (BU_LIST_FOR(lu_p, loopuse, &(s_p->lu_hd)) )
02490                 isect_ray_loopuse(rd, lu_p);
02491 
02492         for (BU_LIST_FOR(eu_p, edgeuse, &(s_p->eu_hd)) )
02493                 isect_ray_edgeuse(rd, eu_p);
02494 
02495         if (s_p->vu_p)
02496                 (void)isect_ray_vertexuse(rd, s_p->vu_p);
02497 
02498         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02499                 bu_log("nmg_isect_ray_shell( done )\n", rd, s_p);
02500 }
02501 
02502 
02503 /**     N M G _ I S E C T _ R A Y _ M O D E L
02504  *
02505  */
02506 void
02507 nmg_isect_ray_model(struct ray_data *rd)
02508 {
02509         struct nmgregion *r_p;
02510         struct shell *s_p;
02511 
02512 
02513         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02514                 bu_log("isect_ray_nmg: Pnt(%g %g %g) Dir(%g %g %g)\n",
02515                         rd->rp->r_pt[0],
02516                         rd->rp->r_pt[1],
02517                         rd->rp->r_pt[2],
02518                         rd->rp->r_dir[0],
02519                         rd->rp->r_dir[1],
02520                         rd->rp->r_dir[2]);
02521 
02522         NMG_CK_MODEL(rd->rd_m);
02523 
02524         /* Caller has assured us that the ray intersects the nmg model,
02525          * check ray for intersecion with rpp's of nmgregion's
02526          */
02527         for ( BU_LIST_FOR(r_p, nmgregion, &rd->rd_m->r_hd) ) {
02528                 NMG_CK_REGION(r_p);
02529                 NMG_CK_REGION_A(r_p->ra_p);
02530 
02531                 /* does ray intersect nmgregion rpp? */
02532                 if (! rt_in_rpp(rd->rp, rd->rd_invdir,
02533                     r_p->ra_p->min_pt, r_p->ra_p->max_pt) )
02534                         continue;
02535 
02536                 /* ray intersects region, check shell intersection */
02537                 for (BU_LIST_FOR(s_p, shell, &r_p->s_hd)) {
02538                         nmg_isect_ray_shell(rd, s_p);
02539                 }
02540         }
02541 
02542 #ifndef FAST_NMG
02543         NMG_CK_HITMISS_LISTS(a_hit, rd);
02544 #endif
02545 
02546         if (rt_g.NMG_debug & DEBUG_RT_ISECT) {
02547                 if (BU_LIST_IS_EMPTY(&rd->rd_hit)) {
02548                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02549                                 bu_log("ray missed NMG\n");
02550                 } else {
02551                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02552                                 nmg_rt_print_hitlist((struct hitmiss*)&rd->rd_hit);
02553                 }
02554         }
02555 }
02556 
02557 /**
02558  *                              N M G _ P L _ H I T M I S S _ L I S T
02559  */
02560 void
02561 nmg_pl_hitmiss_list(const char *str, int num, const struct bu_list *hd, const struct xray *rp)
02562 {
02563         FILE            *fp;
02564         char            buf[128];
02565         struct hitmiss  *hmp;
02566         int             count = 0;
02567 
02568         sprintf( buf, "%s%d.pl", str, num );
02569 
02570         if( bu_list_len(hd) <= 0 )  {
02571                 bu_log("nmg_pl_hitmiss_list(): empty list, no %s written\n", buf);
02572                 return;
02573         }
02574 
02575         if( (fp = fopen(buf, "w")) == (FILE *)NULL )  {
02576                 perror(buf);
02577                 return;
02578         }
02579 
02580         pl_color( fp, 210, 210, 210 );          /* grey ray */
02581 
02582         for( BU_LIST_FOR( hmp, hitmiss, hd ) )  {
02583                 point_t         pt;
02584 #ifndef FAST_NMG
02585                 NMG_CK_HITMISS(hmp);
02586 #endif
02587                 VJOIN1( pt, rp->r_pt, hmp->hit.hit_dist, rp->r_dir );
02588                 if( count++ )
02589                         pdv_3cont( fp, pt );
02590                 else
02591                         pdv_3move( fp, pt );
02592         }
02593         fclose(fp);
02594         bu_log("overlay %s\n", buf);
02595 }
02596 
02597 static int
02598 guess_class_from_hitlist_max(struct ray_data *rd, int *hari_kari, int in_or_out_only)
02599 {
02600         struct hitmiss *a_hit;
02601         struct hitmiss *plus_hit = (struct hitmiss *)NULL;
02602         int pt_class;
02603 
02604         *hari_kari = 0;
02605 
02606         NMG_CK_RD(rd);
02607         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02608                 bu_log("plus guessing\n");
02609 
02610         for (BU_LIST_FOR(a_hit, hitmiss, &rd->rd_hit)) {
02611 
02612 #ifndef FAST_NMG
02613                 NMG_CK_HITMISS(a_hit);
02614 #endif
02615                 if( !in_or_out_only )
02616                 {
02617                         /* if we've got a zero distance hit, that clinches it */
02618                         if (a_hit->hit.hit_dist <= rd->tol->dist &&
02619                             a_hit->hit.hit_dist >= -rd->tol->dist) {
02620                                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02621                                         bu_log("guess_class_from_hitlist_min() returns NMG_CLASS_AonBshared for 0 dist hit\n");
02622 
02623                                 return NMG_CLASS_AonBshared;
02624                         }
02625 
02626                         if (a_hit->hit.hit_dist < -rd->tol->dist)
02627                                 continue;
02628                 }
02629                 else if (a_hit->hit.hit_dist < 0.0)
02630                         continue;
02631 
02632                 if (a_hit->in_out == HMG_HIT_ANY_ANY)
02633                         continue;
02634 
02635                 if (plus_hit == (struct hitmiss *)NULL) {
02636                         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02637                                 bu_log("plus hit = %g (%s)\n", a_hit->hit.hit_dist,
02638                                         nmg_rt_inout_str(a_hit->in_out));
02639                         plus_hit = a_hit;
02640                         *hari_kari = 0;
02641                 } else if (a_hit->hit.hit_dist < plus_hit->hit.hit_dist) {
02642                         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02643                                 bu_log("plus hit = %g (%s)\n", a_hit->hit.hit_dist,
02644                                         nmg_rt_inout_str(a_hit->in_out));
02645                         plus_hit = a_hit;
02646                         *hari_kari = 0;
02647                 } else if ( a_hit->hit.hit_dist == plus_hit->hit.hit_dist) {
02648                         *hari_kari = 1;
02649                 }
02650         }
02651 
02652         /* XXX This needs to be resolved with parity */
02653         if (*hari_kari) {
02654                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02655                         bu_log("Contemplating Hari Kari\n");
02656                 return NMG_CLASS_Unknown;
02657         }
02658         /* figure out what the status is from plus_hit */
02659         if (plus_hit) {
02660                 switch (plus_hit->in_out) {
02661                 case HMG_HIT_IN_IN:
02662                 case HMG_HIT_IN_OUT:
02663                 case HMG_HIT_IN_ON:
02664                         pt_class = NMG_CLASS_AinB;
02665                         break;
02666                 case HMG_HIT_OUT_IN:
02667                 case HMG_HIT_OUT_OUT:
02668                 case HMG_HIT_OUT_ON:
02669                         pt_class = NMG_CLASS_AoutB;
02670                         break;
02671                 case HMG_HIT_ON_IN:
02672                 case HMG_HIT_ON_ON:
02673                 case HMG_HIT_ON_OUT:
02674                         pt_class = NMG_CLASS_AonBshared;
02675                         break;
02676                 default:
02677                         rt_bomb("guess_class_from_hitlist_max() no-class hitpoint\n");
02678                         pt_class = 0; /* shuts up compiler warning */
02679                         break;
02680                 }
02681         } else {
02682                 /* since we didn't hit anything in the positive direction,
02683                  * we've got to be outside, since we don't allow infinite
02684                  * NMG's
02685                  */
02686                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02687                         bu_log("Nothing in the plus direction\n");
02688                 pt_class = NMG_CLASS_AoutB;
02689         }
02690 
02691         return pt_class;
02692 }
02693 
02694 static int
02695 guess_class_from_hitlist_min(struct ray_data *rd, int *hari_kari, int in_or_out_only)
02696 {
02697         struct hitmiss *a_hit;
02698         struct hitmiss *minus_hit = (struct hitmiss *)NULL;
02699         int pt_class;
02700 
02701         *hari_kari = 0;
02702 
02703         NMG_CK_RD(rd);
02704         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02705                 bu_log("minus guessing\n");
02706 
02707         for (BU_LIST_FOR(a_hit, hitmiss, &rd->rd_hit)) {
02708 
02709 #ifndef FAST_NMG
02710                 NMG_CK_HITMISS(a_hit);
02711 #endif
02712                 if( !in_or_out_only )
02713                 {
02714                         /* if we've got a zero distance hit, that clinches it */
02715                         if (a_hit->hit.hit_dist <= rd->tol->dist &&
02716                             a_hit->hit.hit_dist >= -rd->tol->dist) {
02717                                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02718                                         bu_log("guess_class_from_hitlist_min() returns NMG_CLASS_AonBshared for 0 dist hit\n");
02719 
02720                                 return NMG_CLASS_AonBshared;
02721                         }
02722 
02723                         if (a_hit->hit.hit_dist > rd->tol->dist)
02724                                 continue;
02725                 }
02726                 else if (a_hit->hit.hit_dist > 0.0)
02727                         continue;
02728 
02729                 if (a_hit->in_out == HMG_HIT_ANY_ANY)
02730                         continue;
02731 
02732                 if (minus_hit == (struct hitmiss *)NULL) {
02733                         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02734                                 bu_log("minus hit = %g (%s)\n", a_hit->hit.hit_dist,
02735                                         nmg_rt_inout_str(a_hit->in_out));
02736                         minus_hit = a_hit;
02737                         *hari_kari = 0;
02738                 } else if (a_hit->hit.hit_dist > minus_hit->hit.hit_dist) {
02739                         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02740                                 bu_log("minus hit = %g (%s)\n", a_hit->hit.hit_dist,
02741                                         nmg_rt_inout_str(a_hit->in_out));
02742                         minus_hit = a_hit;
02743                         *hari_kari = 0;
02744                 } else if ( a_hit->hit.hit_dist == minus_hit->hit.hit_dist) {
02745                         *hari_kari = 1;
02746                 }
02747         }
02748 
02749         /* XXX This needs to be resolved with parity */
02750         if (*hari_kari) {
02751                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02752                         bu_log("Contemplating Hari Kari\n");
02753                 return NMG_CLASS_Unknown;
02754         }
02755 
02756         /* figure out what the status is from plus_hit */
02757         if (minus_hit) {
02758                 switch (minus_hit->in_out) {
02759                 case HMG_HIT_IN_IN:
02760                 case HMG_HIT_OUT_IN:
02761                 case HMG_HIT_ON_IN:
02762                         pt_class = NMG_CLASS_AinB;
02763                         break;
02764                 case HMG_HIT_OUT_OUT:
02765                 case HMG_HIT_IN_OUT:
02766                 case HMG_HIT_ON_OUT:
02767                         pt_class = NMG_CLASS_AoutB;
02768                         break;
02769                 case HMG_HIT_ON_ON:
02770                 case HMG_HIT_OUT_ON:
02771                 case HMG_HIT_IN_ON:
02772                         pt_class = NMG_CLASS_AonBshared;
02773                         break;
02774                 default:
02775                         rt_bomb("guess_class_from_hitlist_min() no-class hitpoint\n");
02776                         pt_class = 0; /* shuts up compiler warning */
02777                         break;
02778                 }
02779         } else {
02780                 /* since we didn't hit anything in this direction,
02781                  * we've got to be outside, since we don't allow infinite
02782                  * NMG's
02783                  */
02784                 if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02785                         bu_log("Nothing in the minus direction\n");
02786                 pt_class = NMG_CLASS_AoutB;
02787         }
02788 
02789         return pt_class;
02790 }
02791 
02792 
02793 /**
02794  *      N M G _ R A Y _ I S E C T _ S H E L L
02795  *
02796  *      Intended as a support routine for nmg_class_pt_s() in nmg_class.c
02797  *
02798  *      Intersect a ray with a shell, and return whether the ray start
02799  *      point is inside or outside or ON the shell.
02800  *      Count the number of crossings (hit points) along the ray,
02801  *      both in the negative and positive directions.
02802  *      If an even number, point is outside, if an odd number point is inside.
02803  *      If the negative-going and positive-going assessments don't agree,
02804  *      this is a problem.
02805  *
02806  *      If "in_or_out_only" is non-zero, then we will not look for a
02807  *      classification of "ON" the shell.
02808  *
02809  *      The caller must be prepared for a return of NMG_CLASS_Unknown,
02810  *      in which case it should pick a less difficult ray direction to fire
02811  *      and try again.
02812  *
02813  *  Returns -
02814  *      NMG_CLASS_Unknown       Can't tell
02815  *      NMG_CLASS_xxx           Classification of the pt w.r.t. the shell.
02816  */
02817 int
02818 nmg_class_ray_vs_shell(struct xray *rp, const struct shell *s, const int in_or_out_only, const struct bn_tol *tol)
02819 {
02820         struct ray_data rd;
02821         struct application ap;
02822         struct hitmiss *a_hit;
02823         int minus_class, plus_class;
02824         int hari_kari_minus, hari_kari_plus;
02825 
02826 
02827         VUNITIZE(rp->r_dir);
02828         NMG_CK_SHELL(s);
02829         BN_CK_TOL(tol);
02830 
02831         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT)) {
02832                 bu_log("nmg_class_ray_vs_shell(pt(%g %g %g), dir(%g %g %g))\n",
02833                         V3ARGS(rp->r_pt), V3ARGS(rp->r_dir));
02834         }
02835 
02836         RT_APPLICATION_INIT(&ap);
02837 
02838         ap.a_resource = &rt_uniresource;
02839         rt_uniresource.re_magic = RESOURCE_MAGIC;
02840 
02841         if( BU_LIST_UNINITIALIZED( &rt_uniresource.re_nmgfree ) )
02842                 BU_LIST_INIT( &rt_uniresource.re_nmgfree );
02843 
02844         rd.rd_m = nmg_find_model( &s->l.magic );
02845         rd.manifolds = nmg_manifolds(rd.rd_m);
02846 
02847         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
02848         {
02849                 struct faceuse *fu;
02850 
02851                 bu_log( "Manifoldness for shell FU's\n" );
02852                 for( BU_LIST_FOR( fu, faceuse, &s->fu_hd ) )
02853                 {
02854                         if( fu->orientation != OT_SAME )
02855                                 continue;
02856 
02857                         bu_log( "fu x%x: %d\n", fu, NMG_MANIFOLDS( rd.manifolds, fu ) );
02858                 }
02859         }
02860 
02861         /* Compute the inverse of the direction cosines */
02862         if( !NEAR_ZERO( rp->r_dir[X], SQRT_SMALL_FASTF ) )  {
02863                 rd.rd_invdir[X]=1.0/rp->r_dir[X];
02864         } else {
02865                 rd.rd_invdir[X] = INFINITY;
02866                 rp->r_dir[X] = 0.0;
02867         }
02868         if( !NEAR_ZERO( rp->r_dir[Y], SQRT_SMALL_FASTF ) )  {
02869                 rd.rd_invdir[Y]=1.0/rp->r_dir[Y];
02870         } else {
02871                 rd.rd_invdir[Y] = INFINITY;
02872                 rp->r_dir[Y] = 0.0;
02873         }
02874         if( !NEAR_ZERO( rp->r_dir[Z], SQRT_SMALL_FASTF ) )  {
02875                 rd.rd_invdir[Z]=1.0/rp->r_dir[Z];
02876         } else {
02877                 rd.rd_invdir[Z] = INFINITY;
02878                 rp->r_dir[Z] = 0.0;
02879         }
02880 
02881         rd.rp = rp;
02882         rd.tol = tol;
02883         rd.ap = &ap;
02884         rd.stp = (struct soltab *)NULL;
02885         rd.seghead = (struct seg *)NULL;
02886         rd.magic = NMG_RAY_DATA_MAGIC;
02887         rd.hitmiss = (struct hitmiss **)bu_calloc( rd.rd_m->maxindex,
02888                 sizeof(struct hitmiss *), "nmg geom hit list");
02889         rd.classifying_ray = 1;
02890 
02891         /* initialize the lists of things that have been hit/missed */
02892         BU_LIST_INIT(&rd.rd_hit);
02893         BU_LIST_INIT(&rd.rd_miss);
02894 
02895         nmg_isect_ray_shell( &rd, s);
02896 #ifndef FAST_NMG
02897         while (BU_LIST_WHILE(a_hit, hitmiss, &rd.rd_miss)) {
02898                 BU_LIST_DEQUEUE( &a_hit->l );
02899                 bu_free( (char *)a_hit, "Miss list hitmiss struct" );
02900         }
02901 #else
02902         NMG_FREE_HITLIST( &rd.rd_miss, &ap );
02903 #endif
02904         /* count the number of hits */
02905         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT)) {
02906                 bu_log("%s[%d]: shell Hits:\n", __FILE__, __LINE__);
02907                 for( BU_LIST_FOR( a_hit, hitmiss, &rd.rd_hit ) )  {
02908                         if (a_hit->hit.hit_dist >= 0.0)
02909                                 bu_log("Positive dist hit\n");
02910                         else
02911                                 bu_log("Negative dist hit\n");
02912                         nmg_rt_print_hitmiss(a_hit);
02913                 }
02914         }
02915 
02916         if( (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT)) &&
02917             (rt_g.NMG_debug & (DEBUG_PLOTEM)) )  {
02918                 static int      num=0;
02919                 nmg_pl_hitmiss_list( "shell-ray", num++, &rd.rd_hit, rp );
02920         }
02921 
02922         minus_class = guess_class_from_hitlist_min(&rd, &hari_kari_minus, in_or_out_only);
02923         plus_class = guess_class_from_hitlist_max(&rd, &hari_kari_plus, in_or_out_only);
02924 
02925         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT)) {
02926                 bu_log("minus_class = (%s)\n", nmg_class_name(minus_class));
02927                 bu_log("plus_class = (%s)\n", nmg_class_name(plus_class));
02928         }
02929 
02930 
02931 #if 0
02932         /* XXX This should be fixed in the guess_* routines
02933          * instead of being fudged here.
02934          */
02935         if (hari_kari_minus) {
02936                 if(hari_kari_plus)
02937                         rt_bomb("double hari kari");
02938                 if (plus_class == NMG_CLASS_Unknown) {
02939                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02940                                 nmg_rt_print_hitlist(&rd.rd_hit);
02941                         rt_bomb("minus hari kari & plus unknown");
02942                 }
02943                 minus_class = plus_class;
02944 
02945         } else if (minus_class == NMG_CLASS_Unknown) {
02946                 if (plus_class == NMG_CLASS_Unknown) {
02947                         if (rt_g.NMG_debug & DEBUG_RT_ISECT)
02948                                 nmg_rt_print_hitlist(&rd.rd_hit);
02949                         rt_bomb("minus unknown & plus unknown");
02950                 }
02951                 minus_class = plus_class;
02952         } else if (plus_class == NMG_CLASS_Unknown || hari_kari_plus) {
02953                 plus_class = minus_class;
02954         }
02955 
02956 
02957         if (plus_class != minus_class) {
02958                 bu_log("%s:%d plus_class (%s) != minus_class (%s)\n",
02959                         __FILE__, __LINE__,
02960                         nmg_class_name(plus_class),
02961                         nmg_class_name(minus_class) );
02962 
02963                 nmg_rt_print_hitlist(&rd.rd_hit);
02964                 rt_bomb("");
02965         }
02966 #else
02967         /*
02968          *  Rather than blowing up, or guessing, just report that
02969          *  it didn't work, and let caller try another direction.
02970          */
02971         if( hari_kari_minus || hari_kari_plus )  {
02972                 if(rt_g.NMG_debug)
02973                         bu_log("hari_kari = %d, %d\n", hari_kari_minus, hari_kari_plus);
02974                 plus_class = NMG_CLASS_Unknown;
02975                 goto out;
02976         }
02977         if (plus_class != minus_class || plus_class == NMG_CLASS_Unknown ||
02978             minus_class == NMG_CLASS_Unknown ) {
02979 #if 0
02980                 nmg_rt_print_hitlist(&rd.rd_hit);
02981                 bu_log("minus_class = (%s) %d, hari_kari=%d\n", nmg_class_name(minus_class), minus_class, hari_kari_minus);
02982                 bu_log("plus_class = (%s)\n", nmg_class_name(plus_class));
02983                 bu_log("nmg_class_ray_vs_shell() -- can't tell\n");
02984 #endif
02985                 plus_class = NMG_CLASS_Unknown;
02986         }
02987 out:
02988 #endif
02989 
02990 #ifndef FAST_NMG
02991         while (BU_LIST_WHILE(a_hit, hitmiss, &rd.rd_hit)) {
02992                 BU_LIST_DEQUEUE( &a_hit->l );
02993                 bu_free( (char *)a_hit, "hit list hitmiss struct" );
02994         }
02995 #else
02996         NMG_FREE_HITLIST( &rd.rd_hit, &ap );
02997 #endif
02998 
02999         /* free the hitmiss freelist */
03000         if( !BU_LIST_UNINITIALIZED( &rt_uniresource.re_nmgfree ) )  {
03001                 struct hitmiss *hitp;
03002                 while( BU_LIST_WHILE( hitp, hitmiss, &rt_uniresource.re_nmgfree ) )  {
03003                         NMG_CK_HITMISS(hitp);
03004                         BU_LIST_DEQUEUE( (struct bu_list *)hitp );
03005                         bu_free( (genptr_t)hitp, "struct hitmiss" );
03006                 }
03007         }
03008 
03009         /* free the hitmiss table */
03010         bu_free( (char *)rd.hitmiss, "free nmg geom hit list");
03011 
03012         /* free the manifold table */
03013         bu_free( (char *)rd.manifolds, "free manifolds table" );
03014 
03015         if (rt_g.NMG_debug & (DEBUG_CLASSIFY|DEBUG_RT_ISECT))
03016                 bu_log("nmg_class_ray_vs_shell() returns %s(%d)\n",
03017                         nmg_class_name(plus_class), plus_class);
03018 
03019 
03020         return plus_class;
03021 }
03022 
03023 /*
03024  * Local Variables:
03025  * mode: C
03026  * tab-width: 8
03027  * c-basic-offset: 4
03028  * indent-tabs-mode: t
03029  * End:
03030  * ex: shiftwidth=4 tabstop=8
03031  */

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