[postgis-users] postgis 2.0.x segfault crashes
Mario Jurcevic
mario at geosar.ch
Fri Jul 5 08:21:19 PDT 2013
The new functions are not working well with geometries that mix linear
and circular parts.
With the new semantics they are not usable on MULTICURVE and MULTISURFACE.
For example the st_linetocurve will not try to build circular parts by
looking at vertexes position
( angles ) but only outputs circular parts that could be built with
vertex triples.
The semantics in 1.5.3 version:
st_curvetoline: linearize the curved part of geometries and maintain the
linear part
st_linetocurve: discover the part that could be really converted to
curved ( preserving the original shape ) and maintain the linear part.
Take a look at the source of postgis until 2.0 for this functions.
For the st_curvetoline this part in attachment is missing.
Mario
On 07/05/2013 12:52 PM, Sandro Santilli wrote:
> On Fri, Jul 05, 2013 at 11:24:39AM +0200, Mario Jurcevic wrote:
>> The is problem, is with the table "Bodenbedeckung__BBNachfuehrung"
>> which has big polygon inside.
>>
>> Without using custom functions the system does not crashes.
>>
>> The custom functions wraps the postgis functions.
>>
>> Exception is the st_curvetoline and st_linetocurve functions that
>> during migration of postgis 1.5 to 2.0 were changed.
>>
>> Personal opinion: the new behavior is useless.
> Uhm, documentation doesn't say anything about the semantic change:
> http://postgis.net/docs/ST_CurveToLine.html
> http://postgis.net/docs/ST_LineToCurve.html
>
> How did they change ?
>
>> We have modified the 2.0 source to include the old st_cur* st_line*
>> behavior, that seems on 64 bit inside 2.0 source are not working.
> So this sounds like a bug in your patches, right ?
>
> I'll be useful to know more about the semantic change that made you
> patch those functions.
>
> --strk;
>
>> Thanks for taking time for looking,
>> Mario Jurcevic
>>
>>
>> On 07/05/2013 10:08 AM, Sandro Santilli wrote:
>>> On Thu, Jul 04, 2013 at 09:27:11PM +0200, Mario Jurcevic wrote:
>>>> Sorry this one is the right:
>>>>
>>>>
>>>> http://ubuntuone.com/202BqXg3Lvtt9zGuu3sDcx
>>> The SQL file contains calls to custom functions that
>>> are not available here. I handled to get the import
>>> started by creating an "import_cadastre" schema, but
>>> after 2.5k lines of sql the process ends with:
>>>
>>> psql:interlis_import_51d5cb94e16e41.40947245.sql:2661: ERROR: function geo_curvedtolinear(geometry) does not exist
>>> LINE 1: ...o_LinearToCurved(ST_Multi(ST_BuildArea(ST_Collect(Geo_Curved...
>>>
>>> Can you get the segfault crash without using custom
>>> functions ?
>>>
>>> --strk;
> _______________________________________________
> postgis-users mailing list
> postgis-users at lists.osgeo.org
> http://lists.osgeo.org/cgi-bin/mailman/listinfo/postgis-users
-------------- next part --------------
# implement old st_linetocurve behavior
# add these lines to postgis/lwgeom_sqlmm.c
/*******************************************************************************
* ST_CurveToLine 1.5.3 BEHAVIOR
******************************************************************************/
PG_FUNCTION_INFO_V1(LWGEOM_curve_segmentize_old);
Datum LWGEOM_curve_segmentize_old(PG_FUNCTION_ARGS)
{
GSERIALIZED *geom = (GSERIALIZED *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
uint32 perQuad = PG_GETARG_INT32(1);
GSERIALIZED *ret;
LWGEOM *igeom = NULL, *ogeom = NULL;
POSTGIS_DEBUG(2, "LWGEOM_curve_segmentize_old called.");
if (perQuad < 0)
{
elog(ERROR, "2nd argument must be positive.");
PG_RETURN_NULL();
}
POSTGIS_DEBUGF(3, "perQuad = %d", perQuad);
igeom = lwgeom_from_gserialized(geom);
/*
removed to always return a line if curve doesn't have an arc
if ( ! has_arc(igeom) )
{
PG_RETURN_POINTER(geom);
}
*/
ogeom = lwgeom_segmentize_old(igeom, 32);
lwgeom_free(igeom);
if (ogeom == NULL)
PG_RETURN_NULL();
ret = geometry_serialize(ogeom);
lwgeom_free(ogeom);
PG_FREE_IF_COPY(geom, 0);
PG_RETURN_POINTER(ret);
}
/*******************************************************************************
* END ST_CurveToLine 1.5.3 BEHAVIOR
******************************************************************************/
# add these lines to liblwgeom/lwsegmentize.c
/*******************************************************************************
* ST_CurveToLine 1.5.3 BEHAVIOR
******************************************************************************/
uint32_t has_arc_old(LWGEOM *geom);
double lwcircle_center_old(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result);
double interpolate_arc_old(double angle, double zm1, double a1, double zm2, double a2);
POINTARRAY *lwcircle_segmentize_old(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad);
LWLINE *lwcurve_segmentize_old(LWCIRCSTRING *icurve, uint32_t perQuad);
LWLINE *lwcompound_segmentize_old(LWCOMPOUND *icompound, uint32_t perQuad);
LWPOLY *lwcurvepoly_segmentize_old(LWCURVEPOLY *curvepoly, uint32_t perQuad);
LWMLINE *lwmcurve_segmentize_old(LWMCURVE *mcurve, uint32_t perQuad);
LWMPOLY *lwmsurface_segmentize_old(LWMSURFACE *msurface, uint32_t perQuad);
LWCOLLECTION *lwcollection_segmentize_old(LWCOLLECTION *collection, uint32_t perQuad);
LWGEOM *lwgeom_segmentize_old(LWGEOM *geom, uint32_t perQuad);
/*
* Determines (recursively in the case of collections) whether the geometry
* contains at least on arc geometry or segment.
*/
uint32_t
has_arc_old(LWGEOM *geom)
{
LWCOLLECTION *col;
int i;
LWDEBUG(2, "has_arc called.");
switch (geom->type)
{
case POINTTYPE:
case LINETYPE:
case POLYGONTYPE:
case MULTIPOINTTYPE:
case MULTILINETYPE:
case MULTIPOLYGONTYPE:
return 0;
case CIRCSTRINGTYPE:
return 1;
/* It's a collection that MAY contain an arc */
default:
col = (LWCOLLECTION *)geom;
for (i=0; i<col->ngeoms; i++)
{
if (has_arc_old(col->geoms[i]) == 1) return 1;
}
return 0;
}
}
/*
* Determines the center of the circle defined by the three given points.
* In the event the circle is complete, the midpoint of the segment defined
* by the first and second points is returned. If the points are colinear,
* as determined by equal slopes, then NULL is returned. If the interior
* point is coincident with either end point, they are taken as colinear.
*/
double
lwcircle_center_old(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result)
{
POINT4D *c;
double cx, cy, cr;
double temp, bc, cd, det;
LWDEBUGF(2, "lwcircle_center called (%.16f, %.16f), (%.16f, %.16f), (%.16f, %.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);
/* Closed circle */
if (fabs(p1->x - p3->x) < EPSILON_SQLMM
&& fabs(p1->y - p3->y) < EPSILON_SQLMM)
{
cx = p1->x + (p2->x - p1->x) / 2.0;
cy = p1->y + (p2->y - p1->y) / 2.0;
c = lwalloc(sizeof(POINT2D));
c->x = cx;
c->y = cy;
*result = c;
cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
return cr;
}
temp = p2->x*p2->x + p2->y*p2->y;
bc = (p1->x*p1->x + p1->y*p1->y - temp) / 2.0;
cd = (temp - p3->x*p3->x - p3->y*p3->y) / 2.0;
det = (p1->x - p2->x)*(p2->y - p3->y)-(p2->x - p3->x)*(p1->y - p2->y);
/* Check colinearity */
if (fabs(det) < EPSILON_SQLMM)
{
*result = NULL;
return -1.0;
}
det = 1.0 / det;
cx = (bc*(p2->y - p3->y)-cd*(p1->y - p2->y))*det;
cy = ((p1->x - p2->x)*cd-(p2->x - p3->x)*bc)*det;
c = lwalloc(sizeof(POINT4D));
c->x = cx;
c->y = cy;
*result = c;
cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));
return cr;
}
double
interpolate_arc_old(double angle, double zm1, double a1, double zm2, double a2)
{
double frac = fabs((angle - a1) / (a2 - a1));
double result = frac * (zm2 - zm1) + zm1;
LWDEBUG(2, "interpolate_arc called.");
LWDEBUGF(3, "interpolate_arc: angle=%.16f, a1=%.16f, a2=%.16f, z1=%.16f, z2=%.16f, frac=%.16f, result=%.16f", angle, a1, a2, zm1, zm2, frac, result);
return result;
}
POINTARRAY *
lwcircle_segmentize_old(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32_t perQuad)
{
POINTARRAY *result;
POINT4D pbuf;
size_t ptsize = sizeof(POINT4D);
unsigned int ptcount;
POINT4D *pt;
POINT4D *center;
double radius = 0.0,
sweep = 0.0,
angle = 0.0,
increment = 0.0;
double a1, a2, a3, i;
LWDEBUG(2, "lwcircle_segmentize_old called. ");
radius = lwcircle_center_old(p1, p2, p3, ¢er);
if (radius < 0)
{
/* No center because points are colinear */
LWDEBUGF(3, "lwcircle_segmentize_old, (NULL) radius=%.16f", radius);
return NULL;
}
LWDEBUGF(3, "lwcircle_segmentize_old, (%.16f, %.16f) radius=%.16f", center->x, center->y, radius);
a1 = atan2(p1->y - center->y, p1->x - center->x);
a2 = atan2(p2->y - center->y, p2->x - center->x);
a3 = atan2(p3->y - center->y, p3->x - center->x);
LWDEBUGF(3, "a1 = %.16f, a2 = %.16f, a3 = %.16f", a1, a2, a3);
if (fabs(p1->x - p3->x) < EPSILON_SQLMM
&& fabs(p1->y - p3->y) < EPSILON_SQLMM)
{
sweep = 2*M_PI;
}
/* Clockwise */
else if (a1 > a2 && a2 > a3)
{
sweep = a3 - a1;
}
/* Counter-clockwise */
else if (a1 < a2 && a2 < a3)
{
sweep = a3 - a1;
}
/* Clockwise, wrap */
else if ((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3))
{
sweep = a3 - a1 + 2*M_PI;
}
/* Counter-clockwise, wrap */
else if ((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3))
{
sweep = a3 - a1 - 2*M_PI;
}
else
{
sweep = 0.0;
}
ptcount = ceil(fabs(perQuad * sweep / M_PI_2));
result = ptarray_construct(1, 1, ptcount);
increment = M_PI_2 / perQuad;
if (sweep < 0) increment *= -1.0;
angle = a1;
LWDEBUGF(3, "ptcount: %d, perQuad: %d, sweep: %.16f, increment: %.16f", ptcount, perQuad, sweep, increment);
for (i = 0; i < ptcount - 1; i++)
{
pt = getPoint_internal(result, i);
angle += increment;
if (increment > 0.0 && angle > M_PI) angle -= 2*M_PI;
if (increment < 0.0 && angle < -1*M_PI) angle -= 2*M_PI;
pbuf.x = center->x + radius*cos(angle);
pbuf.y = center->y + radius*sin(angle);
if ((sweep > 0 && angle < a2) || (sweep < 0 && angle > a2))
{
if ((sweep > 0 && a1 < a2) || (sweep < 0 && a1 > a2))
{
pbuf.z = interpolate_arc_old(angle, p1->z, a1, p2->z, a2);
pbuf.m = interpolate_arc_old(angle, p1->m, a1, p2->m, a2);
}
else
{
if (sweep > 0)
{
pbuf.z = interpolate_arc_old(angle, p1->z, a1-(2*M_PI), p2->z, a2);
pbuf.m = interpolate_arc_old(angle, p1->m, a1-(2*M_PI), p2->m, a2);
}
else
{
pbuf.z = interpolate_arc_old(angle, p1->z, a1+(2*M_PI), p2->z, a2);
pbuf.m = interpolate_arc_old(angle, p1->m, a1+(2*M_PI), p2->m, a2);
}
}
}
else
{
if ((sweep > 0 && a2 < a3) || (sweep < 0 && a2 > a3))
{
pbuf.z = interpolate_arc_old(angle, p2->z, a2, p3->z, a3);
pbuf.m = interpolate_arc_old(angle, p2->m, a2, p3->m, a3);
}
else
{
if (sweep > 0)
{
pbuf.z = interpolate_arc_old(angle, p2->z, a2-(2*M_PI), p3->z, a3);
pbuf.m = interpolate_arc_old(angle, p2->m, a2-(2*M_PI), p3->m, a3);
}
else
{
pbuf.z = interpolate_arc_old(angle, p2->z, a2+(2*M_PI), p3->z, a3);
pbuf.m = interpolate_arc_old(angle, p2->m, a2+(2*M_PI), p3->m, a3);
}
}
}
memcpy(pt, (POINT4D *)&pbuf, ptsize);
}
pt = getPoint_internal(result, ptcount - 1);
memcpy(pt, (POINT4D *)p3, ptsize);
lwfree(center);
return result;
}
LWLINE *
lwcurve_segmentize_old(LWCIRCSTRING *icurve, uint32_t perQuad)
{
LWLINE *oline;
POINTARRAY *ptarray;
POINTARRAY *tmp;
uint32_t i, j;
POINT4D *p1 = lwalloc(sizeof(POINT4D));
POINT4D *p2 = lwalloc(sizeof(POINT4D));
POINT4D *p3 = lwalloc(sizeof(POINT4D));
POINT4D *p4 = lwalloc(sizeof(POINT4D));
LWDEBUGF(2, "lwcurve_segmentize_old called., dim = %d", icurve->points->dims);
ptarray = ptarray_construct_empty(FLAGS_GET_Z(icurve->points->flags), FLAGS_GET_M(icurve->points->flags), icurve->points->npoints);
if (!getPoint4d_p(icurve->points, 0, p4))
{
lwerror("lwcurve_segmentize_old: Cannot extract point.");
}
ptarray_append_point(ptarray, p4, 1);
for (i = 2; i < icurve->points->npoints; i+=2)
{
LWDEBUGF(3, "lwcurve_segmentize_old: arc ending at point %d", i);
getPoint4d_p(icurve->points, i - 2, p1);
getPoint4d_p(icurve->points, i - 1, p2);
getPoint4d_p(icurve->points, i, p3);
tmp = lwcircle_segmentize_old(p1, p2, p3, perQuad);
if (tmp)
{
LWDEBUGF(3, "lwcurve_segmentize_old: generated %d points", tmp->npoints);
for (j = 0; j < tmp->npoints; j++)
{
getPoint4d_p(tmp, j, p4);
ptarray_append_point(ptarray, p4, 1);
}
lwfree(tmp);
}
else
{
LWDEBUG(3, "lwcurve_segmentize_old: points are colinear, returning curve points as line");
for (j = i - 1 ; j <= i ; j++)
{
getPoint4d_p(icurve->points, j, p4);
ptarray_append_point(ptarray, p4, 1);
}
}
}
oline = lwline_construct(icurve->srid, NULL, ptarray_clone(ptarray));
lwfree(p1);
lwfree(p2);
lwfree(p3);
lwfree(p4);
lwfree(ptarray);
return oline;
}
LWLINE *
lwcompound_segmentize_old(LWCOMPOUND *icompound, uint32_t perQuad)
{
LWLINE *oline;
LWGEOM *geom;
POINTARRAY *ptarray = NULL;
LWLINE *tmp = NULL;
uint32_t i, j;
POINT4D *p = NULL;
LWDEBUG(2, "lwcompound_segmentize_old called.");
p = lwalloc(sizeof(POINT4D));
ptarray = ptarray_construct_empty(FLAGS_GET_Z(((POINTARRAY *)icompound->geoms[0]->data)->flags), FLAGS_GET_M(((POINTARRAY *)icompound->geoms[0]->data)->flags), 2);
for (i = 0; i < icompound->ngeoms; i++)
{
geom = icompound->geoms[i];
if (geom->type == CIRCSTRINGTYPE)
{
tmp = lwcurve_segmentize_old((LWCIRCSTRING *)geom, perQuad);
for (j = 0; j < tmp->points->npoints; j++)
{
getPoint4d_p(tmp->points, j, p);
ptarray_append_point(ptarray, p, 0);
}
lwfree(tmp);
}
else if (geom->type == LINETYPE)
{
tmp = (LWLINE *)geom;
for (j = 0; j < tmp->points->npoints; j++)
{
getPoint4d_p(tmp->points, j, p);
ptarray_append_point(ptarray, p, 0);
}
}
else
{
lwerror("Unsupported geometry type %d found.", geom->type);
return NULL;
}
}
oline = lwline_construct(icompound->srid, NULL, ptarray_clone(ptarray));
lwfree(ptarray);
lwfree(p);
return oline;
}
LWPOLY *
lwcurvepoly_segmentize_old(LWCURVEPOLY *curvepoly, uint32_t perQuad)
{
LWPOLY *ogeom;
LWGEOM *tmp;
LWLINE *line;
POINTARRAY **ptarray;
int i;
LWDEBUG(2, "lwcurvepoly_segmentize_old called.");
ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings);
for (i = 0; i < curvepoly->nrings; i++)
{
tmp = curvepoly->rings[i];
if (tmp->type == CIRCSTRINGTYPE)
{
line = lwcurve_segmentize_old((LWCIRCSTRING *)tmp, perQuad);
ptarray[i] = ptarray_clone(line->points);
lwfree(line);
}
else if (tmp->type == LINETYPE)
{
line = (LWLINE *)tmp;
ptarray[i] = ptarray_clone(line->points);
}
else if (tmp->type == COMPOUNDTYPE)
{
line = lwcompound_segmentize_old((LWCOMPOUND *)tmp, perQuad);
ptarray[i] = ptarray_clone(line->points);
lwfree(line);
}
else
{
lwerror("Invalid ring type found in CurvePoly.");
return NULL;
}
}
ogeom = lwpoly_construct(curvepoly->srid, NULL, curvepoly->nrings, ptarray);
return ogeom;
}
LWMLINE *
lwmcurve_segmentize_old(LWMCURVE *mcurve, uint32_t perQuad)
{
LWMLINE *ogeom;
LWGEOM *tmp;
LWGEOM **lines;
int i;
LWDEBUGF(2, "lwmcurve_segmentize_old called, geoms=%d, dim=%d.", mcurve->ngeoms, TYPE_NDIMS(mcurve->type));
lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms);
for (i = 0; i < mcurve->ngeoms; i++)
{
tmp = mcurve->geoms[i];
if (tmp->type == CIRCSTRINGTYPE)
{
lines[i] = (LWGEOM *)lwcurve_segmentize_old((LWCIRCSTRING *)tmp, perQuad);
}
else if (tmp->type == LINETYPE)
{
lines[i] = (LWGEOM *)lwline_construct(mcurve->srid, NULL, ptarray_clone(((LWLINE *)tmp)->points));
}
else
{
lwerror("Unsupported geometry found in MultiCurve.");
return NULL;
}
}
ogeom = (LWMLINE *)lwcollection_construct(MULTILINETYPE, mcurve->srid, NULL, mcurve->ngeoms, lines);
return ogeom;
}
LWMPOLY *
lwmsurface_segmentize_old(LWMSURFACE *msurface, uint32_t perQuad)
{
LWMPOLY *ogeom;
LWGEOM *tmp;
LWPOLY *poly;
LWGEOM **polys;
POINTARRAY **ptarray;
int i, j;
LWDEBUG(2, "lwmsurface_segmentize_old called.");
polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms);
for (i = 0; i < msurface->ngeoms; i++)
{
tmp = msurface->geoms[i];
if (tmp->type == CURVEPOLYTYPE)
{
polys[i] = (LWGEOM *)lwcurvepoly_segmentize_old((LWCURVEPOLY *)tmp, perQuad);
}
else if (tmp->type == POLYGONTYPE)
{
poly = (LWPOLY *)tmp;
ptarray = lwalloc(sizeof(POINTARRAY *)*poly->nrings);
for (j = 0; j < poly->nrings; j++)
{
ptarray[j] = ptarray_clone(poly->rings[j]);
}
polys[i] = (LWGEOM *)lwpoly_construct(msurface->srid, NULL, poly->nrings, ptarray);
}
}
ogeom = (LWMPOLY *)lwcollection_construct(MULTIPOLYGONTYPE, msurface->srid, NULL, msurface->ngeoms, polys);
return ogeom;
}
LWCOLLECTION *
lwcollection_segmentize_old(LWCOLLECTION *collection, uint32_t perQuad)
{
LWCOLLECTION *ocol;
LWGEOM *tmp;
LWGEOM **geoms;
int i;
LWDEBUG(2, "lwcollection_segmentize_old called.");
geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms);
for (i=0; i<collection->ngeoms; i++)
{
tmp = collection->geoms[i];
switch (tmp->type)
{
case CIRCSTRINGTYPE:
geoms[i] = (LWGEOM *)lwcurve_segmentize_old((LWCIRCSTRING *)tmp, perQuad);
break;
case COMPOUNDTYPE:
geoms[i] = (LWGEOM *)lwcompound_segmentize_old((LWCOMPOUND *)tmp, perQuad);
break;
case CURVEPOLYTYPE:
geoms[i] = (LWGEOM *)lwcurvepoly_segmentize_old((LWCURVEPOLY *)tmp, perQuad);
break;
case COLLECTIONTYPE:
geoms[i] = (LWGEOM *)lwcollection_segmentize_old((LWCOLLECTION *)tmp, perQuad);
break;
default:
geoms[i] = lwgeom_clone(tmp);
break;
}
}
ocol = lwcollection_construct(COLLECTIONTYPE, collection->srid, NULL, collection->ngeoms, geoms);
return ocol;
}
LWGEOM *
lwgeom_segmentize_old(LWGEOM *geom, uint32_t perQuad)
{
LWGEOM * ogeom = NULL;
switch (geom->type)
{
case CIRCSTRINGTYPE:
ogeom = (LWGEOM *)lwcurve_segmentize_old((LWCIRCSTRING *)geom, perQuad);
break;
case COMPOUNDTYPE:
ogeom = (LWGEOM *)lwcompound_segmentize_old((LWCOMPOUND *)geom, perQuad);
break;
case CURVEPOLYTYPE:
ogeom = (LWGEOM *)lwcurvepoly_segmentize_old((LWCURVEPOLY *)geom, perQuad);
break;
case MULTICURVETYPE:
ogeom = (LWGEOM *)lwmcurve_segmentize_old((LWMCURVE *)geom, perQuad);
break;
case MULTISURFACETYPE:
ogeom = (LWGEOM *)lwmsurface_segmentize_old((LWMSURFACE *)geom, perQuad);
break;
case COLLECTIONTYPE:
ogeom = (LWGEOM *)lwcollection_segmentize_old((LWCOLLECTION *)geom, perQuad);
break;
default:
ogeom = lwgeom_clone(geom);
}
return ogeom;
}
/*******************************************************************************
* END ST_CurveToLine 1.5.3 BEHAVIOR
******************************************************************************/
# create function in database as
CREATE OR REPLACE FUNCTION ST_CurveToLine_1_5_3(Geometry geometry)
RETURNS geometry
AS '$libdir/postgis-2.0', 'LWGEOM_curve_segmentize_old'
LANGUAGE 'c' IMMUTABLE STRICT;
More information about the postgis-users
mailing list