BRL-CAD
chull.c
Go to the documentation of this file.
1 /* C H U L L . 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  * Copyright 2001 softSurfer, 2012 Dan Sunday
8  * This code may be freely used and modified for any purpose
9  * providing that this copyright notice is included with it.
10  * SoftSurfer makes no warranty for this code, and cannot be held
11  * liable for any real or imagined damage resulting from its use.
12  * Users of this code must verify correctness for their application.
13  *
14  */
15 /** @file chull.c
16  *
17  * This file implements various algorithms for finding convex hull
18  * of point sets in 2D and 3D.
19  */
20 
21 #include "common.h"
22 
23 #include <stdlib.h>
24 
25 #include "bu/malloc.h"
26 #include "bu/sort.h"
27 #include "bn/chull.h"
28 
29 #include "./bn_private.h"
30 
31 
32 /* isLeft(): test if a point is Left|On|Right of an infinite line.
33  * Input: three points L0, L1, and p
34  * Return: >0 for p left of the line through L0 and L1
35  * =0 for p on the line
36  * <0 for p right of the line
37  */
38 #define isLeft(L0, L1, p) ((L1[X] - L0[X])*(p[Y] - L0[Y]) - (p[X] - L0[X])*(L1[Y] - L0[Y]))
39 
40 
41 /* The implementation of Melkman's algorithm for convex hulls of simple
42  * polylines is a translation of softSurfer's implementation:
43  * http://geomalgorithms.com/a12-_hull-3.html
44  */
45 int
46 bn_polyline_2d_chull(point2d_t** hull, const point2d_t* polyline, int n)
47 {
48  int i;
49 
50  /* initialize a deque D[] from bottom to top so that the
51  1st three vertices of P[] are a ccw triangle */
52  point2d_t* D = (point2d_t *)bu_calloc(2*n+1, sizeof(fastf_t)*3, "dequeue");
53 
54  /* hull vertex counter */
55  int h;
56 
57  /* initial bottom and top deque indices */
58  int bot = n-2;
59  int top = bot+3;
60 
61  /* 3rd vertex is a both bottom and top */
62  V2MOVE(D[top], polyline[2]);
63  V2MOVE(D[bot], D[top]);
64  if (isLeft(polyline[0], polyline[1], polyline[2]) > 0) {
65  V2MOVE(D[bot+1],polyline[0]);
66  V2MOVE(D[bot+2],polyline[1]); /* ccw vertices are: 2,0,1,2 */
67  } else {
68  V2MOVE(D[bot+1],polyline[1]);
69  V2MOVE(D[bot+2],polyline[0]); /* ccw vertices are: 2,1,0,2 */
70  }
71 
72  /* compute the hull on the deque D[] */
73  for (i = 3; i < n; i++) { /* process the rest of vertices */
74  /* test if next vertex is inside the deque hull */
75  if ((isLeft(D[bot], D[bot+1], polyline[i]) > 0) &&
76  (isLeft(D[top-1], D[top], polyline[i]) > 0) )
77  continue; /* skip an interior vertex */
78 
79  /* incrementally add an exterior vertex to the deque hull
80  get the rightmost tangent at the deque bot */
81  while (isLeft(D[bot], D[bot+1], polyline[i]) <= 0)
82  bot = bot + 1; /* remove bot of deque */
83  V2MOVE(D[bot-1],polyline[i]); /* insert P[i] at bot of deque */
84  bot = bot - 1;
85 
86  /* get the leftmost tangent at the deque top */
87  while (isLeft(D[top-1], D[top], polyline[i]) <= 0)
88  top = top - 1; /* pop top of deque */
89  V2MOVE(D[top+1],polyline[i]); /* push P[i] onto top of deque */
90  top = top + 1;
91  }
92 
93  /* transcribe deque D[] to the output hull array hull[] */
94 
95  (*hull) = (point2d_t *)bu_calloc(top - bot + 2, sizeof(fastf_t)*3, "hull");
96  for (h=0; h <= (top-bot); h++)
97  V2MOVE((*hull)[h],D[bot + h]);
98 
99  bu_free(D, "free queue");
100  return h-1;
101 }
102 
103 
104 /* bu_sort functions for points */
105 HIDDEN int
106 pnt_compare_2d(const void *pnt1, const void *pnt2, void *UNUSED(arg))
107 {
108  point2d_t *p1 = (point2d_t *)pnt1;
109  point2d_t *p2 = (point2d_t *)pnt2;
110 
111  if (UNLIKELY(NEAR_ZERO((*p2)[0] - (*p1)[0], SMALL_FASTF) && NEAR_ZERO((*p2)[1] - (*p1)[1], SMALL_FASTF))) return 0;
112  if ((*p1)[0] < (*p2)[0]) return 1;
113  if ((*p1)[0] > (*p2)[0]) return -1;
114  if ((*p1)[1] < (*p2)[1]) return 1;
115  if ((*p1)[1] > (*p2)[1]) return -1;
116 
117  /* should never get here */
118 
119  return 0;
120 }
121 
122 
123 int
124 bn_2d_chull(point2d_t **hull, const point2d_t *points_2d, int n)
125 {
126  int i = 0;
127  int retval = 0;
128  point2d_t *points = (point2d_t *)bu_calloc(n + 1, sizeof(point2d_t), "sorted points_2d");
129  const point2d_t *const_points;
130 
131  /* copy points_2d array to something
132  that can be sorted and sort it */
133  for (i = 0; i < n; i++) {
134  V2MOVE(points[i], points_2d[i]);
135  }
136 
137  bu_sort((void *)points, n, sizeof(point2d_t), pnt_compare_2d, NULL);
138 
139  /* Once sorted, the points can be viewed as describing a simple polyline
140  * and the Melkman algorithm works for a simple polyline even if it
141  * isn't closed. */
142  const_points = (const point2d_t *)points; /* quell */
143  retval = bn_polyline_2d_chull(hull, const_points, n);
144 
145  bu_free(points, "free sorted points");
146 
147  return retval;
148 }
149 
150 int
151 bn_3d_coplanar_chull(point_t **hull, const point_t *points_3d, int n)
152 {
153  int ret = 0;
154  int hull_cnt = 0;
155  point_t origin_pnt;
156  vect_t u_axis, v_axis;
157  point2d_t *hull_2d = (point2d_t *)bu_malloc(sizeof(point2d_t *), "hull pointer");
158  point2d_t *points_tmp = (point2d_t *)bu_calloc(n + 1, sizeof(point2d_t), "points_2d");
159 
160  const point2d_t *const_points_tmp;
161 
162  ret += coplanar_2d_coord_sys(&origin_pnt, &u_axis, &v_axis, points_3d, n);
163  ret += coplanar_3d_to_2d(&points_tmp, (const point_t *)&origin_pnt, (const vect_t *)&u_axis, (const vect_t *)&v_axis, points_3d, n);
164 
165  if (ret)
166  return 0;
167 
168  const_points_tmp = (const point2d_t *)points_tmp;
169 
170  hull_cnt = bn_2d_chull(&hull_2d, const_points_tmp, n);
171  (*hull) = (point_t *)bu_calloc(hull_cnt + 1, sizeof(point_t), "hull array");
172 
173  ret = coplanar_2d_to_3d(hull, (const point_t *)&origin_pnt, (const vect_t *)&u_axis, (const vect_t *)&v_axis, (const point2d_t *)hull_2d, hull_cnt);
174 
175  if (ret)
176  return 0;
177 
178  return hull_cnt;
179 }
180 
181 /*
182  * Local Variables:
183  * mode: C
184  * tab-width: 8
185  * indent-tabs-mode: t
186  * c-file-style: "stroustrup"
187  * End:
188  * ex: shiftwidth=4 tabstop=8
189  */
#define SMALL_FASTF
Definition: defines.h:342
Header file for the BRL-CAD common definitions.
#define HIDDEN
Definition: common.h:86
#define isLeft(L0, L1, p)
Definition: chull.c:38
void * bu_malloc(size_t siz, const char *str)
Definition: malloc.c:314
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
static void top()
int coplanar_2d_coord_sys(point_t *origin_pnt, vect_t *u_axis, vect_t *v_axis, const point_t *points_3d, int n)
Find a 2D coordinate system for a set of co-planar 3D points.
Definition: util.c:29
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 UNUSED(parameter)
Definition: common.h:239
HIDDEN int pnt_compare_2d(const void *pnt1, const void *pnt2, void *arg)
Definition: chull.c:106
int bn_3d_coplanar_chull(point_t **hull, const point_t *points_3d, int n)
Find 3D coplanar point convex hull for unordered co-planar point sets.
Definition: chull.c:151
int coplanar_3d_to_2d(point2d_t **points_2d, const point_t *origin_pnt, const vect_t *u_axis, const vect_t *v_axis, const point_t *points_3d, int n)
Find 2D coordinates for a set of co-planar 3D points.
Definition: util.c:93
void bu_free(void *ptr, const char *str)
Definition: malloc.c:328
int coplanar_2d_to_3d(point_t **points_3d, const point_t *origin_pnt, const vect_t *u_axis, const vect_t *v_axis, const point2d_t *points_2d, int n)
Find 3D coordinates for a set of 2D points given a coordinate system.
Definition: util.c:113
int bn_polyline_2d_chull(point2d_t **hull, const point2d_t *polyline, int n)
Routines for the computation of convex hulls in 2D and 3D.
Definition: chull.c:46
double fastf_t
Definition: defines.h:300
int bn_2d_chull(point2d_t **hull, const point2d_t *points_2d, int n)
Find 2D convex hull for unordered co-planar point sets.
Definition: chull.c:124
#define UNLIKELY(expression)
Definition: common.h:282