BRL-CAD
bezier.c
Go to the documentation of this file.
1 /* B E Z I E R . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2004-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 ray */
21 /** @{ */
22 /** @file librt/bezier.c
23  *
24  * The following routines are for 2D Bezier curves.
25  *
26  * The following routines are borrowed from Graphics Gems I, Academic
27  * Press, Inc., 1990, Andrew S. Glassner (editor), "A Bezier
28  * Curve-based Root-finder", Philip J. Schneider.
29  *
30  * Modifications have been made for inclusion in BRL-CAD and to
31  * generalize the codes for finding intersections with any 2D line
32  * rather than just the X-axis.
33  */
34 
35 #include "common.h"
36 
37 #include "bio.h"
38 
39 #include "vmath.h"
40 #include "raytrace.h"
41 #include "nurb.h"
42 
43 #include "./librt_private.h"
44 
45 
46 #define SGN(_x) (((_x)<0) ? -1 : 1)
47 #define MAXDEPTH 64
48 
49 /*
50  * Count the number of times a Bezier control polygon
51  * crosses the ray. This number is >= the number of roots.
52  */
53 HIDDEN int
55  point2d_t *V, /* 2D Control pts of Bezier curve */
56  int degree, /* Degree of Bezier curve */
57  point2d_t ray_start, /* starting point for ray */
58  point2d_t ray_perp) /* unit vector perpendicular to ray direction */
59 {
60  int i;
61  int n_crossings = 0; /* Number of crossings */
62  int sign, old_sign; /* Sign of coefficients */
63  point2d_t to_pt; /* vector from ray start to a control point */
64 
65  V2SUB2(to_pt, ray_start, V[0]);
66  sign = old_sign = SGN(V2DOT(to_pt, ray_perp));
67  for (i = 1; i <= degree; i++) {
68  V2SUB2(to_pt, ray_start, V[i]);
69  sign = SGN(V2DOT(to_pt, ray_perp));
70  if (sign != old_sign) n_crossings++;
71  old_sign = sign;
72  }
73 
74  return n_crossings;
75 }
76 
77 
78 /*
79  * Check if the control polygon of a Bezier curve is flat enough for
80  * recursive subdivision to bottom out.
81  */
82 HIDDEN int
84  point2d_t *V, /* Control points */
85  int degree, /* Degree of polynomial */
86  fastf_t epsilon) /* Maximum allowable error */
87 {
88  int i; /* Index variable */
89  double *distance; /* Distances from pts to line */
90  double max_distance_above; /* maximum of these */
91  double max_distance_below;
92  double error; /* Precision of root */
93  double intercept_1,
94  intercept_2,
95  left_intercept,
96  right_intercept;
97  double a, b, c; /* Coefficients of implicit */
98  /* eqn for line from V[0]-V[deg]*/
99 
100  /* Find the perpendicular distance */
101  /* from each interior control point to */
102  /* line connecting V[0] and V[degree] */
103  distance = (double *)bu_malloc((unsigned)(degree + 1) * sizeof(double), "control_polygon_flat_enough");
104  {
105  double abSquared;
106 
107  /* Derive the implicit equation for line connecting first */
108  /* and last control points */
109  a = V[0][Y] - V[degree][Y];
110  b = V[degree][X] - V[0][X];
111  c = V[0][X] * V[degree][Y] - V[degree][X] * V[0][Y];
112 
113  abSquared = 1.0 / ((a * a) + (b * b));
114 
115  for (i = 1; i < degree; i++) {
116  /* Compute distance from each of the points to that line */
117  distance[i] = a * V[i][X] + b * V[i][Y] + c;
118  if (distance[i] > 0.0) {
119  distance[i] = (distance[i] * distance[i]) * abSquared;
120  }
121  if (distance[i] < 0.0) {
122  distance[i] = -((distance[i] * distance[i]) * abSquared);
123  }
124  }
125  }
126 
127 
128  /* Find the largest distance */
129  max_distance_above = 0.0;
130  max_distance_below = 0.0;
131  for (i = 1; i < degree; i++) {
132  if (distance[i] < 0.0) {
133  max_distance_below = FMIN(max_distance_below, distance[i]);
134  };
135  if (distance[i] > 0.0) {
136  max_distance_above = FMAX(max_distance_above, distance[i]);
137  }
138  }
139  bu_free((char *)distance, "control_polygon_flat_enough");
140 
141  {
142  double det, dInv;
143  double a1, b1, c1, a2, b2, c2;
144 
145  if (NEAR_ZERO(a, VUNITIZE_TOL)) {
146  a1 = 1.0;
147  b1 = 1.0;
148  c1 = 0.0;
149  } else {
150  /* Implicit equation for zero line */
151  a1 = 0.0;
152  b1 = 1.0;
153  c1 = 0.0;
154  }
155 
156  /* Implicit equation for "above" line */
157  a2 = a;
158  b2 = b;
159  c2 = c + max_distance_above;
160 
161  det = a1 * b2 - a2 * b1;
162  dInv = 1.0/det;
163 
164  intercept_1 = (b1 * c2 - b2 * c1) * dInv;
165 
166  /* Implicit equation for "below" line */
167  a2 = a;
168  b2 = b;
169  c2 = c + max_distance_below;
170 
171  det = a1 * b2 - a2 * b1;
172  dInv = 1.0/det;
173 
174  intercept_2 = (b1 * c2 - b2 * c1) * dInv;
175  }
176 
177  /* Compute intercepts of bounding box */
178  left_intercept = FMIN(intercept_1, intercept_2);
179  right_intercept = FMAX(intercept_1, intercept_2);
180 
181  error = 0.5 * (right_intercept-left_intercept);
182 
183  if (error < epsilon) {
184  return 1;
185  } else {
186  return 0;
187  }
188 }
189 
190 
191 void
193  point2d_t *V, /* Control pts */
194  int degree, /* Degree of bezier curve */
195  double t, /* Parameter value [0..1] */
196  point2d_t *Left, /* RETURN left half ctl pts */
197  point2d_t *Right, /* RETURN right half ctl pts */
198  point2d_t eval_pt, /* RETURN evaluated point */
199  point2d_t normal) /* RETURN unit normal at evaluated pt (may be NULL) */
200 {
201  int i, j; /* Index variables */
202  fastf_t len;
203  point2d_t tangent;
204  point2d_t **Vtemp;
205 
206 
207  Vtemp = (point2d_t **)bu_calloc(degree+1, sizeof(point2d_t *), "bezier: Vtemp array");
208  for (i=0; i<=degree; i++)
209  Vtemp[i] = (point2d_t *)bu_calloc(degree+1, sizeof(point2d_t),
210  "bezier: Vtemp[i] array");
211  /* Copy control points */
212  for (j =0; j <= degree; j++) {
213  V2MOVE(Vtemp[0][j], V[j]);
214  }
215 
216  /* Triangle computation */
217  for (i = 1; i <= degree; i++) {
218  for (j =0; j <= degree - i; j++) {
219  Vtemp[i][j][X] =
220  (1.0 - t) * Vtemp[i-1][j][X] + t * Vtemp[i-1][j+1][X];
221  Vtemp[i][j][Y] =
222  (1.0 - t) * Vtemp[i-1][j][Y] + t * Vtemp[i-1][j+1][Y];
223  }
224  }
225 
226  if (Left != NULL) {
227  for (j = 0; j <= degree; j++) {
228  V2MOVE(Left[j], Vtemp[j][0]);
229  }
230  }
231  if (Right != NULL) {
232  for (j = 0; j <= degree; j++) {
233  V2MOVE(Right[j], Vtemp[degree-j][j]);
234  }
235  }
236 
237  V2MOVE(eval_pt, Vtemp[degree][0]);
238 
239  if (normal) {
240  V2SUB2(tangent, Vtemp[degree-1][1], Vtemp[degree-1][0]);
241  normal[X] = tangent[Y];
242  normal[Y] = -tangent[X];
243  len = sqrt(MAG2SQ(normal));
244  normal[X] /= len;
245  normal[Y] /= len;
246  }
247 
248  for (i=0; i<=degree; i++)
249  bu_free((char *)Vtemp[i], "bezier: Vtemp[i]");
250  bu_free((char *)Vtemp, "bezier: Vtemp");
251 
252  return;
253 }
254 
255 
256 /*
257  * Compute intersection of chord from first control point to last
258  * with ray.
259  *
260  * Returns :
261  *
262  * 0 - no intersection
263  * 1 - found an intersection
264  *
265  * intercept - contains calculated intercept
266  */
267 HIDDEN int
269  point2d_t *V, /* Control points */
270  int degree, /* Degree of curve */
271  point2d_t ray_start, /* starting point of ray */
272  point2d_t ray_dir, /* unit ray direction */
273  point2d_t intercept, /* calculated intercept point */
274  point2d_t normal) /* calculated unit normal at intercept */
275 {
276  fastf_t beta;
277  fastf_t denom;
278  fastf_t len;
279  point2d_t seg_line;
280 
281  denom = (V[degree][X] - V[0][X]) * ray_dir[Y] -
282  (V[degree][Y] - V[0][Y]) * ray_dir[X];
283 
284  if (ZERO(denom))
285  return 0;
286 
287  beta = (V[0][Y] * ray_dir[X] - V[0][X] * ray_dir[Y] +
288  ray_start[X] * ray_dir[Y] - ray_start[Y] * ray_dir[X]) / denom;
289 
290  if (beta < 0.0 || beta > 1.0)
291  return 0;
292 
293  V2SUB2(seg_line, V[degree], V[0]);
294  V2JOIN1(intercept, V[0], beta, seg_line);
295 
296  /* calculate normal */
297  normal[X] = seg_line[Y];
298  normal[Y] = -seg_line[X];
299  len = sqrt(MAG2SQ(seg_line));
300  normal[X] /= len;
301  normal[Y] /= len;
302 
303  return 1;
304 }
305 
306 
307 int
309  point2d_t *w, /* The control points */
310  int degree, /* The degree of the polynomial */
311  point2d_t **intercept, /* list of intersections found */
312  point2d_t **normal, /* corresponding normals */
313  point2d_t ray_start, /* starting point of ray */
314  point2d_t ray_dir, /* Unit direction for ray */
315  point2d_t ray_perp, /* Unit vector normal to ray_dir */
316  int depth, /* The depth of the recursion */
317  fastf_t epsilon) /* maximum allowable error */
318 {
319  int i;
320  point2d_t *Left, /* New left and right */
321  *Right; /* control polygons */
322  int left_count, /* Solution count from */
323  right_count; /* children */
324  point2d_t *left_t, /* Solutions from kids */
325  *right_t;
326  point2d_t *left_n, /* normals from kids */
327  *right_n;
328  int total_count;
329  point2d_t eval_pt;
330 
331  switch (crossing_count(w, degree, ray_start, ray_perp)) {
332  case 0 : {
333  /* No solutions here */
334  return 0;
335  }
336  case 1 : {
337  /* Unique solution */
338  /* Stop recursion when the tree is deep enough */
339  /* if deep enough, return 1 solution at midpoint */
340  if (depth >= MAXDEPTH) {
341  BU_ALLOC(*intercept, point2d_t);
342  BU_ALLOC(*normal, point2d_t);
343  bezier(w, degree, 0.5, NULL, NULL, *intercept[0], *normal[0]);
344  return 1;
345  }
346  if (control_polygon_flat_enough(w, degree, epsilon)) {
347  BU_ALLOC(*intercept, point2d_t);
348  BU_ALLOC(*normal, point2d_t);
349  if (!compute_x_intercept(w, degree, ray_start, ray_dir, *intercept[0], *normal[0])) {
350  bu_free((char *)(*intercept), "bezier_roots: no solution");
351  bu_free((char *)(*normal), "bezier_roots: no solution");
352  return 0;
353  }
354  return 1;
355  }
356  break;
357  }
358  }
359 
360  /* Otherwise, solve recursively after */
361  /* subdividing control polygon */
362  Left = (point2d_t *)bu_calloc(degree+1, sizeof(point2d_t), "bezier_roots: Left");
363  Right = (point2d_t *)bu_calloc(degree+1, sizeof(point2d_t), "bezier_roots: Right");
364  bezier(w, degree, 0.5, Left, Right, eval_pt, NULL);
365 
366  left_count = bezier_roots(Left, degree, &left_t, &left_n, ray_start, ray_dir, ray_perp, depth+1, epsilon);
367  right_count = bezier_roots(Right, degree, &right_t, &right_n, ray_start, ray_dir, ray_perp, depth+1, epsilon);
368 
369  total_count = left_count + right_count;
370 
371  bu_free((char *)Left, "bezier_roots: Left");
372  bu_free((char *)Right, "bezier_roots: Right");
373  if (total_count) {
374  *intercept = (point2d_t *)bu_calloc(total_count, sizeof(point2d_t),
375  "bezier_roots: roots compilation");
376  *normal = (point2d_t *)bu_calloc(total_count, sizeof(point2d_t),
377  "bezier_roots: normal compilation");
378  }
379 
380  /* Gather solutions together */
381  for (i = 0; i < left_count; i++) {
382  V2MOVE((*intercept)[i], left_t[i]);
383  V2MOVE((*normal)[i], left_n[i]);
384  }
385  for (i = 0; i < right_count; i++) {
386  V2MOVE((*intercept)[i+left_count], right_t[i]);
387  V2MOVE((*normal)[i+left_count], right_n[i]);
388  }
389 
390  if (left_count) {
391  bu_free((char *)left_t, "Left roots");
392  bu_free((char *)left_n, "Left normals");
393  }
394  if (right_count) {
395  bu_free((char *)right_t, "Right roots");
396  bu_free((char *)right_n, "Right normals");
397  }
398 
399  /* Send back total number of solutions */
400  return total_count;
401 }
402 
403 
404 struct bezier_2d_list *
405 bezier_subdivide(struct bezier_2d_list *bezier_in, int degree, fastf_t epsilon, int depth)
406 {
407  struct bezier_2d_list *bz_l, *bz_r, *new_head;
408  struct bezier_2d_list *left_rtrn, *rt_rtrn;
409  point2d_t pt;
410 
411  /* create a new head */
412  BU_ALLOC(new_head, struct bezier_2d_list);
413 
414  BU_LIST_INIT(&new_head->l);
415  if (depth >= MAXDEPTH) {
416  BU_LIST_APPEND(&new_head->l, &bezier_in->l);
417  return new_head;
418  }
419 
420  if (control_polygon_flat_enough(bezier_in->ctl, degree, epsilon)) {
421  BU_LIST_APPEND(&new_head->l, &bezier_in->l);
422  return new_head;
423  }
424 
425  /* allocate memory for left and right curves */
426  BU_ALLOC(bz_l, struct bezier_2d_list);
427  BU_LIST_INIT(&bz_l->l);
428  BU_ALLOC(bz_r, struct bezier_2d_list);
429  BU_LIST_INIT(&bz_r->l);
430  bz_l->ctl = (point2d_t *)bu_calloc(degree + 1, sizeof(point2d_t),
431  "bezier_subdivide: bz_l->ctl");
432  bz_r->ctl = (point2d_t *)bu_calloc(degree + 1, sizeof(point2d_t),
433  "bezier_subdivide: bz_r->ctl");
434 
435  /* subdivide at t = 0.5 */
436  bezier(bezier_in->ctl, degree, 0.5, bz_l->ctl, bz_r->ctl, pt, NULL);
437 
438  /* eliminate original */
439  BU_LIST_DEQUEUE(&bezier_in->l);
440  bu_free((char *)bezier_in->ctl, "bezier_subdivide: bezier_in->ctl");
441  bu_free((char *)bezier_in, "bezier_subdivide: bezier_in");
442 
443  /* recurse on left curve */
444  left_rtrn = bezier_subdivide(bz_l, degree, epsilon, depth+1);
445 
446  /* recurse on right curve */
447  rt_rtrn = bezier_subdivide(bz_r, degree, epsilon, depth+1);
448 
449  BU_LIST_APPEND_LIST(&new_head->l, &left_rtrn->l);
450  BU_LIST_APPEND_LIST(&new_head->l, &rt_rtrn->l);
451  bu_free((char *)left_rtrn, "bezier_subdivide: left_rtrn (head)");
452  bu_free((char *)rt_rtrn, "bezier_subdivide: rt_rtrn (head)");
453 
454  return new_head;
455 }
456 
457 
458 /** @} */
459 /*
460  * Local Variables:
461  * mode: C
462  * tab-width: 8
463  * indent-tabs-mode: t
464  * c-file-style: "stroustrup"
465  * End:
466  * ex: shiftwidth=4 tabstop=8
467  */
HIDDEN int sign(double val)
Definition: brep.cpp:1145
int bezier_roots(point2d_t *w, int degree, point2d_t **intercept, point2d_t **normal, point2d_t ray_start, point2d_t ray_dir, point2d_t ray_perp, int depth, fastf_t epsilon)
Definition: bezier.c:308
HIDDEN int compute_x_intercept(point2d_t *V, int degree, point2d_t ray_start, point2d_t ray_dir, point2d_t intercept, point2d_t normal)
Definition: bezier.c:268
Header file for the BRL-CAD common definitions.
#define BU_LIST_APPEND(old, new)
Definition: list.h:197
#define HIDDEN
Definition: common.h:86
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
Definition: color.c:49
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
void * bu_calloc(size_t nelem, size_t elsize, const char *str)
Definition: malloc.c:321
#define NEAR_ZERO(val, epsilon)
Definition: color.c:55
#define SGN(_x)
Definition: bezier.c:46
HIDDEN int control_polygon_flat_enough(point2d_t *V, int degree, fastf_t epsilon)
Definition: bezier.c:83
#define ZERO(val)
Definition: units.c:38
#define BU_LIST_INIT(_hp)
Definition: list.h:148
#define FMAX(a, b)
Definition: common.h:102
void bezier(point2d_t *V, int degree, double t, point2d_t *Left, point2d_t *Right, point2d_t eval_pt, point2d_t normal)
Definition: bezier.c:192
#define MAXDEPTH
Definition: bezier.c:47
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
HIDDEN int crossing_count(point2d_t *V, int degree, point2d_t ray_start, point2d_t ray_perp)
Definition: bezier.c:54
#define BU_LIST_DEQUEUE(cur)
Definition: list.h:209
double fastf_t
Definition: defines.h:300
struct bezier_2d_list * bezier_subdivide(struct bezier_2d_list *bezier_in, int degree, fastf_t epsilon, int depth)
Definition: bezier.c:405
Definition: color.c:50
#define BU_LIST_APPEND_LIST(dest_hp, src_hp)
Definition: list.h:281
#define FMIN(a, b)
Definition: common.h:105