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 (C) 2017 Danny Götte <danny.goette@fem.tu-ilmenau.de>
22 *
23 **********************************************************************/
24
25 #include "postgres.h"
26
27 #include "../postgis_config.h"
28
29 #include <math.h>
30
31 #include "liblwgeom.h" /* For standard geometry types. */
32 #include "lwgeom_pg.h" /* For pg macros. */
33 #include "lwgeom_transform.h" /* For SRID functions */
34
35 Datum geography_centroid(PG_FUNCTION_ARGS);
36
37 /* internal functions */
38 LWPOINT *geography_centroid_from_wpoints(const int32_t srid, const POINT3DM *points, const uint32_t size);
39 LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s);
40 LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s);
41 LWPOINT *cart_to_lwpoint(const double_t x_sum,
42 const double_t y_sum,
43 const double_t z_sum,
44 const double_t weight_sum,
45 const int32_t srid);
46 POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat);
47
48 /**
49 * geography_centroid(GSERIALIZED *g)
50 * returns centroid as point
51 */
52 PG_FUNCTION_INFO_V1(geography_centroid);
geography_centroid(PG_FUNCTION_ARGS)53 Datum geography_centroid(PG_FUNCTION_ARGS)
54 {
55 LWGEOM *lwgeom = NULL;
56 LWGEOM *lwgeom_out = NULL;
57 LWPOINT *lwpoint_out = NULL;
58 GSERIALIZED *g = NULL;
59 GSERIALIZED *g_out = NULL;
60 int32_t srid;
61 bool use_spheroid = true;
62 SPHEROID s;
63
64 /* Get our geometry object loaded into memory. */
65 g = PG_GETARG_GSERIALIZED_P(0);
66 lwgeom = lwgeom_from_gserialized(g);
67
68 if (g == NULL)
69 {
70 PG_RETURN_NULL();
71 }
72
73 srid = lwgeom_get_srid(lwgeom);
74
75 /* on empty input, return empty output */
76 if (gserialized_is_empty(g))
77 {
78 LWCOLLECTION* empty = lwcollection_construct_empty(COLLECTIONTYPE, srid, 0, 0);
79 lwgeom_out = lwcollection_as_lwgeom(empty);
80 g_out = geography_serialize(lwgeom_out);
81 PG_RETURN_POINTER(g_out);
82 }
83
84 /* Initialize spheroid */
85 spheroid_init_from_srid(srid, &s);
86
87 /* Set to sphere if requested */
88 use_spheroid = PG_GETARG_BOOL(1);
89 if ( ! use_spheroid )
90 s.a = s.b = s.radius;
91
92 switch (lwgeom_get_type(lwgeom))
93 {
94
95 case POINTTYPE:
96 {
97 /* centroid of a point is itself */
98 PG_RETURN_POINTER(g);
99 break;
100 }
101
102 case MULTIPOINTTYPE:
103 {
104 LWMPOINT* mpoints = lwgeom_as_lwmpoint(lwgeom);
105
106 /* average between all points */
107 uint32_t size = mpoints->ngeoms;
108 POINT3DM* points = palloc(size*sizeof(POINT3DM));
109
110 uint32_t i;
111 for (i = 0; i < size; i++) {
112 points[i].x = lwpoint_get_x(mpoints->geoms[i]);
113 points[i].y = lwpoint_get_y(mpoints->geoms[i]);
114 points[i].m = 1;
115 }
116
117 lwpoint_out = geography_centroid_from_wpoints(srid, points, size);
118 pfree(points);
119 break;
120 }
121
122 case LINETYPE:
123 {
124 LWLINE* line = lwgeom_as_lwline(lwgeom);
125
126 /* reuse mline function */
127 LWMLINE* mline = lwmline_construct_empty(srid, 0, 0);
128 lwmline_add_lwline(mline, line);
129
130 lwpoint_out = geography_centroid_from_mline(mline, &s);
131 lwmline_free(mline);
132 break;
133 }
134
135 case MULTILINETYPE:
136 {
137 LWMLINE* mline = lwgeom_as_lwmline(lwgeom);
138 lwpoint_out = geography_centroid_from_mline(mline, &s);
139 break;
140 }
141
142 case POLYGONTYPE:
143 {
144 LWPOLY* poly = lwgeom_as_lwpoly(lwgeom);
145
146 /* reuse mpoly function */
147 LWMPOLY* mpoly = lwmpoly_construct_empty(srid, 0, 0);
148 lwmpoly_add_lwpoly(mpoly, poly);
149
150 lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
151 lwmpoly_free(mpoly);
152 break;
153 }
154
155 case MULTIPOLYGONTYPE:
156 {
157 LWMPOLY* mpoly = lwgeom_as_lwmpoly(lwgeom);
158 lwpoint_out = geography_centroid_from_mpoly(mpoly, use_spheroid, &s);
159 break;
160 }
161 default:
162 elog(ERROR, "ST_Centroid(geography) unhandled geography type");
163 PG_RETURN_NULL();
164 }
165
166 PG_FREE_IF_COPY(g, 0);
167
168 lwgeom_out = lwpoint_as_lwgeom(lwpoint_out);
169 g_out = geography_serialize(lwgeom_out);
170
171 PG_RETURN_POINTER(g_out);
172 }
173
174
175 /**
176 * Convert lat-lon-points to x-y-z-coordinates, calculate a weighted average
177 * point and return lat-lon-coordinated
178 */
179 LWPOINT *
geography_centroid_from_wpoints(const int32_t srid,const POINT3DM * points,const uint32_t size)180 geography_centroid_from_wpoints(const int32_t srid, const POINT3DM *points, const uint32_t size)
181 {
182 double_t x_sum = 0;
183 double_t y_sum = 0;
184 double_t z_sum = 0;
185 double_t weight_sum = 0;
186
187 double_t weight = 1;
188 POINT3D* point;
189
190 uint32_t i;
191 for (i = 0; i < size; i++ )
192 {
193 point = lonlat_to_cart(points[i].x, points[i].y);
194 weight = points[i].m;
195
196 x_sum += point->x * weight;
197 y_sum += point->y * weight;
198 z_sum += point->z * weight;
199
200 weight_sum += weight;
201
202 lwfree(point);
203 }
204
205 return cart_to_lwpoint(x_sum, y_sum, z_sum, weight_sum, srid);
206 }
207
lonlat_to_cart(const double_t raw_lon,const double_t raw_lat)208 POINT3D* lonlat_to_cart(const double_t raw_lon, const double_t raw_lat)
209 {
210 double_t lat, lon;
211 double_t sin_lat;
212
213 POINT3D* point = lwalloc(sizeof(POINT3D));;
214
215 // prepare coordinate for trigonometric functions from [-90, 90] -> [0, pi]
216 lat = (raw_lat + 90) / 180 * M_PI;
217
218 // prepare coordinate for trigonometric functions from [-180, 180] -> [-pi, pi]
219 lon = raw_lon / 180 * M_PI;
220
221 /* calculate value only once */
222 sin_lat = sinl(lat);
223
224 /* convert to 3D cartesian coordinates */
225 point->x = sin_lat * cosl(lon);
226 point->y = sin_lat * sinl(lon);
227 point->z = cosl(lat);
228
229 return point;
230 }
231
232 LWPOINT *
cart_to_lwpoint(const double_t x_sum,const double_t y_sum,const double_t z_sum,const double_t weight_sum,const int32_t srid)233 cart_to_lwpoint(const double_t x_sum,
234 const double_t y_sum,
235 const double_t z_sum,
236 const double_t weight_sum,
237 const int32_t srid)
238 {
239 double_t x = x_sum / weight_sum;
240 double_t y = y_sum / weight_sum;
241 double_t z = z_sum / weight_sum;
242
243 /* x-y-z vector length */
244 double_t r = sqrtl(powl(x, 2) + powl(y, 2) + powl(z, 2));
245
246 double_t lon = atan2l(y, x) * 180 / M_PI;
247 double_t lat = acosl(z / r) * 180 / M_PI - 90;
248
249 return lwpoint_make2d(srid, lon, lat);
250 }
251
252 /**
253 * Split lines into segments and calculate with middle of segment as weighted
254 * point.
255 */
geography_centroid_from_mline(const LWMLINE * mline,SPHEROID * s)256 LWPOINT* geography_centroid_from_mline(const LWMLINE* mline, SPHEROID* s)
257 {
258 double_t tolerance = 0.0;
259 uint32_t size = 0;
260 uint32_t i, k, j = 0;
261 POINT3DM* points;
262 LWPOINT* result;
263
264 /* get total number of points */
265 for (i = 0; i < mline->ngeoms; i++) {
266 size += (mline->geoms[i]->points->npoints - 1) * 2;
267 }
268
269 points = palloc(size*sizeof(POINT3DM));
270
271 for (i = 0; i < mline->ngeoms; i++) {
272 LWLINE* line = mline->geoms[i];
273
274 /* add both points of line segment as weighted point */
275 for (k = 0; k < line->points->npoints - 1; k++) {
276 const POINT2D* p1 = getPoint2d_cp(line->points, k);
277 const POINT2D* p2 = getPoint2d_cp(line->points, k+1);
278 double_t weight;
279
280 /* use line-segment length as weight */
281 LWPOINT* lwp1 = lwpoint_make2d(mline->srid, p1->x, p1->y);
282 LWPOINT* lwp2 = lwpoint_make2d(mline->srid, p2->x, p2->y);
283 LWGEOM* lwgeom1 = lwpoint_as_lwgeom(lwp1);
284 LWGEOM* lwgeom2 = lwpoint_as_lwgeom(lwp2);
285 lwgeom_set_geodetic(lwgeom1, LW_TRUE);
286 lwgeom_set_geodetic(lwgeom2, LW_TRUE);
287
288 /* use point distance as weight */
289 weight = lwgeom_distance_spheroid(lwgeom1, lwgeom2, s, tolerance);
290
291 points[j].x = p1->x;
292 points[j].y = p1->y;
293 points[j].m = weight;
294 j++;
295
296 points[j].x = p2->x;
297 points[j].y = p2->y;
298 points[j].m = weight;
299 j++;
300
301 lwgeom_free(lwgeom1);
302 lwgeom_free(lwgeom2);
303 }
304 }
305
306 result = geography_centroid_from_wpoints(mline->srid, points, size);
307 pfree(points);
308 return result;
309 }
310
311
312 /**
313 * Split polygons into triangles and use centroid of the triangle with the
314 * triangle area as weight to calculate the centroid of a (multi)polygon.
315 */
geography_centroid_from_mpoly(const LWMPOLY * mpoly,bool use_spheroid,SPHEROID * s)316 LWPOINT* geography_centroid_from_mpoly(const LWMPOLY* mpoly, bool use_spheroid, SPHEROID* s)
317 {
318 uint32_t size = 0;
319 uint32_t i, ir, ip, j = 0;
320 POINT3DM* points;
321 POINT4D* reference_point = NULL;
322 LWPOINT* result = NULL;
323
324 for (ip = 0; ip < mpoly->ngeoms; ip++) {
325 for (ir = 0; ir < mpoly->geoms[ip]->nrings; ir++) {
326 size += mpoly->geoms[ip]->rings[ir]->npoints - 1;
327 }
328 }
329
330 points = palloc(size*sizeof(POINT3DM));
331
332
333 /* use first point as reference to create triangles */
334 reference_point = (POINT4D*) getPoint2d_cp(mpoly->geoms[0]->rings[0], 0);
335
336 for (ip = 0; ip < mpoly->ngeoms; ip++) {
337 LWPOLY* poly = mpoly->geoms[ip];
338
339 for (ir = 0; ir < poly->nrings; ir++) {
340 POINTARRAY* ring = poly->rings[ir];
341
342 /* split into triangles (two points + reference point) */
343 for (i = 0; i < ring->npoints - 1; i++) {
344 const POINT4D* p1 = (const POINT4D*) getPoint2d_cp(ring, i);
345 const POINT4D* p2 = (const POINT4D*) getPoint2d_cp(ring, i+1);
346 LWPOLY* poly_tri;
347 LWGEOM* geom_tri;
348 double_t weight;
349 POINT3DM triangle[3];
350 LWPOINT* tri_centroid;
351
352 POINTARRAY* pa = ptarray_construct_empty(0, 0, 4);
353 ptarray_insert_point(pa, p1, 0);
354 ptarray_insert_point(pa, p2, 1);
355 ptarray_insert_point(pa, reference_point, 2);
356 ptarray_insert_point(pa, p1, 3);
357
358 poly_tri = lwpoly_construct_empty(mpoly->srid, 0, 0);
359 lwpoly_add_ring(poly_tri, pa);
360
361 geom_tri = lwpoly_as_lwgeom(poly_tri);
362 lwgeom_set_geodetic(geom_tri, LW_TRUE);
363
364 /* Calculate the weight of the triangle. If counter clockwise,
365 * the weight is negative (e.g. for holes in polygons)
366 */
367
368 if ( use_spheroid )
369 weight = lwgeom_area_spheroid(geom_tri, s);
370 else
371 weight = lwgeom_area_sphere(geom_tri, s);
372
373
374 triangle[0].x = p1->x;
375 triangle[0].y = p1->y;
376 triangle[0].m = 1;
377
378 triangle[1].x = p2->x;
379 triangle[1].y = p2->y;
380 triangle[1].m = 1;
381
382 triangle[2].x = reference_point->x;
383 triangle[2].y = reference_point->y;
384 triangle[2].m = 1;
385
386 /* get center of triangle */
387 tri_centroid = geography_centroid_from_wpoints(mpoly->srid, triangle, 3);
388
389 points[j].x = lwpoint_get_x(tri_centroid);
390 points[j].y = lwpoint_get_y(tri_centroid);
391 points[j].m = weight;
392 j++;
393
394 lwpoint_free(tri_centroid);
395 lwgeom_free(geom_tri);
396 }
397 }
398 }
399 result = geography_centroid_from_wpoints(mpoly->srid, points, size);
400 pfree(points);
401 return result;
402 }
403