BRL-CAD
nurb_trim.c
Go to the documentation of this file.
1 /* N U R B _ T R I M . C
2  * BRL-CAD
3  *
4  * Copyright (c) 1990-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 */
21 /** @{ */
22 /** @file primitives/bspline/nurb_trim.c
23  *
24  * Trimming curve routines.
25  *
26  */
27 /** @} */
28 
29 #include "common.h"
30 
31 #include <math.h>
32 #include "bio.h"
33 
34 #include "vmath.h"
35 #include "raytrace.h"
36 #include "nurb.h"
37 
38 extern void rt_clip_cnurb(struct bu_list *plist, struct edge_g_cnurb *crv, fastf_t u, fastf_t v);
39 
41  int axis;
43 };
44 
45 
46 #define QUAD1 0
47 #define QUAD2 1
48 #define QUAD3 2
49 #define QUAD4 3
50 
51 
52 #define TRIM_OUT 0
53 #define TRIM_IN 1
54 #define TRIM_ON 2
55 
56 /* The following defines need to be 0, 2, 3 in order to drive the quad
57  * table and determine the appropriate case for processing the
58  * trimming curve.
59  */
60 
61 #define CASE_A 0
62 #define CASE_B 2
63 #define CASE_C 3
64 
65 static int quad_table[16] = {
66  /* A = 0, B = 2, C = 3 */
67  0, 0, 0, 0, 0, 3, 0, 3, 0, 2, 3, 3, 0, 3, 3, 3
68 };
69 
70 
71 /**
72  * This routine determines what quadrants the trimming curves lies in,
73  * It then uses a table look up to determine the whether its CASE{A,
74  * B, C}, One difference from the paper is the fact that if any of the
75  * points lie on the axis of the u, v quadrant system then the axis is
76  * only in either Quadrant 1 or Quadrant 2 and not q3 or q4. This
77  * handles the case of endpoint problems correctly.
78  */
79 int
80 rt_trim_case(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
81 {
82  int quadrant;
83  int qstats;
84  fastf_t * pts;
85  int coords, rat;
86  int i;
87 
88  qstats = 0;
89 
90  coords = RT_NURB_EXTRACT_COORDS(trim->pt_type);
91  pts = trim->ctl_points;
92  rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);
93 
94  /* Handle rational specially since we need to divide the rational
95  * portion.
96  */
97 
98  if (rat) {
99  for (i = 0; i < trim->c_size; i++) {
100  if (pts[0]/pts[2] > u)
101  quadrant = (pts[1]/pts[2] >= v)? QUAD1:QUAD4;
102  else
103  quadrant = (pts[1]/pts[2] >= v)? QUAD2:QUAD3;
104 
105  qstats |= (1 << quadrant);
106  pts += coords;
107  }
108  } else {
109  for (i = 0; i < trim->c_size; i++) {
110  if (pts[0] > u)
111  quadrant = (pts[1] >= v)? QUAD1:QUAD4;
112  else
113  quadrant = (pts[1] >= v)? QUAD2:QUAD3;
114 
115  qstats |= (1 << quadrant);
116  pts += coords;
117  }
118  }
119 
120  return quad_table[qstats]; /* return the special case of the curve */
121 }
122 
123 
124 /**
125  * Process Case B curves.
126  *
127  * If the two endpoints of the curve lie in different quadrants than
128  * the axis crosses the curve an odd number of times
129  * (TRIM_IN). Otherwise the curve crosses the u, v axis a even number
130  * of times (TRIM_OUT). No further processing is required.
131  */
132 int
133 rt_process_caseb(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
134 {
135  int q1, q2;
136  fastf_t * pts;
137  int rat;
138 
139  rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);
140 
141  pts = trim->ctl_points;
142 
143  if (rat) {
144  if (pts[0]/pts[2] > u) q1 = (pts[1]/pts[2] >= v)?QUAD1:QUAD4;
145  else q1 = (pts[1]/pts[2] >= v)?QUAD2:QUAD3;
146 
147 
148  pts = trim->ctl_points + RT_NURB_EXTRACT_COORDS(trim->pt_type) *
149  (trim->c_size -1);
150  if (pts[0]/pts[2] > u) q2 = (pts[1]/pts[2] >= v)?QUAD1:QUAD4;
151  else q2 = (pts[1]/pts[2] >= v)?QUAD2:QUAD3;
152 
153  } else {
154  if (pts[0] > u) q1 = (pts[1] >= v)?QUAD1:QUAD4;
155  else q1 = (pts[1] >= v)?QUAD2:QUAD3;
156 
157 
158  pts = trim->ctl_points +
159  RT_NURB_EXTRACT_COORDS(trim->pt_type) *
160  (trim->c_size -1);
161  if (pts[0] > u) q2 = (pts[1] >= v)?QUAD1:QUAD4;
162  else q2 = (pts[1] >= v)?QUAD2:QUAD3;
163  }
164 
165  if (q1 != q2)
166  return TRIM_IN;
167  else
168  return TRIM_OUT;
169 
170 }
171 
172 
173 /**
174  * Only check end points of the curve
175  */
176 int
177 rt_nurb_uv_dist(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
178 {
179 
180  fastf_t dist;
181  fastf_t * ptr;
182  int coords;
183  int rat;
184  fastf_t u2, v2;
185 
186  ptr = trim->ctl_points;
187  coords = RT_NURB_EXTRACT_COORDS(trim->pt_type);
188  rat = RT_NURB_IS_PT_RATIONAL(trim->pt_type);
189 
190  u2 = 0.0;
191  v2 = 0.0;
192 
193  if (rat) {
194  u2 = ptr[0]/ptr[2] - u; u2 *= u2;
195  v2 = ptr[1]/ptr[2] - v; v2 *= v2;
196  } else {
197  u2 = ptr[0] - u; u2 *= u2;
198  v2 = ptr[1] - v; v2 *= v2;
199  }
200 
201  dist = sqrt(u2 + v2);
202  if (NEAR_ZERO(dist, 1.0e-4))
203  return TRIM_ON;
204 
205  ptr = trim->ctl_points + coords * (trim->c_size -1);
206 
207  u2 = 0.0;
208  v2 = 0.0;
209 
210  if (rat) {
211  u2 = ptr[0]/ptr[2] - u; u2 *= u2;
212  v2 = ptr[1]/ptr[2] - v; v2 *= v2;
213  } else {
214  u2 = ptr[0] - u; u2 *= u2;
215  v2 = ptr[1] - v; v2 *= v2;
216  }
217 
218  dist = sqrt(u2 + v2);
219  if (NEAR_ZERO(dist, 1.0e-4))
220  return TRIM_ON;
221 
222  return TRIM_OUT;
223 
224 }
225 
226 
227 /**
228  * Process Case C curves;
229  *
230  * A check is placed here to determine if the u, v is on the curve.
231  * Determine how many times the curve will cross the u, v axis. If the
232  * curve crosses an odd number of times than the point is IN, else the
233  * point is OUT. Since a Case C curve need processing a call to clip
234  * the curve so that it becomes either Case B or Case A is required
235  * to determine the number of crossings accurately. Thus we need to keep
236  * the original curve and expect the calling routine to free the
237  * storage. Additional curves are generated in this routine, each of
238  * these new curves are processed, and then are deleted before exiting
239  * this procedure.
240  */
241 int
242 rt_process_casec(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
243 {
244 
245  struct edge_g_cnurb * clip;
246  int jordan_hit;
247  struct bu_list plist;
248  int trim_flag = 0;
249  int caset;
250 
251  /* determine if the u, v values are on the curve */
252 
253  if (rt_nurb_uv_dist(trim, u, v) == TRIM_ON) return TRIM_IN;
254 
255  jordan_hit = 0;
256 
257  BU_LIST_INIT(&plist);
258 
259  if (nurb_crv_is_bezier(trim))
260  rt_clip_cnurb(&plist, trim, u, v);
261  else
262  nurb_c_to_bezier(&plist, trim);
263 
264  while (BU_LIST_WHILE(clip, edge_g_cnurb, &plist)) {
265  BU_LIST_DEQUEUE(&clip->l);
266 
267  caset = rt_trim_case(clip, u, v);
268 
269  trim_flag = 0;
270 
271  if (caset == CASE_B)
272  trim_flag = rt_process_caseb(clip, u, v);
273  if (caset == CASE_C)
274  trim_flag = rt_process_casec(clip, u, v);
275 
276  rt_nurb_free_cnurb(clip);
277 
278  if (trim_flag == TRIM_IN) jordan_hit++;
279  if (trim_flag == TRIM_ON) break;
280  }
281 
282  while (BU_LIST_WHILE(clip, edge_g_cnurb, &plist)) {
283  BU_LIST_DEQUEUE(&clip->l);
284  rt_nurb_free_cnurb(clip);
285  }
286 
287  if (trim_flag == TRIM_ON)
288  return TRIM_ON;
289 
290  else if (jordan_hit & 01)
291  return TRIM_IN;
292  else
293  return TRIM_OUT;
294 }
295 
296 
297 /**
298  * This routine will be called several times, once for each portion of
299  * the trimming curve. It returns whether a line extended from the
300  * <u, v> point will cross the trimming curve an even or odd number of
301  * times. Or the <u, v> point could be on the curve in which case
302  * TRIM_ON will be returned. The algorithm uses the approach taken Tom
303  * Sederburge and uses bezier clipping to produce caseA and caseB
304  * curves. If the original trimming curve is a CASE C curve then
305  * further processing is required.
306  */
307 int
308 rt_uv_in_trim(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
309 {
310 
311  int quad_case;
312 
313  quad_case = rt_trim_case(trim, u, v); /* determine quadrants */
314 
315  /* CASE A */
316  if (quad_case == CASE_A)
317  return TRIM_OUT;
318  if (quad_case == CASE_B) /* CASE B */
319  return rt_process_caseb(trim, u, v);
320  if (quad_case == CASE_C) /* CASE C */
321  return rt_process_casec(trim, u, v);
322 
323  bu_log("rt_uv_in_trim: rt_trim_case() returned illegal value %d\n", quad_case);
324  return -1;
325 }
326 
327 
328 /**
329  * This routines is used to determine how far a point is from the u, v
330  * quadrant axes.
331  *
332  * Equations 3, 4, 5 in Sederberg '90 paper
333  */
334 fastf_t
335 rt_trim_line_pt_dist(struct _interior_line *l, fastf_t *pt, int pt_type)
336 {
337  fastf_t h;
338  int h_flag;
339 
340  h_flag = RT_NURB_IS_PT_RATIONAL(pt_type);
341 
342  if (l->axis == 0) {
343  if (h_flag) h = (pt[1] / pt[2] - l->o_dist) * pt[2]; /* pt[2] is weight */
344  else h = pt[1] - l->o_dist;
345 
346  } else {
347  if (h_flag) h = (pt[0] / pt[2] - l->o_dist) * pt[2];
348  else h = pt[0] - l->o_dist;
349 
350  }
351 
352  return h;
353 }
354 
355 
356 /**
357  * Return the SIGN of the value
358  */
359 int
361 {
362  if (f < 0.0)
363  return -1;
364  else
365  return 1;
366 
367 }
368 
369 
370 /**
371  * We try to clip a curve so that it can be either Case A or Case C.
372  * Sometimes one of the curves is still Case C though, but it is much
373  * smaller than the original, and further clipping will either show that
374  * it is on the curve or provide all Case B or Case A curves. We try
375  * to pick the best axis to clip against, but this may not always
376  * work. One extra step that was included, that is not in the paper
377  * for curves but is for surfaces, is the fact that sometimes the
378  * curve is not clipped enough, if the maximum clip is less than .2
379  * then we subdivide the curve in three equal parts, at .3 and .6 .
380  * Subdivision is done using the Oslo Algorithm, rather than the other
381  * methods which were prossed.
382  */
383 void
384 rt_clip_cnurb(struct bu_list *plist, struct edge_g_cnurb *crv, fastf_t u, fastf_t v)
385 {
386  fastf_t ds1, dt1;
387  struct _interior_line s_line, t_line;
388  int axis, i;
389  fastf_t umin, umax;
390  int coords;
391  struct edge_g_cnurb * c1, *c2, *tmp;
392  fastf_t m1, m2;
393  int zero_changed;
394  fastf_t *ptr;
395  fastf_t dist[10];
396 
397  coords = RT_NURB_EXTRACT_COORDS(crv->pt_type);
398 
399  s_line.axis = 0; s_line.o_dist = v;
400  t_line.axis = 1; t_line.o_dist = u;
401 
402  ds1 = 0.0;
403  dt1 = 0.0;
404 
405  ptr = crv->ctl_points;
406 
407 
408  /* determine what axis to clip against */
409 
410  for (i = 0; i < crv->c_size; i++, ptr += coords) {
411  ds1 +=
412  fabs(rt_trim_line_pt_dist(&s_line, ptr, crv->pt_type));
413  dt1 +=
414  fabs(rt_trim_line_pt_dist(&t_line, ptr, crv->pt_type));
415  }
416 
417  if (ds1 >= dt1) axis = 0; else axis = 1;
418 
419  ptr = crv->ctl_points;
420 
421  for (i = 0; i < crv->c_size; i++) {
422  if (axis == 1)
423  dist[i] = rt_trim_line_pt_dist(&t_line, ptr, crv->pt_type);
424  else
425  dist[i] = rt_trim_line_pt_dist(&s_line, ptr, crv->pt_type);
426 
427  ptr += coords;
428  }
429 
430  /* Find the convex hull of the distances and determine the minimum
431  * and maximum distance to clip against. See the paper for details
432  * about this step
433  */
434 
435  umin = 10e40;
436  umax = -10e40;
437  zero_changed = 0;
438 
439  for (i = 0; i < crv->c_size; i++) {
440  fastf_t d1, d2;
441  fastf_t x0, x1, zero;
442 
443  if (i == (crv->c_size -1)) {
444  d1 = dist[i];
445  d2 = dist[0];
446  x0 = (fastf_t) i / (fastf_t) (crv->c_size - 1);
447  x1 = 0.0;
448  } else {
449  d1 = dist[i];
450  d2 = dist[i+1];
451  x0 = (fastf_t) i / (fastf_t) (crv->c_size - 1);
452  x1 = (i+1.0) / (crv->c_size - 1);
453  }
454 
455  if (_SIGN(d1) != _SIGN(d2)) {
456  zero = x0 - d1 * (x1 - x0)/ (d2-d1);
457  if (zero <= umin)
458  umin = zero * .99;
459  if (zero >= umax)
460  umax = zero * .99 + .01;
461  zero_changed = 1;
462  }
463  }
464 
465  if (!zero_changed)
466  return;
467 
468  /* Clip is not large enough, split in thirds and try again */
469 
470  if (umax - umin < .2) {
471  umin = .3; umax = .6;
472  }
473 
474  /* Translate the 0.0-->1.09 clipping against the real knots */
475 
476  m1 = (crv->k.knots[0] * (1 - umin)) +
477  crv->k.knots[crv->k.k_size -1] * umin;
478 
479  m2 = (crv->k.knots[0] * (1-umax)) +
480  crv->k.knots[crv->k.k_size -1] * umax;
481 
482  /* subdivide the curve */
483  c1 = (struct edge_g_cnurb *) rt_nurb_c_xsplit(crv, m1);
484  c2 = rt_nurb_c_xsplit((struct edge_g_cnurb *) c1->l.forw, m2);
485 
486  tmp = (struct edge_g_cnurb *) c1->l.forw;
487  BU_LIST_DEQUEUE(&tmp->l);
488  rt_nurb_free_cnurb(tmp);
489 
490  BU_LIST_INIT(plist);
491  BU_LIST_INSERT(&c2->l, plist);
492  BU_LIST_APPEND(plist, &c1->l);
493 }
494 
495 
496 int
497 nmg_uv_in_lu(const fastf_t u, const fastf_t v, const struct loopuse *lu)
498 {
499  struct edgeuse *eu;
500  int crossings=0;
501 
502  NMG_CK_LOOPUSE(lu);
503 
504  if (BU_LIST_FIRST_MAGIC(&lu->down_hd) != NMG_EDGEUSE_MAGIC)
505  return 0;
506 
507  for (BU_LIST_FOR(eu, edgeuse, &lu->down_hd)) {
508  struct edge_g_cnurb *eg;
509 
510  if (!eu->g.magic_p) {
511  bu_log("nmg_uv_in_lu: eu (%p) has no geometry!!!\n", (void *)eu);
512  bu_bomb("nmg_uv_in_lu: eu has no geometry!!!\n");
513  }
514 
515  if (*eu->g.magic_p != NMG_EDGE_G_CNURB_MAGIC) {
516  bu_log("nmg_uv_in_lu: Called with lu (%p) containing eu (%p) that is not CNURB!!!!\n",
517  (void *)lu, (void *)eu);
518  bu_bomb("nmg_uv_in_lu: Called with lu containing eu that is not CNURB!!!\n");
519  }
520 
521  eg = eu->g.cnurb_p;
522 
523  if (eg->order <= 0) {
524  struct vertexuse *vu1, *vu2;
525  struct vertexuse_a_cnurb *vua1, *vua2;
526  point_t uv1, uv2;
527  fastf_t slope, intersept;
528  fastf_t u_on_curve;
529 
530  vu1 = eu->vu_p;
531  vu2 = eu->eumate_p->vu_p;
532 
533  if (!vu1->a.magic_p || !vu2->a.magic_p) {
534  bu_log("nmg_uv_in_lu: Called with lu (%p) containing vu with no attribute!!!!\n",
535  (void *)lu);
536  bu_bomb("nmg_uv_in_lu: Called with lu containing vu with no attribute!!!\n");
537  }
538 
539  if (*vu1->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC
540  || *vu2->a.magic_p != NMG_VERTEXUSE_A_CNURB_MAGIC)
541  {
542  bu_log("nmg_uv_in_lu: Called with lu (%p) containing vu that is not CNURB!!!!\n",
543  (void *)lu);
544  bu_bomb("nmg_uv_in_lu: Called with lu containing vu that is not CNURB!!!\n");
545  }
546 
547  vua1 = vu1->a.cnurb_p;
548  vua2 = vu2->a.cnurb_p;
549 
550  VMOVE(uv1, vua1->param);
551  VMOVE(uv2, vua2->param);
552 
553  if (RT_NURB_IS_PT_RATIONAL(eg->pt_type)) {
554  uv1[0] /= uv1[2];
555  uv1[1] /= uv1[2];
556  uv2[0] /= uv2[2];
557  uv2[1] /= uv2[2];
558  }
559 
560  if (uv1[1] < v && uv2[1] < v)
561  continue;
562  if (uv1[1] > v && uv2[1] > v)
563  continue;
564  if (uv1[0] <= u && uv2[0] <= u)
565  continue;
566  if (ZERO(uv1[0] - uv2[0])) {
567  if ((uv1[1] <= v && uv2[1] >= v) ||
568  (uv2[1] <= v && uv1[1] >= v))
569  crossings++;
570 
571  continue;
572  }
573 
574  /* need to calculate intersection */
575  slope = (uv1[1] - uv2[1])/(uv1[0] - uv2[0]);
576  intersept = uv1[1] - slope * uv1[0];
577  u_on_curve = (v - intersept)/slope;
578  if (u_on_curve > u)
579  crossings++;
580  } else
581  crossings += rt_uv_in_trim(eg, u, v);
582  }
583 
584  if (crossings & 01)
585  return 1;
586  else
587  return 0;
588 }
589 
590 
591 /*
592  * Local Variables:
593  * mode: C
594  * tab-width: 8
595  * indent-tabs-mode: t
596  * c-file-style: "stroustrup"
597  * End:
598  * ex: shiftwidth=4 tabstop=8
599  */
#define BU_LIST_FOR(p, structure, hp)
Definition: list.h:365
#define NMG_EDGEUSE_MAGIC
Definition: magic.h:120
void bu_log(const char *,...) _BU_ATTR_PRINTF12
Definition: log.c:176
#define BU_LIST_INSERT(old, new)
Definition: list.h:183
int rt_trim_case(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
Definition: nurb_trim.c:80
#define QUAD4
Definition: nurb_trim.c:49
Definition: list.h:118
#define TRIM_OUT
Definition: nurb_trim.c:52
#define TRIM_ON
Definition: nurb_trim.c:54
int rt_process_caseb(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
Definition: nurb_trim.c:133
lu
Definition: nmg_mod.c:3855
#define QUAD1
Definition: nurb_trim.c:46
#define CASE_A
Definition: nurb_trim.c:61
int rt_nurb_uv_dist(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
Definition: nurb_trim.c:177
int rt_process_casec(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
Definition: nurb_trim.c:242
Header file for the BRL-CAD common definitions.
#define BU_LIST_APPEND(old, new)
Definition: list.h:197
#define QUAD2
Definition: nurb_trim.c:47
int _SIGN(fastf_t f)
Definition: nurb_trim.c:360
#define TRIM_IN
Definition: nurb_trim.c:53
#define QUAD3
Definition: nurb_trim.c:48
NMG_CK_LOOPUSE(lu)
void nurb_c_to_bezier(struct bu_list *clist, struct edge_g_cnurb *crv)
Definition: nurb_bezier.c:134
int rt_uv_in_trim(struct edge_g_cnurb *trim, fastf_t u, fastf_t v)
Definition: nurb_trim.c:308
#define CASE_B
Definition: nurb_trim.c:62
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define CASE_C
Definition: nurb_trim.c:63
int nurb_crv_is_bezier(const struct edge_g_cnurb *crv)
Definition: nurb_bezier.c:108
#define BU_LIST_FIRST_MAGIC(hp)
Definition: list.h:416
#define BU_LIST_WHILE(p, structure, hp)
Definition: list.h:410
#define NMG_VERTEXUSE_A_CNURB_MAGIC
Definition: magic.h:143
int nmg_uv_in_lu(const fastf_t u, const fastf_t v, const struct loopuse *lu)
Definition: nurb_trim.c:497
#define ZERO(val)
Definition: units.c:38
#define BU_LIST_INIT(_hp)
Definition: list.h:148
#define NMG_EDGE_G_CNURB_MAGIC
Definition: magic.h:121
axis
Definition: color.c:48
int clip(fastf_t *, fastf_t *, fastf_t *, fastf_t *)
Definition: clip.c:66
void rt_nurb_free_cnurb(struct edge_g_cnurb *crv)
Definition: nurb_util.c:170
struct edge_g_cnurb * rt_nurb_c_xsplit(struct edge_g_cnurb *crv, fastf_t param)
Definition: nurb_xsplit.c:226
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
void bu_bomb(const char *str) _BU_ATTR_NORETURN
Definition: bomb.c:91
double fastf_t
Definition: defines.h:300
void rt_clip_cnurb(struct bu_list *plist, struct edge_g_cnurb *crv, fastf_t u, fastf_t v)
Definition: nurb_trim.c:384
fastf_t rt_trim_line_pt_dist(struct _interior_line *l, fastf_t *pt, int pt_type)
Definition: nurb_trim.c:335
fastf_t o_dist
Definition: nurb_trim.c:42