[postgis-tickets] r17821 - Speed up ST_Simplify
Raul
raul at rmr.ninja
Tue Sep 24 03:58:30 PDT 2019
Author: algunenano
Date: 2019-09-24 03:58:30 -0700 (Tue, 24 Sep 2019)
New Revision: 17821
Modified:
trunk/NEWS
trunk/liblwgeom/cunit/cu_misc.c
trunk/liblwgeom/liblwgeom.h.in
trunk/liblwgeom/liblwgeom_internal.h
trunk/liblwgeom/lwgeom.c
trunk/liblwgeom/lwgeom_geos_split.c
trunk/liblwgeom/measures.c
trunk/liblwgeom/ptarray.c
trunk/postgis/lwgeom_functions_analytic.c
Log:
Speed up ST_Simplify
Closes #4510
Closes https://github.com/postgis/postgis/pull/480
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/NEWS 2019-09-24 10:58:30 UTC (rev 17821)
@@ -16,6 +16,7 @@
- #4505, Speed up conversion of geometries to/from GEOS (Dan Baston)
- #4507, Use GEOSMakeValid and GEOSBuildArea for GEOS 3.8+ (Dan Baston)
- #4491, Speed up ST_RemoveRepeatedPoints (Raúl Marín)
+ - #4510, Speed up ST_Simplify (Raúl Marín)
PostGIS 3.0.0alpha4
2019/08/10
Modified: trunk/liblwgeom/cunit/cu_misc.c
===================================================================
--- trunk/liblwgeom/cunit/cu_misc.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/cunit/cu_misc.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -62,6 +62,14 @@
lwgeom_free(geom);
lwgeom_free(geom2d);
lwfree(wkt_out);
+
+ geom = lwgeom_from_wkt("POLYGON((0 0,1 1,1 3,0 4,-2 3,-1 1,0 0))", LW_PARSER_CHECK_NONE);
+ geom2d = lwgeom_simplify(geom, 1, LW_FALSE);
+ wkt_out = lwgeom_to_ewkt(geom2d);
+ CU_ASSERT_STRING_EQUAL("POLYGON((0 0,0 4,-2 3,0 0))", wkt_out);
+ lwgeom_free(geom);
+ lwgeom_free(geom2d);
+ lwfree(wkt_out);
}
static void test_misc_count_vertices(void)
Modified: trunk/liblwgeom/liblwgeom.h.in
===================================================================
--- trunk/liblwgeom/liblwgeom.h.in 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/liblwgeom.h.in 2019-09-24 10:58:30 UTC (rev 17821)
@@ -1265,7 +1265,6 @@
/* general utilities 2D */
extern double distance2d_pt_pt(const POINT2D *p1, const POINT2D *p2);
-extern double distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B);
extern double distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B);
extern LWGEOM* lwgeom_closest_line(const LWGEOM *lw1, const LWGEOM *lw2);
extern LWGEOM* lwgeom_furthest_line(const LWGEOM *lw1, const LWGEOM *lw2);
@@ -1362,7 +1361,7 @@
extern void lwgeom_reverse_in_place(LWGEOM *lwgeom);
extern void lwgeom_force_clockwise(LWGEOM *lwgeom);
extern void lwgeom_longitude_shift(LWGEOM *lwgeom);
-extern void lwgeom_simplify_in_place(LWGEOM *igeom, double dist, int preserve_collapsed);
+extern int lwgeom_simplify_in_place(LWGEOM *igeom, double dist, int preserve_collapsed);
extern void lwgeom_affine(LWGEOM *geom, const AFFINE *affine);
extern void lwgeom_scale(LWGEOM *geom, const POINT4D *factors);
extern int lwgeom_remove_repeated_points_in_place(LWGEOM *in, double tolerance);
Modified: trunk/liblwgeom/liblwgeom_internal.h
===================================================================
--- trunk/liblwgeom/liblwgeom_internal.h 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/liblwgeom_internal.h 2019-09-24 10:58:30 UTC (rev 17821)
@@ -210,7 +210,7 @@
/**
* @param minpts minimum number of points to retain, if possible.
*/
-void ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, uint32_t minpts);
+void ptarray_simplify_in_place(POINTARRAY *pa, double tolerance, uint32_t minpts);
/*
* The possible ways a pair of segments can interact. Returned by lw_segment_intersects
Modified: trunk/liblwgeom/lwgeom.c
===================================================================
--- trunk/liblwgeom/lwgeom.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/lwgeom.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -1711,18 +1711,19 @@
/**************************************************************/
-void
+int
lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
{
+ int modified = LW_FALSE;
switch (geom->type)
{
/* No-op! Cannot simplify points or triangles */
case POINTTYPE:
- return;
+ return modified;
case TRIANGLETYPE:
{
if (preserve_collapsed)
- return;
+ return modified;
LWTRIANGLE *t = lwgeom_as_lwtriangle(geom);
POINTARRAY *pa = t->points;
ptarray_simplify_in_place(pa, epsilon, 0);
@@ -1729,13 +1730,17 @@
if (pa->npoints < 3)
{
pa->npoints = 0;
+ modified = LW_TRUE;
}
+ break;
}
case LINETYPE:
{
LWLINE *g = (LWLINE*)(geom);
POINTARRAY *pa = g->points;
+ uint32_t in_npoints = pa->npoints;
ptarray_simplify_in_place(pa, epsilon, 2);
+ modified = in_npoints != pa->npoints;
/* Invalid output */
if (pa->npoints == 1 && pa->maxpoints > 1)
{
@@ -1771,12 +1776,20 @@
/* Skip zero'ed out rings */
if(!pa)
continue;
+ uint32_t in_npoints = pa->npoints;
ptarray_simplify_in_place(pa, epsilon, minpoints);
+ modified |= in_npoints != pa->npoints;
/* Drop collapsed rings */
if(pa->npoints < 4)
{
- ptarray_free(pa);
- continue;
+ /* Any ring deeper than this one is, by OGR definition, smaller
+ * so we drop them all */
+ for (; i < g->nrings; i++)
+ {
+ pa = g->rings[i];
+ ptarray_free(pa);
+ }
+ break;
}
g->rings[j++] = pa;
}
@@ -1797,7 +1810,7 @@
{
LWGEOM *g = col->geoms[i];
if (!g) continue;
- lwgeom_simplify_in_place(g, epsilon, preserve_collapsed);
+ modified |= lwgeom_simplify_in_place(g, epsilon, preserve_collapsed);
/* Drop zero'ed out geometries */
if(lwgeom_is_empty(g))
{
@@ -1816,7 +1829,12 @@
break;
}
}
- return;
+
+ if (modified)
+ {
+ lwgeom_drop_bbox(geom);
+ }
+ return modified;
}
LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
Modified: trunk/liblwgeom/lwgeom_geos_split.c
===================================================================
--- trunk/liblwgeom/lwgeom_geos_split.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/lwgeom_geos_split.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -209,7 +209,7 @@
lwline_split_by_point_to(const LWLINE* lwline_in, const LWPOINT* blade_in,
LWMLINE* v)
{
- double mindist = -1;
+ double mindist_sqr = -1;
POINT4D pt, pt_projected;
POINT4D p1, p2;
POINTARRAY *ipa = lwline_in->points;
@@ -238,23 +238,28 @@
for ( i = 0; i < nsegs; i++ )
{
getPoint4d_p(ipa, i+1, &p2);
- double dist;
- dist = distance2d_pt_seg((POINT2D*)&pt, (POINT2D*)&p1, (POINT2D*)&p2);
- LWDEBUGF(4, " Distance of point %g %g to segment %g %g, %g %g: %g", pt.x, pt.y, p1.x, p1.y, p2.x, p2.y, dist);
- if (i==0 || dist < mindist )
+ double dist_sqr = distance2d_sqr_pt_seg((POINT2D *)&pt, (POINT2D *)&p1, (POINT2D *)&p2);
+ LWDEBUGF(4, "Distance (squared) of point %g %g to segment %g %g, %g %g: %g",
+ pt.x, pt.y,
+ p1.x, p1.y,
+ p2.x, p2.y,
+ dist_sqr);
+ if (i == 0 || dist_sqr < mindist_sqr)
{
- mindist = dist;
+ mindist_sqr = dist_sqr;
seg=i;
- if ( mindist == 0.0 ) break; /* can't be closer than ON line */
+ if (mindist_sqr == 0.0)
+ break; /* can't be closer than ON line */
}
p1 = p2;
}
LWDEBUGF(3, "Closest segment: %d", seg);
- LWDEBUGF(3, "mindist: %g", mindist);
+ LWDEBUGF(3, "mindist: %g", mindist_sqr);
/* No intersection */
- if ( mindist > 0 ) return 0;
+ if (mindist_sqr > 0)
+ return 0;
/* empty or single-point line, intersection on boundary */
if ( seg == UINT32_MAX ) return 1;
Modified: trunk/liblwgeom/measures.c
===================================================================
--- trunk/liblwgeom/measures.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/measures.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -2403,22 +2403,21 @@
return hypot(hside, vside);
}
+/* return distance squared, useful to avoid sqrt calculations */
double
-distance2d_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
+distance2d_sqr_pt_seg(const POINT2D *C, const POINT2D *A, const POINT2D *B)
{
- double r, s;
-
/*if start==end, then use pt distance */
if ((A->x == B->x) && (A->y == B->y))
- return distance2d_pt_pt(p, A);
+ return distance2d_sqr_pt_pt(C, A);
/*
* otherwise, we use comp.graphics.algorithms
* Frequently Asked Questions method
*
- * (1) AC dot AB
+ * (1) AC dot AB
* r = ---------
- * ||AB||^2
+ * ||AB||^2
* r has the following meaning:
* r=0 P = A
* r=1 P = B
@@ -2427,13 +2426,17 @@
* 0<r<1 P is interior to AB
*/
- r = ((p->x - A->x) * (B->x - A->x) + (p->y - A->y) * (B->y - A->y)) /
- ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
+ double ba_x = (B->x - A->x);
+ double ba_y = (B->y - A->y);
+ double ab_length_sqr = (ba_x * ba_x + ba_y * ba_y);
+ double ca_x = (C->x - A->x);
+ double ca_y = (C->y - A->y);
+ double dot_ac_ab = (ca_x * ba_x + ca_y * ba_y);
- if (r < 0)
- return distance2d_pt_pt(p, A);
- if (r > 1)
- return distance2d_pt_pt(p, B);
+ if (dot_ac_ab <= 0)
+ return distance2d_sqr_pt_pt(C, A);
+ if (dot_ac_ab >= ab_length_sqr)
+ return distance2d_sqr_pt_pt(C, B);
/*
* (2)
@@ -2445,45 +2448,12 @@
*
*/
- s = ((A->y - p->y) * (B->x - A->x) - (A->x - p->x) * (B->y - A->y)) /
- ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
+ double s_numerator = ca_x * ba_y - ca_y * ba_x;
- return FP_ABS(s) * sqrt((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
+ /* Distance = (s_num / ab) * (s_num / ab) * ab == s_num * s_num / ab) */
+ return s_numerator * s_numerator / ab_length_sqr;
}
-/* return distance squared, useful to avoid sqrt calculations */
-double
-distance2d_sqr_pt_seg(const POINT2D *p, const POINT2D *A, const POINT2D *B)
-{
- double r, s;
-
- if ((A->x == B->x) && (A->y == B->y))
- return distance2d_sqr_pt_pt(p, A);
-
- r = ((p->x - A->x) * (B->x - A->x) + (p->y - A->y) * (B->y - A->y)) /
- ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
-
- if (r < 0)
- return distance2d_sqr_pt_pt(p, A);
- if (r > 1)
- return distance2d_sqr_pt_pt(p, B);
-
- /*
- * (2)
- * (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
- * s = -----------------------------
- * L^2
- *
- * Then the distance from C to P = |s|*L.
- *
- */
-
- s = ((A->y - p->y) * (B->x - A->x) - (A->x - p->x) * (B->y - A->y)) /
- ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
-
- return s * s * ((B->x - A->x) * (B->x - A->x) + (B->y - A->y) * (B->y - A->y));
-}
-
/**
* Compute the azimuth of segment AB in radians.
* Return 0 on exception (same point), 1 otherwise.
Modified: trunk/liblwgeom/ptarray.c
===================================================================
--- trunk/liblwgeom/ptarray.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/liblwgeom/ptarray.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -1331,13 +1331,13 @@
/* Loop through pointarray looking for nearest segment */
for (t=1; t<pa->npoints; t++)
{
- double dist;
+ double dist_sqr;
end = getPoint2d_cp(pa, t);
- dist = distance2d_pt_seg(&p, start, end);
+ dist_sqr = distance2d_sqr_pt_seg(&p, start, end);
- if ( dist < mindist )
+ if (dist_sqr < mindist)
{
- mindist = dist;
+ mindist = dist_sqr;
seg=t-1;
if ( mindist == 0 )
{
@@ -1348,6 +1348,7 @@
start = end;
}
+ mindist = sqrt(mindist);
if ( mindistout ) *mindistout = mindist;
@@ -1514,143 +1515,151 @@
return;
}
-
-/************************************************************************/
-
-static void
-ptarray_dp_findsplit_in_place(const POINTARRAY *pts, int p1, int p2, int *split, double *dist)
+/* Out of the points in pa [itfist .. itlast], finds the one that's farthest away from
+ * the segment determined by pts[itfist] and pts[itlast].
+ * Returns itfirst if no point was found futher away than max_distance_sqr
+ */
+static uint32_t
+ptarray_dp_findsplit_in_place(const POINTARRAY *pts, uint32_t it_first, uint32_t it_last, double max_distance_sqr)
{
- int k;
- const POINT2D *pk, *pa, *pb;
- double tmp, d;
+ uint32_t split = it_first;
+ if ((it_first - it_last) < 2)
+ return it_first;
- LWDEBUG(4, "function called");
+ const POINT2D *A = getPoint2d_cp(pts, it_first);
+ const POINT2D *B = getPoint2d_cp(pts, it_last);
- *split = p1;
- d = -1;
-
- if (p1 + 1 < p2)
+ if (distance2d_sqr_pt_pt(A, B) < DBL_EPSILON)
{
-
- pa = getPoint2d_cp(pts, p1);
- pb = getPoint2d_cp(pts, p2);
-
- LWDEBUGF(4, "P%d(%f,%f) to P%d(%f,%f)",
- p1, pa->x, pa->y, p2, pb->x, pb->y);
-
- for (k=p1+1; k<p2; k++)
+ /* If p1 == p2, we can just calculate the distance from each point to A */
+ for (uint32_t itk = it_first + 1; itk < it_last; itk++)
{
- pk = getPoint2d_cp(pts, k);
-
- LWDEBUGF(4, "P%d(%f,%f)", k, pk->x, pk->y);
-
- /* distance computation */
- tmp = distance2d_sqr_pt_seg(pk, pa, pb);
-
- if (tmp > d)
+ const POINT2D *pk = getPoint2d_cp(pts, itk);
+ double distance_sqr = distance2d_sqr_pt_pt(pk, A);
+ if (distance_sqr > max_distance_sqr)
{
- d = tmp; /* record the maximum */
- *split = k;
-
- LWDEBUGF(4, "P%d is farthest (%g)", k, d);
+ split = itk;
+ max_distance_sqr = distance_sqr;
}
}
- *dist = d;
+ return split;
}
- else
+
+ /* This is based on distance2d_sqr_pt_seg, but heavily inlined here to avoid recalculations */
+ double ba_x = (B->x - A->x);
+ double ba_y = (B->y - A->y);
+ double ab_length_sqr = (ba_x * ba_x + ba_y * ba_y);
+ /* To avoid the division by ab_length_sqr in the 3rd path, we normalize here
+ * and multiply in the first two paths [(dot_ac_ab < 0) and (> ab_length_sqr)] */
+ max_distance_sqr *= ab_length_sqr;
+ for (uint32_t itk = it_first + 1; itk < it_last; itk++)
{
- LWDEBUG(3, "segment too short, no split/no dist");
- *dist = -1;
+ const POINT2D *C = getPoint2d_cp(pts, itk);
+ double distance_sqr;
+ double ca_x = (C->x - A->x);
+ double ca_y = (C->y - A->y);
+ double dot_ac_ab = (ca_x * ba_x + ca_y * ba_y);
+
+ if (dot_ac_ab <= 0.0)
+ {
+ distance_sqr = distance2d_sqr_pt_pt(C, A) * ab_length_sqr;
+ }
+ else if (dot_ac_ab >= ab_length_sqr)
+ {
+ distance_sqr = distance2d_sqr_pt_pt(C, B) * ab_length_sqr;
+ }
+ else
+ {
+ double s_numerator = ca_x * ba_y - ca_y * ba_x;
+ distance_sqr = s_numerator * s_numerator; /* Missing division by ab_length_sqr on purpose */
+ }
+
+ if (distance_sqr > max_distance_sqr)
+ {
+ split = itk;
+ max_distance_sqr = distance_sqr;
+ }
}
+ return split;
}
-static int
-int_cmp(const void *a, const void *b)
-{
- /* casting pointer types */
- const int *ia = (const int *)a;
- const int *ib = (const int *)b;
- /* returns negative if b > a and positive if a > b */
- return *ia - *ib;
-}
-
void
-ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, uint32_t minpts)
+ptarray_simplify_in_place(POINTARRAY *pa, double tolerance, uint32_t minpts)
{
- static size_t stack_size = 256;
- int *stack, *outlist; /* recursion stack */
- int stack_static[stack_size];
- int outlist_static[stack_size];
- int sp = -1; /* recursion stack pointer */
- int p1, split;
- uint32_t outn = 0;
- int pai = 0;
- uint32_t i;
- double dist;
- double eps_sqr = epsilon * epsilon;
-
/* Do not try to simplify really short things */
- if (pa->npoints < 3) return;
+ if (pa->npoints < 3 || pa->npoints <= minpts)
+ return;
- /* Only heap allocate book-keeping arrays if necessary */
- if (pa->npoints > stack_size)
- {
- stack = lwalloc(sizeof(int) * pa->npoints);
- outlist = lwalloc(sizeof(int) * pa->npoints);
- }
- else
- {
- stack = stack_static;
- outlist = outlist_static;
- }
+ /* We use this array to keep track of the points we are keeping, so
+ * we store just TRUE / FALSE in their position */
+ uint8_t *kept_points = lwalloc(sizeof(uint8_t) * pa->npoints);
+ memset(kept_points, LW_FALSE, sizeof(uint8_t) * pa->npoints);
+ kept_points[0] = LW_TRUE;
+ kept_points[pa->npoints - 1] = LW_TRUE;
+ uint32_t keptn = 2;
- p1 = 0;
- stack[++sp] = pa->npoints-1;
+ /* We use this array as a stack to store the iterators that we are going to need
+ * in the following steps.
+ * This is ~10% faster than iterating over @kept_points looking for them
+ */
+ uint32_t *iterator_stack = lwalloc(sizeof(uint32_t) * pa->npoints);
+ iterator_stack[0] = 0;
+ uint32_t iterator_stack_size = 1;
- /* Add first point to output list */
- outlist[outn++] = 0;
- do
+ uint32_t it_first = 0;
+ uint32_t it_last = pa->npoints - 1;
+
+ const double tolerance_sqr = tolerance * tolerance;
+ /* For the first @minpts points we ignore the tolerance */
+ double it_tol = keptn >= minpts ? tolerance_sqr : -1.0;
+
+ while (iterator_stack_size)
{
- ptarray_dp_findsplit_in_place(pa, p1, stack[sp], &split, &dist);
-
- if ((dist > eps_sqr) || ((outn + sp+1 < minpts) && (dist >= 0)))
+ uint32_t split = ptarray_dp_findsplit_in_place(pa, it_first, it_last, it_tol);
+ if (split == it_first)
{
- stack[++sp] = split;
+ it_first = it_last;
+ it_last = iterator_stack[--iterator_stack_size];
}
else
{
- outlist[outn++] = stack[sp];
- p1 = stack[sp--];
+ kept_points[split] = LW_TRUE;
+ keptn++;
+
+ iterator_stack[iterator_stack_size++] = it_last;
+ it_last = split;
+ it_tol = keptn >= minpts ? tolerance_sqr : -1.0;
}
}
- while (!(sp<0));
- /* Put list of retained points into order */
- qsort(outlist, outn, sizeof(int), int_cmp);
- /* Copy retained points to front of array */
- for (i = 0; i < outn; i++)
+ const size_t pt_size = ptarray_point_size(pa);
+ /* The first point is already in place, so we don't need to copy it */
+ size_t kept_it = 1;
+ if (keptn == 2)
{
- int j = outlist[i];
- /* Indexes the same, means no copy required */
- if (j == pai)
+ /* If there are 2 points remaining, it has to be first and last as
+ * we added those at the start */
+ memcpy(pa->serialized_pointlist + pt_size * kept_it,
+ pa->serialized_pointlist + pt_size * (pa->npoints - 1),
+ pt_size);
+ }
+ else
+ {
+ for (uint32_t i = 1; i < pa->npoints; i++)
{
- pai++;
- continue;
+ if (kept_points[i])
+ {
+ memcpy(pa->serialized_pointlist + pt_size * kept_it,
+ pa->serialized_pointlist + pt_size * i,
+ pt_size);
+ kept_it++;
+ }
}
- /* Indexes different, copy value down */
- ptarray_copy_point(pa, j, pai++);
}
+ pa->npoints = keptn;
- /* Adjust point count on array */
- pa->npoints = outn;
-
- /* Only free if arrays are on heap */
- if (stack != stack_static)
- lwfree(stack);
- if (outlist != outlist_static)
- lwfree(outlist);
-
- return;
+ lwfree(kept_points);
+ lwfree(iterator_stack);
}
/************************************************************************/
Modified: trunk/postgis/lwgeom_functions_analytic.c
===================================================================
--- trunk/postgis/lwgeom_functions_analytic.c 2019-09-24 10:45:43 UTC (rev 17820)
+++ trunk/postgis/lwgeom_functions_analytic.c 2019-09-24 10:58:30 UTC (rev 17821)
@@ -63,33 +63,32 @@
PG_FUNCTION_INFO_V1(LWGEOM_simplify2d);
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS)
{
- GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P(0);
+ GSERIALIZED *geom = PG_GETARG_GSERIALIZED_P_COPY(0);
double dist = PG_GETARG_FLOAT8(1);
GSERIALIZED *result;
int type = gserialized_get_type(geom);
- LWGEOM *in, *out;
+ LWGEOM *in;
bool preserve_collapsed = false;
+ int modified = LW_FALSE;
+ /* Can't simplify points! */
+ if ( type == POINTTYPE || type == MULTIPOINTTYPE )
+ PG_RETURN_POINTER(geom);
+
/* Handle optional argument to preserve collapsed features */
if ((PG_NARGS() > 2) && (!PG_ARGISNULL(2)))
preserve_collapsed = PG_GETARG_BOOL(2);
- /* Can't simplify points! */
- if ( type == POINTTYPE || type == MULTIPOINTTYPE )
- PG_RETURN_POINTER(geom);
-
in = lwgeom_from_gserialized(geom);
- out = lwgeom_simplify(in, dist, preserve_collapsed);
- if ( ! out ) PG_RETURN_NULL();
+ modified = lwgeom_simplify_in_place(in, dist, preserve_collapsed);
+ if (!modified)
+ PG_RETURN_POINTER(geom);
- /* COMPUTE_BBOX TAINTING */
- if (in->bbox)
- lwgeom_refresh_bbox(out);
+ if (!in || lwgeom_is_empty(in))
+ PG_RETURN_NULL();
- result = geometry_serialize(out);
- lwgeom_free(out);
- PG_FREE_IF_COPY(geom, 0);
+ result = geometry_serialize(in);
PG_RETURN_POINTER(result);
}
More information about the postgis-tickets
mailing list