BRL-CAD
polygon.c
Go to the documentation of this file.
1 /* P O L Y G O N . C
2  * BRL-CAD
3  *
4  * Copyright (c) 2013-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 
21 #include "bu/sort.h"
22 #include "bn/polygon.h"
23 #include "bn/plane_struct.h"
24 #include "bn/plane_calc.h"
25 
26 
27 int
28 bn_polygon_area(fastf_t *area, size_t npts, const point_t *pts)
29 {
30  size_t i;
31  vect_t v1, v2, tmp, tot = VINIT_ZERO;
32  plane_t plane_eqn;
33  struct bn_tol tol;
34 
35  if (!pts || !area || npts < 3)
36  return 1;
37  BN_TOL_INIT(&tol);
39  if (bn_mk_plane_3pts(plane_eqn, pts[0], pts[1], pts[2], &tol) == -1)
40  return 1;
41 
42  switch (npts) {
43  case 3:
44  /* Triangular Face - for triangular face T:V0, V1, V2,
45  * area = 0.5 * [(V2 - V0) x (V1 - V0)] */
46  VSUB2(v1, pts[1], pts[0]);
47  VSUB2(v2, pts[2], pts[0]);
48  VCROSS(tot, v2, v1);
49  break;
50  case 4:
51  /* Quadrilateral Face - for planar quadrilateral
52  * Q:V0, V1, V2, V3 with unit normal N,
53  * area = N/2 ⋅ [(V2 - V0) x (V3 - V1)] */
54  VSUB2(v1, pts[2], pts[0]);
55  VSUB2(v2, pts[3], pts[1]);
56  VCROSS(tot, v2, v1);
57  break;
58  default:
59  /* N-Sided Face - compute area using Green's Theorem */
60  for (i = 0; i < npts; i++) {
61  VCROSS(tmp, pts[i], pts[i + 1 == npts ? 0 : i + 1]);
62  VADD2(tot, tot, tmp);
63  }
64  break;
65  }
66  *area = fabs(VDOT(plane_eqn, tot)) * 0.5;
67  return 0;
68 }
69 
70 
71 int
72 bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts)
73 {
74  size_t i;
75  fastf_t x_0 = 0.0;
76  fastf_t x_1 = 0.0;
77  fastf_t y_0 = 0.0;
78  fastf_t y_1 = 0.0;
79  fastf_t z_0 = 0.0;
80  fastf_t z_1 = 0.0;
81  fastf_t a = 0.0;
82  fastf_t signedArea = 0.0;
83 
84  if (!pts || !cent || npts < 3)
85  return 1;
86  /* Calculate Centroid projection for face for x-y-plane */
87  for (i = 0; i < npts-1; i++) {
88  x_0 = pts[i][0];
89  y_0 = pts[i][1];
90  x_1 = pts[i+1][0];
91  y_1 = pts[i+1][1];
92  a = x_0 *y_1 - x_1*y_0;
93  signedArea += a;
94  *cent[0] += (x_0 + x_1)*a;
95  *cent[1] += (y_0 + y_1)*a;
96  }
97  x_0 = pts[i][0];
98  y_0 = pts[i][1];
99  x_1 = pts[0][0];
100  y_1 = pts[0][1];
101  a = x_0 *y_1 - x_1*y_0;
102  signedArea += a;
103  *cent[0] += (x_0 + x_1)*a;
104  *cent[1] += (y_0 + y_1)*a;
105 
106  signedArea *= 0.5;
107  *cent[0] /= (6.0*signedArea);
108  *cent[1] /= (6.0*signedArea);
109 
110  /* calculate Centroid projection for face for x-z-plane */
111 
112  signedArea = 0.0;
113  for (i = 0; i < npts-1; i++) {
114  x_0 = pts[i][0];
115  z_0 = pts[i][2];
116  x_1 = pts[i+1][0];
117  z_1 = pts[i+1][2];
118  a = x_0 *z_1 - x_1*z_0;
119  signedArea += a;
120  *cent[2] += (z_0 + z_1)*a;
121  }
122  x_0 = pts[i][0];
123  z_0 = pts[i][2];
124  x_1 = pts[0][0];
125  z_0 = pts[0][2];
126  a = x_0 *z_1 - x_1*z_0;
127  signedArea += a;
128  *cent[2] += (z_0 + z_1)*a;
129 
130  signedArea *= 0.5;
131  *cent[2] /= (6.0*signedArea);
132  return 0;
133 }
134 
135 
136 int
137 bn_polygon_mk_pts_planes(size_t *npts, point_t **pts, size_t neqs, const plane_t *eqs)
138 {
139  size_t i, j, k, l;
140  if (!npts || !pts || neqs < 4 || !eqs)
141  return 1;
142  /* find all vertices */
143  for (i = 0; i < neqs - 2; i++) {
144  for (j = i + 1; j < neqs - 1; j++) {
145  for (k = j + 1; k < neqs; k++) {
146  point_t pt;
147  int keep_point = 1;
148  if (bn_mkpoint_3planes(pt, eqs[i], eqs[j], eqs[k]) < 0)
149  continue;
150  /* discard pt if it is outside the polyhedron */
151  for (l = 0; l < neqs; l++) {
152  if (l == i || l == j || l == k)
153  continue;
154  if (DIST_PT_PLANE(pt, eqs[l]) > BN_TOL_DIST) {
155  keep_point = 0;
156  break;
157  }
158  }
159  /* found a good point, add it to each of the intersecting faces */
160  if (keep_point) {
161  VMOVE(pts[i][npts[i]], (pt)); npts[i]++;
162  VMOVE(pts[j][npts[j]], (pt)); npts[j]++;
163  VMOVE(pts[k][npts[k]], (pt)); npts[k]++;
164  }
165  }
166  }
167  }
168  return 0;
169 }
170 
171 
172 HIDDEN int
173 sort_ccw(const void *x, const void *y, void *cmp)
174 {
175  vect_t tmp;
176  VCROSS(tmp, ((fastf_t *)x), ((fastf_t *)y));
177  return VDOT(*((point_t *)cmp), tmp);
178 }
179 
180 
181 int
182 bn_polygon_sort_ccw(size_t npts, point_t *pts, plane_t cmp)
183 {
184  if (!pts || npts < 3)
185  return 1;
186  bu_sort(pts, npts, sizeof(point_t), sort_ccw, &cmp);
187  return 0;
188 }
189 
190 
191 /*
192  * Local Variables:
193  * mode: C
194  * tab-width: 8
195  * indent-tabs-mode: t
196  * c-file-style: "stroustrup"
197  * End:
198  * ex: shiftwidth=4 tabstop=8
199  */
#define BN_TOL_INIT(_p)
Definition: tol.h:87
HIDDEN int sort_ccw(const void *x, const void *y, void *cmp)
Definition: polygon.c:173
double dist_sq
dist * dist
Definition: tol.h:74
int bn_mkpoint_3planes(point_t pt, const plane_t a, const plane_t b, const plane_t c)
Given the description of three planes, compute the point of intersection, if any. The direction vecto...
#define HIDDEN
Definition: common.h:86
int bn_polygon_mk_pts_planes(size_t *npts, point_t **pts, size_t neqs, const plane_t *eqs)
Calculate for an array of plane_eqs, which build a polyhedron, the point_t's for each face...
Definition: polygon.c:137
void bu_sort(void *array, size_t nummemb, size_t sizememb, int(*compare)(const void *, const void *, void *), void *context)
Definition: sort.c:110
#define BN_TOL_DIST
Definition: tol.h:109
int bn_polygon_area(fastf_t *area, size_t npts, const point_t *pts)
Functions for working with polygons.
Definition: polygon.c:28
Support for uniform tolerances.
Definition: tol.h:71
int bn_polygon_centroid(point_t *cent, size_t npts, const point_t *pts)
Calculate the centroid of a non self-intersecting polygon.
Definition: polygon.c:72
double fastf_t
Definition: defines.h:300
int bn_mk_plane_3pts(plane_t plane, const point_t a, const point_t b, const point_t c, const struct bn_tol *tol)
int bn_polygon_sort_ccw(size_t npts, point_t *pts, plane_t cmp)
Sort an array of point_ts, building a convex polygon, counter-clockwise.
Definition: polygon.c:182