PostGIS  2.3.8dev-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  *
23  **********************************************************************/
24 
25 #include <float.h>
26 
27 #include "liblwgeom_internal.h"
28 #include "lwgeom_log.h"
29 
30 static void
31 calc_distances_3d(const POINT3D* curr, const POINT3D* points, uint32_t npoints, double* distances)
32 {
33  uint32_t i;
34  for (i = 0; i < npoints; i++)
35  {
36  distances[i] = distance3d_pt_pt(curr, &points[i]);
37  }
38 }
39 
40 static double
41 iterate_3d(POINT3D* curr, const POINT3D* points, uint32_t npoints, double* distances)
42 {
43  uint32_t i;
44  POINT3D next = { 0, 0, 0 };
45  double delta;
46  double denom = 0;
47  char hit = LW_FALSE;
48 
49  calc_distances_3d(curr, points, npoints, distances);
50 
51  for (i = 0; i < npoints; i++)
52  {
53  if (distances[i] == 0)
54  hit = LW_TRUE;
55  else
56  denom += 1.0 / distances[i];
57  }
58 
59  for (i = 0; i < npoints; i++)
60  {
61  if (distances[i] > 0)
62  {
63  next.x += (points[i].x / distances[i]) / denom;
64  next.y += (points[i].y / distances[i]) / denom;
65  next.z += (points[i].z / distances[i]) / denom;
66  }
67  }
68 
69  /* If any of the intermediate points in the calculation is found in the
70  * set of input points, the standard Weiszfeld method gets stuck with a
71  * divide-by-zero.
72  *
73  * To get ourselves out of the hole, we follow an alternate procedure to
74  * get the next iteration, as described in:
75  *
76  * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the
77  * Fermat-Weber location problem." Math. Program., Ser. A 90: 559-566.
78  * DOI 10.1007/s101070100222
79  *
80  * Available online at the time of this writing at
81  * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf
82  */
83  if (hit)
84  {
85  double dx = 0;
86  double dy = 0;
87  double dz = 0;
88  double r_inv;
89  POINT3D alt;
90  for (i = 0; i < npoints; i++)
91  {
92  if (distances[i] > 0)
93  {
94  dx += (points[i].x - curr->x) / distances[i];
95  dy += (points[i].y - curr->y) / distances[i];
96  dz += (points[i].y - curr->z) / distances[i];
97  }
98  }
99 
100  r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz );
101 
102  alt.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x;
103  alt.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y;
104  alt.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z;
105 
106  next = alt;
107  }
108 
109  delta = distance3d_pt_pt(curr, &next);
110 
111  curr->x = next.x;
112  curr->y = next.y;
113  curr->z = next.z;
114 
115  return delta;
116 }
117 
118 static POINT3D
119 init_guess(const POINT3D* points, uint32_t npoints)
120 {
121  POINT3D guess = { 0, 0, 0 };
122  uint32_t i;
123  for (i = 0; i < npoints; i++)
124  {
125  guess.x += points[i].x / npoints;
126  guess.y += points[i].y / npoints;
127  guess.z += points[i].z / npoints;
128  }
129 
130  return guess;
131 }
132 
133 static POINT3D*
134 lwmpoint_extract_points_3d(const LWMPOINT* g, uint32_t* ngeoms)
135 {
136  uint32_t i;
137  uint32_t n = 0;
138  int is_3d = lwgeom_has_z((LWGEOM*) g);
139 
140  POINT3D* points = lwalloc(g->ngeoms * sizeof(POINT3D));
141  for (i = 0; i < g->ngeoms; i++)
142  {
144  if (!lwgeom_is_empty(subg))
145  {
146  getPoint3dz_p(((LWPOINT*) subg)->point, 0, (POINT3DZ*) &points[n++]);
147  if (!is_3d)
148  points[n-1].z = 0.0; /* in case the getPoint functions return NaN in the future for 2d */
149  }
150  }
151 
152  if (ngeoms != NULL)
153  *ngeoms = n;
154 
155  return points;
156 }
157 
158 LWPOINT*
159 lwmpoint_median(const LWMPOINT* g, double tol, uint32_t max_iter, char fail_if_not_converged)
160 {
161  uint32_t npoints; /* we need to count this ourselves so we can exclude empties */
162  uint32_t i;
163  double delta = DBL_MAX;
164  double* distances;
165  POINT3D* points = lwmpoint_extract_points_3d(g, &npoints);
166  POINT3D median;
167 
168  if (npoints == 0)
169  {
170  lwfree(points);
171  return lwpoint_construct_empty(g->srid, 0, 0);
172  }
173 
174  median = init_guess(points, npoints);
175 
176  distances = lwalloc(npoints * sizeof(double));
177 
178  for (i = 0; i < max_iter && delta > tol; i++)
179  {
180  delta = iterate_3d(&median, points, npoints, distances);
181  }
182 
183  lwfree(points);
184  lwfree(distances);
185 
186  if (fail_if_not_converged && delta > tol)
187  {
188  lwerror("Median failed to converge within %g after %d iterations.", tol, max_iter);
189  return NULL;
190  }
191 
192  if (lwgeom_has_z((LWGEOM*) g))
193  {
194  return lwpoint_make3dz(g->srid, median.x, median.y, median.z);
195  }
196  else
197  {
198  return lwpoint_make2d(g->srid, median.x, median.y);
199  }
200 }
201 
202 LWPOINT*
203 lwgeom_median(const LWGEOM* g, double tol, uint32_t max_iter, char fail_if_not_converged)
204 {
205  switch( lwgeom_get_type(g) )
206  {
207  case POINTTYPE:
208  return lwpoint_clone(lwgeom_as_lwpoint(g));
209  case MULTIPOINTTYPE:
210  return lwmpoint_median(lwgeom_as_lwmpoint(g), tol, max_iter, fail_if_not_converged);
211  default:
212  lwerror("Unsupported geometry type in lwgeom_median");
213  return NULL;
214  }
215 }
216 
LWPOINT * lwgeom_median(const LWGEOM *g, double tol, uint32_t max_iter, char fail_if_not_converged)
void lwfree(void *mem)
Definition: lwutil.c:242
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:842
LWPOINT * lwpoint_make2d(int srid, double x, double y)
Definition: lwpoint.c:145
#define MULTIPOINTTYPE
Definition: liblwgeom.h:87
double x
Definition: liblwgeom.h:339
LWPOINT * lwpoint_construct_empty(int srid, char hasz, char hasm)
Definition: lwpoint.c:133
LWPOINT * lwgeom_as_lwpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:93
#define FP_MIN(A, B)
int lwgeom_has_z(const LWGEOM *geom)
Return LW_TRUE if geometry has Z ordinates.
Definition: lwgeom.c:849
int32_t srid
Definition: liblwgeom.h:398
double z
Definition: liblwgeom.h:339
static double iterate_3d(POINT3D *curr, const POINT3D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:41
int getPoint3dz_p(const POINTARRAY *pa, int n, POINT3DZ *point)
Definition: lwgeom_api.c:332
#define LW_FALSE
Definition: liblwgeom.h:76
LWPOINT * lwpoint_clone(const LWPOINT *lwgeom)
Definition: lwpoint.c:221
#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:155
LWPOINT * lwmpoint_median(const LWMPOINT *g, double tol, uint32_t max_iter, char fail_if_not_converged)
static POINT3D * lwmpoint_extract_points_3d(const LWMPOINT *g, uint32_t *ngeoms)
double distance3d_pt_pt(const POINT3D *p1, const POINT3D *p2)
Definition: measures3d.c:795
#define POINTTYPE
LWTYPE numbers, used internally by PostGIS.
Definition: liblwgeom.h:84
LWMPOINT * lwgeom_as_lwmpoint(const LWGEOM *lwgeom)
Definition: lwgeom.c:165
void * lwalloc(size_t size)
Definition: lwutil.c:227
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:1310
int32_t srid
Definition: liblwgeom.h:466
void lwerror(const char *fmt,...)
Write a notice out to the error handler.
Definition: lwutil.c:102
static POINT3D init_guess(const POINT3D *points, uint32_t npoints)
static void calc_distances_3d(const POINT3D *curr, const POINT3D *points, uint32_t npoints, double *distances)
Definition: lwgeom_median.c:31
#define FP_MAX(A, B)