PostGIS  2.5.0dev-r@@SVN_REVISION@@
lwgeom_median.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  * Copyright 2017 Darafei Praliaskouski <me@komzpa.net>
23  *
24  **********************************************************************/
25 
26 #include <assert.h>
27 #include <float.h>
28 
29 #include "liblwgeom_internal.h"
30 #include "lwgeom_log.h"
31 
32 static void
33 calc_weighted_distances_3d(const POINT3D* curr, const POINT4D* points, uint32_t npoints, double* distances)
34 {
35  uint32_t i;
36  for (i = 0; i < npoints; i++)
37  {
38  distances[i] = distance3d_pt_pt(curr, (POINT3D*)&points[i]) / points[i].m;
39  }
40 }
41 
42 static double
43 iterate_4d(POINT3D* curr, const POINT4D* points, uint32_t npoints, double* distances)
44 {
45  uint32_t i;
46  POINT3D next = { 0, 0, 0 };
47  double delta = 0;
48  double denom = 0;
49  char hit = LW_FALSE;
50 
51  calc_weighted_distances_3d(curr, points, npoints, distances);
52 
53  for (i = 0; i < npoints; i++)
54  {
55  /* we need to use lower epsilon than in FP_IS_ZERO in the loop for calculation to converge */
56  if (distances[i] > DBL_EPSILON)
57  {
58  next.x += points[i].x / distances[i];
59  next.y += points[i].y / distances[i];
60  next.z += points[i].z / distances[i];
61  denom += 1.0 / distances[i];
62  }
63  else
64  {
65  hit = LW_TRUE;
66  }
67  }
68  /* negative weight shouldn't get here */
69  assert(denom >= 0);
70 
71  /* denom is zero in case of multipoint of single point when we've converged perfectly */
72  if (denom > DBL_EPSILON)
73  {
74  next.x /= denom;
75  next.y /= denom;
76  next.z /= denom;
77 
78  /* If any of the intermediate points in the calculation is found in the
79  * set of input points, the standard Weiszfeld method gets stuck with a
80  * divide-by-zero.
81  *
82  * To get ourselves out of the hole, we follow an alternate procedure to
83  * get the next iteration, as described in:
84  *
85  * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the
86  * Fermat-Weber location problem." Math. Program., Ser. A 90: 559-566.
87  * DOI 10.1007/s101070100222
88  *
89  * Available online at the time of this writing at
90  * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf
91  */
92  if (hit)
93  {
94  double dx = 0;
95  double dy = 0;
96  double dz = 0;
97  double r_inv;
98  for (i = 0; i < npoints; i++)
99  {
100  if (distances[i] > DBL_EPSILON)
101  {
102  dx += (points[i].x - curr->x) / distances[i];
103  dy += (points[i].y - curr->y) / distances[i];
104  dz += (points[i].z - curr->z) / distances[i];
105  }
106  }
107 
108  r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
109 
110  next.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
111  next.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
112  next.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
113  }
114 
115  delta = distance3d_pt_pt(curr, &next);
116 
117  curr->x = next.x;
118  curr->y = next.y;
119  curr->z = next.z;
120  }
121 
122  return delta;
123 }
124 
125 static POINT3D
126 init_guess(const POINT4D* points, uint32_t npoints)
127 {
128  assert(npoints > 0);
129  POINT3D guess = { 0, 0, 0 };
130  double mass = 0;
131  uint32_t i;
132  for (i = 0; i < npoints; i++)
133  {
134  guess.x += points[i].x * points[i].m;
135  guess.y += points[i].y * points[i].m;
136  guess.z += points[i].z * points[i].m;
137  mass += points[i].m;
138  }
139  guess.x /= mass;
140  guess.y /= mass;
141  guess.z /= mass;
142  return guess;
143 }
144 
145 static POINT4D*
146 lwmpoint_extract_points_4d(const LWMPOINT* g, uint32_t* npoints, int* input_empty)
147 {
148  uint32_t i;
149  uint32_t n = 0;
150  POINT4D* points = lwalloc(g->ngeoms * sizeof(POINT4D));
151  int has_m = lwgeom_has_m((LWGEOM*) g);
152 
153  for (i = 0; i < g->ngeoms; i++)
154  {
156  if (!lwgeom_is_empty(subg))
157  {
158  *input_empty = LW_FALSE;
159  if (!getPoint4d_p(((LWPOINT*) subg)->point, 0, &points[n]))
160  {
161  lwerror("Geometric median: getPoint4d_p reported failure on point (POINT(%g %g %g %g), number %d of %d in input).", points[n].x, points[n].y, points[n].z, points[n].m, i, g->ngeoms);
162  lwfree(points);
163  return NULL;
164  }
165  if (has_m)
166  {
167  /* This limitation on non-negativity can be lifted
168  * by replacing Weiszfeld algorithm with different one.
169  * Possible option described in:
170  *
171  * Drezner, Zvi & O. Wesolowsky, George. (1991).
172  * The Weber Problem On The Plane With Some Negative Weights.
173  * INFOR. Information Systems and Operational Research.
174  * 29. 10.1080/03155986.1991.11732158.
175  */
176  if (points[n].m < 0)
177  {
178  lwerror("Geometric median input contains points with negative weights (POINT(%g %g %g %g), number %d of %d in input). Implementation can't guarantee global minimum convergence.", points[n].x, points[n].y, points[n].z, points[n].m, i, g->ngeoms);
179  lwfree(points);
180  return NULL;
181  }
182 
183  /* points with zero weight are not going to affect calculation, drop them early */
184  if (points[n].m > DBL_EPSILON) n++;
185  }
186  else
187  {
188  points[n].m = 1.0;
189  n++;
190  }
191  }
192  }
193 
194 #if PARANOIA_LEVEL > 0
195  /* check Z=0 for 2D inputs*/
196  if (!lwgeom_has_z((LWGEOM*) g)) for (i = 0; i < n; i++) assert(points[i].z == 0);
197 #endif
198 
199  *npoints = n;
200  return points;
201 }
202 
203 
204 LWPOINT*
205 lwmpoint_median(const LWMPOINT* g, double tol, uint32_t max_iter, char fail_if_not_converged)
206 {
207  /* m ordinate is considered weight, if defined */
208  uint32_t npoints = 0; /* we need to count this ourselves so we can exclude empties and weightless points */
209  uint32_t i;
210  int input_empty = LW_TRUE;
211  double delta = DBL_MAX;
212  POINT3D median;
213  double* distances = NULL;
214  POINT4D* points = lwmpoint_extract_points_4d(g, &npoints, &input_empty);
215 
216  /* input validation failed, error reported already */
217  if (points == NULL) return NULL;
218 
219  if (npoints == 0)
220  {
221  lwfree(points);
222  if (input_empty)
223  {
224  return lwpoint_construct_empty(g->srid, 0, 0);
225  }
226  else
227  {
228  lwerror("Median failed to find non-empty input points with positive weight.");
229  return NULL;
230  }
231  }
232 
233  median = init_guess(points, npoints);
234 
235  distances = lwalloc(npoints * sizeof(double));
236 
237  for (i = 0; i < max_iter && delta > tol; i++)
238  {
239  delta = iterate_4d(&median, points, npoints, distances);
240  }
241 
242  lwfree(distances);
243  lwfree(points);
244 
245  if (fail_if_not_converged && delta > tol)
246  {
247  lwerror("Median failed to converge within %g after %d iterations.", tol, max_iter);
248  return NULL;
249  }
250 
251  if (lwgeom_has_z((LWGEOM*) g))
252  {
253  return lwpoint_make3dz(g->srid, median.x, median.y, median.z);
254  }
255  else
256  {
257  return lwpoint_make2d(g->srid, median.x, median.y);
258  }
259 }
260 
261 LWPOINT*
262 lwgeom_median(const LWGEOM* g, double tol, uint32_t max_iter, char fail_if_not_converged)
263 {
264  switch( lwgeom_get_type(g) )
265  {
266  case POINTTYPE:
267  return lwpoint_clone(lwgeom_as_lwpoint(g));
268  case MULTIPOINTTYPE:
269  return lwmpoint_median(lwgeom_as_lwmpoint(g), tol, max_iter, fail_if_not_converged);
270  default:
271  lwerror("Unsupported geometry type in lwgeom_median");
272  return NULL;
273  }
274 }
275 
double x
Definition: liblwgeom.h:351
LWPOINT * lwgeom_median(const LWGEOM *g, double tol, uint32_t max_iter, char fail_if_not_converged)
double m
Definition: liblwgeom.h:351
uint32_t ngeoms
Definition: liblwgeom.h:467
void lwfree(void *mem)
Definition: lwutil.c:244
double y
Definition: liblwgeom.h:339
LWGEOM * lwcollection_getsubgeom(LWCOLLECTION *col, int gnum)
Definition: lwcollection.c:113
uint32_t lwgeom_get_type(const LWGEOM *geom)
Return LWTYPE number.
Definition: lwgeom.c:916
static POINT3D init_guess(const POINT4D *points, uint32_t npoints)
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:163
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:239
double x
Definition: liblwgeom.h:339
static double iterate_4d(POINT3D *curr, const POINT4D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:43
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:151
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:160
#define FP_MIN(A, B)
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:923
double z
Definition: liblwgeom.h:339
unsigned int uint32_t
Definition: uthash.h:78
#define LW_FALSE
Definition: liblwgeom.h:76
static void calc_weighted_distances_3d(const POINT3D *curr, const POINT4D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:33
#define LW_TRUE
Return types for functions with status returns.
Definition: liblwgeom.h:75
LWPOINT * lwpoint_make3dz(int srid, double x, double y, double z)
Definition: lwpoint.c:173
LWPOINT * lwmpoint_median(const LWMPOINT *g, double tol, uint32_t max_iter, char fail_if_not_converged)
double z
Definition: liblwgeom.h:351
tuple x
Definition: pixval.py:53
int getPoint4d_p(const POINTARRAY *pa, uint32_t n, POINT4D *point)
Definition: lwgeom_api.c:113
static POINT4D * lwmpoint_extract_points_4d(const LWMPOINT *g, uint32_t *npoints, int *input_empty)
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:792
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:232
void * lwalloc(size_t size)
Definition: lwutil.c:229
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:1386
double y
Definition: liblwgeom.h:351
int32_t srid
Definition: liblwgeom.h:466
int lwgeom_has_m(const LWGEOM *geom)
Return LW_TRUE if geometry has M ordinates.
Definition: lwgeom.c:930
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:190
tuple y
Definition: pixval.py:54
#define FP_MAX(A, B)