BRL-CAD
nmg_inter.c
Go to the documentation of this file.
1 /* N M G _ I N T E R . 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_inter.c
23  *
24  * Routines to intersect two NMG regions. When complete, all loops in
25  * each region have a single classification w.r.t. the other region,
26  * i.e. all geometric intersections of the two regions have explicit
27  * topological representations.
28  *
29  * The intersector makes sure that all geometric intersections gets
30  * recorded with explicit geometry and topology that is shared between
31  * both regions. Primary examples of this are (a) the line of
32  * intersection between two planes (faces), and (b) the point of
33  * intersection where two edges cross.
34  *
35  * Entities of one region that are INSIDE, but not ON the other region
36  * do not become shared during the intersection process.
37  *
38  * All point -vs- point comparisons should be done in 3D, for
39  * consistency.
40  *
41  * Method -
42  *
43  * Find all the points of intersection between the two regions, and
44  * insert vertices at those points, breaking edges on those new
45  * vertices as appropriate.
46  *
47  * Call the face cutter to construct and delete edges and loops along
48  * the line of intersection, as appropriate.
49  *
50  * There are no "user interface" routines in here.
51  *
52  */
53 /** @} */
54 
55 #include "common.h"
56 
57 #include <stddef.h>
58 #include <math.h>
59 #include <string.h>
60 #include "bio.h"
61 
62 #include "vmath.h"
63 #include "nmg.h"
64 #include "raytrace.h"
65 #include "plot3.h"
66 
67 
68 #define ISECT_NONE 0
69 #define ISECT_SHARED_V 1
70 #define ISECT_SPLIT1 2
71 #define ISECT_SPLIT2 4
72 
73 struct ee_2d_state {
75  struct edgeuse *eu;
76  point_t start;
77  point_t end;
78  vect_t dir;
79 };
80 
81 
82 static int nmg_isect_edge2p_face2p(struct nmg_inter_struct *is,
83  struct edgeuse *eu, struct faceuse *fu,
84  struct faceuse *eu_fu);
85 
86 
87 static struct nmg_inter_struct *nmg_hack_last_is; /* see nmg_isect2d_final_cleanup() */
88 
89 struct vertexuse *
90 nmg_make_dualvu(struct vertex *v, struct faceuse *fu, const struct bn_tol *tol)
91 {
92  struct loopuse *lu;
93  struct vertexuse *dualvu;
94  struct edgeuse *new_eu;
95 
96  NMG_CK_VERTEX(v);
97  NMG_CK_FACEUSE(fu);
98  BN_CK_TOL(tol);
99 
100  if (RTG.NMG_debug & DEBUG_POLYSECT)
101  bu_log("nmg_make_dualvu(v=%p, fu=%p)\n", (void *)v, (void *)fu);
102 
103  /* check for existing vu */
104  dualvu=nmg_find_v_in_face(v, fu);
105  if (dualvu) {
106  if (RTG.NMG_debug & DEBUG_POLYSECT)
107  bu_log("\tdualvu already exists (%p)\n", (void *)dualvu);
108  return dualvu;
109  }
110 
111  new_eu = (struct edgeuse *)NULL;
112 
113  /* check if v lies within tolerance of an edge in face */
114  if (RTG.NMG_debug & DEBUG_POLYSECT)
115  bu_log("\tLooking for an edge to split\n");
116  for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
117  struct edgeuse *eu;
118 
119  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
120  continue;
121 
122  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
123  int code;
124  fastf_t dist;
125  point_t pca;
126 
127  if (RTG.NMG_debug & DEBUG_POLYSECT)
128  bu_log("\tChecking eu %p (%f %f %f) <-> (%f %f %f)\n",
129  (void *)eu,
130  V3ARGS(eu->vu_p->v_p->vg_p->coord),
131  V3ARGS(eu->eumate_p->vu_p->v_p->vg_p->coord));
132 
133  code = bn_dist_pt3_lseg3(&dist, pca,
134  eu->vu_p->v_p->vg_p->coord,
135  eu->eumate_p->vu_p->v_p->vg_p->coord,
136  v->vg_p->coord, tol);
137 
138  if (RTG.NMG_debug & DEBUG_POLYSECT)
139  bu_log("bn_dist_pt3_lseg3 returns %d, dist=%f\n", code, dist);
140 
141  if (code > 2)
142  continue;
143 
144  /* v is within tolerance of eu */
145  if (code > 0)
146  continue;
147 
148  /* split edge */
149  if (RTG.NMG_debug & DEBUG_POLYSECT)
150  bu_log("nmg_make_dualvu is splitting eu %p at v %p\n", (void *)eu, (void *)v);
151  new_eu = nmg_esplit(v, eu, 1);
152  }
153  }
154 
155  if (new_eu)
156  return new_eu->vu_p;
157 
158  /* need a self loop */
159  lu = nmg_mlv(&fu->l.magic, v, OT_BOOLPLACE);
160  if (RTG.NMG_debug & DEBUG_POLYSECT)
161  bu_log("nmg_make_dualvu is making a self_loop (lu=%p, vu=%p) for v=%p\n",
162  (void *)lu, (void *)BU_LIST_FIRST(vertexuse, &lu->down_hd),
163  (void *)v);
164  nmg_loop_g(lu->l_p, tol);
165  return BU_LIST_FIRST(vertexuse, &lu->down_hd);
166 }
167 
168 
169 /**
170  * Given a vu which represents a point of intersection between shells
171  * s1 and s2, insert it and its dual into lists l1 and l2.
172  * First, determine whether the vu came from s1 or s2, and insert in
173  * the corresponding list.
174  *
175  * Second, try and find a dual of that vertex
176  * in the other shell's faceuse (fu1 or fu2)
177  * (if the entity in the other shell is not a wire), and enlist the dual.
178  * If there is no dual, make a self-loop over there, and enlist that.
179  *
180  * If 'dualvu' is provided, don't search, just use that.
181  *
182  * While it is true that in most cases the calling routine will know
183  * which shell the vu came from, it's cheap to re-determine it here.
184  * This "all in one" packaging, which handles both lists automatically
185  * is *vastly* superior to the previous version, which pushed 10-20
186  * lines of bookkeeping up into *every* place an intersection vu was
187  * created.
188  *
189  * Returns a pointer to vu's dual.
190  *
191  * "Join the Army, young vertexuse".
192  */
193 struct vertexuse *
194 nmg_enlist_vu(struct nmg_inter_struct *is, const struct vertexuse *vu, struct vertexuse *dualvu, fastf_t dist)
195 
196 
197 /* vu's dual in other shell. May be NULL */
198 /* distance along intersect ray for this vu */
199 {
200  struct shell *sv; /* shell of vu */
201  struct loopuse *lu; /* lu of new self-loop */
202  struct faceuse *dualfu = (struct faceuse *)NULL; /* faceuse of vu's dual */
203  struct shell *duals = (struct shell *)NULL; /* shell of vu's dual */
204  struct faceuse *fuv; /* faceuse of vu */
205 
207  NMG_CK_VERTEXUSE(vu);
208  if (dualvu) {
209  NMG_CK_VERTEXUSE(dualvu);
210  if (vu == dualvu) bu_bomb("nmg_enlist_vu() vu == dualvu\n");
211  }
212 
213  if (is->mag_len <= BU_PTBL_END(is->l1) || is->mag_len <= BU_PTBL_END(is->l2))
214  bu_log("Array for distances to vertexuses is too small (%d)\n", is->mag_len);
215 
216  sv = nmg_find_s_of_vu(vu);
217  fuv = nmg_find_fu_of_vu(vu);
218 
219  /* First step: add vu to corresponding list */
220  if (sv == is->s1) {
221  bu_ptbl_ins_unique(is->l1, (long *)&vu->l.magic);
222  if (is->mag_len <= BU_PTBL_END(is->l1)) {
223  if (is->mag_len) {
224  is->mag_len *= 2;
225  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
226  "is->mag1");
227  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
228  "is->mag2");
229  } else {
230  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
231  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
232  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
233  }
234 
235  }
236  if (dist < MAX_FASTF)
237  is->mag1[bu_ptbl_locate(is->l1, (long *)&vu->l.magic)] = dist;
238  duals = is->s2; /* other shell */
239  dualfu = is->fu2;
240  if (is->fu1 && is->fu1->s_p != is->s1) bu_bomb("nmg_enlist_vu() fu1/s1 mismatch\n");
241  if (fuv != is->fu1) {
242  bu_log("fuv=%p, fu1=%p, fu2=%p\n", (void *)fuv, (void *)is->fu1, (void *)is->fu2);
243  bu_log("\tvu=%p (%p)\n", (void *)vu, (void *)vu->v_p);
244  bu_bomb("nmg_enlist_vu() vu/fu1 mis-match\n");
245  }
246  } else if (sv == is->s2) {
247  bu_ptbl_ins_unique(is->l2, (long *)&vu->l.magic);
248  if (is->mag_len <= BU_PTBL_END(is->l2)) {
249  if (is->mag_len) {
250  is->mag_len *= 2;
251  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
252  "is->mag1");
253  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
254  "is->mag2");
255  } else {
256  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
257  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
258  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
259  }
260 
261  }
262  if (dist < MAX_FASTF)
263  is->mag2[bu_ptbl_locate(is->l2, (long *)&vu->l.magic)] = dist;
264  duals = is->s1; /* other shell */
265  dualfu = is->fu1;
266  if (is->fu2 && is->fu2->s_p != is->s2) bu_bomb("nmg_enlist_vu() fu2/s2 mismatch\n");
267  if (fuv != is->fu2) {
268  bu_log("fuv=%p, fu1=%p, fu2=%p\n", (void *)fuv, (void *)is->fu1, (void *)is->fu2);
269  bu_log("\tvu=%p (%p)\n", (void *)vu, (void *)vu->v_p);
270  bu_bomb("nmg_enlist_vu() vu/fu2 mis-match\n");
271  }
272  } else {
273  bu_log("nmg_enlist_vu(vu=%p, dv=%p) sv=%p, s1=%p, s2=%p\n",
274  (void *)vu, (void *)dualvu, (void *)sv, (void *)is->s1, (void *)is->s2);
275  bu_bomb("nmg_enlist_vu: vu is not in s1 or s2\n");
276  }
277 
278  if (dualvu) {
279  if (vu->v_p != dualvu->v_p) bu_bomb("nmg_enlist_vu() dual vu has different vertex\n");
280  if (nmg_find_s_of_vu(dualvu) != duals) {
281  bu_log("nmg_enlist_vu(vu=%p, dv=%p) sv=%p, s1=%p, s2=%p, sdual=%p\n",
282  (void *)vu, (void *)dualvu,
283  (void *)sv, (void *)is->s1, (void *)is->s2, (void *)nmg_find_s_of_vu(dualvu));
284  bu_bomb("nmg_enlist_vu() dual vu shell mis-match\n");
285  }
286  if (dualfu && nmg_find_fu_of_vu(dualvu) != dualfu) bu_bomb("nmg_enlist_vu() dual vu has wrong fu\n");
287  }
288 
289  /* Second, search for vu's dual */
290  if (dualfu) {
291  NMG_CK_FACEUSE(dualfu);
292  if (dualfu->s_p != duals) bu_bomb("nmg_enlist_vu() dual fu's shell is not dual's shell?\n");
293  if (!dualvu)
294  dualvu = nmg_make_dualvu(vu->v_p, dualfu, &(is->tol));
295  else {
296  if (RTG.NMG_debug & DEBUG_POLYSECT) {
297  bu_log("nmg_enlist_vu(vu=%p, dv=%p) re-using dualvu=%p from dualfu=%p\n",
298  (void *)vu, (void *)dualvu,
299  (void *)dualvu, (void *)dualfu);
300  }
301  }
302  } else {
303  /* Must have come from a wire in other shell, make wire loop */
304  bu_log("\tvu=%p, %s, fu1=%p, fu2=%p\n",
305  (void *)vu, (sv==is->s1) ? "shell 1" : "shell 2",
306  (void *)is->fu1, (void *)is->fu2);
307  bu_log("nmg_enlist_vu(): QUESTION: What do I search for wire intersections? Making self-loop\n");
308  if (!dualvu) {
309  dualvu = nmg_find_v_in_shell(vu->v_p, duals, 0);
310  if (!dualvu) {
311  /* Not found, make self-loop in dual shell */
312  lu = nmg_mlv(&duals->l.magic, vu->v_p, OT_BOOLPLACE);
313  nmg_loop_g(lu->l_p, &(is->tol));
314  dualvu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
315  }
316  }
317  }
318  NMG_CK_VERTEXUSE(dualvu);
319 
320  /* Enlist the dual onto the other list */
321  if (sv == is->s1) {
322  bu_ptbl_ins_unique(is->l2, (long *)&dualvu->l.magic);
323  if (is->mag_len <= BU_PTBL_END(is->l2)) {
324  if (is->mag_len) {
325  is->mag_len *= 2;
326  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
327  "is->mag1");
328  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
329  "is->mag2");
330  } else {
331  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
332  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
333  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
334  }
335 
336  }
337  if (dist < MAX_FASTF)
338  is->mag2[bu_ptbl_locate(is->l2, (long *)&dualvu->l.magic)] = dist;
339  } else {
340  bu_ptbl_ins_unique(is->l1, (long *)&dualvu->l.magic);
341  if (is->mag_len <= BU_PTBL_END(is->l1)) {
342  if (is->mag_len) {
343  is->mag_len *= 2;
344  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
345  "is->mag1");
346  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
347  "is->mag2");
348  } else {
349  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
350  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
351  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
352  }
353 
354  }
355  if (dist < MAX_FASTF)
356  is->mag1[bu_ptbl_locate(is->l1, (long *)&dualvu->l.magic)] = dist;
357  }
358 
359  if (RTG.NMG_debug & DEBUG_POLYSECT) {
360  bu_log("nmg_enlist_vu(vu=%p, dv=%p) v=%p, dist=%g (%s) ret=%p\n",
361  (void *)vu, (void *)dualvu, (void *)vu->v_p, dist,
362  (sv == is->s1) ? "shell 1" : "shell 2",
363  (void *)dualvu);
364  }
365 
366  /* Some (expensive) centralized sanity checking */
367  if ((RTG.NMG_debug & DEBUG_VERIFY) && is->fu1 && is->fu2) {
368  nmg_ck_v_in_2fus(vu->v_p, is->fu1, is->fu2, &(is->tol));
369  }
370  return dualvu;
371 }
372 
373 
374 /**
375  * A "lazy evaluator" to obtain the 2D projection of a vertex.
376  * The lazy approach is not a luxury, since new (3D) vertices are created
377  * as the edge/edge intersection proceeds, and their 2D coordinates may
378  * be needed later on in the calculation.
379  * The alternative would be to store the 2D projection each time a
380  * new vertex is created, but that is likely to be a lot of bothersome
381  * code, where one omission would be deadly.
382  *
383  * The return is a 3-tuple, with the Z coordinate set to 0.0 for safety.
384  * This is especially useful when the projected value is printed using
385  * one of the 3D print routines.
386  *
387  * 'assoc_use' is either a pointer to a faceuse, or an edgeuse.
388  */
389 static void
390 nmg_get_2d_vertex(fastf_t *v2d, struct vertex *v, struct nmg_inter_struct *is, const uint32_t *assoc_use)
391 /* a 3-tuple */
392 
393 
394 /* ptr to faceuse/edgeuse associated w/2d projection */
395 {
396  register fastf_t *pt2d;
397  point_t pt;
398  struct vertex_g *vg;
399  uint32_t *this_obj;
400 
402  NMG_CK_VERTEX(v);
403 
404  /* If 2D preparations have not been made yet, do it now */
405  if (!is->vert2d) {
406  nmg_isect2d_prep(is, assoc_use);
407  }
408 
409  if (*assoc_use == NMG_FACEUSE_MAGIC) {
410  this_obj = &((struct faceuse *)assoc_use)->f_p->l.magic;
411  if (this_obj != is->twod)
412  goto bad;
413  } else if (*assoc_use == NMG_EDGEUSE_MAGIC) {
414  this_obj = &((struct edgeuse *)assoc_use)->e_p->magic;
415  if (this_obj != is->twod)
416  goto bad;
417  } else {
418  this_obj = NULL;
419  bad:
420  if (this_obj) {
421  bu_log("nmg_get_2d_vertex(, assoc_use=%p %s) this_obj=%p %s, is->twod=%p %s\n",
422  (void *)assoc_use, bu_identify_magic(*assoc_use),
423  (void *)this_obj, bu_identify_magic(*this_obj),
424  (void *)is->twod, bu_identify_magic(*(is->twod)));
425  } else {
426  bu_log("nmg_get_2d_vertex - this is NULL\n");
427  }
428  bu_bomb("nmg_get_2d_vertex: 2d association mis-match\n");
429  }
430 
431  if (!v->vg_p) {
432  bu_log("nmg_get_2d_vertex: v=%p, assoc_use=%p, null vg_p\n",
433  (void *)v, (void *)assoc_use);
434  bu_bomb("nmg_get_2d_vertex: vertex with no geometry!\n");
435  }
436  vg = v->vg_p;
437  NMG_CK_VERTEX_G(vg);
438  if (v->index >= is->maxindex) {
439  struct model *m;
440  int oldmax;
441  register int i;
442 
443  oldmax = is->maxindex;
444  m = nmg_find_model(&v->magic);
445  NMG_CK_MODEL(m);
446  bu_log("nmg_get_2d_vertex: v=%p, v->index=%ld, is->maxindex=%d, m->maxindex=%ld\n",
447  (void *)v, v->index, is->maxindex, m->maxindex);
448  if (v->index >= m->maxindex) {
449  /* Really off the end */
450  VPRINT("3d vertex", vg->coord);
451  bu_bomb("nmg_get_2d_vertex: array overrun\n");
452  }
453  /* Need to extend array, it's grown. */
454  is->maxindex = m->maxindex * 4;
455  if (RTG.NMG_debug & DEBUG_POLYSECT) {
456  bu_log("nmg_get_2d_vertex() extending vert2d array from %d to %d points (m max=%ld)\n",
457  oldmax, is->maxindex, m->maxindex);
458  }
459  is->vert2d = (fastf_t *)bu_realloc((char *)is->vert2d,
460  is->maxindex * 3 * sizeof(fastf_t), "vert2d[]");
461 
462  /* Clear out the new part of the 2D vertex array, setting flag in [2] to -1 */
463  for (i = (3*is->maxindex)-1-2; i >= oldmax*3; i -= 3) {
464  VSET(&is->vert2d[i], 0, 0, -1);
465  }
466  }
467  pt2d = &is->vert2d[v->index*3];
468  if (ZERO(pt2d[2])) {
469  /* Flag set. Conversion is done. Been here before */
470  v2d[0] = pt2d[0];
471  v2d[1] = pt2d[1];
472  v2d[2] = 0;
473  return;
474  }
475 
476  MAT4X3PNT(pt, is->proj, vg->coord);
477  v2d[0] = pt2d[0] = pt[0];
478  v2d[1] = pt2d[1] = pt[1];
479  v2d[2] = pt2d[2] = 0; /* flag */
480 
481  if (!NEAR_ZERO(pt[2], is->tol.dist)) {
482  struct faceuse *fu = (struct faceuse *)assoc_use;
483  plane_t n;
484  fastf_t dist;
485  NMG_GET_FU_PLANE(n, fu);
486  dist = DIST_PT_PLANE(vg->coord, n);
487  bu_log("nmg_get_2d_vertex ERROR #%ld (%g %g %g) becomes (%g, %g)\n\t%g != zero, dist3d=%g, %g*tol\n",
488  v->index, V3ARGS(vg->coord), V3ARGS(pt),
489  dist, dist/is->tol.dist);
490  if (!NEAR_ZERO(dist, is->tol.dist) &&
491  !NEAR_ZERO(pt[2], 10*is->tol.dist)) {
492  bu_log("nmg_get_2d_vertex(, assoc_use=%p) f=%p, is->twod=%p\n",
493  (void *)assoc_use, (void *)fu->f_p, (void *)is->twod);
494  PLPRINT("fu->f_p N", n);
495  bu_bomb("3D->2D point projection error\n");
496  }
497  }
498 
499  if (RTG.NMG_debug & DEBUG_POLYSECT) {
500  bu_log("2d #%ld (%g %g %g) becomes (%g, %g) %g\n",
501  v->index, V3ARGS(vg->coord), V3ARGS(pt));
502  }
503 }
504 
505 
506 /**
507  * To intersect two co-planar faces, project all vertices from those
508  * faces into 2D.
509  * At the moment, use a memory intensive strategy which allocates a
510  * (3d) point_t for each "index" item, and subscripts the resulting
511  * array by the vertices index number.
512  * Since additional vertices can be created as the intersection process
513  * operates, 2*maxindex items are originally allocated, as a (generous)
514  * upper bound on the amount of intersecting that might happen.
515  *
516  * In the array, the third double of each projected vertex is set to -1 when
517  * that slot has not been filled yet, and 0 when it has been.
518  */
519 /* XXX Set this up so that it can take either an edge pointer
520  * or a face pointer. In case of edge, make edge_g_lseg->dir unit, and
521  * rotate that to +X axis. Make edge_g_lseg->pt be the origin.
522  * This will allow the 2D routines to operate on wires.
523  */
524 void
525 nmg_isect2d_prep(struct nmg_inter_struct *is, const uint32_t *assoc_use)
526 {
527  struct model *m;
528  struct face_g_plane *fg;
529  vect_t to;
530  point_t centroid;
531  point_t centroid_proj;
532  plane_t n;
533  register int i;
534 
536 
537  if (*assoc_use == NMG_FACEUSE_MAGIC) {
538  if (&((struct faceuse *)assoc_use)->f_p->l.magic == is->twod)
539  return; /* Already prepped */
540  } else if (*assoc_use == NMG_EDGEUSE_MAGIC) {
541  if (&((struct edgeuse *)assoc_use)->e_p->magic == is->twod)
542  return; /* Already prepped */
543  } else {
544  bu_bomb("nmg_isect2d_prep() bad assoc_use magic\n");
545  }
546 
548  nmg_hack_last_is = is;
549 
550  m = nmg_find_model(assoc_use);
551 
552  is->maxindex = (2 * m->maxindex);
553  is->vert2d = (fastf_t *)bu_malloc(is->maxindex * 3 * sizeof(fastf_t), "vert2d[]");
554 
555  if (*assoc_use == NMG_FACEUSE_MAGIC) {
556  struct faceuse *fu1 = (struct faceuse *)assoc_use;
557  struct face *f1;
558 
559  f1 = fu1->f_p;
560  fg = f1->g.plane_p;
561  NMG_CK_FACE_G_PLANE(fg);
562  is->twod = &f1->l.magic;
563  if (f1->flip) {
564  VREVERSE(n, fg->N);
565  n[W] = -fg->N[W];
566  } else {
567  HMOVE(n, fg->N);
568  }
569  if (RTG.NMG_debug & DEBUG_POLYSECT) {
570  bu_log("nmg_isect2d_prep(f=%p) flip=%d\n", (void *)f1, f1->flip);
571  PLPRINT("N", n);
572  }
573 
574  /*
575  * Rotate so that f1's N vector points up +Z.
576  * This places all 2D calculations in the XY plane.
577  * Translate so that f1's centroid becomes the 2D origin.
578  * Reasoning: no vertex should be favored by putting it at
579  * the origin. The "desirable" floating point space in the
580  * vicinity of the origin should be used to best advantage,
581  * by centering calculations around it.
582  */
583  VSET(to, 0, 0, 1);
584  bn_mat_fromto(is->proj, n, to, &is->tol);
585  VADD2SCALE(centroid, f1->max_pt, f1->min_pt, 0.5);
586  MAT4X3PNT(centroid_proj, is->proj, centroid);
587  centroid_proj[Z] = n[W]; /* pull dist from origin off newZ */
588  MAT_DELTAS_VEC_NEG(is->proj, centroid_proj);
589  } else if (*assoc_use == NMG_EDGEUSE_MAGIC) {
590  struct edgeuse *eu1 = (struct edgeuse *)assoc_use;
591  struct edge *e1;
592  struct edge_g_lseg *eg;
593 
594  bu_log("2d prep for edgeuse\n");
595  e1 = eu1->e_p;
596  NMG_CK_EDGE(e1);
597  eg = eu1->g.lseg_p;
598  NMG_CK_EDGE_G_LSEG(eg);
599  is->twod = &e1->magic;
600 
601  /*
602  * Rotate so that eg's eg_dir vector points up +X.
603  * The choice of the other axes is arbitrary.
604  * This ensures that all calculations happen on the XY plane.
605  * Translate the edge start point to the origin.
606  */
607  VSET(to, 1, 0, 0);
608  bn_mat_fromto(is->proj, eg->e_dir, to, &is->tol);
609  MAT_DELTAS_VEC_NEG(is->proj, eg->e_pt);
610  } else {
611  bu_bomb("nmg_isect2d_prep() bad assoc_use magic\n");
612  }
613 
614  /* Clear out the 2D vertex array, setting flag in [2] to -1 */
615  for (i = (3*is->maxindex)-1-2; i >= 0; i -= 3) {
616  VSET(&is->vert2d[i], 0, 0, -1);
617  }
618 }
619 
620 
621 /**
622  * Common routine to zap 2d vertex cache, and release dynamic storage.
623  */
624 void
626 {
628 
629  nmg_hack_last_is = (struct nmg_inter_struct *)NULL;
630 
631  if (!is->vert2d) return;
632  bu_free((char *)is->vert2d, "vert2d");
633  is->vert2d = NULL;
634  is->twod = NULL;
635 }
636 
637 
638 /**
639  * XXX Hack routine used for storage reclamation by G-JACK for
640  * XXX calculation of the reportcard without gobbling lots of memory
641  * XXX on bu_bomb() longjmp()s.
642  * Can be called by the longjmp handler with impunity.
643  * If a pointer to busy dynamic memory is still handy, it will be freed.
644  * If not, no harm done.
645  */
646 void
648 {
649  if (nmg_hack_last_is && nmg_hack_last_is->magic == NMG_INTER_STRUCT_MAGIC)
650  nmg_isect2d_cleanup(nmg_hack_last_is);
651 }
652 
653 
654 /**
655  * Handle the complete intersection of a vertex which lies on the
656  * plane of a face. *every* intersection is performed.
657  *
658  * If already part of the topology of the face, do nothing more.
659  * If it intersects one of the edges of the face, break the edge there.
660  * Otherwise, add a self-loop into the face as a marker.
661  *
662  * All vertexuse pairs are enlisted on the intersection line.
663  * Assuming that there is one (is->l1 non null).
664  *
665  * Called by -
666  * nmg_isect_3vertex_3face()
667  * nmg_isect_two_face2p()
668  */
669 void
670 nmg_isect_vert2p_face2p(struct nmg_inter_struct *is, struct vertexuse *vu1, struct faceuse *fu2)
671 {
672  struct vertexuse *vu2;
673  struct loopuse *lu2;
674  pointp_t pt;
675  int ret = 0;
676 
677  if (RTG.NMG_debug & DEBUG_POLYSECT)
678  bu_log("nmg_isect_vert2p_face2p(, vu1=%p, fu2=%p)\n", (void *)vu1, (void *)fu2);
680  NMG_CK_VERTEXUSE(vu1);
681  NMG_CK_FACEUSE(fu2);
682 
683  pt = vu1->v_p->vg_p->coord;
684 
685  /* Prep the 2D cache, if the face changed */
686  nmg_isect2d_prep(is, &fu2->l.magic);
687 
688  /* For every edge and vert, check topo AND geometric intersection */
689  for (BU_LIST_FOR(lu2, loopuse, &fu2->lu_hd)) {
690  struct edgeuse *eu2;
691 
692  NMG_CK_LOOPUSE(lu2);
693  if (BU_LIST_FIRST_MAGIC(&lu2->down_hd) == NMG_VERTEXUSE_MAGIC) {
694  vu2 = BU_LIST_FIRST(vertexuse, &lu2->down_hd);
695  if (vu1->v_p == vu2->v_p) {
696 
697  if (is->l1) nmg_enlist_vu(is, vu1, vu2, MAX_FASTF);
698  ret++;
699  continue;
700  }
701  /* Use 3D comparisons for uniformity */
702  if (bn_pt3_pt3_equal(pt, vu2->v_p->vg_p->coord, &is->tol)) {
703  /* Fuse the two verts together */
704  nmg_jv(vu1->v_p, vu2->v_p);
705  if (is->l1) nmg_enlist_vu(is, vu1, vu2, MAX_FASTF);
706  ret++;
707  continue;
708  }
709  continue;
710  }
711  for (BU_LIST_FOR(eu2, edgeuse, &lu2->down_hd)) {
712  struct edgeuse *new_eu;
713 
714  if (eu2->vu_p->v_p == vu1->v_p) {
715  if (is->l1) nmg_enlist_vu(is, vu1, eu2->vu_p, MAX_FASTF);
716  ret++;
717  continue;
718  }
719 
720  new_eu = nmg_break_eu_on_v(eu2, vu1->v_p, fu2, is);
721  if (new_eu) {
722  if (is->l1) nmg_enlist_vu(is, vu1, new_eu->vu_p, MAX_FASTF);
723  ret++;
724  continue;
725  }
726  }
727  }
728 
729  if (ret == 0) {
730  /* The vertex lies in the face, but touches nothing. Place marker */
731  if (RTG.NMG_debug & DEBUG_POLYSECT)
732  VPRINT("Making vertexloop", pt);
733 
734  lu2 = nmg_mlv(&fu2->l.magic, vu1->v_p, OT_BOOLPLACE);
735  nmg_loop_g(lu2->l_p, &is->tol);
736  vu2 = BU_LIST_FIRST(vertexuse, &lu2->down_hd);
737  if (is->l1) nmg_enlist_vu(is, vu1, vu2, MAX_FASTF);
738  }
739 }
740 
741 
742 /**
743  * intersect a vertex with a face (primarily for intersecting
744  * loops of a single vertex with a face).
745  *
746  * XXX It would be useful to have one of the new vu's in fu returned
747  * XXX as a flag, so that nmg_find_v_in_face() wouldn't have to be called
748  * XXX to re-determine what was just done.
749  */
750 static void
751 nmg_isect_3vertex_3face(struct nmg_inter_struct *is, struct vertexuse *vu, struct faceuse *fu)
752 {
753  struct vertexuse *vup;
754  pointp_t pt;
755  fastf_t dist;
756  plane_t n;
757 
759  NMG_CK_VERTEXUSE(vu);
760  NMG_CK_VERTEX(vu->v_p);
761  NMG_CK_FACEUSE(fu);
762 
763  if (RTG.NMG_debug & DEBUG_POLYSECT)
764  bu_log("nmg_isect_3vertex_3face(, vu=%p, fu=%p) v=%p\n", (void *)vu, (void *)fu, (void *)vu->v_p);
765 
766  /* check the topology first */
767  vup=nmg_find_v_in_face(vu->v_p, fu);
768  if (vup) {
769  if (RTG.NMG_debug & DEBUG_POLYSECT) bu_log("\tvu lies in face (topology 1)\n");
770  (void)bu_ptbl_ins_unique(is->l1, (long *)&vu->l.magic);
771  (void)bu_ptbl_ins_unique(is->l2, (long *)&vup->l.magic);
772  return;
773  }
774 
775 
776  /* since the topology didn't tell us anything, we need to check with
777  * the geometry
778  */
779  pt = vu->v_p->vg_p->coord;
780  NMG_GET_FU_PLANE(n, fu);
781  dist = DIST_PT_PLANE(pt, n);
782 
783  if (!NEAR_ZERO(dist, is->tol.dist)) {
784  if (RTG.NMG_debug & DEBUG_POLYSECT) bu_log("\tvu not on face (geometry)\n");
785  return;
786  }
787 
788  /*
789  * The point lies on the plane of the face, by geometry.
790  * This is now a 2-D problem.
791  */
792  (void)nmg_isect_vert2p_face2p(is, vu, fu);
793 }
794 
795 
796 /**
797  * Having decided that an edge(use) crosses a plane of intersection,
798  * stick a vertex at the point of intersection along the edge.
799  *
800  * vu1_final in fu1 is BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p after return.
801  * vu2_final is the returned value, and is in fu2.
802  *
803  */
804 static struct vertexuse *
805 nmg_break_3edge_at_plane(const fastf_t *hit_pt, struct faceuse *fu2, struct nmg_inter_struct *is, struct edgeuse *eu1)
806 
807 /* The face that eu intersects */
808 
809 /* Edge to be broken (in fu1) */
810 {
811  struct vertexuse *vu1_final;
812  struct vertexuse *vu2_final; /* hit_pt's vu in fu2 */
813  struct vertex *v2;
814  struct loopuse *plu2; /* "point" loopuse */
815  struct edgeuse *eu1forw; /* New eu, after break, forw of eu1 */
816  struct vertex *v1;
817  struct vertex *v1mate;
818  fastf_t dist;
819 
821  NMG_CK_EDGEUSE(eu1);
822 
823  v1 = eu1->vu_p->v_p;
824  NMG_CK_VERTEX(v1);
825  v1mate = eu1->eumate_p->vu_p->v_p;
826  NMG_CK_VERTEX(v1mate);
827 
828  /* Intersection is between first and second vertex points.
829  * Insert new vertex at intersection point.
830  */
831  if (RTG.NMG_debug & DEBUG_POLYSECT) {
832  bu_log("nmg_break_3edge_at_plane() Splitting %g, %g, %g <-> %g, %g, %g\n",
833  V3ARGS(v1->vg_p->coord),
834  V3ARGS(v1mate->vg_p->coord));
835  VPRINT("\tAt point of intersection", hit_pt);
836  }
837 
838  /* Double check for bad behavior */
839  if (bn_pt3_pt3_equal(hit_pt, v1->vg_p->coord, &is->tol))
840  bu_bomb("nmg_break_3edge_at_plane() hit_pt equal to v1\n");
841  if (bn_pt3_pt3_equal(hit_pt, v1mate->vg_p->coord, &is->tol))
842  bu_bomb("nmg_break_3edge_at_plane() hit_pt equal to v1mate\n");
843 
844  {
845  vect_t va, vb;
846  VSUB2(va, hit_pt, eu1->vu_p->v_p->vg_p->coord);
847  VSUB2(vb, eu1->eumate_p->vu_p->v_p->vg_p->coord, hit_pt);
848  VUNITIZE(va);
849  VUNITIZE(vb);
850  if (VDOT(va, vb) < M_SQRT1_2) {
851  bu_bomb("nmg_break_3edge_at_plane() eu1 changes direction?\n");
852  }
853  }
854  {
855  struct bn_tol t2;
856  t2 = is->tol; /* Struct copy */
857 
858  t2.dist = is->tol.dist * 4;
859  t2.dist_sq = t2.dist * t2.dist;
860  dist = DIST_PT_PT(hit_pt, v1->vg_p->coord);
861  if (bn_pt3_pt3_equal(hit_pt, v1->vg_p->coord, &t2))
862  bu_log("NOTICE: nmg_break_3edge_at_plane() hit_pt nearly equal to v1 %g*tol\n", dist/is->tol.dist);
863  dist = DIST_PT_PT(hit_pt, v1mate->vg_p->coord);
864  if (bn_pt3_pt3_equal(hit_pt, v1mate->vg_p->coord, &t2))
865  bu_log("NOTICE: nmg_break_3edge_at_plane() hit_pt nearly equal to v1mate %g*tol\n", dist/is->tol.dist);
866  }
867 
868  /* if we can't find the appropriate vertex in the
869  * other face by a geometry search, build a new vertex.
870  * Otherwise, re-use the existing one.
871  * Can't just search other face, might miss relevant vert.
872  */
873  v2 = nmg_find_pt_in_model(fu2->s_p->r_p->m_p, hit_pt, &(is->tol));
874  if (v2) {
875  /* the other face has a convenient vertex for us */
876  if (RTG.NMG_debug & DEBUG_POLYSECT)
877  bu_log("re-using vertex v=%p from other shell\n", (void *)v2);
878 
879  eu1forw = nmg_ebreaker(v2, eu1, &(is->tol));
880  vu1_final = eu1forw->vu_p;
881  vu2_final = nmg_enlist_vu(is, vu1_final, 0, MAX_FASTF);
882  } else {
883  /* The other face has no vertex in this vicinity */
884  /* If hit_pt falls outside all the loops in fu2,
885  * then there is no need to break this edge.
886  * XXX It is probably cheaper to call nmg_isect_3vertex_3face()
887  * XXX here first, causing any ON cases to be resolved into
888  * XXX shared topology first (and also cutting fu2 edges NOW),
889  * XXX and then run the classifier to answer IN/OUT.
890  * This is expensive. For getting started, tolerate it.
891  */
892  int nmg_class;
893  nmg_class = nmg_class_pt_fu_except(hit_pt, fu2,
894  (struct loopuse *)NULL,
895  (void (*)(struct edgeuse *, point_t, const char *))NULL,
896  (void (*)(struct vertexuse *, point_t, const char *))NULL,
897  (const char *)NULL, 0,
898  0, &is->tol);
899 
900  eu1forw = nmg_ebreaker((struct vertex *)NULL, eu1, &is->tol);
901  vu1_final = eu1forw->vu_p;
902  nmg_vertex_gv(vu1_final->v_p, hit_pt);
903  if (RTG.NMG_debug & DEBUG_POLYSECT)
904  bu_log("Made new vertex vu=%p, v=%p\n",
905  (void *)vu1_final, (void *)vu1_final->v_p);
906 
907  NMG_CK_VERTEX_G(eu1->vu_p->v_p->vg_p);
908  NMG_CK_VERTEX_G(eu1->eumate_p->vu_p->v_p->vg_p);
909  NMG_CK_VERTEX_G(eu1forw->vu_p->v_p->vg_p);
910  NMG_CK_VERTEX_G(eu1forw->eumate_p->vu_p->v_p->vg_p);
911 
912  if (RTG.NMG_debug & DEBUG_POLYSECT) {
913  register pointp_t p1 = eu1->vu_p->v_p->vg_p->coord;
914  register pointp_t p2 = eu1->eumate_p->vu_p->v_p->vg_p->coord;
915 
916  bu_log("After split eu1 %p= %g, %g, %g -> %g, %g, %g\n",
917  (void *)eu1,
918  V3ARGS(p1), V3ARGS(p2));
919  p1 = eu1forw->vu_p->v_p->vg_p->coord;
920  p2 = eu1forw->eumate_p->vu_p->v_p->vg_p->coord;
921  bu_log("\teu1forw %p = %g, %g, %g -> %g, %g, %g\n",
922  (void *)eu1forw,
923  V3ARGS(p1), V3ARGS(p2));
924  }
925 
926  switch (nmg_class) {
927  case NMG_CLASS_AinB:
928  /* point inside a face loop, break edge */
929  break;
930  case NMG_CLASS_AonBshared:
931  /* point is on a loop boundary. Break fu2 loop too? */
932  if (RTG.NMG_debug & DEBUG_POLYSECT)
933  bu_log("%%%%%% point is on loop boundary. Break fu2 loop too?\n");
934  nmg_isect_3vertex_3face(is, vu1_final, fu2);
935  /* XXX should get new vu2 from isect_3vertex_3face! */
936  vu2_final = nmg_find_v_in_face(vu1_final->v_p, fu2);
937  if (!vu2_final) bu_bomb("%%%%%% missed!\n");
938  NMG_CK_VERTEXUSE(vu2_final);
939  nmg_enlist_vu(is, vu1_final, vu2_final, MAX_FASTF);
940  return vu2_final;
941  case NMG_CLASS_AoutB:
942  /* Can't optimize this, break edge anyway. */
943  break;
944  default:
945  bu_bomb("nmg_break_3edge_at_plane() bad classification return from nmg_class_pt_f()\n");
946  }
947 
948  /* stick this vertex in the other shell
949  * and make sure it is in the other shell's
950  * list of vertices on the intersect line
951  */
952  plu2 = nmg_mlv(&fu2->l.magic, vu1_final->v_p, OT_BOOLPLACE);
953  vu2_final = BU_LIST_FIRST(vertexuse, &plu2->down_hd);
954  NMG_CK_VERTEXUSE(vu2_final);
955  nmg_loop_g(plu2->l_p, &is->tol);
956 
957  if (RTG.NMG_debug & DEBUG_POLYSECT) {
958  bu_log("Made vertexloop in other face. lu=%p vu=%p on v=%p\n",
959  (void *)plu2,
960  (void *)vu2_final, (void *)vu2_final->v_p);
961  }
962  vu2_final = nmg_enlist_vu(is, vu1_final, vu2_final, MAX_FASTF);
963  }
964 
965  if (RTG.NMG_debug & DEBUG_POLYSECT) {
966  register pointp_t p1, p2;
967  p1 = eu1->vu_p->v_p->vg_p->coord;
968  p2 = eu1->eumate_p->vu_p->v_p->vg_p->coord;
969  bu_log("\tNow %g, %g, %g <-> %g, %g, %g\n",
970  V3ARGS(p1), V3ARGS(p2));
971  p1 = eu1forw->vu_p->v_p->vg_p->coord;
972  p2 = eu1forw->eumate_p->vu_p->v_p->vg_p->coord;
973  bu_log("\tand %g, %g, %g <-> %g, %g, %g\n\n",
974  V3ARGS(p1), V3ARGS(p2));
975  }
976  return vu2_final;
977 }
978 
979 
980 /**
981  * The vertex 'v2' is known to lie in the plane of eu1's face.
982  * If v2 lies between the two endpoints of eu1, break eu1 and
983  * return the new edgeuse pointer.
984  *
985  * If an edgeuse vertex is joined with v2, v2 remains as the survivor,
986  * as the caller is working on it explicitly, and the edgeuse vertices
987  * are dealt with implicitly (by dereferencing the eu pointers).
988  * Otherwise, we will invalidate our caller's v2 pointer.
989  *
990  * Note that no "intersection line" stuff is done, the goal here is
991  * just to get the edge appropriately broken.
992  *
993  * Either faceuse can be passed in, but it needs to be consistent with the
994  * faceuse used to establish the 2d vertex cache.
995  *
996  * Returns -
997  * new_eu if edge is broken
998  * 0 otherwise
999  */
1000 struct edgeuse *
1001 nmg_break_eu_on_v(struct edgeuse *eu1, struct vertex *v2, struct faceuse *fu, struct nmg_inter_struct *is)
1002 
1003 
1004 /* for plane equation of (either) face */
1005 
1006 {
1007  point_t a;
1008  point_t b;
1009  point_t p;
1010  int code;
1011  fastf_t dist;
1012  struct vertex *v1a;
1013  struct vertex *v1b;
1014  struct edgeuse *new_eu = (struct edgeuse *)0;
1015 
1016  NMG_CK_EDGEUSE(eu1);
1017  NMG_CK_VERTEX(v2);
1018  NMG_CK_FACEUSE(fu);
1019  NMG_CK_INTER_STRUCT(is);
1020 
1021  v1a = eu1->vu_p->v_p;
1022  v1b = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p->v_p;
1023 
1024  /* Check for already shared topology */
1025  if (v1a == v2 || v1b == v2) {
1026  goto out;
1027  }
1028 
1029  /* Map to 2d */
1030  nmg_get_2d_vertex(a, v1a, is, &fu->l.magic);
1031  nmg_get_2d_vertex(b, v1b, is, &fu->l.magic);
1032  nmg_get_2d_vertex(p, v2, is, &fu->l.magic);
1033 
1034  dist = -INFINITY;
1035  code = bn_isect_pt2_lseg2(&dist, a, b, p, &(is->tol));
1036 
1037  switch (code) {
1038  case -2:
1039  /* P outside AB */
1040  break;
1041  default:
1042  case -1:
1043  /* P not on line */
1044  /* This can happen when v2 is a long way from the lseg */
1045  break;
1046  case 1:
1047  /* P is at A */
1048  nmg_jv(v2, v1a); /* v2 must be surviving vertex */
1049  break;
1050  case 2:
1051  /* P is at B */
1052  nmg_jv(v2, v1b); /* v2 must be surviving vertex */
1053  break;
1054  case 3:
1055  /* P is in the middle, break edge */
1056  new_eu = nmg_ebreaker(v2, eu1, &is->tol);
1057  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1058  bu_log("nmg_break_eu_on_v() breaking eu=%p on v=%p, new_eu=%p\n",
1059  (void *)eu1, (void *)v2, (void *)new_eu);
1060  }
1061  break;
1062  }
1063 
1064 out:
1065  return new_eu;
1066 }
1067 
1068 
1069 /**
1070  * Given a vertex 'v' which is already known to have geometry that lies
1071  * on the line defined by 'eg', break all the edgeuses along 'eg'
1072  * which cross 'v'.
1073  *
1074  * Calculation is done in 1 dimension: parametric distance along 'eg'.
1075  * Edge direction vector needs to be made unit length so that tol->dist
1076  * makes sense.
1077  */
1078 void
1079 nmg_break_eg_on_v(const struct edge_g_lseg *eg, struct vertex *v, const struct bn_tol *tol)
1080 {
1081  register struct edgeuse **eup;
1082  struct bu_ptbl eutab;
1083  vect_t dir;
1084  double vdist;
1085 
1086  NMG_CK_EDGE_G_LSEG(eg);
1087  NMG_CK_VERTEX(v);
1088  BN_CK_TOL(tol);
1089 
1090  VMOVE(dir, eg->e_dir);
1091  VUNITIZE(dir);
1092  vdist = bn_dist_pt3_along_line3(eg->e_pt, dir, v->vg_p->coord);
1093 
1094  /* This has to be a table, because nmg_ebreaker() will
1095  * change the list on the fly, otherwise.
1096  */
1097  nmg_edgeuse_with_eg_tabulate(&eutab, eg);
1098 
1099  for (eup = (struct edgeuse **)BU_PTBL_LASTADDR(&eutab);
1100  eup >= (struct edgeuse **)BU_PTBL_BASEADDR(&eutab);
1101  eup--
1102  ) {
1103  struct vertex *va;
1104  struct vertex *vb;
1105  double a;
1106  double b;
1107  struct edgeuse *new_eu;
1108 
1109  NMG_CK_EDGEUSE(*eup);
1110  if ((*eup)->g.lseg_p != eg) bu_bomb("nmg_break_eg_on_v() eu disowns eg\n");
1111 
1112  va = (*eup)->vu_p->v_p;
1113  vb = (*eup)->eumate_p->vu_p->v_p;
1114  if (v == va || v == vb) continue;
1115  if (bn_pt3_pt3_equal(v->vg_p->coord, va->vg_p->coord, tol)) {
1116  nmg_jv(v, va);
1117  continue;
1118  }
1119  if (bn_pt3_pt3_equal(v->vg_p->coord, vb->vg_p->coord, tol)) {
1120  nmg_jv(v, vb);
1121  continue;
1122  }
1123  a = bn_dist_pt3_along_line3(eg->e_pt, dir, va->vg_p->coord);
1124  b = bn_dist_pt3_along_line3(eg->e_pt, dir, vb->vg_p->coord);
1125  if (NEAR_EQUAL(a, vdist, tol->dist)) continue;
1126  if (NEAR_EQUAL(b, vdist, tol->dist)) continue;
1127  if (!bn_between(a, vdist, b, tol)) continue;
1128  new_eu = nmg_ebreaker(v, *eup, tol);
1129  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1130  bu_log("nmg_break_eg_on_v(eg=%p, v=%p) new_eu=%p\n",
1131  (void *)eg, (void *)v, (void *)new_eu);
1132  }
1133  }
1134  bu_ptbl_free(&eutab);
1135 }
1136 
1137 
1138 /**
1139  * Perform edge mutual breaking only on two collinear edgeuses.
1140  * This can result in 2 new edgeuses showing up in either loop (case A & D).
1141  * The vertexuse lists are updated to have all participating vu's and
1142  * their duals.
1143  *
1144  * Two collinear line segments (eu1 and eu2, or just "1" and "2" in the
1145  * diagram) can overlap each other in one of 9 configurations,
1146  * labeled A through I:
1147  *
1148  * A B C D E F G H I
1149  *
1150  * vu1b, vu2b
1151  * * * * * * * * * *=*
1152  * 1 1 2 2 1 2 1 2 1 2
1153  * 1=* 1 2 *=2 1=* *=2 * * 1 2
1154  * 1 2 *=* *=* 1 2 1 2 1 2 1 2
1155  * 1 2 2 1 1 2 1 2 1 2 * * 1 2
1156  * 1=* 2 1 *=2 *=2 1=* 2 1 1 2
1157  * 1 * * 2 2 1 * * 1 2
1158  * * * * * *=*
1159  * vu1a, vu2a
1160  *
1161  * To ensure nothing is missed, break every edgeuse on all 4 vertices.
1162  * If a new edgeuse is created, add it to the list of edgeuses still to be
1163  * broken.
1164  * Brute force, but *certain* not to miss anything.
1165  *
1166  * There is nothing to prevent eu1 and eu2 from being edgeuses in the same
1167  * loop. This creates interesting patterns if one is NEXT of the other,
1168  * such as vu[1] == vu[2]. Just handle it gracefully.
1169  *
1170  * Returns the number of edgeuses that resulted,
1171  * which is always at least the original 2.
1172  *
1173  */
1174 int
1175 nmg_isect_2colinear_edge2p(struct edgeuse *eu1, struct edgeuse *eu2, struct faceuse *fu, struct nmg_inter_struct *is, struct bu_ptbl *l1, struct bu_ptbl *l2)
1176 
1177 
1178 /* for plane equation of (either) face */
1179 
1180 /* optional: list of new eu1 pieces */
1181 /* optional: list of new eu2 pieces */
1182 {
1183  struct edgeuse *eu[10];
1184  struct vertexuse *vu[4];
1185  register int i;
1186  register int j;
1187  int neu; /* Number of edgeuses */
1188 
1189  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1190  bu_log("nmg_isect_2colinear_edge2p(eu1=%p, eu2=%p) START\n",
1191  (void *)eu1, (void *)eu2);
1192  }
1193 
1194  NMG_CK_EDGEUSE(eu1);
1195  NMG_CK_EDGEUSE(eu2);
1196  NMG_CK_FACEUSE(fu); /* Don't check it, just pass it on down. */
1197  NMG_CK_INTER_STRUCT(is);
1198  if (l1) BU_CK_PTBL(l1);
1199  if (l2) BU_CK_PTBL(l2);
1200 
1201  vu[0] = eu1->vu_p;
1202  vu[1] = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
1203  vu[2] = eu2->vu_p;
1204  vu[3] = BU_LIST_PNEXT_CIRC(edgeuse, eu2)->vu_p;
1205 
1206  eu[0] = eu1;
1207  eu[1] = eu2;
1208  neu = 2;
1209 
1210  for (i = 0; i < neu; i++) {
1211  for (j = 0; j < 4; j++) {
1212  eu[neu] = nmg_break_eu_on_v(eu[i], vu[j]->v_p, fu, is);
1213  if (eu[neu]) {
1214  nmg_enlist_vu(is, eu[neu]->vu_p, vu[j], MAX_FASTF);
1215  if (l1 && eu[neu]->e_p == eu1->e_p)
1216  bu_ptbl_ins_unique(l1, (long *)&eu[neu]->l.magic);
1217  else if (l2 && eu[neu]->e_p == eu2->e_p)
1218  bu_ptbl_ins_unique(l2, (long *)&eu[neu]->l.magic);
1219  neu++;
1220  }
1221  }
1222  }
1223 
1224  /* Now join 'em up */
1225  /* This step should no longer be necessary, as nmg_ebreaker()
1226  * from nmg_break_eu_on_v() should have already handled this. */
1227  for (i=0; i < neu-1; i++) {
1228  for (j=i+1; j < neu; j++) {
1229  if (!NMG_ARE_EUS_ADJACENT(eu[i], eu[j])) continue;
1230  nmg_radial_join_eu(eu[i], eu[j], &(is->tol));
1231  }
1232  }
1233 
1234  /* Enlist all four of the original endpoints */
1235  for (i=0; i < 4; i++) {
1236  for (j=0; j < 4; j++) {
1237  if (i==j) continue;
1238  if (vu[i] == vu[j]) continue; /* Happens if eu2 follows eu1 in loop */
1239  if (vu[i]->v_p == vu[j]->v_p) {
1240  nmg_enlist_vu(is, vu[i], vu[j], MAX_FASTF);
1241  goto next_i;
1242  }
1243  }
1244  /* No match, let subroutine hunt for dual */
1245  nmg_enlist_vu(is, vu[i], 0, MAX_FASTF);
1246  next_i: ;
1247  }
1248 
1249  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1250  bu_log("nmg_isect_2colinear_edge2p(eu1=%p, eu2=%p) ret #eu=%d\n",
1251  (void *)eu1, (void *)eu2, neu);
1252  }
1253  return neu;
1254 }
1255 
1256 
1257 /**
1258  * Actual 2d edge/edge intersector
1259  *
1260  * One or both of the edges may be wire edges, i.e.
1261  * either or both of the fu1 and fu2 args may be null.
1262  * If so, the vert_list's are unimport4ant.
1263  *
1264  * Returns a bit vector -
1265  * ISECT_NONE no intersection
1266  * ISECT_SHARED_V intersection was at (at least one) shared vertex
1267  * ISECT_SPLIT1 eu1 was split at (geometric) intersection.
1268  * ISECT_SPLIT2 eu2 was split at (geometric) intersection.
1269  */
1270 int
1271 nmg_isect_edge2p_edge2p(struct nmg_inter_struct *is, struct edgeuse *eu1, struct edgeuse *eu2, struct faceuse *fu1, struct faceuse *fu2)
1272 
1273 
1274 /* fu of eu1, for plane equation */
1275 /* fu of eu2, for error checks */
1276 {
1277  point_t eu1_start;
1278  point_t eu1_end;
1279  vect_t eu1_dir;
1280  point_t eu2_start;
1281  point_t eu2_end;
1282  vect_t eu2_dir;
1283  vect_t dir3d;
1284  fastf_t dist[2];
1285  int status;
1286  point_t hit_pt;
1287  struct vertexuse *vu;
1288  struct vertexuse *vu1a, *vu1b;
1289  struct vertexuse *vu2a, *vu2b;
1290  struct model *m;
1291  int ret = 0;
1292 
1293  NMG_CK_INTER_STRUCT(is);
1294  NMG_CK_EDGEUSE(eu1);
1295  NMG_CK_EDGEUSE(eu2);
1296  m = nmg_find_model(&eu1->l.magic);
1297  NMG_CK_MODEL(m);
1298  /*
1299  * Important note: don't use eu1->eumate_p->vu_p here,
1300  * because that vu is in the opposite orientation faceuse.
1301  * Putting those vu's on the intersection line makes for big trouble.
1302  */
1303  vu1a = eu1->vu_p;
1304  vu1b = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
1305  vu2a = eu2->vu_p;
1306  vu2b = BU_LIST_PNEXT_CIRC(edgeuse, eu2)->vu_p;
1307  NMG_CK_VERTEXUSE(vu1a);
1308  NMG_CK_VERTEXUSE(vu1b);
1309  NMG_CK_VERTEXUSE(vu2a);
1310  NMG_CK_VERTEXUSE(vu2b);
1311 
1312  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1313  bu_log("nmg_isect_edge2p_edge2p(eu1=%p, eu2=%p) START\n\tfu1=%p, fu2=%p\n\tvu1a=%p vu1b=%p, vu2a=%p vu2b=%p\n\tv1a=%p v1b=%p, v2a=%p v2b=%p\n",
1314  (void *)eu1, (void *)eu2,
1315  (void *)fu1, (void *)fu2,
1316  (void *)vu1a, (void *)vu1b, (void *)vu2a, (void *)vu2b,
1317  (void *)vu1a->v_p, (void *)vu1b->v_p, (void *)vu2a->v_p, (void *)vu2b->v_p);
1318  }
1319 
1320  /*
1321  * Topology check.
1322  * If both endpoints of both edges match, this is a trivial accept.
1323  */
1324  if (vu1a->v_p == vu2a->v_p && vu1b->v_p == vu2b->v_p) {
1325  if (RTG.NMG_debug & DEBUG_POLYSECT)
1326  bu_log("nmg_isect_edge2p_edge2p: shared edge topology, both ends\n");
1327  nmg_radial_join_eu(eu1, eu2, &is->tol);
1328  nmg_enlist_vu(is, vu1a, vu2a, MAX_FASTF);
1329  nmg_enlist_vu(is, vu1b, vu2b, MAX_FASTF);
1330  ret = ISECT_SHARED_V;
1331  goto out; /* vu1a, vu1b already listed */
1332  }
1333  if (vu1a->v_p == vu2b->v_p && vu1b->v_p == vu2a->v_p) {
1334  if (RTG.NMG_debug & DEBUG_POLYSECT)
1335  bu_log("nmg_isect_edge2p_edge2p: shared edge topology, both ends, reversed.\n");
1336  nmg_radial_join_eu(eu1, eu2, &is->tol);
1337  nmg_enlist_vu(is, vu1a, vu2b, MAX_FASTF);
1338  nmg_enlist_vu(is, vu1b, vu2a, MAX_FASTF);
1339  ret = ISECT_SHARED_V;
1340  goto out; /* vu1a, vu1b already listed */
1341  }
1342 
1343  /*
1344  * The 3D line in is->pt and is->dir is prepared by the caller.
1345  * is->pt is *not* one of the endpoints of this edge.
1346  *
1347  * IMPORTANT NOTE: The edge-ray used for the edge intersection
1348  * calculations is collinear with the "intersection line",
1349  * but the edge-ray starts at vu1a and points to vu1b,
1350  * while the intersection line has to satisfy different constraints.
1351  * Don't confuse the two!
1352  */
1353  nmg_get_2d_vertex(eu1_start, vu1a->v_p, is, &fu2->l.magic); /* 2D line */
1354  nmg_get_2d_vertex(eu1_end, vu1b->v_p, is, &fu2->l.magic);
1355  V2SUB2(eu1_dir, eu1_end, eu1_start);
1356 
1357  nmg_get_2d_vertex(eu2_start, vu2a->v_p, is, &fu2->l.magic);
1358  nmg_get_2d_vertex(eu2_end, vu2b->v_p, is, &fu2->l.magic);
1359  V2SUB2(eu2_dir, eu2_end, eu2_start);
1360 
1361  dist[0] = dist[1] = 0; /* for clean prints, below */
1362 
1363  /* The "proper" thing to do is intersect two line segments.
1364  * However, this means that none of the intersections of edge "line"
1365  * with the exterior of the loop are computed, and that
1366  * violates the strategy assumptions of the face-cutter.
1367  */
1368  /* To pick up ALL intersection points, the source edge is a line */
1369  status = bn_isect_line2_lseg2(dist, eu1_start, eu1_dir,
1370  eu2_start, eu2_dir, &is->tol);
1371 
1372  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1373  bu_log("\tbn_isect_line2_lseg2()=%d, dist: %g, %g\n",
1374  status, dist[0], dist[1]);
1375  }
1376 
1377  /*
1378  * Whether geometry hits or misses, as long as not collinear, check topo.
1379  * If one endpoint matches, and edges are not collinear,
1380  * then accept the one shared vertex as the intersection point.
1381  * Can't do this before geometry check, or we might miss the
1382  * collinear condition, and not do the mutual intersection.
1383  */
1384  if (status != 0 &&
1385  (vu1a->v_p == vu2a->v_p || vu1a->v_p == vu2b->v_p ||
1386  vu1b->v_p == vu2a->v_p || vu1b->v_p == vu2b->v_p)
1387  ) {
1388  if (RTG.NMG_debug & DEBUG_POLYSECT)
1389  bu_log("edge2p_edge2p: non-colinear edges share one vertex (topology)\n");
1390  if (vu1a->v_p == vu2a->v_p)
1391  nmg_enlist_vu(is, vu1a, vu2a, MAX_FASTF);
1392  else if (vu1a->v_p == vu2b->v_p)
1393  nmg_enlist_vu(is, vu1a, vu2b, MAX_FASTF);
1394 
1395  if (vu1b->v_p == vu2a->v_p)
1396  nmg_enlist_vu(is, vu1b, vu2a, MAX_FASTF);
1397  else if (vu1b->v_p == vu2b->v_p)
1398  nmg_enlist_vu(is, vu1b, vu2b, MAX_FASTF);
1399 
1400  ret = ISECT_SHARED_V;
1401  goto out; /* vu1a, vu1b already listed */
1402  }
1403 
1404  if (status < 0) {
1405  ret = ISECT_NONE; /* No geometric intersection */
1406  goto topo; /* Still need to list vu1a, vu2b */
1407  }
1408 
1409  if (status == 0) {
1410  /* Lines are co-linear and on line of intersection. */
1411  /* Perform full mutual intersection, and vu enlisting. */
1412  if (nmg_isect_2colinear_edge2p(eu1, eu2, fu2, is, (struct bu_ptbl *)0, (struct bu_ptbl *)0) > 2) {
1413  /* Can't tell which edgeuse(s) got split */
1414  ret = ISECT_SPLIT1 | ISECT_SPLIT2;
1415  } else {
1416  /* XXX Can't tell if some sharing ensued. Does it matter? */
1417  /* No, not for the one place we are called. */
1418  ret = ISECT_NONE;
1419  }
1420  goto out; /* vu1a, vu1b listed by nmg_isect_2colinear_edge2p */
1421  }
1422 
1423  /* There is only one intersect point. Break one or both edges. */
1424 
1425 
1426  /* The ray defined by the edgeuse line eu1 intersects the lseg eu2.
1427  * Tolerances have already been factored in.
1428  * The edge exists over values of 0 <= dist <= 1.
1429  */
1430  VSUB2(dir3d, vu1b->v_p->vg_p->coord, vu1a->v_p->vg_p->coord);
1431  VJOIN1(hit_pt, vu1a->v_p->vg_p->coord, dist[0], dir3d);
1432 
1433  if (ZERO(dist[0])) {
1434  /* First point of eu1 is on eu2, by geometry */
1435  if (RTG.NMG_debug & DEBUG_POLYSECT)
1436  bu_log("\tvu=%p vu1a is intersect point\n", (void *)vu1a);
1437  if (dist[1] < 0 || dist[1] > 1) {
1438  if (RTG.NMG_debug & DEBUG_POLYSECT)
1439  bu_log("\teu1 line intersects eu2 outside vu2a...vu2b range, ignore.\n");
1440  ret = ISECT_NONE;
1441  goto topo;
1442  }
1443 
1444  /* Edges not collinear. Either join up with a matching vertex,
1445  * or break eu2 on our vert.
1446  */
1447  if (ZERO(dist[1])) {
1448  if (RTG.NMG_debug & DEBUG_POLYSECT)
1449  bu_log("\tvu2a matches vu1a\n");
1450  nmg_jv(vu1a->v_p, vu2a->v_p);
1451  nmg_enlist_vu(is, vu1a, vu2a, MAX_FASTF);
1452  ret = ISECT_SHARED_V;
1453  goto topo;
1454  }
1455  if (ZERO(dist[1] - 1.0)) {
1456  if (RTG.NMG_debug & DEBUG_POLYSECT)
1457  bu_log("\tsecond point of eu2 matches vu1a\n");
1458  nmg_jv(vu1a->v_p, vu2b->v_p);
1459  nmg_enlist_vu(is, vu1a, vu2b, MAX_FASTF);
1460  ret = ISECT_SHARED_V;
1461  goto topo;
1462  }
1463  /* Break eu2 on our first vertex */
1464  if (RTG.NMG_debug & DEBUG_POLYSECT)
1465  bu_log("\tbreaking eu2 on vu1a\n");
1466  vu = nmg_ebreaker(vu1a->v_p, eu2, &is->tol)->vu_p;
1467  nmg_enlist_vu(is, vu1a, vu, MAX_FASTF);
1468  ret = ISECT_SPLIT2; /* eu1 not broken, just touched */
1469  goto topo;
1470  }
1471 
1472  if (ZERO(dist[0] - 1.0)) {
1473  /* Second point of eu1 is on eu2, by geometry */
1474  if (RTG.NMG_debug & DEBUG_POLYSECT)
1475  bu_log("\tvu=%p vu1b is intersect point\n", (void *)vu1b);
1476  if (dist[1] < 0 || dist[1] > 1) {
1477  if (RTG.NMG_debug & DEBUG_POLYSECT)
1478  bu_log("\teu1 line intersects eu2 outside vu2a...vu2b range, ignore.\n");
1479  ret = ISECT_NONE;
1480  goto topo;
1481  }
1482 
1483  /* Edges not collinear. Either join up with a matching vertex,
1484  * or break eu2 on our vert.
1485  */
1486  if (ZERO(dist[1])) {
1487  if (RTG.NMG_debug & DEBUG_POLYSECT)
1488  bu_log("\tvu2a matches vu1b\n");
1489  nmg_jv(vu1b->v_p, vu2a->v_p);
1490  nmg_enlist_vu(is, vu1b, vu2a, MAX_FASTF);
1491  ret = ISECT_SHARED_V;
1492  goto topo;
1493  }
1494  if (ZERO(dist[1] - 1.0)) {
1495  if (RTG.NMG_debug & DEBUG_POLYSECT)
1496  bu_log("\tsecond point of eu2 matches vu1b\n");
1497  nmg_jv(vu1b->v_p, vu2b->v_p);
1498  nmg_enlist_vu(is, vu1b, vu2b, MAX_FASTF);
1499  ret = ISECT_SHARED_V;
1500  goto topo;
1501  }
1502  /* Break eu2 on our second vertex */
1503  if (RTG.NMG_debug & DEBUG_POLYSECT)
1504  bu_log("\tbreaking eu2 on vu1b\n");
1505  vu = nmg_ebreaker(vu1b->v_p, eu2, &is->tol)->vu_p;
1506  nmg_enlist_vu(is, vu1b, vu, MAX_FASTF);
1507  ret = ISECT_SPLIT2; /* eu1 not broken, just touched */
1508  goto topo;
1509  }
1510 
1511  /* eu2 intersect point is on eu1 line, but not between vertices.
1512  * Since it crosses the line of intersection, it must be broken.
1513  */
1514  if (dist[0] < 0 || dist[0] > 1) {
1515  if (RTG.NMG_debug & DEBUG_POLYSECT)
1516  bu_log("\tIntersect point on eu2 is outside vu1a...vu1b. Break eu2 anyway.\n");
1517 
1518  if (ZERO(dist[1])) {
1519  nmg_enlist_vu(is, vu2a, 0, MAX_FASTF);
1520  ret = ISECT_SHARED_V; /* eu1 was not broken */
1521  goto topo;
1522  } else if (ZERO(dist[1] - 1.0)) {
1523  nmg_enlist_vu(is, vu2b, 0, MAX_FASTF);
1524  ret = ISECT_SHARED_V; /* eu1 was not broken */
1525  goto topo;
1526  } else if (dist[1] > 0 && dist[1] < 1) {
1527  /* Break eu2 somewhere in the middle */
1528  struct vertexuse *new_vu2;
1529  struct vertex *new_v2;
1530  if (RTG.NMG_debug & DEBUG_POLYSECT)
1531  VPRINT("\t\tBreaking eu2 at intersect point", hit_pt);
1532  new_v2 = nmg_find_pt_in_model(m, hit_pt, &(is->tol));
1533  new_vu2 = nmg_ebreaker(new_v2, eu2, &is->tol)->vu_p;
1534  if (!new_v2) {
1535  /* A new vertex was created, assign geom */
1536  nmg_vertex_gv(new_vu2->v_p, hit_pt); /* 3d geom */
1537  }
1538  nmg_enlist_vu(is, new_vu2, 0, MAX_FASTF);
1539  ret = ISECT_SPLIT2; /* eu1 was not broken */
1540  goto topo;
1541  }
1542 
1543  /* Hit point not on either eu1 or eu2, nothing to do */
1544  ret = ISECT_NONE;
1545  goto topo;
1546  }
1547 
1548  /* Intersection is in the middle of the reference edge (eu1) */
1549  /* dist[0] >= 0 && dist[0] <= 1) */
1550  if (RTG.NMG_debug & DEBUG_POLYSECT)
1551  bu_log("\tintersect is in middle of eu1, breaking it\n");
1552 
1553  /* Edges not collinear. Either join up with a matching vertex,
1554  * or break eu2 on our vert.
1555  */
1556  if (ZERO(dist[1])) {
1557  if (RTG.NMG_debug & DEBUG_POLYSECT)
1558  bu_log("\t\tintersect point is vu2a\n");
1559  vu = nmg_ebreaker(vu2a->v_p, eu1, &is->tol)->vu_p;
1560  nmg_enlist_vu(is, vu2a, vu, MAX_FASTF);
1561  ret |= ISECT_SPLIT1;
1562  goto topo;
1563  } else if (ZERO(dist[1] - 1.0)) {
1564  if (RTG.NMG_debug & DEBUG_POLYSECT)
1565  bu_log("\t\tintersect point is vu2b\n");
1566  vu = nmg_ebreaker(vu2b->v_p, eu1, &is->tol)->vu_p;
1567  nmg_enlist_vu(is, vu2b, vu, MAX_FASTF);
1568  ret |= ISECT_SPLIT1;
1569  goto topo;
1570  } else if (dist[1] > 0 && dist[1] < 1) {
1571  /* Intersection is in the middle of both, split edge */
1572  struct vertex *new_v;
1573  if (RTG.NMG_debug & DEBUG_POLYSECT)
1574  VPRINT("\t\tBreaking both edges at intersect point", hit_pt);
1575  ret = ISECT_SPLIT1 | ISECT_SPLIT2;
1576  new_v = nmg_e2break(eu1, eu2);
1577  nmg_vertex_gv(new_v, hit_pt); /* 3d geometry */
1578 
1579  /* new_v is at far end of eu1 and eu2 */
1580  if (eu1->eumate_p->vu_p->v_p != new_v) bu_bomb("new_v 1\n");
1581  if (eu2->eumate_p->vu_p->v_p != new_v) bu_bomb("new_v 2\n");
1582  /* Can't use eumate_p here, it's in wrong orientation face */
1583  nmg_enlist_vu(is, BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p,
1584  BU_LIST_PNEXT_CIRC(edgeuse, eu2)->vu_p, MAX_FASTF);
1585  goto topo;
1586  } else {
1587  /* Intersection is in middle of eu1, which lies on the
1588  * line of intersection being computed, but is outside
1589  * the endpoints of eu2. There is no point in breaking
1590  * eu1 here -- it does not connect up with anything.
1591  */
1592  ret = ISECT_NONE;
1593  goto topo;
1594  }
1595 
1596 topo:
1597  /*
1598  * Listing of any vu's from eu2 will have been done above.
1599  *
1600  * The *original* vu1a and vu1b (and their duals) MUST be
1601  * forcibly listed on
1602  * the intersection line, since eu1 lies ON the line!
1603  *
1604  * This is done last, so that the intersection code (above) has
1605  * the opportunity to create the duals.
1606  * vu1a and vu1b don't have to have anything to do with eu2,
1607  * hence the 2nd vu argument is unspecified (0).
1608  * For our purposes here, we will be satisfied with *any* use
1609  * of the same vertex in the other face.
1610  */
1611  nmg_enlist_vu(is, vu1a, 0, MAX_FASTF);
1612  nmg_enlist_vu(is, vu1b, 0, MAX_FASTF);
1613 out:
1614  /* By here, vu1a and vu1b MUST have been enlisted */
1615  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1616  bu_log("nmg_isect_edge2p_edge2p(eu1=%p, eu2=%p) END, ret=%d %s%s%s\n",
1617  (void *)eu1, (void *)eu2, ret,
1618  (ret&ISECT_SHARED_V)? "SHARED_V|" :
1619  ((ret==0) ? "NONE" : ""),
1620  (ret&ISECT_SPLIT1)? "SPLIT1|" : "",
1621  (ret&ISECT_SPLIT2)? "SPLIT2" : ""
1622  );
1623  }
1624 
1625  return ret;
1626 }
1627 
1628 
1629 /**
1630  * Intersect an edge eu1 with a faceuse fu2.
1631  * eu1 may belong to fu1, or it may be a wire edge.
1632  *
1633  * XXX It is not clear whether we need the caller to provide the
1634  * line equation, or if we should just create it here.
1635  * If done here, the start pt needs to be outside fu2 (fu1 also?)
1636  *
1637  * Returns -
1638  * 0 If everything went well
1639  * 1 If vu[] list along the intersection line needs to be re-done.
1640  */
1641 static int
1642 nmg_isect_wireedge3p_face3p(struct nmg_inter_struct *is, struct edgeuse *eu1, struct faceuse *fu2)
1643 {
1644  struct vertexuse *vu1_final = (struct vertexuse *)NULL;
1645  struct vertexuse *vu2_final = (struct vertexuse *)NULL;
1646  struct vertex *v1a; /* vertex at start of eu1 */
1647  struct vertex *v1b; /* vertex at end of eu1 */
1648  point_t hit_pt;
1649  vect_t edge_vect;
1650  fastf_t edge_len; /* MAGNITUDE(edge_vect) */
1651  fastf_t dist; /* parametric dist to hit point */
1652  fastf_t dist_to_plane; /* distance to hit point, in mm */
1653  int status;
1654  vect_t start_pt;
1655  struct edgeuse *eunext;
1656  struct faceuse *fu1; /* fu that contains eu1 */
1657  plane_t n2;
1658  int ret = 0;
1659 
1660  if (RTG.NMG_debug & DEBUG_POLYSECT)
1661  bu_log("nmg_isect_wireedge3p_face3p(, eu1=%p, fu2=%p) START\n", (void *)eu1, (void *)fu2);
1662 
1663  NMG_CK_INTER_STRUCT(is);
1664  NMG_CK_EDGEUSE(eu1);
1665  NMG_CK_VERTEXUSE(eu1->vu_p);
1666  v1a = eu1->vu_p->v_p;
1667  NMG_CK_VERTEX(v1a);
1668  NMG_CK_VERTEX_G(v1a->vg_p);
1669 
1670  NMG_CK_EDGEUSE(eu1->eumate_p);
1671  NMG_CK_VERTEXUSE(eu1->eumate_p->vu_p);
1672  v1b = eu1->eumate_p->vu_p->v_p;
1673  NMG_CK_VERTEX(v1b);
1674  NMG_CK_VERTEX_G(v1b->vg_p);
1675 
1676  NMG_CK_FACEUSE(fu2);
1677  if (fu2->orientation != OT_SAME) bu_bomb("nmg_isect_wireedge3p_face3p() fu2 not OT_SAME\n");
1678  fu1 = nmg_find_fu_of_eu(eu1); /* May be NULL */
1679 
1680  /*
1681  * Form a ray that starts at one vertex of the edgeuse
1682  * and points to the other vertex.
1683  */
1684  VSUB2(edge_vect, v1b->vg_p->coord, v1a->vg_p->coord);
1685  edge_len = MAGNITUDE(edge_vect);
1686 
1687  VMOVE(start_pt, v1a->vg_p->coord);
1688 
1689  {
1690  /* XXX HACK */
1691  double dot;
1692  dot = fabs(VDOT(is->dir, edge_vect) / edge_len) - 1;
1693  if (!NEAR_ZERO(dot, .01)) {
1694  bu_log("HACK HACK cough cough. Resetting is->pt, is->dir\n");
1695  VPRINT("old is->pt ", is->pt);
1696  VPRINT("old is->dir", is->dir);
1697  VMOVE(is->pt, start_pt);
1698  VMOVE(is->dir, edge_vect);
1699  VUNITIZE(is->dir);
1700  VPRINT("new is->pt ", is->pt);
1701  VPRINT("new is->dir", is->dir);
1702  }
1703  }
1704 
1705  NMG_GET_FU_PLANE(n2, fu2);
1706  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1707  bu_log("Testing (%g, %g, %g) -> (%g, %g, %g) dir=(%g, %g, %g)\n",
1708  V3ARGS(start_pt),
1709  V3ARGS(v1b->vg_p->coord),
1710  V3ARGS(edge_vect));
1711  PLPRINT("\t", n2);
1712  }
1713 
1714  status = bn_isect_line3_plane(&dist, start_pt, edge_vect,
1715  n2, &is->tol);
1716 
1717  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1718  if (status >= 0)
1719  bu_log("\tHit. bn_isect_line3_plane=%d, dist=%g (%e)\n",
1720  status, dist, dist);
1721  else
1722  bu_log("\tMiss. Boring status of bn_isect_line3_plane: %d\n",
1723  status);
1724  }
1725  if (status == 0) {
1726  struct nmg_inter_struct is2;
1727 
1728  /*
1729  * Edge (ray) lies in the plane of the other face,
1730  * by geometry. Drop into 2D code to handle all
1731  * possible intersections (there may be many),
1732  * and any cut/joins, then resume with the previous work.
1733  */
1734  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1735  bu_log("nmg_isect_wireedge3p_face3p: edge lies ON face, using 2D code\n@ @ @ @ @ @ @ @ @ @ 2D CODE, START\n");
1736  bu_log(" The status of the face/face intersect line, before 2d:\n");
1737  nmg_pr_ptbl_vert_list("l1", is->l1, is->mag1);
1738  nmg_pr_ptbl_vert_list("l2", is->l2, is->mag2);
1739  }
1740 
1741  is2 = *is; /* make private copy */
1742  is2.vert2d = 0; /* Don't use previously initialized stuff */
1743 
1744  ret = nmg_isect_edge2p_face2p(&is2, eu1, fu2, fu1);
1745 
1746  nmg_isect2d_cleanup(&is2);
1747 
1748  /*
1749  * Because nmg_isect_edge2p_face2p() calls the face cutter,
1750  * vu's in lone lu's that are listed in the current l1 or
1751  * l2 lists may have been destroyed. Its ret is ours.
1752  */
1753 
1754  /* Only do this if list is still OK */
1755  if (RTG.NMG_debug & DEBUG_POLYSECT && ret == 0) {
1756  bu_log("nmg_isect_wireedge3p_face3p: @ @ @ @ @ @ @ @ @ @ 2D CODE, END, resume 3d problem.\n");
1757  bu_log(" The status of the face/face intersect line, so far:\n");
1758  nmg_pr_ptbl_vert_list("l1", is->l1, is->mag1);
1759  nmg_pr_ptbl_vert_list("l2", is->l2, is->mag2);
1760  }
1761 
1762  /* See if start vertex is now shared */
1763  vu2_final=nmg_find_v_in_face(eu1->vu_p->v_p, fu2);
1764  if (vu2_final) {
1765  if (RTG.NMG_debug & DEBUG_POLYSECT)
1766  bu_log("\tEdge start vertex lies on other face (2d topology).\n");
1767  vu1_final = eu1->vu_p;
1768  (void)bu_ptbl_ins_unique(is->l1, (long *)&vu1_final->l.magic);
1769  (void)bu_ptbl_ins_unique(is->l2, (long *)&vu2_final->l.magic);
1770  }
1771  /* XXX HACK HACK -- shut off error checking */
1772  vu1_final = vu2_final = (struct vertexuse *)NULL;
1773  goto out;
1774  }
1775 
1776  /*
1777  * We now know that the edge does not lie +in+ the other face,
1778  * so it will intersect the face in at most one point.
1779  * Before looking at the results of the geometric calculation,
1780  * check the topology. If the topology says that starting vertex
1781  * of this edgeuse is on the other face, that is the hit point.
1782  * Enter the two vertexuses of that starting vertex in the list,
1783  * and return.
1784  *
1785  * XXX Lee wonders if there might be a benefit to violating the
1786  * XXX "only ask geom question once" rule, and doing a geom
1787  * XXX calculation here before the topology check.
1788  */
1789  vu2_final=nmg_find_v_in_face(v1a, fu2);
1790  if (vu2_final) {
1791  vu1_final = eu1->vu_p;
1792  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1793  bu_log("\tEdge start vertex lies on other face (topology).\n\tAdding vu1_final=%p (v=%p), vu2_final=%p (v=%p)\n",
1794  (void *)vu1_final, (void *)vu1_final->v_p,
1795  (void *)vu2_final, (void *)vu2_final->v_p);
1796  }
1797  (void)bu_ptbl_ins_unique(is->l1, (long *)&vu1_final->l.magic);
1798  (void)bu_ptbl_ins_unique(is->l2, (long *)&vu2_final->l.magic);
1799  goto out;
1800  }
1801 
1802  if (status < 0) {
1803  /* Ray does not strike plane.
1804  * See if start point lies on plane.
1805  */
1806  dist = VDOT(start_pt, n2) - n2[W];
1807  if (!NEAR_ZERO(dist, is->tol.dist))
1808  goto out; /* No geometric intersection */
1809 
1810  /* XXX Does this ever happen, now that geom calc is done
1811  * XXX above, and there is 2D handling as well? Lets find out.
1812  */
1813  bu_bomb("nmg_isect_wireedge3p_face3p: Edge start vertex lies on other face (geometry)\n");
1814 
1815  /* Start point lies on plane of other face */
1816  if (RTG.NMG_debug & DEBUG_POLYSECT)
1817  bu_log("\tEdge start vertex lies on other face (geometry)\n");
1818  dist = VSUB2DOT(v1a->vg_p->coord, start_pt, edge_vect)
1819  / edge_len;
1820  }
1821 
1822  /* The ray defined by the edgeuse intersects the plane
1823  * of the other face. Check to see if the distance to
1824  * intersection is between limits of the endpoints of
1825  * this edge(use).
1826  * The edge exists over values of 0 <= dist <= 1, i.e.,
1827  * over values of 0 <= dist_to_plane <= edge_len.
1828  * The tolerance, an absolute distance, can only be compared
1829  * to other absolute distances like dist_to_plane & edge_len.
1830  * The vertices are "fattened" by +/- is->tol units.
1831  */
1832  dist_to_plane = edge_len * dist;
1833 
1834  if (RTG.NMG_debug & DEBUG_POLYSECT)
1835  bu_log("\tedge_len=%g, dist=%g, dist_to_plane=%g\n",
1836  edge_len, dist, dist_to_plane);
1837 
1838  if (dist_to_plane < -is->tol.dist) {
1839  /* Hit is behind first point */
1840  if (RTG.NMG_debug & DEBUG_POLYSECT)
1841  bu_log("\tplane behind first point\n");
1842  goto out;
1843  }
1844 
1845  if (dist_to_plane > edge_len + is->tol.dist) {
1846  if (RTG.NMG_debug & DEBUG_POLYSECT)
1847  bu_log("\tplane beyond second point\n");
1848  goto out;
1849  }
1850 
1851  VJOIN1(hit_pt, start_pt, dist, edge_vect);
1852 
1853  /* Check hit_pt against face/face intersection line */
1854  {
1855  fastf_t ff_dist;
1856  ff_dist = bn_dist_line3_pt3(is->pt, is->dir, hit_pt);
1857  if (ff_dist > is->tol.dist) {
1858  bu_log("WARNING nmg_isect_wireedge3p_face3p() hit_pt off f/f line %g*tol (%e, tol=%e)\n",
1859  ff_dist/is->tol.dist,
1860  ff_dist, is->tol.dist);
1861  /* XXX now what? */
1862  }
1863  }
1864 
1865  /*
1866  * If the vertex on the other end of this edgeuse is on the face,
1867  * then make a linkage to an existing face vertex (if found),
1868  * and give up on this edge, knowing that we'll pick up the
1869  * intersection of the next edgeuse with the face later.
1870  */
1871  if (dist_to_plane < is->tol.dist) {
1872  /* First point is on plane of face, by geometry */
1873  if (RTG.NMG_debug & DEBUG_POLYSECT)
1874  bu_log("\tedge starts at plane intersect\n");
1875  vu1_final = eu1->vu_p;
1876  vu2_final = nmg_enlist_vu(is, vu1_final, 0, MAX_FASTF);
1877  goto out;
1878  }
1879 
1880  if (dist_to_plane < edge_len - is->tol.dist) {
1881  /* Intersection is between first and second vertex points.
1882  * Insert new vertex at intersection point.
1883  */
1884  vu2_final = nmg_break_3edge_at_plane(hit_pt, fu2, is, eu1);
1885  if (vu2_final)
1886  vu1_final = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
1887  goto out;
1888  }
1889 
1890  /* Second point is on plane of face, by geometry */
1891  if (RTG.NMG_debug & DEBUG_POLYSECT)
1892  bu_log("\tedge ends at plane intersect\n");
1893 
1894  eunext = BU_LIST_PNEXT_CIRC(edgeuse, eu1);
1895  NMG_CK_EDGEUSE(eunext);
1896  if (eunext->vu_p->v_p != v1b)
1897  bu_bomb("nmg_isect_wireedge3p_face3p: discontinuous eu loop\n");
1898 
1899  vu1_final = eunext->vu_p;
1900  vu2_final = nmg_enlist_vu(is, vu1_final, 0, MAX_FASTF);
1901 
1902 out:
1903  /* If vu's were added to list, run some quick checks here */
1904  if (vu1_final && vu2_final) {
1905  if (vu1_final->v_p != vu2_final->v_p) bu_bomb("nmg_isect_wireedge3p_face3p() vertex mis-match\n");
1906 
1907  dist = bn_dist_line3_pt3(is->pt, is->dir,
1908  vu1_final->v_p->vg_p->coord);
1909  if (dist > 100*is->tol.dist) {
1910  bu_log("ERROR nmg_isect_wireedge3p_face3p() vu1=%p point off line by %g > 100*dist_tol (%g)\n",
1911  (void *)vu1_final, dist, 100*is->tol.dist);
1912  VPRINT("is->pt|", is->pt);
1913  VPRINT("is->dir", is->dir);
1914  VPRINT(" coord ", vu1_final->v_p->vg_p->coord);
1915  bu_bomb("nmg_isect_wireedge3p_face3p()\n");
1916  }
1917  if (dist > is->tol.dist) {
1918  bu_log("WARNING nmg_isect_wireedge3p_face3p() vu1=%p pt off line %g*tol (%e, tol=%e)\n",
1919  (void *)vu1_final, dist/is->tol.dist,
1920  dist, is->tol.dist);
1921  }
1922  }
1923 
1924  if (RTG.NMG_debug & DEBUG_POLYSECT)
1925  bu_log("nmg_isect_wireedge3p_face3p(, eu1=%p, fu2=%p) ret=%d END\n",
1926  (void *)eu1, (void *)fu2, ret);
1927  return ret;
1928 }
1929 
1930 
1931 /**
1932  * Intersect a single loop with another face.
1933  * Note that it may be a wire loop.
1934  *
1935  * Returns -
1936  * 0 everything is ok
1937  * >0 vu[] list along intersection line needs to be re-done.
1938  */
1939 static int
1940 nmg_isect_wireloop3p_face3p(struct nmg_inter_struct *bs, struct loopuse *lu, struct faceuse *fu)
1941 {
1942  struct edgeuse *eu;
1943  uint32_t magic1;
1944  int discards = 0;
1945 
1946  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1947  plane_t n;
1948  bu_log("nmg_isect_wireloop3p_face3p(, lu=%p, fu=%p) START\n", (void *)lu, (void *)fu);
1949  NMG_GET_FU_PLANE(n, fu);
1950  HPRINT(" fg N", n);
1951  }
1952 
1953  NMG_CK_INTER_STRUCT(bs);
1954  NMG_CK_LOOPUSE(lu);
1955  NMG_CK_LOOP(lu->l_p);
1956  NMG_CK_LOOP_G(lu->l_p->lg_p);
1957 
1958  NMG_CK_FACEUSE(fu);
1959 
1960  magic1 = BU_LIST_FIRST_MAGIC(&lu->down_hd);
1961  if (magic1 == NMG_VERTEXUSE_MAGIC) {
1962  struct vertexuse *vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
1963  /* this is most likely a loop inserted when we split
1964  * up fu2 wrt fu1 (we're now splitting fu1 wrt fu2)
1965  */
1966  nmg_isect_3vertex_3face(bs, vu, fu);
1967  return 0;
1968  } else if (magic1 != NMG_EDGEUSE_MAGIC) {
1969  bu_bomb("nmg_isect_wireloop3p_face3p() Unknown type of NMG loopuse\n");
1970  }
1971 
1972  /* Process loop consisting of a list of edgeuses.
1973  *
1974  * By going backwards around the list we avoid
1975  * re-processing an edgeuse that was just created
1976  * by nmg_isect_wireedge3p_face3p. This is because the edgeuses
1977  * point in the "next" direction, and when one of
1978  * them is split, it inserts a new edge AHEAD or
1979  * "nextward" of the current edgeuse.
1980  */
1981  for (eu = BU_LIST_LAST(edgeuse, &lu->down_hd);
1982  BU_LIST_NOT_HEAD(eu, &lu->down_hd);
1983  eu = BU_LIST_PLAST(edgeuse, eu)) {
1984  NMG_CK_EDGEUSE(eu);
1985 
1986  if (eu->up.magic_p != &lu->l.magic) {
1987  bu_bomb("nmg_isect_wireloop3p_face3p: edge does not share loop\n");
1988  }
1989 
1990  discards += nmg_isect_wireedge3p_face3p(bs, eu, fu);
1991 
1992  nmg_ck_lueu(lu, "nmg_isect_wireloop3p_face3p");
1993  }
1994 
1995  if (RTG.NMG_debug & DEBUG_POLYSECT) {
1996  bu_log("nmg_isect_wireloop3p_face3p(, lu=%p, fu=%p) END, discards=%d\n",
1997  (void *)lu, (void *)fu, discards);
1998  }
1999  return discards;
2000 }
2001 
2002 
2003 /**
2004  * Construct a nice ray for is->pt, is->dir
2005  * which contains the line of intersection, is->on_eg.
2006  *
2007  * See the comment in nmg_isect_two_generics_faces() for details
2008  * on the constraints on this ray, and the algorithm.
2009  *
2010  * XXX Danger?
2011  * The ray -vs- RPP check is being done in 3D.
2012  * It really ought to be done in 2D, to ensure that
2013  * long edge lines on nearly axis-aligned faces don't
2014  * get discarded prematurely!
2015  * XXX Can't just comment out the code, I think the selection
2016  * XXX of is->pt is significant:
2017  * 1) All intersections are at positive distances on the ray,
2018  * 2) dir cross N will point "left".
2019  *
2020  * Returns -
2021  * 0 OK
2022  * 1 ray misses fu2 bounding box
2023  */
2024 int
2025 nmg_isect_construct_nice_ray(struct nmg_inter_struct *is, struct faceuse *fu2)
2026 {
2027  struct xray line;
2028  vect_t invdir;
2029 
2030  NMG_CK_INTER_STRUCT(is);
2031  NMG_CK_FACEUSE(fu2);
2032 
2033  VMOVE(line.r_pt, is->on_eg->e_pt); /* 3D line */
2034  VMOVE(line.r_dir, is->on_eg->e_dir);
2035  VUNITIZE(line.r_dir);
2036  VINVDIR(invdir, line.r_dir);
2037 
2038  /* nmg_loop_g() makes sure there are no 0-thickness faces */
2039  if (!rt_in_rpp(&line, invdir, fu2->f_p->min_pt, fu2->f_p->max_pt)) {
2040  /* The edge ray missed the face RPP, nothing to do. */
2041  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2042  VPRINT("r_pt ", line.r_pt);
2043  VPRINT("r_dir", line.r_dir);
2044  VPRINT("fu2 min", fu2->f_p->min_pt);
2045  VPRINT("fu2 max", fu2->f_p->max_pt);
2046  bu_log("r_min=%g, r_max=%g\n", line.r_min, line.r_max);
2047  bu_log("nmg_isect_construct_nice_ray() edge ray missed face bounding RPP, ret=1\n");
2048  }
2049  return 1; /* Missed */
2050  }
2051  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2052  VPRINT("fu2 min", fu2->f_p->min_pt);
2053  VPRINT("fu2 max", fu2->f_p->max_pt);
2054  bu_log("r_min=%g, r_max=%g\n", line.r_min, line.r_max);
2055  }
2056  /* Start point will lie at min or max dist, outside of face RPP */
2057  VJOIN1(is->pt, line.r_pt, line.r_min, line.r_dir);
2058  if (line.r_min > line.r_max) {
2059  /* Direction is heading the wrong way, flip it */
2060  VREVERSE(is->dir, line.r_dir);
2061  if (RTG.NMG_debug & DEBUG_POLYSECT)
2062  bu_log("flipping dir\n");
2063  } else {
2064  VMOVE(is->dir, line.r_dir);
2065  }
2066  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2067  VPRINT("r_pt ", line.r_pt);
2068  VPRINT("r_dir", line.r_dir);
2069  VPRINT("->pt ", is->pt);
2070  VPRINT("->dir", is->dir);
2071  bu_log("nmg_isect_construct_nice_ray() ret=0\n");
2072  }
2073  return 0;
2074 }
2075 
2076 
2077 /**
2078  * Given one (2D) edge (eu1) lying in the plane of another face (fu2),
2079  * intersect with all the other edges of that face.
2080  * The line of intersection is defined by the geometry of this edgeuse.
2081  * Therefore, all edgeuses in fu1 which share edge geometry are,
2082  * by definition, ON the intersection line. We process all edgeuses
2083  * which share geometry at once, followed by cutjoin operation.
2084  * It is up to the caller not to recall for the other edgeuses of this edge_g.
2085  *
2086  * XXX eu1 may be a wire edge, in which case there is no fu1 face!
2087  *
2088  * Note that this routine completely conducts the
2089  * intersection operation, so that edges may come and go, loops
2090  * may join or split, each time it is called.
2091  * This imposes special requirements on handling the march through
2092  * the linked lists in this routine.
2093  *
2094  * This also means that much of argument "is" is changed each call.
2095  *
2096  * It further means that vu's in lone lu's found along the edge
2097  * "intersection line" here may get merged in, causing the lu to
2098  * be killed, and the vu, which is listed in the 3D (calling)
2099  * routine's l1/l2 list, is now invalid.
2100  *
2101  * NOTE-
2102  * Since this routine calls the face cutter, *all* points of intersection
2103  * along the line, for *both* faces, need to be found.
2104  * Otherwise, the parity requirements of the face cutter will be violated.
2105  * This means that eu1 needs to be intersected with all of fu1 also,
2106  * including itself (so that the vu's at the ends of eu1 are listed).
2107  *
2108  * Returns -
2109  * 0 Topology is completely shared (or no sharing). l1/l2 valid.
2110  * >0 Caller needs to invalidate his l1/l2 list.
2111  */
2112 static int
2113 nmg_isect_edge2p_face2p(struct nmg_inter_struct *is, struct edgeuse *eu1, struct faceuse *fu2, struct faceuse *fu1)
2114 
2115 /* edge to be intersected w/fu2 */
2116 /* face to be intersected w/eu1 */
2117 /* fu that eu1 is from */
2118 {
2119  struct bu_ptbl vert_list1, vert_list2;
2120  fastf_t *mag1, *mag2;
2121  struct vertexuse *vu1;
2122  struct vertexuse *vu2;
2123  struct edgeuse *fu2_eu; /* use of edge in fu2 */
2124  int ret = 0;
2125  struct bu_ptbl eu1_list;
2126  struct bu_ptbl eu2_list;
2127 
2128  NMG_CK_INTER_STRUCT(is);
2129  NMG_CK_EDGEUSE(eu1);
2130  NMG_CK_FACEUSE(fu2);
2131  if (fu1) NMG_CK_FACEUSE(fu1); /* fu1 may be null */
2132 
2133  if (RTG.NMG_debug & DEBUG_POLYSECT)
2134  bu_log("nmg_isect_edge2p_face2p(eu1=%p, fu2=%p, fu1=%p) START\n",
2135  (void *)eu1, (void *)fu2, (void *)fu1);
2136 
2137  if (fu2->orientation != OT_SAME) bu_bomb("nmg_isect_edge2p_face2p() fu2 not OT_SAME\n");
2138  if (fu1 && fu1->orientation != OT_SAME) bu_bomb("nmg_isect_edge2p_face2p() fu1 not OT_SAME\n");
2139 
2140  mag1 = (fastf_t *)NULL;
2141  mag2 = (fastf_t *)NULL;
2142 
2143  /* See if an edge exists in other face that connects these 2 verts */
2144  fu2_eu = nmg_find_eu_in_face(eu1->vu_p->v_p, eu1->eumate_p->vu_p->v_p,
2145  fu2, (const struct edgeuse *)NULL, 0);
2146  if (fu2_eu != (struct edgeuse *)NULL) {
2147  /* There is an edge in other face that joins these 2 verts. */
2148  NMG_CK_EDGEUSE(fu2_eu);
2149  if (fu2_eu->e_p != eu1->e_p) {
2150  /* Not the same edge, fuse! */
2151  bu_log("nmg_isect_edge2p_face2p() fusing unshared shared edge\n");
2152  nmg_radial_join_eu(eu1, fu2_eu, &is->tol);
2153  }
2154  /* Topology is completely shared */
2155  if (RTG.NMG_debug & DEBUG_POLYSECT)
2156  bu_log("nmg_isect_edge2p_face2p() topology is shared\n");
2157  ret = 0;
2158  goto do_ret;
2159  }
2160 
2161  bu_ptbl_init(&vert_list1, 64, "vert_list1 buffer");
2162  bu_ptbl_init(&vert_list2, 64, "vert_list1 buffer");
2163  bu_ptbl_init(&eu1_list, 64, "eu1_list1 buffer");
2164  bu_ptbl_init(&eu2_list, 64, "eu2_list1 buffer");
2165 
2166  NMG_CK_EDGE_G_LSEG(eu1->g.lseg_p);
2167  is->on_eg = eu1->g.lseg_p;
2168  is->l1 = &vert_list1;
2169  is->l2 = &vert_list2;
2170  is->s1 = nmg_find_s_of_eu(eu1); /* may be wire edge */
2171  is->s2 = fu2->s_p;
2172  is->fu1 = fu1;
2173  is->fu2 = fu2;
2174 
2175  if (fu1 && RTG.NMG_debug & (DEBUG_POLYSECT|DEBUG_FCUT|DEBUG_MESH)
2176  && RTG.NMG_debug & DEBUG_PLOTEM) {
2177  nmg_pl_2fu("Iface%d.plot3", fu2, fu1, 0);
2178  }
2179 
2180  vu1 = eu1->vu_p;
2181  vu2 = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
2182  if (vu1->v_p == vu2->v_p) {
2183  bu_log("nmg_isect_edge2p_face2p(eu1=%p) skipping 0-len edge (topology)\n", (void *)eu1);
2184  /* Call nmg_k0eu() ? */
2185  goto out;
2186  }
2187 
2188  /*
2189  * Construct the ray which contains the line of intersection,
2190  * i.e. the line that contains the edge "eu1" (is->on_eg).
2191  */
2192  if (nmg_isect_construct_nice_ray(is, fu2)) goto out;
2193 
2194  if (RTG.NMG_debug & DEBUG_VERIFY) {
2195  nmg_fu_touchingloops(fu2);
2196  if (fu1)nmg_fu_touchingloops(fu1);
2197  nmg_region_v_unique(is->s1->r_p, &is->tol);
2198  nmg_region_v_unique(is->s2->r_p, &is->tol);
2199  }
2200 
2201  /* Build list of all edgeuses in eu1/fu1 and fu2 */
2202  if (fu1) {
2203  nmg_edgeuse_tabulate(&eu1_list, &fu1->l.magic);
2204  } else {
2205  nmg_edgeuse_tabulate(&eu1_list, &eu1->l.magic);
2206  }
2207  nmg_edgeuse_tabulate(&eu2_list, &fu2->l.magic);
2208 
2209  is->mag_len = 2 * (BU_PTBL_END(&eu1_list) + BU_PTBL_END(&eu2_list));
2210  mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag1");
2211  mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag2");
2212 
2213  is->mag1 = mag1;
2214  is->mag2 = mag2;
2215 
2216  /* Run infinite line containing eu1 through fu2 */
2217  nmg_isect_line2_face2pNEW(is, fu2, fu1, &eu2_list, &eu1_list);
2218 
2219  /* If eu1 is a wire, there is no fu1 to run line through. */
2220  if (fu1) {
2221  /* We are intersecting with ourselves */
2222  nmg_isect_line2_face2pNEW(is, fu1, fu2, &eu1_list, &eu2_list);
2223  }
2224 
2225  if (RTG.NMG_debug & DEBUG_FCUT) {
2226  bu_log("nmg_isect_edge2p_face2p(eu1=%p, fu2=%p) vert_lists C:\n", (void *)eu1, (void *)fu2);
2227  nmg_pr_ptbl_vert_list("vert_list1", &vert_list1, mag1);
2228  nmg_pr_ptbl_vert_list("vert_list2", &vert_list2, mag2);
2229  }
2230 
2231  if (RTG.NMG_debug & DEBUG_FCUT) {
2232  bu_log("nmg_isect_edge2p_face2p(eu1=%p, fu2=%p) vert_lists D:\n", (void *)eu1, (void *)fu2);
2233  nmg_pr_ptbl_vert_list("vert_list1", &vert_list1, mag1);
2234  nmg_pr_ptbl_vert_list("vert_list2", &vert_list2, mag2);
2235  }
2236 
2237  if (vert_list1.end == 0 && vert_list2.end == 0) goto out;
2238 
2239  /* Invoke the face cutter to snip and join loops along isect line */
2240  if (fu1 && fu2) {
2241  is->on_eg = nmg_face_cutjoin(&vert_list1, &vert_list2, mag1, mag2, fu1, fu2, is->pt, is->dir, is->on_eg, &is->tol);
2242  ret = 1; /* face cutter was called. */
2243  }
2244 
2245 out:
2246  (void)bu_ptbl_free(&vert_list1);
2247  (void)bu_ptbl_free(&vert_list2);
2248  (void)bu_ptbl_free(&eu1_list);
2249  (void)bu_ptbl_free(&eu2_list);
2250  if (mag1)
2251  bu_free((char *)mag1, "nmg_isect_edge2p_face2p: mag1");
2252  if (mag2)
2253  bu_free((char *)mag2, "nmg_isect_edge2p_face2p: mag2");
2254 
2255 do_ret:
2256  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2257  bu_log("nmg_isect_edge2p_face2p(eu1=%p, fu2=%p) ret=%d\n",
2258  (void *)eu1, (void *)fu2, ret);
2259  }
2260  return ret;
2261 }
2262 
2263 
2264 void
2265 nmg_enlist_one_vu(struct nmg_inter_struct *is, const struct vertexuse *vu, fastf_t dist)
2266 
2267 
2268 /* distance along intersect ray for this vu */
2269 {
2270  struct shell *sv; /* shell of vu */
2271 
2272  NMG_CK_INTER_STRUCT(is);
2273  NMG_CK_VERTEXUSE(vu);
2274 
2275  if (is->mag_len <= BU_PTBL_END(is->l1) || is->mag_len <= BU_PTBL_END(is->l2))
2276  bu_log("Array for distances to vertexuses is too small (%d)\n", is->mag_len);
2277 
2278  sv = nmg_find_s_of_vu(vu);
2279 
2280  /* First step: add vu to corresponding list */
2281  if (sv == is->s1) {
2282  bu_ptbl_ins_unique(is->l1, (long *)&vu->l.magic);
2283  if (is->mag_len <= BU_PTBL_END(is->l1)) {
2284  if (is->mag_len) {
2285  is->mag_len *= 2;
2286  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
2287  "is->mag1");
2288  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
2289  "is->mag2");
2290  } else {
2291  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
2292  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
2293  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
2294  }
2295 
2296  }
2297  if (dist < MAX_FASTF)
2298  is->mag1[bu_ptbl_locate(is->l1, (long *)&vu->l.magic)] = dist;
2299  } else if (sv == is->s2) {
2300  bu_ptbl_ins_unique(is->l2, (long *)&vu->l.magic);
2301  if (is->mag_len <= BU_PTBL_END(is->l2)) {
2302  if (is->mag_len) {
2303  is->mag_len *= 2;
2304  is->mag1 = (fastf_t *)bu_realloc((char *)is->mag1, is->mag_len*sizeof(fastf_t),
2305  "is->mag1");
2306  is->mag2 = (fastf_t *)bu_realloc((char *)is->mag2, is->mag_len*sizeof(fastf_t),
2307  "is->mag2");
2308  } else {
2309  is->mag_len = 2*(BU_PTBL_END(is->l1) + BU_PTBL_END(is->l2));
2310  is->mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag1");
2311  is->mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "is->mag2");
2312  }
2313 
2314  }
2315  if (dist < MAX_FASTF)
2316  is->mag2[bu_ptbl_locate(is->l2, (long *)&vu->l.magic)] = dist;
2317  } else {
2318  bu_log("nmg_enlist_one_vu(vu=%p) sv=%p, s1=%p, s2=%p\n",
2319  (void *)vu, (void *)sv, (void *)is->s1, (void *)is->s2);
2320  bu_bomb("nmg_enlist_one_vu: vu is not in s1 or s2\n");
2321  }
2322 
2323  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2324  bu_log("nmg_enlist_one_vu(vu=%p) v=%p, dist=%g (%s)\n",
2325  (void *)vu, (void *)vu->v_p, dist,
2326  (sv == is->s1) ? "shell 1" : "shell 2");
2327  }
2328 
2329  /* Some (expensive) centralized sanity checking */
2330  if ((RTG.NMG_debug & DEBUG_VERIFY) && is->fu1 && is->fu2) {
2331  nmg_ck_v_in_2fus(vu->v_p, is->fu1, is->fu2, &(is->tol));
2332  }
2333 }
2334 
2335 
2336 static void
2337 nmg_coplanar_face_vertex_fuse(struct faceuse *fu1, struct faceuse *fu2, struct bn_tol *tol)
2338 {
2339  struct bu_ptbl verts;
2340  struct faceuse *faces[4];
2341  struct loopuse *lu;
2342  struct edgeuse *eu;
2343  struct vertexuse *vu;
2344  int i;
2345 
2346  bu_ptbl_init(&verts, 128, "&verts");
2347 
2348  faces[0] = fu1;
2349  faces[1] = fu1->fumate_p;
2350  faces[2] = fu2;
2351  faces[3] = fu2->fumate_p;
2352 
2353  for (i = 0 ; i < 4 ; i++) {
2354  NMG_CK_FACEUSE(faces[i]);
2355  for (BU_LIST_FOR(lu, loopuse, &faces[i]->lu_hd)) {
2356  NMG_CK_LOOPUSE(lu);
2357  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_EDGEUSE_MAGIC) {
2358  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
2359  NMG_CK_EDGEUSE(eu);
2360  NMG_CK_VERTEXUSE(eu->vu_p);
2361  NMG_CK_VERTEX(eu->vu_p->v_p);
2362  (void)bu_ptbl_ins_unique(&verts, (long *)eu->vu_p->v_p);
2363  }
2364  } else if (BU_LIST_FIRST_MAGIC(&lu->down_hd) == NMG_VERTEXUSE_MAGIC) {
2365  vu = BU_LIST_FIRST(vertexuse, &lu->down_hd);
2366  NMG_CK_VERTEXUSE(vu);
2367  NMG_CK_VERTEX(vu->v_p);
2368  (void)bu_ptbl_ins_unique(&verts, (long *)vu->v_p);
2369  } else {
2370  bu_bomb("nmg_coplanar_face_vertex_fuse(): invalid loopuse");
2371  }
2372  }
2373  }
2374 
2375  (void)nmg_vertex_fuse((const uint32_t *)&verts, tol);
2376 
2377  bu_ptbl_free(&verts);
2378 }
2379 
2380 
2381 static void
2382 nmg_isect_two_face2p_jra(struct nmg_inter_struct *is, struct faceuse *fu1, struct faceuse *fu2)
2383 {
2384  struct model *m;
2385  struct loopuse *lu;
2386  struct bu_ptbl eu1_list;
2387  struct bu_ptbl eu2_list;
2388  struct bu_ptbl v_list;
2389  struct bu_ptbl vert_list1, vert_list2;
2390  fastf_t *mag1, *mag2;
2391  int i, j;
2392  const fastf_t opsff = (1.0 + SMALL_FASTF);
2393  const fastf_t omsff = (1.0 - SMALL_FASTF);
2394 
2395  NMG_CK_FACEUSE(fu1);
2396  NMG_CK_FACEUSE(fu2);
2397  NMG_CK_INTER_STRUCT(is);
2398 
2399  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2400  bu_log("nmg_isect_two)face2p_jra: fu1=%p, fu2=%p\n", (void *)fu1, (void *)fu2);
2401  }
2402 
2403  nmg_coplanar_face_vertex_fuse(fu1, fu2, &is->tol);
2404 
2405  m = nmg_find_model(&fu1->l.magic);
2406  NMG_CK_MODEL(m);
2407 
2408  nmg_edgeuse_tabulate(&eu1_list, &fu1->l.magic);
2409  nmg_edgeuse_tabulate(&eu2_list, &fu2->l.magic);
2410 
2411  is->mag_len = 2 * (BU_PTBL_END(&eu1_list) + BU_PTBL_END(&eu2_list));
2412  mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag1");
2413  mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag2");
2414 
2415  for (i = 0; i < is->mag_len; i++) {
2416  mag1[i] = MAX_FASTF;
2417  mag2[i] = MAX_FASTF;
2418  }
2419 
2420  is->s1 = fu1->s_p;
2421  is->s2 = fu2->s_p;
2422  is->fu1 = fu1;
2423  is->fu2 = fu2;
2424  is->mag1 = mag1;
2425  is->mag2 = mag2;
2426 
2427  /* First split all edgeuses that intersect */
2428  for (i = 0; i < BU_PTBL_END(&eu1_list); i++) {
2429  struct edgeuse *eu1;
2430  struct vertex_g *vg1a, *vg1b;
2431  vect_t vt1_3d;
2432 
2433  eu1 = (struct edgeuse *)BU_PTBL_GET(&eu1_list, i);
2434  NMG_CK_EDGEUSE(eu1);
2435 
2436  vg1a = eu1->vu_p->v_p->vg_p;
2437  NMG_CK_VERTEX_G(vg1a);
2438  vg1b = eu1->eumate_p->vu_p->v_p->vg_p;
2439  NMG_CK_VERTEX_G(vg1b);
2440 
2441  VSUB2(vt1_3d, vg1b->coord, vg1a->coord);
2442 
2443  for (j = 0; j < BU_PTBL_END(&eu2_list); j++) {
2444  struct edgeuse *eu2 = NULL;
2445  struct vertex_g *vg2a = NULL;
2446  struct vertex_g *vg2b = NULL;
2447  int code = 0;
2448  vect_t vt2_3d = VINIT_ZERO;
2449  fastf_t dist[2] = V2INIT_ZERO;
2450  point_t hit_pt = VINIT_ZERO;
2451  int hit_no = 0;
2452  int hit_count = 0;
2453 
2454  eu2 = (struct edgeuse *)BU_PTBL_GET(&eu2_list, j);
2455  NMG_CK_EDGEUSE(eu2);
2456 
2457  vg2a = eu2->vu_p->v_p->vg_p;
2458  vg2b = eu2->eumate_p->vu_p->v_p->vg_p;
2459  VSUB2(vt2_3d, vg2b->coord, vg2a->coord);
2460 
2461  code = bn_isect_lseg3_lseg3(dist, vg1a->coord, vt1_3d,
2462  vg2a->coord, vt2_3d, &is->tol);
2463 
2464  if (code < 0)
2465  continue;
2466 
2467  if (code == 0) {
2468  /* When there is only one hit, place the distance in
2469  * dist[0] otherwise use both dist[0] and dist[1].
2470  * In the description below the line p0->p1 is eu1 where
2471  * p0 is the start of eu1 and p1 is the end of eu1.
2472  * The same is done for eu2 with line q0->q1.
2473  * When eu1 and eu2 are collinear, the value of dist[0]
2474  * returned from function 'bn_isect_lseg3_lseg3' is the
2475  * scaled distance from p0->q0 and dist[1] is the scaled
2476  * distance from p0->q1.
2477  */
2478  if ((dist[0] < -SMALL_FASTF && ((dist[1] > SMALL_FASTF) &&
2479  (dist[1] < omsff))) ||
2480  ((dist[1] > SMALL_FASTF) && (dist[1] < omsff)
2481  && (dist[0] > opsff))) {
2482  /* true when q1 is within p0->p1 */
2483  hit_count = 1;
2484  dist[0] = dist[1];
2485  dist[1] = MAX_FASTF; /* sanity */
2486  } else if ((dist[0] > SMALL_FASTF) && (dist[0] < omsff) &&
2487  (dist[1] > SMALL_FASTF) && (dist[1] < omsff)) {
2488  /* true when both q0 and q1 is within p0->p1 */
2489  hit_count = 2;
2490  /* dist[0] and dist[1] are already the correct values */
2491  } else if (((dist[0] > SMALL_FASTF) && (dist[0] < omsff) && (dist[1] < -SMALL_FASTF)) ||
2492  ((dist[0] > SMALL_FASTF) && (dist[0] < omsff) && (dist[1] > opsff))) {
2493  /* true when q0 is within p0->p1 */
2494  hit_count = 1;
2495  dist[1] = MAX_FASTF; /* sanity */
2496  /* dist[0] is already the correct value */
2497  } else if (((dist[0] < -SMALL_FASTF) && (dist[1] > opsff)) ||
2498  ((dist[0] > opsff) && (dist[1] < -SMALL_FASTF))) {
2499  /* true when both p0 and p1 is within q0->q1 */
2500  /* eu1 is not cut */
2501  hit_count = 0; /* sanity */
2502  dist[0] = dist[1] = MAX_FASTF; /* sanity */
2503  continue;
2504 
2505  } else if ((ZERO(dist[0]) && EQUAL(dist[1], 1.0)) ||
2506  (ZERO(dist[1]) && EQUAL(dist[0], 1.0))) {
2507  /* true when eu1 and eu2 shared the same vertices */
2508  /* eu1 is not cut */
2509  hit_count = 0; /* sanity */
2510  dist[0] = dist[1] = MAX_FASTF; /* sanity */
2511  continue;
2512  } else if (((dist[0] < -SMALL_FASTF) && ZERO(dist[1])) ||
2513  (ZERO(dist[0]) && (dist[1] < -SMALL_FASTF)) ||
2514  (EQUAL(dist[0], 1.0) && (dist[1] > opsff)) ||
2515  (EQUAL(dist[1], 1.0) && (dist[0] > opsff))) {
2516  /* true when eu2 shares one of eu1 vertices and the
2517  * other eu2 vertex is outside p0->p1 (i.e. eu1)
2518  */
2519  /* eu1 is not cut */
2520  hit_count = 0; /* sanity */
2521  dist[0] = dist[1] = MAX_FASTF; /* sanity */
2522  continue;
2523  } else if ((ZERO(dist[0]) && (dist[1] > SMALL_FASTF) &&
2524  (dist[1] < omsff)) ||
2525  ((dist[1] > SMALL_FASTF) &&
2526  (dist[1] < omsff) &&
2527  EQUAL(dist[0], 1.0))) {
2528  /* true when q1 is within p0->p1 and q0 = p0 or q0 = p1 */
2529  hit_count = 1;
2530  dist[0] = dist[1];
2531  dist[1] = MAX_FASTF; /* sanity */
2532  } else if ((ZERO(dist[1]) && (dist[0] > SMALL_FASTF) &&
2533  (dist[0] < omsff)) ||
2534  ((dist[0] > SMALL_FASTF) &&
2535  (dist[0] < omsff) &&
2536  EQUAL(dist[1], 1.0))) {
2537  /* true when q0 is within p0->p1 and q1 = p0 or q1 = p1 */
2538  hit_count = 1;
2539  dist[1] = MAX_FASTF; /* sanity */
2540  /* dist[0] is already the correct value */
2541  } else if ((ZERO(dist[0]) && (dist[1] > opsff)) ||
2542  (ZERO(dist[1]) && (dist[0] > opsff)) ||
2543  (EQUAL(dist[0], 1.0) && (dist[1] < -SMALL_FASTF)) ||
2544  ((dist[0] < -SMALL_FASTF) && EQUAL(dist[1], 1.0))) {
2545  /* true when eu2 shares one of the vertices of eu1 and
2546  * the other vertex in eu2 is on the far side of eu1
2547  * outside eu1 (i.e. p0->p1).
2548  */
2549  /* eu1 is not cut */
2550  hit_count = 0; /* sanity */
2551  dist[0] = dist[1] = MAX_FASTF; /* sanity */
2552  continue;
2553  } else {
2554  /* should never get here */
2555  bu_log("nmg_isect_two_face2p_jra(): dist[0] = %f dist[1] = %f\n",
2556  dist[0], dist[1]);
2557  bu_bomb("nmg_isect_two_face2p_jra(): unexpected condition\n");
2558  }
2559  } else {
2560  if (ZERO(dist[0]) || EQUAL(dist[0], 1.0)) {
2561  /* eu1 was hit on a vertex, nothing to cut */
2562  continue;
2563  } else if (UNLIKELY((dist[0] < -SMALL_FASTF) || (dist[0] > opsff))) {
2564  bu_bomb("nmg_isect_two_face2p_jra(): dist[0] not within 0-1\n");
2565  } else {
2566  hit_count = 1;
2567  }
2568  }
2569 
2570  for (hit_no=0; hit_no < hit_count; hit_no++) {
2571  struct edgeuse *new_eu;
2572  struct vertex *hitv;
2573  struct vertexuse *hit_vu;
2574 
2575  if (dist[hit_no] < -SMALL_FASTF || dist[hit_no] > opsff)
2576  continue;
2577 
2578  hitv = (struct vertex *)NULL;
2579  if (ZERO(dist[hit_no])) {
2580  hitv = eu1->vu_p->v_p;
2581  VMOVE(hit_pt, hitv->vg_p->coord);
2582  } else if (EQUAL(dist[hit_no], 1.0)) {
2583  hitv = eu1->eumate_p->vu_p->v_p;
2584  VMOVE(hit_pt, hitv->vg_p->coord);
2585  } else {
2586  VJOIN1(hit_pt, vg1a->coord, dist[hit_no], vt1_3d);
2587  if ((hit_vu = nmg_find_pt_in_face(fu2, hit_pt, &is->tol))) {
2588  hitv = hit_vu->v_p;
2589  } else if ((hit_vu = nmg_find_pt_in_face(fu1, hit_pt, &is->tol))) {
2590  hitv = hit_vu->v_p;
2591  } else {
2592  hitv = nmg_find_pt_in_model(nmg_find_model(&fu1->l.magic), hit_pt, &is->tol);
2593  }
2594  if (UNLIKELY(!hitv && (code == 0))) {
2595  bu_log("dist[0] = %f dist[1] = %f hit_no = %d\n", dist[0], dist[1], hit_no);
2596  bu_bomb("nmg_isect_two_face2p_jra(): Can not find vertexuse");
2597  }
2598  }
2599  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2600  bu_log("eus %p and %p intersect #%d at (%f %f %f)\n",
2601  (void *)eu1, (void *)eu2, hit_no, V3ARGS(hit_pt));
2602  }
2603  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT && hitv)) {
2604  bu_log("Found vertex (%p) at hit_pt\n", (void *)hitv);
2605  }
2606 
2607  if (hitv != eu1->vu_p->v_p && hitv != eu1->eumate_p->vu_p->v_p) {
2608  struct edgeuse *next_eu, *prev_eu;
2609 
2610  next_eu = BU_LIST_PNEXT_CIRC(edgeuse, &eu1->l);
2611  prev_eu = BU_LIST_PPREV_CIRC(edgeuse, &eu1->l);
2612 
2613  if (hitv != prev_eu->vu_p->v_p && hitv != next_eu->eumate_p->vu_p->v_p) {
2614  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2615  bu_log("Splitting eu1 %p\n", (void *)eu1);
2616  }
2617  new_eu = nmg_esplit(hitv, eu1, 1);
2618  hitv = new_eu->vu_p->v_p;
2619  if (!hitv->vg_p)
2620  nmg_vertex_gv(hitv, hit_pt);
2621  vg1b = eu1->eumate_p->vu_p->v_p->vg_p;
2622  VSUB2(vt1_3d, vg1b->coord, vg1a->coord);
2623  }
2624  }
2625  if (code == 1 && hitv != eu2->vu_p->v_p && hitv != eu2->eumate_p->vu_p->v_p) {
2626  struct edgeuse *next_eu, *prev_eu;
2627 
2628  next_eu = BU_LIST_PNEXT_CIRC(edgeuse, &eu2->l);
2629  prev_eu = BU_LIST_PPREV_CIRC(edgeuse, &eu2->l);
2630 
2631  if (hitv != prev_eu->vu_p->v_p && hitv != next_eu->eumate_p->vu_p->v_p) {
2632  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2633  vect_t tmp1, tmp2;
2634  VSUB2(tmp1, hit_pt, eu2->vu_p->v_p->vg_p->coord);
2635  VSUB2(tmp2, hit_pt, eu2->eumate_p->vu_p->v_p->vg_p->coord);
2636  bu_log("Splitting eu2 %p\n", (void *)eu2);
2637  bu_log("Distance to hit_pt = %g from vu1, %g from vu2\n",
2638  MAGNITUDE(tmp1), MAGNITUDE(tmp2));
2639  }
2640  new_eu = nmg_esplit(hitv, eu2, 1);
2641  hitv = new_eu->vu_p->v_p;
2642  if (!hitv->vg_p)
2643  nmg_vertex_gv(hitv, hit_pt);
2644  bu_ptbl_ins(&eu2_list, (long *)new_eu);
2645  }
2646  }
2647 
2648  if (hitv) {
2649  (void)nmg_break_all_es_on_v(&fu1->l.magic, hitv, &is->tol);
2650  (void)nmg_break_all_es_on_v(&fu2->l.magic, hitv, &is->tol);
2651  }
2652  }
2653  }
2654  }
2655 
2656  bu_ptbl_free(&eu1_list);
2657  bu_ptbl_free(&eu2_list);
2658 
2659  /* Make sure every vertex in fu1 has dual in fu2
2660  * (if they overlap)
2661  */
2662  nmg_vertex_tabulate(&v_list, &fu1->l.magic);
2663 
2664  for (i=0; i<BU_PTBL_END(&v_list); i++) {
2665  struct vertex *v;
2666  int nmg_class;
2667 
2668  v = (struct vertex *)BU_PTBL_GET(&v_list, i);
2669  NMG_CK_VERTEX(v);
2670 
2671  if (nmg_find_v_in_face(v, fu2))
2672  continue;
2673 
2674  /* Check if this vertex is within other FU */
2675  nmg_class = nmg_class_pt_fu_except(v->vg_p->coord, fu2, NULL, NULL, NULL,
2676  (char *)NULL, 0, 0, &is->tol);
2677 
2678  if (nmg_class == NMG_CLASS_AinB) {
2679  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2680  bu_log("Making dualvu of vertex %p in fu2 %p\n", (void *)v, (void *)fu2);
2681  }
2682  (void)nmg_make_dualvu(v, fu2, &is->tol);
2683  }
2684  }
2685  bu_ptbl_reset(&v_list);
2686 
2687  /* same for fu2 */
2688  nmg_vertex_tabulate(&v_list, &fu2->l.magic);
2689 
2690  for (i=0; i<BU_PTBL_END(&v_list); i++) {
2691  struct vertex *v;
2692  int nmg_class;
2693 
2694  v = (struct vertex *)BU_PTBL_GET(&v_list, i);
2695  NMG_CK_VERTEX(v);
2696 
2697  if (nmg_find_v_in_face(v, fu1))
2698  continue;
2699 
2700  /* Check if this vertex is within other FU */
2701  nmg_class = nmg_class_pt_fu_except(v->vg_p->coord, fu1, NULL, NULL, NULL,
2702  (char *)NULL, 0, 0, &is->tol);
2703 
2704  if (nmg_class == NMG_CLASS_AinB) {
2705  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2706  bu_log("Making dualvu of vertex %p in fu1 %p\n", (void *)v, (void *)fu1);
2707  }
2708  (void)nmg_make_dualvu(v, fu1, &is->tol);
2709  }
2710  }
2711 
2712  bu_ptbl_free(&v_list);
2713 
2714  bu_ptbl_init(&vert_list1, 64, " &vert_list1");
2715  bu_ptbl_init(&vert_list2, 64, " &vert_list2");
2716  is->l1 = &vert_list1;
2717  is->l2 = &vert_list2;
2718 
2719  for (BU_LIST_FOR(lu, loopuse, &fu1->lu_hd)) {
2720  struct edgeuse *eu;
2721 
2722  NMG_CK_LOOPUSE(lu);
2723 
2724  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
2725  continue;
2726 
2727  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
2728  struct vertexuse *vu;
2729  struct vertex *v1, *v2;
2730 
2731  NMG_CK_EDGEUSE(eu);
2732 
2733  v1 = eu->vu_p->v_p;
2734  NMG_CK_VERTEX(v1);
2735  v2 = eu->eumate_p->vu_p->v_p;
2736  NMG_CK_VERTEX(v2);
2737 
2738  if (!nmg_find_v_in_face(v1, fu2) ||
2739  !nmg_find_v_in_face(v2, fu2))
2740  continue;
2741 
2742  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2743  bu_log("Making EU %p an intersect line for face cutting\n", (void *)eu);
2744  }
2745 
2746  for (BU_LIST_FOR(vu, vertexuse, &v1->vu_hd)) {
2747  struct faceuse *fu;
2748 
2749  fu = nmg_find_fu_of_vu(vu);
2750 
2751  if (fu == fu2)
2752  nmg_enlist_one_vu(is, vu, 0.0);
2753  }
2754 
2755  for (BU_LIST_FOR(vu, vertexuse, &v2->vu_hd)) {
2756  struct faceuse *fu;
2757 
2758  fu = nmg_find_fu_of_vu(vu);
2759 
2760  if (fu == fu2)
2761  nmg_enlist_one_vu(is, vu, 1.0);
2762  }
2763 
2764  /* Now do face cutting */
2765  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2766  bu_log("Calling face cutter for fu2 %p\n", (void *)fu2);
2767  }
2768  nmg_fcut_face_2d(is->l2, is->mag2, fu2, fu1, &is->tol);
2769 
2770  bu_ptbl_reset(is->l1);
2771  bu_ptbl_reset(is->l2);
2772  }
2773  }
2774 
2775  for (BU_LIST_FOR(lu, loopuse, &fu2->lu_hd)) {
2776  struct edgeuse *eu;
2777 
2778  NMG_CK_LOOPUSE(lu);
2779 
2780  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
2781  continue;
2782 
2783  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
2784  struct vertexuse *vu;
2785  struct vertex *v1, *v2;
2786 
2787  NMG_CK_EDGEUSE(eu);
2788 
2789  v1 = eu->vu_p->v_p;
2790  NMG_CK_VERTEX(v1);
2791  v2 = eu->eumate_p->vu_p->v_p;
2792  NMG_CK_VERTEX(v2);
2793 
2794  if (!nmg_find_v_in_face(v1, fu1) ||
2795  !nmg_find_v_in_face(v2, fu1))
2796  continue;
2797 
2798  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2799  bu_log("Making EU %p an intersect line for face cutting\n", (void *)eu);
2800  }
2801 
2802  for (BU_LIST_FOR(vu, vertexuse, &v1->vu_hd)) {
2803  struct faceuse *fu;
2804 
2805  fu = nmg_find_fu_of_vu(vu);
2806 
2807  if (fu == fu1)
2808  nmg_enlist_one_vu(is, vu, 0.0);
2809  }
2810 
2811  for (BU_LIST_FOR(vu, vertexuse, &v2->vu_hd)) {
2812  struct faceuse *fu;
2813 
2814  fu = nmg_find_fu_of_vu(vu);
2815 
2816  if (fu == fu1)
2817  nmg_enlist_one_vu(is, vu, 1.0);
2818  }
2819 
2820  /* Now do face cutting */
2821  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
2822  bu_log("Calling face cutter for fu1 %p\n", (void *)fu1);
2823  }
2824  nmg_fcut_face_2d(is->l1, is->mag1, fu1, fu2, &is->tol);
2825 
2826  bu_ptbl_reset(is->l1);
2827  bu_ptbl_reset(is->l2);
2828  }
2829  }
2830  if (mag1)
2831  bu_free((char *)mag1, "mag1");
2832  if (mag2)
2833  bu_free((char *)mag2, "mag2");
2834 
2835  bu_ptbl_free(is->l1);
2836  bu_ptbl_free(is->l2);
2837 }
2838 
2839 
2840 /**
2841  * A parallel to nmg_isect_edge2p_edge2p().
2842  *
2843  * Intersect the line with eu1, from fu1.
2844  * The resulting vu's are added to "list", not is->l1 or is->l2.
2845  * fu2 is the "other" face on this intersect line, and is used only
2846  * when searching for existing vertex structs suitable for re-use.
2847  *
2848  * Returns -
2849  * Number of times edge is broken (0 or 1).
2850  */
2851 int
2852 nmg_isect_line2_edge2p(struct nmg_inter_struct *is, struct bu_ptbl *list, struct edgeuse *eu1, struct faceuse *fu1, struct faceuse *fu2)
2853 {
2854  point_t eu1_start; /* 2D */
2855  point_t eu1_end; /* 2D */
2856  vect_t eu1_dir; /* 2D */
2857  fastf_t dist[2];
2858  int status;
2859  point_t hit_pt; /* 3D */
2860  struct vertexuse *vu1a, *vu1b;
2861  int ret = 0;
2862 
2863  NMG_CK_INTER_STRUCT(is);
2864  BU_CK_PTBL(list);
2865  NMG_CK_EDGEUSE(eu1);
2866  NMG_CK_FACEUSE(fu1);
2867  NMG_CK_FACEUSE(fu2);
2868 
2869  /*
2870  * Important note: don't use eu1->eumate_p->vu_p here,
2871  * because that vu is in the opposite orientation faceuse.
2872  * Putting those vu's on the intersection line makes for big trouble.
2873  */
2874  vu1a = eu1->vu_p;
2875  vu1b = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
2876  NMG_CK_VERTEXUSE(vu1a);
2877  NMG_CK_VERTEXUSE(vu1b);
2878 
2879  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2880  bu_log("nmg_isect_line2_edge2p(eu1=%p, fu1=%p)\n\tvu1a=%p vu1b=%p\n\tv2a=%p v2b=%p\n",
2881  (void *)eu1, (void *)fu1,
2882  (void *)vu1a, (void *)vu1b,
2883  (void *)vu1a->v_p, (void *)vu1b->v_p);
2884  }
2885 
2886  /*
2887  * The 3D line in is->pt and is->dir is prepared by the caller.
2888  */
2889  nmg_get_2d_vertex(eu1_start, vu1a->v_p, is, &fu1->l.magic);
2890  nmg_get_2d_vertex(eu1_end, vu1b->v_p, is, &fu1->l.magic);
2891  V2SUB2(eu1_dir, eu1_end, eu1_start);
2892 
2893  dist[0] = dist[1] = 0; /* for clean prints, below */
2894 
2895  /* Intersect the line with the edge, in 2D */
2896  status = bn_isect_line2_lseg2(dist, is->pt2d, is->dir2d,
2897  eu1_start, eu1_dir, &is->tol);
2898 
2899  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2900  bu_log("\tbn_isect_line2_lseg2()=%d, dist: %g, %g\n",
2901  status, dist[0], dist[1]);
2902  }
2903 
2904  if (status < 0) goto out; /* No geometric intersection */
2905 
2906  if (status == 0) {
2907  /*
2908  * The edge is collinear with the line.
2909  * List both vertexuse structures, and return.
2910  */
2911  if (RTG.NMG_debug & DEBUG_POLYSECT)
2912  bu_log("\t\tedge collinear with isect line. Listing vu1a, vu1b\n");
2913  nmg_enlist_vu(is, vu1a, 0, MAX_FASTF);
2914  nmg_enlist_vu(is, vu1b, 0, MAX_FASTF);
2915  ret = 0;
2916  goto out;
2917  }
2918 
2919  /* There is only one intersect point. Break the edge there. */
2920 
2921  VJOIN1(hit_pt, is->pt, dist[0], is->dir); /* 3D hit */
2922 
2923  /* Edges not collinear. Either list a vertex,
2924  * or break eu1.
2925  */
2926  if (status == 1 || ZERO(dist[1])) {
2927  if (RTG.NMG_debug & DEBUG_POLYSECT)
2928  bu_log("\t\tintersect point is vu1a\n");
2929  if (!bn_pt3_pt3_equal(hit_pt, vu1a->v_p->vg_p->coord, &(is->tol)))
2930  bu_bomb("vu1a does not match calculated point\n");
2931  nmg_enlist_vu(is, vu1a, 0, MAX_FASTF);
2932  ret = 0;
2933  } else if (status == 2 || ZERO(dist[1] - 1.0)) {
2934  if (RTG.NMG_debug & DEBUG_POLYSECT)
2935  bu_log("\t\tintersect point is vu1b\n");
2936  if (!bn_pt3_pt3_equal(hit_pt, vu1b->v_p->vg_p->coord, &(is->tol)))
2937  bu_bomb("vu1b does not match calculated point\n");
2938  nmg_enlist_vu(is, vu1b, 0, MAX_FASTF);
2939  ret = 0;
2940  } else {
2941  /* Intersection is in the middle of eu1, split edge */
2942  fastf_t distance;
2943  struct vertexuse *vu1_final;
2944  struct vertex *new_v;
2945  if (RTG.NMG_debug & DEBUG_POLYSECT) {
2946  int code;
2947  bu_log("\t2D: pt2d=(%g, %g), dir2d=(%g, %g)\n",
2948  is->pt2d[X], is->pt2d[Y],
2949  is->dir2d[X], is->dir2d[Y]);
2950  bu_log("\t2D: eu1_start=(%g, %g), eu1_dir=(%g, %g)\n",
2951  eu1_start[X], eu1_start[Y],
2952  eu1_dir[X], eu1_dir[Y]);
2953  VPRINT("\t3D: is->pt ", is->pt);
2954  VPRINT("\t3D: is->dir", is->dir);
2955  bu_log("\t2d: Breaking eu1 at isect.\n");
2956  VPRINT(" vu1a", vu1a->v_p->vg_p->coord);
2957  VPRINT("hit_pt", hit_pt);
2958  VPRINT(" vu1b", vu1b->v_p->vg_p->coord);
2959  /* XXX Perform a (not-so) quick check */
2960  code = bn_isect_pt_lseg(&distance, vu1a->v_p->vg_p->coord,
2961  vu1b->v_p->vg_p->coord,
2962  hit_pt, &(is->tol));
2963  bu_log("\tbn_isect_pt_lseg() dist=%g, ret=%d\n", distance, code);
2964  if (code < 0) bu_bomb("3D point not on 3D lseg\n");
2965 
2966  /* Ensure that the 3D hit_pt is between the end pts */
2967  if (!bn_between(vu1a->v_p->vg_p->coord[X], hit_pt[X], vu1b->v_p->vg_p->coord[X], &(is->tol)) ||
2968  !bn_between(vu1a->v_p->vg_p->coord[Y], hit_pt[Y], vu1b->v_p->vg_p->coord[Y], &(is->tol)) ||
2969  !bn_between(vu1a->v_p->vg_p->coord[Z], hit_pt[Z], vu1b->v_p->vg_p->coord[Z], &(is->tol))) {
2970  VPRINT("vu1a", vu1a->v_p->vg_p->coord);
2971  VPRINT("hitp", hit_pt);
2972  VPRINT("vu1b", vu1b->v_p->vg_p->coord);
2973  bu_bomb("nmg_isect_line2_edge2p() hit point not between edge verts!\n");
2974  }
2975 
2976  }
2977 
2978  /* if we can't find the appropriate vertex
2979  * by a geometry search, build a new vertex.
2980  * Otherwise, re-use the existing one.
2981  * Can't just search other face, might miss relevant vert.
2982  */
2983  new_v = nmg_find_pt_in_model(fu2->s_p->r_p->m_p, hit_pt, &(is->tol));
2984  vu1_final = nmg_ebreaker(new_v, eu1, &is->tol)->vu_p;
2985  ret = 1;
2986  if (!new_v) {
2987  nmg_vertex_gv(vu1_final->v_p, hit_pt); /* 3d geom */
2988  if (RTG.NMG_debug & DEBUG_POLYSECT)
2989  bu_log("\t\tmaking new vertex vu=%p v=%p\n",
2990  (void *)vu1_final, (void *)vu1_final->v_p);
2991  } else {
2992  if (RTG.NMG_debug & DEBUG_POLYSECT)
2993  bu_log("\t\tre-using vertex v=%p vu=%p\n", (void *)new_v, (void *)vu1_final);
2994  }
2995  nmg_enlist_vu(is, vu1_final, 0, MAX_FASTF);
2996 
2998  }
2999 
3000 out:
3001  if (RTG.NMG_debug & DEBUG_POLYSECT)
3002  bu_log("nmg_isect_line2_edge2p(eu1=%p, fu1=%p) END ret=%d\n",
3003  (void *)eu1, (void *)fu1, ret);
3004  return ret;
3005 }
3006 
3007 
3008 /**
3009  * If this lone vertex lies along the intersect line, then add it to
3010  * the lists.
3011  *
3012  * Called from nmg_isect_line2_face2p().
3013  */
3014 void
3015 nmg_isect_line2_vertex2(struct nmg_inter_struct *is, struct vertexuse *vu1, struct faceuse *fu1)
3016 {
3017 
3018  NMG_CK_INTER_STRUCT(is);
3019  NMG_CK_VERTEXUSE(vu1);
3020  NMG_CK_FACEUSE(fu1);
3021 
3022  if (RTG.NMG_debug & DEBUG_POLYSECT)
3023  bu_log("nmg_isect_line2_vertex2(vu=%p)\n", (void *)vu1);
3024 
3025  /* Needs to be a 3D comparison */
3026  if (bn_distsq_line3_pt3(is->pt, is->dir, vu1->v_p->vg_p->coord) > is->tol.dist_sq)
3027  return;
3028 
3029  if (RTG.NMG_debug & DEBUG_POLYSECT)
3030  bu_log("nmg_isect_line2_vertex2(vu=%p) line hits vertex v=%p\n",
3031  (void *)vu1, (void *)vu1->v_p);
3032 
3033  nmg_enlist_vu(is, vu1, 0, MAX_FASTF);
3034 }
3035 
3036 
3037 /**
3038  *
3039  * Given two pointer tables filled with edgeuses representing two different
3040  * edge geometry lines, see if there is a common vertex of intersection.
3041  * If so, enlist the intersection.
3042  *
3043  * Returns -
3044  * 1 intersection found
3045  * 0 no intersection
3046  */
3047 int
3048 nmg_isect_two_ptbls(struct nmg_inter_struct *is, const struct bu_ptbl *t1, const struct bu_ptbl *t2)
3049 {
3050  const struct edgeuse **eu1;
3051  const struct edgeuse **eu2;
3052  struct vertexuse *vu1a;
3053  struct vertexuse *vu1b;
3054 
3055  NMG_CK_INTER_STRUCT(is);
3056  BU_CK_PTBL(t1);
3057 
3058  for (eu1 = (const struct edgeuse **)BU_PTBL_LASTADDR(t1);
3059  eu1 >= (const struct edgeuse **)BU_PTBL_BASEADDR(t1); eu1--
3060  ) {
3061  struct vertex *v1a;
3062  struct vertex *v1b;
3063 
3064  vu1a = (*eu1)->vu_p;
3065  vu1b = BU_LIST_PNEXT_CIRC(edgeuse, (*eu1))->vu_p;
3066  NMG_CK_VERTEXUSE(vu1a);
3067  NMG_CK_VERTEXUSE(vu1b);
3068  v1a = vu1a->v_p;
3069  v1b = vu1b->v_p;
3070 
3071  for (eu2 = (const struct edgeuse **)BU_PTBL_LASTADDR(t2);
3072  eu2 >= (const struct edgeuse **)BU_PTBL_BASEADDR(t2); eu2--
3073  ) {
3074  register struct vertexuse *vu2a;
3075  register struct vertexuse *vu2b;
3076 
3077  vu2a = (*eu2)->vu_p;
3078  vu2b = BU_LIST_PNEXT_CIRC(edgeuse, (*eu2))->vu_p;
3079  NMG_CK_VERTEXUSE(vu2a);
3080  NMG_CK_VERTEXUSE(vu2b);
3081 
3082  if (v1a == vu2a->v_p) {
3083  vu1b = vu2a;
3084  goto enlist;
3085  }
3086  if (v1a == vu2b->v_p) {
3087  vu1b = vu2b;
3088  goto enlist;
3089  }
3090  if (v1b == vu2a->v_p) {
3091  vu1a = vu1b;
3092  vu1b = vu2a;
3093  goto enlist;
3094  }
3095  if (v1b == vu2b->v_p) {
3096  vu1a = vu1b;
3097  vu1b = vu2b;
3098  goto enlist;
3099  }
3100  }
3101  }
3102  return 0;
3103 enlist:
3104  /* Two vu's are now vu1a, vu1b */
3105  if (nmg_find_s_of_vu(vu1a) == nmg_find_s_of_vu(vu1b)) {
3106  vu1b = 0;
3107  }
3108 
3109  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3110  bu_log("nmg_isect_two_ptbls() intersection! vu=%p, vu_dual=%p\n",
3111  (void *)vu1a, (void *)vu1b);
3112  }
3113  nmg_enlist_vu(is, vu1a, vu1b, MAX_FASTF);
3114  return 1;
3115 }
3116 
3117 
3118 /**
3119  * Do a geometric search to find an edge_g_lseg on the given line.
3120  * If the fuser did its job, there should be only one.
3121  */
3122 struct edge_g_lseg *
3123 nmg_find_eg_on_line(const uint32_t *magic_p, const fastf_t *pt, const fastf_t *dir, const struct bn_tol *tol)
3124 {
3125  struct bu_ptbl eutab;
3126  struct edgeuse **eup;
3127  struct edge_g_lseg *ret = (struct edge_g_lseg *)NULL;
3128  vect_t dir1, dir2;
3129 
3130  BN_CK_TOL(tol);
3131 
3132  nmg_edgeuse_on_line_tabulate(&eutab, magic_p, pt, dir, tol);
3133 
3134  for (eup = (struct edgeuse **)BU_PTBL_LASTADDR(&eutab);
3135  eup >= (struct edgeuse **)BU_PTBL_BASEADDR(&eutab); eup--
3136  ) {
3137  if (!ret) {
3138  /* No edge_g_lseg found yet, use this one. */
3139  ret = (*eup)->g.lseg_p;
3140  continue;
3141  }
3142  if ((*eup)->g.lseg_p == ret) continue; /* OK */
3143 
3144  /* Found 2 different edge_g_lseg, pick best one */
3145  VMOVE(dir1, ret->e_dir);
3146  VUNITIZE(dir1);
3147  VMOVE(dir2, (*eup)->g.lseg_p->e_dir);
3148  VUNITIZE(dir2);
3149  if (fabs(VDOT(dir1, dir)) > fabs(VDOT(dir2, dir))) {
3150  /* ret is better, do nothing */
3151  } else {
3152  /* *eup is better, take it instead */
3153  ret = (*eup)->g.lseg_p;
3154  }
3155  bu_log("nmg_find_eg_on_line() 2 different eg's, taking better one.\n");
3156  }
3157  (void)bu_ptbl_free(&eutab);
3158  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3159  bu_log("rt_find_eg_on_line(%p) ret=%p\n", (void *)magic_p, (void *)ret);
3160  }
3161  return ret;
3162 }
3163 
3164 
3165 /**
3166  * Kill all 0-length edgeuses that start and end on this vertex.
3167  *
3168  * Returns -
3169  * 0 If none were found
3170  * count Number of 0-length edgeuses killed (not counting mates)
3171  */
3172 int
3173 nmg_k0eu(struct vertex *v)
3174 {
3175  struct edgeuse *eu;
3176  struct vertexuse *vu;
3177  int count = 0;
3178 
3179  NMG_CK_VERTEX(v);
3180 top:
3181  for (BU_LIST_FOR(vu, vertexuse, &v->vu_hd)) {
3182  NMG_CK_VERTEXUSE(vu);
3183  if (*vu->up.magic_p != NMG_EDGEUSE_MAGIC) continue;
3184  eu = vu->up.eu_p;
3185  NMG_CK_EDGEUSE(eu);
3186  if (eu->eumate_p->vu_p->v_p != v) continue;
3187  bu_log("nmg_k0eu(v=%p) killing 0-len eu=%p, mate=%p\n",
3188  (void *)v, (void *)eu, (void *)eu->eumate_p);
3189  nmg_pr_eu_briefly(eu, 0);
3190  nmg_pr_eu_briefly(eu->eumate_p, 0);
3191  if (nmg_keu(eu)) {
3192  bu_log("nmg_k0eu() WARNING: parent now has no edgeuses\n");
3193  /* XXX Now what? */
3194  }
3195  count++;
3196  goto top; /* vu_hd list is altered by nmg_keu() */
3197  }
3198  return count;
3199 }
3200 
3201 
3202 /**
3203  * Attempt to join two vertices which both claim to be the intersection
3204  * of two lines. If they are close enough, repair the damage.
3205  *
3206  * Returns -
3207  * hit_v If repair succeeds. vertex 'v' is now invalid.
3208  * NULL If repair fails.
3209  * If 'bomb' is non-zero, bu_bomb() is called.
3210  */
3211 struct vertex *
3212 nmg_repair_v_near_v(struct vertex *hit_v, struct vertex *v, const struct edge_g_lseg *eg1, const struct edge_g_lseg *eg2, int bomb, const struct bn_tol *tol)
3213 
3214 
3215 /* edge_g_lseg of hit_v */
3216 /* edge_g_lseg of v */
3217 
3218 
3219 {
3220  NMG_CK_VERTEX(hit_v);
3221  NMG_CK_VERTEX(v);
3222  if (eg1) NMG_CK_EDGE_G_LSEG(eg1); /* eg1 may be NULL */
3223  NMG_CK_EDGE_G_LSEG(eg2);
3224  BN_CK_TOL(tol);
3225 
3226  bu_log("nmg_repair_v_near_v(hit_v=%p, v=%p)\n", (void *)hit_v, (void *)v);
3227 
3228  VPRINT("v ", v->vg_p->coord);
3229  VPRINT("hit", hit_v->vg_p->coord);
3230  bu_log("dist v-hit=%g, equal=%d\n",
3231  bn_dist_pt3_pt3(v->vg_p->coord, hit_v->vg_p->coord),
3232  bn_pt3_pt3_equal(v->vg_p->coord, hit_v->vg_p->coord, tol)
3233  );
3234  if (eg1) {
3235  if (bn_2line3_colinear(eg1->e_pt, eg1->e_dir, eg2->e_pt, eg2->e_dir, 1e5, tol))
3236  bu_bomb("ERROR: nmg_repair_v_near_v() eg1 and eg2 are collinear!\n");
3237  bu_log("eg1: line/ vu dist=%g, hit dist=%g\n",
3238  bn_dist_line3_pt3(eg1->e_pt, eg1->e_dir, v->vg_p->coord),
3239  bn_dist_line3_pt3(eg1->e_pt, eg1->e_dir, hit_v->vg_p->coord));
3240  bu_log("eg2: line/ vu dist=%g, hit dist=%g\n",
3241  bn_dist_line3_pt3(eg2->e_pt, eg2->e_dir, v->vg_p->coord),
3242  bn_dist_line3_pt3(eg2->e_pt, eg2->e_dir, hit_v->vg_p->coord));
3243  nmg_pr_eg(&eg1->l.magic, 0);
3244  nmg_pr_eg(&eg2->l.magic, 0);
3245  }
3246 
3247  if (bn_dist_pt3_pt3(v->vg_p->coord,
3248  hit_v->vg_p->coord) < 10 * tol->dist) {
3249  struct edgeuse *eu0;
3250  bu_log("NOTICE: The intersection of two lines has resulted in 2 different intersect points\n");
3251  bu_log(" Since the two points are 'close', they are being fused.\n");
3252 
3253  /* See if there is an edge between them */
3254  eu0 = nmg_findeu(hit_v, v, (struct shell *)NULL,
3255  (struct edgeuse *)NULL, 0);
3256  if (eu0) {
3257  bu_log("DANGER: a 0-length edge is being created eu0=%p\n", (void *)eu0);
3258  }
3259 
3260  nmg_jv(hit_v, v);
3261  (void)nmg_k0eu(hit_v);
3262  goto out;
3263  }
3264  /* Separation is too great */
3265  if (bomb)
3266  bu_bomb("nmg_repair_v_near_v() separation is too great to repair.\n");
3267  hit_v = (struct vertex *)NULL;
3268 out:
3269  bu_log("nmg_repair_v_near_v(v=%p) ret=%p\n", (void *)v, (void *)hit_v);
3270  return hit_v;
3271 }
3272 
3273 
3274 /**
3275  * Search all edgeuses referring to this vu's vertex.
3276  * If the vertex is used by edges on both eg1 and eg2, then it's a "hit"
3277  * between the two edge geometries.
3278  * If a new hit happens at a different vertex from a previous hit,
3279  * that is a fatal error.
3280  *
3281  * This routine exists only as a support routine for nmg_common_v_2eg().
3282  *
3283  * XXX This is a lame name.
3284  */
3285 struct vertex *
3286 nmg_search_v_eg(const struct edgeuse *eu, int second, const struct edge_g_lseg *eg1, const struct edge_g_lseg *eg2, register struct vertex *hit_v, const struct bn_tol *tol)
3287 
3288 /* 2nd vu on eu, not 1st */
3289 
3290 
3291 /* often will be NULL */
3292 
3293 {
3294  struct vertex *v;
3295  register struct vertexuse *vu1;
3296  register struct edgeuse *seen1 = (struct edgeuse *)NULL;
3297  register struct edgeuse *seen2 = (struct edgeuse *)NULL;
3298 
3299  NMG_CK_EDGEUSE(eu);
3300  NMG_CK_EDGE_G_LSEG(eg1);
3301  NMG_CK_EDGE_G_LSEG(eg2);
3302  BN_CK_TOL(tol);
3303 
3304  if (second) {
3305  v = BU_LIST_PNEXT_CIRC(edgeuse, eu)->vu_p->v_p;
3306  if (v != eu->eumate_p->vu_p->v_p)
3307  bu_bomb("nmg_search_v_eg() next vu not mate's vu?\n");
3308  } else {
3309  v = eu->vu_p->v_p;
3310  }
3311  NMG_CK_VERTEX(v);
3312 
3313  if (eu->g.lseg_p != eg1) bu_bomb("nmg_search_v_eg() eu not on eg1\n");
3314 
3315  /* vu lies on eg1 by topology. Check this assertion. */
3316  if (bn_distsq_line3_pt3(eg1->e_pt, eg1->e_dir, v->vg_p->coord) > tol->dist_sq) {
3317  VPRINT("v", v->vg_p->coord);
3318  nmg_pr_eu(eu, (char *)NULL);
3319  nmg_pr_eg(&eg1->l.magic, 0);
3320  bu_bomb("nmg_search_v_eg() eu vertex not on eg line\n");
3321  }
3322 
3323  /* This loop accounts for 30% of the runtime of a boolean! */
3324  for (BU_LIST_FOR(vu1, vertexuse, &v->vu_hd)) {
3325  register struct edgeuse *eu1;
3326 
3327  if (*vu1->up.magic_p != NMG_EDGEUSE_MAGIC) continue;
3328  eu1 = vu1->up.eu_p;
3329  if (eu1->g.lseg_p == eg1) seen1 = eu1;
3330  if (eu1->g.lseg_p == eg2) seen2 = eu1;
3331  if (!seen1 || !seen2) continue;
3332 
3333  /* Both edge_g's have been seen at 'v', this is a hit. */
3334  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3335  bu_log(" seen1=%p, seen2=%p, hit_v=%p, v=%p\n",
3336  (void *)seen1, (void *)seen2, (void *)hit_v, (void *)v);
3337  }
3338  if (!hit_v) {
3339  hit_v = v;
3340  break;
3341  }
3342 
3343  /* Is it a different vertex than hit_v? */
3344  if (hit_v == v) break;
3345 
3346  /* Different vertices, this "can't happen" */
3347  bu_log("ERROR seen1=%p, seen2=%p, hit_v=%p != v=%p\n",
3348  (void *)seen1, (void *)seen2, (void *)hit_v, (void *)v);
3349  if (nmg_repair_v_near_v(hit_v, v, eg1, eg2, 0, tol))
3350  break;
3351 
3352  bu_bomb("nmg_search_v_eg() two different vertices for intersect point?\n");
3353  }
3354  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3355  bu_log("nmg_search_v_eg(eu=%p, %d, eg1=%p, eg2=%p) ret=%p\n",
3356  (void *)eu, second, (void *)eg1, (void *)eg2, (void *)hit_v);
3357  }
3358  return hit_v;
3359 }
3360 
3361 
3362 /**
3363  * Perform a topology search for a common vertex between two edge geometry
3364  * lines.
3365  */
3366 struct vertex *
3367 nmg_common_v_2eg(struct edge_g_lseg *eg1, struct edge_g_lseg *eg2, const struct bn_tol *tol)
3368 {
3369  struct edgeuse *eu1;
3370  struct vertex *hit_v = (struct vertex *)NULL;
3371  struct bu_list *midway; /* &eu->l2, midway into edgeuse */
3372 
3373  NMG_CK_EDGE_G_LSEG(eg1);
3374  NMG_CK_EDGE_G_LSEG(eg2);
3375  BN_CK_TOL(tol);
3376 
3377  if (eg1 == eg2)
3378  bu_bomb("nmg_common_v_2eg() eg1 and eg2 are collinear\n");
3379 
3380  /* Scan all edgeuses in the model that use eg1 */
3381  for (BU_LIST_FOR(midway, bu_list, &eg1->eu_hd2)) {
3382  NMG_CKMAG(midway, NMG_EDGEUSE2_MAGIC, "edgeuse2 [l2]");
3383  eu1 = BU_LIST_MAIN_PTR(edgeuse, midway, l2);
3384  NMG_CK_EDGEUSE(eu1);
3385  if (eu1->g.lseg_p != eg1) bu_bomb("nmg_common_v_2eg() eu disavows eg\n");
3386  /* Both verts of eu1 lie on line eg1 */
3387  hit_v = nmg_search_v_eg(eu1, 0, eg1, eg2, hit_v, tol);
3388  hit_v = nmg_search_v_eg(eu1, 1, eg1, eg2, hit_v, tol);
3389  }
3390  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3391  bu_log("nmg_common_v_2eg(eg1=%p, eg2=%p) hit_v=%p\n",
3392  (void *)eg1, (void *)eg2, (void *)hit_v);
3393  }
3394  return hit_v;
3395 }
3396 
3397 
3398 #define VDIST(a, b) sqrt((a[X]-b[X])*(a[X]-b[X]) + (a[Y]-b[Y])*(a[Y]-b[Y]) + (a[Z]-b[Z])*(a[Z]-b[Z]))
3399 #define VDIST_SQ(a, b) ((a[X]-b[X])*(a[X]-b[X]) + (a[Y]-b[Y])*(a[Y]-b[Y]) + (a[Z]-b[Z])*(a[Z]-b[Z]))
3400 
3401 int
3402 nmg_is_vertex_on_inter(struct vertex *v, struct faceuse *fu1, struct faceuse *fu2, struct nmg_inter_struct *is)
3403 {
3404  struct vertex_g *vg;
3405  plane_t pl1, pl2;
3406  int code;
3407  fastf_t dist;
3408 
3409  NMG_CK_VERTEX(v);
3410  NMG_CK_FACEUSE(fu1);
3411  NMG_CK_FACEUSE(fu2);
3412  NMG_CK_INTER_STRUCT(is);
3413 
3414  if (nmg_find_v_in_face(v, fu1) && nmg_find_v_in_face(v, fu2))
3415  return 1;
3416 
3417  NMG_GET_FU_PLANE(pl1, fu1);
3418  NMG_GET_FU_PLANE(pl2, fu2);
3419 
3420  vg = v->vg_p;
3421  NMG_CK_VERTEX_G(vg);
3422 
3423  /* check if vertex is in plane of fu's */
3424  dist = fabs(DIST_PT_PLANE(vg->coord, pl1));
3425  if (dist > is->tol.dist)
3426  return 0;
3427  dist = fabs(DIST_PT_PLANE(vg->coord, pl2));
3428  if (dist > is->tol.dist)
3429  return 0;
3430 
3431  /* check if it is on intersection line */
3432  if (bn_distsq_line3_pt3(is->pt, is->dir, vg->coord) > is->tol.dist_sq)
3433  return 0;
3434 
3435  /* check if it is within fu's */
3436  code = nmg_class_pt_fu_except(vg->coord, fu1, (struct loopuse *)NULL,
3437  (void (*)(struct edgeuse *, point_t, const char *))NULL,
3438  (void (*)(struct vertexuse *, point_t, const char *))NULL,
3439  (const char *)NULL, 0, 0, &is->tol);
3440  if (code != NMG_CLASS_AinB)
3441  return 0;
3442 
3443  code = nmg_class_pt_fu_except(vg->coord, fu2, (struct loopuse *)NULL,
3444  (void (*)(struct edgeuse *, point_t, const char *))NULL,
3445  (void (*)(struct vertexuse *, point_t, const char *))NULL,
3446  (const char *)NULL, 0, 0, &is->tol);
3447  if (code != NMG_CLASS_AinB)
3448  return 0;
3449 
3450  return 1;
3451 }
3452 
3453 
3454 void
3455 nmg_isect_eu_verts(struct edgeuse *eu, struct vertex_g *vg1, struct vertex_g *vg2, struct bu_ptbl *verts, struct bu_ptbl *inters, const struct bn_tol *tol)
3456 {
3457  int i;
3458  struct vertex *v1, *v2;
3459 
3460  NMG_CK_EDGEUSE(eu);
3461  NMG_CK_VERTEX_G(vg1);
3462  NMG_CK_VERTEX_G(vg2);
3463  BU_CK_PTBL(verts);
3464  BU_CK_PTBL(inters);
3465  BN_CK_TOL(tol);
3466 
3467  v1 = eu->vu_p->v_p;
3468  v2 = eu->eumate_p->vu_p->v_p;
3469 
3470  for (i=0; i<BU_PTBL_END(verts); i++) {
3471  struct vertex *v;
3472  fastf_t dist;
3473  point_t pca;
3474  int code;
3475 
3476  v = (struct vertex *)BU_PTBL_GET(verts, i);
3477  if (v == v1 || v == v2) {
3478  bu_ptbl_ins_unique(inters, (long *)v);
3479  continue;
3480  }
3481 
3482  code = bn_dist_pt3_lseg3(&dist, pca, vg1->coord,
3483  vg2->coord, v->vg_p->coord, tol);
3484 
3485  if (code)
3486  continue;
3487 
3488  bu_ptbl_ins_unique(inters, (long *)v);
3489  }
3490 
3491  return;
3492 }
3493 
3494 
3495 void
3496 nmg_isect_eu_eu(struct edgeuse *eu1, struct vertex_g *vg1a, struct vertex_g *vg1b, fastf_t *dir1, struct edgeuse *eu2, struct bu_ptbl *verts, struct bu_ptbl *inters, const struct bn_tol *tol)
3497 {
3498  struct model *m;
3499  struct vertex_g *vg2a, *vg2b;
3500  vect_t dir2, diff;
3501  fastf_t dist[2];
3502  int code;
3503  point_t hit_pt;
3504  vect_t e1_min_pt, e1_max_pt, e2_min_pt, e2_max_pt;
3505 
3506  if (RTG.NMG_debug & DEBUG_POLYSECT)
3507  bu_log("nmg_isect_eu_eu(eu1=%p, eu2=%p)\n", (void *)eu1, (void *)eu2);
3508 
3509  m = nmg_find_model(&eu1->l.magic);
3510 
3511  vg2a = eu2->vu_p->v_p->vg_p;
3512  vg2b = eu2->eumate_p->vu_p->v_p->vg_p;
3513 
3514  /* compute bounding box for eu1 */
3515  VMOVE(e1_min_pt, vg1a->coord);
3516  VMIN(e1_min_pt, vg1b->coord);
3517  VMOVE(e1_max_pt, vg1a->coord);
3518  VMAX(e1_max_pt, vg1b->coord);
3519 
3520  /* compute bounding box for eu2 */
3521  VMOVE(e2_min_pt, vg2a->coord);
3522  VMIN(e2_min_pt, vg2b->coord);
3523  VMOVE(e2_max_pt, vg2a->coord);
3524  VMAX(e2_max_pt, vg2b->coord);
3525 
3526  if (V3RPP_DISJOINT_TOL(e1_min_pt, e1_max_pt, e2_min_pt, e2_max_pt, tol->dist)) {
3527  return;
3528  }
3529 
3530  VSUB2(dir2, vg2b->coord, vg2a->coord);
3531 
3532  code = bn_isect_lseg3_lseg3(dist, vg1a->coord, dir1, vg2a->coord, dir2, tol);
3533 
3534  if (code < 0) {
3535  if (RTG.NMG_debug & DEBUG_POLYSECT)
3536  bu_log("\tnmg_isect_eu_eu: No intersection\n");
3537  return;
3538  }
3539 
3540  if (code == 1) {
3541  point_t hit_pt1, hit_pt2;
3542  struct vertex *v=(struct vertex *)NULL;
3543  struct edgeuse *new_eu;
3544 
3545  /* normal intersection (one point) */
3546 
3547  if (eu1->vu_p->v_p == eu2->vu_p->v_p ||
3548  eu1->vu_p->v_p == eu2->eumate_p->vu_p->v_p ||
3549  eu1->eumate_p->vu_p->v_p == eu2->vu_p->v_p ||
3550  eu1->eumate_p->vu_p->v_p == eu2->eumate_p->vu_p->v_p)
3551  return;
3552 
3553  if (ZERO(dist[0]) || ZERO(dist[0] - 1.0))
3554  return;
3555 
3556  if (ZERO(dist[1])) {
3557  bu_ptbl_ins_unique(inters, (long *)eu2->vu_p->v_p);
3558  return;
3559  }
3560  if (ZERO(dist[1] - 1.0)) {
3561  bu_ptbl_ins_unique(inters, (long *)eu2->eumate_p->vu_p->v_p);
3562  return;
3563  }
3564 
3565  VJOIN1(hit_pt1, vg1a->coord, dist[0], dir1);
3566  VJOIN1(hit_pt2, vg2a->coord, dist[1], dir2);
3567 
3568  VBLEND2(hit_pt, 0.5, hit_pt1, 0.5, hit_pt2);
3569 
3570  v = nmg_find_pt_in_model(m, hit_pt, tol);
3571 
3572  if (RTG.NMG_debug & DEBUG_POLYSECT) {
3573  bu_log("nmg_isect_eu_eu: intersection at (%g %g %g)\n", V3ARGS(hit_pt));
3574  bu_log("splitting eu %p at v=%p\n", (void *)eu2, (void *)v);
3575  }
3576  new_eu = nmg_esplit(v, eu2, 1);
3577  if (!v) {
3578  v = new_eu->vu_p->v_p;
3579  nmg_vertex_gv(v, hit_pt);
3580  if (RTG.NMG_debug & DEBUG_POLYSECT)
3581  bu_log("\tcreated new vertex %p\n", (void *)v);
3582  }
3583  bu_ptbl_ins_unique(inters, (long *)v);
3584  bu_ptbl_ins_unique(verts, (long *)v);
3585  return;
3586  }
3587 
3588  /* code == 0, could be two intersection points
3589  * But there should be no vertex creation here
3590  */
3591 
3592  VSUB2(diff, vg2a->coord, vg1a->coord);
3593  if (VDOT(diff, dir1) > SMALL_FASTF) {
3594  VSUB2(diff, vg1b->coord, vg2a->coord);
3595  if (VDOT(diff, dir1) > SMALL_FASTF)
3596  bu_ptbl_ins_unique(inters, (long *)eu2->vu_p->v_p);
3597  }
3598 
3599  VSUB2(diff, vg2b->coord, vg1a->coord);
3600  if (VDOT(diff, dir1) > SMALL_FASTF) {
3601  VSUB2(diff, vg1b->coord, vg2b->coord);
3602  if (VDOT(diff, dir1) > SMALL_FASTF)
3603  bu_ptbl_ins_unique(inters, (long *)eu2->eumate_p->vu_p->v_p);
3604  }
3605 }
3606 
3607 
3608 void
3609 nmg_isect_eu_fu(struct nmg_inter_struct *is, struct bu_ptbl *verts, struct edgeuse *eu, struct faceuse *fu)
3610 {
3611  struct model *m;
3612  struct vertex_g *vg1, *vg2;
3613  struct loopuse *lu;
3614  plane_t pl;
3615  fastf_t dist;
3616  fastf_t eu_len;
3617  fastf_t one_over_len;
3618  point_t hit_pt;
3619  vect_t edir;
3620  vect_t dir;
3621  struct bu_ptbl inters;
3622  fastf_t *inter_dist;
3623  int i;
3624 
3625  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3626  bu_log("nmg_isect_eu_fu: eu=%p, fu=%p START\n", (void *)eu, (void *)fu);
3627 
3628  NMG_CK_INTER_STRUCT(is);
3629  NMG_CK_FACEUSE(fu);
3630  NMG_CK_EDGEUSE(eu);
3631  BU_CK_PTBL(verts);
3632 
3633  if (UNLIKELY(nmg_find_fu_of_eu(eu) == fu)) {
3634  bu_log("nmg_isect_eu_fu() called with eu (%p) from its own fu (%p)\n",
3635  (void *)eu, (void *)fu);
3636  bu_bomb("nmg_isect_eu_fu() called with eu from its own fu");
3637  }
3638 
3639  m = nmg_find_model(&fu->l.magic);
3640  NMG_CK_MODEL(m);
3641  if (UNLIKELY(nmg_find_model(&eu->l.magic) != m)) {
3642  bu_log("nmg_isect_eu_fu() called with EU (%p) from model (%p)\n",
3643  (void *)eu, (void *)nmg_find_model(&eu->l.magic));
3644  bu_log("\tand FU (%p) from model (%p)\n", (void *)fu, (void *)m);
3645  bu_bomb("nmg_isect_eu_fu() called with EU and FU from different models");
3646  }
3647 
3648  vg1 = eu->vu_p->v_p->vg_p;
3649  NMG_CK_VERTEX_G(vg1);
3650  vg2 = eu->eumate_p->vu_p->v_p->vg_p;
3651  NMG_CK_VERTEX_G(vg2);
3652 
3653  VSUB2(dir, vg2->coord, vg1->coord);
3654  VMOVE(edir, dir);
3655  eu_len = MAGNITUDE(dir);
3656  if (eu_len < is->tol.dist) {
3657  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3658  bu_log("\tnmg_isec_eu_fu: 0 length edge\n");
3659  return;
3660  }
3661 
3662  one_over_len = 1.0/eu_len;
3663  VSCALE(dir, dir, one_over_len);
3664 
3665  NMG_GET_FU_PLANE(pl, fu);
3666  /* check if edge line intersects plane of fu */
3667  if (bn_isect_line3_plane(&dist, vg1->coord, dir, pl, &is->tol) < 1) {
3668  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3669  bu_log("\tnmg_isec_eu_fu: no intersection\n");
3670  return;
3671  }
3672 
3673  VJOIN1(hit_pt, vg1->coord, dist, dir);
3674 
3675  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3676  bu_log("\tintersection point at (%g %g %g)\n", V3ARGS(hit_pt));
3677 
3678  /* create a list of intersection vertices */
3679  bu_ptbl_init(&inters, 64, " &inters");
3680 
3681  /* add vertices from fu to list */
3682  nmg_isect_eu_verts(eu, vg1, vg2, verts, &inters, &is->tol);
3683 
3684  /* break FU EU's that intersect our eu, and put vertices on list */
3685  for (BU_LIST_FOR(lu, loopuse, &fu->lu_hd)) {
3686  struct edgeuse *eu_fu;
3687 
3688  NMG_CK_LOOPUSE(lu);
3689 
3690  /* vertices of FU are handled above */
3691  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
3692  continue;
3693 
3694  for (BU_LIST_FOR(eu_fu, edgeuse, &lu->down_hd)) {
3695  NMG_CK_EDGEUSE(eu_fu);
3696 
3697  nmg_isect_eu_eu(eu, vg1, vg2, edir, eu_fu, verts, &inters, &is->tol);
3698 
3699  }
3700  }
3701 
3702  /* Now break eu at every vertex in the "inters" list */
3703 
3704  if (BU_PTBL_END(&inters) == 0) {
3705  struct vertex *v=(struct vertex *)NULL;
3706  int nmg_class;
3707  fastf_t dist_to_plane;
3708 
3709  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3710  bu_log("\tNo intersection points found\n");
3711 
3712  /* check if EU endpoints are within tolerance of FU
3713  * If so, the endpoint is the intersection and nothing to do
3714  */
3715  dist_to_plane = fabs(DIST_PT_PLANE(vg1->coord, pl));
3716  if (dist_to_plane < is->tol.dist) {
3717  /* check if hit point is within fu */
3718  nmg_class = nmg_class_pt_fu_except(vg1->coord, fu, (struct loopuse *)NULL,
3719  0, 0, (char *)NULL, 0, 0, &is->tol);
3720  if (nmg_class != NMG_CLASS_AinB)
3721  goto out;
3722 
3723  v = eu->vu_p->v_p;
3724  if (v && !nmg_find_v_in_face(v, fu)) {
3725  struct vertexuse *new_vu;
3726 
3727  new_vu = nmg_make_dualvu(v, fu, &is->tol);
3728  bu_ptbl_ins_unique(verts, (long *)new_vu->v_p);
3729  }
3730  goto out;
3731  }
3732 
3733  dist_to_plane = fabs(DIST_PT_PLANE(vg2->coord, pl));
3734  if (dist_to_plane < is->tol.dist) {
3735  /* check if hit point is within fu */
3736  nmg_class = nmg_class_pt_fu_except(vg2->coord, fu, (struct loopuse *)NULL,
3737  0, 0, (char *)NULL, 0, 0, &is->tol);
3738  if (nmg_class != NMG_CLASS_AinB)
3739  goto out;
3740 
3741  v = eu->eumate_p->vu_p->v_p;
3742  if (v && !nmg_find_v_in_face(v, fu)) {
3743  struct vertexuse *new_vu;
3744 
3745  new_vu = nmg_make_dualvu(v, fu, &is->tol);
3746  bu_ptbl_ins_unique(verts, (long *)new_vu->v_p);
3747  }
3748  goto out;
3749  }
3750 
3751  /* no vertices found, we might need to create a self loop */
3752 
3753  /* make sure intersection is within limits of eu */
3754  if (dist < (-is->tol.dist) || dist > eu_len+is->tol.dist) {
3755  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3756  bu_log("\tnmg_isec_eu_fu: intersection beyond ends of EU\n");
3757  goto out;
3758  }
3759 
3760  /* check if hit_point is within tolerance of an end of eu */
3761  if (bn_pt3_pt3_equal(hit_pt, vg1->coord, &is->tol)) {
3762  v = eu->vu_p->v_p;
3763  VMOVE(hit_pt, vg1->coord);
3764  } else if (bn_pt3_pt3_equal(hit_pt, vg2->coord, &is->tol)) {
3765  v = eu->eumate_p->vu_p->v_p;
3766  VMOVE(hit_pt, vg2->coord);
3767  }
3768 
3769  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
3770  bu_log("\tHit point is not within tolerance of eu endpoints\n");
3771  bu_log("\t\thit_pt=(%g %g %g), eu=(%g %g %g)<->(%g %g %g)\n",
3772  V3ARGS(hit_pt), V3ARGS(vg1->coord), V3ARGS(vg2->coord));
3773  }
3774 
3775  /* check if hit point is within fu */
3776  nmg_class = nmg_class_pt_fu_except(hit_pt, fu, (struct loopuse *)NULL,
3777  0, 0, (char *)NULL, 0, 0, &is->tol);
3778 
3779  if (nmg_class == NMG_CLASS_AinB) {
3780  struct edgeuse *new_eu;
3781 
3782  /* may need to split eu */
3783  if (!v)
3784  v = nmg_find_pt_in_model(m, hit_pt, &is->tol);
3785  if (v != eu->vu_p->v_p && v != eu->eumate_p->vu_p->v_p) {
3786  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3787  bu_log("\tsplitting eu (%p) at hit_pt (v=%p)\n", (void *)eu, (void *)v);
3788 
3789  new_eu = nmg_esplit(v, eu, 1);
3790  if (!v) {
3791  v = new_eu->vu_p->v_p;
3792  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3793  bu_log("\tnew vertex at hit point is %p\n", (void *)v);
3794  nmg_vertex_gv(v, hit_pt);
3795  }
3796  }
3797 
3798  if (v && !nmg_find_v_in_face(v, fu)) {
3799  struct vertexuse *new_vu;
3800 
3801  new_vu = nmg_make_dualvu(v, fu, &is->tol);
3802  bu_ptbl_ins_unique(verts, (long *)new_vu->v_p);
3803  }
3804 
3805  }
3806  goto out;
3807  }
3808 
3809  if (BU_PTBL_END(&inters) == 1) {
3810  struct vertex *v;
3811 
3812  /* only one vertex, just split */
3813  v = (struct vertex *)BU_PTBL_GET(&inters, 0);
3814  NMG_CK_VERTEX(v);
3815 
3816  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3817  bu_log("Only one intersect vertex (%p), just split all EU's at (%p)\n", (void *)v, (void *)eu);
3818 
3819  if (v == eu->vu_p->v_p || v == eu->eumate_p->vu_p->v_p)
3820  goto out;
3821 
3822  (void)nmg_break_all_es_on_v(&fu->l.magic, v, &is->tol);
3823  (void)nmg_break_all_es_on_v(&eu->l.magic, v, &is->tol);
3824 
3825  goto out;
3826  }
3827 
3828  /* must do them in order from furthest to nearest */
3829  inter_dist = (fastf_t *)bu_calloc(BU_PTBL_END(&inters), sizeof(fastf_t),
3830  "nmg_isect_eu_fu: inter_dist");
3831 
3832  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3833  bu_log("%ld intersect vertices along eu (%p)\n", BU_PTBL_END(&inters), (void *)eu);
3834 
3835  for (i=0; i<BU_PTBL_END(&inters); i++) {
3836  struct vertex *v;
3837  struct vertex_g *vg;
3838  vect_t diff;
3839 
3840  v = (struct vertex *)BU_PTBL_GET(&inters, i);
3841  NMG_CK_VERTEX(v);
3842  vg = v->vg_p;
3843  NMG_CK_VERTEX_G(vg);
3844 
3845  VSUB2(diff, vg->coord, vg1->coord);
3846  if (UNLIKELY(VDOT(diff, dir) < -SMALL_FASTF)) {
3847  bu_bomb("nmg_isect_eu_fu: intersection point not on eu\n");
3848  }
3849  inter_dist[i] = MAGSQ(diff);
3850  }
3851 
3852  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT)) {
3853  bu_log("Intersect vertices along eu %p:\n", (void *)eu);
3854  for (i = 0; i < BU_PTBL_END(&inters); i++)
3855  bu_log("%d %p %g\n", i+1, (void *)BU_PTBL_GET(&inters, i), inter_dist[i]);
3856  }
3857 
3858  while (1) {
3859  struct vertex *v;
3860  fastf_t max_dist;
3861  int index_at_max;
3862 
3863  max_dist = (-1.0);
3864  index_at_max = (-1);
3865 
3866  for (i=0; i<BU_PTBL_END(&inters); i++) {
3867  if (inter_dist[i] > max_dist) {
3868  max_dist = inter_dist[i];
3869  index_at_max = i;
3870  }
3871  }
3872 
3873  if (index_at_max < 0)
3874  break;
3875 
3876  v = (struct vertex *)BU_PTBL_GET(&inters, index_at_max);
3877  NMG_CK_VERTEX(v);
3878 
3879  if (v != eu->vu_p->v_p && v != eu->eumate_p->vu_p->v_p) {
3880  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3881  bu_log("Breaking edges at vertex #%d, dist=%g, v=%p\n", i+1, inter_dist[i], (void *)v);
3882  (void)nmg_break_all_es_on_v(&fu->l.magic, v, &is->tol);
3883  (void)nmg_break_all_es_on_v(&eu->l.magic, v, &is->tol);
3884  }
3885 
3886  inter_dist[index_at_max] = (-10.0);
3887  }
3888 
3889  bu_free((char *)inter_dist, "nmg_isect_eu_fu: inter_dist");
3890 
3891 out:
3892  bu_ptbl_free(&inters);
3893 
3894  if (UNLIKELY(RTG.NMG_debug & DEBUG_POLYSECT))
3895  bu_log("nmg_isect_eu_fu: eu=%p, fu=%p END\n", (void *)eu, (void *)fu);
3896 
3897 }
3898 
3899 
3900 void
3901 nmg_isect_fu_jra(struct nmg_inter_struct *is, struct faceuse *fu1, struct faceuse *fu2, struct bu_ptbl *eu1_list, struct bu_ptbl *eu2_list)
3902 {
3903  struct model *m;
3904  struct bu_ptbl verts1, verts2, eus;
3905  ssize_t i;
3906 
3907  if (RTG.NMG_debug & DEBUG_POLYSECT)
3908  bu_log("nmg_isect_fu_jra(fu1=%p, fu2=%p) START\n", (void *)fu1, (void *)fu2);
3909 
3910  NMG_CK_INTER_STRUCT(is);
3911  NMG_CK_FACEUSE(fu1);
3912  NMG_CK_FACEUSE(fu2);
3913  BU_CK_PTBL(eu1_list);
3914  BU_CK_PTBL(eu2_list);
3915 
3916  m = nmg_find_model(&fu1->l.magic);
3917  NMG_CK_MODEL(m);
3918 
3919  /* Intersect fu1 edgeuses */
3920  nmg_vertex_tabulate(&verts2, &fu2->l.magic);
3921  nmg_edgeuse_tabulate(&eus, &fu1->l.magic);
3922  for (i = 0; i < BU_PTBL_END(&eus); i++) {
3923  struct edgeuse *eu;
3924 
3925  eu = (struct edgeuse *)BU_PTBL_GET(&eus, i);
3926  NMG_CK_EDGEUSE(eu);
3927  if (eu->g.magic_p && *eu->g.magic_p == NMG_EDGE_G_CNURB_MAGIC) {
3928  continue;
3929  }
3930  nmg_isect_eu_fu(is, &verts2, eu, fu2);
3931  }
3932  bu_ptbl_free(&verts2);
3933  bu_ptbl_free(&eus);
3934 
3935  /* Intersect fu2 edgeuses */
3936  nmg_vertex_tabulate(&verts1, &fu1->l.magic);
3937  nmg_edgeuse_tabulate(&eus, &fu2->l.magic);
3938  for (i = 0; i < BU_PTBL_END(&eus); i++) {
3939  struct edgeuse *eu;
3940 
3941  eu = (struct edgeuse *)BU_PTBL_GET(&eus, i);
3942  NMG_CK_EDGEUSE(eu);
3943  if (eu->g.magic_p && *eu->g.magic_p == NMG_EDGE_G_CNURB_MAGIC) {
3944  continue;
3945  }
3946  nmg_isect_eu_fu(is, &verts1, eu, fu1);
3947  }
3948  bu_ptbl_free(&verts1);
3949  bu_ptbl_free(&eus);
3950 
3951  /* check for existing vertices along intersection */
3952  /* XXXX this is the second time this tabulate is being done,
3953  * but for now it's safer this way
3954  */
3955  nmg_vertex_tabulate(&verts1, &fu1->l.magic);
3956  nmg_vertex_tabulate(&verts2, &fu2->l.magic);
3957 
3958  /* merge the two lists */
3959  for (i=0; i<BU_PTBL_END(&verts2); i++) {
3960  struct vertex *v;
3961 
3962  v = (struct vertex *)BU_PTBL_GET(&verts2, i);
3963  NMG_CK_VERTEX(v);
3964 
3965  bu_ptbl_ins_unique(&verts1, (long *)v);
3966  }
3967  bu_ptbl_free(&verts2);
3968 
3969  for (i=0; i< BU_PTBL_END(&verts1); i++) {
3970  struct vertex *v;
3971  struct vertexuse *vu;
3972  fastf_t dist;
3973 
3974  v = (struct vertex *)BU_PTBL_GET(&verts1, i);
3975  NMG_CK_VERTEX(v);
3976 
3977  if (!nmg_is_vertex_on_inter(v, fu1, fu2, is))
3978  continue;
3979 
3980  /* calculate distance along intersect ray */
3981  dist = VDIST(is->pt, v->vg_p->coord);
3982 
3983  /* this vertex is on the intersection line
3984  * add all uses from fu1 and fu2 to intersection list
3985  */
3986  for (BU_LIST_FOR(vu, vertexuse, &v->vu_hd)) {
3987  struct faceuse *fu_tmp;
3988 
3989  fu_tmp = nmg_find_fu_of_vu(vu);
3990  if (fu_tmp == fu1 || fu_tmp == fu2) {
3991  if (RTG.NMG_debug & DEBUG_POLYSECT)
3992  bu_log("\tenlisting vu %p (%p) from fu (%p)\n",
3993  (void *)vu, (void *)v, (void *)fu_tmp);
3994  nmg_enlist_one_vu(is, vu, dist);
3995  }
3996  }
3997  }
3998 
3999  bu_ptbl_free(&verts1);
4000 
4001  if (RTG.NMG_debug & DEBUG_POLYSECT)
4002  bu_log("nmg_isect_fu_jra(fu1=%p, fu2=%p) END\n", (void *)fu1, (void *)fu2);
4003 }
4004 
4005 
4006 /**
4007  * HEART
4008  *
4009  * For each distinct edge_g_lseg LINE on the face (composed of potentially many
4010  * edgeuses and many different edges), intersect with the edge_g_lseg LINE
4011  * which represents the face/face intersection line.
4012  *
4013  * Note that the geometric intersection of the two faces is
4014  * stored in is->pt and is->dir.
4015  * If is->on_eg is set, it is the callers' responsibility to make sure
4016  * it is not much different than the original geometric one.
4017  *
4018  * Go to great pains to ensure that two non-collinear lines intersect
4019  * at either 0 or 1 points, and no more.
4020  *
4021  * Called from -
4022  * nmg_isect_edge2p_face2p()
4023  * nmg_isect_two_face3p()
4024  */
4025 void
4026 nmg_isect_line2_face2pNEW(struct nmg_inter_struct *is, struct faceuse *fu1, struct faceuse *fu2, struct bu_ptbl *eu1_list, struct bu_ptbl *eu2_list)
4027 {
4028  struct bu_ptbl eg_list;
4029  struct edge_g_lseg **eg1;
4030  struct edgeuse *eu1;
4031  struct edgeuse *eu2;
4032  fastf_t dist[2];
4033  int code = 0;
4034  point_t eg_pt2d; /* 2D */
4035  vect_t eg_dir2d; /* 2D */
4036  struct loopuse *lu1;
4037  point_t hit3d;
4038  point_t hit2d; /* 2D */
4039  struct edgeuse *new_eu;
4040  int eu1_index;
4041  int eu2_index;
4042  int nmg_class;
4043  fastf_t distance;
4044 
4045  NMG_CK_INTER_STRUCT(is);
4046  NMG_CK_FACEUSE(fu1);
4047  NMG_CK_FACEUSE(fu2);
4048  BU_CK_PTBL(eu1_list);
4049  BU_CK_PTBL(eu2_list);
4050 
4051  if (RTG.NMG_debug & DEBUG_POLYSECT)
4052  bu_log("nmg_isect_line2_face2pNEW(, fu1=%p, fu2=%p) on_eg=%p\n",
4053  (void *)fu1, (void *)fu2, (void *)is->on_eg);
4054 
4055  /* Project the intersect line into 2D. Build matrix first. */
4056  nmg_isect2d_prep(is, &fu1->l.magic);
4057  /* XXX Need subroutine for this!! */
4058  MAT4X3PNT(is->pt2d, is->proj, is->pt);
4059  MAT4X3VEC(is->dir2d, is->proj, is->dir);
4060 
4061 re_tabulate:
4062  /* Build list of all edge_g_lseg's in fu1 */
4063  /* XXX This could be more cheaply done by cooking down eu1_list */
4064  nmg_edge_g_tabulate(&eg_list, &fu1->l.magic);
4065 
4066  /* Process each distinct line in the face fu1 */
4067  for (eg1 = (struct edge_g_lseg **)BU_PTBL_LASTADDR(&eg_list);
4068  eg1 >= (struct edge_g_lseg **)BU_PTBL_BASEADDR(&eg_list); eg1--
4069  ) {
4070  struct vertex *hit_v = (struct vertex *)NULL;
4071 
4072  NMG_CK_EDGE_G_LSEG(*eg1);
4073 
4074  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4075  bu_log("\tChecking eg=%p\n", (void *)*eg1);
4076  }
4077 
4078  if (*eg1 == is->on_eg) {
4079  point_t pca;
4080  struct edgeuse *eu_end;
4081  struct vertex_g *vg;
4082  plane_t pl1, pl2;
4083 
4084  colinear:
4085  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4086  bu_log("\tThis edge_geom generated the line. Enlisting.\n");
4087  }
4088 
4089  NMG_GET_FU_PLANE(pl1, is->fu1);
4090  NMG_GET_FU_PLANE(pl2, is->fu2);
4091  for (eu1_index=0; eu1_index < BU_PTBL_END(eu1_list); eu1_index++) {
4092  eu1 = (struct edgeuse *)BU_PTBL_GET(eu1_list, eu1_index);
4093  NMG_CK_EDGEUSE(eu1);
4094  if (eu1->g.lseg_p != is->on_eg) continue;
4095  /* eu1 is from fu1 */
4096 
4097  vg = eu1->vu_p->v_p->vg_p;
4098  NMG_CK_VERTEX_G(vg);
4099  (void)bn_dist_pt3_line3(&distance, pca, is->pt, is->dir, vg->coord, &(is->tol));
4100  if (distance <= is->tol.dist) {
4101  /* vertex is on intersection line */
4102 
4103  if (DIST_PT_PLANE(vg->coord, pl2) <= is->tol.dist) {
4104  /* and in plane of fu2 */
4105  if (nmg_class_pt_fu_except(vg->coord, fu2, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol)) != NMG_CLASS_AoutB) {
4106  /* and within fu2 */
4107  distance = VDIST(is->pt, vg->coord);
4108  nmg_enlist_vu(is, eu1->vu_p, 0, distance);
4109  }
4110  }
4111  }
4112  eu_end = BU_LIST_PNEXT_CIRC(edgeuse, &eu1->l);
4113  vg = eu_end->vu_p->v_p->vg_p;
4114  NMG_CK_VERTEX_G(vg);
4115  code = bn_dist_pt3_line3(&distance, pca, is->pt, is->dir, eu_end->vu_p->v_p->vg_p->coord, &(is->tol));
4116  if (distance <= is->tol.dist) {
4117  /* vertex is on intersection line */
4118 
4119  if (DIST_PT_PLANE(vg->coord, pl2) <= is->tol.dist) {
4120  /* and in plane of fu2 */
4121  if (nmg_class_pt_fu_except(vg->coord, fu2, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol)) != NMG_CLASS_AoutB) {
4122  /* and within fu2 */
4123  distance = VDIST(is->pt, vg->coord);
4124  nmg_enlist_vu(is, eu_end->vu_p, 0, distance);
4125  }
4126  }
4127  }
4128  }
4129 
4130  for (eu2_index=0; eu2_index < BU_PTBL_END(eu2_list); eu2_index++) {
4131  eu2 = (struct edgeuse *)BU_PTBL_GET(eu2_list, eu2_index);
4132  NMG_CK_EDGEUSE(eu2);
4133 
4134  if (eu2->g.lseg_p != is->on_eg) continue;
4135  /* eu2 is from fu2 */
4136 
4137  vg = eu2->vu_p->v_p->vg_p;
4138  NMG_CK_VERTEX_G(vg);
4139  code = bn_dist_pt3_line3(&distance, pca, is->pt, is->dir, vg->coord, &(is->tol));
4140  if (distance <= is->tol.dist) {
4141  /* vertex is on intersection line */
4142 
4143  if (DIST_PT_PLANE(vg->coord, pl1) <= is->tol.dist) {
4144  /* and in plane of fu1 */
4145  if (nmg_class_pt_fu_except(vg->coord, fu1, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol)) != NMG_CLASS_AoutB) {
4146  /* and within fu1 */
4147  distance = VDIST(is->pt, vg->coord);
4148  nmg_enlist_vu(is, eu2->vu_p, 0, distance);
4149  }
4150  }
4151  }
4152  eu_end = BU_LIST_PNEXT_CIRC(edgeuse, &eu2->l);
4153  vg = eu_end->vu_p->v_p->vg_p;
4154  NMG_CK_VERTEX_G(vg);
4155  code = bn_dist_pt3_line3(&distance, pca, is->pt, is->dir, eu_end->vu_p->v_p->vg_p->coord, &(is->tol));
4156  if (distance <= is->tol.dist) {
4157  /* vertex is on intersection line */
4158 
4159  if (DIST_PT_PLANE(vg->coord, pl1) <= is->tol.dist) {
4160  /* and in plane of fu1 */
4161  if (nmg_class_pt_fu_except(vg->coord, fu1, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol)) != NMG_CLASS_AoutB) {
4162  /* and within fu1 */
4163  distance = VDIST(is->pt, vg->coord);
4164  nmg_enlist_vu(is, eu_end->vu_p, 0, distance);
4165  }
4166  }
4167  }
4168  }
4169  continue;
4170  }
4171 
4172  /*
4173  * eg1 is now known to be NOT collinear with on_eg.
4174  * From here on, only 0 or 1 points of intersection are possible.
4175  */
4176 
4177  /* The 3D line in is->pt and is->dir is prepared by our caller. */
4178  MAT4X3PNT(eg_pt2d, is->proj, (*eg1)->e_pt);
4179  MAT4X3VEC(eg_dir2d, is->proj, (*eg1)->e_dir);
4180 
4181  /* Calculate 2D geometric intersection, but don't look at answer yet */
4182  dist[0] = dist[1] = -INFINITY;
4183  code = bn_isect_line2_line2(dist, is->pt2d, is->dir2d,
4184  eg_pt2d, eg_dir2d, &(is->tol));
4185 
4186  /* Do this check before topology search */
4187  if (code == 0) {
4188  /* Geometry says lines are collinear. Egads! This can't be! */
4189  if (is->on_eg) {
4190  bu_log("nmg_isect_line2_face2pNEW() edge_g not shared, geometry says lines are collinear.\n");
4191  goto fixup;
4192  }
4193  /* on_eg wasn't set, use it and continue on */
4194  if (RTG.NMG_debug & DEBUG_POLYSECT)
4195  bu_log("NOTICE: setting on_eg to eg1 and continuing with collinear case.\n");
4196  is->on_eg = (*eg1);
4197  goto colinear;
4198  }
4199 
4200  /* Double check */
4201  if (is->on_eg && bn_2line3_colinear(
4202  (*eg1)->e_pt, (*eg1)->e_dir,
4203  is->on_eg->e_pt, is->on_eg->e_dir, 1e5, &(is->tol))) {
4204  fixup:
4205  nmg_pr_eg(&(*eg1)->l.magic, 0);
4206  nmg_pr_eg(&is->on_eg->l.magic, 0);
4207  bu_log("nmg_isect_line2_face2pNEW() eg1 collinear to on_eg?\n");
4208 
4209  /* fuse eg1 with on_eg, handle as collinear */
4210  bu_log("fusing eg1 with on_eg, handling as collinear\n");
4211  nmg_jeg(is->on_eg, *eg1);
4212  goto colinear;
4213  }
4214 
4215  /* If on_eg was specified, do a search for topology intersection */
4216  if (is->on_eg && !hit_v) {
4217  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4218  bu_log("non-colinear. Searching for topological intersection between on_eg and eg1\n");
4219  }
4220  /* See if any vu along eg1 is used by edge from on_eg */
4221  hit_v = nmg_common_v_2eg(*eg1, is->on_eg, &(is->tol));
4222  /*
4223  * Note that while eg1 contains an eu from fu1,
4224  * the intersection vertex just found may occur
4225  * well outside *either* faceuse, but lie on some
4226  * other face that shares this face geometry.
4227  */
4228  if (hit_v) {
4229  fastf_t dist1, dist2;
4230  plane_t n1, n2;
4231 
4232  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4233  static int num=0;
4234  char buf[128];
4235  FILE *fp;
4236  sprintf(buf, "Itopo%d.plot3", num++);
4237  if ((fp=fopen(buf, "wb")) != NULL) {
4238  pl_color(fp, 255, 0, 0);
4239  pdv_3ray(fp, is->on_eg->e_pt, is->on_eg->e_dir, 1.0);
4240  pdv_3cont(fp, hit_v->vg_p->coord);
4241  pl_color(fp, 255, 255, 0);
4242  pdv_3ray(fp, (*eg1)->e_pt, (*eg1)->e_dir, 1.0);
4243  pdv_3cont(fp, hit_v->vg_p->coord);
4244  fclose(fp);
4245  bu_log("overlay %s red is on_eg, yellow is eg1\n", buf);
4246  } else perror(buf);
4247  bu_log("\tTopology intersection. hit_v=%p (%g, %g, %g)\n",
4248  (void *)hit_v,
4249  V3ARGS(hit_v->vg_p->coord));
4250  }
4251 
4252  /*
4253  * If the hit point is outside BOTH faces
4254  * bounding RPP, then it can be ignored.
4255  * Otherwise, it is needed for generating
4256  * accurate state transitions in the
4257  * face cutter.
4258  */
4259  if (!V3PT_IN_RPP(hit_v->vg_p->coord, fu1->f_p->min_pt, fu1->f_p->max_pt) &&
4260  !V3PT_IN_RPP(hit_v->vg_p->coord, fu2->f_p->min_pt, fu2->f_p->max_pt)
4261  ) {
4262  /* Lines intersect outside bounds of both faces. */
4263  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4264  VPRINT("\t\tisect pt outside fu1 & fu2 RPP:", hit_v->vg_p->coord);
4265  bu_log("\t\tfu1 RPP: (%g %g %g) <-> (%g %g %g)\n",
4266  V3ARGS(fu1->f_p->min_pt), V3ARGS(fu1->f_p->max_pt));
4267  bu_log("\t\tfu2 RPP: (%g %g %g) <-> (%g %g %g)\n",
4268  V3ARGS(fu2->f_p->min_pt), V3ARGS(fu2->f_p->max_pt));
4269  }
4270  continue;
4271  }
4272 
4273  /*
4274  * Geometry check.
4275  * Since on_eg is the line of intersection
4276  * between fu1 and fu2, and since eg1
4277  * lies in the plane of fu1,
4278  * and since hit_v lies on eg1
4279  * at the point where eg1 intersects on_eg,
4280  * it's pretty hard to see how hit_v
4281  * could fail to lie in the plane of fu2.
4282  *
4283  * However, in BigWedge r1 (m14.r) this
4284  * case occurs: hit_v is off f2 by -60mm!
4285  * These two fu's don't overlap at all.
4286  *
4287  * Rather than letting nmg_ck_v_in_2fus()
4288  * blow its mind over this, catch it here
4289  * and discard the point.
4290  * (Hopefully the RPP checks above will have
4291  * already discarded it).
4292  */
4293  NMG_GET_FU_PLANE(n1, fu1);
4294  NMG_GET_FU_PLANE(n2, fu2);
4295  dist1 = DIST_PT_PLANE(hit_v->vg_p->coord, n1);
4296  dist2 = DIST_PT_PLANE(hit_v->vg_p->coord, n2);
4297 
4298  if (!NEAR_ZERO(dist1, is->tol.dist) || !NEAR_ZERO(dist2, is->tol.dist)) {
4299  continue;
4300  }
4301  if ((nmg_class=nmg_class_pt_fu_except(hit_v->vg_p->coord, fu1, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol))) == NMG_CLASS_AoutB) {
4302  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4303  VPRINT("\t\tisect pt outside face fu1 (nmg_class_pt_fu_except):", hit_v->vg_p->coord);
4304  }
4305  continue;
4306  } else if ((nmg_class=nmg_class_pt_fu_except(hit_v->vg_p->coord, fu2, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol))) == NMG_CLASS_AoutB) {
4307  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4308  VPRINT("\t\tisect pt outside face fu2 (nmg_class_pt_fu_except):", hit_v->vg_p->coord);
4309  }
4310  continue;
4311  }
4312  /* Ignore further geometry checking, isect is good.
4313  * Head right to searching for vertexuses
4314  */
4315  goto force_isect;
4316  } else if (code < 0) {
4317  /* geometry says lines are parallel, so lack of intersection is expected */
4318  continue;
4319  }
4320  } else if (code < 0) {
4321  /* geometry says lines are parallel, but we have an intersection */
4322  bu_log("NOTICE: geom/topo mis-match, enlisting topo vu, hit_v=%p\n", (void *)hit_v);
4323  VPRINT("hit_v", hit_v->vg_p->coord);
4324  nmg_pr_eg(&(*eg1)->l.magic, 0);
4325  nmg_pr_eg(&is->on_eg->l.magic, 0);
4326  bu_log(" dist to eg1=%e, dist to on_eg=%e\n",
4327  bn_dist_line3_pt3((*eg1)->e_pt, (*eg1)->e_dir, hit_v->vg_p->coord),
4328  bn_dist_line3_pt3(is->on_eg->e_pt, is->on_eg->e_dir, hit_v->vg_p->coord));
4329  VPRINT("is->pt2d ", is->pt2d);
4330  VPRINT("is->dir2d", is->dir2d);
4331  VPRINT("eg_pt2d ", eg_pt2d);
4332  VPRINT("eg_dir2d ", eg_dir2d);
4333  bu_log(" 3d line isect, code=%d\n",
4334  bn_isect_line3_line3(&dist[0], &dist[1],
4335  is->pt, is->dir,
4336  (*eg1)->e_pt,
4337  (*eg1)->e_dir,
4338  &(is->tol)));
4339  goto force_isect;
4340  }
4341 
4342  /* Geometry says 2 lines intersect at a point */
4343  VJOIN1(hit3d, is->pt, dist[0], is->dir);
4344  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4345  VPRINT("\t2 lines intersect at", hit3d);
4346  }
4347 
4348  if (!V3PT_IN_RPP(hit3d, fu1->f_p->min_pt, fu1->f_p->max_pt)) {
4349  /* Lines intersect outside bounds of this face. */
4350  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4351  VPRINT("\t\tisect pt outside fu1 face RPP:", hit3d);
4352  bu_log("\t\tface RPP: (%g %g %g) <-> (%g %g %g)\n",
4353  V3ARGS(fu1->f_p->min_pt), V3ARGS(fu1->f_p->max_pt));
4354  }
4355  continue;
4356  }
4357 
4358  if ((nmg_class=nmg_class_pt_fu_except(hit3d, fu1, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol))) == NMG_CLASS_AoutB) {
4359  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4360  VPRINT("\t\tisect pt outside face fu1 (nmg_class_pt_fu_except):", hit3d);
4361  }
4362  continue;
4363  } else if ((nmg_class=nmg_class_pt_fu_except(hit3d, fu2, (struct loopuse *)NULL, NULL, NULL, (char *)NULL, 0, 0, &(is->tol))) == NMG_CLASS_AoutB) {
4364  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4365  VPRINT("\t\tisect pt outside face fu2 (nmg_class_pt_fu_except):", hit3d);
4366  }
4367  continue;
4368  } else if (RTG.NMG_debug & DEBUG_POLYSECT) {
4369  bu_log("\t\tnmg_class_pt_fu_except(fu1) returns %s\n", nmg_class_name(nmg_class));
4370  }
4371 
4372  V2JOIN1(hit2d, is->pt2d, dist[0], is->dir2d);
4373 
4374  /* Consistency check between geometry, and hit_v. */
4375  if (hit_v) {
4376  force_isect:
4377  /* Force things to be consistent, use geom from hit_v */
4378  VMOVE(hit3d, hit_v->vg_p->coord);
4379  nmg_get_2d_vertex(hit2d, hit_v, is, &fu1->l.magic);
4380  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4381  bu_log("hit_v=%p\n", (void *)hit_v);
4382  VPRINT("hit3d", hit3d);
4383  V2PRINT("hit2d", hit2d);
4384  }
4385  }
4386 
4387  eu_search:
4388  /* Search all eu's on eg1 for vu's to enlist.
4389  * There may be many, and the list may grow.
4390  */
4391  for (eu1_index=0; eu1_index < BU_PTBL_END(eu1_list); eu1_index++) {
4392  struct vertexuse *vu1a, *vu1b;
4393  struct vertexuse *vu1_midpt;
4394  fastf_t ldist;
4395  point_t eu1_pt2d; /* 2D */
4396  point_t eu1_end2d; /* 2D */
4397  double tmp_dist_sq;
4398 
4399  eu1 = (struct edgeuse *)BU_PTBL_GET(eu1_list, eu1_index);
4400 
4401  /* This EU may have been killed by nmg_repair_v_near_v() */
4402  if (eu1->l.magic != NMG_EDGEUSE_MAGIC)
4403  continue;
4404 
4405  NMG_CK_EDGEUSE(eu1);
4406 
4407  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4408  bu_log("\tChecking eu %p\n", (void *)eu1);
4409  }
4410 
4411  if (eu1->g.lseg_p != *eg1) {
4412  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4413  bu_log("\t\teg %p is not eg1=%p\n", (void *)eu1->g.lseg_p, (void *)*eg1);
4414  }
4415  continue;
4416  }
4417  vu1a = eu1->vu_p;
4418  vu1b = BU_LIST_PNEXT_CIRC(edgeuse, eu1)->vu_p;
4419 
4420  /* First, a topology check of both endpoints */
4421  if (vu1a->v_p == hit_v) {
4422  hit_a:
4423  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4424  bu_log("\tlisting intersect point at vu1a=%p\n", (void *)vu1a);
4425  }
4426  /* Note that the distance dist[0] may not actually
4427  * be the exact distance to vu1a, but it is the distance
4428  * to the actual intersection point. This is an
4429  * attempt to get the correct ordering of vertices
4430  * on the intersection list, since using the
4431  * actual distance can get them reversed when
4432  * a VU is chosen over the actual intersection
4433  * point.
4434  */
4435  nmg_enlist_vu(is, vu1a, 0, dist[0]);
4436  }
4437  if (vu1b->v_p == hit_v) {
4438  hit_b:
4439  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4440  bu_log("\tlisting intersect point at vu1b=%p\n", (void *)vu1b);
4441  }
4442  /* see above note about dist[0] */
4443  nmg_enlist_vu(is, vu1b, 0, dist[0]);
4444  }
4445  if (vu1a->v_p == hit_v || vu1b->v_p == hit_v) continue;
4446 
4447  /* Second, a geometry check on the edgeuse ENDPOINTS
4448  * -vs- the line segment. This is 3D, for consistency
4449  * with comparisons elsewhere.
4450  */
4451  tmp_dist_sq = bn_distsq_line3_pt3(is->pt, is->dir, vu1a->v_p->vg_p->coord);
4452  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4453  bu_log("\tvu1a is sqrt(%g) from the intersect line\n", tmp_dist_sq);
4454  }
4455  if (tmp_dist_sq <= is->tol.dist_sq) {
4456  if (!hit_v) {
4457  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4458  bu_log("\tnmg_isect_line2_face2pNEW: using nearby vu1a vertex %p from eu %p\n",
4459  (void *)vu1a->v_p, (void *)eu1);
4460  }
4461  hit_v = vu1a->v_p;
4462  goto hit_a;
4463  }
4464  if (hit_v == vu1a->v_p) goto hit_a;
4465 
4466  /* Fall through to bn_isect_pt2_lseg2() */
4467  }
4468  tmp_dist_sq = bn_distsq_line3_pt3(is->pt, is->dir, vu1b->v_p->vg_p->coord);
4469  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4470  bu_log("\tvu1b is sqrt(%g) from the intersect line\n", tmp_dist_sq);
4471  }
4472  if (tmp_dist_sq <= is->tol.dist_sq) {
4473  if (!hit_v) {
4474  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4475  bu_log("\tnmg_isect_line2_face2pNEW: using nearby vu1b vertex %p from eu %p\n",
4476  (void *)vu1b->v_p, (void *)eu1);
4477  }
4478  hit_v = vu1b->v_p;
4479  goto hit_b;
4480  }
4481  if (hit_v == vu1b->v_p) goto hit_b;
4482 
4483  /* Fall through to bn_isect_pt2_lseg2() */
4484  }
4485 
4486  /* Third, a geometry check of the HITPT -vs- the line segment */
4487  nmg_get_2d_vertex(eu1_pt2d, vu1a->v_p, is, &fu1->l.magic);
4488  nmg_get_2d_vertex(eu1_end2d, vu1b->v_p, is, &fu1->l.magic);
4489  ldist = 0;
4490  code = bn_isect_pt2_lseg2(&ldist, eu1_pt2d, eu1_end2d, hit2d, &(is->tol));
4491  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4492  bu_log("\tbn_isect_pt2_lseg2() returned %d, ldist=%g\n", code, ldist);
4493  }
4494  switch (code) {
4495  case -2:
4496  continue; /* outside lseg AB pts */
4497  default:
4498  case -1:
4499  continue; /* Point not on lseg */
4500  case 1:
4501  /* Point is at A (vu1a) by geometry */
4502  if (hit_v && hit_v != vu1a->v_p) {
4503  nmg_repair_v_near_v(hit_v, vu1a->v_p,
4504  is->on_eg, *eg1, 1, &(is->tol));
4505  bu_ptbl_free(&eg_list);
4506  goto re_tabulate;
4507  }
4508  hit_v = vu1a->v_p;
4509  if (RTG.NMG_debug & DEBUG_POLYSECT)
4510  bu_log("\thit_v = %p (vu1a)\n", (void *)hit_v);
4511  goto hit_a;
4512  case 2:
4513  /* Point is at B (vu1b) by geometry */
4514  if (hit_v && hit_v != vu1b->v_p) {
4515  nmg_repair_v_near_v(hit_v, vu1b->v_p,
4516  is->on_eg, *eg1, 1, &(is->tol));
4517  bu_ptbl_free(&eg_list);
4518  goto re_tabulate;
4519  }
4520  hit_v = vu1b->v_p;
4521  if (RTG.NMG_debug & DEBUG_POLYSECT)
4522  bu_log("\thit_v = %p (vu1b)\n", (void *)hit_v);
4523  goto hit_b;
4524  case 3:
4525  /* Point hits the line segment amidships! Split edge!
4526  * If we don't have a hit vertex yet,
4527  * search for one in whole model.
4528  */
4529  if (!hit_v) {
4530  hit_v = nmg_find_pt_in_model(fu1->s_p->r_p->m_p,
4531  hit3d, &(is->tol));
4532  if (hit_v == vu1a->v_p || hit_v == vu1b->v_p)
4533  bu_bomb("About to make 0-length edge!\n");
4534  }
4535  new_eu = nmg_ebreaker(hit_v, eu1, &is->tol);
4536  bu_ptbl_ins_unique(eu1_list, (long *)&new_eu->l.magic);
4537  /* "eu1" must now be considered invalid */
4538  vu1_midpt = new_eu->vu_p;
4539  if (!hit_v) {
4540  hit_v = vu1_midpt->v_p;
4541  nmg_vertex_gv(hit_v, hit3d);
4542  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4543  bu_log("\tmaking new vertex vu=%p hit_v=%p\n",
4544  (void *)vu1_midpt, (void *)hit_v);
4545  }
4546  /* Before we loose track of the fact
4547  * that this vertex lies on *both*
4548  * lines, break any edges in the
4549  * intersection line that cross it.
4550  */
4551  if (is->on_eg) {
4553  hit_v, &(is->tol));
4554  }
4555  } else {
4556  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4557  bu_log("\tre-using hit_v=%p, vu=%p\n", (void *)hit_v, (void *)vu1_midpt);
4558  }
4559  if (hit_v != vu1_midpt->v_p) bu_bomb("hit_v changed?\n");
4560  }
4561  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4562  bu_log("Faceuses after nmg_ebreaker() call\n");
4563  bu_log("fu1:\n");
4564  nmg_pr_fu_briefly(fu1, "\t");
4565  bu_log("fu2:\n");
4566  nmg_pr_fu_briefly(fu2, "\t");
4567  }
4568  nmg_enlist_vu(is, vu1_midpt, 0, dist[0]);
4569  /* Neither old nor new edgeuse need further handling */
4570  /* Because "eu1" is now invalid, restart loop. */
4571  goto eu_search;
4572  }
4573  }
4574  }
4575 
4576  /* Don't forget to do self loops with no edges */
4577  for (BU_LIST_FOR(lu1, loopuse, &fu1->lu_hd)) {
4578  struct vertexuse *vu1;
4579 
4580  if (BU_LIST_FIRST_MAGIC(&lu1->down_hd) != NMG_VERTEXUSE_MAGIC) continue;
4581 
4582  /* Intersect line with lone vertex vu1 */
4583  vu1 = BU_LIST_FIRST(vertexuse, &lu1->down_hd);
4584  NMG_CK_VERTEXUSE(vu1);
4585 
4586  /* Needs to be a 3D comparison */
4587  if (bn_distsq_line3_pt3(is->pt, is->dir,
4588  vu1->v_p->vg_p->coord) > is->tol.dist_sq)
4589  continue;
4590 
4591  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4592  bu_log("\tself-loop vu=%p lies on line of intersection\n", (void *)vu1);
4593  }
4594 
4595  /* Break all edgeuses in the model on line on_eg, at vu1 */
4596  if (is->on_eg) {
4597  nmg_break_eg_on_v(is->on_eg, vu1->v_p, &(is->tol));
4598  nmg_enlist_vu(is, vu1, 0, MAX_FASTF);
4599  }
4600 
4601  /* FIXME: This can probably be removed */
4602  /* recent OLD WAY: */
4603  /* To prevent an OT_BOOLPLACE from being overlooked,
4604  * break *both* sets of eu's
4605  */
4606 
4607  /* Break all edges from fu1 list */
4608  for (eu1_index=0; eu1_index < BU_PTBL_END(eu1_list); eu1_index++) {
4609  eu1 = (struct edgeuse *)BU_PTBL_GET(eu1_list, eu1_index);
4610 
4611  /* This EU may have been killed by nmg_repair_v_near_v() */
4612  if (eu1->l.magic != NMG_EDGEUSE_MAGIC)
4613  continue;
4614 
4615  NMG_CK_EDGEUSE(eu1);
4616  if (eu1->g.lseg_p != is->on_eg) continue;
4617  /* eu1 is from fu1 and on the intersection line */
4618  new_eu = nmg_break_eu_on_v(eu1, vu1->v_p, fu1, is);
4619  if (!new_eu) continue;
4620  bu_ptbl_ins_unique(eu1_list, (long *)&new_eu->l.magic);
4621  nmg_enlist_vu(is, new_eu->vu_p, 0, MAX_FASTF);
4622  }
4623 
4624  /* Break all edges from fu2 list */
4625  for (eu2_index=0; eu2_index < BU_PTBL_END(eu2_list); eu2_index++) {
4626  eu2 = (struct edgeuse *)BU_PTBL_GET(eu2_list, eu2_index);
4627 
4628  /* This EU may have been killed by nmg_repair_v_near_v() */
4629  if (eu2->l.magic != NMG_EDGEUSE_MAGIC)
4630  continue;
4631 
4632  NMG_CK_EDGEUSE(eu2);
4633  if (eu2->g.lseg_p != is->on_eg) continue;
4634  /* eu2 is from fu2 and on the intersection line */
4635  new_eu = nmg_break_eu_on_v(eu2, vu1->v_p, fu1, is);
4636  if (!new_eu) continue;
4637  bu_ptbl_ins_unique(eu2_list, (long *)&new_eu->l.magic);
4638  nmg_enlist_vu(is, new_eu->vu_p, 0, MAX_FASTF);
4639  }
4640 
4641  }
4642 
4643  bu_ptbl_free(&eg_list);
4644 }
4645 
4646 
4647 /* XXX move to nmg_info.c */
4648 int
4649 nmg_is_eu_on_line3(const struct edgeuse *eu, const fastf_t *UNUSED(pt), const fastf_t *dir, const struct bn_tol *tol)
4650 {
4651  struct edge_g_lseg *eg;
4652 
4653  NMG_CK_EDGEUSE(eu);
4654  BN_CK_TOL(tol);
4655 
4656  eg = eu->g.lseg_p;
4657  NMG_CK_EDGE_G_LSEG(eg);
4658 
4659  /* Ensure direction vectors are generally parallel */
4660  /* These are not unit vectors */
4661  /* tol->para and RT_DOT_TOL are too tight a tolerance. 0.1 is 5 degrees */
4662  if (fabs(VDOT(eg->e_dir, dir)) <
4663  0.9 * MAGNITUDE(eg->e_dir) * MAGNITUDE(dir)) return 0;
4664 
4665  /* XXX: this does not take pt+dir into account, a bug? */
4666 
4667  /* Ensure that vertices on edge are within tol of line */
4668  if (bn_distsq_line3_pt3(eg->e_pt, eg->e_dir,
4669  eu->vu_p->v_p->vg_p->coord) > tol->dist_sq) return 0;
4670  if (bn_distsq_line3_pt3(eg->e_pt, eg->e_dir,
4671  eu->eumate_p->vu_p->v_p->vg_p->coord) > tol->dist_sq) return 0;
4672 
4673 
4674  return 1;
4675 }
4676 
4677 
4678 /* XXX Move to nmg_info.c */
4679 /**
4680  * Perform a topology search to determine if two face geometries (specified
4681  * by their faceuses) share an edge geometry in common.
4682  * The edge_g is returned, even if there are no existing uses of it
4683  * in *either* fu1 or fu2. It still represents the intersection line
4684  * between the two face geometries, found topologically.
4685  *
4686  * If there are multiple edgeuses in common, ensure that they all refer
4687  * to the same edge_g geometry structure. The intersection of two planes
4688  * (non-coplanar) must be a single line.
4689  *
4690  * Calling this routine when the two faces share face geometry
4691  * is illegal.
4692  *
4693  * NULL is returned if no common edge geometry could be found.
4694  */
4695 struct edge_g_lseg *
4696 nmg_find_eg_between_2fg(const struct faceuse *ofu1, const struct faceuse *fu2, const struct bn_tol *tol)
4697 {
4698  const struct faceuse *fu1;
4699  const struct loopuse *lu1;
4700  const struct face_g_plane *fg1;
4701  const struct face_g_plane *fg2;
4702  const struct face *f1;
4703  struct edgeuse *ret = (struct edgeuse *)NULL;
4704  int coincident;
4705 
4706  NMG_CK_FACEUSE(ofu1);
4707  NMG_CK_FACEUSE(fu2);
4708  BN_CK_TOL(tol);
4709 
4710  fg1 = ofu1->f_p->g.plane_p;
4711  fg2 = fu2->f_p->g.plane_p;
4712  NMG_CK_FACE_G_PLANE(fg1);
4713  NMG_CK_FACE_G_PLANE(fg2);
4714 
4715  if (fg1 == fg2) bu_bomb("nmg_find_eg_between_2fg() face_g_plane shared, infinitely many results\n");
4716 
4717  if (RTG.NMG_debug & DEBUG_BASIC) {
4718  nmg_pr_fus_in_fg(&fg1->magic);
4719  nmg_pr_fus_in_fg(&fg2->magic);
4720  }
4721 
4722  /* For all faces using fg1 */
4723  for (BU_LIST_FOR(f1, face, &fg1->f_hd)) {
4724  NMG_CK_FACE(f1);
4725 
4726  /* Arbitrarily pick one of the two fu's using f1 */
4727  fu1 = f1->fu_p;
4728  NMG_CK_FACEUSE(fu1);
4729 
4730  for (BU_LIST_FOR(lu1, loopuse, &fu1->lu_hd)) {
4731  const struct edgeuse *eu1;
4732  NMG_CK_LOOPUSE(lu1);
4733  if (BU_LIST_FIRST_MAGIC(&lu1->down_hd) == NMG_VERTEXUSE_MAGIC)
4734  continue;
4735  if (RTG.NMG_debug & DEBUG_BASIC) {
4736  bu_log(" visiting lu1=%p, fu1=%p, fg1=%p\n",
4737  (void *)lu1, (void *)fu1, (void *)fg1);
4738  }
4739  restart:
4740  for (BU_LIST_FOR(eu1, edgeuse, &lu1->down_hd)) {
4741  struct edgeuse *eur;
4742 
4743  NMG_CK_EDGEUSE(eu1);
4744  /* Walk radially around the edge */
4745  for (
4746  eur = eu1->radial_p;
4747  eur != eu1->eumate_p;
4748  eur = eur->eumate_p->radial_p
4749  ) {
4750  const struct faceuse *tfu;
4751 
4752  if (*eur->up.magic_p != NMG_LOOPUSE_MAGIC) continue;
4753  if (*eur->up.lu_p->up.magic_p != NMG_FACEUSE_MAGIC) continue;
4754  tfu = eur->up.lu_p->up.fu_p;
4755  if (tfu->f_p->g.plane_p != fg2) continue;
4756  NMG_CK_EDGE_G_EITHER(eur->g.lseg_p);
4757 
4758  /* Found the other face on this edge! */
4759  if (RTG.NMG_debug & DEBUG_BASIC) {
4760  bu_log(" Found shared edge, eur=%p, eg=%p\n", (void *)eur, (void *)eur->g.lseg_p);
4761  nmg_pr_eu_briefly(eur, (char *)NULL);
4762  nmg_pr_eu_briefly(eur->eumate_p, (char *)NULL);
4763  nmg_pr_eg(eur->g.magic_p, 0);
4764  nmg_pr_fu_around_eu(eur, tol);
4765  }
4766 
4767  if (!ret) {
4768  /* First common edge found */
4769  if (RTG.NMG_debug & DEBUG_BASIC) {
4770  nmg_pl_lu_around_eu(eur);
4771  }
4772  ret = eur;
4773  continue;
4774  }
4775 
4776  /* Previous edge found, check edge_g */
4777  if (eur->g.lseg_p == ret->g.lseg_p) continue;
4778 
4779  /* Edge geometry differs. vu's same? */
4780  if (NMG_ARE_EUS_ADJACENT(eur, ret)) {
4781  if (RTG.NMG_debug & DEBUG_BASIC) {
4782  bu_log("nmg_find_eg_between_2fg() joining edges eur=%p, ret=%p\n",
4783  (void *)eur, (void *)ret);
4784  }
4785  nmg_radial_join_eu(ret, eur, tol);
4786  goto restart;
4787  }
4788 
4789  /* This condition "shouldn't happen" */
4790  bu_log("eur=%p, eg_p=%p; ret=%p, eg_p=%p\n",
4791  (void *)eur, (void *)eur->g.lseg_p,
4792  (void *)ret, (void *)ret->g.lseg_p);
4793  nmg_pr_eg(eur->g.magic_p, 0);
4794  nmg_pr_eg(ret->g.magic_p, 0);
4795  nmg_pr_eu_endpoints(eur, 0);
4796  nmg_pr_eu_endpoints(ret, 0);
4797 
4798  coincident = nmg_2edgeuse_g_coincident(eur, ret, tol);
4799  if (coincident) {
4800  /* Change eur to use ret's eg */
4801  bu_log("nmg_find_eg_between_2fg() belatedly fusing e1=%p, eg1=%p, e2=%p, eg2=%p\n",
4802  (void *)eur->e_p, (void *)eur->g.lseg_p,
4803  (void *)ret->e_p, (void *)ret->g.lseg_p);
4804  nmg_jeg(ret->g.lseg_p, eur->g.lseg_p);
4805  /* See if there are any others. */
4806  nmg_model_fuse(nmg_find_model(&eur->l.magic), tol);
4807  } else {
4808  bu_bomb("nmg_find_eg_between_2fg() 2 faces intersect with differing edge geometries?\n");
4809  }
4810  goto restart;
4811  }
4812  }
4813  }
4814  }
4815  if (RTG.NMG_debug & DEBUG_BASIC) {
4816  bu_log("nmg_find_eg_between_2fg(fu1=%p, fu2=%p) edge_g=%p\n",
4817  (void *)ofu1, (void *)fu2, ret ? (void *)ret->g.lseg_p : (void *)0);
4818  }
4819  if (ret)
4820  return ret->g.lseg_p;
4821  return (struct edge_g_lseg *)NULL;
4822 }
4823 
4824 
4825 /* XXX Move to nmg_info.c */
4826 /**
4827  * See if any edgeuse in the given faceuse
4828  * lies on the indicated edge geometry (edge_g).
4829  * This is a topology check only.
4830  *
4831  * Returns -
4832  * NULL No
4833  * eu Yes, here is one edgeuse that does. There may be more.
4834  */
4835 struct edgeuse *
4836 nmg_does_fu_use_eg(const struct faceuse *fu1, const uint32_t *eg)
4837 {
4838  const struct loopuse *lu1;
4839  register struct edgeuse *eu1;
4840 
4841  NMG_CK_FACEUSE(fu1);
4842  NMG_CK_EDGE_G_EITHER(eg);
4843 
4844  for (BU_LIST_FOR(lu1, loopuse, &fu1->lu_hd)) {
4845  NMG_CK_LOOPUSE(lu1);
4846  if (BU_LIST_FIRST_MAGIC(&lu1->down_hd) == NMG_VERTEXUSE_MAGIC)
4847  continue;
4848  if (RTG.NMG_debug & DEBUG_BASIC) {
4849  bu_log(" visiting lu1=%p, fu1=%p\n",
4850  (void *)lu1, (void *)fu1);
4851  }
4852  for (BU_LIST_FOR(eu1, edgeuse, &lu1->down_hd)) {
4853  if (eu1->g.magic_p == eg) goto out;
4854  }
4855  }
4856  eu1 = (struct edgeuse *)NULL;
4857 out:
4858  if (RTG.NMG_debug & DEBUG_BASIC) {
4859  bu_log("nmg_does_fu_use_eg(fu1=%p, eg=%p) eu1=%p\n",
4860  (void *)fu1, (void *)eg, (void *)eu1);
4861  }
4862  return eu1;
4863 }
4864 
4865 
4866 /* XXX move to plane.c */
4867 /**
4868  * Returns -
4869  * 1 line is on plane, within tol
4870  * 0 line does not lie on the plane
4871  */
4872 int
4873 rt_line_on_plane(const fastf_t *pt, const fastf_t *dir, const fastf_t *plane, const struct bn_tol *tol)
4874 {
4875  vect_t unitdir;
4876  fastf_t dist;
4877 
4878  BN_CK_TOL(tol);
4879 
4880  dist = DIST_PT_PLANE(pt, plane);
4881  if (!NEAR_ZERO(dist, tol->dist)) return 0;
4882 
4883  VMOVE(unitdir, dir);
4884  VUNITIZE(unitdir);
4885 /* XXX This is *way* too tight TOO_STRICT */
4886  if (fabs(VDOT(unitdir, plane)) >= tol->para) {
4887  /* Vectors are parallel */
4888  /* ray parallel to plane, and point is on it */
4889  return 1;
4890  }
4891  return 0;
4892 }
4893 
4894 
4895 /**
4896  * Handle the complete mutual intersection of
4897  * two 3-D non-coplanar planar faces,
4898  * including cutjoin and meshing.
4899  *
4900  * The line of intersection has already been computed.
4901  * Handle as two 2-D line/face intersection problems
4902  *
4903  * This is the HEART of the intersection code.
4904  */
4905 static void
4906 nmg_isect_two_face3p(struct nmg_inter_struct *is, struct faceuse *fu1, struct faceuse *fu2)
4907 {
4908  struct bu_ptbl vert_list1, vert_list2;
4909  struct bu_ptbl eu1_list; /* all eu's in fu1 */
4910  struct bu_ptbl eu2_list; /* all eu's in fu2 */
4911  fastf_t *mag1=(fastf_t *)NULL;
4912  fastf_t *mag2=(fastf_t *)NULL;
4913  int i;
4914 
4915  NMG_CK_INTER_STRUCT(is);
4916  NMG_CK_FACEUSE(fu1);
4917  NMG_CK_FACEUSE(fu2);
4918 
4919  if (RTG.NMG_debug & DEBUG_POLYSECT) {
4920  bu_log("nmg_isect_two_face3p(fu1=%p, fu2=%p) START12\n", (void *)fu1, (void *)fu2);
4921  VPRINT("isect ray is->pt ", is->pt);
4922  VPRINT("isect ray is->dir", is->dir);
4923  }
4924 
4925  if (RTG.NMG_debug & DEBUG_VERIFY) {
4926  nmg_vfu(&fu1->s_p->fu_hd, fu1->s_p);
4927  nmg_vfu(&fu2->s_p->fu_hd, fu2->s_p);
4928  nmg_fu_touchingloops(fu1);
4929  nmg_fu_touchingloops(fu2);
4930  }
4931 
4932  /* Verify that line is within tolerance of both planes */
4933 #ifdef TOO_STRICT
4934  NMG_GET_FU_PLANE(n1, fu1);
4935  if (!rt_line_on_plane(is->pt, is->dir, n1, &(is->tol)))
4936  bu_log("WARNING: intersect line not on plane of fu1\n");
4937  NMG_GET_FU_PLANE(n2, fu2);
4938  if (!rt_line_on_plane(is->pt, is->dir, n2, &(is->tol)))
4939  bu_log("WARNING: intersect line not on plane of fu2\n");
4940 #endif
4941 
4942  if (RTG.NMG_debug & (DEBUG_POLYSECT|DEBUG_FCUT|DEBUG_MESH)
4943  && RTG.NMG_debug & DEBUG_PLOTEM) {
4944  nmg_pl_2fu("Iface%d.plot3", fu1, fu2, 0);
4945  }
4946 
4947  bu_ptbl_init(&vert_list1, 64, "vert_list1 buffer");
4948  bu_ptbl_init(&vert_list2, 64, "vert_list2 buffer");
4949 
4950  /* Build list of all edgeuses in fu1 and fu2 */
4951  nmg_edgeuse_tabulate(&eu1_list, &fu1->l.magic);
4952  nmg_edgeuse_tabulate(&eu2_list, &fu2->l.magic);
4953 
4954  is->mag_len = 2 * (BU_PTBL_END(&eu1_list) + BU_PTBL_END(&eu2_list));
4955  mag1 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag1");
4956  mag2 = (fastf_t *)bu_calloc(is->mag_len, sizeof(fastf_t), "mag2");
4957 
4958  for (i=0; i<is->mag_len; i++) {
4959  mag1[i] = MAX_FASTF;
4960  mag2[i] = MAX_FASTF;
4961  }
4962 
4963 
4964  is->l1 = &vert_list1;
4965  is->l2 = &vert_list2;
4966  is->s1 = fu1->s_p;
4967  is->s2 = fu2->s_p;
4968  is->fu1 = fu1;
4969  is->fu2 = fu2;
4970  is->mag1 = mag1;
4971  is->mag2 = mag2;
4972 
4973  /*
4974  * Now do the intersection with the other face.
4975  */
4976  is->l2 = &vert_list1;
4977  is->l1 = &vert_list2;
4978  is->s2 = fu1->s_p;
4979  is->s1 = fu2->s_p;
4980  is->fu2 = fu1;
4981  is->fu1 = fu2;
4982  is->mag1 = mag2;
4983  is->mag2 = mag1;
4984 /* nmg_isect_line2_face2pNEW(is, fu2, fu1, &eu2_list, &eu1_list); */
4985  is->on_eg = (struct edge_g_lseg *)NULL;
4986  nmg_isect_fu_jra(is, fu1, fu2, &eu1_list, &eu2_list);
4987 
4988  if (RTG.NMG_debug & DEBUG_VERIFY) {
4989  nmg_fu_touchingloops(fu1);
4990  nmg_fu_touchingloops(fu2);
4991  nmg_vfu(&fu1->s_p->fu_hd, fu1->s_p);
4992  nmg_vfu(&fu2->s_p->fu_hd, fu2->s_p);
4993  }
4994 
4995  if (RTG.NMG_debug & DEBUG_FCUT) {
4996  bu_log("nmg_isect_two_face3p(fu1=%p, fu2=%p) vert_lists B:\n", (void *)fu1, (void *)fu2);
4997  nmg_pr_ptbl_vert_list("vert_list1", &vert_list1, mag1);
4998  nmg_pr_ptbl_vert_list("vert_list2", &vert_list2, mag2);
4999  }
5000 
5001  if (vert_list1.end == 0 && vert_list2.end == 0) {
5002  /* there were no intersections */
5003  goto out;
5004  }
5005 
5006  if (RTG.NMG_debug & DEBUG_POLYSECT) {
5007  bu_log("nmg_isect_two_face3p(fu1=%p, fu2=%p) MIDDLE\n", (void *)fu1, (void *)fu2);
5008  }
5009 
5010  is->on_eg = nmg_face_cutjoin(&vert_list1, &vert_list2, mag1, mag2, fu1, fu2, is->pt, is->dir, is->on_eg, &is->tol);
5011 
5012  if (RTG.NMG_debug & DEBUG_VERIFY) {
5013  nmg_fu_touchingloops(fu1);
5014  nmg_fu_touchingloops(fu2);
5015  nmg_region_v_unique(fu1->s_p->r_p, &is->tol);
5016  nmg_region_v_unique(fu2->s_p->r_p, &is->tol);
5017  nmg_vfu(&fu1->s_p->fu_hd, fu1->s_p);
5018  nmg_vfu(&fu2->s_p->fu_hd, fu2->s_p);
5019  }
5020 
5021  nmg_mesh_faces(fu1, fu2, &is->tol);
5022  if (RTG.NMG_debug & DEBUG_VERIFY) {
5023  nmg_fu_touchingloops(fu1);
5024  nmg_fu_touchingloops(fu2);
5025  }
5026 
5027 out:
5028  (void)bu_ptbl_free(&vert_list1);
5029  (void)bu_ptbl_free(&vert_list2);
5030  (void)bu_ptbl_free(&eu1_list);
5031  (void)bu_ptbl_free(&eu2_list);
5032  if (mag1)
5033  bu_free((char *)mag1, "nmg_isect_two_face3p: mag1");
5034  if (mag2)
5035  bu_free((char *)mag2, "nmg_isect_two_face3p: mag2");
5036 
5037 
5038  if (RTG.NMG_debug & DEBUG_VERIFY) {
5039  nmg_vfu(&fu1->s_p->fu_hd, fu1->s_p);
5040  nmg_vfu(&fu2->s_p->fu_hd, fu2->s_p);
5041  }
5042  if (RTG.NMG_debug & DEBUG_POLYSECT) {
5043  bu_log("nmg_isect_two_face3p(fu1=%p, fu2=%p) END\n", (void *)fu1, (void *)fu2);
5044  VPRINT("isect ray is->pt ", is->pt);
5045  VPRINT("isect ray is->dir", is->dir);
5046  }
5047 }
5048 
5049 
5050 void
5051 nmg_cut_lu_into_coplanar_and_non(struct loopuse *lu, fastf_t *pl, struct nmg_inter_struct *is)
5052 {
5053  struct model *m;
5054  struct edgeuse *eu;
5055  struct vertex_g *vg2;
5056  struct vertex *vcut1, *vcut2;
5057  struct bu_ptbl cut_list;
5058  fastf_t dist, dist2;
5059  int class1, class2;
5060  int in=0;
5061  int on=0;
5062  int out=0;
5063  int i;
5064 
5065  if (RTG.NMG_debug & DEBUG_POLYSECT)
5066  bu_log("nmg_cut_lu_into_coplanar_and_non(lu=%p, pl=%g %g %g %g)\n", (void *)lu, V4ARGS(pl));
5067 
5068  NMG_CK_LOOPUSE(lu);
5069  NMG_CK_INTER_STRUCT(is);
5070 
5071  m = nmg_find_model(&is->fu1->l.magic);
5072 
5073  /* check if this loop even needs to be considered */
5074  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
5075  return;
5076 
5077  bu_ptbl_init(&cut_list, 64, " &cut_list");
5078 
5079  eu = BU_LIST_FIRST(edgeuse, &lu->down_hd);
5080  vg2 = eu->vu_p->v_p->vg_p;
5081  NMG_CK_VERTEX_G(vg2);
5082  dist2 = DIST_PT_PLANE(vg2->coord, pl);
5083  if (dist2 > is->tol.dist) {
5084  class2 = NMG_CLASS_AoutB;
5085  out++;
5086  } else if (dist2 < (-is->tol.dist)) {
5087  class2 = NMG_CLASS_AinB;
5088  in++;
5089  } else {
5090  class2 = NMG_CLASS_AonBshared;
5091  on++;
5092  }
5093 
5094  vcut1 = (struct vertex *)NULL;
5095  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
5096  class1 = class2;
5097 
5098  vg2 = eu->eumate_p->vu_p->v_p->vg_p;
5099  dist2 = DIST_PT_PLANE(vg2->coord, pl);
5100 
5101  if (dist2 > is->tol.dist) {
5102  class2 = NMG_CLASS_AoutB;
5103  out++;
5104  } else if (dist2 < (-is->tol.dist)) {
5105  class2 = NMG_CLASS_AinB;
5106  in++;
5107  } else {
5108  class2 = NMG_CLASS_AonBshared;
5109  on++;
5110  }
5111 
5112  if (class1 == NMG_CLASS_AonBshared && class2 != class1) {
5113  if (!vcut1) {
5114  vcut1 = eu->vu_p->v_p;
5115  } else if (vcut1 != eu->vu_p->v_p) {
5116  bu_ptbl_ins(&cut_list, (long *)vcut1);
5117  bu_ptbl_ins(&cut_list, (long *)eu->vu_p->v_p);
5118  vcut1 = (struct vertex *)NULL;
5119  }
5120  } else if (class2 == NMG_CLASS_AonBshared && class1 != class2) {
5121  if (!vcut1) {
5122  vcut1 = eu->eumate_p->vu_p->v_p;
5123  } else if (vcut1 != eu->eumate_p->vu_p->v_p) {
5124  bu_ptbl_ins(&cut_list, (long *)vcut1);
5125  bu_ptbl_ins(&cut_list, (long *)eu->eumate_p->vu_p->v_p);
5126  vcut1 = (struct vertex *)NULL;
5127  }
5128  }
5129  }
5130 
5131  if (RTG.NMG_debug & DEBUG_POLYSECT) {
5132  bu_log("\t pl=(%g %g %g %g)\n", V4ARGS(pl));
5133  bu_log("\tcut_lists=%ld, on=%d, in=%d, out=%d\n", BU_PTBL_END(&cut_list), on, in, out);
5134  if (BU_PTBL_END(&cut_list)) {
5135  bu_log("\tcut_lists:\n");
5136  for (i=0; i<BU_PTBL_END(&cut_list); i++) {
5137  struct vertex *v;
5138 
5139  v = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5140  bu_log("\t\t%d, %p\n", i+1, (void *)v);
5141  }
5142  }
5143  }
5144 
5145  if (!on)
5146  return;
5147 
5148  if (BU_PTBL_END(&cut_list) < 2) {
5149  bu_ptbl_free(&cut_list);
5150 
5151  if (RTG.NMG_debug & DEBUG_POLYSECT)
5152  bu_log("No loops need cutting\n");
5153  return;
5154  }
5155 
5156  if (nmg_loop_is_a_crack(lu)) {
5157  struct bu_ptbl lus;
5158 
5159  if (RTG.NMG_debug & DEBUG_POLYSECT)
5160  bu_log("Loop is a crack\n");
5161 
5162  i = 0;
5163  while (i < BU_PTBL_END(&cut_list)) {
5164  struct vertexuse *vu;
5165 
5166  vcut1 = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5167  for (BU_LIST_FOR(vu, vertexuse, &vcut1->vu_hd)) {
5168  if (nmg_find_lu_of_vu(vu) != lu)
5169  continue;
5170 
5171  eu = vu->up.eu_p;
5172  if (NMG_ARE_EUS_ADJACENT(eu, BU_LIST_PNEXT_CIRC(edgeuse, &eu->l))) {
5173  i--;
5174  bu_ptbl_rm(&cut_list, (long *)vcut1);
5175  } else if (NMG_ARE_EUS_ADJACENT(eu, BU_LIST_PPREV_CIRC(edgeuse, &eu->l))) {
5176  i--;
5177  bu_ptbl_rm(&cut_list, (long *)vcut1);
5178  }
5179  }
5180  i++;
5181  }
5182 
5183  if (BU_PTBL_END(&cut_list) == 0) {
5184  bu_ptbl_free(&cut_list);
5185 
5186  if (RTG.NMG_debug & DEBUG_POLYSECT)
5187  bu_log("no loops need cutting\n");
5188  return;
5189  }
5190 
5191  bu_ptbl_init(&lus, 64, " &lus");
5192  bu_ptbl_ins(&lus, (long *)lu);
5193  for (i=0; i<BU_PTBL_END(&cut_list); i++) {
5194  int j;
5195 
5196  vcut1 = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5197 
5198  for (j=0; j<BU_PTBL_END(&lus); j++) {
5199  int did_split=0;
5200  struct loopuse *lu1;
5201  struct vertexuse *vu1;
5202  struct loopuse *new_lu;
5203 
5204  lu1 = (struct loopuse *)BU_PTBL_GET(&lus, j);
5205 
5206  for (BU_LIST_FOR(vu1, vertexuse, &vcut1->vu_hd)) {
5207  if (nmg_find_lu_of_vu(vu1) == lu1) {
5208  if (RTG.NMG_debug & DEBUG_POLYSECT)
5209  bu_log("Splitting lu %p at vu %p\n", (void *)lu1, (void *)vu1);
5210  new_lu = nmg_split_lu_at_vu(lu1, vu1);
5211  nmg_lu_reorient(lu1);
5212  nmg_lu_reorient(new_lu);
5213  nmg_loop_g(new_lu->l_p, &is->tol);
5214  nmg_loop_g(lu1->l_p, &is->tol);
5215  bu_ptbl_ins(&lus, (long *)new_lu);
5216  did_split = 1;
5217  break;
5218  }
5219  }
5220  if (did_split)
5221  break;
5222  }
5223  }
5224  bu_ptbl_free(&lus);
5225  bu_ptbl_free(&cut_list);
5226  return;
5227  }
5228 
5229  if (BU_PTBL_END(&cut_list)%2) {
5230  bu_log("Uneven number (%ld) of vertices on cut list\n", BU_PTBL_END(&cut_list));
5231  bu_bomb("Uneven number of vertices on cut list");
5232  }
5233 
5234  /* Sort vertices on cut list into some order */
5235  if (BU_PTBL_END(&cut_list) > 2) {
5236  struct vertex *v1=(struct vertex *)NULL;
5237  struct vertex *v2=(struct vertex *)NULL;
5238  struct vertex *end1 = NULL, *end2 = NULL;
5239  fastf_t max_dist = 0.0;
5240  vect_t diff;
5241  fastf_t *dist_array;
5242  int done;
5243 
5244  /* find longest distance between two vertices */
5245  for (i = 0; i < BU_PTBL_END(&cut_list); i++) {
5246  int j;
5247 
5248  v1 = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5249 
5250  for (j=i; j<BU_PTBL_END(&cut_list); j++) {
5251  fastf_t tmp_dist;
5252 
5253  v2 = (struct vertex *)BU_PTBL_GET(&cut_list, j);
5254  VSUB2(diff, v1->vg_p->coord, v2->vg_p->coord);
5255  tmp_dist = MAGSQ(diff);
5256  if (tmp_dist > max_dist) {
5257  max_dist = tmp_dist;
5258  end1 = v1;
5259  end2 = v2;
5260  }
5261  }
5262  }
5263  if (!end1 || !end2) {
5264  bu_log("nmg_cut_lu_into_coplanar_and_non: Cannot find endpoints\n");
5265  bu_bomb("nmg_cut_lu_into_coplanar_and_non: Cannot find endpoints\n");
5266  }
5267 
5268  /* create array of distances (SQ) along the line from end1 to end2 */
5269  dist_array = (fastf_t *)bu_calloc(sizeof(fastf_t), BU_PTBL_END(&cut_list), "distance array");
5270  for (i=0; i<BU_PTBL_END(&cut_list); i++) {
5271  v1 = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5272  if (v1 == end1) {
5273  dist_array[i] = 0.0;
5274  continue;
5275  }
5276  if (v1 == end2) {
5277  dist_array[i] = max_dist;
5278  continue;
5279  }
5280 
5281  VSUB2(diff, v1->vg_p->coord, end1->vg_p->coord);
5282  dist_array[i] = MAGSQ(diff);
5283  }
5284 
5285  /* sort vertices according to distance array */
5286  done = 0;
5287  while (!done) {
5288  fastf_t tmp_dist;
5289  long *tmp_v;
5290 
5291  done = 1;
5292  for (i=1; i<BU_PTBL_END(&cut_list); i++) {
5293  if (dist_array[i-1] <= dist_array[i])
5294  continue;
5295 
5296  /* swap distances in array */
5297  tmp_dist = dist_array[i];
5298  dist_array[i] = dist_array[i-1];
5299  dist_array[i-1] = tmp_dist;
5300 
5301  /* swap vertices */
5302  tmp_v = cut_list.buffer[i];
5303  cut_list.buffer[i] = cut_list.buffer[i-1];
5304  cut_list.buffer[i-1] = tmp_v;
5305 
5306  done = 0;
5307  }
5308  }
5309 
5310  if (RTG.NMG_debug & DEBUG_POLYSECT) {
5311  bu_log("After sorting:\n");
5312  for (i = 0; i < BU_PTBL_END(&cut_list); i++)
5313  bu_log("v=%p, dist=%g\n", (void *)BU_PTBL_GET(&cut_list, i), dist_array[i]);
5314  }
5315 
5316  bu_free((char *)dist_array, "distance array");
5317  }
5318 
5319  for (i=0; i<BU_PTBL_END(&cut_list); i += 2) {
5320  struct loopuse *lu1;
5321  struct vertexuse *vu;
5322  vect_t dir;
5323  point_t hit_pt;
5324  struct vertex *hit_v;
5325  struct vertexuse *hit_vu;
5326  struct edgeuse *new_eu;
5327  fastf_t len;
5328  fastf_t inv_len;
5329  int skip=0;
5330 
5331  vcut1 = (struct vertex *)BU_PTBL_GET(&cut_list, i);
5332  vcut2 = (struct vertex *)BU_PTBL_GET(&cut_list, i+1);
5333 
5334  if (vcut1 == vcut2)
5335  continue;
5336 
5337  /* make sure these are not the ends of an edge of lu */
5338  for (BU_LIST_FOR(vu, vertexuse, &vcut1->vu_hd)) {
5339  if (nmg_find_lu_of_vu(vu) != lu)
5340  continue;
5341 
5342  eu = vu->up.eu_p;
5343  if (eu->eumate_p->vu_p->v_p == vcut2) {
5344  /* already an edge here */
5345  skip = 1;
5346  break;
5347  }
5348  }
5349  if (skip)
5350  continue;
5351 
5352  /* need to cut face along line from vcut1 to vcut2 */
5353  VMOVE(is->pt, vcut1->vg_p->coord);
5354  VSUB2(dir, vcut2->vg_p->coord, vcut1->vg_p->coord);
5355  len = MAGNITUDE(dir);
5356  if (len <= is->tol.dist)
5357  continue;
5358 
5359  inv_len = 1.0/len;
5360  VSCALE(is->dir, dir, inv_len);
5361 
5362  /* add vertexuses to intersect list */
5363  bu_ptbl_reset(is->l1);
5364 
5365  /* add uses of vcut1 */
5366  for (BU_LIST_FOR(vu, vertexuse, &vcut1->vu_hd)) {
5367  if (nmg_find_fu_of_vu(vu) == is->fu1)
5368  nmg_enlist_one_vu(is, vu, 0.0);
5369  }
5370 
5371  /* add uses of vcut2 */
5372  for (BU_LIST_FOR(vu, vertexuse, &vcut2->vu_hd)) {
5373  if (nmg_find_fu_of_vu(vu) == is->fu1)
5374  nmg_enlist_one_vu(is, vu, len);
5375  }
5376 
5377  /* look for other edges that may intersect this line */
5378  for (BU_LIST_FOR(lu1, loopuse, &is->fu1->lu_hd)) {
5379  struct edgeuse *eu1;
5380