Commit 3233ccd1 authored by Darafei Praliaskouski's avatar Darafei Praliaskouski

ST_LocateBetween[Elevations]: support POLYGON

Output may be invalid for non-convex, area is preserved.

References #4155
Closes https://github.com/postgis/postgis/pull/319



git-svn-id: http://svn.osgeo.org/postgis/trunk@16952 b70326c6-7e19-0410-871a-916f4a2858ee
parent 6df47db9
Pipeline #34373613 passed with stage
in 34 minutes and 14 seconds
......@@ -27,7 +27,7 @@ PostGIS 3.0.0
Praliaskouski)
- #3457, Fix raster envelope shortcut in ST_Clip (Sai-bot)
- #4215, Use floating point compare in ST_DumpAsPolygons (Darafei Praliaskouski)
- #4155, Support for GEOMETRYCOLLECTION in ST_LocateBetween (Darafei
- #4155, Support for GEOMETRYCOLLECTION and POLYGON in ST_LocateBetween (Darafei
Praliaskouski)
PostGIS 2.5.0
......
......@@ -546,15 +546,15 @@ static void test_lwline_interpolate_points(void)
static void test_lwline_clip(void)
{
LWCOLLECTION *c;
LWLINE *line = NULL;
LWLINE *l51 = NULL;
LWGEOM *line = NULL;
LWGEOM *l51 = NULL;
char *ewkt;
/* Vertical line with vertices at y integers */
l51 = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
l51 = lwgeom_from_wkt("LINESTRING(0 0, 0 1, 0 2, 0 3, 0 4)", LW_PARSER_CHECK_NONE);
/* Clip in the middle, mid-range. */
c = lwline_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.5, 2.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
......@@ -562,7 +562,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Clip off the top. */
c = lwline_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', 3.5, 5.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 3.5,0 4))");
......@@ -570,7 +570,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Clip off the bottom. */
c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 2.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 2.5))" );
......@@ -578,7 +578,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Range holds entire object. */
c = lwline_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.5, 5.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))" );
......@@ -586,7 +586,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Clip on vertices. */
c = lwline_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', 1.0, 2.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1,0 2))" );
......@@ -594,7 +594,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Clip on vertices off the bottom. */
c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 2.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 0,0 1,0 2))" );
......@@ -602,7 +602,7 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* Clip on top. */
c = lwline_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0);
c = lwgeom_clip_to_ordinate_range(l51, 'Y', -1.0, 0.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(0 0))" );
......@@ -610,107 +610,144 @@ static void test_lwline_clip(void)
lwcollection_free(c);
/* ST_LocateBetweenElevations(ST_GeomFromEWKT('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)'), 1, 2)) */
line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
/* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 2)) */
line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 2.0);
line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 2.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("a = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2,1 1 1))" );
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
/* ST_LocateBetweenElevations('LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)', 1, 1)) */
line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
line = lwgeom_from_wkt("LINESTRING(1 2 3, 4 5 6, 6 6 6, 1 1 1)", LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("b = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
/* ST_LocateBetweenElevations('LINESTRING(1 1 1, 1 2 2)', 1,1) */
line = (LWLINE*)lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
c = lwline_clip_to_ordinate_range(line, 'Z', 1.0, 1.0);
line = lwgeom_from_wkt("LINESTRING(1 1 1, 1 2 2)", LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 1.0, 1.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 1 1))" );
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
lwline_free(l51);
lwgeom_free(l51);
}
static void
test_lwpoly_clip(void)
{
LWCOLLECTION *c;
LWGEOM *g = NULL;
char *ewkt;
g = lwgeom_from_wkt(
"POLYGON ((0.51 -0.25, 1.27 -0.14, 1.27 0.25, 0.6 0.3, 0.7 0.7, 1.2 0.7, 0.8 0.5, 1.3 0.4, 1.2 1.2, 0.5 1.2, 0.5 -0.1, 0.3 -0.1, 0.3 1.3, -0.18 1.25, -0.17 -0.25, 0.51 -0.25))",
LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(g, 'X', 0.0, 1.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM *)c);
// printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(
ewkt,
"MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))");
lwfree(ewkt);
lwcollection_free(c);
lwgeom_free(g);
g = lwgeom_from_wkt(
"MULTIPOLYGON(((0.51 -0.25,1 -0.179078947368,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1.2,0.5 1.2,0.5 -0.1,0.3 -0.1,0.3 1.3,0 1.26875,0 -0.25,0.51 -0.25)))",
LW_PARSER_CHECK_NONE);
c = lwgeom_clip_to_ordinate_range(g, 'Y', 0.0, 1.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM *)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(
ewkt,
"MULTIPOLYGON(((1 0,1 0.270149253731,0.6 0.3,0.7 0.7,1 0.7,1 0.6,0.8 0.5,1 0.46,1 1,0.5 1,0.5 0,0.3 0,0.3 1,0 1,0 0,1 0)))");
lwfree(ewkt);
lwcollection_free(c);
lwgeom_free(g);
}
static void test_lwmline_clip(void)
{
LWCOLLECTION *c;
char *ewkt;
LWMLINE *mline = NULL;
LWLINE *line = NULL;
LWGEOM *mline = NULL;
LWGEOM *line = NULL;
/*
** Set up the input line. Trivial one-member case.
*/
mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
mline = lwgeom_from_wkt("MULTILINESTRING((0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
/* Clip in the middle, mid-range. */
c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 1.5, 2.5);
c = lwgeom_clip_to_ordinate_range(mline, 'Y', 1.5, 2.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0 1.5,0 2,0 2.5))");
lwfree(ewkt);
lwcollection_free(c);
lwmline_free(mline);
lwgeom_free(mline);
/*
** Set up the input line. Two-member case.
*/
mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 1,1 2,1 3,1 4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
mline = lwgeom_from_wkt("MULTILINESTRING((1 0,1 1,1 2,1 3,1 4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
/* Clip off the top. */
c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 3.5, 5.5);
c = lwgeom_clip_to_ordinate_range(mline, 'Y', 3.5, 5.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((1 3.5,1 4),(0 3.5,0 4))");
lwfree(ewkt);
lwcollection_free(c);
lwmline_free(mline);
lwgeom_free(mline);
/*
** Set up staggered input line to create multi-type output.
*/
mline = (LWMLINE*)lwgeom_from_wkt("MULTILINESTRING((1 0,1 -1,1 -2,1 -3,1 -4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
mline =
lwgeom_from_wkt("MULTILINESTRING((1 0,1 -1,1 -2,1 -3,1 -4), (0 0,0 1,0 2,0 3,0 4))", LW_PARSER_CHECK_NONE);
/* Clip from 0 upwards.. */
c = lwcollection_clip_to_ordinate_range((LWCOLLECTION *)mline, 'Y', 0.0, 2.5);
c = lwgeom_clip_to_ordinate_range(mline, 'Y', 0.0, 2.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "GEOMETRYCOLLECTION(POINT(1 0),LINESTRING(0 0,0 1,0 2,0 2.5))");
lwfree(ewkt);
lwcollection_free(c);
lwmline_free(mline);
lwgeom_free(mline);
/*
** Set up input line from MAC
*/
line = (LWLINE*)lwgeom_from_wkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)", LW_PARSER_CHECK_NONE);
line = lwgeom_from_wkt("LINESTRING(0 0 0 0,1 1 1 1,2 2 2 2,3 3 3 3,4 4 4 4,3 3 3 5,2 2 2 6,1 1 1 7,0 0 0 8)",
LW_PARSER_CHECK_NONE);
/* Clip from 3 to 3.5 */
c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 3.5);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 3.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5))");
......@@ -718,7 +755,7 @@ static void test_lwmline_clip(void)
lwcollection_free(c);
/* Clip from 2 to 3.5 */
c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.5);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3,3.5 3.5 3.5 3.5),(3.5 3.5 3.5 4.5,3 3 3 5,2 2 2 6))");
......@@ -726,7 +763,7 @@ static void test_lwmline_clip(void)
lwcollection_free(c);
/* Clip from 3 to 4 */
c = lwline_clip_to_ordinate_range(line, 'Z', 3.0, 4.0);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 3.0, 4.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((3 3 3 3,4 4 4 4,3 3 3 5))");
......@@ -734,24 +771,20 @@ static void test_lwmline_clip(void)
lwcollection_free(c);
/* Clip from 2 to 3 */
c = lwline_clip_to_ordinate_range(line, 'Z', 2.0, 3.0);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 2.0, 3.0, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((2 2 2 2,3 3 3 3),(3 3 3 5,2 2 2 6))");
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
}
static void test_lwline_clip_big(void)
{
POINTARRAY *pa = ptarray_construct(1, 0, 3);
LWLINE *line = lwline_construct(SRID_UNKNOWN, NULL, pa);
LWGEOM *line = (LWGEOM *)lwline_construct(SRID_UNKNOWN, NULL, pa);
LWCOLLECTION *c;
char *ewkt;
POINT4D p;
......@@ -771,14 +804,14 @@ static void test_lwline_clip_big(void)
p.z = 2.0;
ptarray_set_point4d(pa, 2, &p);
c = lwline_clip_to_ordinate_range(line, 'Z', 0.5, 1.5);
c = lwgeom_clip_to_ordinate_range(line, 'Z', 0.5, 1.5, 0);
ewkt = lwgeom_to_ewkt((LWGEOM*)c);
//printf("c = %s\n", ewkt);
CU_ASSERT_STRING_EQUAL(ewkt, "MULTILINESTRING((0.5 0.5 0.5,1 1 1,1.5 1.5 1.5))" );
lwfree(ewkt);
lwcollection_free(c);
lwline_free(line);
lwgeom_free(line);
}
static void test_geohash_precision(void)
......@@ -1497,6 +1530,7 @@ void algorithms_suite_setup(void)
PG_ADD_TEST(suite,test_point_interpolate);
PG_ADD_TEST(suite,test_lwline_interpolate_points);
PG_ADD_TEST(suite,test_lwline_clip);
PG_ADD_TEST(suite, test_lwpoly_clip);
PG_ADD_TEST(suite,test_lwline_clip_big);
PG_ADD_TEST(suite,test_lwmline_clip);
PG_ADD_TEST(suite,test_geohash_point);
......
......@@ -234,26 +234,6 @@ void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value);
int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value);
/**
* Clip a line based on the from/to range of one of its ordinates. Use for m- and z- clipping
*/
LWCOLLECTION *lwline_clip_to_ordinate_range(const LWLINE *line, char ordinate, double from, double to);
/**
* Clip collection based on the from/to range of one of its ordinates. Use for m- and z- clipping
*/
LWCOLLECTION *lwcollection_clip_to_ordinate_range(const LWCOLLECTION *col, char ordinate, double from, double to);
/**
* Clip a multi-point based on the from/to range of one of its ordinates. Use for m- and z- clipping
*/
LWCOLLECTION *lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to);
/**
* Clip a point based on the from/to range of one of its ordinates. Use for m- and z- clipping
*/
LWCOLLECTION *lwpoint_clip_to_ordinate_range(const LWPOINT *mpoint, char ordinate, double from, double to);
/*
* Geohash
*/
......
......@@ -248,54 +248,49 @@ lwgeom_locate_along(const LWGEOM *lwin, double m, double offset)
* @param ordinate number (1=x, 2=y, 3=z, 4=m)
* @return d value at that ordinate
*/
double lwpoint_get_ordinate(const POINT4D *p, char ordinate)
inline double
lwpoint_get_ordinate(const POINT4D *p, char ordinate)
{
if ( ! p )
if (!p)
{
lwerror("Null input geometry.");
return 0.0;
}
if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
switch (ordinate)
{
lwerror("Cannot extract %c ordinate.", ordinate);
return 0.0;
}
if ( ordinate == 'X' )
case 'X':
return p->x;
if ( ordinate == 'Y' )
case 'Y':
return p->y;
if ( ordinate == 'Z' )
case 'Z':
return p->z;
if ( ordinate == 'M' )
case 'M':
return p->m;
/* X */
return p->x;
}
lwerror("Cannot extract %c ordinate.", ordinate);
return 0.0;
}
/**
* Given a point, ordinate number and value, set that ordinate on the
* point.
*/
void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
inline void
lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
{
if ( ! p )
if (!p)
{
lwerror("Null input geometry.");
return;
}
if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
if (!(ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M'))
{
lwerror("Cannot set %c ordinate.", ordinate);
return;
}
LWDEBUGF(4, " setting ordinate %c to %g", ordinate, value);
switch ( ordinate )
{
case 'X':
......@@ -318,7 +313,14 @@ void lwpoint_set_ordinate(POINT4D *p, char ordinate, double value)
* generate a new point that is proportionally between the input points,
* using the values in the provided dimension as the scaling factors.
*/
int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz, int hasm, char ordinate, double interpolation_value)
inline int
point_interpolate(const POINT4D *p1,
const POINT4D *p2,
POINT4D *p,
int hasz,
int hasm,
char ordinate,
double interpolation_value)
{
static char* dims = "XYZM";
double p1_value = lwpoint_get_ordinate(p1, ordinate);
......@@ -326,41 +328,51 @@ int point_interpolate(const POINT4D *p1, const POINT4D *p2, POINT4D *p, int hasz
double proportion;
int i = 0;
if ( ! ( ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M' ) )
#if PARANOIA_LEVEL > 0
if (!(ordinate == 'X' || ordinate == 'Y' || ordinate == 'Z' || ordinate == 'M'))
{
lwerror("Cannot set %c ordinate.", ordinate);
return 0;
lwerror("Cannot interpolate over %c ordinate.", ordinate);
return LW_FAILURE;
}
if ( FP_MIN(p1_value, p2_value) > interpolation_value ||
FP_MAX(p1_value, p2_value) < interpolation_value )
if (FP_MIN(p1_value, p2_value) > interpolation_value || FP_MAX(p1_value, p2_value) < interpolation_value)
{
lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).", interpolation_value, p1_value, p2_value);
return 0;
lwerror("Cannot interpolate to a value (%g) not between the input points (%g, %g).",
interpolation_value,
p1_value,
p2_value);
return LW_FAILURE;
}
#endif
proportion = fabs((interpolation_value - p1_value) / (p2_value - p1_value));
for ( i = 0; i < 4; i++ )
for (i = 0; i < 4; i++)
{
double newordinate = 0.0;
if ( dims[i] == 'Z' && ! hasz ) continue;
if ( dims[i] == 'M' && ! hasm ) continue;
p1_value = lwpoint_get_ordinate(p1, dims[i]);
p2_value = lwpoint_get_ordinate(p2, dims[i]);
newordinate = p1_value + proportion * (p2_value - p1_value);
lwpoint_set_ordinate(p, dims[i], newordinate);
LWDEBUGF(4, " clip ordinate(%c) p1_value(%g) p2_value(%g) proportion(%g) newordinate(%g) ", dims[i], p1_value, p2_value, proportion, newordinate );
if (dims[i] == 'Z' && !hasz)
continue;
if (dims[i] == 'M' && !hasm)
continue;
if (dims[i] == ordinate)
lwpoint_set_ordinate(p, dims[i], interpolation_value);
else
{
p1_value = lwpoint_get_ordinate(p1, dims[i]);
p2_value = lwpoint_get_ordinate(p2, dims[i]);
newordinate = p1_value + proportion * (p2_value - p1_value);
lwpoint_set_ordinate(p, dims[i], newordinate);
}
}
return 1;
return LW_SUCCESS;
}
/**
* Clip an input POINT between two values, on any ordinate input.
*/
LWCOLLECTION*
static inline LWCOLLECTION *
lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from, double to)
{
LWCOLLECTION *lwgeom_out = NULL;
......@@ -368,18 +380,6 @@ lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from,
POINT4D p4d;
double ordinate_value;
/* Nothing to do with NULL */
if ( ! point )
lwerror("Null input geometry.");
/* Ensure 'from' is less than 'to'. */
if ( to < from )
{
double t = from;
from = to;
to = t;
}
/* Read Z/M info */
hasz = lwgeom_has_z(lwpoint_as_lwgeom(point));
hasm = lwgeom_has_m(lwpoint_as_lwgeom(point));
......@@ -396,39 +396,19 @@ lwpoint_clip_to_ordinate_range(const LWPOINT *point, char ordinate, double from,
lwcollection_add_lwgeom(lwgeom_out, lwpoint_as_lwgeom(lwp));
}
/* Set the bbox, if necessary */
if ( lwgeom_out->bbox )
{
lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
}
return lwgeom_out;
}
/**
* Clip an input MULTIPOINT between two values, on any ordinate input.
*/
LWCOLLECTION*
static inline LWCOLLECTION *
lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double from, double to)
{
LWCOLLECTION *lwgeom_out = NULL;
char hasz, hasm;
uint32_t i;
/* Nothing to do with NULL */
if ( ! mpoint )
lwerror("Null input geometry.");
/* Ensure 'from' is less than 'to'. */
if ( to < from )
{
double t = from;
from = to;
to = t;
}
/* Read Z/M info */
hasz = lwgeom_has_z(lwmpoint_as_lwgeom(mpoint));
hasm = lwgeom_has_m(lwmpoint_as_lwgeom(mpoint));
......@@ -453,62 +433,112 @@ lwmpoint_clip_to_ordinate_range(const LWMPOINT *mpoint, char ordinate, double fr
}
/* Set the bbox, if necessary */
if ( lwgeom_out->bbox )
{
if (mpoint->bbox)
lwgeom_refresh_bbox((LWGEOM*)lwgeom_out);
}
return lwgeom_out;
}
/**
* Clip an input COLLECTION between two values, on any ordinate input.
*/
LWCOLLECTION *
lwcollection_clip_to_ordinate_range(const LWCOLLECTION *icol, char ordinate, double from, double to)
static inline POINTARRAY *
ptarray_clamp_to_ordinate_range(const POINTARRAY *ipa, char ordinate, double from, double to)
{
LWCOLLECTION *lwgeom_out = NULL;
POINT4D p1, p2;
POINTARRAY *opa;
double ovp1, ovp2;
POINT4D *t;
int8_t p1out, p2out; /* -1 - smaller than from, 0 - in range, 1 - larger than to */
uint32_t i;
uint8_t hasz = FLAGS_GET_Z(ipa->flags);
uint8_t hasm = FLAGS_GET_M(ipa->flags);
if (!icol)
{
lwerror("Null input geometry.");
return NULL;
}
assert(from <= to);
if (icol->ngeoms == 1)
lwgeom_out = lwgeom_clip_to_ordinate_range(icol->geoms[0], ordinate, from, to, 0);
else
t = lwalloc(sizeof(POINT4D));
/* Initial storage */
opa = ptarray_construct_empty(hasz, hasm, ipa->npoints);
/* Add first point */
getPoint4d_p(ipa, 0, &p1);
ovp1 = lwpoint_get_ordinate(&p1, ordinate);
p1out = (ovp1 < from) ? -1 : ((ovp1 > to) ? 1 : 0);
if (from <= ovp1 && ovp1 <= to)
ptarray_append_point(opa, &p1, LW_FALSE);
/* Loop on all other input points */
for (i = 1; i < ipa->npoints; i++)
{
LWCOLLECTION *col;
char hasz = lwgeom_has_z(lwcollection_as_lwgeom(icol));
char hasm = lwgeom_has_m(lwcollection_as_lwgeom(icol));
uint32_t i;
lwgeom_out = lwcollection_construct_empty(icol->type, icol->srid, hasz, hasm);
FLAGS_SET_Z(lwgeom_out->flags, hasz);
FLAGS_SET_M(lwgeom_out->flags, hasm);
for (i = 0; i < icol->ngeoms; i++)
getPoint4d_p(ipa, i, &p2);
ovp2 = lwpoint_get_ordinate(&p2, ordinate);
p2out = (ovp2 < from) ? -1 : ((ovp2 > to) ? 1 : 0);
if (p1out == 0 && p2out == 0) /* both visible */
{
col = lwgeom_clip_to_ordinate_range(icol->geoms[i], ordinate, from, to, 0);
if (col)
{
if (col->type != icol->type)
lwgeom_out->type = COLLECTIONTYPE;
lwgeom_out = lwcollection_concat_in_place(lwgeom_out, col);
lwcollection_release(col);
}
ptarray_append_point(opa, &p2, LW_FALSE);
}
if (lwgeom_out->bbox)
lwgeom_refresh_bbox((LWGEOM *)lwgeom_out);
else if (p1out == p2out && p1out != 0) /* both invisible */
{
/* skip */
}
else if (p1out == -1 && p2out == 0)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
ptarray_append_point(opa, t, LW_FALSE);
ptarray_append_point(opa, &p2, LW_FALSE);
}
else if (p1out == -1 && p2out == 1)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
ptarray_append_point(opa, t, LW_FALSE);
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
ptarray_append_point(opa, t, LW_FALSE);
}
else if (p1out == 0 && p2out == -1)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
ptarray_append_point(opa, t, LW_FALSE);
}
else if (p1out == 0 && p2out == 1)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
ptarray_append_point(opa, t, LW_FALSE);
}
else if (p1out == 1 && p2out == -1)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
ptarray_append_point(opa, t, LW_FALSE);
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, from);
ptarray_append_point(opa, t, LW_FALSE);
}
else if (p1out == 1 && p2out == 0)
{
point_interpolate(&p1, &p2, t, hasz, hasm, ordinate, to);
ptarray_append_point(opa, t, LW_FALSE);
ptarray_append_point(opa, &p2, LW_FALSE);
}
p1 = p2;
p1out = p2out;
LW_ON_INTERRUPT(ptarray_free(opa); return NULL);
}
return lwgeom_out;
}
if (ptarray_is_closed(ipa) && opa->npoints > 0)
{
getPoint4d_p(opa, 0, &p1);
ptarray_append_point(opa, &p1, LW_FALSE);
}
lwfree(t);