From: Daniel Baston Date: Sat, 5 Mar 2016 01:27:32 +0000 (+0000) Subject: Add new files for #3364 X-Git-Tag: 2.3.0beta1~186 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=40dd6f8f66b5f19206d58e8d40772e8910556583;p=postgis Add new files for #3364 git-svn-id: http://svn.osgeo.org/postgis/trunk@14747 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/doc/html/image_src/st_geometricmedian01.wkt b/doc/html/image_src/st_geometricmedian01.wkt new file mode 100644 index 000000000..73e122728 --- /dev/null +++ b/doc/html/image_src/st_geometricmedian01.wkt @@ -0,0 +1,4 @@ +Style6;POINT(25 25) +Style5;POINT(65 65) +Style7-smallpoint;MULTIPOINT((20 20), (20 30), (30 20), (190 190)) + diff --git a/liblwgeom/lwgeom_median.c b/liblwgeom/lwgeom_median.c new file mode 100644 index 000000000..e8c528261 --- /dev/null +++ b/liblwgeom/lwgeom_median.c @@ -0,0 +1,216 @@ +/********************************************************************** + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.net + * + * PostGIS is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 2 of the License, or + * (at your option) any later version. + * + * PostGIS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with PostGIS. If not, see . + * + ********************************************************************** + * + * Copyright 2015 Daniel Baston + * + **********************************************************************/ + +#include + +#include "liblwgeom_internal.h" +#include "lwgeom_log.h" + +static void +calc_distances_3d(const POINT3D* curr, const POINT3D* points, uint32_t npoints, double* distances) +{ + uint32_t i; + for (i = 0; i < npoints; i++) + { + distances[i] = distance3d_pt_pt(curr, &points[i]); + } +} + +static double +iterate_3d(POINT3D* curr, const POINT3D* points, uint32_t npoints, double* distances) +{ + uint32_t i; + POINT3D next = { 0, 0, 0 }; + double delta; + double denom = 0; + char hit = LW_FALSE; + + calc_distances_3d(curr, points, npoints, distances); + + for (i = 0; i < npoints; i++) + { + if (distances[i] == 0) + hit = LW_TRUE; + else + denom += 1.0 / distances[i]; + } + + for (i = 0; i < npoints; i++) + { + if (distances[i] > 0) + { + next.x += (points[i].x / distances[i]) / denom; + next.y += (points[i].y / distances[i]) / denom; + next.z += (points[i].z / distances[i]) / denom; + } + } + + /* If any of the intermediate points in the calculation is found in the + * set of input points, the standard Weiszfeld method gets stuck with a + * divide-by-zero. + * + * To get ourselves out of the hole, we follow an alternate procedure to + * get the next iteration, as described in: + * + * Vardi, Y. and Zhang, C. (2011) "A modified Weiszfeld algorithm for the + * Fermat-Weber location problem." Math. Program., Ser. A 90: 559-566. + * DOI 10.1007/s101070100222 + * + * Available online at the time of this writing at + * http://www.stat.rutgers.edu/home/cunhui/papers/43.pdf + */ + if (hit) + { + double dx = 0; + double dy = 0; + double dz = 0; + double r_inv; + POINT3D alt; + for (i = 0; i < npoints; i++) + { + if (distances[i] > 0) + { + dx += (points[i].x - curr->x) / distances[i]; + dy += (points[i].y - curr->y) / distances[i]; + dz += (points[i].y - curr->z) / distances[i]; + } + } + + r_inv = 1.0 / sqrt ( dx*dx + dy*dy + dz*dz ); + + alt.x = FP_MAX(0, 1.0 - r_inv)*next.x + FP_MIN(1.0, r_inv)*curr->x; + alt.y = FP_MAX(0, 1.0 - r_inv)*next.y + FP_MIN(1.0, r_inv)*curr->y; + alt.z = FP_MAX(0, 1.0 - r_inv)*next.z + FP_MIN(1.0, r_inv)*curr->z; + + next = alt; + } + + delta = distance3d_pt_pt(curr, &next); + + curr->x = next.x; + curr->y = next.y; + curr->z = next.z; + + return delta; +} + +static POINT3D +init_guess(const POINT3D* points, uint32_t npoints) +{ + POINT3D guess = { 0, 0, 0 }; + uint32_t i; + for (i = 0; i < npoints; i++) + { + guess.x += points[i].x / npoints; + guess.y += points[i].y / npoints; + guess.z += points[i].z / npoints; + } + + return guess; +} + +static POINT3D* +lwmpoint_extract_points_3d(const LWMPOINT* g, uint32_t* ngeoms) +{ + uint32_t i; + uint32_t n = 0; + int is_3d = lwgeom_has_z((LWGEOM*) g); + + POINT3D* points = lwalloc(g->ngeoms * sizeof(POINT3D)); + for (i = 0; i < g->ngeoms; i++) + { + LWGEOM* subg = lwcollection_getsubgeom((LWCOLLECTION*) g, i); + if (!lwgeom_is_empty(subg)) + { + getPoint3dz_p(((LWPOINT*) subg)->point, 0, (POINT3DZ*) &points[n++]); + if (!is_3d) + points[n-1].z = 0.0; /* in case the getPoint functions return NaN in the future for 2d */ + } + } + + if (ngeoms != NULL) + *ngeoms = n; + + return points; +} + +LWPOINT* +lwmpoint_median(const LWMPOINT* g, double tol, uint32_t max_iter, char fail_if_not_converged) +{ + uint32_t npoints; /* we need to count this ourselves so we can exclude empties */ + uint32_t i; + double delta = DBL_MAX; + double* distances; + POINT3D* points = lwmpoint_extract_points_3d(g, &npoints); + POINT3D median; + + if (npoints == 0) + { + lwfree(points); + return lwpoint_construct_empty(g->srid, 0, 0); + } + + median = init_guess(points, npoints); + + distances = lwalloc(npoints * sizeof(double)); + + for (i = 0; i < max_iter && delta > tol; i++) + { + delta = iterate_3d(&median, points, npoints, distances); + } + + lwfree(points); + lwfree(distances); + + if (fail_if_not_converged && delta > tol) + { + lwerror("Median failed to converge within %g after %d iterations.", tol, max_iter); + return NULL; + } + + if (lwgeom_has_z((LWGEOM*) g)) + { + return lwpoint_make3dz(g->srid, median.x, median.y, median.z); + } + else + { + return lwpoint_make2d(g->srid, median.x, median.y); + } +} + +LWPOINT* +lwgeom_median(const LWGEOM* g, double tol, uint32_t max_iter, char fail_if_not_converged) +{ + switch( lwgeom_get_type(g) ) + { + case POINTTYPE: + return lwpoint_clone(lwgeom_as_lwpoint(g)); + case MULTIPOINTTYPE: + return lwmpoint_median(lwgeom_as_lwmpoint(g), tol, max_iter, fail_if_not_converged); + default: + lwerror("Unsupported geometry type in lwgeom_median"); + return NULL; + } +} + diff --git a/regress/geometric_median.sql b/regress/geometric_median.sql new file mode 100644 index 000000000..f65ada0f7 --- /dev/null +++ b/regress/geometric_median.sql @@ -0,0 +1,16 @@ +-- postgres + +SELECT 't1', ST_GeometricMedian(null) IS NULL; +SELECT 't2', ST_AsText(ST_GeometricMedian('LINESTRING (1 1, 2 2)')); +SELECT 't3', 'POINT EMPTY' = ST_AsText(ST_GeometricMedian('POINT EMPTY')); +SELECT 't4', 'POINT EMPTY' = ST_AsText(ST_GeometricMedian('MULTIPOINT EMPTY')); +SELECT 't5', ST_SRID(ST_GeometricMedian('SRID=32611;POINT (1 1)')) = 32611; +SELECT 't6', ST_SRID(ST_GeometricMedian('SRID=32611;MULTIPOINT ((1 1), (2 7))')) = 32611; +SELECT 't7', ST_SRID(ST_GeometricMedian('SRID=32611;MULTIPOINT ((1 1))')) = 32611; +SELECT 't8', ST_Equals('POINTZ (15 15 15)', ST_AsText(ST_GeometricMedian('MULTIPOINT ((10 10 10), (10 20 10), (20 10 10), (20 20 10), (10 10 20), (10 20 20), (20 10 20), (20 20 20))'))); +-- Doesn't fail even though we don't let it iterate +SELECT 't9', ST_IsEmpty(ST_GeometricMedian('MULTIPOINT (0 0, 1 1, 0 1, 2 2)', max_iter := 1)); +-- Unless we enfore that it must converge. +SELECT 't10', ST_IsEmpty(ST_GeometricMedian('MULTIPOINT (0 0, 1 1, 0 1, 2 2)', max_iter := 1, fail_if_not_converged := true)); +-- But if we drop the tolerance, it's OK +SELECT 't11', ST_IsEmpty(ST_GeometricMedian('MULTIPOINT (0 0, 1 1, 0 1, 2 2)', max_iter := 1, tolerance := 0.1, fail_if_not_converged := true)); diff --git a/regress/geometric_median_expected b/regress/geometric_median_expected new file mode 100644 index 000000000..ba2e48ff5 --- /dev/null +++ b/regress/geometric_median_expected @@ -0,0 +1,11 @@ +t1|t +ERROR: Unsupported geometry type in lwgeom_median +t3|t +t4|t +t5|t +t6|t +t7|t +t8|t +t9|f +ERROR: Median failed to converge within 2e-6 after 1 iterations. +t11|f