BRL-CAD
primitive_util.c
Go to the documentation of this file.
1 /* P R I M I T I V E _ U T I L . 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 /** @file primitive_util.c
21  *
22  * General utility routines for librt primitives that are not part of
23  * the librt API. Prototypes for these routines are located in
24  * librt_private.h.
25  */
26 
27 #include "../librt_private.h"
28 
29 /**
30  * If only the absolute tolerance is valid (positive), it is returned.
31  *
32  * If only the relative tolerance is valid (in (0.0, 1.0)), the return value
33  * is the relative tolerance multiplied by rel_to_abs.
34  *
35  * If both tolerances are valid, the smaller is used.
36  * If neither is valid, a default of .1 * rel_to_abs is used.
37  */
38 #define DEFAULT_REL_TOL .1
39 fastf_t
41  const struct rt_tess_tol *ttol,
42  fastf_t rel_to_abs)
43 {
44  fastf_t tol, rel_tol, abs_tol;
45  int rel_tol_is_valid;
46 
47  rel_tol = ttol->rel;
48  abs_tol = ttol->abs;
49 
50  rel_tol_is_valid = 0;
51  if (rel_tol > 0.0 && rel_tol < 1.0) {
52  rel_tol_is_valid = 1;
53  }
54 
55  if (abs_tol > 0.0 && (!rel_tol_is_valid || rel_tol > abs_tol)) {
56  tol = abs_tol;
57  } else {
58  if (!rel_tol_is_valid) {
59  rel_tol = DEFAULT_REL_TOL;
60  }
61  tol = rel_tol * rel_to_abs;
62  }
63 
64  return tol;
65 }
66 
67 /**
68  * Gives a rough estimate of the maximum number of times a primitive's bounding
69  * box diagonal will be sampled based on the sample density of the view.
70  *
71  * Practically, it is an estimate of the maximum number of pixels that would be
72  * used if the diagonal line were drawn in the current view window.
73  *
74  * It is currently used in adaptive plot routines to help choose how many
75  * sample points should be used for plotted curves.
76  */
77 fastf_t
79  struct rt_db_internal *ip,
80  const struct rt_view_info *info)
81 {
82  point_t bbox_min, bbox_max;
83  fastf_t primitive_diagonal_mm, samples_per_mm;
84  fastf_t diagonal_samples;
85 
86  ip->idb_meth->ft_bbox(ip, &bbox_min, &bbox_max, info->tol);
87  primitive_diagonal_mm = DIST_PT_PT(bbox_min, bbox_max);
88 
89  samples_per_mm = 1.0 / info->point_spacing;
90  diagonal_samples = samples_per_mm * primitive_diagonal_mm;
91 
92  return diagonal_samples;
93 }
94 
95 /**
96  * Estimate the number of evenly spaced cross-section curves needed to meet a
97  * target curve spacing.
98  */
99 fastf_t
101  struct rt_db_internal *ip,
102  const struct rt_view_info *info)
103 {
104  point_t bbox_min, bbox_max;
105  fastf_t x_len, y_len, z_len, avg_len;
106 
107  if (!ip->idb_meth->ft_bbox) {
108  return 0.0;
109  }
110 
111  ip->idb_meth->ft_bbox(ip, &bbox_min, &bbox_max, info->tol);
112 
113  x_len = fabs(bbox_max[X] - bbox_min[X]);
114  y_len = fabs(bbox_max[Y] - bbox_min[Y]);
115  z_len = fabs(bbox_max[Z] - bbox_min[Z]);
116 
117  avg_len = (x_len + y_len + z_len) / 3.0;
118 
119  return avg_len / info->curve_spacing;
120 }
121 
122 /* Calculate the length of the shortest distance between a point and a line in
123  * the y-z plane.
124  */
125 static fastf_t
126 distance_point_to_line(point_t p, fastf_t slope, fastf_t intercept)
127 {
128  /* input line:
129  * z = slope * y + intercept
130  *
131  * implicit form is:
132  * az + by + c = 0, where a = -slope, b = 1, c = -intercept
133  *
134  * standard 2D point-line distance calculation:
135  * d = |aPy + bPz + c| / sqrt(a^2 + b^2)
136  */
137  return fabs(-slope * p[Y] + p[Z] - intercept) / sqrt(slope * slope + 1);
138 }
139 
140 /* Assume the existence of a line and a parabola at the origin, both in the y-z
141  * plane (x = 0.0). The parabola and line are described by arguments p and m,
142  * from the respective equations (z = y^2 / 4p) and (z = my + b).
143  *
144  * The line is assumed to intersect the parabola at two points.
145  *
146  * The portion of the parabola between these two intersection points takes on a
147  * single maximum value with respect to the intersecting line when its slope is
148  * the same as the line's (i.e. when the tangent line and intersecting line are
149  * parallel).
150  *
151  * The first derivate slope equation z' = y / 2p implies that the parabola has
152  * slope m when y = 2pm. So we calculate y from p and m, and then z from y.
153  */
154 static void
155 parabola_point_farthest_from_line(point_t max_point, fastf_t p, fastf_t m)
156 {
157  fastf_t y = 2.0 * p * m;
158 
159  max_point[X] = 0.0;
160  max_point[Y] = y;
161  max_point[Z] = (y * y) / ( 4.0 * p);
162 }
163 
164 /* Approximate part of the parabola (y - h)^2 = 4p(z - k) lying in the y-z
165  * plane.
166  *
167  * pts->p: the vertex (0, h, k)
168  * pts->next->p: another point on the parabola
169  * pts->next->next: NULL
170  * p: the constant from the above equation
171  *
172  * This routine inserts num_new_points points between the two input points to
173  * better approximate the parabolic curve passing between them.
174  *
175  * Returns number of points successfully added.
176  */
177 int
178 approximate_parabolic_curve(struct rt_pt_node *pts, fastf_t p, int num_new_points)
179 {
180  fastf_t error, max_error, seg_slope, seg_intercept;
181  point_t v, point, p0, p1, new_point = VINIT_ZERO;
182  struct rt_pt_node *node, *worst_node, *new_node;
183  int i;
184 
185  if (pts == NULL || pts->next == NULL || num_new_points < 1) {
186  return 0;
187  }
188 
189  VMOVE(v, pts->p);
190 
191  for (i = 0; i < num_new_points; ++i) {
192  worst_node = node = pts;
193  max_error = 0.0;
194 
195  /* Find the least accurate line segment, and calculate a new parabola
196  * point to insert between its endpoints.
197  */
198  while (node->next != NULL) {
199  VMOVE(p0, node->p);
200  VMOVE(p1, node->next->p);
201 
202  seg_slope = (p1[Z] - p0[Z]) / (p1[Y] - p0[Y]);
203  seg_intercept = p0[Z] - seg_slope * p0[Y];
204 
205  /* get farthest point on the equivalent parabola at the origin */
206  parabola_point_farthest_from_line(point, p, seg_slope);
207 
208  /* translate result to actual parabola position */
209  point[Y] += v[Y];
210  point[Z] += v[Z];
211 
212  /* see if the maximum distance between the parabola and the line
213  * (the error) is larger than the largest recorded error
214  * */
215  error = distance_point_to_line(point, seg_slope, seg_intercept);
216 
217  if (error > max_error) {
218  max_error = error;
219  worst_node = node;
220  VMOVE(new_point, point);
221  }
222 
223  node = node->next;
224  }
225 
226  /* insert new point between endpoints of the least accurate segment */
227  BU_ALLOC(new_node, struct rt_pt_node);
228  VMOVE(new_node->p, new_point);
229  new_node->next = worst_node->next;
230  worst_node->next = new_node;
231  }
232 
233  return num_new_points;
234 }
235 
236 /* Assume the existence of a line and a hyperbola with asymptote origin at the
237  * origin, both in the y-z plane (x = 0.0). The hyperbola and line are
238  * described by arguments a, b, and m, from the respective equations
239  * (z = +-(a/b) * sqrt(y^2 + b^2)) and (z = my + b).
240  *
241  * The line is assumed to intersect the positive half of the hyperbola at two
242  * points.
243  *
244  * The portion of the hyperbola between these two intersection points takes on
245  * a single maximum value with respect to the intersecting line when its slope
246  * is the same as the line's (i.e. when the tangent line and intersecting line
247  * are parallel).
248  *
249  * The first derivate slope equation z' = ay / (b * sqrt(y^2 + b^2)) implies
250  * that the hyperbola has slope m when y = mb^2 / sqrt(a^2 - b^2m^2).
251  *
252  * We calculate y from a, b, and m, and then z from y.
253  */
254 static void
255 hyperbola_point_farthest_from_line(point_t max_point, fastf_t a, fastf_t b, fastf_t m)
256 {
257  fastf_t y = (m * b * b) / sqrt((a * a) - (b * b * m * m));
258 
259  max_point[X] = 0.0;
260  max_point[Y] = y;
261  max_point[Z] = (a / b) * sqrt(y * y + b * b);
262 }
263 
264 /* Approximate part of the hyperbola lying in the Y-Z plane described by:
265  * ((z - k)^2 / a^2) - ((y - h)^2 / b^2) = 1
266  *
267  * pts->p: the asymptote origin (0, h, k)
268  * pts->next->p: another point on the hyperbola
269  * pts->next->next: NULL
270  * a, b: the constants from the above equation
271  *
272  * This routine inserts num_new_points points between the two input points to
273  * better approximate the hyperbolic curve passing between them.
274  *
275  * Returns number of points successfully added.
276  */
277 int
278 approximate_hyperbolic_curve(struct rt_pt_node *pts, fastf_t a, fastf_t b, int num_new_points)
279 {
280  fastf_t error, max_error, seg_slope, seg_intercept;
281  point_t v, point, p0, p1, new_point = VINIT_ZERO;
282  struct rt_pt_node *node, *worst_node, *new_node;
283  int i;
284 
285  if (pts == NULL || pts->next == NULL || num_new_points < 1) {
286  return 0;
287  }
288 
289  VMOVE(v, pts->p);
290 
291  for (i = 0; i < num_new_points; ++i) {
292  worst_node = node = pts;
293  max_error = 0.0;
294 
295  /* Find the least accurate line segment, and calculate a new hyperbola
296  * point to insert between its endpoints.
297  */
298  while (node->next != NULL) {
299  VMOVE(p0, node->p);
300  VMOVE(p1, node->next->p);
301 
302  seg_slope = (p1[Z] - p0[Z]) / (p1[Y] - p0[Y]);
303  seg_intercept = p0[Z] - seg_slope * p0[Y];
304 
305  /* get farthest point on the equivalent hyperbola with asymptote origin at
306  * the origin
307  */
308  hyperbola_point_farthest_from_line(point, a, b, seg_slope);
309 
310  /* translate result to actual hyperbola position */
311  point[Y] += v[Y];
312  point[Z] += v[Z] - a;
313 
314  /* see if the maximum distance between the hyperbola and the line
315  * (the error) is larger than the largest recorded error
316  * */
317  error = distance_point_to_line(point, seg_slope, seg_intercept);
318 
319  if (error > max_error) {
320  max_error = error;
321  worst_node = node;
322  VMOVE(new_point, point);
323  }
324 
325  node = node->next;
326  }
327 
328  /* insert new point between endpoints of the least accurate segment */
329  BU_ALLOC(new_node, struct rt_pt_node);
330  VMOVE(new_node->p, new_point);
331  new_node->next = worst_node->next;
332  worst_node->next = new_node;
333  }
334 
335  return num_new_points;
336 }
337 
338 void
340  point_t result,
341  const vect_t center,
342  const vect_t axis_a,
343  const vect_t axis_b,
344  fastf_t radian)
345 {
346  fastf_t cos_rad, sin_rad;
347 
348  cos_rad = cos(radian);
349  sin_rad = sin(radian);
350 
351  VJOIN2(result, center, cos_rad, axis_a, sin_rad, axis_b);
352 }
353 
354 void
356  struct bu_list *vhead,
357  const vect_t center,
358  const vect_t axis_a,
359  const vect_t axis_b,
360  int num_points)
361 {
362  int i;
363  point_t p;
364  fastf_t radian, radian_step;
365 
366  radian_step = M_2PI / num_points;
367 
368  ellipse_point_at_radian(p, center, axis_a, axis_b,
369  radian_step * (num_points - 1));
370  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_MOVE);
371 
372  radian = 0;
373  for (i = 0; i < num_points; ++i) {
374  ellipse_point_at_radian(p, center, axis_a, axis_b, radian);
375  RT_ADD_VLIST(vhead, p, BN_VLIST_LINE_DRAW);
376 
377  radian += radian_step;
378  }
379 }
int(* ft_bbox)(struct rt_db_internal *, point_t *, point_t *, const struct bn_tol *)
Definition: raytrace.h:2196
Definition: list.h:118
void plot_ellipse(struct bu_list *vhead, const vect_t center, const vect_t axis_a, const vect_t axis_b, int num_points)
const struct bn_tol * tol
Definition: raytrace.h:1927
int approximate_hyperbolic_curve(struct rt_pt_node *pts, fastf_t a, fastf_t b, int num_new_points)
int approximate_parabolic_curve(struct rt_pt_node *pts, fastf_t p, int num_new_points)
fastf_t primitive_curve_count(struct rt_db_internal *ip, const struct rt_view_info *info)
fastf_t curve_spacing
Definition: raytrace.h:1940
double rel
rel dist tol
Definition: raytrace.h:181
Definition: color.c:49
#define RT_ADD_VLIST(hd, pnt, draw)
Definition: raytrace.h:1865
#define BU_ALLOC(_ptr, _type)
Definition: malloc.h:223
#define BN_VLIST_LINE_MOVE
Definition: vlist.h:82
const struct rt_functab * idb_meth
for ft_ifree(), etc.
Definition: raytrace.h:194
#define BN_VLIST_LINE_DRAW
Definition: vlist.h:83
Coord * point
Definition: chull3d.cpp:52
struct rt_pt_node * next
ptr to next pt
Definition: raytrace.h:1912
fastf_t primitive_diagonal_samples(struct rt_db_internal *ip, const struct rt_view_info *info)
fastf_t point_spacing
Definition: raytrace.h:1934
double abs
absolute dist tol
Definition: raytrace.h:180
point_t p
a point
Definition: raytrace.h:1911
Definition: color.c:51
fastf_t primitive_get_absolute_tolerance(const struct rt_tess_tol *ttol, fastf_t rel_to_abs)
#define DEFAULT_REL_TOL
double fastf_t
Definition: defines.h:300
void ellipse_point_at_radian(point_t result, const vect_t center, const vect_t axis_a, const vect_t axis_b, fastf_t radian)
Definition: color.c:50