PostGIS  2.5.7dev-r@@SVN_REVISION@@
lwboundingcircle.c
Go to the documentation of this file.
1 /**********************************************************************
2  *
3  * PostGIS - Spatial Types for PostgreSQL
4  * http://postgis.net
5  *
6  * PostGIS is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU General Public License as published by
8  * the Free Software Foundation, either version 2 of the License, or
9  * (at your option) any later version.
10  *
11  * PostGIS is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU General Public License for more details.
15  *
16  * You should have received a copy of the GNU General Public License
17  * along with PostGIS. If not, see <http://www.gnu.org/licenses/>.
18  *
19  **********************************************************************
20  *
21  * Copyright 2015 Daniel Baston <dbaston@gmail.com>
22  *
23  **********************************************************************/
24 
25 
26 #include <string.h>
27 #include "liblwgeom_internal.h"
28 
29 typedef struct {
30  const POINT2D* p1;
31  const POINT2D* p2;
32  const POINT2D* p3;
34 
35 static SUPPORTING_POINTS*
37 {
39  s->p1 = NULL;
40  s->p2 = NULL;
41  s->p3 = NULL;
42 
43  return s;
44 }
45 
46 static void
48 {
49  lwfree(s);
50 }
51 
52 static uint32_t
54 {
55  uint32_t N = 0;
56 
57  if (support->p1 != NULL)
58  N++;
59  if (support->p2 != NULL)
60  N++;
61  if (support->p3 != NULL)
62  N++;
63 
64  return N;
65 }
66 
67 static int
69 {
70  switch(num_supporting_points(support))
71  {
72  case 0: support->p1 = p;
73  break;
74  case 1: support->p2 = p;
75  break;
76  case 2: support->p3 = p;
77  break;
78  default: return LW_FAILURE;
79  }
80 
81  return LW_SUCCESS;
82 }
83 
84 static int
86 {
87  if (!c)
88  return LW_FALSE;
89 
90  if (distance2d_pt_pt(p, c->center) - c->radius > DBL_EPSILON)
91  return LW_FALSE;
92 
93  return LW_TRUE;
94 }
95 
96 static double
97 det(double m00, double m01, double m10, double m11)
98 {
99  return m00 * m11 - m01 * m10;
100 }
101 
102 static void
103 circumcenter(const POINT2D* a, const POINT2D* b, const POINT2D* c, POINT2D* result)
104 {
105  double cx = c->x;
106  double cy = c->y;
107  double ax = a->x - cx;
108  double ay = a->y - cy;
109  double bx = b->x - cx;
110  double by = b->y - cy;
111 
112  double denom = 2 * det(ax, ay, bx, by);
113  double numx = det(ay, ax * ax + ay * ay, by, bx * bx + by * by);
114  double numy = det(ax, ax * ax + ay * ay, bx, bx * bx + by * by);
115 
116  result->x = cx - numx / denom;
117  result->y = cy + numy / denom;
118 }
119 
120 static void
122 {
123  mbc->radius = 0;
124  mbc->center->x = support->p1->x;
125  mbc->center->y = support->p1->y;
126 }
127 
128 static void
130 {
131  double d1, d2;
132 
133  mbc->center->x = 0.5*(support->p1->x + support->p2->x);
134  mbc->center->y = 0.5*(support->p1->y + support->p2->y);
135 
136  d1 = distance2d_pt_pt(mbc->center, support->p1);
137  d2 = distance2d_pt_pt(mbc->center, support->p2);
138 
139  mbc->radius = FP_MAX(d1, d2);
140 }
141 
142 static void
144 {
145  double d1, d2, d3;
146  circumcenter(support->p1, support->p2, support->p3, mbc->center);
147 
148  d1 = distance2d_pt_pt(mbc->center, support->p1);
149  d2 = distance2d_pt_pt(mbc->center, support->p2);
150  d3 = distance2d_pt_pt(mbc->center, support->p3);
151 
152  mbc->radius = FP_MAX(FP_MAX(d1, d2), d3);
153 }
154 
155 static int
157 {
158  switch(num_supporting_points(support))
159  {
160  case 0: break;
161  case 1: calculate_mbc_1(support, mbc);
162  break;
163  case 2: calculate_mbc_2(support, mbc);
164  break;
165  case 3: calculate_mbc_3(support, mbc);
166  break;
167  default: return LW_FAILURE;
168  }
169 
170  return LW_SUCCESS;
171 }
172 
173 static int
174 calculate_mbc(const POINT2D** points, uint32_t max_n, SUPPORTING_POINTS* support, LWBOUNDINGCIRCLE* mbc)
175 {
176  uint32_t i;
177 
178  if(!calculate_mbc_from_support(support, mbc))
179  {
180  return LW_FAILURE;
181  }
182 
183  if (num_supporting_points(support) == 3)
184  {
185  /* If we're entering the function with three supporting points already, our circle
186  * is already fully constrained - we couldn't add another supporting point if we
187  * needed to. So, there's no point in proceeding further. Welzl (1991) provides
188  * a much better explanation of why this works.
189  * */
190  return LW_SUCCESS;
191  }
192 
193  for (i = 0; i < max_n; i++)
194  {
195  if (!point_inside_circle(points[i], mbc))
196  {
197  /* We've run into a point that isn't inside our circle. To fix this, we'll
198  * go back in time, and re-run the algorithm for each point we've seen so
199  * far, with the constraint that the current point must be on the boundary
200  * of the circle. Then, we'll continue on in this loop with the modified
201  * circle that by definition includes the current point. */
202  SUPPORTING_POINTS next_support;
203  memcpy(&next_support, support, sizeof(SUPPORTING_POINTS));
204 
205  add_supporting_point(&next_support, points[i]);
206  if (!calculate_mbc(points, i, &next_support, mbc))
207  {
208  return LW_FAILURE;
209  }
210  }
211  }
212 
213  return LW_SUCCESS;
214 }
215 
216 static LWBOUNDINGCIRCLE*
218 {
220  c->center = lwalloc(sizeof(POINT2D));
221 
222  c->radius = 0.0;
223  c->center->x = 0.0;
224  c->center->y = 0.0;
225 
226  return c;
227 }
228 
229 void
231 {
232  lwfree(c->center);
233  lwfree(c);
234 }
235 
238 {
239  SUPPORTING_POINTS* support;
240  LWBOUNDINGCIRCLE* result;
241  LWPOINTITERATOR* it;
242  uint32_t num_points;
243  POINT2D** points;
244  POINT4D p;
245  uint32_t i;
246  int success;
247 
248  if(g == NULL || lwgeom_is_empty(g))
249  return LW_FAILURE;
250 
251  num_points = lwgeom_count_vertices(g);
252  it = lwpointiterator_create(g);
253  points = lwalloc(num_points * sizeof(POINT2D*));
254  for (i = 0; i < num_points; i++)
255  {
256  if(!lwpointiterator_next(it, &p))
257  {
258  uint32_t j;
259  for (j = 0; j < i; j++)
260  {
261  lwfree(points[j]);
262  }
264  lwfree(points);
265  return LW_FAILURE;
266  }
267 
268  points[i] = lwalloc(sizeof(POINT2D));
269  points[i]->x = p.x;
270  points[i]->y = p.y;
271  }
273 
274  support = supporting_points_create();
275  result = lwboundingcircle_create();
276  /* Technically, a randomized algorithm would demand that we shuffle the input points
277  * before we call calculate_mbc(). However, we make the (perhaps poor) assumption that
278  * the order we happen to find the points is as good as random, or close enough.
279  * */
280  success = calculate_mbc((const POINT2D**) points, num_points, support, result);
281 
282  for (i = 0; i < num_points; i++)
283  {
284  lwfree(points[i]);
285  }
286  lwfree(points);
287  supporting_points_destroy(support);
288 
289  if (!success)
290  return NULL;
291 
292  return result;
293 }
char * s
Definition: cu_in_wkt.c:23
#define LW_FALSE
Definition: liblwgeom.h:77
double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2)
Definition: measures.c:2314
LWPOINTITERATOR * lwpointiterator_create(const LWGEOM *g)
Create a new LWPOINTITERATOR over supplied LWGEOM*.
Definition: lwiterator.c:244
#define LW_FAILURE
Definition: liblwgeom.h:79
int lwpointiterator_next(LWPOINTITERATOR *s, POINT4D *p)
Attempts to assign the next point in the iterator to p, and advances the iterator to the next point.
Definition: lwiterator.c:212
#define LW_SUCCESS
Definition: liblwgeom.h:80
void lwpointiterator_destroy(LWPOINTITERATOR *s)
Free all memory associated with the iterator.
Definition: lwiterator.c:269
uint32_t lwgeom_count_vertices(const LWGEOM *geom)
Count the total number of vertices in any LWGEOM.
Definition: lwgeom.c:1235
void lwfree(void *mem)
Definition: lwutil.c:244
int lwgeom_is_empty(const LWGEOM *geom)
Return true or false depending on whether a geometry is an "empty" geometry (no vertices members)
Definition: lwgeom.c:1393
void * lwalloc(size_t size)
Definition: lwutil.c:229
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:76
#define FP_MAX(A, B)
static int calculate_mbc_from_support(SUPPORTING_POINTS *support, LWBOUNDINGCIRCLE *mbc)
static double det(double m00, double m01, double m10, double m11)
static void supporting_points_destroy(SUPPORTING_POINTS *s)
static int calculate_mbc(const POINT2D **points, uint32_t max_n, SUPPORTING_POINTS *support, LWBOUNDINGCIRCLE *mbc)
static SUPPORTING_POINTS * supporting_points_create()
static LWBOUNDINGCIRCLE * lwboundingcircle_create()
LWBOUNDINGCIRCLE * lwgeom_calculate_mbc(const LWGEOM *g)
static void calculate_mbc_3(const SUPPORTING_POINTS *support, LWBOUNDINGCIRCLE *mbc)
static int add_supporting_point(SUPPORTING_POINTS *support, const POINT2D *p)
static void calculate_mbc_1(const SUPPORTING_POINTS *support, LWBOUNDINGCIRCLE *mbc)
static uint32_t num_supporting_points(SUPPORTING_POINTS *support)
void lwboundingcircle_destroy(LWBOUNDINGCIRCLE *c)
static void circumcenter(const POINT2D *a, const POINT2D *b, const POINT2D *c, POINT2D *result)
static void calculate_mbc_2(const SUPPORTING_POINTS *support, LWBOUNDINGCIRCLE *mbc)
static int point_inside_circle(const POINT2D *p, const LWBOUNDINGCIRCLE *c)
POINT2D * center
Definition: liblwgeom.h:1650
double y
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:331
double x
Definition: liblwgeom.h:355
double y
Definition: liblwgeom.h:355
const POINT2D * p2
const POINT2D * p3
const POINT2D * p1
unsigned int uint32_t
Definition: uthash.h:78