BRL-CAD
nmg_pt_fu.c
Go to the documentation of this file.
1 /* N M G _ P T _ F U . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1994-2014 United States Government as represented by
5  * the U.S. Army Research Laboratory.
6  *
7  * This library is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public License
9  * version 2.1 as published by the Free Software Foundation.
10  *
11  * This library is distributed in the hope that it will be useful, but
12  * WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14  * Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public
17  * License along with this file; see the file named COPYING for more
18  * information.
19  */
20 /** @addtogroup nmg */
21 /** @{ */
22 /** @file primitives/nmg/nmg_pt_fu.c
23  *
24  * Routines for classifying a point with respect to a faceuse.
25  *
26  */
27 /** @} */
28 
29 #include "common.h"
30 
31 #include <stdlib.h>
32 #include <math.h>
33 #include "bio.h"
34 
35 #include "vmath.h"
36 #include "nmg.h"
37 #include "raytrace.h"
38 #include "plot3.h"
39 
40 
41 /* vertex/edge distance
42  * Each loop geometry element (edge/vertex) has one of these computed.
43  * We keep track of them for the whole face so that this value is computed
44  * only once per geometry element. This way we get consistent answers
45  */
46 struct ve_dist {
47  struct bu_list l;
48  uint32_t *magic_p; /* pointer to edge/vertex structure */
49  fastf_t dist; /* distance squared from point to edge */
50  struct vertex *v1;
51  struct vertex *v2;
52  int status; /* return code from bn_dist_pt3_lseg3 */
53 };
54 #define NMG_VE_DIST_MAGIC 0x102938
55 #define NMG_CK_VED(_p) NMG_CKMAG(_p, NMG_VE_DIST_MAGIC, "vertex/edge_dist")
56 
57 /* The loop builds a list of these things so that it can figure out the point
58  * classification based upon the classification of the pt vs each edge of
59  * the loop. This list is sorted in ascending "ved_p->dist" order.
60  */
61 struct edge_info {
62  struct bu_list l;
63  struct ve_dist *ved_p; /* ptr to ve_dist for this item */
64  struct edgeuse *eu_p; /* edgeuse pointer */
65  int nmg_class; /* pt classification WRT this item use */
66 };
67 #define NMG_EDGE_INFO_MAGIC 0xe100
68 #define NMG_CK_EI(_p) NMG_CKMAG(_p, NMG_EDGE_INFO_MAGIC, "edge_info")
69 
70 struct fpi {
71  uint32_t magic;
72  const struct bn_tol *tol;
73  const struct faceuse *fu_p;
74  struct bu_list ve_dh; /* ve_dist list head */
75  plane_t norm; /* surface normal for face(use) */
76  point_t pt; /* pt in plane of face to classify */
77  void (*eu_func)(struct edgeuse *, point_t, const char *); /* call w/eu when pt on edgeuse */
78  void (*vu_func)(struct vertexuse *, point_t, const char *); /* call w/vu when pt on vertexuse */
79  const char *priv; /* caller's private data */
80  int hits; /* flag PERUSE/PERGEOM */
81 };
82 #define NMG_FPI_MAGIC 12345678 /* fpi\0 */
83 #define NMG_CK_FPI(_fpi) \
84  NMG_CKMAG(_fpi, NMG_FPI_MAGIC, "fu_pt_info"); \
85  BN_CK_TOL(_fpi->tol); \
86  BU_CK_LIST_HEAD(&_fpi->ve_dh)
87 
88 #define NMG_FPI_TOUCHED 27
89 #define NMG_FPI_MISSED 32768
90 
91 static int nmg_class_pt_vu(struct fpi *fpi, struct vertexuse *vu);
92 static struct edge_info *nmg_class_pt_eu(struct fpi *fpi, struct edgeuse *eu, struct edge_info *edge_list, const int in_or_out_only);
93 static int compute_loop_class(struct fpi *fpi, const struct loopuse *lu, struct edge_info *edge_list);
94 static int nmg_class_pt_lu(struct loopuse *lu, struct fpi *fpi, const int in_or_out_only);
95 int nmg_class_pt_fu_except(const point_t pt, const struct faceuse *fu, const struct loopuse *ignore_lu, void (*eu_func)(struct edgeuse *, point_t, const char *), void (*vu_func)(struct vertexuse *, point_t, const char *), const char *priv, const int call_on_hits, const int in_or_out_only, const struct bn_tol *tol);
96 
97 
98 /**
99  * Find the square of the distance from a point P to a line segment described
100  * by the two endpoints A and B.
101  *
102  * P
103  * *
104  * /.
105  * / .
106  * / .
107  * / . (dist)
108  * / .
109  * / .
110  * *------*--------*
111  * A PCA B
112  *
113  * There are six distinct cases, with these return codes -
114  * 0 P is within tolerance of lseg AB. *dist = 0.
115  * 1 P is within tolerance of point A. *dist = 0.
116  * 2 P is within tolerance of point B. *dist = 0.
117  * 3 PCA is within tolerance of A. *dist = |P-A|**2.
118  * 4 PCA is within tolerance of B. *dist = |P-B|**2.
119  * 5 P is "above/below" lseg AB. *dist=|PCA-P|**2.
120  *
121  * This routine was formerly called bn_dist_pt_lseg().
122  * This is a special version that returns the square of the distance
123  * and does not actually calculate PCA.
124  *
125  */
126 int
127 bn_distsq_pt3_lseg3(fastf_t *dist, const fastf_t *a, const fastf_t *b, const fastf_t *p, const struct bn_tol *tol)
128 {
129  vect_t PtoA; /* P-A */
130  vect_t PtoB; /* P-B */
131  vect_t AtoB; /* B-A */
132  fastf_t P_A_sq; /* |P-A|**2 */
133  fastf_t P_B_sq; /* |P-B|**2 */
134  fastf_t B_A_sq; /* |B-A|**2 */
135  fastf_t t; /* distance squared along ray of projection of P */
136  fastf_t dot;
137 
138  BN_CK_TOL(tol);
139 
140  /* Check proximity to endpoint A */
141  VSUB2(PtoA, p, a);
142  P_A_sq = MAGSQ(PtoA);
143  if (P_A_sq < tol->dist_sq) {
144  /* P is within the tol->dist radius circle around A */
145  *dist = 0.0;
146  return 1;
147  }
148 
149  /* Check proximity to endpoint B */
150  VSUB2(PtoB, p, b);
151  P_B_sq = MAGSQ(PtoB);
152  if (P_B_sq < tol->dist_sq) {
153  /* P is within the tol->dist radius circle around B */
154  *dist = 0.0;
155  return 2;
156  }
157 
158  VSUB2(AtoB, b, a);
159  B_A_sq = MAGSQ(AtoB);
160 
161  /* compute distance squared (in actual units) along line to PROJECTION of
162  * point p onto the line: point pca
163  */
164  dot = VDOT(PtoA, AtoB);
165  t = dot * dot / B_A_sq;
166 
167  if (t < tol->dist_sq) {
168  /* PCA is at A */
169  *dist = P_A_sq;
170  return 3;
171  }
172 
173  if (dot < -SMALL_FASTF && t > tol->dist_sq) {
174  /* P is "left" of A */
175  *dist = P_A_sq;
176  return 3;
177  }
178  if (t < (B_A_sq - tol->dist_sq)) {
179  /* PCA falls between A and B */
180  register fastf_t dsq;
181 
182  /* Find distance from PCA to line segment (Pythagorus) */
183  dsq = P_A_sq - t;
184  if (dsq < tol->dist_sq) {
185  /* PCA is on lseg */
186  *dist = 0.0;
187  return 0;
188  }
189 
190  /* P is above or below lseg */
191  *dist = dsq;
192  return 5;
193  }
194 
195  /* P is "right" of B */
196  *dist = P_B_sq;
197  return 4;
198 }
199 
200 
201 /**
202  * Classify a point vs a vertex (touching/missed)
203  */
204 static int
205 nmg_class_pt_vu(struct fpi *fpi, struct vertexuse *vu)
206 {
207  vect_t delta;
208  struct ve_dist *ved;
209 
210  NMG_CK_VERTEXUSE(vu);
211 
212  /* see if we've classified this vertex WRT the point already */
213  for (BU_LIST_FOR(ved, ve_dist, &fpi->ve_dh)) {
214  NMG_CK_VED(ved);
215  if (ved->magic_p == &vu->v_p->magic) {
216  goto found;
217  }
218  }
219 
220  /* if we get here, we didn't find the vertex in the list of
221  * previously classified geometry. Create an entry in the
222  * face's list of processed geometry.
223  */
224  VSUB2(delta, vu->v_p->vg_p->coord, fpi->pt);
225 
226  BU_ALLOC(ved, struct ve_dist);
227  ved->magic_p = &vu->v_p->magic;
228  ved->dist = MAGNITUDE(delta);
229  if (ved->dist < fpi->tol->dist_sq) {
230  ved->status = NMG_FPI_TOUCHED;
231  if (fpi->hits == NMG_FPI_PERGEOM) {
232  /* need to cast vu_func pointer for actual use as a function */
233  void (*cfp)(struct vertexuse *, point_t, const char*);
234  cfp = (void (*)(struct vertexuse *, point_t, const char *))fpi->vu_func;
235  cfp(vu, fpi->pt, fpi->priv);
236  }
237  } else ved->status = NMG_FPI_MISSED;
238 
239  ved->v1 = ved->v2 = vu->v_p;
240 
242  BU_LIST_APPEND(&fpi->ve_dh, &ved->l);
243 found:
244 
245  if (fpi->vu_func &&
246  ved->status == NMG_FPI_TOUCHED &&
247  fpi->hits == NMG_FPI_PERUSE) {
248  /* need to cast vu_func pointer for actual use as a function */
249  void (*cfp)(struct vertexuse *, point_t, const char*);
250  cfp = (void (*)(struct vertexuse *, point_t, const char *))fpi->vu_func;
251  cfp(vu, fpi->pt, fpi->priv);
252  }
253 
254  return ved->status;
255 }
256 
257 
258 static int
259 Quadrant(fastf_t x, fastf_t y)
260 {
261  if (x > -SMALL_FASTF) {
262  if (y > -SMALL_FASTF)
263  return 1;
264  else
265  return 4;
266  } else {
267  if (y > -SMALL_FASTF)
268  return 2;
269  else
270  return 3;
271  }
272 }
273 
274 
275 int
276 nmg_eu_is_part_of_crack(const struct edgeuse *eu)
277 {
278  struct loopuse *lu;
279  struct edgeuse *eu_test;
280 
281  NMG_CK_EDGEUSE(eu);
282 
283  /* must be part of a loop to be a crack */
284  if (*eu->up.magic_p != NMG_LOOPUSE_MAGIC)
285  return 0;
286 
287  lu = eu->up.lu_p;
288  NMG_CK_LOOPUSE(lu);
289 
290  for (BU_LIST_FOR(eu_test, edgeuse, &lu->down_hd)) {
291  if (eu_test == eu)
292  continue;
293 
294  if (eu_test->vu_p->v_p == eu->eumate_p->vu_p->v_p &&
295  eu_test->eumate_p->vu_p->v_p == eu->vu_p->v_p)
296  return 1;
297  }
298 
299  return 0;
300 }
301 
302 
303 /**
304  * Classify a point with respect to an EU where the VU is the
305  * closest to the point. The EU and its left vector form an XY
306  * coordinate system in the face, with EU along the X-axis and
307  * its left vector along the Y-axis. Use these to decompose the
308  * direction of the prev_eu into X and Y coordinates (xo, yo) then
309  * do the same for the vector to the point to be classed (xpt, ypt).
310  * If (xpt, ypt) makes a smaller angle with EU than does (xo, yo),
311  * then PT is in the face, otherwise it is out.
312  *
313  *
314  * It is assumed that eu is from an OT_SAME faceuse, and that
315  * the PCA of PT to EU is at eu->vu_p.
316  */
317 
318 int
319 nmg_class_pt_euvu(const fastf_t *pt, struct edgeuse *eu_in, const struct bn_tol *tol)
320 {
321  struct loopuse *lu;
322  struct edgeuse *prev_eu;
323  struct edgeuse *eu;
324  struct vertex *v0, *v1, *v2;
325  vect_t left;
326  vect_t eu_dir;
327  vect_t other_eudir;
328  vect_t pt_dir;
329  fastf_t xo, yo;
330  fastf_t xpt, ypt;
331  fastf_t len;
332  int quado, quadpt;
333  int nmg_class = NMG_CLASS_Unknown;
334  int eu_is_crack = 0;
335  int prev_eu_is_crack = 0;
336 
337  NMG_CK_EDGEUSE(eu_in);
338  BN_CK_TOL(tol);
339 
340  eu = eu_in;
341 
342  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
343  bu_log("nmg_class_pt_euvu((%g %g %g), eu=%p)\n", V3ARGS(pt), (void *)eu);
344 
345  if (UNLIKELY(*eu->up.magic_p != NMG_LOOPUSE_MAGIC)) {
346  bu_log("nmg_class_pt_euvu() called with eu (%p) that isn't part of a loop\n", (void *)eu);
347  bu_bomb("nmg_class_pt_euvu() called with eu that isn't part of a loop");
348  }
349  lu = eu->up.lu_p;
350  NMG_CK_LOOPUSE(lu);
351 
352  eu_is_crack = nmg_eu_is_part_of_crack(eu);
353 
354  prev_eu = BU_LIST_PPREV_CIRC(edgeuse, &eu->l);
355 
356  prev_eu_is_crack = nmg_eu_is_part_of_crack(prev_eu);
357 
358  /* if both EU's are cracks, we cannot classify */
359  if (eu_is_crack && prev_eu_is_crack)
360  return NMG_CLASS_Unknown;
361 
362  if (eu_is_crack) {
363  struct edgeuse *eu_test;
364  int done = 0;
365 
366  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
367  bu_log("nmg_class_pt_euvu: eu %p is a crack\n", (void *)eu);
368 
369  /* find next eu from this vertex that is not a crack */
370  eu_test = BU_LIST_PNEXT_CIRC(edgeuse, &eu->l);
371  while (!done) {
372  while (eu_test->vu_p->v_p != eu->vu_p->v_p && eu_test != eu)
373  eu_test = BU_LIST_PNEXT_CIRC(edgeuse, &eu_test->l);
374 
375  if (eu_test == eu)
376  done = 1;
377 
378  if (!nmg_eu_is_part_of_crack(eu_test))
379  done = 1;
380 
381  if (!done)
382  eu_test = BU_LIST_PNEXT_CIRC(edgeuse, &eu_test->l);
383  }
384 
385  if (eu_test == eu) /* can't get away from crack */
386  return NMG_CLASS_Unknown;
387  else
388  eu = eu_test;
389 
390  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
391  bu_log("\tUsing eu %p instead\n", (void *)eu);
392  }
393 
394  if (prev_eu_is_crack) {
395  struct edgeuse *eu_test;
396  int done = 0;
397 
398  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
399  bu_log("nmg_class_pt_euvu: prev_eu (%p) is a crack\n", (void *)prev_eu);
400 
401  /* find previous eu ending at this vertex that is not a crack */
402  eu_test = BU_LIST_PPREV_CIRC(edgeuse, &prev_eu->l);
403  while (!done) {
404  while (eu_test->eumate_p->vu_p->v_p != eu->vu_p->v_p && eu_test != prev_eu)
405  eu_test = BU_LIST_PPREV_CIRC(edgeuse, &eu_test->l);
406 
407  if (eu_test == prev_eu)
408  done = 1;
409 
410  if (!nmg_eu_is_part_of_crack(eu_test))
411  done = 1;
412 
413  if (!done)
414  eu_test = BU_LIST_PPREV_CIRC(edgeuse, &eu_test->l);
415  }
416 
417  if (eu_test == prev_eu) /* can't get away from crack */
418  return NMG_CLASS_Unknown;
419  else
420  prev_eu = eu_test;
421 
422  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
423  bu_log("\tUsing prev_eu %p instead\n", (void *)prev_eu);
424  }
425 
426  /* left is the Y-axis of our XY-coordinate system */
427  if (UNLIKELY(nmg_find_eu_leftvec(left, eu))) {
428  bu_log("nmg_class_pt_euvu: nmg_find_eu_leftvec() for eu=%p failed!\n", (void *)eu);
429  bu_bomb("nmg_class_pt_euvu: nmg_find_eu_leftvec() failed!");
430  }
431 
432  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
433  bu_log("\tprev_eu = %p, left = (%g %g %g)\n", (void *)prev_eu, V3ARGS(left));
434 
435  /* v0 is the origin of the XY-coordinate system */
436  v0 = eu->vu_p->v_p;
437  NMG_CK_VERTEX(v0);
438 
439  /* v1 is on the X-axis */
440  v1 = eu->eumate_p->vu_p->v_p;
441  NMG_CK_VERTEX(v1);
442 
443  /* v2 determines angle prev_eu makes with X-axis */
444  v2 = prev_eu->vu_p->v_p;
445  NMG_CK_VERTEX(v2);
446 
447  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
448  bu_log("\tv0=%p, v1=%p, v2=%p\n", (void *)v0, (void *)v1, (void *)v2);
449 
450  /* eu_dir is our X-direction */
451  VSUB2(eu_dir, v1->vg_p->coord, v0->vg_p->coord);
452 
453  /* other_eudir is direction along the previous EU (from origin) */
454  VSUB2(other_eudir, v2->vg_p->coord, v0->vg_p->coord);
455 
456  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
457  bu_log("\teu_dir=(%g %g %g), other_eudir=(%g %g %g)\n", V3ARGS(eu_dir), V3ARGS(other_eudir));
458 
459  /* get X and Y components for other_eu */
460  xo = VDOT(eu_dir, other_eudir);
461  yo = VDOT(left, other_eudir);
462 
463  /* which quadrant does this XY point lie in */
464  quado = Quadrant(xo, yo);
465 
466  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
467  bu_log("\txo=%g, yo=%g, quadrant=%d\n", xo, yo, quado);
468 
469  /* get direction to PT from origin */
470  VSUB2(pt_dir, pt, v0->vg_p->coord);
471 
472  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
473  bu_log("\tpt_dir=(%g %g %g)\n", V3ARGS(pt_dir));
474 
475  /* get X and Y components for PT */
476  xpt = VDOT(eu_dir, pt_dir);
477  ypt = VDOT(left, pt_dir);
478 
479  /* which quadrant does this XY point lie in */
480  quadpt = Quadrant(xpt, ypt);
481 
482  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
483  bu_log("\txpt=%g, ypt=%g, quadrant=%d\n", xpt, ypt, quadpt);
484 
485  /* do a quadrant comparison first (cheap!!!) */
486  if (quadpt < quado)
487  return NMG_CLASS_AinB;
488 
489  if (quadpt > quado)
490  return NMG_CLASS_AoutB;
491 
492  /* both are in the same quadrant, need to normalize the coordinates */
493  len = sqrt(xo*xo + yo*yo);
494  xo = xo/len;
495  yo = yo/len;
496 
497  len = sqrt(xpt*xpt + ypt*ypt);
498  xpt = xpt/len;
499  ypt = ypt/len;
500 
501  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
502  bu_log("\tNormalized xo, yo=(%g %g), xpt, ypt=(%g %g)\n", xo, yo, xpt, ypt);
503 
504  switch (quadpt) {
505  case 1:
506  if (xpt >= xo && ypt <= yo)
507  nmg_class = NMG_CLASS_AinB;
508  else
509  nmg_class = NMG_CLASS_AoutB;
510  break;
511  case 2:
512  if (xpt >= xo && ypt >= yo)
513  nmg_class = NMG_CLASS_AinB;
514  else
515  nmg_class = NMG_CLASS_AoutB;
516  break;
517  case 3:
518  if (xpt <= xo && ypt >= yo)
519  nmg_class = NMG_CLASS_AinB;
520  else
521  nmg_class = NMG_CLASS_AoutB;
522  break;
523  case 4:
524  if (xpt <= xo && ypt <= yo)
525  nmg_class = NMG_CLASS_AinB;
526  else
527  nmg_class = NMG_CLASS_AoutB;
528  break;
529  default:
530  bu_log("This can't happen (illegal quadrant %d)\n", quadpt);
531  bu_bomb("This can't happen (illegal quadrant)\n");
532  break;
533  }
534  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
535  bu_log("returning %s\n", nmg_class_name(nmg_class));
536 
537  return nmg_class;
538 }
539 
540 
541 /**
542  * If there is no ve_dist structure for this edge, compute one and
543  * add it to the list.
544  *
545  * Sort an edge_info structure into the loops list of edgeuse status
546  */
547 static struct edge_info *
548 nmg_class_pt_eu(struct fpi *fpi, struct edgeuse *eu, struct edge_info *edge_list, const int in_or_out_only)
549 {
550  struct bn_tol tmp_tol;
551  struct edgeuse *next_eu;
552  struct ve_dist *ved, *ed;
553  struct edge_info *ei_p;
554  struct edge_info *ei;
555  pointp_t eu_pt;
556  vect_t left;
557  vect_t v_to_pt;
558  int found_data = 0;
559 
560  NMG_CK_FPI(fpi);
561  BN_CK_TOL(fpi->tol);
562 
563  if (RTG.NMG_debug & DEBUG_PT_FU) {
564  bu_log("pt (%g %g %g) vs_edge (%g %g %g) (%g %g %g) (eu=%p)\n",
565  V3ARGS(fpi->pt),
566  V3ARGS(eu->vu_p->v_p->vg_p->coord),
567  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord), (void *)eu);
568  }
569 
570  /* we didn't find a ve_dist structure for this edge, so we'll
571  * have to do the calculations.
572  */
573  tmp_tol = (*fpi->tol);
574  if (in_or_out_only) {
575  tmp_tol.dist = 0.0;
576  tmp_tol.dist_sq = 0.0;
577  }
578 
579  BU_ALLOC(ved, struct ve_dist);
580  ved->magic_p = &eu->e_p->magic;
581  ved->status = bn_distsq_pt3_lseg3(&ved->dist,
582  eu->vu_p->v_p->vg_p->coord,
583  eu->eumate_p->vu_p->v_p->vg_p->coord,
584  fpi->pt,
585  &tmp_tol);
586  ved->v1 = eu->vu_p->v_p;
587  ved->v2 = eu->eumate_p->vu_p->v_p;
589  BU_LIST_APPEND(&fpi->ve_dh, &ved->l);
590  eu_pt = ved->v1->vg_p->coord;
591 
592  if (RTG.NMG_debug & DEBUG_PT_FU) {
593  bu_log("nmg_class_pt_eu: status for eu %p (%g %g %g)<->(%g %g %g) vs pt (%g %g %g) is %d\n",
594  (void *)eu, V3ARGS(eu->vu_p->v_p->vg_p->coord),
595  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord),
596  V3ARGS(fpi->pt), ved->status);
597  bu_log("\tdist = %g\n", ved->dist);
598  }
599 
600 
601  /* Add a struct for this edgeuse to the loop's list of dist-sorted
602  * edgeuses.
603  */
604  BU_ALLOC(ei, struct edge_info);
605  ei->ved_p = ved;
606  ei->eu_p = eu;
608 
609  /* compute the status (ei->status) of the point WRT this edge */
610 
611  switch (ved->status) {
612  case 0: /* pt is on the edge(use) */
613  ei->nmg_class = NMG_CLASS_AonBshared;
614  if (fpi->eu_func &&
615  (fpi->hits == NMG_FPI_PERUSE ||
616  (fpi->hits == NMG_FPI_PERGEOM && !found_data))) {
617  /* need to cast eu_func pointer for actual use as a function */
618  void (*cfp)(struct edgeuse *, point_t, const char*);
619  cfp = (void (*)(struct edgeuse *, point_t, const char *))fpi->eu_func;
620  cfp(eu, fpi->pt, fpi->priv);
621  }
622  break;
623  case 1: /* within tolerance of endpoint at ved->v1 */
624  ei->nmg_class = NMG_CLASS_AonBshared;
625  /* add an entry for the vertex in the edge list so that
626  * other uses of this vertex will claim the point is within
627  * tolerance without re-computing
628  */
629  BU_ALLOC(ed, struct ve_dist);
630  ed->magic_p = &ved->v1->magic;
631  ed->status = ved->status;
632  ed->v1 = ed->v2 = ved->v1;
633 
635  BU_LIST_APPEND(&fpi->ve_dh, &ed->l);
636 
637  if (fpi->vu_func &&
638  (fpi->hits == NMG_FPI_PERUSE ||
639  (fpi->hits == NMG_FPI_PERGEOM && !found_data))) {
640  /* need to cast vu_func pointer for actual use as a function */
641  void (*cfp)(struct vertexuse *, point_t, const char*);
642  cfp = (void (*)(struct vertexuse *, point_t, const char *))fpi->vu_func;
643  cfp(eu->vu_p, fpi->pt, fpi->priv);
644  }
645 
646  break;
647  case 2: /* within tolerance of endpoint at ved->v2 */
648  ei->nmg_class = NMG_CLASS_AonBshared;
649  /* add an entry for the vertex in the edge list so that
650  * other uses of this vertex will claim the point is within
651  * tolerance without re-computing
652  */
653  BU_ALLOC(ed, struct ve_dist);
654  ed->magic_p = &ved->v2->magic;
655  ed->status = ved->status;
656  ed->v1 = ed->v2 = ved->v2;
657 
659  BU_LIST_APPEND(&fpi->ve_dh, &ed->l);
660  if (fpi->vu_func &&
661  (fpi->hits == NMG_FPI_PERUSE ||
662  (fpi->hits == NMG_FPI_PERGEOM && !found_data))) {
663  /* need to cast vu_func pointer for actual use as a function */
664  void (*cfp)(struct vertexuse *, point_t, const char*);
665  cfp = (void (*)(struct vertexuse *, point_t, const char *))fpi->vu_func;
666  cfp(eu->eumate_p->vu_p, fpi->pt, fpi->priv);
667  }
668  break;
669 
670  case 3: /* PCA of pt on line is within tolerance of ved->v1 of segment */
671  ei->nmg_class = nmg_class_pt_euvu(fpi->pt, eu, fpi->tol);
672  if (ei->nmg_class == NMG_CLASS_Unknown)
673  ei->ved_p->dist = MAX_FASTF;
674  break;
675  case 4: /* PCA of pt on line is within tolerance of ved->v2 of segment */
676  next_eu = BU_LIST_PNEXT_CIRC(edgeuse, &eu->l);
677  ei->nmg_class = nmg_class_pt_euvu(fpi->pt, next_eu, fpi->tol);
678  if (ei->nmg_class == NMG_CLASS_Unknown)
679  ei->ved_p->dist = MAX_FASTF;
680  break;
681 
682  case 5: /* PCA is along length of edge, but point is NOT on edge. */
683  if (nmg_find_eu_left_non_unit(left, eu))
684  bu_bomb("can't find left vector\n");
685  /* take dot product of v->pt vector with left to determine
686  * if pt is inside/left of edge
687  */
688  VSUB2(v_to_pt, fpi->pt, eu_pt);
689  if (VDOT(v_to_pt, left) > -SMALL_FASTF)
690  ei->nmg_class = NMG_CLASS_AinB;
691  else
692  ei->nmg_class = NMG_CLASS_AoutB;
693  break;
694  default:
695  bu_log("%s:%d status = %d\n", __FILE__, __LINE__, ved->status);
696  bu_bomb("Why did this happen?");
697  break;
698  }
699 
700 
701  if (RTG.NMG_debug & DEBUG_PT_FU) {
702  bu_log("pt @ dist %g from edge classed %s vs edge\n",
703  ei->ved_p->dist, nmg_class_name(ei->nmg_class));
704 /* pl_pt_e(fpi, ei); */
705  }
706 
707  /* now that it's complete, add ei to the edge list */
708  for (BU_LIST_FOR(ei_p, edge_info, &edge_list->l)) {
709  /* if the distance to this edge is smaller, or
710  * if the distance is the same & the edge is the same
711  * Insert edge_info struct here in list
712  */
713  if (ved->dist < ei_p->ved_p->dist
714  || (ZERO(ved->dist - ei_p->ved_p->dist)
715  && ei_p->ved_p->magic_p == ved->magic_p))
716  {
717  break;
718  }
719  }
720 
721  BU_LIST_INSERT(&ei_p->l, &ei->l);
722  return ei;
723 }
724 
725 
726 /**
727  * Make a list of all edgeuses which are at the same distance as the
728  * first element on the list. Toss out opposing pairs of edgeuses of the
729  * same edge.
730  *
731  */
732 HIDDEN void make_near_list(struct edge_info *edge_list, struct bu_list *near1, const struct bn_tol *tol)
733 {
734  struct edge_info *ei;
735  struct edge_info *ei_p;
736  struct edge_info *tmp;
737  fastf_t dist;
738 
739  BU_CK_LIST_HEAD(&edge_list->l);
740  BU_CK_LIST_HEAD(near1);
741 
742  /* toss opposing pairs of uses of the same edge from the list */
743  ei = BU_LIST_FIRST(edge_info, &edge_list->l);
744  while (BU_LIST_NOT_HEAD(&ei->l, &edge_list->l)) {
745  NMG_CK_EI(ei);
746  ei_p = BU_LIST_FIRST(edge_info, &edge_list->l);
747  while (BU_LIST_NOT_HEAD(&ei_p->l, &edge_list->l)) {
748  NMG_CK_EI(ei_p);
749  NMG_CK_VED(ei_p->ved_p);
750 
751  /* if we've found an opposing use of the same
752  * edge toss the pair of them
753  */
754  if (ei_p->ved_p->magic_p == ei->ved_p->magic_p &&
755  ei_p->eu_p->eumate_p->vu_p->v_p == ei->eu_p->vu_p->v_p &&
756  ei_p->eu_p->vu_p->v_p == ei->eu_p->eumate_p->vu_p->v_p) {
757  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
758  bu_log("tossing edgeuse pair:\n");
759  bu_log("(%g %g %g) -> (%g %g %g)\n",
760  V3ARGS(ei->eu_p->vu_p->v_p->vg_p->coord),
761  V3ARGS(ei->eu_p->eumate_p->vu_p->v_p->vg_p->coord));
762  bu_log("(%g %g %g) -> (%g %g %g)\n",
763  V3ARGS(ei_p->eu_p->vu_p->v_p->vg_p->coord),
764  V3ARGS(ei_p->eu_p->eumate_p->vu_p->v_p->vg_p->coord));
765  }
766 
767  tmp = ei_p;
768  ei_p = BU_LIST_PLAST(edge_info, &ei_p->l);
769  BU_LIST_DEQUEUE(&tmp->l);
770  bu_free((char *)tmp, "edge info struct");
771  tmp = ei;
772  ei = BU_LIST_PLAST(edge_info, &ei->l);
773  BU_LIST_DEQUEUE(&tmp->l);
774  bu_free((char *)tmp, "edge info struct");
775  break;
776  }
777  ei_p = BU_LIST_PNEXT(edge_info, &ei_p->l);
778  }
779  ei = BU_LIST_PNEXT(edge_info, &ei->l);
780  }
781 
782  if (BU_LIST_IS_EMPTY(&edge_list->l))
783  return;
784 
785  ei = BU_LIST_FIRST(edge_info, &edge_list->l);
786  NMG_CK_EI(ei);
787  NMG_CK_VED(ei->ved_p);
788  dist = ei->ved_p->dist;
789 
790  /* create "near" list with all ei's at this dist */
791  for (BU_LIST_FOR(ei, edge_info, &edge_list->l)) {
792  NMG_CK_EI(ei);
793  NMG_CK_VED(ei->ved_p);
794  if (NEAR_EQUAL(ei->ved_p->dist, dist, tol->dist_sq)) {
795  ei_p = BU_LIST_PLAST(edge_info, &ei->l);
796  BU_LIST_DEQUEUE(&ei->l);
797  BU_LIST_APPEND(near1, &ei->l);
798  ei = ei_p;
799  }
800  }
801 
802  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
803  bu_log("dist %g near list\n", dist);
804  for (BU_LIST_FOR(ei, edge_info, near1)) {
805  bu_log("\t(%g %g %g) -> (%g %g %g)\n",
806  V3ARGS(ei->eu_p->vu_p->v_p->vg_p->coord),
807  V3ARGS(ei->eu_p->eumate_p->vu_p->v_p->vg_p->coord));
808  bu_log("\tdist:%g class:%s status:%d\n\t\tv1(%g %g %g) v2(%g %g %g)\n",
809  ei->ved_p->dist, nmg_class_name(ei->nmg_class),
810  ei->ved_p->status,
811  V3ARGS(ei->ved_p->v1->vg_p->coord),
812  V3ARGS(ei->ved_p->v2->vg_p->coord));
813  bu_log("\tei->ved_p->magic_p=%p, ei->eu_p->vu_p=%p, ei->eu_p->eumate_p->vu_p=%p\n",
814  (void *)ei->ved_p->magic_p, (void *)ei->eu_p->vu_p, (void *)ei->eu_p->eumate_p->vu_p);
815  }
816  }
817 }
818 
819 
820 static void
821 pl_pt_lu(struct fpi *fpi, const struct loopuse *lu, struct edge_info *ei)
822 {
823  FILE *fp;
824  char name[25];
825  long *b;
826  static int plot_file_number = 0;
827  int i;
828  point_t p1, p2;
829  point_t pca;
830  fastf_t dist;
831  struct bn_tol tmp_tol;
832 
833  NMG_CK_FPI(fpi);
834  NMG_CK_FACEUSE(fpi->fu_p);
835  NMG_CK_LOOPUSE(lu);
836  NMG_CK_EI(ei);
837 
838  sprintf(name, "pt_lu%02d.plot3", plot_file_number++);
839  fp = fopen(name, "wb");
840  if (fp == (FILE *)NULL) {
841  perror(name);
842  bu_bomb("unable to open file for writing");
843  }
844 
845  bu_log("\toverlay %s\n", name);
846  b = (long *)bu_calloc(fpi->fu_p->s_p->r_p->m_p->maxindex,
847  sizeof(long), "bit vec"),
848 
849  pl_erase(fp);
850  pd_3space(fp,
851  fpi->fu_p->f_p->min_pt[0]-1.0,
852  fpi->fu_p->f_p->min_pt[1]-1.0,
853  fpi->fu_p->f_p->min_pt[2]-1.0,
854  fpi->fu_p->f_p->max_pt[0]+1.0,
855  fpi->fu_p->f_p->max_pt[1]+1.0,
856  fpi->fu_p->f_p->max_pt[2]+1.0);
857 
858  nmg_pl_lu(fp, lu, b, 255, 255, 255);
859 
860  tmp_tol.magic = BN_TOL_MAGIC;
861  tmp_tol.dist = 0.0005;
862  tmp_tol.dist_sq = tmp_tol.dist * tmp_tol.dist;
863  tmp_tol.perp = 1e-6;
864  tmp_tol.para = 1 - tmp_tol.perp;
865 
866  (void)bn_dist_pt3_lseg3(&dist, pca, ei->eu_p->vu_p->v_p->vg_p->coord,
867  ei->eu_p->eumate_p->vu_p->v_p->vg_p->coord, fpi->pt, &tmp_tol);
868 
869  pl_color(fp, 255, 255, 50);
870  pdv_3line(fp, pca, fpi->pt);
871 
872  pl_color(fp, 255, 64, 255);
873 
874  /* make a nice axis-cross at the point in question */
875  for (i = 0; i < 3; i++) {
876  VMOVE(p1, fpi->pt);
877  p1[i] -= 1.0;
878  VMOVE(p2, fpi->pt);
879  p2[i] += 1.0;
880  pdv_3line(fp, p1, p2);
881  }
882 
883  bu_free((char *)b, "bit vec");
884  fclose(fp);
885 }
886 
887 
888 /**
889  * Given a list of edge_info structures for the edges of a loop,
890  * determine what the classification for the loop should be.
891  *
892  * If passed a "crack" loop, will produce random results.
893  */
894 
895 static int
896 compute_loop_class(struct fpi *fpi,
897  const struct loopuse *lu,
898  struct edge_info *edge_list)
899 {
900  struct edge_info *ei;
901  struct bu_list near1;
902  int lu_class = NMG_CLASS_Unknown;
903 
904  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
905  bu_log("compute_loop_class()\n");
906  for (BU_LIST_FOR(ei, edge_info, &edge_list->l)) {
907  bu_log("dist:%g class:%s status:%d\n\tv1(%g %g %g) v2(%g %g %g)\n",
908  ei->ved_p->dist, nmg_class_name(ei->nmg_class),
909  ei->ved_p->status,
910  V3ARGS(ei->ved_p->v1->vg_p->coord),
911  V3ARGS(ei->ved_p->v2->vg_p->coord));
912  }
913  }
914 
915  BU_CK_LIST_HEAD(&edge_list->l);
916  BU_LIST_INIT(&near1);
917 
918  /* get a list of "closest/useful" edges to use in classifying
919  * the pt WRT the loop
920  */
921  while (BU_LIST_IS_EMPTY(&near1) && BU_LIST_NON_EMPTY(&edge_list->l)) {
922  make_near_list(edge_list, &near1, fpi->tol);
923  }
924 
925  if (BU_LIST_IS_EMPTY(&near1)) {
926  /* This was a "crack" or "snake" loop */
927 
928  if (lu->orientation == OT_SAME) {
929  lu_class = NMG_CLASS_AoutB;
930  } else if (lu->orientation == OT_OPPOSITE) {
931  lu_class = NMG_CLASS_AinB;
932  } else
933  bu_bomb("bad lu orientation\n");
934 
935  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
936  bu_log("list was empty, so class is %s\n",
937  nmg_class_name(lu_class));
938  }
939 
940  return lu_class;
941  }
942 
943  for (BU_LIST_FOR(ei, edge_info, &near1)) {
944  int done = 0;
945  NMG_CK_EI(ei);
946  NMG_CK_VED(ei->ved_p);
947  switch (ei->ved_p->status) {
948  case 0: /* pt is on edge */
949  case 1: /* pt is on ei->ved_p->v1 */
950  case 2: /* pt is on ei->ved_p->v2 */
951  lu_class = NMG_CLASS_AonBshared;
952  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
953  pl_pt_lu(fpi, lu, ei);
954  done = 1;
955  break;
956  case 3: /* pt pca is v1 */
957  case 4: /* pt pca is v2 */
958  case 5: /* pt pca between v1 and v2 */
959  lu_class = ei->nmg_class;
960  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
961  bu_log("found status 5 edge, loop class is %s\n",
962  nmg_class_name(lu_class));
963  }
964  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU))
965  pl_pt_lu(fpi, lu, ei);
966  done = 1;
967  break;
968  default:
969  bu_log("%s:%d status = %d\n",
970  __FILE__, __LINE__, ei->ved_p->status);
971  bu_bomb("Why did this happen?");
972  break;
973  }
974  if (done) {
975  break;
976  }
977  }
978 
979  /* the caller will get whatever is left of the edge_list, but
980  * we need to free up the "near" list
981  */
982  while (BU_LIST_WHILE(ei, edge_info, &near1)) {
983  BU_LIST_DEQUEUE(&ei->l);
984  bu_free((char *)ei, "edge_info struct");
985  }
986 
987  if (UNLIKELY(RTG.NMG_debug & DEBUG_PT_FU)) {
988  bu_log("compute_loop_class() returns %s\n",
989  nmg_class_name(lu_class));
990  }
991 
992  return lu_class;
993 }
994 /**
995  *
996  * For each edgeuse, compute an edge_info structure.
997  *
998  * if min_dist == 0 return ON
999  *
1000  * For dist min -> max
1001  * if even # of uses of edge in this loop
1002  * "ignore" all uses of this edge since we can't answer the
1003  * "spike" problem: *---------*
1004  * | . |
1005  * *---* | .
1006  * | *----*
1007  * | |
1008  * *---------*
1009  * else (odd # of uses of edge in this loop)
1010  * "ignore" consecutive uses of the same edge to avoid the
1011  * "accordian pleat" problem: *-------*
1012  * | . |
1013  * *----* |
1014  * *----* |
1015  * *----* |
1016  * *----*--*
1017  * classify pt WRT remaining edgeuse
1018  *
1019  * The "C-clamp" problem:
1020  * *---------------*
1021  * | |
1022  * | *----* |
1023  * | | | . |
1024  * | | *-------*
1025  * | | | |
1026  * | *----* |
1027  * | |
1028  * *---------------*
1029  *
1030  */
1031 static int
1032 nmg_class_pt_lu(struct loopuse *lu, struct fpi *fpi, const int in_or_out_only)
1033 {
1034  int lu_class = NMG_CLASS_Unknown;
1035 
1036  NMG_CK_FPI(fpi);
1037  NMG_CK_LOOPUSE(lu);
1038  NMG_CK_LOOP(lu->l_p);
1039  NMG_CK_LOOP_G(lu->l_p->lg_p);
1040 
1041 
1042  if (V3PT_OUT_RPP_TOL(fpi->pt, lu->l_p->lg_p->min_pt, lu->l_p->lg_p->max_pt, fpi->tol->dist)) {
1043  /* point is not in RPP of loop */
1044 
1045  if (RTG.NMG_debug & DEBUG_PT_FU) {
1046  bu_log("nmg_class_pt_lu(pt(%g %g %g) outside loop RPP\n",
1047  V3ARGS(fpi->pt));
1048  bu_log(" lu RPP: (%g %g %g) <-> (%g %g %g)\n",
1049  V3ARGS(lu->l_p->lg_p->min_pt), V3ARGS(lu->l_p->lg_p->max_pt));
1050  }
1051 
1052  if (lu->orientation == OT_SAME)
1053  return NMG_CLASS_AoutB;
1054  else if (lu->orientation == OT_OPPOSITE)
1055  return NMG_CLASS_AinB;
1056  else if (lu->orientation == OT_UNSPEC)
1057  return NMG_CLASS_Unknown;
1058 
1059  }
1060 
1061  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_EDGEUSE_MAGIC) {
1062  register struct edgeuse *eu;
1063  struct edge_info edge_list;
1064  struct edge_info *ei;
1065  int is_crack;
1066 
1067  is_crack = nmg_loop_is_a_crack(lu);
1068  if (lu->orientation == OT_OPPOSITE && is_crack) {
1069  /* Even if point lies on a crack, if it's an
1070  * OT_OPPOSITE crack loop, it subtracts nothing.
1071  * Just ignore it.
1072  */
1073  lu_class = NMG_CLASS_AinB;
1074  if (RTG.NMG_debug & DEBUG_PT_FU)
1075  bu_log("nmg_class_pt_lu() ignoring OT_OPPOSITE crack loop\n");
1076  goto out;
1077  }
1078 
1079  BU_LIST_INIT(&edge_list.l);
1080 
1081  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1082  ei = nmg_class_pt_eu(fpi, eu, &edge_list, in_or_out_only);
1083  NMG_CK_EI(ei);
1084  NMG_CK_VED(ei->ved_p);
1085  if (!in_or_out_only && ei->ved_p->dist < fpi->tol->dist_sq) {
1086  lu_class = NMG_CLASS_AinB;
1087  break;
1088  }
1089  }
1090  /* */
1091  if (lu_class == NMG_CLASS_Unknown) {
1092  /* pt does not touch any edge or vertex */
1093  if (is_crack) {
1094  /* orientation here is OT_SAME */
1095  lu_class = NMG_CLASS_AoutB;
1096  } else {
1097  lu_class = compute_loop_class(fpi, lu, &edge_list);
1098  }
1099  } else {
1100  /* pt touches edge or vertex */
1101  if (RTG.NMG_debug & DEBUG_PT_FU)
1102  bu_log("loop class already known (pt must touch edge)\n");
1103  }
1104 
1105  /* free up the edge_list elements */
1106  while (BU_LIST_WHILE(ei, edge_info, &edge_list.l)) {
1107  BU_LIST_DEQUEUE(&ei->l);
1108  bu_free((char *)ei, "edge info struct");
1109  }
1110  } else if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
1111  register struct vertexuse *vu;
1112  int v_class;
1113 
1114  vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
1115  v_class = nmg_class_pt_vu(fpi, vu);
1116 
1117  switch (lu->orientation) {
1118  case OT_BOOLPLACE:
1119  case OT_SAME:
1120  if (v_class == NMG_FPI_TOUCHED)
1121  lu_class = NMG_CLASS_AinB;
1122  else
1123  lu_class = NMG_CLASS_AoutB;
1124  break;
1125  case OT_OPPOSITE:
1126  /* Why even call nmg_class_pt_vu() here, if return isn't used? */
1127  lu_class = NMG_CLASS_AinB;
1128  break;
1129  case OT_UNSPEC:
1130  lu_class = NMG_CLASS_Unknown;
1131  break;
1132  default:
1133  bu_log("nmg_class_pt_lu() hit %s loop at vu=%p\n",
1134  nmg_orientation(lu->orientation), (void *)vu);
1135  bu_bomb("nmg_class_pt_lu() Loop orientation error\n");
1136  break;
1137  }
1138  } else {
1139  bu_log("%s:%d bad child of loopuse\n", __FILE__, __LINE__);
1140  bu_bomb("nmg_class_pt_lu() crash and burn\n");
1141  }
1142 
1143 
1144 out:
1145  if (RTG.NMG_debug & DEBUG_PT_FU)
1146  bu_log("nmg_class_pt_lu() pt classed %s vs loop\n", nmg_class_name(lu_class));
1147 
1148  return lu_class;
1149 }
1150 
1151 
1152 static void
1153 plot_parity_error(const struct faceuse *fu, const fastf_t *pt)
1154 {
1155  long *b;
1156  FILE *fp;
1157  point_t p1, p2;
1158  int i;
1159 
1160  NMG_CK_FACEUSE(fu);
1161 
1162  fp=fopen("pt_fu_parity_error.plot3", "wb");
1163  if (!fp)
1164  bu_bomb("error opening pt_fu_parity_error.plot3\n");
1165 
1166  bu_log("overlay pt_fu_parity_error.plot3\n");
1167 
1168  b = (long *)bu_calloc(fu->s_p->r_p->m_p->maxindex,
1169  sizeof(long), "bit vec"),
1170 
1171  pl_erase(fp);
1172  pd_3space(fp,
1173  fu->f_p->min_pt[0]-1.0,
1174  fu->f_p->min_pt[1]-1.0,
1175  fu->f_p->min_pt[2]-1.0,
1176  fu->f_p->max_pt[0]+1.0,
1177  fu->f_p->max_pt[1]+1.0,
1178  fu->f_p->max_pt[2]+1.0);
1179 
1180  nmg_pl_fu(fp, fu, b, 200, 200, 200);
1181 
1182 
1183  /* make a nice axis-cross at the point in question */
1184  for (i = 0; i < 3; i++) {
1185  VMOVE(p1, pt);
1186  p1[i] -= 1.0;
1187  VMOVE(p2, pt);
1188  p2[i] += 1.0;
1189  pdv_3line(fp, p1, p2);
1190  }
1191 
1192  bu_free((char *)b, "plot table");
1193  fclose(fp);
1194 
1195 }
1196 
1197 
1198 /**
1199  *
1200  * Classify a point on a face's plane as being inside/outside the area
1201  * of the face.
1202  *
1203  * For each loopuse, compute IN/ON/OUT
1204  *
1205  * if any loop has pt classified as "ON" return "ON" (actually returns "IN" -jra)
1206  *
1207  * ignore all OT_SAME loops w/pt classified as "OUT"
1208  * ignore all OT_OPPOSITE loops w/pt classified as "IN"
1209  * If (# OT_SAME loops == # OT_OPPOSITE loops)
1210  * pt is "OUT"
1211  * else if (# OT_SAME loops - 1 == # OT_OPPOSITE loops)
1212  * pt is "IN"
1213  * else
1214  * Error! panic!
1215  *
1216  *
1217  * Values for "call_on_hits"
1218  * 1 find all elements pt touches, call user routine for each geom.
1219  * 2 find all elements pt touches, call user routine for each use
1220  *
1221  * in_or_out_only:
1222  * non-zero pt is known NOT to be on an EU of FU
1223  * 0 pt may be on an EU of FU
1224  *
1225  * Returns -
1226  * NMG_CLASS_AonB, etc...
1227  */
1228 int
1229 nmg_class_pt_fu_except(const fastf_t *pt, const struct faceuse *fu, const struct loopuse *ignore_lu,
1230  void (*eu_func) (struct edgeuse *, point_t, const char *), void (*vu_func) (struct vertexuse *, point_t, const char *), const char *priv,
1231  const int call_on_hits, const int in_or_out_only, const struct bn_tol *tol)
1232 
1233 /* func to call when pt on edgeuse */
1234 /* func to call when pt on vertexuse*/
1235 /* private data for [ev]u_func */
1236 
1237 
1238 {
1239  struct fpi fpi;
1240  struct loopuse *lu;
1241  int ot_same_in = 0;
1242  int ot_opposite_out = 0;
1243  int ot_same[4];
1244  int ot_opposite[4];
1245  int lu_class;
1246  int fu_class = NMG_CLASS_Unknown;
1247  double dist;
1248  struct ve_dist *ved_p;
1249  int i;
1250 
1251  if (RTG.NMG_debug & DEBUG_PT_FU)
1252  bu_log("nmg_class_pt_fu_except(pt=(%g %g %g), fu=%p)\n", V3ARGS(pt), (void *)fu);
1253 
1254  if (fu->orientation != OT_SAME) bu_bomb("nmg_class_pt_fu_except() not OT_SAME\n");
1255 
1256  NMG_CK_FACEUSE(fu);
1257  NMG_CK_FACE(fu->f_p);
1258  NMG_CK_FACE_G_PLANE(fu->f_p->g.plane_p);
1259  if (ignore_lu) NMG_CK_LOOPUSE(ignore_lu);
1260  BN_CK_TOL(tol);
1261 
1262  /* Validate distance from point to plane */
1263  NMG_GET_FU_PLANE(fpi.norm, fu);
1264  if ((dist=fabs(DIST_PT_PLANE(pt, fpi.norm))) > tol->dist) {
1265  bu_log("nmg_class_pt_fu_except() ERROR, point (%g, %g, %g)\nnot on face %g %g %g %g, \ndist=%g\n",
1266  V3ARGS(pt), V4ARGS(fpi.norm), dist);
1267  }
1268 
1269  if (V3PT_OUT_RPP_TOL(pt, fu->f_p->min_pt, fu->f_p->max_pt, tol->dist)) {
1270  /* point is not in RPP of face, so there's NO WAY this point
1271  * is anything but OUTSIDE
1272  */
1273  if (RTG.NMG_debug & DEBUG_PT_FU)
1274  bu_log("nmg_class_pt_fu_except((%g %g %g) outside face RPP\n",
1275  V3ARGS(pt));
1276 
1277  return NMG_CLASS_AoutB;
1278  }
1279 
1280  for (i = 0; i < 4; i++) {
1281  ot_same[i] = ot_opposite[i] = 0;
1282  }
1283 
1284  fpi.fu_p = fu;
1285  fpi.tol = tol;
1286  BU_LIST_INIT(&fpi.ve_dh);
1287  VMOVE(fpi.pt, pt);
1288  fpi.eu_func = eu_func;
1289  fpi.vu_func = vu_func;
1290  fpi.priv = priv;
1291  fpi.hits = call_on_hits;
1292  fpi.magic = NMG_FPI_MAGIC;
1293 
1294  for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
1295  if (ignore_lu && (ignore_lu==lu || ignore_lu==lu->lumate_p))
1296  continue;
1297 
1298  /* Ignore OT_BOOLPLACE, etc. */
1299  if (lu->orientation != OT_SAME && lu->orientation != OT_OPPOSITE) {
1300  if (lu->orientation != OT_BOOLPLACE) {
1301  bu_bomb("nmg_class_pt_fu_except() lu orientation is not OT_SAME, OT_OPPOSITE or OT_BOOLPLACE\n");
1302  }
1303  continue;
1304  }
1305 
1306  lu_class = nmg_class_pt_lu(lu, &fpi, in_or_out_only);
1307  if (RTG.NMG_debug & DEBUG_PT_FU)
1308  bu_log("loop %s says pt is %s\n",
1309  nmg_orientation(lu->orientation),
1310  nmg_class_name(lu_class));
1311 
1312  if (lu_class < 0 || lu_class > 3) {
1313  bu_log("nmg_class_pt_fu_except() lu_class=%s %d\n",
1314  nmg_class_name(lu_class), lu_class);
1315  bu_bomb("nmg_class_pt_fu_except() bad lu_class\n");
1316  }
1317 
1318  if (lu->orientation == OT_OPPOSITE) {
1319  ot_opposite[lu_class]++;
1320  if (lu_class == NMG_CLASS_AoutB) {
1321  ot_opposite_out++;
1322  }
1323  } else if (lu->orientation == OT_SAME) {
1324  ot_same[lu_class]++;
1325  if (lu_class == NMG_CLASS_AinB
1326  || lu_class == NMG_CLASS_AonBshared
1327  || lu_class == NMG_CLASS_AonBanti)
1328  {
1329  ot_same_in++;
1330  }
1331  }
1332  }
1333 
1334  if (RTG.NMG_debug & DEBUG_PT_FU) {
1335  bu_log("loops ot_same_in:%d ot_opposite_out:%d\n",
1336  ot_same_in, ot_opposite_out);
1337  bu_log("loops in/onS/onA/out ot_same=%d/%d/%d/%d ot_opp=%d/%d/%d/%d\n",
1338  ot_same[0], ot_same[1], ot_same[2], ot_same[3],
1339  ot_opposite[0], ot_opposite[1], ot_opposite[2], ot_opposite[3]);
1340  }
1341 
1342  if (ot_same_in == ot_opposite_out) {
1343  /* All the holes cancel out the solid loops */
1344  fu_class = NMG_CLASS_AoutB;
1345  } else if (ot_same_in > ot_opposite_out) {
1346  /* XXX How can this difference be > 1 ? */
1347  fu_class = NMG_CLASS_AinB;
1348  } else {
1349  /* Panic time! How did I get a parity mis-match? */
1350  bu_log("loops in/onS/onA/out ot_same=%d/%d/%d/%d ot_opp=%d/%d/%d/%d\n",
1351  ot_same[0], ot_same[1], ot_same[2], ot_same[3],
1352  ot_opposite[0], ot_opposite[1], ot_opposite[2], ot_opposite[3]);
1353  bu_log("nmg_class_pt_fu_except(%g %g %g)\nParity error @ %s:%d ot_same_in:%d ot_opposite_out:%d\n",
1354  V3ARGS(pt), __FILE__, __LINE__,
1355  ot_same_in, ot_opposite_out);
1356  bu_log("fu=%p\n", (void *)fu);
1357  nmg_pr_fu_briefly(fu, "");
1358 
1359  plot_parity_error(fu, pt);
1360 
1361  bu_bomb("nmg_class_pt_fu_except() loop classification parity error\n");
1362  }
1363 
1364  while (BU_LIST_WHILE(ved_p, ve_dist, &fpi.ve_dh)) {
1365  BU_LIST_DEQUEUE(&ved_p->l);
1366  bu_free((char *)ved_p, "ve_dist struct");
1367  }
1368 
1369 
1370  if (RTG.NMG_debug & DEBUG_PT_FU)
1371  bu_log("nmg_class_pt_fu_except() returns %s\n",
1372  nmg_class_name(fu_class));
1373 
1374  return fu_class;
1375 }
1376 
1377 
1378 /**
1379  * Classify a point as being in/on/out of the area bounded by a loop,
1380  * ignoring any uses of a particular edge in the loop.
1381  *
1382  * This routine must be called with a face-loop of edges!
1383  * It will not work properly on crack loops.
1384  */
1385 int
1386 nmg_class_pt_lu_except(fastf_t *pt, const struct loopuse *lu, const struct edge *e_p, const struct bn_tol *tol)
1387 {
1388  register struct edgeuse *eu;
1389  struct edge_info edge_list;
1390  struct edge_info *ei;
1391  struct fpi fpi;
1392  int lu_class = NMG_CLASS_Unknown;
1393  struct ve_dist *ved_p;
1394  double dist;
1395 
1396  if (RTG.NMG_debug & DEBUG_PT_FU) {
1397  bu_log("nmg_class_pt_lu_except((%g %g %g) %p ", V3ARGS(pt), (void *)e_p);
1398  if (e_p)
1399  bu_log(" e_p=(%g %g %g) <-> (%g %g %g))\n",
1400  V3ARGS(e_p->eu_p->vu_p->v_p->vg_p->coord),
1401  V3ARGS(e_p->eu_p->eumate_p->vu_p->v_p->vg_p->coord));
1402  else
1403  bu_log(" e_p=(NULL))\n");
1404  }
1405 
1406  NMG_CK_LOOPUSE(lu);
1407 
1408  if (e_p) NMG_CK_EDGE(e_p);
1409 
1410  NMG_CK_FACEUSE(lu->up.fu_p);
1411 
1412  /* Validate distance from point to plane */
1413  NMG_GET_FU_PLANE(fpi.norm, lu->up.fu_p);
1414  if ((dist=fabs(DIST_PT_PLANE(pt, fpi.norm))) > tol->dist) {
1415  bu_log("nmg_class_pt_lu_except() ERROR, point (%g, %g, %g)\nnot on face %g %g %g %g, \ndist=%g\n",
1416  V3ARGS(pt), V4ARGS(fpi.norm), dist);
1417  }
1418 
1419  if (V3PT_OUT_RPP_TOL(pt, lu->l_p->lg_p->min_pt, lu->l_p->lg_p->max_pt, tol->dist)) {
1420  /* point is not in RPP of loop */
1421 
1422  if (RTG.NMG_debug & DEBUG_PT_FU)
1423  bu_log("nmg_class_pt_lu_except(pt(%g %g %g) outside loop RPP\n",
1424  V3ARGS(pt));
1425 
1426  if (lu->orientation == OT_SAME) return NMG_CLASS_AoutB;
1427  else if (lu->orientation == OT_OPPOSITE) return NMG_CLASS_AinB;
1428  else {
1429  bu_log("What kind of loop is this anyway? %s?\n",
1430  nmg_orientation(lu->orientation));
1431  bu_bomb("");
1432  }
1433  }
1434 
1435  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
1436  bu_log("%s:%d Improper use of nmg_class_pt_lu_except(pt(%g %g %g), vu)\n",
1437  __FILE__, __LINE__, V3ARGS(pt));
1438  bu_bomb("giving up\n");
1439  }
1440 
1441  BU_LIST_INIT(&edge_list.l);
1442  fpi.fu_p = lu->up.fu_p;
1443 
1444  fpi.tol = tol;
1445  BU_LIST_INIT(&fpi.ve_dh);
1446  VMOVE(fpi.pt, pt);
1447  fpi.eu_func = (void (*)(struct edgeuse *, point_t, const char *))NULL;
1448  fpi.vu_func = (void (*)(struct vertexuse *, point_t, const char *))NULL;
1449  fpi.priv = (char *)NULL;
1450  fpi.hits = 0;
1451  fpi.magic = NMG_FPI_MAGIC;
1452 
1453  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
1454  if (eu->e_p == e_p) {
1455  if (RTG.NMG_debug & DEBUG_PT_FU)
1456  bu_log("skipping edgeuse (%g %g %g) -> (%g %g %g) on \"except\" edge\n",
1457  V3ARGS(eu->vu_p->v_p->vg_p->coord),
1458  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
1459 
1460  continue;
1461  }
1462 
1463  ei = nmg_class_pt_eu(&fpi, eu, &edge_list, 0);
1464  NMG_CK_EI(ei);
1465  NMG_CK_VED(ei->ved_p);
1466  if (ei->ved_p->dist < tol->dist_sq) {
1467  lu_class = NMG_CLASS_AonBshared;
1468  break;
1469  }
1470  }
1471  if (lu_class == NMG_CLASS_Unknown)
1472  lu_class = compute_loop_class(&fpi, lu, &edge_list);
1473  else if (RTG.NMG_debug & DEBUG_PT_FU)
1474  bu_log("loop class already known (pt must touch edge)\n");
1475 
1476  /* free up the edge_list elements */
1477  while (BU_LIST_WHILE(ei, edge_info, &edge_list.l)) {
1478  BU_LIST_DEQUEUE(&ei->l);
1479  bu_free((char *)ei, "edge info struct");
1480  }
1481 
1482  while (BU_LIST_WHILE(ved_p, ve_dist, &fpi.ve_dh)) {
1483  BU_LIST_DEQUEUE(&ved_p->l);
1484  bu_free((char *)ved_p, "ve_dist struct");
1485  }
1486 
1487  if (RTG.NMG_debug & DEBUG_PT_FU)
1488  bu_log("nmg_class_pt_lu_except() returns %s\n",
1489  nmg_class_name(lu_class));
1490 
1491  return lu_class;
1492 }
1493 
1494 
1495 /*
1496  * Local Variables:
1497  * mode: C
1498  * tab-width: 8
1499  * indent-tabs-mode: t
1500  * c-file-style: "stroustrup"
1501  * End:
1502  * ex: shiftwidth=4 tabstop=8
1503  */
#define BU_LIST_PNEXT_CIRC(structure, p)
Definition: list.h:442
#define NMG_CK_FPI(_fpi)
Definition: nmg_pt_fu.c:83
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
#define NMG_EDGEUSE_MAGIC
Definition: magic.h:120
void(* vu_func)(struct vertexuse *, point_t, const char *)
Definition: nmg_pt_fu.c:78
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
Definition: list.h:118
int hits
Definition: nmg_pt_fu.c:80
struct edgeuse * eu_p
Definition: nmg_pt_fu.c:64
int nmg_class_pt_fu_except(const point_t pt, const struct faceuse *fu, const struct loopuse *ignore_lu, void(*eu_func)(struct edgeuse *, point_t, const char *), void(*vu_func)(struct vertexuse *, point_t, const char *), const char *priv, const int call_on_hits, const int in_or_out_only, const struct bn_tol *tol)
double dist
>= 0
Definition: tol.h:73
int nmg_find_eu_left_non_unit(fastf_t *left, const struct edgeuse *eu)
Definition: nmg_info.c:1408
Definition: clone.c:90
#define NMG_FPI_TOUCHED
Definition: nmg_pt_fu.c:88
lu
Definition: nmg_mod.c:3855
struct vertex * v1
Definition: nmg_pt_fu.c:50
#define BU_LIST_IS_EMPTY(hp)
Definition: list.h:295
struct bu_list l
Definition: nmg_pt_fu.c:47
double dist_sq
dist * dist
Definition: tol.h:74
int nmg_class_pt_lu_except(fastf_t *pt, const struct loopuse *lu, const struct edge *e_p, const struct bn_tol *tol)
Definition: nmg_pt_fu.c:1386
#define BU_LIST_MAGIC_SET(_l, _magic)
Definition: list.h:129
int nmg_class_pt_euvu(const fastf_t *pt, struct edgeuse *eu_in, const struct bn_tol *tol)
Definition: nmg_pt_fu.c:319
HIDDEN void make_near_list(struct edge_info *edge_list, struct bu_list *near1, const struct bn_tol *tol)
Definition: nmg_pt_fu.c:732
#define BN_TOL_MAGIC
Definition: magic.h:74
#define SMALL_FASTF
Definition: defines.h:342
struct vertex * v2
Definition: nmg_pt_fu.c:51
Header file for the BRL-CAD common definitions.
#define BU_LIST_APPEND(old, new)
Definition: list.h:197
#define BU_LIST_NON_EMPTY(hp)
Definition: list.h:296
#define MAX_FASTF
Definition: defines.h:340
#define NMG_LOOPUSE_MAGIC
Definition: magic.h:130
void pl_erase(register FILE *plotfp)
Definition: plot3.c:271
#define HIDDEN
Definition: common.h:86
NMG_CK_LOOPUSE(lu)
struct bu_list ve_dh
Definition: nmg_pt_fu.c:74
#define NMG_CK_EI(_p)
Definition: nmg_pt_fu.c:68
void nmg_pl_fu(FILE *fp, const struct faceuse *fu, long *b, int red, int green, int blue)
Definition: nmg_plot.c:722
#define BU_LIST_PLAST(structure, p)
Definition: list.h:424
int nmg_loop_is_a_crack(const struct loopuse *lu)
Definition: nmg_info.c:449
void pd_3space(register FILE *plotfp, double px1, double py1, double pz1, double px2, double py2, double pz2)
Definition: plot3.c:580
uint32_t NMG_debug
debug bits for NMG's see nmg.h
Definition: raytrace.h:1699
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
#define NMG_VE_DIST_MAGIC
Definition: nmg_pt_fu.c:54
#define V3ARGS(a)
Definition: color.c:56
const struct bn_tol * tol
Definition: nmg_pt_fu.c:72
#define BU_LIST_PNEXT(structure, p)
Definition: list.h:422
uint32_t * magic_p
Definition: nmg_pt_fu.c:48
void(* eu_func)(struct edgeuse *, point_t, const char *)
Definition: nmg_pt_fu.c:77
struct bu_list l
Definition: nmg_pt_fu.c:62
goto out
Definition: nmg_mod.c:3846
Support for uniform tolerances.
Definition: tol.h:71
#define BN_CK_TOL(_p)
Definition: tol.h:82
#define BU_LIST_FIRST_MAGIC(hp)
Definition: list.h:416
oldeumate e_p
Definition: nmg_mod.c:3936
int nmg_find_eu_leftvec(fastf_t *left, const struct edgeuse *eu)
Definition: nmg_info.c:1274
uint32_t magic
Definition: nmg_pt_fu.c:71
#define NMG_VERTEXUSE_MAGIC
Definition: magic.h:145
#define BU_LIST_WHILE(p, structure, hp)
Definition: list.h:410
fastf_t dist
Definition: nmg_pt_fu.c:49
void pl_color(register FILE *plotfp, int r, int g, int b)
Definition: plot3.c:325
void nmg_pl_lu(FILE *fp, const struct loopuse *lu, long *b, int red, int green, int blue)
Definition: nmg_plot.c:710
void pdv_3line(register FILE *plotfp, const fastf_t *a, const fastf_t *b)
Definition: plot3.c:642
#define BU_LIST_PPREV_CIRC(structure, p)
Definition: list.h:450
#define ZERO(val)
Definition: units.c:38
#define BU_LIST_INIT(_hp)
Definition: list.h:148
plane_t norm
Definition: nmg_pt_fu.c:75
#define NMG_FPI_MAGIC
Definition: nmg_pt_fu.c:82
int bn_dist_pt3_lseg3(fastf_t *dist, point_t pca, const point_t a, const point_t b, const point_t p, const struct bn_tol *tol)
Find the distance from a point P to a line segment described by the two endpoints A and B...
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
#define BU_CK_LIST_HEAD(_p)
Definition: list.h:142
const char * nmg_class_name(int nmg_class)
Definition: nmg_eval.c:122
void nmg_pr_fu_briefly(const struct faceuse *fu, char *h)
Definition: nmg_pr.c:359
const struct faceuse * fu_p
Definition: nmg_pt_fu.c:73
int nmg_eu_is_part_of_crack(const struct edgeuse *eu)
Definition: nmg_pt_fu.c:276
int bn_distsq_pt3_lseg3(fastf_t *dist, const fastf_t *a, const fastf_t *b, const fastf_t *p, const struct bn_tol *tol)
Definition: nmg_pt_fu.c:127
int status
Definition: nmg_pt_fu.c:52
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
#define NMG_CK_VED(_p)
Definition: nmg_pt_fu.c:55
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
int nmg_class
Definition: nmg_pt_fu.c:65
HIDDEN const point_t delta
Definition: sh_prj.c:618
double fastf_t
Definition: defines.h:300
char * nmg_orientation(int orientation)
Definition: nmg_pr.c:50
#define BU_LIST_NOT_HEAD(p, hp)
Definition: list.h:324
#define NMG_FPI_MISSED
Definition: nmg_pt_fu.c:89
#define NMG_EDGE_INFO_MAGIC
Definition: nmg_pt_fu.c:67
struct ve_dist * ved_p
Definition: nmg_pt_fu.c:63
const char * priv
Definition: nmg_pt_fu.c:79
#define BU_LIST_FIRST(structure, hp)
Definition: list.h:312
Definition: nmg_pt_fu.c:70
point_t pt
Definition: nmg_pt_fu.c:76
#define UNLIKELY(expression)
Definition: common.h:282
struct rt_g RTG
Definition: globals.c:39