Skip to content

Commit

Permalink
Add ST_LineClipZ(geometry, from, to) SQL and C functions.
Browse files Browse the repository at this point in the history 
git-svn-id: http://svn.osgeo.org/postgis/trunk@3456 b70326c6-7e19-0410-871a-916f4a2858ee
  • Loading branch information
@pramsey
pramsey committed Dec 21, 2008
1 parent 2d6e86f commit f95ff43
Show file tree
Hide file tree
Showing 5 changed files with 114 additions and 17 deletions.
36 changes: 35 additions & 1 deletion liblwgeom/cunit/cu_algorithm.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,8 @@ CU_pSuite register_cg_suite(void)
(NULL == CU_add_test(pSuite, "test_lwpoint_set_ordinate()", test_lwpoint_set_ordinate)) ||
(NULL == CU_add_test(pSuite, "test_lwpoint_get_ordinate()", test_lwpoint_get_ordinate)) ||
(NULL == CU_add_test(pSuite, "test_lwpoint_interpolate()", test_lwpoint_interpolate)) ||
(NULL == CU_add_test(pSuite, "test_lwline_clip()", test_lwline_clip))
(NULL == CU_add_test(pSuite, "test_lwline_clip()", test_lwline_clip)) ||
(NULL == CU_add_test(pSuite, "test_lwline_clip_big()", test_lwline_clip_big))
)
{
CU_cleanup_registry();
Expand Down Expand Up @@ -674,3 +675,36 @@ void test_lwline_clip(void)

}

void test_lwline_clip_big(void)
{
POINTARRAY *pa = ptarray_construct(1, 0, 3);
LWLINE *line = lwline_construct(-1, NULL, pa);
LWCOLLECTION *c;
char *ewkt;
int rv = 0;

p->x = 0.0;
p->y = 0.0;
p->z = 0.0;
setPoint4d(pa, 0, p);

p->x = 1.0;
p->y = 1.0;
p->z = 1.0;
setPoint4d(pa, 1, p);

p->x = 2.0;
p->y = 2.0;
p->z = 2.0;
setPoint4d(pa, 2, p);

c = lwline_clip_to_ordinate_range(line, 2, 0.5, 1.5);
ewkt = lwgeom_to_ewkt((LWGEOM*)c,0);
//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);
lwfree_collection(c);
lwfree_line(line);
}

2 changes: 2 additions & 0 deletions liblwgeom/cunit/cu_algorithm.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,3 +35,5 @@ void test_lwpoint_set_ordinate(void);
void test_lwpoint_get_ordinate(void);
void test_lwpoint_interpolate(void);
void test_lwline_clip(void);
void test_lwline_clip_big(void);

47 changes: 31 additions & 16 deletions liblwgeom/lwalgorithm.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -403,9 +403,9 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
double ordinate_value_p = 0.0, ordinate_value_q = 0.0;
char hasz = TYPE_HASZ(line->type);
char hasm = TYPE_HASM(line->type);
char dims = TYPE_NDIMS(line->type);
char hassrid = TYPE_HASSRID(line->type);


LWDEBUGF(5, "hassrid = %d", hassrid);

/* Null input, nothing we can do. */
Expand All @@ -423,12 +423,10 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
to = t;
}

int ndims = TYPE_NDIMS(line->type);

/* Asking for an ordinate we don't have. Error. */
if( ordinate >= ndims )
if( ordinate >= dims )
{
lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, ndims);
lwerror("Cannot clip on ordinate %d in a %d-d geometry.", ordinate, dims);
return NULL;
}

Expand Down Expand Up @@ -471,7 +469,7 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
if( ! added_last_point )
{
if( dp ) lwfree(dp);
dp = dynptarray_create(64, ndims);
dp = dynptarray_create(64, line->type);
}
rv = dynptarray_addPoint4d(dp, p, 0);
added_last_point = 2; /* Added on a boundary. */
Expand All @@ -488,21 +486,21 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
/* We didn't add the previous point, so this is a new segment.
* Make a new point array. */
if( dp ) lwfree(dp);
dp = dynptarray_create(64, ndims);
dp = dynptarray_create(64, line->type);

/* We're transiting into the range so add an interpolated
* point at the range boundary. */
if ( ordinate_value_p > from && ordinate_value_p < to && i > 0 )
{
double interpolation_value;
(ordinate_value_q > to) ? (interpolation_value = to) : (interpolation_value = from);
rv = lwpoint_interpolate(q, p, r, ndims, ordinate, interpolation_value);
rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
rv = dynptarray_addPoint4d(dp, r, 0);
}
}
/* Add the current vertex to the point array. */
rv = dynptarray_addPoint4d(dp, p, 0);
added_last_point = 1;
added_last_point = 1; /* Added inside range. */
}
/* Is this point inside the ordinate range? No. */
else
Expand All @@ -513,7 +511,7 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
* to the point array at the range boundary. */
double interpolation_value;
(ordinate_value_p > to) ? (interpolation_value = to) : (interpolation_value = from);
rv = lwpoint_interpolate(q, p, r, ndims, ordinate, interpolation_value);
rv = lwpoint_interpolate(q, p, r, dims, ordinate, interpolation_value);
rv = dynptarray_addPoint4d(dp, r, 0);
}
else if ( added_last_point == 2 )
Expand All @@ -526,21 +524,21 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
* so we need to add *two* interpolated points! */
pa_out = ptarray_construct(hasz, hasm, 2);
/* Interpolate lower point. */
rv = lwpoint_interpolate(p, q, r, ndims, ordinate, from);
rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
setPoint4d(pa_out, 0, r);
/* Interpolate upper point. */
rv = lwpoint_interpolate(p, q, r, ndims, ordinate, to);
rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
setPoint4d(pa_out, 1, r);
}
else if ( ordinate_value_q > to && ordinate_value_p < from ) {
/* We just hopped over the whole range, from top to bottom,
* so we need to add *two* interpolated points! */
pa_out = ptarray_construct(hasz, hasm, 2);
/* Interpolate upper point. */
rv = lwpoint_interpolate(p, q, r, ndims, ordinate, to);
rv = lwpoint_interpolate(p, q, r, dims, ordinate, to);
setPoint4d(pa_out, 0, r);
/* Interpolate lower point. */
rv = lwpoint_interpolate(p, q, r, ndims, ordinate, from);
rv = lwpoint_interpolate(p, q, r, dims, ordinate, from);
setPoint4d(pa_out, 1, r);
}
/* We have an extant point-array, save it out to a multi-line. */
Expand All @@ -553,19 +551,28 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
if( dp->pa->npoints == 1 )
{
oline = (LWGEOM*)lwpoint_construct(line->SRID, NULL, dp->pa);
oline->type = lwgeom_makeType(hasz, hasm, hassrid, POINTTYPE);
lwgeom_out->type = lwgeom_makeType(hasz, hasm, hassrid, COLLECTIONTYPE);
}
else
{
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
}
}
else
{
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, pa_out);
}
lwgeom_out->ngeoms++;
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
if( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
{
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
}
else
{
lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
}
lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
lwgeom_dropBBOX((LWGEOM*)lwgeom_out);
lwgeom_addBBOX((LWGEOM*)lwgeom_out);
Expand All @@ -582,8 +589,16 @@ LWCOLLECTION *lwline_clip_to_ordinate_range(LWLINE *line, int ordinate, double f
if( dp && dp->pa->npoints > 0 ) {
LWGEOM *oline;
oline = (LWGEOM*)lwline_construct(line->SRID, NULL, dp->pa);
oline->type = lwgeom_makeType(hasz, hasm, hassrid, LINETYPE);
lwgeom_out->ngeoms++;
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
if( lwgeom_out->geoms ) /* We can't just realloc, since repalloc chokes on a starting null ptr. */
{
lwgeom_out->geoms = lwrealloc(lwgeom_out->geoms, sizeof(LWGEOM*) * lwgeom_out->ngeoms);
}
else
{
lwgeom_out->geoms = lwalloc(sizeof(LWGEOM*) * lwgeom_out->ngeoms);
}
lwgeom_out->geoms[lwgeom_out->ngeoms - 1] = oline;
lwgeom_dropBBOX((LWGEOM*)lwgeom_out);
lwgeom_addBBOX((LWGEOM*)lwgeom_out);
Expand Down
40 changes: 40 additions & 0 deletions lwgeom/lwgeom_functions_analytic.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@ LWCOLLECTION *simplify2d_collection(const LWCOLLECTION *igeom, double dist);
LWGEOM *simplify2d_lwgeom(const LWGEOM *igeom, double dist);
Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS);
Datum crossingDirection(PG_FUNCTION_ARGS);
Datum ST_LineClipZ(PG_FUNCTION_ARGS);

double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point);
Expand Down Expand Up @@ -993,6 +994,45 @@ Datum crossingDirection(PG_FUNCTION_ARGS)

}

PG_FUNCTION_INFO_V1(ST_LineClipZ);
Datum ST_LineClipZ(PG_FUNCTION_ARGS)
{
PG_LWGEOM *geom_in = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
double from = PG_GETARG_FLOAT8(1);
double to = PG_GETARG_FLOAT8(2);
LWCOLLECTION *geom_out = NULL;
LWLINE *line_in = NULL;
uchar type = (uchar)SERIALIZED_FORM(geom_in)[0];
char geomtype = TYPE_GETTYPE(type);
char hasz = TYPE_HASZ(type);
static int ordinate = 2; /* Z */

if ( geomtype != LINETYPE )
{
elog(ERROR,"This function only accepts LINESTRING as arguments.");
PG_RETURN_NULL();
}

if( ! hasz )
{
elog(ERROR,"This function only accepts LINESTRING with Z values as arguments.");
PG_RETURN_NULL();
}

line_in = lwline_deserialize(SERIALIZED_FORM(geom_in));
geom_out = lwline_clip_to_ordinate_range(line_in, ordinate, from, to);
lwfree_line(line_in);

if( ! geom_out ) {
elog(ERROR,"The lwline_clip_to_ordinate_range returned null.");
PG_RETURN_NULL();
}

PG_FREE_IF_COPY(geom_in, 0);
PG_RETURN_POINTER(pglwgeom_serialize((LWGEOM*)geom_out));

}

/***********************************************************************
* [email protected]
***********************************************************************/
Expand Down
6 changes: 6 additions & 0 deletions lwgeom/lwpostgis.sql.in.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3972,6 +3972,12 @@ CREATEFUNCTION ST_CrossingDirection(geometry, geometry)
AS 'MODULE_PATHNAME', 'crossingDirection'
LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable);

-- Only accepts LINESTRING as parameters.
-- Availability: 1.4.0
CREATEFUNCTION ST_LineClipZ(geometry, float8, float8)
RETURNS geometry
AS 'MODULE_PATHNAME', 'ST_LineClipZ'
LANGUAGE 'C' _IMMUTABLE_STRICT; -- WITH (isstrict,iscachable);

#if POSTGIS_GEOS_VERSION >= 30
-- Requires GEOS >= 3.0.0
Expand Down

0 comments on commit f95ff43

Please sign in to comment.