BRL-CAD
sketch_tess.cpp
Go to the documentation of this file.
1 /* S K E T C H _ T E S S . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2012-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 primitives */
21 /** @{ */
22 /** @file primitives/sketch/sketch_tess.c
23  *
24  * An extension of sketch.c for functions using the openNURBS library. Includes
25  * support for approximation of Bezier curves using circular arcs.
26  */
27 /** @} */
28 
29 /* NOTE: Current Bezier approximation routine is meant to handle cubic and lower
30  * degree Bezier curves, and may not work properly for higher degrees.
31  *
32  * TODO:
33  * rt_sketch_surf_area: add logic to determine if part of a sketch
34  * "removes" area from the total area, e.g. a sketch of a circumscribed square
35  * has area less than the sum of the areas of the square and circle.
36  *
37  * rt_sketch_tess: implement this.
38  */
39 
40 #include "common.h"
41 
42 #include <vector>
43 
44 #include "raytrace.h"
45 #include "rtgeom.h"
46 #include "vmath.h"
47 #include "brep.h"
48 
49 
50 /* returns the incenter of the inscribed circle inside the triangle defined by
51  * points a, b, c
52  */
53 HIDDEN ON_2dPoint
54 incenter(const ON_2dPoint a, const ON_2dPoint b, const ON_2dPoint c)
55 {
56  fastf_t a_b, a_c, b_c, sum;
57  ON_2dPoint incenter;
58 
59  a_b = a.DistanceTo(b);
60  a_c = a.DistanceTo(c);
61  b_c = b.DistanceTo(c);
62  sum = a_b + a_c + b_c;
63 
64  incenter.x = (b_c * a.x + a_c * b.x + a_b * c.x) / sum;
65  incenter.y = (b_c * a.y + a_c * b.y + a_b * c.y) / sum;
66  return incenter;
67 }
68 
69 
70 /* create a biarc for a bezier curve.
71  *
72  * extends the tangent lines to the bezier curve at its first and last control
73  * points, and intersects them to find a third point.
74  * the biarc passes through the first and last control points, and the incenter
75  * of the circle defined by the first, last and intersection points.
76  */
77 HIDDEN ON_Arc
78 make_biarc(const ON_BezierCurve& bezier)
79 {
80  ON_2dPoint isect, arc_pt;
81  ON_2dPoint p_start(bezier.PointAt(0)), p_end(bezier.PointAt(1.0));
82  ON_2dVector t_start(bezier.TangentAt(0)), t_end(bezier.TangentAt(1.0));
83  ON_Ray r_start(p_start, t_start), r_end(p_end, t_end);
84 
85  r_start.IntersectRay(r_end, isect);
86  arc_pt = incenter(p_start, p_end, isect);
87 
88  return ON_Arc(p_start, arc_pt, p_end);
89 }
90 
91 
92 /* NOTE: MINSTEP and MAXSTEP were determined by experimentation. If MINSTEP is
93  * much smaller (i.e. 0.00001), approx_bezier() slows significantly on curves with
94  * high curvature over a large part of its domain. MAXSTEP represents a step
95  * size of 1/10th the domain of a Bezier curve, and MINSTEP represents 1/10000th.
96  */
97 #define MINSTEP 0.0001
98 #define MAXSTEP 0.1
99 
100 #define GETSTEPSIZE(step) ((step) > MAXSTEP ? MAXSTEP : ((step) < MINSTEP ? MINSTEP : (step)))
101 
102 #define POW3(x) ((x)*(x)*(x))
103 #define SIGN(x) ((x) >= 0 ? 1 : -1)
104 #define CURVATURE(d1, d2) (V2CROSS((d1), (d2)) / POW3((d1).Length()))
105 
106 /* find a point of inflection on a bezier curve, if it exists, by finding the
107  * value of parameter 't' where the signed curvature of the bezier changes
108  * signs. Returns true if an inflection point is found.
109  */
110 HIDDEN bool
111 bezier_inflection(const ON_BezierCurve& bezier, fastf_t& inflection_pt)
112 {
113  int sign;
114  fastf_t t, step, crv;
115  ON_3dVector d1, d2; // first derivative, second derivative
116  ON_3dPoint tmp; // not used, but needed by Ev2Der
117 
118  // calculate curvature at t=0
119  bezier.Ev2Der(0, tmp, d1, d2);
120  crv = CURVATURE(d1, d2);
121  // step size decreases as |crv| -> 0
122  step = GETSTEPSIZE(fabs(crv));
123 
124  sign = SIGN(crv);
125 
126  for (t = step; t <= 1.0; t += step) {
127  bezier.Ev2Der(t, tmp, d1, d2);
128  crv = CURVATURE(d1, d2);
129  // if sign changes, t is an inflection point
130  if (sign != SIGN(crv)) {
131  inflection_pt = t;
132  return true;
133  }
134  step = GETSTEPSIZE(fabs(crv));
135  }
136  return false;
137 }
138 
139 
140 /* approximates a bezier curve with a set of circular arcs by dividing where
141  * the bezier's deviation from its approximating biarc is at a maximum, then
142  * recursively calling on the subsections until it is approximated to
143  * tolerance by the biarc
144  */
145 HIDDEN void
146 approx_bezier(const ON_BezierCurve& bezier, const ON_Arc& biarc, const struct bn_tol *tol, std::vector<ON_Arc>& approx)
147 {
148  fastf_t t = 0.0, step = 0.0;
149  fastf_t crv = 0.0, err = 0.0, max_t = 0.0, max_err = 0.0;
150  ON_3dPoint test;
151  ON_3dVector d1, d2;
152 
153  // walk the bezier curve at interval given by step
154  for (t = 0; t <= 1.0; t += step) {
155  bezier.Ev2Der(t, test, d1, d2);
156  err = fabs((test - biarc.Center()).Length() - biarc.Radius());
157  // find the maximum point of deviation
158  if (err > max_err) {
159  max_t = t;
160  max_err = err;
161  }
162  crv = CURVATURE(d1, d2);
163  // step size decreases as |crv| -> 1
164  step = GETSTEPSIZE(1.0 - fabs(crv));
165  }
166 
167  if (max_err + VDIVIDE_TOL < tol->dist) {
168  // max deviation is less than the given tolerance, add the biarc approximation
169  approx.push_back(biarc);
170  } else {
171  ON_BezierCurve head, tail;
172  // split bezier at point of maximum deviation and recurse on the new subsections
173  bezier.Split(max_t, head, tail);
174  approx_bezier(head, make_biarc(head), tol, approx);
175  approx_bezier(tail, make_biarc(tail), tol, approx);
176  }
177 }
178 
179 
180 /* approximates a bezier curve with a set of circular arcs.
181  * returns approximation in carcs
182  */
183 HIDDEN void
184 bezier_to_carcs(const ON_BezierCurve& bezier, const struct bn_tol *tol, std::vector<ON_Arc>& carcs)
185 {
186  bool skip_while = true, curvature_changed = false;
187  fastf_t inflection_pt, biarc_angle;
188  ON_Arc biarc;
189  ON_BezierCurve current, next;
190  std::vector<ON_BezierCurve> rest;
191 
192  // find inflection point, if it exists
193  if (bezier_inflection(bezier, inflection_pt)) {
194  curvature_changed = true;
195  bezier.Split(inflection_pt, current, next);
196  rest.push_back(next);
197  } else {
198  current = bezier;
199  }
200 
201  while (skip_while || !rest.empty()) {
202  if (skip_while) skip_while = false;
203  biarc = make_biarc(current);
204  if ((biarc_angle = biarc.AngleRadians()) <= M_PI_2) {
205  // approximate the current bezier segment and add its biarc
206  // approximation to carcs
207  approx_bezier(current, biarc, tol, carcs);
208  } else if (biarc_angle <= M_PI) {
209  // divide the current bezier segment in half
210  current.Split(0.5, current, next);
211  // approximate first bezier segment
212  approx_bezier(current, biarc, tol, carcs);
213  // approximate second bezier segment
214  approx_bezier(next, biarc, tol, carcs);
215  } else {
216  fastf_t t = 1.0;
217  ON_Arc test_biarc;
218  ON_BezierCurve test_bezier;
219  // divide the current bezier such that the first curve segment would
220  // have an approximating biarc segment <=90 degrees
221  do {
222  t *= 0.5;
223  current.Split(t, test_bezier, next);
224  test_biarc = make_biarc(test_bezier);
225  } while (test_biarc.AngleRadians() > M_PI_2);
226 
227  approx_bezier(test_bezier, test_biarc, tol, carcs);
228  current = next;
229  skip_while = true;
230  continue;
231  }
232 
233  if (curvature_changed) {
234  curvature_changed = false;
235  current = rest.back();
236  rest.pop_back();
237  // continue even if we just popped the last element
238  skip_while = true;
239  }
240  }
241 }
242 
243 
244 #define SKETCH_PT(idx) sketch_ip->verts[(idx)]
245 
246 #define DIST_PT2D_PT2D_SQ(_p1, _p2) \
247  (((_p2)[X] - (_p1)[X]) * ((_p2)[X] - (_p1)[X]) + \
248  ((_p2)[Y] - (_p1)[Y]) * ((_p2)[Y] - (_p1)[Y]))
249 
250 #define DIST_PT2D_PT2D(_p1, _p2) sqrt(DIST_PT2D_PT2D_SQ(_p1, _p2))
251 
252 /**
253  * calculate approximate surface area for a sketch primitive by iterating through
254  * each curve segment in the sketch, calculating the area of the polygon
255  * created by the start and end points of each curve segment, as well as the
256  * additional areas for circular segments.
257  *
258  * line_seg: calculate the area for the polygon edge Start->End
259  * carc_seg: if the segment is a full circle, calculate its area.
260  * else, calculate the area for the polygon edge Start->End, and the area
261  * of the circular segment
262  * bezier_seg: approximate the bezier using the bezier_to_carcs() function. for
263  * each carc_seg, calculate the area for the polygon edges Start->End and
264  * the area of the circular segment
265  */
266 extern "C" void
268 {
269  int j;
270  size_t i;
271  fastf_t poly_area = 0, carc_area = 0;
272 
273  struct bn_tol tol;
274  struct rt_sketch_internal *sketch_ip = (struct rt_sketch_internal *)ip->idb_ptr;
275  RT_SKETCH_CK_MAGIC(sketch_ip);
276  struct rt_curve crv = sketch_ip->curve;
277 
278  // a sketch with no curves has no area
279  if (UNLIKELY(crv.count == 0)) {
280  return;
281  }
282 
283  tol.magic = BN_TOL_MAGIC;
284  tol.dist = RT_DOT_TOL;
285  tol.dist_sq = RT_DOT_TOL * RT_DOT_TOL;
286 
287  // iterate through each curve segment in the sketch
288  for (i = 0; i < crv.count; i++) {
289  const struct line_seg *lsg;
290  const struct carc_seg *csg;
291  const struct bezier_seg *bsg;
292 
293  const uint32_t *lng = (uint32_t *)crv.segment[i];
294 
295  switch (*lng) {
296  case CURVE_LSEG_MAGIC:
297  lsg = (struct line_seg *)lng;
298  // calculate area for polygon edge
299  poly_area += V2CROSS(SKETCH_PT(lsg->start), SKETCH_PT(lsg->end));
300  break;
301  case CURVE_CARC_MAGIC:
302  csg = (struct carc_seg *)lng;
303  if (csg->radius < 0) {
304  // calculate full circle area
305  carc_area += M_PI * DIST_PT2D_PT2D_SQ(SKETCH_PT(csg->start), SKETCH_PT(csg->end));
306  } else {
307  fastf_t theta, side_ratio;
308  // calculate area for polygon edge
309  poly_area += V2CROSS(SKETCH_PT(csg->start), SKETCH_PT(csg->end));
310  // calculate area for circular segment
311  side_ratio = DIST_PT2D_PT2D(SKETCH_PT(csg->start), SKETCH_PT(csg->end)) / (2.0 * csg->radius);
312  theta = asin(side_ratio);
313  carc_area += 0.5 * csg->radius * csg->radius * (theta - side_ratio);
314  }
315  break;
316  case CURVE_BEZIER_MAGIC: {
317  bsg = (struct bezier_seg *)lng;
318  ON_2dPointArray bez_pts(bsg->degree + 1);
319  std::vector<ON_Arc> carcs;
320  // convert struct bezier_seg into ON_BezierCurve
321  for (j = 0; j < bsg->degree + 1; j++) {
322  bez_pts.Append(SKETCH_PT(bsg->ctl_points[j]));
323  }
324  // approximate bezier curve by a set of circular arcs
325  bezier_to_carcs(ON_BezierCurve(bez_pts), &tol, carcs);
326  for (std::vector<ON_Arc>::iterator it = carcs.begin(); it != carcs.end(); ++it) {
327  // calculate area for each polygon edge
328  poly_area += V2CROSS(it->StartPoint(), it->EndPoint());
329  // calculate area for each circular segment
330  carc_area += 0.5 * it->Radius() * it->Radius() * (it->AngleRadians() - sin(it->AngleRadians()));
331  }
332  }
333  break;
334  case CURVE_NURB_MAGIC:
335  default:
336  break;
337  }
338  }
339  *area = 0.5 * fabs(poly_area) + carc_area;
340 }
341 
342 
343 /*
344  * Local Variables:
345  * mode: C
346  * tab-width: 8
347  * indent-tabs-mode: t
348  * c-file-style: "stroustrup"
349  * End:
350  * ex: shiftwidth=4 tabstop=8
351  */
HIDDEN int sign(double val)
Definition: brep.cpp:1145
#define DIST_PT2D_PT2D_SQ(_p1, _p2)
#define RT_DOT_TOL
Definition: raytrace.h:170
HIDDEN ON_2dPoint incenter(const ON_2dPoint a, const ON_2dPoint b, const ON_2dPoint c)
Definition: sketch_tess.cpp:54
double dist
>= 0
Definition: tol.h:73
HIDDEN bool bezier_inflection(const ON_BezierCurve &bezier, fastf_t &inflection_pt)
HIDDEN void approx_bezier(const ON_BezierCurve &bezier, const ON_Arc &biarc, const struct bn_tol *tol, std::vector< ON_Arc > &approx)
void rt_sketch_surf_area(fastf_t *area, const struct rt_db_internal *ip)
#define M_PI
Definition: fft.h:35
double dist_sq
dist * dist
Definition: tol.h:74
#define CURVE_CARC_MAGIC
Definition: magic.h:198
#define BN_TOL_MAGIC
Definition: magic.h:74
Header file for the BRL-CAD common definitions.
#define SKETCH_PT(idx)
#define CURVE_NURB_MAGIC
Definition: magic.h:200
#define HIDDEN
Definition: common.h:86
#define CURVE_BEZIER_MAGIC
Definition: magic.h:197
HIDDEN void bezier_to_carcs(const ON_BezierCurve &bezier, const struct bn_tol *tol, std::vector< ON_Arc > &carcs)
Support for uniform tolerances.
Definition: tol.h:71
#define SIGN(x)
HIDDEN ON_Arc make_biarc(const ON_BezierCurve &bezier)
Definition: sketch_tess.cpp:78
void * idb_ptr
Definition: raytrace.h:195
uint32_t magic
Definition: tol.h:72
#define GETSTEPSIZE(step)
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 CURVATURE(d1, d2)
#define DIST_PT2D_PT2D(_p1, _p2)
#define CURVE_LSEG_MAGIC
Definition: magic.h:199
double fastf_t
Definition: defines.h:300
#define UNLIKELY(expression)
Definition: common.h:282