From: Mark Cave-Ayland Date: Thu, 22 May 2008 20:43:00 +0000 (+0000) Subject: Since PGXS compiles libraries with -Wall, attempt to remove as many warnings as possi... X-Git-Tag: 1.4.0b1~904 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=5603cdfaf285a6db62e9038a5e5c5371fdd79cb3;p=postgis Since PGXS compiles libraries with -Wall, attempt to remove as many warnings as possible. Most of these are missing function prototypes at the top of each file. git-svn-id: http://svn.osgeo.org/postgis/trunk@2781 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/lwgeom/liblwgeom.h b/lwgeom/liblwgeom.h index 8d2b495aa..0ae0c0406 100644 --- a/lwgeom/liblwgeom.h +++ b/lwgeom/liblwgeom.h @@ -912,6 +912,7 @@ extern float nextUp_f(double d); extern double nextDown_d(float d); extern double nextUp_d(float d); +extern float nextafterf_custom(float x, float y); #define LW_MAX(a,b) ((a) > (b) ? (a) : (b)) @@ -1070,6 +1071,7 @@ extern int32 lwgeom_nrings_recursive(uchar *serialized); extern void ptarray_reverse(POINTARRAY *pa); extern POINTARRAY *ptarray_substring(POINTARRAY *, double, double); extern double ptarray_locate_point(POINTARRAY *, POINT2D *); +extern void closest_point_on_segment(POINT2D *p, POINT2D *A, POINT2D *B, POINT2D *ret); /* * Ensure every segment is at most 'dist' long. @@ -1163,6 +1165,10 @@ typedef struct * LWCURVE functions ******************************************************************/ +/* Casts LWGEOM->LW* (return NULL if cast is illegal) */ +extern LWCURVE *lwgeom_as_lwcurve(LWGEOM *lwgeom); + + LWCURVE *lwcurve_construct(int SRID, BOX2DFLOAT4 *bbox, POINTARRAY *points); /* diff --git a/lwgeom/long_xact.c b/lwgeom/long_xact.c index 4570154a8..96f9ee208 100644 --- a/lwgeom/long_xact.c +++ b/lwgeom/long_xact.c @@ -8,6 +8,7 @@ #define ABORT_ON_AUTH_FAILURE 1 Datum check_authorization(PG_FUNCTION_ARGS); +Datum getTransactionID(PG_FUNCTION_ARGS); /* * This trigger will check for authorization before diff --git a/lwgeom/lwcurve.c b/lwgeom/lwcurve.c index 501107174..a5dc40093 100644 --- a/lwgeom/lwcurve.c +++ b/lwgeom/lwcurve.c @@ -1,750 +1,762 @@ -/********************************************************************** - * $Id$ - * - * PostGIS - Spatial Types for PostgreSQL - * http://postgis.refractions.net - * Copyright 2001-2006 Refractions Research Inc. - * - * This is free software; you can redistribute and/or modify it under - * the terms of the GNU General Public Licence. See the COPYING file. - * - **********************************************************************/ - -/* basic LWCURVE functions */ - -#include -#include -#include -#include -#include "liblwgeom.h" - -/*#define PGIS_DEBUG_CALLS 1 */ -/*#define PGIS_DEBUG 1 */ - -#ifndef MAXFLOAT - #define MAXFLOAT 3.402823466e+38F -#endif - -/* - * Construct a new LWCURVE. points will *NOT* be copied - * use SRID=-1 for unknown SRID (will have 8bit type's S = 0) - */ -LWCURVE * -lwcurve_construct(int SRID, BOX2DFLOAT4 *bbox, POINTARRAY *points) -{ - LWCURVE *result; - - /* - * The first arc requires three points. Each additional - * arc requires two more points. Thus the minimum point count - * is three, and the count must be odd. - */ - if(points->npoints % 2 != 1 || points->npoints < 3) - { - lwerror("lwcurve_construct: invalid point count %d", points->npoints); - return NULL; - } - - result = (LWCURVE*) lwalloc(sizeof(LWCURVE)); - - result->type = lwgeom_makeType_full( - TYPE_HASZ(points->dims), - TYPE_HASM(points->dims), - (SRID!=-1), CURVETYPE, 0); - result->SRID = SRID; - result->points = points; - result->bbox = bbox; - - return result; -} - -/* - * given the LWGEOM serialized form (or a point into a multi* one) - * construct a propert LWCURVE. - * serialized_form should point to the 8bit type format (with type = 8) - * See serialized form doc - */ -LWCURVE * -lwcurve_deserialize(uchar *serialized_form) -{ - uchar type; - LWCURVE *result; - uchar *loc=NULL; - uint32 npoints; - POINTARRAY *pa; - - type = (uchar)serialized_form[0]; - if(lwgeom_getType(type) != CURVETYPE) - { - lwerror("lwcurve_deserialize: attempt to deserialize a curve which is really a %s", lwgeom_typename(type)); - return NULL; - } - - result = (LWCURVE*) lwalloc(sizeof(LWCURVE)); - result->type = type; - - loc = serialized_form + 1; - - if(lwgeom_hasBBOX(type)) - { -#ifdef PGIS_DEBUG - lwnotice("lwcurve_deserialize: input has bbox"); -#endif - result->bbox = lwalloc(sizeof(BOX2DFLOAT4)); - memcpy(result->bbox, loc, sizeof(BOX2DFLOAT4)); - loc += sizeof(BOX2DFLOAT4); - } - else - { -#ifdef PGIS_DEBUG - lwnotice("lwcurve_deserialize: input lacks bbox"); -#endif - result->bbox = NULL; - } - - if(lwgeom_hasSRID(type)) - { -#ifdef PGIS_DEBUG - lwnotice("lwcurve_deserialize: input has srid"); -#endif - result->SRID = lw_get_int32(loc); - loc += 4; /* type + SRID */ - } - else - { -#ifdef PGIS_DEBUG - lwnotice("lwcurve_deserialize: input lacks srid"); -#endif - result->SRID = -1; - } - - /* we've read the type (1 byte) and SRID (4 bytes, if present) */ - - npoints = lw_get_uint32(loc); -#ifdef PGIS_DEBUG - lwnotice("curve npoints = %d", npoints); -#endif - loc += 4; - pa = pointArray_construct(loc, TYPE_HASZ(type), TYPE_HASM(type), npoints); - result->points = pa; - return result; -} - -/* - * convert this curve into its serialized form - * result's first char will be the 8bit type. See serialized form doc - */ -uchar * -lwcurve_serialize(LWCURVE *curve) -{ - size_t size, retsize; - uchar * result; - - if(curve == NULL) { - lwerror("lwcurve_serialize:: given null curve"); - return NULL; - } - - size = lwcurve_serialize_size(curve); - result = lwalloc(size); - lwcurve_serialize_buf(curve, result, &retsize); - if(retsize != size) - lwerror("lwcurve_serialize_size returned %d, ..selialize_buf returned %d", size, retsize); - return result; -} - -/* - * convert this curve into its serialized form writing it into - * the given buffer, and returning number of bytes written into - * the given int pointer. - * result's first char will be the 8bit type. See serialized form doc - */ -void lwcurve_serialize_buf(LWCURVE *curve, uchar *buf, size_t *retsize) -{ - char hasSRID; - uchar *loc; - int ptsize; - size_t size; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_serialize_buf(%p, %p, %p) called", - curve, buf, retsize); -#endif - - if(curve == NULL) - { - lwerror("lwcurve_serialize:: given null curve"); - return; - } - - if(TYPE_GETZM(curve->type) != TYPE_GETZM(curve->points->dims)) - { - lwerror("Dimensions mismatch in lwcurve"); - return; - } - - ptsize = pointArray_ptsize(curve->points); - - hasSRID = (curve->SRID != -1); - - buf[0] = (uchar)lwgeom_makeType_full( - TYPE_HASZ(curve->type), TYPE_HASM(curve->type), - hasSRID, CURVETYPE, curve->bbox ? 1 : 0); - loc = buf+1; - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_serialize_buf added type (%d)", curve->type); -#endif - - if(curve->bbox) - { - memcpy(loc, curve->bbox, sizeof(BOX2DFLOAT4)); - loc += sizeof(BOX2DFLOAT4); - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_serialize_buf added BBOX"); -#endif - - } - - if(hasSRID) - { - memcpy(loc, &curve->SRID, sizeof(int32)); - loc += sizeof(int32); - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_serialize_buf added SRID"); -#endif - - } - - memcpy(loc, &curve->points->npoints, sizeof(uint32)); - loc += sizeof(uint32); - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_serialize_buf added npoints (%d)", - curve->points->npoints); -#endif - - /* copy in points */ - size = curve->points->npoints * ptsize; - memcpy(loc, getPoint_internal(curve->points, 0), size); - loc += size; - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_serialize_buf copied serialized_pointlist (%d bytes)", - ptsize * curve->points->npoints); -#endif - - if(retsize) *retsize = loc-buf; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_serialize_buf returning (loc: %p, size: %d)", - loc, loc-buf); -#endif -} - -/* find length of this deserialized curve */ -size_t -lwcurve_serialize_size(LWCURVE *curve) -{ - size_t size = 1; /* type */ - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_serialize_size called"); -#endif - - if(curve->SRID != -1) size += 4; /* SRID */ - if(curve->bbox) size += sizeof(BOX2DFLOAT4); - - size += 4; /* npoints */ - size += pointArray_ptsize(curve->points) * curve->points->npoints; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_serialize_size returning %d", size); -#endif - - return size; -} - -BOX3D * -lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3) -{ - double x1, x2, y1, y2, z1, z2; - double angle, radius, sweep; - /* - double top, left; - */ - double a1, a2, a3; - double xe = 0.0, ye = 0.0; - POINT4D *center; - int i; - BOX3D *box; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcircle_compute_box3d called."); -#endif - - center = lwalloc(sizeof(POINT4D)); - radius = lwcircle_center(p1, p2, p3, ¢er); - if(radius < 0.0) return NULL; - - /* - top = center->y + radius; - left = center->x - radius; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcircle_compute_box3d: top=%.16f, left=%.16f", top, left); -#endif - */ - - x1 = MAXFLOAT; - x2 = -1 * MAXFLOAT; - y1 = MAXFLOAT; - y2 = -1 * MAXFLOAT; - - 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); - - /* Determine sweep angle */ - 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; - } - -#ifdef PGIS_DEBUG - lwnotice("a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep); -#endif - - angle = 0.0; - for(i=0; i < 6; i++) - { - switch(i) { - case 0: - angle = 0.0; - xe = center->x + radius; - ye = center->y; - break; - case 1: - angle = M_PI_2; - xe = center->x; - ye = center->y + radius; - break; - case 2: - angle = M_PI; - xe = center->x - radius; - ye = center->y; - break; - case 3: - angle = -1 * M_PI_2; - xe = center->x; - ye = center->y - radius; - break; - case 4: - angle = a1; - xe = p1->x; - ye = p1->y; - break; - case 5: - angle = a3; - xe = p3->x; - ye = p3->y; - break; - } - if(i < 4) - { - if(sweep > 0.0 && (angle > a3 || angle < a1)) continue; - if(sweep < 0.0 && (angle < a3 || angle > a1)) continue; - } - -#ifdef PGIS_DEBUG - lwnotice("lwcircle_compute_box3d: potential extreame %d (%.16f, %.16f)", i, xe, ye); -#endif - x1 = (x1 < xe) ? x1 : xe; - y1 = (y1 < ye) ? y1 : ye; - x2 = (x2 > xe) ? x2 : xe; - y2 = (y2 > ye) ? y2 : ye; - } -#ifdef PGIS_DEBUG - lwnotice("lwcircle_compute_box3d: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2); -#endif - - /* - x1 = center->x + x1 * radius; - x2 = center->x + x2 * radius; - y1 = center->y + y1 * radius; - y2 = center->y + y2 * radius; - */ - z1 = (p1->z < p2->z) ? p1->z : p2->z; - z1 = (z1 < p3->z) ? z1 : p3->z; - z2 = (p1->z > p2->z) ? p1->z : p2->z; - z2 = (z2 > p3->z) ? z2 : p3->z; - - box = lwalloc(sizeof(BOX3D)); - box->xmin = x1; box->xmax = x2; - box->ymin = y1; box->ymax = y2; - box->zmin = z1; box->zmax = z2; - - lwfree(center); - - return box; -} - -/* - * Find bounding box (standard one) - * zmin=zmax=NO_Z_VALUE if 2d - * TODO: This ignores curvature, which should be taken into account. - */ -BOX3D * -lwcurve_compute_box3d(LWCURVE *curve) -{ - BOX3D *box, *tmp; - int i; - POINT4D *p1 = lwalloc(sizeof(POINT4D)); - POINT4D *p2 = lwalloc(sizeof(POINT4D)); - POINT4D *p3 = lwalloc(sizeof(POINT4D)); - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_compute_box3d called."); -#endif - - /* initialize box values */ - box = lwalloc(sizeof(BOX3D)); - box->xmin = MAXFLOAT; box->xmax = -1 * MAXFLOAT; - box->ymin = MAXFLOAT; box->ymax = -1 * MAXFLOAT; - box->zmin = MAXFLOAT; box->zmax = -1 * MAXFLOAT; - - for(i = 2; i < curve->points->npoints; i+=2) - { - getPoint4d_p(curve->points, i-2, p1); - getPoint4d_p(curve->points, i-1, p2); - getPoint4d_p(curve->points, i, p3); - tmp = lwcircle_compute_box3d(p1, p2, p3); - if(tmp == NULL) continue; - box->xmin = (box->xmin < tmp->xmin) ? box->xmin : tmp->xmin; - box->xmax = (box->xmax > tmp->xmax) ? box->xmax : tmp->xmax; - box->ymin = (box->ymin < tmp->ymin) ? box->ymin : tmp->ymin; - box->ymax = (box->ymax > tmp->ymax) ? box->ymax : tmp->ymax; - box->zmin = (box->zmin < tmp->zmin) ? box->zmin : tmp->zmin; - box->zmax = (box->zmax > tmp->zmax) ? box->zmax : tmp->zmax; -#ifdef PGIS_DEBUG_CALLS - lwnotice("curve %d x=(%.16f,%.16f) y=(%.16f,%.16f) z=(%.16f,%.16f)", i/2, box->xmin, box->xmax, box->ymin, box->ymax, box->zmin, box->zmax); -#endif - } - - - return box; -} - -int -lwcurve_compute_box2d_p(LWCURVE *curve, BOX2DFLOAT4 *result) -{ - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurve_compute_box2d_p called."); -#endif - - BOX3D *box = lwcurve_compute_box3d(curve); - if(box == NULL) return 0; - box3d_to_box2df_p(box, result); - return 1; -} - -void pfree_curve(LWCURVE *curve) -{ - lwfree(curve->points); - lwfree(curve); -} - -/* find length of this serialized curve */ -size_t -lwgeom_size_curve(const uchar *serialized_curve) -{ - int type = (uchar)serialized_curve[0]; - uint32 result = 1; /* type */ - const uchar *loc; - uint32 npoints; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwgeom_size_curve called"); -#endif - if(lwgeom_getType(type) != CURVETYPE) - lwerror("lwgeom_size_curve::attempt to find the length of a non-curve"); - - loc = serialized_curve + 1; - if(lwgeom_hasBBOX(type)) - { - loc += sizeof(BOX2DFLOAT4); - result += sizeof(BOX2DFLOAT4); - } - - if(lwgeom_hasSRID(type)) - { - loc += 4; /* type + SRID */ - result += 4; - } - - /* we've read the type (1 byte) and SRID (4 bytes, if present) */ - npoints = lw_get_uint32(loc); - result += sizeof(uint32); /* npoints */ - - result += TYPE_NDIMS(type) * sizeof(double) * npoints; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwgeom_size_curve returning %d", result); -#endif - - return result; -} - -void printLWCURVE(LWCURVE *curve) -{ - lwnotice("LWCURVE {"); - lwnotice(" ndims = %i", (int)TYPE_NDIMS(curve->type)); - lwnotice(" SRID = %i", (int)curve->SRID); - printPA(curve->points); - lwnotice("}"); -} - -/* Clone LWCURVE object. POINTARRAY is not copied. */ -LWCURVE * -lwcurve_clone(const LWCURVE *g) -{ - LWCURVE *ret = lwalloc(sizeof(LWCURVE)); - memcpy(ret, g, sizeof(LWCURVE)); - if(g->bbox) ret->bbox = box2d_clone(g->bbox); - return ret; -} - -/* - * Add 'what' to this curve at position 'where'. - * where=0 == prepend - * where=-1 == append - * Returns a MULTICURVE or a GEOMETRYCOLLECTION - */ -LWGEOM * -lwcurve_add(const LWCURVE *to, uint32 where, const LWGEOM *what) -{ - LWCOLLECTION *col; - LWGEOM **geoms; - int newtype; - - if(where != -1 && where != 0) - { - lwerror("lwcurve_add only supports 0 or -1 as second argument %d", where); - return NULL; - } - - /* dimensions compatibility are checked by caller */ - - /* Construct geoms array */ - geoms = lwalloc(sizeof(LWGEOM *)*2); - if(where == -1) /* append */ - { - geoms[0] = lwgeom_clone((LWGEOM *)to); - geoms[1] = lwgeom_clone(what); - } - else /* prepend */ - { - geoms[0] = lwgeom_clone(what); - geoms[1] = lwgeom_clone((LWGEOM *)to); - } - - /* reset SRID and wantbbox flag from component types */ - geoms[0]->SRID = geoms[1]->SRID = -1; - TYPE_SETHASSRID(geoms[0]->type, 0); - TYPE_SETHASSRID(geoms[1]->type, 0); - TYPE_SETHASBBOX(geoms[0]->type, 0); - TYPE_SETHASBBOX(geoms[1]->type, 0); - - /* Find appropriate geom type */ - if(TYPE_GETTYPE(what->type) == CURVETYPE || TYPE_GETTYPE(what->type) == LINETYPE) newtype = MULTICURVETYPE; - else newtype = COLLECTIONTYPE; - - col = lwcollection_construct(newtype, - to->SRID, NULL, - 2, geoms); - - return (LWGEOM *)col; -} - -void lwcurve_reverse(LWCURVE *curve) -{ - ptarray_reverse(curve->points); -} - -/* - * TODO: Invalid segmentization. This should accomodate the curvature. - */ -LWCURVE * -lwcurve_segmentize2d(LWCURVE *curve, double dist) -{ - return lwcurve_construct(curve->SRID, NULL, - ptarray_segmentize2d(curve->points, dist)); -} - -/* check coordinate equality */ -char -lwcurve_same(const LWCURVE *me, const LWCURVE *you) -{ - return ptarray_same(me->points, you->points); -} - -/* - * Construct a LWCURVE from an array of LWPOINTs - * LWCURVE dimensions are large enough to host all input dimensions. - */ -LWCURVE * -lwcurve_from_lwpointarray(int SRID, unsigned int npoints, LWPOINT **points) -{ - int zmflag=0; - unsigned int i; - POINTARRAY *pa; - uchar *newpoints, *ptr; - size_t ptsize, size; - - /* - * Find output dimensions, check integrity - */ - for(i = 0; i < npoints; i++) - { - if(TYPE_GETTYPE(points[i]->type) != POINTTYPE) - { - lwerror("lwcurve_from_lwpointarray: invalid input type: %s", - lwgeom_typename(TYPE_GETTYPE(points[i]->type))); - return NULL; - } - if(TYPE_HASZ(points[i]->type)) zmflag |= 2; - if(TYPE_HASM(points[i]->type)) zmflag |=1; - if(zmflag == 3) break; - } - - if(zmflag == 0) ptsize = 2 * sizeof(double); - else if(zmflag == 3) ptsize = 4 * sizeof(double); - else ptsize = 3 * sizeof(double); - - /* - * Allocate output points array - */ - size = ptsize * npoints; - newpoints = lwalloc(size); - memset(newpoints, 0, size); - - ptr = newpoints; - for(i = 0; i < npoints; i++) - { - size = pointArray_ptsize(points[i]->point); - memcpy(ptr, getPoint_internal(points[i]->point, 0), size); - ptr += ptsize; - } - pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, npoints); - - return lwcurve_construct(SRID, NULL, pa); -} - -/* - * Construct a LWCURVE from a LWMPOINT - */ -LWCURVE * -lwcurve_from_lwmpoint(int SRID, LWMPOINT *mpoint) -{ - unsigned int i; - POINTARRAY *pa; - char zmflag = TYPE_GETZM(mpoint->type); - size_t ptsize, size; - uchar *newpoints, *ptr; - - if(zmflag == 0) ptsize = 2 * sizeof(double); - else if(zmflag == 3) ptsize = 4 * sizeof(double); - else ptsize = 3 * sizeof(double); - - /* Allocate space for output points */ - size = ptsize * mpoint->ngeoms; - newpoints = lwalloc(size); - memset(newpoints, 0, size); - - ptr = newpoints; - for(i = 0; i < mpoint->ngeoms; i++) - { - memcpy(ptr, - getPoint_internal(mpoint->geoms[i]->point, 0), - ptsize); - ptr += ptsize; - } - - pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, - mpoint->ngeoms); - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_from_lwmpoint: constructed pointarray for %d points, %d zmflag", mpoint->ngeoms, zmflag); -#endif - - return lwcurve_construct(SRID, NULL, pa); -} - -LWCURVE * -lwcurve_addpoint(LWCURVE *curve, LWPOINT *point, unsigned int where) -{ - POINTARRAY *newpa; - LWCURVE *ret; - - newpa = ptarray_addPoint(curve->points, - getPoint_internal(point->point, 0), - TYPE_NDIMS(point->type), where); - ret = lwcurve_construct(curve->SRID, NULL, newpa); - - return ret; -} - -LWCURVE * -lwcurve_removepoint(LWCURVE *curve, unsigned int index) -{ - POINTARRAY *newpa; - LWCURVE *ret; - - newpa = ptarray_removePoint(curve->points, index); - ret = lwcurve_construct(curve->SRID, NULL, newpa); - - return ret; -} - -/* - * Note: input will be changed, make sure you have permissions for this. - * */ -void -lwcurve_setPoint4d(LWCURVE *curve, unsigned int index, POINT4D *newpoint) -{ - setPoint4d(curve->points, index, newpoint); -} - - - - - - - - - - - - +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2001-2006 Refractions Research Inc. + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +/* basic LWCURVE functions */ + +#include +#include +#include +#include +#include "liblwgeom.h" + +BOX3D *lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3); +void printLWCURVE(LWCURVE *curve); +void lwcurve_reverse(LWCURVE *curve); +LWCURVE *lwcurve_segmentize2d(LWCURVE *curve, double dist); +char lwcurve_same(const LWCURVE *me, const LWCURVE *you); +LWCURVE *lwcurve_from_lwpointarray(int SRID, unsigned int npoints, LWPOINT **points); +LWCURVE *lwcurve_from_lwmpoint(int SRID, LWMPOINT *mpoint); +LWCURVE *lwcurve_addpoint(LWCURVE *curve, LWPOINT *point, unsigned int where); +LWCURVE *lwcurve_removepoint(LWCURVE *curve, unsigned int index); +void lwcurve_setPoint4d(LWCURVE *curve, unsigned int index, POINT4D *newpoint); + + +/*#define PGIS_DEBUG_CALLS 1 */ +/*#define PGIS_DEBUG 1 */ + +#ifndef MAXFLOAT + #define MAXFLOAT 3.402823466e+38F +#endif + +/* + * Construct a new LWCURVE. points will *NOT* be copied + * use SRID=-1 for unknown SRID (will have 8bit type's S = 0) + */ +LWCURVE * +lwcurve_construct(int SRID, BOX2DFLOAT4 *bbox, POINTARRAY *points) +{ + LWCURVE *result; + + /* + * The first arc requires three points. Each additional + * arc requires two more points. Thus the minimum point count + * is three, and the count must be odd. + */ + if(points->npoints % 2 != 1 || points->npoints < 3) + { + lwerror("lwcurve_construct: invalid point count %d", points->npoints); + return NULL; + } + + result = (LWCURVE*) lwalloc(sizeof(LWCURVE)); + + result->type = lwgeom_makeType_full( + TYPE_HASZ(points->dims), + TYPE_HASM(points->dims), + (SRID!=-1), CURVETYPE, 0); + result->SRID = SRID; + result->points = points; + result->bbox = bbox; + + return result; +} + +/* + * given the LWGEOM serialized form (or a point into a multi* one) + * construct a propert LWCURVE. + * serialized_form should point to the 8bit type format (with type = 8) + * See serialized form doc + */ +LWCURVE * +lwcurve_deserialize(uchar *serialized_form) +{ + uchar type; + LWCURVE *result; + uchar *loc=NULL; + uint32 npoints; + POINTARRAY *pa; + + type = (uchar)serialized_form[0]; + if(lwgeom_getType(type) != CURVETYPE) + { + lwerror("lwcurve_deserialize: attempt to deserialize a curve which is really a %s", lwgeom_typename(type)); + return NULL; + } + + result = (LWCURVE*) lwalloc(sizeof(LWCURVE)); + result->type = type; + + loc = serialized_form + 1; + + if(lwgeom_hasBBOX(type)) + { +#ifdef PGIS_DEBUG + lwnotice("lwcurve_deserialize: input has bbox"); +#endif + result->bbox = lwalloc(sizeof(BOX2DFLOAT4)); + memcpy(result->bbox, loc, sizeof(BOX2DFLOAT4)); + loc += sizeof(BOX2DFLOAT4); + } + else + { +#ifdef PGIS_DEBUG + lwnotice("lwcurve_deserialize: input lacks bbox"); +#endif + result->bbox = NULL; + } + + if(lwgeom_hasSRID(type)) + { +#ifdef PGIS_DEBUG + lwnotice("lwcurve_deserialize: input has srid"); +#endif + result->SRID = lw_get_int32(loc); + loc += 4; /* type + SRID */ + } + else + { +#ifdef PGIS_DEBUG + lwnotice("lwcurve_deserialize: input lacks srid"); +#endif + result->SRID = -1; + } + + /* we've read the type (1 byte) and SRID (4 bytes, if present) */ + + npoints = lw_get_uint32(loc); +#ifdef PGIS_DEBUG + lwnotice("curve npoints = %d", npoints); +#endif + loc += 4; + pa = pointArray_construct(loc, TYPE_HASZ(type), TYPE_HASM(type), npoints); + result->points = pa; + return result; +} + +/* + * convert this curve into its serialized form + * result's first char will be the 8bit type. See serialized form doc + */ +uchar * +lwcurve_serialize(LWCURVE *curve) +{ + size_t size, retsize; + uchar * result; + + if(curve == NULL) { + lwerror("lwcurve_serialize:: given null curve"); + return NULL; + } + + size = lwcurve_serialize_size(curve); + result = lwalloc(size); + lwcurve_serialize_buf(curve, result, &retsize); + if(retsize != size) + lwerror("lwcurve_serialize_size returned %d, ..selialize_buf returned %d", size, retsize); + return result; +} + +/* + * convert this curve into its serialized form writing it into + * the given buffer, and returning number of bytes written into + * the given int pointer. + * result's first char will be the 8bit type. See serialized form doc + */ +void lwcurve_serialize_buf(LWCURVE *curve, uchar *buf, size_t *retsize) +{ + char hasSRID; + uchar *loc; + int ptsize; + size_t size; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_serialize_buf(%p, %p, %p) called", + curve, buf, retsize); +#endif + + if(curve == NULL) + { + lwerror("lwcurve_serialize:: given null curve"); + return; + } + + if(TYPE_GETZM(curve->type) != TYPE_GETZM(curve->points->dims)) + { + lwerror("Dimensions mismatch in lwcurve"); + return; + } + + ptsize = pointArray_ptsize(curve->points); + + hasSRID = (curve->SRID != -1); + + buf[0] = (uchar)lwgeom_makeType_full( + TYPE_HASZ(curve->type), TYPE_HASM(curve->type), + hasSRID, CURVETYPE, curve->bbox ? 1 : 0); + loc = buf+1; + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_serialize_buf added type (%d)", curve->type); +#endif + + if(curve->bbox) + { + memcpy(loc, curve->bbox, sizeof(BOX2DFLOAT4)); + loc += sizeof(BOX2DFLOAT4); + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_serialize_buf added BBOX"); +#endif + + } + + if(hasSRID) + { + memcpy(loc, &curve->SRID, sizeof(int32)); + loc += sizeof(int32); + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_serialize_buf added SRID"); +#endif + + } + + memcpy(loc, &curve->points->npoints, sizeof(uint32)); + loc += sizeof(uint32); + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_serialize_buf added npoints (%d)", + curve->points->npoints); +#endif + + /* copy in points */ + size = curve->points->npoints * ptsize; + memcpy(loc, getPoint_internal(curve->points, 0), size); + loc += size; + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_serialize_buf copied serialized_pointlist (%d bytes)", + ptsize * curve->points->npoints); +#endif + + if(retsize) *retsize = loc-buf; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_serialize_buf returning (loc: %p, size: %d)", + loc, loc-buf); +#endif +} + +/* find length of this deserialized curve */ +size_t +lwcurve_serialize_size(LWCURVE *curve) +{ + size_t size = 1; /* type */ + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_serialize_size called"); +#endif + + if(curve->SRID != -1) size += 4; /* SRID */ + if(curve->bbox) size += sizeof(BOX2DFLOAT4); + + size += 4; /* npoints */ + size += pointArray_ptsize(curve->points) * curve->points->npoints; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_serialize_size returning %d", size); +#endif + + return size; +} + +BOX3D * +lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3) +{ + double x1, x2, y1, y2, z1, z2; + double angle, radius, sweep; + /* + double top, left; + */ + double a1, a2, a3; + double xe = 0.0, ye = 0.0; + POINT4D *center; + int i; + BOX3D *box; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcircle_compute_box3d called."); +#endif + + center = lwalloc(sizeof(POINT4D)); + radius = lwcircle_center(p1, p2, p3, ¢er); + if(radius < 0.0) return NULL; + + /* + top = center->y + radius; + left = center->x - radius; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcircle_compute_box3d: top=%.16f, left=%.16f", top, left); +#endif + */ + + x1 = MAXFLOAT; + x2 = -1 * MAXFLOAT; + y1 = MAXFLOAT; + y2 = -1 * MAXFLOAT; + + 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); + + /* Determine sweep angle */ + 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; + } + +#ifdef PGIS_DEBUG + lwnotice("a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep); +#endif + + angle = 0.0; + for(i=0; i < 6; i++) + { + switch(i) { + case 0: + angle = 0.0; + xe = center->x + radius; + ye = center->y; + break; + case 1: + angle = M_PI_2; + xe = center->x; + ye = center->y + radius; + break; + case 2: + angle = M_PI; + xe = center->x - radius; + ye = center->y; + break; + case 3: + angle = -1 * M_PI_2; + xe = center->x; + ye = center->y - radius; + break; + case 4: + angle = a1; + xe = p1->x; + ye = p1->y; + break; + case 5: + angle = a3; + xe = p3->x; + ye = p3->y; + break; + } + if(i < 4) + { + if(sweep > 0.0 && (angle > a3 || angle < a1)) continue; + if(sweep < 0.0 && (angle < a3 || angle > a1)) continue; + } + +#ifdef PGIS_DEBUG + lwnotice("lwcircle_compute_box3d: potential extreame %d (%.16f, %.16f)", i, xe, ye); +#endif + x1 = (x1 < xe) ? x1 : xe; + y1 = (y1 < ye) ? y1 : ye; + x2 = (x2 > xe) ? x2 : xe; + y2 = (y2 > ye) ? y2 : ye; + } +#ifdef PGIS_DEBUG + lwnotice("lwcircle_compute_box3d: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2); +#endif + + /* + x1 = center->x + x1 * radius; + x2 = center->x + x2 * radius; + y1 = center->y + y1 * radius; + y2 = center->y + y2 * radius; + */ + z1 = (p1->z < p2->z) ? p1->z : p2->z; + z1 = (z1 < p3->z) ? z1 : p3->z; + z2 = (p1->z > p2->z) ? p1->z : p2->z; + z2 = (z2 > p3->z) ? z2 : p3->z; + + box = lwalloc(sizeof(BOX3D)); + box->xmin = x1; box->xmax = x2; + box->ymin = y1; box->ymax = y2; + box->zmin = z1; box->zmax = z2; + + lwfree(center); + + return box; +} + +/* + * Find bounding box (standard one) + * zmin=zmax=NO_Z_VALUE if 2d + * TODO: This ignores curvature, which should be taken into account. + */ +BOX3D * +lwcurve_compute_box3d(LWCURVE *curve) +{ + BOX3D *box, *tmp; + int i; + POINT4D *p1 = lwalloc(sizeof(POINT4D)); + POINT4D *p2 = lwalloc(sizeof(POINT4D)); + POINT4D *p3 = lwalloc(sizeof(POINT4D)); + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_compute_box3d called."); +#endif + + /* initialize box values */ + box = lwalloc(sizeof(BOX3D)); + box->xmin = MAXFLOAT; box->xmax = -1 * MAXFLOAT; + box->ymin = MAXFLOAT; box->ymax = -1 * MAXFLOAT; + box->zmin = MAXFLOAT; box->zmax = -1 * MAXFLOAT; + + for(i = 2; i < curve->points->npoints; i+=2) + { + getPoint4d_p(curve->points, i-2, p1); + getPoint4d_p(curve->points, i-1, p2); + getPoint4d_p(curve->points, i, p3); + tmp = lwcircle_compute_box3d(p1, p2, p3); + if(tmp == NULL) continue; + box->xmin = (box->xmin < tmp->xmin) ? box->xmin : tmp->xmin; + box->xmax = (box->xmax > tmp->xmax) ? box->xmax : tmp->xmax; + box->ymin = (box->ymin < tmp->ymin) ? box->ymin : tmp->ymin; + box->ymax = (box->ymax > tmp->ymax) ? box->ymax : tmp->ymax; + box->zmin = (box->zmin < tmp->zmin) ? box->zmin : tmp->zmin; + box->zmax = (box->zmax > tmp->zmax) ? box->zmax : tmp->zmax; +#ifdef PGIS_DEBUG_CALLS + lwnotice("curve %d x=(%.16f,%.16f) y=(%.16f,%.16f) z=(%.16f,%.16f)", i/2, box->xmin, box->xmax, box->ymin, box->ymax, box->zmin, box->zmax); +#endif + } + + + return box; +} + +int +lwcurve_compute_box2d_p(LWCURVE *curve, BOX2DFLOAT4 *result) +{ + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurve_compute_box2d_p called."); +#endif + + BOX3D *box = lwcurve_compute_box3d(curve); + if(box == NULL) return 0; + box3d_to_box2df_p(box, result); + return 1; +} + +void pfree_curve(LWCURVE *curve) +{ + lwfree(curve->points); + lwfree(curve); +} + +/* find length of this serialized curve */ +size_t +lwgeom_size_curve(const uchar *serialized_curve) +{ + int type = (uchar)serialized_curve[0]; + uint32 result = 1; /* type */ + const uchar *loc; + uint32 npoints; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwgeom_size_curve called"); +#endif + if(lwgeom_getType(type) != CURVETYPE) + lwerror("lwgeom_size_curve::attempt to find the length of a non-curve"); + + loc = serialized_curve + 1; + if(lwgeom_hasBBOX(type)) + { + loc += sizeof(BOX2DFLOAT4); + result += sizeof(BOX2DFLOAT4); + } + + if(lwgeom_hasSRID(type)) + { + loc += 4; /* type + SRID */ + result += 4; + } + + /* we've read the type (1 byte) and SRID (4 bytes, if present) */ + npoints = lw_get_uint32(loc); + result += sizeof(uint32); /* npoints */ + + result += TYPE_NDIMS(type) * sizeof(double) * npoints; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwgeom_size_curve returning %d", result); +#endif + + return result; +} + +void printLWCURVE(LWCURVE *curve) +{ + lwnotice("LWCURVE {"); + lwnotice(" ndims = %i", (int)TYPE_NDIMS(curve->type)); + lwnotice(" SRID = %i", (int)curve->SRID); + printPA(curve->points); + lwnotice("}"); +} + +/* Clone LWCURVE object. POINTARRAY is not copied. */ +LWCURVE * +lwcurve_clone(const LWCURVE *g) +{ + LWCURVE *ret = lwalloc(sizeof(LWCURVE)); + memcpy(ret, g, sizeof(LWCURVE)); + if(g->bbox) ret->bbox = box2d_clone(g->bbox); + return ret; +} + +/* + * Add 'what' to this curve at position 'where'. + * where=0 == prepend + * where=-1 == append + * Returns a MULTICURVE or a GEOMETRYCOLLECTION + */ +LWGEOM * +lwcurve_add(const LWCURVE *to, uint32 where, const LWGEOM *what) +{ + LWCOLLECTION *col; + LWGEOM **geoms; + int newtype; + + if(where != -1 && where != 0) + { + lwerror("lwcurve_add only supports 0 or -1 as second argument %d", where); + return NULL; + } + + /* dimensions compatibility are checked by caller */ + + /* Construct geoms array */ + geoms = lwalloc(sizeof(LWGEOM *)*2); + if(where == -1) /* append */ + { + geoms[0] = lwgeom_clone((LWGEOM *)to); + geoms[1] = lwgeom_clone(what); + } + else /* prepend */ + { + geoms[0] = lwgeom_clone(what); + geoms[1] = lwgeom_clone((LWGEOM *)to); + } + + /* reset SRID and wantbbox flag from component types */ + geoms[0]->SRID = geoms[1]->SRID = -1; + TYPE_SETHASSRID(geoms[0]->type, 0); + TYPE_SETHASSRID(geoms[1]->type, 0); + TYPE_SETHASBBOX(geoms[0]->type, 0); + TYPE_SETHASBBOX(geoms[1]->type, 0); + + /* Find appropriate geom type */ + if(TYPE_GETTYPE(what->type) == CURVETYPE || TYPE_GETTYPE(what->type) == LINETYPE) newtype = MULTICURVETYPE; + else newtype = COLLECTIONTYPE; + + col = lwcollection_construct(newtype, + to->SRID, NULL, + 2, geoms); + + return (LWGEOM *)col; +} + +void lwcurve_reverse(LWCURVE *curve) +{ + ptarray_reverse(curve->points); +} + +/* + * TODO: Invalid segmentization. This should accomodate the curvature. + */ +LWCURVE * +lwcurve_segmentize2d(LWCURVE *curve, double dist) +{ + return lwcurve_construct(curve->SRID, NULL, + ptarray_segmentize2d(curve->points, dist)); +} + +/* check coordinate equality */ +char +lwcurve_same(const LWCURVE *me, const LWCURVE *you) +{ + return ptarray_same(me->points, you->points); +} + +/* + * Construct a LWCURVE from an array of LWPOINTs + * LWCURVE dimensions are large enough to host all input dimensions. + */ +LWCURVE * +lwcurve_from_lwpointarray(int SRID, unsigned int npoints, LWPOINT **points) +{ + int zmflag=0; + unsigned int i; + POINTARRAY *pa; + uchar *newpoints, *ptr; + size_t ptsize, size; + + /* + * Find output dimensions, check integrity + */ + for(i = 0; i < npoints; i++) + { + if(TYPE_GETTYPE(points[i]->type) != POINTTYPE) + { + lwerror("lwcurve_from_lwpointarray: invalid input type: %s", + lwgeom_typename(TYPE_GETTYPE(points[i]->type))); + return NULL; + } + if(TYPE_HASZ(points[i]->type)) zmflag |= 2; + if(TYPE_HASM(points[i]->type)) zmflag |=1; + if(zmflag == 3) break; + } + + if(zmflag == 0) ptsize = 2 * sizeof(double); + else if(zmflag == 3) ptsize = 4 * sizeof(double); + else ptsize = 3 * sizeof(double); + + /* + * Allocate output points array + */ + size = ptsize * npoints; + newpoints = lwalloc(size); + memset(newpoints, 0, size); + + ptr = newpoints; + for(i = 0; i < npoints; i++) + { + size = pointArray_ptsize(points[i]->point); + memcpy(ptr, getPoint_internal(points[i]->point, 0), size); + ptr += ptsize; + } + pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, npoints); + + return lwcurve_construct(SRID, NULL, pa); +} + +/* + * Construct a LWCURVE from a LWMPOINT + */ +LWCURVE * +lwcurve_from_lwmpoint(int SRID, LWMPOINT *mpoint) +{ + unsigned int i; + POINTARRAY *pa; + char zmflag = TYPE_GETZM(mpoint->type); + size_t ptsize, size; + uchar *newpoints, *ptr; + + if(zmflag == 0) ptsize = 2 * sizeof(double); + else if(zmflag == 3) ptsize = 4 * sizeof(double); + else ptsize = 3 * sizeof(double); + + /* Allocate space for output points */ + size = ptsize * mpoint->ngeoms; + newpoints = lwalloc(size); + memset(newpoints, 0, size); + + ptr = newpoints; + for(i = 0; i < mpoint->ngeoms; i++) + { + memcpy(ptr, + getPoint_internal(mpoint->geoms[i]->point, 0), + ptsize); + ptr += ptsize; + } + + pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, + mpoint->ngeoms); + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_from_lwmpoint: constructed pointarray for %d points, %d zmflag", mpoint->ngeoms, zmflag); +#endif + + return lwcurve_construct(SRID, NULL, pa); +} + +LWCURVE * +lwcurve_addpoint(LWCURVE *curve, LWPOINT *point, unsigned int where) +{ + POINTARRAY *newpa; + LWCURVE *ret; + + newpa = ptarray_addPoint(curve->points, + getPoint_internal(point->point, 0), + TYPE_NDIMS(point->type), where); + ret = lwcurve_construct(curve->SRID, NULL, newpa); + + return ret; +} + +LWCURVE * +lwcurve_removepoint(LWCURVE *curve, unsigned int index) +{ + POINTARRAY *newpa; + LWCURVE *ret; + + newpa = ptarray_removePoint(curve->points, index); + ret = lwcurve_construct(curve->SRID, NULL, newpa); + + return ret; +} + +/* + * Note: input will be changed, make sure you have permissions for this. + * */ +void +lwcurve_setPoint4d(LWCURVE *curve, unsigned int index, POINT4D *newpoint) +{ + setPoint4d(curve->points, index, newpoint); +} + + + + + + + + + + + + diff --git a/lwgeom/lwgeom_box.c b/lwgeom/lwgeom_box.c index c3b3bdce4..2babe68be 100644 --- a/lwgeom/lwgeom_box.c +++ b/lwgeom/lwgeom_box.c @@ -5,6 +5,9 @@ #include "liblwgeom.h" +void box_to_box2df(BOX *box, BOX2DFLOAT4 *out); + + /* convert postgresql BOX to BOX2D */ void diff --git a/lwgeom/lwgeom_chip.c b/lwgeom/lwgeom_chip.c index 7fe271b79..22c29c6da 100644 --- a/lwgeom/lwgeom_chip.c +++ b/lwgeom/lwgeom_chip.c @@ -18,10 +18,47 @@ /* Define this to debug CHIP ops */ /*#define DEBUG_CHIP 1*/ +typedef unsigned short int UINT16; +typedef float FLOAT32; + +typedef struct PIXEL_T { + int type; /* 1=float32, 5=int24, 6=int16 */ + uchar val[4]; +} PIXEL; + +typedef struct RGB_T { + uchar red; + uchar green; + uchar blue; +} RGB; + + /* Internal funcs */ void swap_char(char *a,char *b); void flip_endian_double(char *d); void flip_endian_int32(char *i); +const char* pixelOpName(int op); +const char* pixelHEX(PIXEL* p); +UINT16 pixel_readUINT16(PIXEL *p); +void pixel_writeUINT16(PIXEL *p, UINT16 i); +PIXEL pixel_readval(char *buf); +void pixel_writeval(PIXEL *p, char *buf, size_t maxlen); +void pixel_add_float32(PIXEL *where, PIXEL *what); +void pixel_add_int24(PIXEL *where, PIXEL *what); +void pixel_add_int16(PIXEL *where, PIXEL *what); +void pixel_add(PIXEL *where, PIXEL *what); +size_t chip_xy_off(CHIP *c, size_t x, size_t y); +void chip_setPixel(CHIP *c, int x, int y, PIXEL *p); +PIXEL chip_getPixel(CHIP *c, int x, int y); +void chip_draw_pixel(CHIP *chip, int x, int y, PIXEL *pixel, int op); +void chip_draw_segment(CHIP *chip, int x1, int y1, int x2, int y2, PIXEL *pixel, int op); +void chip_fill(CHIP *chip, PIXEL *pixel, int op); +CHIP * pgchip_construct(BOX3D *bvol, int SRID, int width, int height, int datatype, PIXEL *initvalue); +void chip_draw_ptarray(CHIP *chip, POINTARRAY *pa, PIXEL *pixel, int op); +void chip_draw_lwpoint(CHIP *chip, LWPOINT *lwpoint, PIXEL* pixel, int op); +void chip_draw_lwline(CHIP *chip, LWLINE *lwline, PIXEL* pixel, int op); +void chip_draw_lwgeom(CHIP *chip, LWGEOM *lwgeom, PIXEL *pixel, int op); +char * text_to_cstring(text *t); /* Prototypes */ @@ -36,6 +73,13 @@ Datum CHIP_getCompression(PG_FUNCTION_ARGS); Datum CHIP_getHeight(PG_FUNCTION_ARGS); Datum CHIP_getWidth(PG_FUNCTION_ARGS); Datum CHIP_setSRID(PG_FUNCTION_ARGS); +Datum CHIP_send(PG_FUNCTION_ARGS); +Datum CHIP_dump(PG_FUNCTION_ARGS); +Datum CHIP_construct(PG_FUNCTION_ARGS); +Datum CHIP_getpixel(PG_FUNCTION_ARGS); +Datum CHIP_setpixel(PG_FUNCTION_ARGS); +Datum CHIP_draw(PG_FUNCTION_ARGS); +Datum CHIP_fill(PG_FUNCTION_ARGS); /* @@ -355,21 +399,6 @@ pixelOpName(int op) return pixelop_name[op]; } - -typedef struct RGB_T { - uchar red; - uchar green; - uchar blue; -} RGB; - -typedef unsigned short int UINT16; -typedef float FLOAT32; - -typedef struct PIXEL_T { - int type; /* 1=float32, 5=int24, 6=int16 */ - uchar val[4]; -} PIXEL; - const char* pixelHEX(PIXEL* p) { diff --git a/lwgeom/lwgeom_functions_analytic.c b/lwgeom/lwgeom_functions_analytic.c index 43d573947..94bda2dbd 100644 --- a/lwgeom/lwgeom_functions_analytic.c +++ b/lwgeom/lwgeom_functions_analytic.c @@ -43,6 +43,15 @@ LWCOLLECTION *simplify2d_collection(const LWCOLLECTION *igeom, double dist); LWGEOM *simplify2d_lwgeom(const LWGEOM *igeom, double dist); Datum LWGEOM_simplify2d(PG_FUNCTION_ARGS); +double determineSide(POINT2D *seg1, POINT2D *seg2, POINT2D *point); +int isOnSegment(POINT2D *seg1, POINT2D *seg2, POINT2D *point); +int point_in_ring(RTREE_NODE *root, POINT2D *point); +int point_in_ring_deprecated(POINTARRAY *pts, POINT2D *point); +int point_in_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point); +int point_in_polygon_deprecated(LWPOLY *polygon, LWPOINT *point); +int point_outside_polygon(RTREE_NODE **root, int ringCount, LWPOINT *point); +int point_outside_polygon_deprecated(LWPOLY *polygon, LWPOINT *point); + /* * Search farthest point from segment p1-p2 diff --git a/lwgeom/lwgeom_functions_basic.c b/lwgeom/lwgeom_functions_basic.c index a356a37b4..228e4e105 100644 --- a/lwgeom/lwgeom_functions_basic.c +++ b/lwgeom/lwgeom_functions_basic.c @@ -73,8 +73,15 @@ Datum LWGEOM_asEWKT(PG_FUNCTION_ARGS); Datum LWGEOM_hasBBOX(PG_FUNCTION_ARGS); Datum LWGEOM_azimuth(PG_FUNCTION_ARGS); Datum LWGEOM_affine(PG_FUNCTION_ARGS); - +Datum LWGEOM_longitude_shift(PG_FUNCTION_ARGS); Datum optimistic_overlap(PG_FUNCTION_ARGS); + +void lwgeom_affine_ptarray(POINTARRAY *pa, double afac, double bfac, double cfac, + double dfac, double efac, double ffac, double gfac, double hfac, double ifac, double xoff, double yoff, double zoff); + +void lwgeom_affine_recursive(uchar *serialized, double afac, double bfac, double cfac, + double dfac, double efac, double ffac, double gfac, double hfac, double ifac, double xoff, double yoff, double zoff); + /*------------------------------------------------------------------*/ /*find the size of geometry */ diff --git a/lwgeom/lwgeom_inout.c b/lwgeom/lwgeom_inout.c index 36a613d34..a2e0aecaf 100644 --- a/lwgeom/lwgeom_inout.c +++ b/lwgeom/lwgeom_inout.c @@ -38,6 +38,7 @@ Datum LWGEOM_out(PG_FUNCTION_ARGS); Datum LWGEOM_to_text(PG_FUNCTION_ARGS); Datum LWGEOM_to_bytea(PG_FUNCTION_ARGS); Datum LWGEOM_from_bytea(PG_FUNCTION_ARGS); +Datum LWGEOM_asHEXEWKB(PG_FUNCTION_ARGS); Datum parse_WKT_lwgeom(PG_FUNCTION_ARGS); #if POSTGIS_PGSQL_VERSION > 73 Datum LWGEOM_recv(PG_FUNCTION_ARGS); diff --git a/lwgeom/lwgeom_ogc.c b/lwgeom/lwgeom_ogc.c index 87873d33e..2b181b4c9 100644 --- a/lwgeom/lwgeom_ogc.c +++ b/lwgeom/lwgeom_ogc.c @@ -81,6 +81,8 @@ Datum LWGEOM_isclosed_linestring(PG_FUNCTION_ARGS); static int32 lwgeom_numpoints_linestring_recursive(const uchar *serialized); static int32 lwgeom_dimension_recursive(const uchar *serialized); char line_is_closed(LWLINE *line); +char curve_is_closed(LWCURVE *curve); +char compound_is_closed(LWCOMPOUND *compound); /*------------------------------------------------------------------*/ diff --git a/lwgeom/lwgeom_rtree.c b/lwgeom/lwgeom_rtree.c index 2bf3b0e8c..52e038d07 100644 --- a/lwgeom/lwgeom_rtree.c +++ b/lwgeom/lwgeom_rtree.c @@ -1,475 +1,480 @@ -/********************************************************************** - * $Id$ - * - * PostGIS - Spatial Types for PostgreSQL - * http://postgis.refractions.net - * Copyright 2001-2005 Refractions Research Inc. - * - * This is free software; you can redistribute and/or modify it under - * the terms of the GNU General Public Licence. See the COPYING file. - * - **********************************************************************/ - -#include "lwgeom_pg.h" -#include "liblwgeom.h" -#include "lwgeom_rtree.h" - -/* - * Creates an rtree given a pointer to the point array. - * Must copy the point array. - */ -RTREE_NODE *createTree(POINTARRAY *pointArray) -{ - RTREE_NODE *root; - RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints); - int i, nodeCount; - int childNodes, parentNodes; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createTree called with pointarray %p", pointArray); -#endif - - nodeCount = pointArray->npoints - 1; - -#ifdef PGIS_DEBUG - lwnotice("Total leaf nodes: %d", nodeCount); -#endif - - /* - * Create a leaf node for every line segment. - */ - for(i = 0; i < nodeCount; i++) - { - nodes[i] = createLeafNode(pointArray, i); - } - - /* - * Next we group nodes by pairs. If there's an odd number of nodes, - * we bring the last node up a level as is. Continue until we have - * a single top node. - */ - childNodes = nodeCount; - parentNodes = nodeCount / 2; - while(parentNodes > 0) - { -#ifdef PGIS_DEBUG - lwnotice("Merging %d children into %d parents.", childNodes, parentNodes); -#endif - i = 0; - while(i < parentNodes) - { - nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]); - i++; - } - /* - * Check for an odd numbered final node. - */ - if(parentNodes * 2 < childNodes) - { -#ifdef PGIS_DEBUG - lwnotice("Shuffling child %d to parent %d", childNodes - 1, i); -#endif - nodes[i] = nodes[childNodes - 1]; - parentNodes++; - } - childNodes = parentNodes; - parentNodes = parentNodes / 2; - } - - root = nodes[0]; - -#ifdef PGIS_DEBUG - lwnotice("createTree returning %p", root); -#endif - return root; -} - -/* - * Creates an interior node given the children. - */ -RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){ - RTREE_NODE *parent; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createInteriorNode called for children %p, %p", left, right); -#endif - parent = lwalloc(sizeof(RTREE_NODE)); - parent->leftNode = left; - parent->rightNode = right; - parent->interval = mergeIntervals(left->interval, right->interval); - parent->segment = NULL; -#ifdef PGIS_DEBUG_CALLS - lwnotice("createInteriorNode returning %p", parent); -#endif - return parent; -} - -/* - * Creates a leaf node given the pointer to the start point of the segment. - */ -RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint) -{ - RTREE_NODE *parent; - LWLINE *line; - double value1; - double value2; - POINT4D tmp; - POINTARRAY *npa; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createLeafNode called for point %d of %p", startPoint, pa); -#endif - if(pa->npoints < startPoint + 2) - { - lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint); - } - - /* - * The given point array will be part of a geometry that will be freed - * independently of the index. Since we may want to cache the index, - * we must create independent arrays. - */ - npa = lwalloc(sizeof(POINTARRAY)); - npa->dims = 0; - npa->npoints = 2; - TYPE_SETZM(npa->dims, 0, 0); - npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2); - - getPoint4d_p(pa, startPoint, &tmp); - - setPoint4d(npa, 0, &tmp); - value1 = tmp.y; - - getPoint4d_p(pa, startPoint + 1, &tmp); - - setPoint4d(npa, 1, &tmp); - value2 = tmp.y; - - line = lwline_construct(-1, NULL, npa); - - parent = lwalloc(sizeof(RTREE_NODE)); - parent->interval = createInterval(value1, value2); - parent->segment = line; - parent->leftNode = NULL; - parent->rightNode = NULL; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createLeafNode returning %p", parent); -#endif - return parent; -} - -/* - * Creates an interval with the total extents of the two given intervals. - */ -INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2) -{ - INTERVAL *interval; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("mergeIntervals called with %p, %p", inter1, inter2); -#endif - interval = lwalloc(sizeof(INTERVAL)); - interval->max = FP_MAX(inter1->max, inter2->max); - interval->min = FP_MIN(inter1->min, inter2->min); -#ifdef PGIS_DEBUG - lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max); -#endif - return interval; -} - -/* - * Creates an interval given the min and max values, in arbitrary order. - */ -INTERVAL *createInterval(double value1, double value2) -{ - INTERVAL *interval; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createInterval called with %8.3f, %8.3f", value1, value2); -#endif - interval = lwalloc(sizeof(INTERVAL)); - interval->max = FP_MAX(value1, value2); - interval->min = FP_MIN(value1, value2); -#ifdef PGIS_DEBUG - lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max); -#endif - return interval; -} - -/* - * Recursively frees the child nodes, the interval and the line before - * freeing the root node. - */ -void freeTree(RTREE_NODE *root) -{ - -#ifdef PGIS_DEBUG_CALLS - lwnotice("freeTree called for %p", root); -#endif - if(root->leftNode) - freeTree(root->leftNode); - if(root->rightNode) - freeTree(root->rightNode); - lwfree(root->interval); - if(root->segment) - lwfree(root->segment); - lwfree(root); -} - -/* - * Retrieves a collection of line segments given the root and crossing value. - * The collection is a multilinestring consisting of two point lines - * representing the segments of the ring that may be crossed by the - * horizontal projection line at the given y value. - */ -LWMLINE *findLineSegments(RTREE_NODE *root, double value) -{ - LWMLINE *tmp, *result; - LWGEOM **lwgeoms; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("findLineSegments called for tree %p and value %8.3f", root, value); -#endif - result = NULL; - - if(!isContained(root->interval, value)) - { -#ifdef PGIS_DEBUG - lwnotice("findLineSegments %p: not contained.", root); -#endif - return NULL; - } - - /* If there is a segment defined for this node, include it. */ - if(root->segment) - { -#ifdef PGIS_DEBUG - lwnotice("findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type)); -#endif - lwgeoms = lwalloc(sizeof(LWGEOM *)); - lwgeoms[0] = (LWGEOM *)root->segment; - -#ifdef PGIS_DEBUG - lwnotice("Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type)); -#endif - result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms); - } - - /* If there is a left child node, recursively include its results. */ - if(root->leftNode) - { -#ifdef PGIS_DEBUG - lwnotice("findLineSegments %p: recursing left.", root); -#endif - tmp = findLineSegments(root->leftNode, value); - if(tmp) - { -#ifdef PGIS_DEBUG - lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); -#endif - if(result) - result = mergeMultiLines(result, tmp); - else - result = tmp; - } - } - - /* Same for any right child. */ - if(root->rightNode) - { -#ifdef PGIS_DEBUG - lwnotice("findLineSegments %p: recursing right.", root); -#endif - tmp = findLineSegments(root->rightNode, value); - if(tmp) - { -#ifdef PGIS_DEBUG - lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); -#endif - if(result) - result = mergeMultiLines(result, tmp); - else - result = tmp; - } - } - - return result; -} - -/* Merges two multilinestrings into a single multilinestring. */ -LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2) -{ - LWGEOM **geoms; - LWCOLLECTION *col; - int i, j, ngeoms; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type)); -#endif - ngeoms = line1->ngeoms + line2->ngeoms; - geoms = lwalloc(sizeof(LWGEOM *) * ngeoms); - - j = 0; - for(i = 0; i < line1->ngeoms; i++, j++) - { - geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]); - } - for(i = 0; i < line2->ngeoms; i++, j++) - { - geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]); - } - col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms); - -#ifdef PGIS_DEBUG_CALLS - lwnotice("mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type)); -#endif - - return (LWMLINE *)col; -} - -/* - * Returns 1 if min < value <= max, 0 otherwise. */ -uint32 isContained(INTERVAL *interval, double value) -{ - return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0; -} - -PG_FUNCTION_INFO_V1(LWGEOM_polygon_index); -Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS) -{ - PG_LWGEOM *igeom, *result; - LWGEOM *geom; - LWPOLY *poly; - LWMLINE *mline; - RTREE_NODE *root; - double yval; - -#ifdef PGIS_DEBUG_CALLS - int i; - lwnotice("polygon_index called."); -#endif - result = NULL; - igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - yval = PG_GETARG_FLOAT8(1); - geom = lwgeom_deserialize(SERIALIZED_FORM(igeom)); - if(TYPE_GETTYPE(geom->type) != POLYGONTYPE) - { - lwgeom_release(geom); - PG_FREE_IF_COPY(igeom, 0); - PG_RETURN_NULL(); - } - poly = (LWPOLY *)geom; - root = createTree(poly->rings[0]); - - mline = findLineSegments(root, yval); - -#ifdef PGIS_DEBUG - lwnotice("mline returned %p %d", mline, TYPE_GETTYPE(mline->type)); - for(i = 0; i < mline->ngeoms; i++) - { - lwnotice("geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type)); - } -#endif - - if(mline) - result = pglwgeom_serialize((LWGEOM *)mline); - -#ifdef PGIS_DEBUG - lwnotice("returning result %p", result); -#endif - lwfree(root); - - PG_FREE_IF_COPY(igeom, 0); - lwgeom_release((LWGEOM *)poly); - lwgeom_release((LWGEOM *)mline); - PG_RETURN_POINTER(result); - -} - -RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly) -{ - RTREE_POLY_CACHE *result; - int i, length; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("createNewCache called with %p", poly); -#endif - result = lwalloc(sizeof(RTREE_POLY_CACHE)); - result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings); - result->ringCount = poly->nrings; - length = lwgeom_size_poly(serializedPoly); - result->poly = lwalloc(length); - memcpy(result->poly, serializedPoly, length); - for(i = 0; i < result->ringCount; i++) - { - result->ringIndices[i] = createTree(poly->rings[i]); - } -#ifdef PGIS_DEBUG - lwnotice("createNewCache returning %p", result); -#endif - return result; -} - -/* - * Creates a new cachable index if needed, or returns the current cache if - * it is applicable to the current polygon. - * The memory context must be changed to function scope before calling this - * method. The method will allocate memory for the cache it creates, - * as well as freeing the memory of any cache that is no longer applicable. - */ -RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, - RTREE_POLY_CACHE *currentCache) -{ - int i, length; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache); -#endif - if(!currentCache) - { -#ifdef PGIS_DEBUG - lwnotice("No existing cache, create one."); -#endif - return createNewCache(poly, serializedPoly); - } - if(!(currentCache->poly)) - { -#ifdef PGIS_DEBUG - lwnotice("Cache contains no polygon, creating new cache."); -#endif - return createNewCache(poly, serializedPoly); - } - - length = lwgeom_size_poly(serializedPoly); - - if(lwgeom_size_poly(currentCache->poly) != length) - { -#ifdef PGIS_DEBUG - lwnotice("Polygon size mismatch, creating new cache."); -#endif - lwfree(currentCache->poly); - lwfree(currentCache); - return createNewCache(poly, serializedPoly); - } - for(i = 0; i < length; i++) - { - uchar a = serializedPoly[i]; - uchar b = currentCache->poly[i]; - if(a != b) - { -#ifdef PGIS_DEBUG - lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b); -#endif - lwfree(currentCache->poly); - lwfree(currentCache); - return createNewCache(poly, serializedPoly); - } - } - -#ifdef PGIS_DEBUG - lwnotice("Polygon match, retaining current cache, %p.", currentCache); -#endif - return currentCache; -} - +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2001-2005 Refractions Research Inc. + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include "lwgeom_pg.h" +#include "liblwgeom.h" +#include "lwgeom_rtree.h" + + +Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS); +RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly); + + +/* + * Creates an rtree given a pointer to the point array. + * Must copy the point array. + */ +RTREE_NODE *createTree(POINTARRAY *pointArray) +{ + RTREE_NODE *root; + RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints); + int i, nodeCount; + int childNodes, parentNodes; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createTree called with pointarray %p", pointArray); +#endif + + nodeCount = pointArray->npoints - 1; + +#ifdef PGIS_DEBUG + lwnotice("Total leaf nodes: %d", nodeCount); +#endif + + /* + * Create a leaf node for every line segment. + */ + for(i = 0; i < nodeCount; i++) + { + nodes[i] = createLeafNode(pointArray, i); + } + + /* + * Next we group nodes by pairs. If there's an odd number of nodes, + * we bring the last node up a level as is. Continue until we have + * a single top node. + */ + childNodes = nodeCount; + parentNodes = nodeCount / 2; + while(parentNodes > 0) + { +#ifdef PGIS_DEBUG + lwnotice("Merging %d children into %d parents.", childNodes, parentNodes); +#endif + i = 0; + while(i < parentNodes) + { + nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]); + i++; + } + /* + * Check for an odd numbered final node. + */ + if(parentNodes * 2 < childNodes) + { +#ifdef PGIS_DEBUG + lwnotice("Shuffling child %d to parent %d", childNodes - 1, i); +#endif + nodes[i] = nodes[childNodes - 1]; + parentNodes++; + } + childNodes = parentNodes; + parentNodes = parentNodes / 2; + } + + root = nodes[0]; + +#ifdef PGIS_DEBUG + lwnotice("createTree returning %p", root); +#endif + return root; +} + +/* + * Creates an interior node given the children. + */ +RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){ + RTREE_NODE *parent; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createInteriorNode called for children %p, %p", left, right); +#endif + parent = lwalloc(sizeof(RTREE_NODE)); + parent->leftNode = left; + parent->rightNode = right; + parent->interval = mergeIntervals(left->interval, right->interval); + parent->segment = NULL; +#ifdef PGIS_DEBUG_CALLS + lwnotice("createInteriorNode returning %p", parent); +#endif + return parent; +} + +/* + * Creates a leaf node given the pointer to the start point of the segment. + */ +RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint) +{ + RTREE_NODE *parent; + LWLINE *line; + double value1; + double value2; + POINT4D tmp; + POINTARRAY *npa; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createLeafNode called for point %d of %p", startPoint, pa); +#endif + if(pa->npoints < startPoint + 2) + { + lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint); + } + + /* + * The given point array will be part of a geometry that will be freed + * independently of the index. Since we may want to cache the index, + * we must create independent arrays. + */ + npa = lwalloc(sizeof(POINTARRAY)); + npa->dims = 0; + npa->npoints = 2; + TYPE_SETZM(npa->dims, 0, 0); + npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2); + + getPoint4d_p(pa, startPoint, &tmp); + + setPoint4d(npa, 0, &tmp); + value1 = tmp.y; + + getPoint4d_p(pa, startPoint + 1, &tmp); + + setPoint4d(npa, 1, &tmp); + value2 = tmp.y; + + line = lwline_construct(-1, NULL, npa); + + parent = lwalloc(sizeof(RTREE_NODE)); + parent->interval = createInterval(value1, value2); + parent->segment = line; + parent->leftNode = NULL; + parent->rightNode = NULL; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createLeafNode returning %p", parent); +#endif + return parent; +} + +/* + * Creates an interval with the total extents of the two given intervals. + */ +INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2) +{ + INTERVAL *interval; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("mergeIntervals called with %p, %p", inter1, inter2); +#endif + interval = lwalloc(sizeof(INTERVAL)); + interval->max = FP_MAX(inter1->max, inter2->max); + interval->min = FP_MIN(inter1->min, inter2->min); +#ifdef PGIS_DEBUG + lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max); +#endif + return interval; +} + +/* + * Creates an interval given the min and max values, in arbitrary order. + */ +INTERVAL *createInterval(double value1, double value2) +{ + INTERVAL *interval; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createInterval called with %8.3f, %8.3f", value1, value2); +#endif + interval = lwalloc(sizeof(INTERVAL)); + interval->max = FP_MAX(value1, value2); + interval->min = FP_MIN(value1, value2); +#ifdef PGIS_DEBUG + lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max); +#endif + return interval; +} + +/* + * Recursively frees the child nodes, the interval and the line before + * freeing the root node. + */ +void freeTree(RTREE_NODE *root) +{ + +#ifdef PGIS_DEBUG_CALLS + lwnotice("freeTree called for %p", root); +#endif + if(root->leftNode) + freeTree(root->leftNode); + if(root->rightNode) + freeTree(root->rightNode); + lwfree(root->interval); + if(root->segment) + lwfree(root->segment); + lwfree(root); +} + +/* + * Retrieves a collection of line segments given the root and crossing value. + * The collection is a multilinestring consisting of two point lines + * representing the segments of the ring that may be crossed by the + * horizontal projection line at the given y value. + */ +LWMLINE *findLineSegments(RTREE_NODE *root, double value) +{ + LWMLINE *tmp, *result; + LWGEOM **lwgeoms; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("findLineSegments called for tree %p and value %8.3f", root, value); +#endif + result = NULL; + + if(!isContained(root->interval, value)) + { +#ifdef PGIS_DEBUG + lwnotice("findLineSegments %p: not contained.", root); +#endif + return NULL; + } + + /* If there is a segment defined for this node, include it. */ + if(root->segment) + { +#ifdef PGIS_DEBUG + lwnotice("findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type)); +#endif + lwgeoms = lwalloc(sizeof(LWGEOM *)); + lwgeoms[0] = (LWGEOM *)root->segment; + +#ifdef PGIS_DEBUG + lwnotice("Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type)); +#endif + result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms); + } + + /* If there is a left child node, recursively include its results. */ + if(root->leftNode) + { +#ifdef PGIS_DEBUG + lwnotice("findLineSegments %p: recursing left.", root); +#endif + tmp = findLineSegments(root->leftNode, value); + if(tmp) + { +#ifdef PGIS_DEBUG + lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); +#endif + if(result) + result = mergeMultiLines(result, tmp); + else + result = tmp; + } + } + + /* Same for any right child. */ + if(root->rightNode) + { +#ifdef PGIS_DEBUG + lwnotice("findLineSegments %p: recursing right.", root); +#endif + tmp = findLineSegments(root->rightNode, value); + if(tmp) + { +#ifdef PGIS_DEBUG + lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type)); +#endif + if(result) + result = mergeMultiLines(result, tmp); + else + result = tmp; + } + } + + return result; +} + +/* Merges two multilinestrings into a single multilinestring. */ +LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2) +{ + LWGEOM **geoms; + LWCOLLECTION *col; + int i, j, ngeoms; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type)); +#endif + ngeoms = line1->ngeoms + line2->ngeoms; + geoms = lwalloc(sizeof(LWGEOM *) * ngeoms); + + j = 0; + for(i = 0; i < line1->ngeoms; i++, j++) + { + geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]); + } + for(i = 0; i < line2->ngeoms; i++, j++) + { + geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]); + } + col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms); + +#ifdef PGIS_DEBUG_CALLS + lwnotice("mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type)); +#endif + + return (LWMLINE *)col; +} + +/* + * Returns 1 if min < value <= max, 0 otherwise. */ +uint32 isContained(INTERVAL *interval, double value) +{ + return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0; +} + +PG_FUNCTION_INFO_V1(LWGEOM_polygon_index); +Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *igeom, *result; + LWGEOM *geom; + LWPOLY *poly; + LWMLINE *mline; + RTREE_NODE *root; + double yval; + +#ifdef PGIS_DEBUG_CALLS + int i; + lwnotice("polygon_index called."); +#endif + result = NULL; + igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + yval = PG_GETARG_FLOAT8(1); + geom = lwgeom_deserialize(SERIALIZED_FORM(igeom)); + if(TYPE_GETTYPE(geom->type) != POLYGONTYPE) + { + lwgeom_release(geom); + PG_FREE_IF_COPY(igeom, 0); + PG_RETURN_NULL(); + } + poly = (LWPOLY *)geom; + root = createTree(poly->rings[0]); + + mline = findLineSegments(root, yval); + +#ifdef PGIS_DEBUG + lwnotice("mline returned %p %d", mline, TYPE_GETTYPE(mline->type)); + for(i = 0; i < mline->ngeoms; i++) + { + lwnotice("geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type)); + } +#endif + + if(mline) + result = pglwgeom_serialize((LWGEOM *)mline); + +#ifdef PGIS_DEBUG + lwnotice("returning result %p", result); +#endif + lwfree(root); + + PG_FREE_IF_COPY(igeom, 0); + lwgeom_release((LWGEOM *)poly); + lwgeom_release((LWGEOM *)mline); + PG_RETURN_POINTER(result); + +} + +RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly) +{ + RTREE_POLY_CACHE *result; + int i, length; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("createNewCache called with %p", poly); +#endif + result = lwalloc(sizeof(RTREE_POLY_CACHE)); + result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings); + result->ringCount = poly->nrings; + length = lwgeom_size_poly(serializedPoly); + result->poly = lwalloc(length); + memcpy(result->poly, serializedPoly, length); + for(i = 0; i < result->ringCount; i++) + { + result->ringIndices[i] = createTree(poly->rings[i]); + } +#ifdef PGIS_DEBUG + lwnotice("createNewCache returning %p", result); +#endif + return result; +} + +/* + * Creates a new cachable index if needed, or returns the current cache if + * it is applicable to the current polygon. + * The memory context must be changed to function scope before calling this + * method. The method will allocate memory for the cache it creates, + * as well as freeing the memory of any cache that is no longer applicable. + */ +RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, + RTREE_POLY_CACHE *currentCache) +{ + int i, length; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache); +#endif + if(!currentCache) + { +#ifdef PGIS_DEBUG + lwnotice("No existing cache, create one."); +#endif + return createNewCache(poly, serializedPoly); + } + if(!(currentCache->poly)) + { +#ifdef PGIS_DEBUG + lwnotice("Cache contains no polygon, creating new cache."); +#endif + return createNewCache(poly, serializedPoly); + } + + length = lwgeom_size_poly(serializedPoly); + + if(lwgeom_size_poly(currentCache->poly) != length) + { +#ifdef PGIS_DEBUG + lwnotice("Polygon size mismatch, creating new cache."); +#endif + lwfree(currentCache->poly); + lwfree(currentCache); + return createNewCache(poly, serializedPoly); + } + for(i = 0; i < length; i++) + { + uchar a = serializedPoly[i]; + uchar b = currentCache->poly[i]; + if(a != b) + { +#ifdef PGIS_DEBUG + lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b); +#endif + lwfree(currentCache->poly); + lwfree(currentCache); + return createNewCache(poly, serializedPoly); + } + } + +#ifdef PGIS_DEBUG + lwnotice("Polygon match, retaining current cache, %p.", currentCache); +#endif + return currentCache; +} + diff --git a/lwgeom/lwgeom_sqlmm.c b/lwgeom/lwgeom_sqlmm.c index ccdda9d1d..172d852ac 100644 --- a/lwgeom/lwgeom_sqlmm.c +++ b/lwgeom/lwgeom_sqlmm.c @@ -1,1176 +1,1197 @@ -/********************************************************************** - * $Id$ - * - * PostGIS - Spatial Types for PostgreSQL - * http://postgis.refractions.net - * Copyright 2001-2006 Refractions Research Inc. - * - * This is free software; you can redistribute and/or modify it under - * the terms of the GNU General Public Licence. See the COPYING file. - * - **********************************************************************/ - -#include -#include -#include -#include -#include - -#include "postgres.h" -#include "liblwgeom.h" -#include "fmgr.h" -#include "wktparse.h" -#include "lwgeom_pg.h" - -/* - * Tolerance used to determine equality. - */ -#define EPSILON_SQLMM 1e-8 - -/* - * Determines (recursively in the case of collections) whether the geometry - * contains at least on arc geometry or segment. - */ -uint32 -has_arc(LWGEOM *geom) -{ - LWCOLLECTION *col; - int i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("has_arc called."); -#endif - - switch(lwgeom_getType(geom->type)) - { - case POINTTYPE: - case LINETYPE: - case POLYGONTYPE: - case MULTIPOINTTYPE: - case MULTILINETYPE: - case MULTIPOLYGONTYPE: - return 0; - case CURVETYPE: - return 1; - /* It's a collection that MAY contain an arc */ - default: - col = (LWCOLLECTION *)geom; - for(i=0; ingeoms; i++) - { - if(has_arc(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(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result) -{ - POINT4D *c; - double cx, cy, cr; - double temp, bc, cd, det; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcircle_center called (%.16f, %.16f), (%.16f, %.16f), (%.16f, %.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y); -#endif - - /* 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(double angle, double zm1, double a1, double zm2, double a2) -{ - -#ifdef PGIS_DEBUG_CALLS - lwnotice("interpolate_arc called."); -#endif - - double frac = fabs((angle - a1) / (a2 - a1)); - double result = frac * (zm2 - zm1) + zm1; - -#ifdef PGIS_DEBUG - lwnotice("interpolate_arc: angle=%.16f, a1=%.16f, a2=%.16f, z1=%.16f, z2=%.16f, frac=%.16f, result=%.16f", angle, a1, a2, zm1, zm2, frac, result); -#endif - - return result; -} - -/******************************************************************************* - * Begin curve segmentize functions - ******************************************************************************/ - -POINTARRAY * -lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad) -{ - POINTARRAY *result; - POINT4D pbuf; - size_t ptsize = sizeof(POINT4D); - unsigned int ptcount; - uchar *pt; - - POINT4D *center; - double radius = 0.0, - sweep = 0.0, - angle = 0.0, - increment = 0.0; - double a1, a2, a3, i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcircle_segmentize called. "); -#endif - - radius = lwcircle_center(p1, p2, p3, ¢er); -#ifdef PGIS_DEBUG - lwnotice("lwcircle_segmentize, (%.16f, %.16f) radius=%.16f", center->x, center->y, radius); -#endif - if(radius < 0) - { - return NULL; - } - - 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); - -#ifdef PGIS_DEBUG - lwnotice("a1 = %.16f, a2 = %.16f, a3 = %.16f", a1, a2, a3); -#endif - - 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; - -#ifdef PGIS_DEBUG - lwnotice("ptcount: %d, perQuad: %d, sweep: %.16f, increment: %.16f", ptcount, perQuad, sweep, increment); -#endif - - 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(angle, p1->z, a1, p2->z, a2); - pbuf.m = interpolate_arc(angle, p1->m, a1, p2->m, a2); - } - else - { - if(sweep > 0) - { - pbuf.z = interpolate_arc(angle, p1->z, a1-(2*M_PI), p2->z, a2); - pbuf.m = interpolate_arc(angle, p1->m, a1-(2*M_PI), p2->m, a2); - } - else - { - pbuf.z = interpolate_arc(angle, p1->z, a1+(2*M_PI), p2->z, a2); - pbuf.m = interpolate_arc(angle, p1->m, a1+(2*M_PI), p2->m, a2); - } - } - } - else - { - if((sweep > 0 && a2 < a3) || (sweep < 0 && a2 > a3)) - { - pbuf.z = interpolate_arc(angle, p2->z, a2, p3->z, a3); - pbuf.m = interpolate_arc(angle, p2->m, a2, p3->m, a3); - } - else - { - if(sweep > 0) - { - pbuf.z = interpolate_arc(angle, p2->z, a2-(2*M_PI), p3->z, a3); - pbuf.m = interpolate_arc(angle, p2->m, a2-(2*M_PI), p3->m, a3); - } - else - { - pbuf.z = interpolate_arc(angle, p2->z, a2+(2*M_PI), p3->z, a3); - pbuf.m = interpolate_arc(angle, p2->m, a2+(2*M_PI), p3->m, a3); - } - } - } - memcpy(pt, (uchar *)&pbuf, ptsize); - } - - pt = getPoint_internal(result, ptcount - 1); - memcpy(pt, (uchar *)p3, ptsize); - - lwfree(center); - - return result; -} - -LWLINE * -lwcurve_segmentize(LWCURVE *icurve, uint32 perQuad) -{ - LWLINE *oline; - DYNPTARRAY *ptarray; - POINTARRAY *tmp; - uint32 i, j; - POINT4D *p1 = lwalloc(sizeof(POINT4D)); - POINT4D *p2 = lwalloc(sizeof(POINT4D)); - POINT4D *p3 = lwalloc(sizeof(POINT4D)); - POINT4D *p4 = lwalloc(sizeof(POINT4D)); - - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_segmentize called., dim = %d", icurve->points->dims); -#endif - - ptarray = dynptarray_create(icurve->points->npoints, icurve->points->dims); - if(!getPoint4d_p(icurve->points, 0, p4)) - { - elog(ERROR, "curve_segmentize: Cannot extract point."); - } - dynptarray_addPoint4d(ptarray, p4, 1); - - for(i = 2; i < icurve->points->npoints; i+=2) - { - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_segmentize: arc ending at point %d", i); -#endif - - getPoint4d_p(icurve->points, i - 2, p1); - getPoint4d_p(icurve->points, i - 1, p2); - getPoint4d_p(icurve->points, i, p3); - tmp = lwcircle_segmentize(p1, p2, p3, perQuad); - -#ifdef PGIS_DEBUG - lwnotice("lwcurve_segmentize: generated %d points", tmp->npoints); -#endif - - for(j = 0; j < tmp->npoints; j++) - { - getPoint4d_p(tmp, j, p4); - dynptarray_addPoint4d(ptarray, p4, 1); - } - lwfree(tmp); - } - oline = lwline_construct(icurve->SRID, NULL, ptarray_clone(ptarray->pa)); - - lwfree(p1); - lwfree(p2); - lwfree(p3); - lwfree(p4); - lwfree(ptarray); - return oline; -} - -LWLINE * -lwcompound_segmentize(LWCOMPOUND *icompound, uint32 perQuad) -{ - LWLINE *oline; - LWGEOM *geom; - DYNPTARRAY *ptarray = NULL; - LWLINE *tmp = NULL; - uint32 i, j; - POINT4D *p = NULL; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcompound_segmentize called."); -#endif - p = lwalloc(sizeof(POINT4D)); - - ptarray = dynptarray_create(2, ((POINTARRAY *)icompound->geoms[0]->data)->dims); - - for(i = 0; i < icompound->ngeoms; i++) - { - geom = icompound->geoms[i]; - if(lwgeom_getType(geom->type) == CURVETYPE) - { - tmp = lwcurve_segmentize((LWCURVE *)geom, perQuad); - for(j = 0; j < tmp->points->npoints; j++) - { - getPoint4d_p(tmp->points, j, p); - dynptarray_addPoint4d(ptarray, p, 0); - } - lwfree(tmp); - } - else if(lwgeom_getType(geom->type) == LINETYPE) - { - tmp = (LWLINE *)geom; - for(j = 0; j < tmp->points->npoints; j++) - { - getPoint4d_p(tmp->points, j, p); - dynptarray_addPoint4d(ptarray, p, 0); - } - } - else - { - lwerror("Unsupported geometry type %d found.", lwgeom_getType(geom->type)); - return NULL; - } - } - oline = lwline_construct(icompound->SRID, NULL, ptarray_clone(ptarray->pa)); - lwfree(ptarray); - lwfree(p); - return oline; -} - -LWPOLY * -lwcurvepoly_segmentize(LWCURVEPOLY *curvepoly, uint32 perQuad) -{ - LWPOLY *ogeom; - LWGEOM *tmp; - LWLINE *line; - POINTARRAY **ptarray; - int i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcurvepoly_segmentize called."); -#endif - - ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings); - - for(i = 0; i < curvepoly->nrings; i++) - { - tmp = curvepoly->rings[i]; - if(lwgeom_getType(tmp->type) == CURVETYPE) - { - line = lwcurve_segmentize((LWCURVE *)tmp, perQuad); - ptarray[i] = ptarray_clone(line->points); - lwfree(line); - } - else if(lwgeom_getType(tmp->type) == LINETYPE) - { - line = (LWLINE *)tmp; - ptarray[i] = ptarray_clone(line->points); - } - else - { - lwerror("Invalid ring type found in CurvePoly."); - return NULL; - } - } - - ogeom = lwpoly_construct(curvepoly->SRID, NULL, curvepoly->nrings, ptarray); - return ogeom; -} - -LWMLINE * -lwmcurve_segmentize(LWMCURVE *mcurve, uint32 perQuad) -{ - LWMLINE *ogeom; - LWGEOM *tmp; - LWGEOM **lines; - int i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwmcurve_segmentize called, geoms=%d, dim=%d.", mcurve->ngeoms, TYPE_NDIMS(mcurve->type)); -#endif - - lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms); - - for(i = 0; i < mcurve->ngeoms; i++) - { - tmp = mcurve->geoms[i]; - if(lwgeom_getType(tmp->type) == CURVETYPE) - { - lines[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad); - } - else if(lwgeom_getType(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(LWMSURFACE *msurface, uint32 perQuad) -{ - LWMPOLY *ogeom; - LWGEOM *tmp; - LWPOLY *poly; - LWGEOM **polys; - POINTARRAY **ptarray; - int i, j; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwmsurface_segmentize called."); -#endif - - polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms); - - for(i = 0; i < msurface->ngeoms; i++) - { - tmp = msurface->geoms[i]; - if(lwgeom_getType(tmp->type) == CURVEPOLYTYPE) - { - polys[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad); - } - else if(lwgeom_getType(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(LWCOLLECTION *collection, uint32 perQuad) -{ - LWCOLLECTION *ocol; - LWGEOM *tmp; - LWGEOM **geoms; - int i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwcollection_segmentize called."); -#endif - - if(has_arc((LWGEOM *)collection) == 0) - { - return collection; - } - - geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms); - - for(i=0; ingeoms; i++) - { - tmp = collection->geoms[i]; - switch(lwgeom_getType(tmp->type)) { - case CURVETYPE: - geoms[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad); - break; - case COMPOUNDTYPE: - geoms[i] = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad); - break; - case CURVEPOLYTYPE: - geoms[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)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(LWGEOM *geom, uint32 perQuad) -{ - LWGEOM * ogeom = NULL; - switch(lwgeom_getType(geom->type)) { - case CURVETYPE: - ogeom = (LWGEOM *)lwcurve_segmentize((LWCURVE *)geom, perQuad); - break; - case COMPOUNDTYPE: - ogeom = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)geom, perQuad); - break; - case CURVEPOLYTYPE: - ogeom = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)geom, perQuad); - break; - case MULTICURVETYPE: - ogeom = (LWGEOM *)lwmcurve_segmentize((LWMCURVE *)geom, perQuad); - break; - case MULTISURFACETYPE: - ogeom = (LWGEOM *)lwmsurface_segmentize((LWMSURFACE *)geom, perQuad); - break; - case COLLECTIONTYPE: - ogeom = (LWGEOM *)lwcollection_segmentize((LWCOLLECTION *)geom, perQuad); - break; - default: - ogeom = lwgeom_clone(geom); - } - return ogeom; -} - -/******************************************************************************* - * End curve segmentize functions - ******************************************************************************/ -LWGEOM * -append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int SRID) -{ - LWGEOM *result; - int currentType, i; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("append_segment called %p, %p, %d, %d", geom, pts, type, SRID); -#endif - - if(geom == NULL) - { - if(type == LINETYPE) - { -#ifdef PGIS_DEBUG - lwnotice("append_segment: line to NULL"); -#endif - return (LWGEOM *)lwline_construct(SRID, NULL, pts); - } - else if(type == CURVETYPE) - { -#ifdef PGIS_DEBUG - lwnotice("append_segment: curve to NULL %d", pts->npoints); - POINT4D tmp; - for(i=0; inpoints; i++) - { - getPoint4d_p(pts, i, &tmp); - lwnotice("new point: (%.16f,%.16f)",tmp.x,tmp.y); - } -#endif - return (LWGEOM *)lwcurve_construct(SRID, NULL, pts); - } - else - { - lwerror("Invalid segment type %d.", type); - } - } - - currentType = lwgeom_getType(geom->type); - if(currentType == LINETYPE && type == LINETYPE) - { - POINTARRAY *newPoints; - POINT4D pt; - LWLINE *line = (LWLINE *)geom; -#ifdef PGIS_DEBUG - lwnotice("append_segment: line to line"); -#endif - newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + line->points->npoints - 1); - for(i=0; ipoints->npoints; i++) - { - getPoint4d_p(pts, i, &pt); - setPoint4d(newPoints, i, &pt); - } - for(i=1; inpoints; i++) - { - getPoint4d_p(line->points, i, &pt); - setPoint4d(newPoints, i + pts->npoints, &pt); - } - result = (LWGEOM *)lwline_construct(SRID, NULL, newPoints); - lwgeom_release(geom); - return result; - } - else if(currentType == CURVETYPE && type == CURVETYPE) - { - POINTARRAY *newPoints; - POINT4D pt; - LWCURVE *curve = (LWCURVE *)geom; -#ifdef PGIS_DEBUG - lwnotice("append_segment: curve to curve"); -#endif - newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + curve->points->npoints - 1); -#ifdef PGIS_DEBUG - lwnotice("New array length: %d", pts->npoints + curve->points->npoints - 1); -#endif - for(i=0; ipoints->npoints; i++) - { - getPoint4d_p(curve->points, i, &pt); -#ifdef PGIS_DEBUG - lwnotice("orig point %d: (%.16f,%.16f)", i, pt.x, pt.y); -#endif - setPoint4d(newPoints, i, &pt); - } - for(i=1; inpoints;i++) - { - getPoint4d_p(pts, i, &pt); -#ifdef PGIS_DEBUG - lwnotice("new point %d: (%.16f,%.16f)", i + curve->points->npoints - 1, pt.x, pt.y); -#endif - setPoint4d(newPoints, i + curve->points->npoints - 1, &pt); - } - result = (LWGEOM *)lwcurve_construct(SRID, NULL, newPoints); - lwgeom_release(geom); - return result; - } - else if(currentType == CURVETYPE && type == LINETYPE) - { - LWLINE *line; - LWGEOM **geomArray; -#ifdef PGIS_DEBUG - lwnotice("append_segment: line to curve"); -#endif - geomArray = lwalloc(sizeof(LWGEOM *)*2); - geomArray[0] = lwgeom_clone(geom); - - line = lwline_construct(SRID, NULL, pts); - geomArray[1] = lwgeom_clone((LWGEOM *)line); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray); - lwfree((LWGEOM *)line); - lwgeom_release(geom); - return result; - } - else if(currentType == LINETYPE && type == CURVETYPE) - { - LWCURVE *curve; - LWGEOM **geomArray; -#ifdef PGIS_DEBUG - lwnotice("append_segment: curve to line"); -#endif - geomArray = lwalloc(sizeof(LWGEOM *)*2); - geomArray[0] = lwgeom_clone(geom); - - curve = lwcurve_construct(SRID, NULL, pts); - geomArray[1] = lwgeom_clone((LWGEOM *)curve); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray); - lwfree((LWGEOM *)curve); - lwgeom_release(geom); - return result; - } - else if(currentType == COMPOUNDTYPE) - { - LWGEOM *newGeom; - LWCOMPOUND *compound; - int count; - LWGEOM **geomArray; - - compound = (LWCOMPOUND *)geom; - count = compound->ngeoms + 1; - geomArray = lwalloc(sizeof(LWGEOM *)*count); - for(i=0; ingeoms; i++) - { - geomArray[i] = lwgeom_clone(compound->geoms[i]); - } - if(type == LINETYPE) - { -#ifdef PGIS_DEBUG - lwnotice("append_segment: line to compound"); -#endif - newGeom = (LWGEOM *)lwline_construct(SRID, NULL, pts); - } - else if(type == CURVETYPE) - { -#ifdef PGIS_DEBUG - lwnotice("append_segment: curve to compound"); -#endif - newGeom = (LWGEOM *)lwcurve_construct(SRID, NULL, pts); - } - else - { - lwerror("Invalid segment type %d.", type); - return NULL; - } - geomArray[compound->ngeoms] = lwgeom_clone(newGeom); - - result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, count, geomArray); - lwfree(newGeom); - lwgeom_release(geom); - return result; - } - lwerror("Invalid state %d-%d", currentType, type); - return NULL; -} - -LWGEOM * -pta_desegmentize(POINTARRAY *points, int type, int SRID) -{ - int i, j, commit, isline, count; - double last_angle, last_length; - double dxab, dyab, dxbc, dybc, theta, length; - POINT4D a, b, c, tmp; - POINTARRAY *pts; - LWGEOM *geom = NULL; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("pta_desegmentize called."); -#endif - - getPoint4d_p(points, 0, &a); - getPoint4d_p(points, 1, &b); - getPoint4d_p(points, 2, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - - theta = atan2(dyab, dxab); - last_angle = theta - atan2(dybc, dxbc); - last_length = sqrt(dxbc*dxbc+dybc*dybc); - length = sqrt(dxab*dxab+dyab*dyab); - if((last_length - length) < EPSILON_SQLMM) - { - isline = -1; -#ifdef PGIS_DEBUG - lwnotice("Starting as unknown."); -#endif - } - else - { - isline = 1; -#ifdef PGIS_DEBUG - lwnotice("Starting as line."); -#endif - } - - commit = 0; - for(i=3; inpoints; i++) - { - getPoint4d_p(points, i-2, &a); - getPoint4d_p(points, i-1, &b); - getPoint4d_p(points, i, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - -#ifdef PGIS_DEBUG - lwnotice("(dxab, dyab, dxbc, dybc) (%.16f, %.16f, %.16f, %.16f)", dxab, dyab, dxbc, dybc); -#endif - - theta = atan2(dyab, dxab); - theta = theta - atan2(dybc, dxbc); - length = sqrt(dxbc*dxbc+dybc*dybc); -#ifdef PGIS_DEBUG - lwnotice("Last/current length and angle %.16f/%.16f, %.16f/%.16f", last_angle, theta, last_length, length); -#endif - /* Found a line segment */ - if(fabs(length - last_length) > EPSILON_SQLMM || - fabs(theta - last_angle) > EPSILON_SQLMM) - { - last_length = length; - last_angle = theta; - /* We are tracking a line, keep going */ - if(isline > 0) - { - } - /* We were tracking a curve, commit it and start line*/ - else if(isline == 0) - { -#ifdef PGIS_DEBUG - lwnotice("Building curve, %d - %d", commit, i); -#endif - count = i - commit; - pts = ptarray_construct( - TYPE_HASZ(type), - TYPE_HASM(type), - 3); - getPoint4d_p(points, commit, &tmp); - setPoint4d(pts, 0, &tmp); - getPoint4d_p(points, - commit + count/2, &tmp); - setPoint4d(pts, 1, &tmp); - getPoint4d_p(points, i - 1, &tmp); - setPoint4d(pts, 2, &tmp); - - commit = i-1; - geom = append_segment(geom, pts, CURVETYPE, SRID); - isline = -1; - - /* - * We now need to move ahead one point to - * determine if it's a potential new curve, - * since the last_angle value is corrupt. - */ - i++; - getPoint4d_p(points, i-2, &a); - getPoint4d_p(points, i-1, &b); - getPoint4d_p(points, i, &c); - - dxab = b.x - a.x; - dyab = b.y - a.y; - dxbc = c.x - b.x; - dybc = c.y - b.y; - - theta = atan2(dyab, dxab); - last_angle = theta - atan2(dybc, dxbc); - last_length = sqrt(dxbc*dxbc+dybc*dybc); - length = sqrt(dxab*dxab+dyab*dyab); - if((last_length - length) < EPSILON_SQLMM) - { - isline = -1; -#ifdef PGIS_DEBUG - lwnotice("Restarting as unknown."); -#endif - } - else - { - isline = 1; -#ifdef PGIS_DEBUG - lwnotice("Restarting as line."); -#endif - } - - - } - /* We didn't know what we were tracking, now we do. */ - else - { -#ifdef PGIS_DEBUG - lwnotice("It's a line"); -#endif - isline = 1; - } - } - /* Found a curve segment */ - else - { - /* We were tracking a curve, commit it and start line */ - if(isline > 0) - { -#ifdef PGIS_DEBUG - lwnotice("Building line, %d - %d", commit, i-2); -#endif - count = i - commit - 2; - - pts = ptarray_construct( - TYPE_HASZ(type), - TYPE_HASM(type), - count); - for(j=commit;j 2) - { -#ifdef PGIS_DEBUG - lwnotice("Finishing curve %d,%d.", commit, i); -#endif - pts = ptarray_construct( - TYPE_HASZ(type), - TYPE_HASM(type), - 3); - getPoint4d_p(points, commit, &tmp); - setPoint4d(pts, 0, &tmp); - getPoint4d_p(points, commit + count/2, &tmp); - setPoint4d(pts, 1, &tmp); - getPoint4d_p(points, i - 1, &tmp); - setPoint4d(pts, 2, &tmp); - - geom = append_segment(geom, pts, CURVETYPE, SRID); - } - else - { -#ifdef PGIS_DEBUG - lwnotice("Finishing line %d,%d.", commit, i); -#endif - pts = ptarray_construct( - TYPE_HASZ(type), - TYPE_HASM(type), - count); - for(j=commit;jpoints, line->type, line->SRID); -} - -LWGEOM * -lwpolygon_desegmentize(LWPOLY *poly) -{ - LWGEOM **geoms; - int i, hascurve = 0; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwpolygon_desegmentize called."); -#endif - - geoms = lwalloc(sizeof(LWGEOM *)*poly->nrings); - for(i=0; inrings; i++) - { - geoms[i] = pta_desegmentize(poly->rings[i], poly->type, poly->SRID); - if(lwgeom_getType(geoms[i]->type) == CURVETYPE || - lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE) - { - hascurve = 1; - } - } - if(hascurve == 0) - { - for(i=0; inrings; i++) - { - lwfree(geoms[i]); - } - return lwgeom_clone((LWGEOM *)poly); - } - - return (LWGEOM *)lwcollection_construct(CURVEPOLYTYPE, poly->SRID, NULL, poly->nrings, geoms); -} - -LWGEOM * -lwmline_desegmentize(LWMLINE *mline) -{ - LWGEOM **geoms; - int i, hascurve = 0; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwmline_desegmentize called."); -#endif - - geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms); - for(i=0; ingeoms; i++) - { - geoms[i] = lwline_desegmentize((LWLINE *)mline->geoms[i]); - if(lwgeom_getType(geoms[i]->type) == CURVETYPE || - lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE) - { - hascurve = 1; - } - } - if(hascurve == 0) - { - for(i=0; ingeoms; i++) - { - lwfree(geoms[i]); - } - return lwgeom_clone((LWGEOM *)mline); - } - return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->SRID, NULL, mline->ngeoms, geoms); -} - -LWGEOM * -lwmpolygon_desegmentize(LWMPOLY *mpoly) -{ - LWGEOM **geoms; - int i, hascurve = 0; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwmpoly_desegmentize called."); -#endif - - geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms); - for(i=0; ingeoms; i++) - { - geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]); - if(lwgeom_getType(geoms[i]->type) == CURVEPOLYTYPE) - { - hascurve = 1; - } - } - if(hascurve == 0) - { - for(i=0; ingeoms; i++) - { - lwfree(geoms[i]); - } - return lwgeom_clone((LWGEOM *)mpoly); - } - return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->SRID, NULL, mpoly->ngeoms, geoms); -} - -LWGEOM * -lwgeom_desegmentize(LWGEOM *geom) -{ - int type = lwgeom_getType(geom->type); - -#ifdef PGIS_DEBUG_CALLS - lwnotice("lwgeom_desegmentize called."); -#endif - - switch(type) { - case LINETYPE: - return lwline_desegmentize((LWLINE *)geom); - case POLYGONTYPE: - return lwpolygon_desegmentize((LWPOLY *)geom); - case MULTILINETYPE: - return lwmline_desegmentize((LWMLINE *)geom); - case MULTIPOLYGONTYPE: - return lwmpolygon_desegmentize((LWMPOLY *)geom); - default: - return lwgeom_clone(geom); - } -} - -/******************************************************************************* - * Begin PG_FUNCTIONs - ******************************************************************************/ - -PG_FUNCTION_INFO_V1(LWGEOM_has_arc); -Datum LWGEOM_has_arc(PG_FUNCTION_ARGS) -{ - PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - uint32 result = has_arc(lwgeom_deserialize(SERIALIZED_FORM(geom))); - PG_RETURN_BOOL(result == 1); -} - -/* - * Converts any curve segments of the geometry into a linear approximation. - * Curve centers are determined by projecting the defining points into the 2d - * plane. Z and M values are assigned by linear interpolation between - * defining points. - */ -PG_FUNCTION_INFO_V1(LWGEOM_curve_segmentize); -Datum LWGEOM_curve_segmentize(PG_FUNCTION_ARGS) -{ - PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - uint32 perQuad = PG_GETARG_INT32(1); - PG_LWGEOM *ret; - LWGEOM *igeom = NULL, *ogeom = NULL; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("LWGEOM_curve_segmentize called."); -#endif - - if(perQuad < 0) - { - elog(ERROR, "2nd argument must be positive."); - PG_RETURN_NULL(); - } -#ifdef PGIS_DEBUG - else - { - lwnotice("perQuad = %d", perQuad); - } -#endif - igeom = lwgeom_deserialize(SERIALIZED_FORM(geom)); - ogeom = lwgeom_segmentize(igeom, perQuad); - if(ogeom == NULL) PG_RETURN_NULL(); - ret = pglwgeom_serialize(ogeom); - lwgeom_release(igeom); - lwgeom_release(ogeom); - PG_FREE_IF_COPY(geom, 0); - PG_RETURN_POINTER(ret); -} - -PG_FUNCTION_INFO_V1(LWGEOM_line_desegmentize); -Datum LWGEOM_line_desegmentize(PG_FUNCTION_ARGS) -{ - PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); - PG_LWGEOM *ret; - LWGEOM *igeom = NULL, *ogeom = NULL; - -#ifdef PGIS_DEBUG_CALLS - lwnotice("LWGEOM_line_desegmentize."); -#endif - - igeom = lwgeom_deserialize(SERIALIZED_FORM(geom)); - ogeom = lwgeom_desegmentize(igeom); - if(ogeom == NULL) - { - lwgeom_release(igeom); - PG_RETURN_NULL(); - } - ret = pglwgeom_serialize(ogeom); - lwgeom_release(igeom); - lwgeom_release(ogeom); - PG_FREE_IF_COPY(geom, 0); - PG_RETURN_POINTER(ret); -} - -/******************************************************************************* - * End PG_FUNCTIONs - ******************************************************************************/ +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2001-2006 Refractions Research Inc. + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU General Public Licence. See the COPYING file. + * + **********************************************************************/ + +#include +#include +#include +#include +#include + +#include "postgres.h" +#include "liblwgeom.h" +#include "fmgr.h" +#include "wktparse.h" +#include "lwgeom_pg.h" + + +double interpolate_arc(double angle, double zm1, double a1, double zm2, double a2); +POINTARRAY *lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad); +LWLINE *lwcurve_segmentize(LWCURVE *icurve, uint32 perQuad); +LWLINE *lwcompound_segmentize(LWCOMPOUND *icompound, uint32 perQuad); +LWPOLY *lwcurvepoly_segmentize(LWCURVEPOLY *curvepoly, uint32 perQuad); +LWMLINE *lwmcurve_segmentize(LWMCURVE *mcurve, uint32 perQuad); +LWMPOLY *lwmsurface_segmentize(LWMSURFACE *msurface, uint32 perQuad); +LWCOLLECTION *lwcollection_segmentize(LWCOLLECTION *collection, uint32 perQuad); +LWGEOM *append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int SRID); +LWGEOM *pta_desegmentize(POINTARRAY *points, int type, int SRID); +LWGEOM *lwline_desegmentize(LWLINE *line); +LWGEOM *lwpolygon_desegmentize(LWPOLY *poly); +LWGEOM *lwmline_desegmentize(LWMLINE *mline); +LWGEOM *lwmpolygon_desegmentize(LWMPOLY *mpoly); +LWGEOM *lwgeom_desegmentize(LWGEOM *geom); +Datum LWGEOM_has_arc(PG_FUNCTION_ARGS); +Datum LWGEOM_curve_segmentize(PG_FUNCTION_ARGS); +Datum LWGEOM_line_desegmentize(PG_FUNCTION_ARGS); + + +/* + * Tolerance used to determine equality. + */ +#define EPSILON_SQLMM 1e-8 + +/* + * Determines (recursively in the case of collections) whether the geometry + * contains at least on arc geometry or segment. + */ +uint32 +has_arc(LWGEOM *geom) +{ + LWCOLLECTION *col; + int i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("has_arc called."); +#endif + + switch(lwgeom_getType(geom->type)) + { + case POINTTYPE: + case LINETYPE: + case POLYGONTYPE: + case MULTIPOINTTYPE: + case MULTILINETYPE: + case MULTIPOLYGONTYPE: + return 0; + case CURVETYPE: + return 1; + /* It's a collection that MAY contain an arc */ + default: + col = (LWCOLLECTION *)geom; + for(i=0; ingeoms; i++) + { + if(has_arc(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(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result) +{ + POINT4D *c; + double cx, cy, cr; + double temp, bc, cd, det; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcircle_center called (%.16f, %.16f), (%.16f, %.16f), (%.16f, %.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y); +#endif + + /* 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(double angle, double zm1, double a1, double zm2, double a2) +{ + +#ifdef PGIS_DEBUG_CALLS + lwnotice("interpolate_arc called."); +#endif + + double frac = fabs((angle - a1) / (a2 - a1)); + double result = frac * (zm2 - zm1) + zm1; + +#ifdef PGIS_DEBUG + lwnotice("interpolate_arc: angle=%.16f, a1=%.16f, a2=%.16f, z1=%.16f, z2=%.16f, frac=%.16f, result=%.16f", angle, a1, a2, zm1, zm2, frac, result); +#endif + + return result; +} + +/******************************************************************************* + * Begin curve segmentize functions + ******************************************************************************/ + +POINTARRAY * +lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad) +{ + POINTARRAY *result; + POINT4D pbuf; + size_t ptsize = sizeof(POINT4D); + unsigned int ptcount; + uchar *pt; + + POINT4D *center; + double radius = 0.0, + sweep = 0.0, + angle = 0.0, + increment = 0.0; + double a1, a2, a3, i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcircle_segmentize called. "); +#endif + + radius = lwcircle_center(p1, p2, p3, ¢er); +#ifdef PGIS_DEBUG + lwnotice("lwcircle_segmentize, (%.16f, %.16f) radius=%.16f", center->x, center->y, radius); +#endif + if(radius < 0) + { + return NULL; + } + + 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); + +#ifdef PGIS_DEBUG + lwnotice("a1 = %.16f, a2 = %.16f, a3 = %.16f", a1, a2, a3); +#endif + + 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; + +#ifdef PGIS_DEBUG + lwnotice("ptcount: %d, perQuad: %d, sweep: %.16f, increment: %.16f", ptcount, perQuad, sweep, increment); +#endif + + 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(angle, p1->z, a1, p2->z, a2); + pbuf.m = interpolate_arc(angle, p1->m, a1, p2->m, a2); + } + else + { + if(sweep > 0) + { + pbuf.z = interpolate_arc(angle, p1->z, a1-(2*M_PI), p2->z, a2); + pbuf.m = interpolate_arc(angle, p1->m, a1-(2*M_PI), p2->m, a2); + } + else + { + pbuf.z = interpolate_arc(angle, p1->z, a1+(2*M_PI), p2->z, a2); + pbuf.m = interpolate_arc(angle, p1->m, a1+(2*M_PI), p2->m, a2); + } + } + } + else + { + if((sweep > 0 && a2 < a3) || (sweep < 0 && a2 > a3)) + { + pbuf.z = interpolate_arc(angle, p2->z, a2, p3->z, a3); + pbuf.m = interpolate_arc(angle, p2->m, a2, p3->m, a3); + } + else + { + if(sweep > 0) + { + pbuf.z = interpolate_arc(angle, p2->z, a2-(2*M_PI), p3->z, a3); + pbuf.m = interpolate_arc(angle, p2->m, a2-(2*M_PI), p3->m, a3); + } + else + { + pbuf.z = interpolate_arc(angle, p2->z, a2+(2*M_PI), p3->z, a3); + pbuf.m = interpolate_arc(angle, p2->m, a2+(2*M_PI), p3->m, a3); + } + } + } + memcpy(pt, (uchar *)&pbuf, ptsize); + } + + pt = getPoint_internal(result, ptcount - 1); + memcpy(pt, (uchar *)p3, ptsize); + + lwfree(center); + + return result; +} + +LWLINE * +lwcurve_segmentize(LWCURVE *icurve, uint32 perQuad) +{ + LWLINE *oline; + DYNPTARRAY *ptarray; + POINTARRAY *tmp; + uint32 i, j; + POINT4D *p1 = lwalloc(sizeof(POINT4D)); + POINT4D *p2 = lwalloc(sizeof(POINT4D)); + POINT4D *p3 = lwalloc(sizeof(POINT4D)); + POINT4D *p4 = lwalloc(sizeof(POINT4D)); + + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_segmentize called., dim = %d", icurve->points->dims); +#endif + + ptarray = dynptarray_create(icurve->points->npoints, icurve->points->dims); + if(!getPoint4d_p(icurve->points, 0, p4)) + { + elog(ERROR, "curve_segmentize: Cannot extract point."); + } + dynptarray_addPoint4d(ptarray, p4, 1); + + for(i = 2; i < icurve->points->npoints; i+=2) + { + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_segmentize: arc ending at point %d", i); +#endif + + getPoint4d_p(icurve->points, i - 2, p1); + getPoint4d_p(icurve->points, i - 1, p2); + getPoint4d_p(icurve->points, i, p3); + tmp = lwcircle_segmentize(p1, p2, p3, perQuad); + +#ifdef PGIS_DEBUG + lwnotice("lwcurve_segmentize: generated %d points", tmp->npoints); +#endif + + for(j = 0; j < tmp->npoints; j++) + { + getPoint4d_p(tmp, j, p4); + dynptarray_addPoint4d(ptarray, p4, 1); + } + lwfree(tmp); + } + oline = lwline_construct(icurve->SRID, NULL, ptarray_clone(ptarray->pa)); + + lwfree(p1); + lwfree(p2); + lwfree(p3); + lwfree(p4); + lwfree(ptarray); + return oline; +} + +LWLINE * +lwcompound_segmentize(LWCOMPOUND *icompound, uint32 perQuad) +{ + LWLINE *oline; + LWGEOM *geom; + DYNPTARRAY *ptarray = NULL; + LWLINE *tmp = NULL; + uint32 i, j; + POINT4D *p = NULL; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcompound_segmentize called."); +#endif + p = lwalloc(sizeof(POINT4D)); + + ptarray = dynptarray_create(2, ((POINTARRAY *)icompound->geoms[0]->data)->dims); + + for(i = 0; i < icompound->ngeoms; i++) + { + geom = icompound->geoms[i]; + if(lwgeom_getType(geom->type) == CURVETYPE) + { + tmp = lwcurve_segmentize((LWCURVE *)geom, perQuad); + for(j = 0; j < tmp->points->npoints; j++) + { + getPoint4d_p(tmp->points, j, p); + dynptarray_addPoint4d(ptarray, p, 0); + } + lwfree(tmp); + } + else if(lwgeom_getType(geom->type) == LINETYPE) + { + tmp = (LWLINE *)geom; + for(j = 0; j < tmp->points->npoints; j++) + { + getPoint4d_p(tmp->points, j, p); + dynptarray_addPoint4d(ptarray, p, 0); + } + } + else + { + lwerror("Unsupported geometry type %d found.", lwgeom_getType(geom->type)); + return NULL; + } + } + oline = lwline_construct(icompound->SRID, NULL, ptarray_clone(ptarray->pa)); + lwfree(ptarray); + lwfree(p); + return oline; +} + +LWPOLY * +lwcurvepoly_segmentize(LWCURVEPOLY *curvepoly, uint32 perQuad) +{ + LWPOLY *ogeom; + LWGEOM *tmp; + LWLINE *line; + POINTARRAY **ptarray; + int i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcurvepoly_segmentize called."); +#endif + + ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings); + + for(i = 0; i < curvepoly->nrings; i++) + { + tmp = curvepoly->rings[i]; + if(lwgeom_getType(tmp->type) == CURVETYPE) + { + line = lwcurve_segmentize((LWCURVE *)tmp, perQuad); + ptarray[i] = ptarray_clone(line->points); + lwfree(line); + } + else if(lwgeom_getType(tmp->type) == LINETYPE) + { + line = (LWLINE *)tmp; + ptarray[i] = ptarray_clone(line->points); + } + else + { + lwerror("Invalid ring type found in CurvePoly."); + return NULL; + } + } + + ogeom = lwpoly_construct(curvepoly->SRID, NULL, curvepoly->nrings, ptarray); + return ogeom; +} + +LWMLINE * +lwmcurve_segmentize(LWMCURVE *mcurve, uint32 perQuad) +{ + LWMLINE *ogeom; + LWGEOM *tmp; + LWGEOM **lines; + int i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwmcurve_segmentize called, geoms=%d, dim=%d.", mcurve->ngeoms, TYPE_NDIMS(mcurve->type)); +#endif + + lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms); + + for(i = 0; i < mcurve->ngeoms; i++) + { + tmp = mcurve->geoms[i]; + if(lwgeom_getType(tmp->type) == CURVETYPE) + { + lines[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad); + } + else if(lwgeom_getType(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(LWMSURFACE *msurface, uint32 perQuad) +{ + LWMPOLY *ogeom; + LWGEOM *tmp; + LWPOLY *poly; + LWGEOM **polys; + POINTARRAY **ptarray; + int i, j; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwmsurface_segmentize called."); +#endif + + polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms); + + for(i = 0; i < msurface->ngeoms; i++) + { + tmp = msurface->geoms[i]; + if(lwgeom_getType(tmp->type) == CURVEPOLYTYPE) + { + polys[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad); + } + else if(lwgeom_getType(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(LWCOLLECTION *collection, uint32 perQuad) +{ + LWCOLLECTION *ocol; + LWGEOM *tmp; + LWGEOM **geoms; + int i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwcollection_segmentize called."); +#endif + + if(has_arc((LWGEOM *)collection) == 0) + { + return collection; + } + + geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms); + + for(i=0; ingeoms; i++) + { + tmp = collection->geoms[i]; + switch(lwgeom_getType(tmp->type)) { + case CURVETYPE: + geoms[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad); + break; + case COMPOUNDTYPE: + geoms[i] = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad); + break; + case CURVEPOLYTYPE: + geoms[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)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(LWGEOM *geom, uint32 perQuad) +{ + LWGEOM * ogeom = NULL; + switch(lwgeom_getType(geom->type)) { + case CURVETYPE: + ogeom = (LWGEOM *)lwcurve_segmentize((LWCURVE *)geom, perQuad); + break; + case COMPOUNDTYPE: + ogeom = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)geom, perQuad); + break; + case CURVEPOLYTYPE: + ogeom = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)geom, perQuad); + break; + case MULTICURVETYPE: + ogeom = (LWGEOM *)lwmcurve_segmentize((LWMCURVE *)geom, perQuad); + break; + case MULTISURFACETYPE: + ogeom = (LWGEOM *)lwmsurface_segmentize((LWMSURFACE *)geom, perQuad); + break; + case COLLECTIONTYPE: + ogeom = (LWGEOM *)lwcollection_segmentize((LWCOLLECTION *)geom, perQuad); + break; + default: + ogeom = lwgeom_clone(geom); + } + return ogeom; +} + +/******************************************************************************* + * End curve segmentize functions + ******************************************************************************/ +LWGEOM * +append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int SRID) +{ + LWGEOM *result; + int currentType, i; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("append_segment called %p, %p, %d, %d", geom, pts, type, SRID); +#endif + + if(geom == NULL) + { + if(type == LINETYPE) + { +#ifdef PGIS_DEBUG + lwnotice("append_segment: line to NULL"); +#endif + return (LWGEOM *)lwline_construct(SRID, NULL, pts); + } + else if(type == CURVETYPE) + { +#ifdef PGIS_DEBUG + lwnotice("append_segment: curve to NULL %d", pts->npoints); + POINT4D tmp; + for(i=0; inpoints; i++) + { + getPoint4d_p(pts, i, &tmp); + lwnotice("new point: (%.16f,%.16f)",tmp.x,tmp.y); + } +#endif + return (LWGEOM *)lwcurve_construct(SRID, NULL, pts); + } + else + { + lwerror("Invalid segment type %d.", type); + } + } + + currentType = lwgeom_getType(geom->type); + if(currentType == LINETYPE && type == LINETYPE) + { + POINTARRAY *newPoints; + POINT4D pt; + LWLINE *line = (LWLINE *)geom; +#ifdef PGIS_DEBUG + lwnotice("append_segment: line to line"); +#endif + newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + line->points->npoints - 1); + for(i=0; ipoints->npoints; i++) + { + getPoint4d_p(pts, i, &pt); + setPoint4d(newPoints, i, &pt); + } + for(i=1; inpoints; i++) + { + getPoint4d_p(line->points, i, &pt); + setPoint4d(newPoints, i + pts->npoints, &pt); + } + result = (LWGEOM *)lwline_construct(SRID, NULL, newPoints); + lwgeom_release(geom); + return result; + } + else if(currentType == CURVETYPE && type == CURVETYPE) + { + POINTARRAY *newPoints; + POINT4D pt; + LWCURVE *curve = (LWCURVE *)geom; +#ifdef PGIS_DEBUG + lwnotice("append_segment: curve to curve"); +#endif + newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + curve->points->npoints - 1); +#ifdef PGIS_DEBUG + lwnotice("New array length: %d", pts->npoints + curve->points->npoints - 1); +#endif + for(i=0; ipoints->npoints; i++) + { + getPoint4d_p(curve->points, i, &pt); +#ifdef PGIS_DEBUG + lwnotice("orig point %d: (%.16f,%.16f)", i, pt.x, pt.y); +#endif + setPoint4d(newPoints, i, &pt); + } + for(i=1; inpoints;i++) + { + getPoint4d_p(pts, i, &pt); +#ifdef PGIS_DEBUG + lwnotice("new point %d: (%.16f,%.16f)", i + curve->points->npoints - 1, pt.x, pt.y); +#endif + setPoint4d(newPoints, i + curve->points->npoints - 1, &pt); + } + result = (LWGEOM *)lwcurve_construct(SRID, NULL, newPoints); + lwgeom_release(geom); + return result; + } + else if(currentType == CURVETYPE && type == LINETYPE) + { + LWLINE *line; + LWGEOM **geomArray; +#ifdef PGIS_DEBUG + lwnotice("append_segment: line to curve"); +#endif + geomArray = lwalloc(sizeof(LWGEOM *)*2); + geomArray[0] = lwgeom_clone(geom); + + line = lwline_construct(SRID, NULL, pts); + geomArray[1] = lwgeom_clone((LWGEOM *)line); + + result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray); + lwfree((LWGEOM *)line); + lwgeom_release(geom); + return result; + } + else if(currentType == LINETYPE && type == CURVETYPE) + { + LWCURVE *curve; + LWGEOM **geomArray; +#ifdef PGIS_DEBUG + lwnotice("append_segment: curve to line"); +#endif + geomArray = lwalloc(sizeof(LWGEOM *)*2); + geomArray[0] = lwgeom_clone(geom); + + curve = lwcurve_construct(SRID, NULL, pts); + geomArray[1] = lwgeom_clone((LWGEOM *)curve); + + result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray); + lwfree((LWGEOM *)curve); + lwgeom_release(geom); + return result; + } + else if(currentType == COMPOUNDTYPE) + { + LWGEOM *newGeom; + LWCOMPOUND *compound; + int count; + LWGEOM **geomArray; + + compound = (LWCOMPOUND *)geom; + count = compound->ngeoms + 1; + geomArray = lwalloc(sizeof(LWGEOM *)*count); + for(i=0; ingeoms; i++) + { + geomArray[i] = lwgeom_clone(compound->geoms[i]); + } + if(type == LINETYPE) + { +#ifdef PGIS_DEBUG + lwnotice("append_segment: line to compound"); +#endif + newGeom = (LWGEOM *)lwline_construct(SRID, NULL, pts); + } + else if(type == CURVETYPE) + { +#ifdef PGIS_DEBUG + lwnotice("append_segment: curve to compound"); +#endif + newGeom = (LWGEOM *)lwcurve_construct(SRID, NULL, pts); + } + else + { + lwerror("Invalid segment type %d.", type); + return NULL; + } + geomArray[compound->ngeoms] = lwgeom_clone(newGeom); + + result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, count, geomArray); + lwfree(newGeom); + lwgeom_release(geom); + return result; + } + lwerror("Invalid state %d-%d", currentType, type); + return NULL; +} + +LWGEOM * +pta_desegmentize(POINTARRAY *points, int type, int SRID) +{ + int i, j, commit, isline, count; + double last_angle, last_length; + double dxab, dyab, dxbc, dybc, theta, length; + POINT4D a, b, c, tmp; + POINTARRAY *pts; + LWGEOM *geom = NULL; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("pta_desegmentize called."); +#endif + + getPoint4d_p(points, 0, &a); + getPoint4d_p(points, 1, &b); + getPoint4d_p(points, 2, &c); + + dxab = b.x - a.x; + dyab = b.y - a.y; + dxbc = c.x - b.x; + dybc = c.y - b.y; + + theta = atan2(dyab, dxab); + last_angle = theta - atan2(dybc, dxbc); + last_length = sqrt(dxbc*dxbc+dybc*dybc); + length = sqrt(dxab*dxab+dyab*dyab); + if((last_length - length) < EPSILON_SQLMM) + { + isline = -1; +#ifdef PGIS_DEBUG + lwnotice("Starting as unknown."); +#endif + } + else + { + isline = 1; +#ifdef PGIS_DEBUG + lwnotice("Starting as line."); +#endif + } + + commit = 0; + for(i=3; inpoints; i++) + { + getPoint4d_p(points, i-2, &a); + getPoint4d_p(points, i-1, &b); + getPoint4d_p(points, i, &c); + + dxab = b.x - a.x; + dyab = b.y - a.y; + dxbc = c.x - b.x; + dybc = c.y - b.y; + +#ifdef PGIS_DEBUG + lwnotice("(dxab, dyab, dxbc, dybc) (%.16f, %.16f, %.16f, %.16f)", dxab, dyab, dxbc, dybc); +#endif + + theta = atan2(dyab, dxab); + theta = theta - atan2(dybc, dxbc); + length = sqrt(dxbc*dxbc+dybc*dybc); +#ifdef PGIS_DEBUG + lwnotice("Last/current length and angle %.16f/%.16f, %.16f/%.16f", last_angle, theta, last_length, length); +#endif + /* Found a line segment */ + if(fabs(length - last_length) > EPSILON_SQLMM || + fabs(theta - last_angle) > EPSILON_SQLMM) + { + last_length = length; + last_angle = theta; + /* We are tracking a line, keep going */ + if(isline > 0) + { + } + /* We were tracking a curve, commit it and start line*/ + else if(isline == 0) + { +#ifdef PGIS_DEBUG + lwnotice("Building curve, %d - %d", commit, i); +#endif + count = i - commit; + pts = ptarray_construct( + TYPE_HASZ(type), + TYPE_HASM(type), + 3); + getPoint4d_p(points, commit, &tmp); + setPoint4d(pts, 0, &tmp); + getPoint4d_p(points, + commit + count/2, &tmp); + setPoint4d(pts, 1, &tmp); + getPoint4d_p(points, i - 1, &tmp); + setPoint4d(pts, 2, &tmp); + + commit = i-1; + geom = append_segment(geom, pts, CURVETYPE, SRID); + isline = -1; + + /* + * We now need to move ahead one point to + * determine if it's a potential new curve, + * since the last_angle value is corrupt. + */ + i++; + getPoint4d_p(points, i-2, &a); + getPoint4d_p(points, i-1, &b); + getPoint4d_p(points, i, &c); + + dxab = b.x - a.x; + dyab = b.y - a.y; + dxbc = c.x - b.x; + dybc = c.y - b.y; + + theta = atan2(dyab, dxab); + last_angle = theta - atan2(dybc, dxbc); + last_length = sqrt(dxbc*dxbc+dybc*dybc); + length = sqrt(dxab*dxab+dyab*dyab); + if((last_length - length) < EPSILON_SQLMM) + { + isline = -1; +#ifdef PGIS_DEBUG + lwnotice("Restarting as unknown."); +#endif + } + else + { + isline = 1; +#ifdef PGIS_DEBUG + lwnotice("Restarting as line."); +#endif + } + + + } + /* We didn't know what we were tracking, now we do. */ + else + { +#ifdef PGIS_DEBUG + lwnotice("It's a line"); +#endif + isline = 1; + } + } + /* Found a curve segment */ + else + { + /* We were tracking a curve, commit it and start line */ + if(isline > 0) + { +#ifdef PGIS_DEBUG + lwnotice("Building line, %d - %d", commit, i-2); +#endif + count = i - commit - 2; + + pts = ptarray_construct( + TYPE_HASZ(type), + TYPE_HASM(type), + count); + for(j=commit;j 2) + { +#ifdef PGIS_DEBUG + lwnotice("Finishing curve %d,%d.", commit, i); +#endif + pts = ptarray_construct( + TYPE_HASZ(type), + TYPE_HASM(type), + 3); + getPoint4d_p(points, commit, &tmp); + setPoint4d(pts, 0, &tmp); + getPoint4d_p(points, commit + count/2, &tmp); + setPoint4d(pts, 1, &tmp); + getPoint4d_p(points, i - 1, &tmp); + setPoint4d(pts, 2, &tmp); + + geom = append_segment(geom, pts, CURVETYPE, SRID); + } + else + { +#ifdef PGIS_DEBUG + lwnotice("Finishing line %d,%d.", commit, i); +#endif + pts = ptarray_construct( + TYPE_HASZ(type), + TYPE_HASM(type), + count); + for(j=commit;jpoints, line->type, line->SRID); +} + +LWGEOM * +lwpolygon_desegmentize(LWPOLY *poly) +{ + LWGEOM **geoms; + int i, hascurve = 0; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwpolygon_desegmentize called."); +#endif + + geoms = lwalloc(sizeof(LWGEOM *)*poly->nrings); + for(i=0; inrings; i++) + { + geoms[i] = pta_desegmentize(poly->rings[i], poly->type, poly->SRID); + if(lwgeom_getType(geoms[i]->type) == CURVETYPE || + lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE) + { + hascurve = 1; + } + } + if(hascurve == 0) + { + for(i=0; inrings; i++) + { + lwfree(geoms[i]); + } + return lwgeom_clone((LWGEOM *)poly); + } + + return (LWGEOM *)lwcollection_construct(CURVEPOLYTYPE, poly->SRID, NULL, poly->nrings, geoms); +} + +LWGEOM * +lwmline_desegmentize(LWMLINE *mline) +{ + LWGEOM **geoms; + int i, hascurve = 0; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwmline_desegmentize called."); +#endif + + geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms); + for(i=0; ingeoms; i++) + { + geoms[i] = lwline_desegmentize((LWLINE *)mline->geoms[i]); + if(lwgeom_getType(geoms[i]->type) == CURVETYPE || + lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE) + { + hascurve = 1; + } + } + if(hascurve == 0) + { + for(i=0; ingeoms; i++) + { + lwfree(geoms[i]); + } + return lwgeom_clone((LWGEOM *)mline); + } + return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->SRID, NULL, mline->ngeoms, geoms); +} + +LWGEOM * +lwmpolygon_desegmentize(LWMPOLY *mpoly) +{ + LWGEOM **geoms; + int i, hascurve = 0; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwmpoly_desegmentize called."); +#endif + + geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms); + for(i=0; ingeoms; i++) + { + geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]); + if(lwgeom_getType(geoms[i]->type) == CURVEPOLYTYPE) + { + hascurve = 1; + } + } + if(hascurve == 0) + { + for(i=0; ingeoms; i++) + { + lwfree(geoms[i]); + } + return lwgeom_clone((LWGEOM *)mpoly); + } + return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->SRID, NULL, mpoly->ngeoms, geoms); +} + +LWGEOM * +lwgeom_desegmentize(LWGEOM *geom) +{ + int type = lwgeom_getType(geom->type); + +#ifdef PGIS_DEBUG_CALLS + lwnotice("lwgeom_desegmentize called."); +#endif + + switch(type) { + case LINETYPE: + return lwline_desegmentize((LWLINE *)geom); + case POLYGONTYPE: + return lwpolygon_desegmentize((LWPOLY *)geom); + case MULTILINETYPE: + return lwmline_desegmentize((LWMLINE *)geom); + case MULTIPOLYGONTYPE: + return lwmpolygon_desegmentize((LWMPOLY *)geom); + default: + return lwgeom_clone(geom); + } +} + +/******************************************************************************* + * Begin PG_FUNCTIONs + ******************************************************************************/ + +PG_FUNCTION_INFO_V1(LWGEOM_has_arc); +Datum LWGEOM_has_arc(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + uint32 result = has_arc(lwgeom_deserialize(SERIALIZED_FORM(geom))); + PG_RETURN_BOOL(result == 1); +} + +/* + * Converts any curve segments of the geometry into a linear approximation. + * Curve centers are determined by projecting the defining points into the 2d + * plane. Z and M values are assigned by linear interpolation between + * defining points. + */ +PG_FUNCTION_INFO_V1(LWGEOM_curve_segmentize); +Datum LWGEOM_curve_segmentize(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + uint32 perQuad = PG_GETARG_INT32(1); + PG_LWGEOM *ret; + LWGEOM *igeom = NULL, *ogeom = NULL; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("LWGEOM_curve_segmentize called."); +#endif + + if(perQuad < 0) + { + elog(ERROR, "2nd argument must be positive."); + PG_RETURN_NULL(); + } +#ifdef PGIS_DEBUG + else + { + lwnotice("perQuad = %d", perQuad); + } +#endif + igeom = lwgeom_deserialize(SERIALIZED_FORM(geom)); + ogeom = lwgeom_segmentize(igeom, perQuad); + if(ogeom == NULL) PG_RETURN_NULL(); + ret = pglwgeom_serialize(ogeom); + lwgeom_release(igeom); + lwgeom_release(ogeom); + PG_FREE_IF_COPY(geom, 0); + PG_RETURN_POINTER(ret); +} + +PG_FUNCTION_INFO_V1(LWGEOM_line_desegmentize); +Datum LWGEOM_line_desegmentize(PG_FUNCTION_ARGS) +{ + PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0)); + PG_LWGEOM *ret; + LWGEOM *igeom = NULL, *ogeom = NULL; + +#ifdef PGIS_DEBUG_CALLS + lwnotice("LWGEOM_line_desegmentize."); +#endif + + igeom = lwgeom_deserialize(SERIALIZED_FORM(geom)); + ogeom = lwgeom_desegmentize(igeom); + if(ogeom == NULL) + { + lwgeom_release(igeom); + PG_RETURN_NULL(); + } + ret = pglwgeom_serialize(ogeom); + lwgeom_release(igeom); + lwgeom_release(ogeom); + PG_FREE_IF_COPY(geom, 0); + PG_RETURN_POINTER(ret); +} + +/******************************************************************************* + * End PG_FUNCTIONs + ******************************************************************************/ diff --git a/lwgeom/lwgeom_transform.c b/lwgeom/lwgeom_transform.c index dd4654cf0..3fac1e8e8 100644 --- a/lwgeom/lwgeom_transform.c +++ b/lwgeom/lwgeom_transform.c @@ -28,6 +28,7 @@ Datum transform(PG_FUNCTION_ARGS); Datum transform_geom(PG_FUNCTION_ARGS); Datum postgis_proj_version(PG_FUNCTION_ARGS); + /* if POSTGIS_PROJ_VERSION undefined, we get a do-nothing transform() function */ #ifndef POSTGIS_PROJ_VERSION @@ -145,7 +146,7 @@ void DeleteFromPROJ4SRSCache(PROJ4PortalCache *PROJ4Cache, int srid); /* Search path for PROJ.4 library */ static bool IsPROJ4LibPathSet = false; -void SetPROJ4LibPath(); +void SetPROJ4LibPath(void); /* Memory context cache functions */ static void PROJ4SRSCacheInit(MemoryContext context); @@ -244,7 +245,7 @@ PROJ4SRSCacheReset(MemoryContext context) { /* * Do nothing, but we must supply a function since this call is mandatory according to tgl - * (see postgis-devel archives July 2007) + * (see postgis-devel archives July 2007) */ } @@ -253,7 +254,7 @@ PROJ4SRSCacheIsEmpty(MemoryContext context) { /* * Always return false since this call is mandatory according to tgl - * (see postgis-devel archives July 2007) + * (see postgis-devel archives July 2007) */ return FALSE; } @@ -263,7 +264,7 @@ PROJ4SRSCacheStats(MemoryContext context) { /* * Simple stats display function - we must supply a function since this call is mandatory according to tgl - * (see postgis-devel archives July 2007) + * (see postgis-devel archives July 2007) */ fprintf(stderr, "%s: PROJ4 context\n", context->name); @@ -580,7 +581,7 @@ void DeleteFromPROJ4SRSCache(PROJ4PortalCache *PROJ4Cache, int srid) * since the method of determining the current installation * path are different on older PostgreSQL versions. */ -void SetPROJ4LibPath() +void SetPROJ4LibPath(void) { #if POSTGIS_PGSQL_VERSION >= 80 char *path; diff --git a/lwgeom/wktunparse.c b/lwgeom/wktunparse.c index 152fa7f49..b6dffed31 100644 --- a/lwgeom/wktunparse.c +++ b/lwgeom/wktunparse.c @@ -42,6 +42,8 @@ uchar* output_single(uchar* geom,int supress); uchar* output_collection(uchar* geom,outfunc func,int supress); uchar* output_collection_2(uchar* geom,int suppress); uchar* output_multipoint(uchar* geom,int suppress); +uchar* output_compound(uchar* geom, int suppress); +uchar* output_multisurface(uchar* geom, int suppress); void write_wkb_hex_bytes(uchar* ptr, unsigned int cnt, size_t size); void write_wkb_bin_bytes(uchar* ptr, unsigned int cnt, size_t size);