static void test_lwgeom_simplify(void)
{
- LWGEOM *l;
- LWGEOM *g;
- char *ewkt;
-
- /* Simplify but only so far... */
- g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 10, LW_TRUE);
- ewkt = lwgeom_to_ewkt(l);
- CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
- lwgeom_free(g);
- lwgeom_free(l);
- lwfree(ewkt);
-
- /* Simplify but only so far... */
- g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 10, LW_TRUE);
- ewkt = lwgeom_to_ewkt(l);
- CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
- lwgeom_free(g);
- lwgeom_free(l);
- lwfree(ewkt);
-
- /* Simplify and collapse */
- g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 10, LW_FALSE);
- CU_ASSERT_EQUAL(l, NULL);
- lwgeom_free(g);
- lwgeom_free(l);
-
- /* Simplify and collapse */
- g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 10, LW_FALSE);
- CU_ASSERT_EQUAL(l, NULL);
- lwgeom_free(g);
- lwgeom_free(l);
-
- /* Not simplifiable */
- g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 1.0, LW_FALSE);
- ewkt = lwgeom_to_ewkt(l);
- CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
- lwgeom_free(g);
- lwgeom_free(l);
- lwfree(ewkt);
-
- /* Simplifiable */
- g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
- l = lwgeom_simplify(g, 1.0, LW_FALSE);
- ewkt = lwgeom_to_ewkt(l);
- CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
- lwgeom_free(g);
- lwgeom_free(l);
- lwfree(ewkt);
+ LWGEOM *l;
+ LWGEOM *g;
+ char *ewkt;
+
+ /* Simplify but only so far... */
+ g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 10, LW_TRUE);
+ ewkt = lwgeom_to_ewkt(l);
+ CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,0 0)");
+ lwgeom_free(g);
+ lwgeom_free(l);
+ lwfree(ewkt);
+
+ /* Simplify but only so far... */
+ g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 10, LW_TRUE);
+ ewkt = lwgeom_to_ewkt(l);
+ CU_ASSERT_STRING_EQUAL(ewkt, "POLYGON((0 0,1 0,1 1,0 0))");
+ lwgeom_free(g);
+ lwgeom_free(l);
+ lwfree(ewkt);
+
+ /* Simplify and collapse */
+ g = lwgeom_from_wkt("LINESTRING(0 0, 1 0, 1 1, 0 1, 0 0)", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 10, LW_FALSE);
+ CU_ASSERT_EQUAL(l, NULL);
+ lwgeom_free(g);
+ lwgeom_free(l);
+
+ /* Simplify and collapse */
+ g = lwgeom_from_wkt("POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 10, LW_FALSE);
+ CU_ASSERT_EQUAL(l, NULL);
+ lwgeom_free(g);
+ lwgeom_free(l);
+
+ /* Not simplifiable */
+ g = lwgeom_from_wkt("LINESTRING(0 0, 50 1.00001, 100 0)", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 1.0, LW_FALSE);
+ ewkt = lwgeom_to_ewkt(l);
+ CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,50 1.00001,100 0)");
+ lwgeom_free(g);
+ lwgeom_free(l);
+ lwfree(ewkt);
+
+ /* Simplifiable */
+ g = lwgeom_from_wkt("LINESTRING(0 0,50 0.99999,100 0)", LW_PARSER_CHECK_NONE);
+ l = lwgeom_simplify(g, 1.0, LW_FALSE);
+ ewkt = lwgeom_to_ewkt(l);
+ CU_ASSERT_STRING_EQUAL(ewkt, "LINESTRING(0 0,100 0)");
+ lwgeom_free(g);
+ lwgeom_free(l);
+ lwfree(ewkt);
}
+extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed);
+
/****************************************************************
* READ/WRITE FUNCTIONS
*
extern void lwgeom_reverse_in_place(LWGEOM *lwgeom);
extern void lwgeom_force_clockwise(LWGEOM *lwgeom);
extern void lwgeom_longitude_shift(LWGEOM *lwgeom);
-extern LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed);
+extern void 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 LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance);
-/****************************************************************/
-
-
/**
* @brief wrap geometry on given cut x value
*
/**
* Remove repeated points!
*/
-extern LWGEOM* lwgeom_remove_repeated_points(const LWGEOM *in, double tolerance);
/**
* Swap ordinate values in every vertex of the geometry.
/**
* @param minpts minimun number of points to retain, if possible.
*/
-POINTARRAY* ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts);
-LWLINE* lwline_simplify(const LWLINE *iline, double dist, int preserve_collapsed);
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed);
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed);
+void ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts);
/*
* The possible ways a pair of segments can interact. Returned by lw_segment_intersects
*/
void ptarray_longitude_shift(POINTARRAY *pa);
+/*
+* Support for in place modification of point arrays, fast
+* function to move coordinate values around
+*/
+void ptarray_copy_point(POINTARRAY *pa, int from, int to);
+
/*
* Reverse
*/
return v;
}
-LWCOLLECTION* lwcollection_simplify(const LWCOLLECTION *igeom, double dist, int preserve_collapsed)
-{
- int i;
- LWCOLLECTION *out = lwcollection_construct_empty(igeom->type, igeom->srid, FLAGS_GET_Z(igeom->flags), FLAGS_GET_M(igeom->flags));
-
- if( lwcollection_is_empty(igeom) )
- return out; /* should we return NULL instead ? */
-
- for( i = 0; i < igeom->ngeoms; i++ )
- {
- LWGEOM *ngeom = lwgeom_simplify(igeom->geoms[i], dist, preserve_collapsed);
- if ( ngeom ) out = lwcollection_add_lwgeom(out, ngeom);
- }
-
- return out;
-}
int lwcollection_allows_subtype(int collectiontype, int subtype)
{
}
}
-LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
+/**************************************************************/
+
+
+void
+lwgeom_simplify_in_place(LWGEOM *geom, double epsilon, int preserve_collapsed)
{
- switch (igeom->type)
+ switch (geom->type)
{
- case POINTTYPE:
- case MULTIPOINTTYPE:
- return lwgeom_clone(igeom);
- case LINETYPE:
- return (LWGEOM*)lwline_simplify((LWLINE*)igeom, dist, preserve_collapsed);
- case POLYGONTYPE:
- return (LWGEOM*)lwpoly_simplify((LWPOLY*)igeom, dist, preserve_collapsed);
- case MULTILINETYPE:
- case MULTIPOLYGONTYPE:
- case COLLECTIONTYPE:
- return (LWGEOM*)lwcollection_simplify((LWCOLLECTION *)igeom, dist, preserve_collapsed);
- default:
- lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(igeom->type));
+ /* No-op! Cannot simplify points */
+ case POINTTYPE:
+ return;
+ case LINETYPE:
+ {
+ LWLINE *g = (LWLINE*)(geom);
+ POINTARRAY *pa = g->points;
+ ptarray_simplify_in_place(pa, epsilon, 2);
+ /* Invalid output */
+ if (pa->npoints == 1 && pa->maxpoints > 1)
+ {
+ /* Use first point as last point */
+ if (preserve_collapsed)
+ {
+ pa->npoints = 2;
+ ptarray_copy_point(pa, 0, 1);
+ }
+ /* Finish the collapse process */
+ else
+ {
+ pa->npoints = 0;
+ }
+ }
+ /* Duped output, force collapse */
+ if (pa->npoints == 2 && !preserve_collapsed)
+ {
+ if (p2d_same(getPoint2d_cp(pa, 0), getPoint2d_cp(pa, 1)))
+ pa->npoints = 0;
+ }
+ break;
+ }
+ case POLYGONTYPE:
+ {
+ int i, j = 0;
+ LWPOLY *g = (LWPOLY*)(geom);
+ for (i = 0; i < g->nrings; i++)
+ {
+ POINTARRAY *pa = g->rings[i];
+ /* Only stop collapse on first ring */
+ int minpoints = (preserve_collapsed && i == 0) ? 4 : 0;
+ /* Skip zero'ed out rings */
+ if(!pa)
+ continue;
+ ptarray_simplify_in_place(pa, epsilon, minpoints);
+ /* Drop collapsed rings */
+ if(pa->npoints < 4)
+ {
+ ptarray_free(pa);
+ continue;
+ }
+ g->rings[j++] = pa;
+ }
+ /* Update ring count */
+ g->nrings = j;
+ break;
+ }
+ /* Can process all multi* types as generic collection */
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ {
+ int i, j = 0;
+ LWCOLLECTION *col = (LWCOLLECTION*)geom;
+ for (i = 0; i < col->ngeoms; i++)
+ {
+ LWGEOM *g = col->geoms[i];
+ if (!g) continue;
+ lwgeom_simplify_in_place(g, epsilon, preserve_collapsed);
+ /* Drop zero'ed out geometries */
+ if(lwgeom_is_empty(g))
+ {
+ lwgeom_free(g);
+ continue;
+ }
+ col->geoms[j++] = g;
+ }
+ /* Update geometry count */
+ col->ngeoms = j;
+ break;
+ }
+ default:
+ {
+ lwerror("%s: unsupported geometry type: %s", __func__, lwtype_name(geom->type));
+ break;
+ }
}
- return NULL;
+ return;
}
+
+/**************************************************************/
+
+LWGEOM* lwgeom_simplify(const LWGEOM *igeom, double dist, int preserve_collapsed)
+{
+ LWGEOM *lwgeom_out = lwgeom_clone_deep(igeom);
+ lwgeom_simplify_in_place(lwgeom_out, dist, preserve_collapsed);
+ if (lwgeom_is_empty(lwgeom_out)) return NULL;
+ return lwgeom_out;
+}
+
+
double lwgeom_area(const LWGEOM *geom)
{
int type = geom->type;
}
}
+void
+ptarray_copy_point(POINTARRAY *pa, int from, int to)
+{
+ int ndims = FLAGS_NDIMS(pa->flags);
+ switch (ndims)
+ {
+ case 2:
+ {
+ POINT2D *p_from = (POINT2D*)(getPoint_internal(pa, from));
+ POINT2D *p_to = (POINT2D*)(getPoint_internal(pa, to));
+ *p_to = *p_from;
+ return;
+ }
+ case 3:
+ {
+ POINT3D *p_from = (POINT3D*)(getPoint_internal(pa, from));
+ POINT3D *p_to = (POINT3D*)(getPoint_internal(pa, to));
+ *p_to = *p_from;
+ return;
+ }
+ case 4:
+ {
+ POINT4D *p_from = (POINT4D*)(getPoint_internal(pa, from));
+ POINT4D *p_to = (POINT4D*)(getPoint_internal(pa, to));
+ *p_to = *p_from;
+ return;
+ }
+ default:
+ {
+ lwerror("%s: unsupported number of dimensions - %d", __func__, ndims);
+ return;
+ }
+ }
+ return;
+}
+
/************************************************
* debugging routines
return line->points->npoints;
}
-LWLINE* lwline_simplify(const LWLINE *iline, double dist, int preserve_collapsed)
-{
- static const int minvertices = 2; /* TODO: allow setting this */
- LWLINE *oline;
- POINTARRAY *pa;
-
- LWDEBUG(2, "function called");
-
- /* Skip empty case */
- if( lwline_is_empty(iline) )
- return NULL;
-
- pa = ptarray_simplify(iline->points, dist, minvertices);
- if ( ! pa ) return NULL;
-
- /* Make sure single-point collapses have two points */
- if ( pa->npoints == 1 )
- {
- /* Make sure single-point collapses have two points */
- if ( preserve_collapsed )
- {
- POINT4D pt;
- getPoint4d_p(pa, 0, &pt);
- ptarray_append_point(pa, &pt, LW_TRUE);
- }
- /* Return null for collapse */
- else
- {
- ptarray_free(pa);
- return NULL;
- }
- }
-
- oline = lwline_construct(iline->srid, NULL, pa);
- oline->type = iline->type;
- return oline;
-}
-
double lwline_length(const LWLINE *line)
{
if ( lwline_is_empty(line) )
return v;
}
-LWPOLY* lwpoly_simplify(const LWPOLY *ipoly, double dist, int preserve_collapsed)
-{
- int i;
- LWPOLY *opoly = lwpoly_construct_empty(ipoly->srid, FLAGS_GET_Z(ipoly->flags), FLAGS_GET_M(ipoly->flags));
-
- LWDEBUGF(2, "%s: simplifying polygon with %d rings", __func__, ipoly->nrings);
-
- if ( lwpoly_is_empty(ipoly) )
- {
- lwpoly_free(opoly);
- return NULL;
- }
-
- for ( i = 0; i < ipoly->nrings; i++ )
- {
- POINTARRAY *opts;
- int minvertices = 0;
-
- /* We'll still let holes collapse, but if we're preserving */
- /* and this is a shell, we ensure it is kept */
- if ( preserve_collapsed && i == 0 )
- minvertices = 4;
-
- opts = ptarray_simplify(ipoly->rings[i], dist, minvertices);
-
- LWDEBUGF(3, "ring%d simplified from %d to %d points", i, ipoly->rings[i]->npoints, opts->npoints);
-
- /* Less points than are needed to form a closed ring, we can't use this */
- if ( opts->npoints < 4 )
- {
- LWDEBUGF(3, "ring%d skipped (% pts)", i, opts->npoints);
- ptarray_free(opts);
- if ( i ) continue;
- else break; /* Don't scan holes if shell is collapsed */
- }
-
- /* Add ring to simplified polygon */
- if( lwpoly_add_ring(opoly, opts) == LW_FAILURE )
- {
- lwpoly_free(opoly);
- return NULL;
- }
- }
-
- LWDEBUGF(3, "simplified polygon with %d rings", ipoly->nrings);
- opoly->type = ipoly->type;
-
- if( lwpoly_is_empty(opoly) )
- {
- lwpoly_free(opoly);
- return NULL;
- }
-
- return opoly;
-}
-
/**
* Find the area of the outer ring - sum (area of inner rings).
*/
return ptarray_remove_repeated_points_minpoints(in, tolerance, 2);
}
+/************************************************************************/
+
static void
-ptarray_dp_findsplit(POINTARRAY *pts, int p1, int p2, int *split, double *dist)
+ptarray_dp_findsplit_in_place(const POINTARRAY *pts, int p1, int p2, int *split, double *dist)
{
int k;
const POINT2D *pk, *pa, *pb;
}
}
*dist = d;
-
- } /* length---should be redone if can == 0 */
+ }
else
{
LWDEBUG(3, "segment too short, no split/no dist");
*dist = -1;
}
+}
+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;
}
-POINTARRAY *
-ptarray_simplify(POINTARRAY *inpts, double epsilon, unsigned int minpts)
+void
+ptarray_simplify_in_place(POINTARRAY *pa, double epsilon, unsigned int minpts)
{
- int *stack; /* recursion stack */
- int sp=-1; /* recursion stack pointer */
+ 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;
+ int outn = 0;
+ int pai = 0;
+ int i;
double dist;
- POINTARRAY *outpts;
- POINT4D pt;
-
double eps_sqr = epsilon * epsilon;
- /* Allocate recursion stack */
- stack = lwalloc(sizeof(int)*inpts->npoints);
-
- p1 = 0;
- stack[++sp] = inpts->npoints-1;
-
- LWDEBUGF(2, "Input has %d pts and %d dims", inpts->npoints,
- FLAGS_NDIMS(inpts->flags));
+ /* Do not try to simplify really short things */
+ if (pa->npoints < 3) return;
- /* Allocate output POINTARRAY, and add first point. */
- outpts = ptarray_construct_empty(FLAGS_GET_Z(inpts->flags), FLAGS_GET_M(inpts->flags), inpts->npoints);
- getPoint4d_p(inpts, 0, &pt);
- ptarray_append_point(outpts, &pt, LW_FALSE);
+ /* 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;
+ }
- LWDEBUG(3, "Added P0 to simplified point array (size 1)");
+ p1 = 0;
+ stack[++sp] = pa->npoints-1;
+ /* Add first point to output list */
+ outlist[outn++] = 0;
do
{
+ ptarray_dp_findsplit_in_place(pa, p1, stack[sp], &split, &dist);
- ptarray_dp_findsplit(inpts, p1, stack[sp], &split, &dist);
-
- LWDEBUGF(3, "Farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, dist);
-
- if (dist > eps_sqr || ( outpts->npoints+sp+1 < minpts && dist >= 0 ) )
+ if ((dist > eps_sqr) || ((outn + sp+1 < minpts) && (dist >= 0)))
{
- LWDEBUGF(4, "Added P%d to stack (outpts:%d)", split, sp);
stack[++sp] = split;
}
else
{
- getPoint4d_p(inpts, stack[sp], &pt);
- LWDEBUGF(4, "npoints , minpoints %d %d", outpts->npoints, minpts);
- ptarray_append_point(outpts, &pt, LW_FALSE);
-
- LWDEBUGF(4, "Added P%d to simplified point array (size: %d)", stack[sp], outpts->npoints);
-
+ outlist[outn++] = stack[sp];
p1 = stack[sp--];
}
+ }
+ while (!(sp<0));
- LWDEBUGF(4, "stack pointer = %d", sp);
+ /* 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++)
+ {
+ int j = outlist[i];
+ /* Indexes the same, means no copy required */
+ if (j == pai)
+ {
+ pai++;
+ continue;
+ }
+ /* Indexes different, copy value down */
+ ptarray_copy_point(pa, j, pai++);
}
- while (! (sp<0) );
- lwfree(stack);
- return outpts;
+ /* 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;
}
+/************************************************************************/
+
/**
* Find the 2d length of the given #POINTARRAY, using circular
* arc interpolation between each coordinate triple.
projects / postgis / commitdiff