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))
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.
* 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);
/*
#define ABORT_ON_AUTH_FAILURE 1
Datum check_authorization(PG_FUNCTION_ARGS);
+Datum getTransactionID(PG_FUNCTION_ARGS);
/*
* This trigger will check for authorization before
-/**********************************************************************\r
- * $Id$\r
- *\r
- * PostGIS - Spatial Types for PostgreSQL\r
- * http://postgis.refractions.net\r
- * Copyright 2001-2006 Refractions Research Inc.\r
- *\r
- * This is free software; you can redistribute and/or modify it under\r
- * the terms of the GNU General Public Licence. See the COPYING file.\r
- * \r
- **********************************************************************/\r
-\r
-/* basic LWCURVE functions */\r
-\r
-#include <stdio.h>\r
-#include <stdlib.h>\r
-#include <string.h>\r
-#include <math.h>\r
-#include "liblwgeom.h"\r
-\r
-/*#define PGIS_DEBUG_CALLS 1 */\r
-/*#define PGIS_DEBUG 1 */\r
-\r
-#ifndef MAXFLOAT\r
- #define MAXFLOAT 3.402823466e+38F\r
-#endif\r
-\r
-/*\r
- * Construct a new LWCURVE. points will *NOT* be copied\r
- * use SRID=-1 for unknown SRID (will have 8bit type's S = 0)\r
- */\r
-LWCURVE *\r
-lwcurve_construct(int SRID, BOX2DFLOAT4 *bbox, POINTARRAY *points)\r
-{\r
- LWCURVE *result;\r
- \r
- /*\r
- * The first arc requires three points. Each additional\r
- * arc requires two more points. Thus the minimum point count\r
- * is three, and the count must be odd.\r
- */\r
- if(points->npoints % 2 != 1 || points->npoints < 3) \r
- {\r
- lwerror("lwcurve_construct: invalid point count %d", points->npoints);\r
- return NULL;\r
- }\r
- \r
- result = (LWCURVE*) lwalloc(sizeof(LWCURVE));\r
-\r
- result->type = lwgeom_makeType_full(\r
- TYPE_HASZ(points->dims),\r
- TYPE_HASM(points->dims),\r
- (SRID!=-1), CURVETYPE, 0);\r
- result->SRID = SRID;\r
- result->points = points;\r
- result->bbox = bbox;\r
- \r
- return result;\r
-}\r
-\r
-/*\r
- * given the LWGEOM serialized form (or a point into a multi* one)\r
- * construct a propert LWCURVE.\r
- * serialized_form should point to the 8bit type format (with type = 8)\r
- * See serialized form doc\r
- */\r
-LWCURVE *\r
-lwcurve_deserialize(uchar *serialized_form)\r
-{\r
- uchar type;\r
- LWCURVE *result;\r
- uchar *loc=NULL;\r
- uint32 npoints;\r
- POINTARRAY *pa;\r
-\r
- type = (uchar)serialized_form[0];\r
- if(lwgeom_getType(type) != CURVETYPE)\r
- {\r
- lwerror("lwcurve_deserialize: attempt to deserialize a curve which is really a %s", lwgeom_typename(type));\r
- return NULL;\r
- }\r
-\r
- result = (LWCURVE*) lwalloc(sizeof(LWCURVE));\r
- result->type = type;\r
-\r
- loc = serialized_form + 1;\r
-\r
- if(lwgeom_hasBBOX(type))\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_deserialize: input has bbox");\r
-#endif\r
- result->bbox = lwalloc(sizeof(BOX2DFLOAT4));\r
- memcpy(result->bbox, loc, sizeof(BOX2DFLOAT4));\r
- loc += sizeof(BOX2DFLOAT4); \r
- }\r
- else\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_deserialize: input lacks bbox"); \r
-#endif \r
- result->bbox = NULL;\r
- }\r
-\r
- if(lwgeom_hasSRID(type))\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_deserialize: input has srid");\r
-#endif\r
- result->SRID = lw_get_int32(loc); \r
- loc += 4; /* type + SRID */\r
- }\r
- else\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_deserialize: input lacks srid");\r
-#endif\r
- result->SRID = -1; \r
- }\r
-\r
- /* we've read the type (1 byte) and SRID (4 bytes, if present) */\r
-\r
- npoints = lw_get_uint32(loc);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("curve npoints = %d", npoints);\r
-#endif \r
- loc += 4;\r
- pa = pointArray_construct(loc, TYPE_HASZ(type), TYPE_HASM(type), npoints);\r
- result->points = pa;\r
- return result;\r
-}\r
-\r
-/*\r
- * convert this curve into its serialized form\r
- * result's first char will be the 8bit type. See serialized form doc\r
- */\r
-uchar *\r
-lwcurve_serialize(LWCURVE *curve)\r
-{\r
- size_t size, retsize;\r
- uchar * result;\r
-\r
- if(curve == NULL) {\r
- lwerror("lwcurve_serialize:: given null curve");\r
- return NULL;\r
- }\r
-\r
- size = lwcurve_serialize_size(curve);\r
- result = lwalloc(size);\r
- lwcurve_serialize_buf(curve, result, &retsize);\r
- if(retsize != size)\r
- lwerror("lwcurve_serialize_size returned %d, ..selialize_buf returned %d", size, retsize);\r
- return result;\r
-}\r
-\r
-/*\r
- * convert this curve into its serialized form writing it into \r
- * the given buffer, and returning number of bytes written into \r
- * the given int pointer.\r
- * result's first char will be the 8bit type. See serialized form doc\r
- */\r
-void lwcurve_serialize_buf(LWCURVE *curve, uchar *buf, size_t *retsize)\r
-{\r
- char hasSRID;\r
- uchar *loc;\r
- int ptsize;\r
- size_t size;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_serialize_buf(%p, %p, %p) called",\r
- curve, buf, retsize);\r
-#endif\r
-\r
- if(curve == NULL) \r
- {\r
- lwerror("lwcurve_serialize:: given null curve");\r
- return;\r
- }\r
-\r
- if(TYPE_GETZM(curve->type) != TYPE_GETZM(curve->points->dims))\r
- {\r
- lwerror("Dimensions mismatch in lwcurve");\r
- return;\r
- }\r
-\r
- ptsize = pointArray_ptsize(curve->points);\r
-\r
- hasSRID = (curve->SRID != -1);\r
-\r
- buf[0] = (uchar)lwgeom_makeType_full(\r
- TYPE_HASZ(curve->type), TYPE_HASM(curve->type),\r
- hasSRID, CURVETYPE, curve->bbox ? 1 : 0);\r
- loc = buf+1;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_serialize_buf added type (%d)", curve->type);\r
-#endif\r
-\r
- if(curve->bbox)\r
- {\r
- memcpy(loc, curve->bbox, sizeof(BOX2DFLOAT4));\r
- loc += sizeof(BOX2DFLOAT4);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_serialize_buf added BBOX");\r
-#endif\r
-\r
- } \r
-\r
- if(hasSRID)\r
- {\r
- memcpy(loc, &curve->SRID, sizeof(int32));\r
- loc += sizeof(int32);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_serialize_buf added SRID");\r
-#endif\r
-\r
- }\r
-\r
- memcpy(loc, &curve->points->npoints, sizeof(uint32));\r
- loc += sizeof(uint32);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_serialize_buf added npoints (%d)",\r
- curve->points->npoints);\r
-#endif\r
-\r
- /* copy in points */\r
- size = curve->points->npoints * ptsize;\r
- memcpy(loc, getPoint_internal(curve->points, 0), size);\r
- loc += size;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_serialize_buf copied serialized_pointlist (%d bytes)",\r
- ptsize * curve->points->npoints); \r
-#endif\r
-\r
- if(retsize) *retsize = loc-buf;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_serialize_buf returning (loc: %p, size: %d)",\r
- loc, loc-buf);\r
-#endif\r
-}\r
-\r
-/* find length of this deserialized curve */\r
-size_t\r
-lwcurve_serialize_size(LWCURVE *curve)\r
-{\r
- size_t size = 1; /* type */\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_serialize_size called"); \r
-#endif\r
-\r
- if(curve->SRID != -1) size += 4; /* SRID */\r
- if(curve->bbox) size += sizeof(BOX2DFLOAT4);\r
-\r
- size += 4; /* npoints */\r
- size += pointArray_ptsize(curve->points) * curve->points->npoints;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_serialize_size returning %d", size);\r
-#endif\r
-\r
- return size;\r
-}\r
-\r
-BOX3D *\r
-lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)\r
-{\r
- double x1, x2, y1, y2, z1, z2;\r
- double angle, radius, sweep;\r
- /*\r
- double top, left;\r
- */\r
- double a1, a2, a3;\r
- double xe = 0.0, ye = 0.0;\r
- POINT4D *center;\r
- int i;\r
- BOX3D *box;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcircle_compute_box3d called.");\r
-#endif\r
-\r
- center = lwalloc(sizeof(POINT4D));\r
- radius = lwcircle_center(p1, p2, p3, ¢er);\r
- if(radius < 0.0) return NULL;\r
-\r
- /*\r
- top = center->y + radius;\r
- left = center->x - radius;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcircle_compute_box3d: top=%.16f, left=%.16f", top, left);\r
-#endif\r
- */\r
-\r
- x1 = MAXFLOAT;\r
- x2 = -1 * MAXFLOAT;\r
- y1 = MAXFLOAT;\r
- y2 = -1 * MAXFLOAT;\r
-\r
- a1 = atan2(p1->y - center->y, p1->x - center->x);\r
- a2 = atan2(p2->y - center->y, p2->x - center->x);\r
- a3 = atan2(p3->y - center->y, p3->x - center->x);\r
-\r
- /* Determine sweep angle */\r
- if(a1 > a2 && a2 > a3) \r
- {\r
- sweep = a3 - a1;\r
- }\r
- /* Counter-clockwise */\r
- else if(a1 < a2 && a2 < a3)\r
- {\r
- sweep = a3 - a1;\r
- }\r
- /* Clockwise, wrap */\r
- else if((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3))\r
- {\r
- sweep = a3 - a1 + 2*M_PI;\r
- }\r
- /* Counter-clockwise, wrap */\r
- else if((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3))\r
- {\r
- sweep = a3 - a1 - 2*M_PI;\r
- } \r
- else \r
- {\r
- sweep = 0.0;\r
- }\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep);\r
-#endif\r
-\r
- angle = 0.0;\r
- for(i=0; i < 6; i++)\r
- {\r
- switch(i) {\r
- case 0:\r
- angle = 0.0;\r
- xe = center->x + radius;\r
- ye = center->y;\r
- break;\r
- case 1:\r
- angle = M_PI_2;\r
- xe = center->x;\r
- ye = center->y + radius;\r
- break;\r
- case 2:\r
- angle = M_PI;\r
- xe = center->x - radius;\r
- ye = center->y;\r
- break;\r
- case 3:\r
- angle = -1 * M_PI_2;\r
- xe = center->x;\r
- ye = center->y - radius;\r
- break;\r
- case 4:\r
- angle = a1;\r
- xe = p1->x;\r
- ye = p1->y;\r
- break;\r
- case 5:\r
- angle = a3;\r
- xe = p3->x;\r
- ye = p3->y;\r
- break;\r
- }\r
- if(i < 4) \r
- {\r
- if(sweep > 0.0 && (angle > a3 || angle < a1)) continue;\r
- if(sweep < 0.0 && (angle < a3 || angle > a1)) continue;\r
- }\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcircle_compute_box3d: potential extreame %d (%.16f, %.16f)", i, xe, ye);\r
-#endif\r
- x1 = (x1 < xe) ? x1 : xe;\r
- y1 = (y1 < ye) ? y1 : ye;\r
- x2 = (x2 > xe) ? x2 : xe;\r
- y2 = (y2 > ye) ? y2 : ye;\r
- }\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcircle_compute_box3d: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2);\r
-#endif\r
-\r
- /*\r
- x1 = center->x + x1 * radius;\r
- x2 = center->x + x2 * radius;\r
- y1 = center->y + y1 * radius;\r
- y2 = center->y + y2 * radius;\r
- */\r
- z1 = (p1->z < p2->z) ? p1->z : p2->z;\r
- z1 = (z1 < p3->z) ? z1 : p3->z;\r
- z2 = (p1->z > p2->z) ? p1->z : p2->z;\r
- z2 = (z2 > p3->z) ? z2 : p3->z;\r
-\r
- box = lwalloc(sizeof(BOX3D));\r
- box->xmin = x1; box->xmax = x2;\r
- box->ymin = y1; box->ymax = y2;\r
- box->zmin = z1; box->zmax = z2;\r
-\r
- lwfree(center);\r
-\r
- return box;\r
-}\r
-\r
-/*\r
- * Find bounding box (standard one)\r
- * zmin=zmax=NO_Z_VALUE if 2d\r
- * TODO: This ignores curvature, which should be taken into account.\r
- */\r
-BOX3D *\r
-lwcurve_compute_box3d(LWCURVE *curve)\r
-{\r
- BOX3D *box, *tmp; \r
- int i;\r
- POINT4D *p1 = lwalloc(sizeof(POINT4D));\r
- POINT4D *p2 = lwalloc(sizeof(POINT4D));\r
- POINT4D *p3 = lwalloc(sizeof(POINT4D));\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_compute_box3d called.");\r
-#endif\r
-\r
- /* initialize box values */\r
- box = lwalloc(sizeof(BOX3D));\r
- box->xmin = MAXFLOAT; box->xmax = -1 * MAXFLOAT;\r
- box->ymin = MAXFLOAT; box->ymax = -1 * MAXFLOAT;\r
- box->zmin = MAXFLOAT; box->zmax = -1 * MAXFLOAT;\r
- \r
- for(i = 2; i < curve->points->npoints; i+=2)\r
- {\r
- getPoint4d_p(curve->points, i-2, p1);\r
- getPoint4d_p(curve->points, i-1, p2);\r
- getPoint4d_p(curve->points, i, p3);\r
- tmp = lwcircle_compute_box3d(p1, p2, p3);\r
- if(tmp == NULL) continue;\r
- box->xmin = (box->xmin < tmp->xmin) ? box->xmin : tmp->xmin;\r
- box->xmax = (box->xmax > tmp->xmax) ? box->xmax : tmp->xmax;\r
- box->ymin = (box->ymin < tmp->ymin) ? box->ymin : tmp->ymin;\r
- box->ymax = (box->ymax > tmp->ymax) ? box->ymax : tmp->ymax;\r
- box->zmin = (box->zmin < tmp->zmin) ? box->zmin : tmp->zmin;\r
- box->zmax = (box->zmax > tmp->zmax) ? box->zmax : tmp->zmax;\r
-#ifdef PGIS_DEBUG_CALLS\r
- 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);\r
-#endif\r
- }\r
-\r
- \r
- return box;\r
-}\r
-\r
-int\r
-lwcurve_compute_box2d_p(LWCURVE *curve, BOX2DFLOAT4 *result)\r
-{\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurve_compute_box2d_p called.");\r
-#endif\r
-\r
- BOX3D *box = lwcurve_compute_box3d(curve);\r
- if(box == NULL) return 0;\r
- box3d_to_box2df_p(box, result);\r
- return 1;\r
-}\r
-\r
-void pfree_curve(LWCURVE *curve)\r
-{\r
- lwfree(curve->points);\r
- lwfree(curve);\r
-}\r
-\r
-/* find length of this serialized curve */\r
-size_t\r
-lwgeom_size_curve(const uchar *serialized_curve)\r
-{\r
- int type = (uchar)serialized_curve[0];\r
- uint32 result = 1; /* type */\r
- const uchar *loc;\r
- uint32 npoints;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwgeom_size_curve called");\r
-#endif\r
- if(lwgeom_getType(type) != CURVETYPE)\r
- lwerror("lwgeom_size_curve::attempt to find the length of a non-curve");\r
-\r
- loc = serialized_curve + 1;\r
- if(lwgeom_hasBBOX(type))\r
- {\r
- loc += sizeof(BOX2DFLOAT4);\r
- result += sizeof(BOX2DFLOAT4);\r
- }\r
-\r
- if(lwgeom_hasSRID(type))\r
- {\r
- loc += 4; /* type + SRID */\r
- result += 4;\r
- }\r
-\r
- /* we've read the type (1 byte) and SRID (4 bytes, if present) */\r
- npoints = lw_get_uint32(loc);\r
- result += sizeof(uint32); /* npoints */\r
-\r
- result += TYPE_NDIMS(type) * sizeof(double) * npoints;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwgeom_size_curve returning %d", result);\r
-#endif \r
-\r
- return result;\r
-}\r
-\r
-void printLWCURVE(LWCURVE *curve)\r
-{\r
- lwnotice("LWCURVE {");\r
- lwnotice(" ndims = %i", (int)TYPE_NDIMS(curve->type));\r
- lwnotice(" SRID = %i", (int)curve->SRID);\r
- printPA(curve->points);\r
- lwnotice("}");\r
-}\r
-\r
-/* Clone LWCURVE object. POINTARRAY is not copied. */\r
-LWCURVE *\r
-lwcurve_clone(const LWCURVE *g)\r
-{\r
- LWCURVE *ret = lwalloc(sizeof(LWCURVE));\r
- memcpy(ret, g, sizeof(LWCURVE));\r
- if(g->bbox) ret->bbox = box2d_clone(g->bbox);\r
- return ret;\r
-}\r
-\r
-/*\r
- * Add 'what' to this curve at position 'where'.\r
- * where=0 == prepend\r
- * where=-1 == append\r
- * Returns a MULTICURVE or a GEOMETRYCOLLECTION\r
- */\r
-LWGEOM *\r
-lwcurve_add(const LWCURVE *to, uint32 where, const LWGEOM *what)\r
-{\r
- LWCOLLECTION *col;\r
- LWGEOM **geoms;\r
- int newtype;\r
- \r
- if(where != -1 && where != 0)\r
- {\r
- lwerror("lwcurve_add only supports 0 or -1 as second argument %d", where);\r
- return NULL;\r
- }\r
-\r
- /* dimensions compatibility are checked by caller */\r
-\r
- /* Construct geoms array */\r
- geoms = lwalloc(sizeof(LWGEOM *)*2);\r
- if(where == -1) /* append */\r
- {\r
- geoms[0] = lwgeom_clone((LWGEOM *)to);\r
- geoms[1] = lwgeom_clone(what);\r
- }\r
- else /* prepend */\r
- {\r
- geoms[0] = lwgeom_clone(what);\r
- geoms[1] = lwgeom_clone((LWGEOM *)to);\r
- }\r
- \r
- /* reset SRID and wantbbox flag from component types */\r
- geoms[0]->SRID = geoms[1]->SRID = -1;\r
- TYPE_SETHASSRID(geoms[0]->type, 0);\r
- TYPE_SETHASSRID(geoms[1]->type, 0);\r
- TYPE_SETHASBBOX(geoms[0]->type, 0);\r
- TYPE_SETHASBBOX(geoms[1]->type, 0);\r
-\r
- /* Find appropriate geom type */\r
- if(TYPE_GETTYPE(what->type) == CURVETYPE || TYPE_GETTYPE(what->type) == LINETYPE) newtype = MULTICURVETYPE;\r
- else newtype = COLLECTIONTYPE;\r
-\r
- col = lwcollection_construct(newtype, \r
- to->SRID, NULL, \r
- 2, geoms);\r
-\r
- return (LWGEOM *)col;\r
-}\r
-\r
-void lwcurve_reverse(LWCURVE *curve)\r
-{\r
- ptarray_reverse(curve->points);\r
-}\r
-\r
-/*\r
- * TODO: Invalid segmentization. This should accomodate the curvature.\r
- */\r
-LWCURVE *\r
-lwcurve_segmentize2d(LWCURVE *curve, double dist)\r
-{\r
- return lwcurve_construct(curve->SRID, NULL,\r
- ptarray_segmentize2d(curve->points, dist));\r
-}\r
- \r
-/* check coordinate equality */\r
-char\r
-lwcurve_same(const LWCURVE *me, const LWCURVE *you)\r
-{\r
- return ptarray_same(me->points, you->points);\r
-}\r
-\r
-/*\r
- * Construct a LWCURVE from an array of LWPOINTs\r
- * LWCURVE dimensions are large enough to host all input dimensions.\r
- */\r
-LWCURVE *\r
-lwcurve_from_lwpointarray(int SRID, unsigned int npoints, LWPOINT **points)\r
-{\r
- int zmflag=0;\r
- unsigned int i;\r
- POINTARRAY *pa;\r
- uchar *newpoints, *ptr;\r
- size_t ptsize, size;\r
-\r
- /*\r
- * Find output dimensions, check integrity\r
- */\r
- for(i = 0; i < npoints; i++)\r
- {\r
- if(TYPE_GETTYPE(points[i]->type) != POINTTYPE)\r
- {\r
- lwerror("lwcurve_from_lwpointarray: invalid input type: %s",\r
- lwgeom_typename(TYPE_GETTYPE(points[i]->type)));\r
- return NULL;\r
- }\r
- if(TYPE_HASZ(points[i]->type)) zmflag |= 2;\r
- if(TYPE_HASM(points[i]->type)) zmflag |=1;\r
- if(zmflag == 3) break;\r
- }\r
-\r
- if(zmflag == 0) ptsize = 2 * sizeof(double);\r
- else if(zmflag == 3) ptsize = 4 * sizeof(double);\r
- else ptsize = 3 * sizeof(double);\r
-\r
- /*\r
- * Allocate output points array\r
- */\r
- size = ptsize * npoints;\r
- newpoints = lwalloc(size);\r
- memset(newpoints, 0, size);\r
-\r
- ptr = newpoints;\r
- for(i = 0; i < npoints; i++)\r
- {\r
- size = pointArray_ptsize(points[i]->point);\r
- memcpy(ptr, getPoint_internal(points[i]->point, 0), size);\r
- ptr += ptsize;\r
- }\r
- pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, npoints);\r
-\r
- return lwcurve_construct(SRID, NULL, pa);\r
-}\r
-\r
-/*\r
- * Construct a LWCURVE from a LWMPOINT\r
- */\r
-LWCURVE *\r
-lwcurve_from_lwmpoint(int SRID, LWMPOINT *mpoint)\r
-{\r
- unsigned int i;\r
- POINTARRAY *pa;\r
- char zmflag = TYPE_GETZM(mpoint->type);\r
- size_t ptsize, size;\r
- uchar *newpoints, *ptr;\r
-\r
- if(zmflag == 0) ptsize = 2 * sizeof(double);\r
- else if(zmflag == 3) ptsize = 4 * sizeof(double);\r
- else ptsize = 3 * sizeof(double);\r
-\r
- /* Allocate space for output points */\r
- size = ptsize * mpoint->ngeoms;\r
- newpoints = lwalloc(size);\r
- memset(newpoints, 0, size);\r
-\r
- ptr = newpoints;\r
- for(i = 0; i < mpoint->ngeoms; i++)\r
- {\r
- memcpy(ptr,\r
- getPoint_internal(mpoint->geoms[i]->point, 0),\r
- ptsize);\r
- ptr += ptsize;\r
- }\r
-\r
- pa = pointArray_construct(newpoints, zmflag&2, zmflag&1,\r
- mpoint->ngeoms);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_from_lwmpoint: constructed pointarray for %d points, %d zmflag", mpoint->ngeoms, zmflag);\r
-#endif\r
- \r
- return lwcurve_construct(SRID, NULL, pa);\r
-}\r
-\r
-LWCURVE *\r
-lwcurve_addpoint(LWCURVE *curve, LWPOINT *point, unsigned int where)\r
-{\r
- POINTARRAY *newpa;\r
- LWCURVE *ret;\r
-\r
- newpa = ptarray_addPoint(curve->points, \r
- getPoint_internal(point->point, 0),\r
- TYPE_NDIMS(point->type), where);\r
- ret = lwcurve_construct(curve->SRID, NULL, newpa);\r
-\r
- return ret;\r
-}\r
-\r
-LWCURVE *\r
-lwcurve_removepoint(LWCURVE *curve, unsigned int index)\r
-{\r
- POINTARRAY *newpa;\r
- LWCURVE *ret;\r
-\r
- newpa = ptarray_removePoint(curve->points, index);\r
- ret = lwcurve_construct(curve->SRID, NULL, newpa);\r
-\r
- return ret;\r
-}\r
-\r
-/*\r
- * Note: input will be changed, make sure you have permissions for this.\r
- * */\r
-void\r
-lwcurve_setPoint4d(LWCURVE *curve, unsigned int index, POINT4D *newpoint)\r
-{\r
- setPoint4d(curve->points, index, newpoint);\r
-}\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
-\r
+/**********************************************************************
+ * $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 <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#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);
+}
+
+
+
+
+
+
+
+
+
+
+
+
#include "liblwgeom.h"
+void box_to_box2df(BOX *box, BOX2DFLOAT4 *out);
+
+
/* convert postgresql BOX to BOX2D */
void
/* 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 */
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);
/*
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)
{
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
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 */
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);
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);
/*------------------------------------------------------------------*/
-/**********************************************************************\r
- * $Id$\r
- *\r
- * PostGIS - Spatial Types for PostgreSQL\r
- * http://postgis.refractions.net\r
- * Copyright 2001-2005 Refractions Research Inc.\r
- *\r
- * This is free software; you can redistribute and/or modify it under\r
- * the terms of the GNU General Public Licence. See the COPYING file.\r
- * \r
- **********************************************************************/\r
-\r
-#include "lwgeom_pg.h"\r
-#include "liblwgeom.h"\r
-#include "lwgeom_rtree.h"\r
-\r
-/*\r
- * Creates an rtree given a pointer to the point array.\r
- * Must copy the point array.\r
- */\r
-RTREE_NODE *createTree(POINTARRAY *pointArray)\r
-{\r
- RTREE_NODE *root;\r
- RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints);\r
- int i, nodeCount;\r
- int childNodes, parentNodes;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createTree called with pointarray %p", pointArray);\r
-#endif\r
-\r
- nodeCount = pointArray->npoints - 1;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Total leaf nodes: %d", nodeCount);\r
-#endif\r
-\r
- /*\r
- * Create a leaf node for every line segment.\r
- */\r
- for(i = 0; i < nodeCount; i++)\r
- {\r
- nodes[i] = createLeafNode(pointArray, i);\r
- }\r
-\r
- /*\r
- * Next we group nodes by pairs. If there's an odd number of nodes,\r
- * we bring the last node up a level as is. Continue until we have\r
- * a single top node.\r
- */\r
- childNodes = nodeCount;\r
- parentNodes = nodeCount / 2;\r
- while(parentNodes > 0)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Merging %d children into %d parents.", childNodes, parentNodes);\r
-#endif\r
- i = 0;\r
- while(i < parentNodes) \r
- {\r
- nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]);\r
- i++;\r
- }\r
- /*\r
- * Check for an odd numbered final node.\r
- */\r
- if(parentNodes * 2 < childNodes)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Shuffling child %d to parent %d", childNodes - 1, i);\r
-#endif\r
- nodes[i] = nodes[childNodes - 1];\r
- parentNodes++;\r
- }\r
- childNodes = parentNodes;\r
- parentNodes = parentNodes / 2;\r
- }\r
-\r
- root = nodes[0];\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("createTree returning %p", root);\r
-#endif\r
- return root;\r
-}\r
-\r
-/* \r
- * Creates an interior node given the children. \r
- */\r
-RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){\r
- RTREE_NODE *parent;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createInteriorNode called for children %p, %p", left, right);\r
-#endif\r
- parent = lwalloc(sizeof(RTREE_NODE));\r
- parent->leftNode = left;\r
- parent->rightNode = right;\r
- parent->interval = mergeIntervals(left->interval, right->interval);\r
- parent->segment = NULL;\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createInteriorNode returning %p", parent);\r
-#endif\r
- return parent;\r
-}\r
-\r
-/*\r
- * Creates a leaf node given the pointer to the start point of the segment.\r
- */\r
-RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint)\r
-{\r
- RTREE_NODE *parent;\r
- LWLINE *line;\r
- double value1;\r
- double value2;\r
- POINT4D tmp;\r
- POINTARRAY *npa;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createLeafNode called for point %d of %p", startPoint, pa);\r
-#endif\r
- if(pa->npoints < startPoint + 2)\r
- {\r
- lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);\r
- }\r
-\r
- /*\r
- * The given point array will be part of a geometry that will be freed\r
- * independently of the index. Since we may want to cache the index,\r
- * we must create independent arrays.\r
- */ \r
- npa = lwalloc(sizeof(POINTARRAY));\r
- npa->dims = 0;\r
- npa->npoints = 2;\r
- TYPE_SETZM(npa->dims, 0, 0);\r
- npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2);\r
-\r
- getPoint4d_p(pa, startPoint, &tmp);\r
-\r
- setPoint4d(npa, 0, &tmp);\r
- value1 = tmp.y;\r
-\r
- getPoint4d_p(pa, startPoint + 1, &tmp);\r
-\r
- setPoint4d(npa, 1, &tmp);\r
- value2 = tmp.y;\r
- \r
- line = lwline_construct(-1, NULL, npa);\r
- \r
- parent = lwalloc(sizeof(RTREE_NODE));\r
- parent->interval = createInterval(value1, value2);\r
- parent->segment = line;\r
- parent->leftNode = NULL;\r
- parent->rightNode = NULL;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createLeafNode returning %p", parent);\r
-#endif\r
- return parent;\r
-}\r
-\r
-/*\r
- * Creates an interval with the total extents of the two given intervals.\r
- */\r
-INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2)\r
-{\r
- INTERVAL *interval;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("mergeIntervals called with %p, %p", inter1, inter2);\r
-#endif\r
- interval = lwalloc(sizeof(INTERVAL));\r
- interval->max = FP_MAX(inter1->max, inter2->max);\r
- interval->min = FP_MIN(inter1->min, inter2->min);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);\r
-#endif\r
- return interval;\r
-}\r
-\r
-/*\r
- * Creates an interval given the min and max values, in arbitrary order.\r
- */\r
-INTERVAL *createInterval(double value1, double value2)\r
-{\r
- INTERVAL *interval;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createInterval called with %8.3f, %8.3f", value1, value2);\r
-#endif\r
- interval = lwalloc(sizeof(INTERVAL));\r
- interval->max = FP_MAX(value1, value2);\r
- interval->min = FP_MIN(value1, value2);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);\r
-#endif\r
- return interval;\r
-}\r
-\r
-/*\r
- * Recursively frees the child nodes, the interval and the line before \r
- * freeing the root node.\r
- */\r
-void freeTree(RTREE_NODE *root)\r
-{\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("freeTree called for %p", root);\r
-#endif\r
- if(root->leftNode)\r
- freeTree(root->leftNode);\r
- if(root->rightNode)\r
- freeTree(root->rightNode);\r
- lwfree(root->interval);\r
- if(root->segment)\r
- lwfree(root->segment);\r
- lwfree(root);\r
-}\r
-\r
-/*\r
- * Retrieves a collection of line segments given the root and crossing value.\r
- * The collection is a multilinestring consisting of two point lines \r
- * representing the segments of the ring that may be crossed by the \r
- * horizontal projection line at the given y value.\r
- */\r
-LWMLINE *findLineSegments(RTREE_NODE *root, double value)\r
-{\r
- LWMLINE *tmp, *result;\r
- LWGEOM **lwgeoms;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("findLineSegments called for tree %p and value %8.3f", root, value);\r
-#endif\r
- result = NULL;\r
-\r
- if(!isContained(root->interval, value))\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("findLineSegments %p: not contained.", root);\r
-#endif\r
- return NULL;\r
- }\r
-\r
- /* If there is a segment defined for this node, include it. */\r
- if(root->segment)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type));\r
-#endif\r
- lwgeoms = lwalloc(sizeof(LWGEOM *));\r
- lwgeoms[0] = (LWGEOM *)root->segment;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type));\r
-#endif\r
- result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms);\r
- }\r
-\r
- /* If there is a left child node, recursively include its results. */\r
- if(root->leftNode)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("findLineSegments %p: recursing left.", root);\r
-#endif\r
- tmp = findLineSegments(root->leftNode, value);\r
- if(tmp)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));\r
-#endif\r
- if(result)\r
- result = mergeMultiLines(result, tmp);\r
- else\r
- result = tmp;\r
- }\r
- }\r
-\r
- /* Same for any right child. */\r
- if(root->rightNode)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("findLineSegments %p: recursing right.", root);\r
-#endif\r
- tmp = findLineSegments(root->rightNode, value);\r
- if(tmp)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));\r
-#endif\r
- if(result)\r
- result = mergeMultiLines(result, tmp);\r
- else\r
- result = tmp;\r
- }\r
- }\r
- \r
- return result;\r
-}\r
-\r
-/* Merges two multilinestrings into a single multilinestring. */\r
-LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2)\r
-{\r
- LWGEOM **geoms;\r
- LWCOLLECTION *col;\r
- int i, j, ngeoms;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type));\r
-#endif\r
- ngeoms = line1->ngeoms + line2->ngeoms;\r
- geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);\r
-\r
- j = 0;\r
- for(i = 0; i < line1->ngeoms; i++, j++)\r
- {\r
- geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);\r
- }\r
- for(i = 0; i < line2->ngeoms; i++, j++)\r
- {\r
- geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);\r
- }\r
- col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms);\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type));\r
-#endif\r
-\r
- return (LWMLINE *)col;\r
-}\r
-\r
-/*\r
- * Returns 1 if min < value <= max, 0 otherwise. */\r
-uint32 isContained(INTERVAL *interval, double value)\r
-{\r
- return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;\r
-}\r
-\r
-PG_FUNCTION_INFO_V1(LWGEOM_polygon_index);\r
-Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS)\r
-{\r
- PG_LWGEOM *igeom, *result;\r
- LWGEOM *geom;\r
- LWPOLY *poly;\r
- LWMLINE *mline;\r
- RTREE_NODE *root;\r
- double yval;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- int i;\r
- lwnotice("polygon_index called.");\r
-#endif\r
- result = NULL;\r
- igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));\r
- yval = PG_GETARG_FLOAT8(1);\r
- geom = lwgeom_deserialize(SERIALIZED_FORM(igeom));\r
- if(TYPE_GETTYPE(geom->type) != POLYGONTYPE)\r
- {\r
- lwgeom_release(geom);\r
- PG_FREE_IF_COPY(igeom, 0);\r
- PG_RETURN_NULL();\r
- }\r
- poly = (LWPOLY *)geom;\r
- root = createTree(poly->rings[0]);\r
-\r
- mline = findLineSegments(root, yval);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("mline returned %p %d", mline, TYPE_GETTYPE(mline->type));\r
- for(i = 0; i < mline->ngeoms; i++)\r
- {\r
- lwnotice("geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type));\r
- }\r
-#endif\r
-\r
- if(mline)\r
- result = pglwgeom_serialize((LWGEOM *)mline);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("returning result %p", result);\r
-#endif\r
- lwfree(root);\r
-\r
- PG_FREE_IF_COPY(igeom, 0);\r
- lwgeom_release((LWGEOM *)poly);\r
- lwgeom_release((LWGEOM *)mline);\r
- PG_RETURN_POINTER(result);\r
- \r
-}\r
-\r
-RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly)\r
-{\r
- RTREE_POLY_CACHE *result;\r
- int i, length;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("createNewCache called with %p", poly);\r
-#endif\r
- result = lwalloc(sizeof(RTREE_POLY_CACHE));\r
- result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);\r
- result->ringCount = poly->nrings;\r
- length = lwgeom_size_poly(serializedPoly);\r
- result->poly = lwalloc(length);\r
- memcpy(result->poly, serializedPoly, length); \r
- for(i = 0; i < result->ringCount; i++)\r
- {\r
- result->ringIndices[i] = createTree(poly->rings[i]);\r
- }\r
-#ifdef PGIS_DEBUG\r
- lwnotice("createNewCache returning %p", result);\r
-#endif\r
- return result;\r
-}\r
-\r
-/* \r
- * Creates a new cachable index if needed, or returns the current cache if\r
- * it is applicable to the current polygon.\r
- * The memory context must be changed to function scope before calling this\r
- * method. The method will allocate memory for the cache it creates,\r
- * as well as freeing the memory of any cache that is no longer applicable.\r
- */\r
-RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly, \r
- RTREE_POLY_CACHE *currentCache)\r
-{\r
- int i, length;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache);\r
-#endif\r
- if(!currentCache)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("No existing cache, create one.");\r
-#endif\r
- return createNewCache(poly, serializedPoly);\r
- }\r
- if(!(currentCache->poly))\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Cache contains no polygon, creating new cache.");\r
-#endif\r
- return createNewCache(poly, serializedPoly);\r
- }\r
-\r
- length = lwgeom_size_poly(serializedPoly);\r
-\r
- if(lwgeom_size_poly(currentCache->poly) != length)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Polygon size mismatch, creating new cache.");\r
-#endif\r
- lwfree(currentCache->poly);\r
- lwfree(currentCache);\r
- return createNewCache(poly, serializedPoly);\r
- }\r
- for(i = 0; i < length; i++) \r
- {\r
- uchar a = serializedPoly[i];\r
- uchar b = currentCache->poly[i];\r
- if(a != b) \r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b);\r
-#endif\r
- lwfree(currentCache->poly);\r
- lwfree(currentCache);\r
- return createNewCache(poly, serializedPoly);\r
- }\r
- }\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Polygon match, retaining current cache, %p.", currentCache);\r
-#endif\r
- return currentCache;\r
-}\r
-\r
+/**********************************************************************
+ * $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;
+}
+
-/**********************************************************************\r
- * $Id$\r
- *\r
- * PostGIS - Spatial Types for PostgreSQL\r
- * http://postgis.refractions.net\r
- * Copyright 2001-2006 Refractions Research Inc.\r
- *\r
- * This is free software; you can redistribute and/or modify it under\r
- * the terms of the GNU General Public Licence. See the COPYING file.\r
- * \r
- **********************************************************************/\r
-\r
-#include <stdio.h>\r
-#include <stdlib.h>\r
-#include <stdarg.h>\r
-#include <string.h>\r
-#include <math.h>\r
-\r
-#include "postgres.h"\r
-#include "liblwgeom.h"\r
-#include "fmgr.h"\r
-#include "wktparse.h"\r
-#include "lwgeom_pg.h"\r
-\r
-/*\r
- * Tolerance used to determine equality.\r
- */\r
-#define EPSILON_SQLMM 1e-8\r
-\r
-/*\r
- * Determines (recursively in the case of collections) whether the geometry\r
- * contains at least on arc geometry or segment.\r
- */\r
-uint32\r
-has_arc(LWGEOM *geom)\r
-{\r
- LWCOLLECTION *col;\r
- int i;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("has_arc called.");\r
-#endif\r
-\r
- switch(lwgeom_getType(geom->type)) \r
- {\r
- case POINTTYPE:\r
- case LINETYPE:\r
- case POLYGONTYPE:\r
- case MULTIPOINTTYPE:\r
- case MULTILINETYPE:\r
- case MULTIPOLYGONTYPE:\r
- return 0;\r
- case CURVETYPE:\r
- return 1;\r
- /* It's a collection that MAY contain an arc */\r
- default:\r
- col = (LWCOLLECTION *)geom;\r
- for(i=0; i<col->ngeoms; i++)\r
- {\r
- if(has_arc(col->geoms[i]) == 1) return 1;\r
- }\r
- return 0;\r
- }\r
-}\r
-\r
-/*\r
- * Determines the center of the circle defined by the three given points.\r
- * In the event the circle is complete, the midpoint of the segment defined\r
- * by the first and second points is returned. If the points are colinear,\r
- * as determined by equal slopes, then NULL is returned. If the interior\r
- * point is coincident with either end point, they are taken as colinear.\r
- */\r
-double\r
-lwcircle_center(POINT4D *p1, POINT4D *p2, POINT4D *p3, POINT4D **result)\r
-{\r
- POINT4D *c;\r
- double cx, cy, cr;\r
- double temp, bc, cd, det;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcircle_center called (%.16f, %.16f), (%.16f, %.16f), (%.16f, %.16f).", p1->x, p1->y, p2->x, p2->y, p3->x, p3->y);\r
-#endif\r
-\r
- /* Closed circle */\r
- if(fabs(p1->x - p3->x) < EPSILON_SQLMM\r
- && fabs(p1->y - p3->y) < EPSILON_SQLMM)\r
- {\r
- cx = p1->x + (p2->x - p1->x) / 2.0;\r
- cy = p1->y + (p2->y - p1->y) / 2.0;\r
- c = lwalloc(sizeof(POINT2D));\r
- c->x = cx;\r
- c->y = cy;\r
- *result = c;\r
- cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));\r
- return cr;\r
- }\r
-\r
- temp = p2->x*p2->x + p2->y*p2->y;\r
- bc = (p1->x*p1->x + p1->y*p1->y - temp) / 2.0;\r
- cd = (temp - p3->x*p3->x - p3->y*p3->y) / 2.0;\r
- det = (p1->x - p2->x)*(p2->y - p3->y)-(p2->x - p3->x)*(p1->y - p2->y);\r
-\r
- /* Check colinearity */\r
- if(fabs(det) < EPSILON_SQLMM)\r
- {\r
- *result = NULL;\r
- return -1.0;\r
- }\r
-\r
- det = 1.0 / det;\r
- cx = (bc*(p2->y - p3->y)-cd*(p1->y - p2->y))*det;\r
- cy = ((p1->x - p2->x)*cd-(p2->x - p3->x)*bc)*det;\r
- c = lwalloc(sizeof(POINT4D));\r
- c->x = cx;\r
- c->y = cy;\r
- *result = c;\r
- cr = sqrt((cx-p1->x)*(cx-p1->x)+(cy-p1->y)*(cy-p1->y));\r
- return cr;\r
-}\r
-\r
-double\r
-interpolate_arc(double angle, double zm1, double a1, double zm2, double a2)\r
-{\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("interpolate_arc called.");\r
-#endif\r
-\r
- double frac = fabs((angle - a1) / (a2 - a1));\r
- double result = frac * (zm2 - zm1) + zm1;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("interpolate_arc: angle=%.16f, a1=%.16f, a2=%.16f, z1=%.16f, z2=%.16f, frac=%.16f, result=%.16f", angle, a1, a2, zm1, zm2, frac, result);\r
-#endif\r
-\r
- return result;\r
-}\r
-\r
-/*******************************************************************************\r
- * Begin curve segmentize functions\r
- ******************************************************************************/\r
-\r
-POINTARRAY *\r
-lwcircle_segmentize(POINT4D *p1, POINT4D *p2, POINT4D *p3, uint32 perQuad)\r
-{\r
- POINTARRAY *result;\r
- POINT4D pbuf;\r
- size_t ptsize = sizeof(POINT4D);\r
- unsigned int ptcount;\r
- uchar *pt;\r
-\r
- POINT4D *center;\r
- double radius = 0.0, \r
- sweep = 0.0,\r
- angle = 0.0,\r
- increment = 0.0;\r
- double a1, a2, a3, i;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcircle_segmentize called. ");\r
-#endif\r
-\r
- radius = lwcircle_center(p1, p2, p3, ¢er);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcircle_segmentize, (%.16f, %.16f) radius=%.16f", center->x, center->y, radius);\r
-#endif\r
- if(radius < 0)\r
- {\r
- return NULL;\r
- }\r
-\r
- a1 = atan2(p1->y - center->y, p1->x - center->x);\r
- a2 = atan2(p2->y - center->y, p2->x - center->x);\r
- a3 = atan2(p3->y - center->y, p3->x - center->x);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("a1 = %.16f, a2 = %.16f, a3 = %.16f", a1, a2, a3);\r
-#endif\r
-\r
- if(fabs(p1->x - p3->x) < EPSILON_SQLMM\r
- && fabs(p1->y - p3->y) < EPSILON_SQLMM)\r
- {\r
- sweep = 2*M_PI;\r
- }\r
- /* Clockwise */\r
- else if(a1 > a2 && a2 > a3) \r
- {\r
- sweep = a3 - a1;\r
- }\r
- /* Counter-clockwise */\r
- else if(a1 < a2 && a2 < a3)\r
- {\r
- sweep = a3 - a1;\r
- }\r
- /* Clockwise, wrap */\r
- else if((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3))\r
- {\r
- sweep = a3 - a1 + 2*M_PI;\r
- }\r
- /* Counter-clockwise, wrap */\r
- else if((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3))\r
- {\r
- sweep = a3 - a1 - 2*M_PI;\r
- } \r
- else \r
- {\r
- sweep = 0.0;\r
- }\r
- \r
- ptcount = ceil(fabs(perQuad * sweep / M_PI_2));\r
- \r
- result = ptarray_construct(1, 1, ptcount);\r
-\r
- increment = M_PI_2 / perQuad;\r
- if(sweep < 0) increment *= -1.0;\r
- angle = a1;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("ptcount: %d, perQuad: %d, sweep: %.16f, increment: %.16f", ptcount, perQuad, sweep, increment);\r
-#endif\r
-\r
- for(i = 0; i < ptcount - 1; i++)\r
- {\r
- pt = getPoint_internal(result, i);\r
- angle += increment;\r
- if(increment > 0.0 && angle > M_PI) angle -= 2*M_PI;\r
- if(increment < 0.0 && angle < -1*M_PI) angle -= 2*M_PI;\r
- pbuf.x = center->x + radius*cos(angle);\r
- pbuf.y = center->y + radius*sin(angle);\r
- if((sweep > 0 && angle < a2) || (sweep < 0 && angle > a2))\r
- {\r
- if((sweep > 0 && a1 < a2) || (sweep < 0 && a1 > a2))\r
- {\r
- pbuf.z = interpolate_arc(angle, p1->z, a1, p2->z, a2);\r
- pbuf.m = interpolate_arc(angle, p1->m, a1, p2->m, a2);\r
- }\r
- else\r
- {\r
- if(sweep > 0)\r
- {\r
- pbuf.z = interpolate_arc(angle, p1->z, a1-(2*M_PI), p2->z, a2);\r
- pbuf.m = interpolate_arc(angle, p1->m, a1-(2*M_PI), p2->m, a2);\r
- }\r
- else\r
- {\r
- pbuf.z = interpolate_arc(angle, p1->z, a1+(2*M_PI), p2->z, a2);\r
- pbuf.m = interpolate_arc(angle, p1->m, a1+(2*M_PI), p2->m, a2);\r
- }\r
- }\r
- }\r
- else\r
- {\r
- if((sweep > 0 && a2 < a3) || (sweep < 0 && a2 > a3))\r
- {\r
- pbuf.z = interpolate_arc(angle, p2->z, a2, p3->z, a3);\r
- pbuf.m = interpolate_arc(angle, p2->m, a2, p3->m, a3);\r
- }\r
- else\r
- {\r
- if(sweep > 0)\r
- {\r
- pbuf.z = interpolate_arc(angle, p2->z, a2-(2*M_PI), p3->z, a3);\r
- pbuf.m = interpolate_arc(angle, p2->m, a2-(2*M_PI), p3->m, a3);\r
- }\r
- else\r
- {\r
- pbuf.z = interpolate_arc(angle, p2->z, a2+(2*M_PI), p3->z, a3);\r
- pbuf.m = interpolate_arc(angle, p2->m, a2+(2*M_PI), p3->m, a3);\r
- }\r
- }\r
- }\r
- memcpy(pt, (uchar *)&pbuf, ptsize);\r
- }\r
-\r
- pt = getPoint_internal(result, ptcount - 1);\r
- memcpy(pt, (uchar *)p3, ptsize);\r
-\r
- lwfree(center);\r
-\r
- return result;\r
-}\r
-\r
-LWLINE *\r
-lwcurve_segmentize(LWCURVE *icurve, uint32 perQuad)\r
-{\r
- LWLINE *oline;\r
- DYNPTARRAY *ptarray;\r
- POINTARRAY *tmp;\r
- uint32 i, j;\r
- POINT4D *p1 = lwalloc(sizeof(POINT4D));\r
- POINT4D *p2 = lwalloc(sizeof(POINT4D));\r
- POINT4D *p3 = lwalloc(sizeof(POINT4D));\r
- POINT4D *p4 = lwalloc(sizeof(POINT4D));\r
-\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_segmentize called., dim = %d", icurve->points->dims);\r
-#endif\r
-\r
- ptarray = dynptarray_create(icurve->points->npoints, icurve->points->dims);\r
- if(!getPoint4d_p(icurve->points, 0, p4))\r
- {\r
- elog(ERROR, "curve_segmentize: Cannot extract point.");\r
- }\r
- dynptarray_addPoint4d(ptarray, p4, 1);\r
-\r
- for(i = 2; i < icurve->points->npoints; i+=2) \r
- {\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_segmentize: arc ending at point %d", i);\r
-#endif\r
-\r
- getPoint4d_p(icurve->points, i - 2, p1);\r
- getPoint4d_p(icurve->points, i - 1, p2);\r
- getPoint4d_p(icurve->points, i, p3);\r
- tmp = lwcircle_segmentize(p1, p2, p3, perQuad);\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("lwcurve_segmentize: generated %d points", tmp->npoints);\r
-#endif\r
-\r
- for(j = 0; j < tmp->npoints; j++)\r
- {\r
- getPoint4d_p(tmp, j, p4);\r
- dynptarray_addPoint4d(ptarray, p4, 1);\r
- }\r
- lwfree(tmp);\r
- }\r
- oline = lwline_construct(icurve->SRID, NULL, ptarray_clone(ptarray->pa));\r
-\r
- lwfree(p1);\r
- lwfree(p2);\r
- lwfree(p3);\r
- lwfree(p4);\r
- lwfree(ptarray);\r
- return oline;\r
-}\r
-\r
-LWLINE *\r
-lwcompound_segmentize(LWCOMPOUND *icompound, uint32 perQuad)\r
-{\r
- LWLINE *oline;\r
- LWGEOM *geom;\r
- DYNPTARRAY *ptarray = NULL;\r
- LWLINE *tmp = NULL;\r
- uint32 i, j;\r
- POINT4D *p = NULL;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcompound_segmentize called.");\r
-#endif \r
- p = lwalloc(sizeof(POINT4D));\r
-\r
- ptarray = dynptarray_create(2, ((POINTARRAY *)icompound->geoms[0]->data)->dims);\r
-\r
- for(i = 0; i < icompound->ngeoms; i++)\r
- {\r
- geom = icompound->geoms[i];\r
- if(lwgeom_getType(geom->type) == CURVETYPE)\r
- {\r
- tmp = lwcurve_segmentize((LWCURVE *)geom, perQuad);\r
- for(j = 0; j < tmp->points->npoints; j++)\r
- {\r
- getPoint4d_p(tmp->points, j, p);\r
- dynptarray_addPoint4d(ptarray, p, 0);\r
- }\r
- lwfree(tmp);\r
- }\r
- else if(lwgeom_getType(geom->type) == LINETYPE)\r
- {\r
- tmp = (LWLINE *)geom;\r
- for(j = 0; j < tmp->points->npoints; j++)\r
- {\r
- getPoint4d_p(tmp->points, j, p);\r
- dynptarray_addPoint4d(ptarray, p, 0);\r
- }\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry type %d found.", lwgeom_getType(geom->type));\r
- return NULL;\r
- }\r
- }\r
- oline = lwline_construct(icompound->SRID, NULL, ptarray_clone(ptarray->pa));\r
- lwfree(ptarray);\r
- lwfree(p);\r
- return oline;\r
-}\r
-\r
-LWPOLY *\r
-lwcurvepoly_segmentize(LWCURVEPOLY *curvepoly, uint32 perQuad)\r
-{\r
- LWPOLY *ogeom;\r
- LWGEOM *tmp;\r
- LWLINE *line;\r
- POINTARRAY **ptarray;\r
- int i;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcurvepoly_segmentize called.");\r
-#endif\r
-\r
- ptarray = lwalloc(sizeof(POINTARRAY *)*curvepoly->nrings);\r
-\r
- for(i = 0; i < curvepoly->nrings; i++)\r
- {\r
- tmp = curvepoly->rings[i];\r
- if(lwgeom_getType(tmp->type) == CURVETYPE)\r
- {\r
- line = lwcurve_segmentize((LWCURVE *)tmp, perQuad);\r
- ptarray[i] = ptarray_clone(line->points);\r
- lwfree(line);\r
- }\r
- else if(lwgeom_getType(tmp->type) == LINETYPE)\r
- {\r
- line = (LWLINE *)tmp;\r
- ptarray[i] = ptarray_clone(line->points);\r
- }\r
- else\r
- {\r
- lwerror("Invalid ring type found in CurvePoly.");\r
- return NULL;\r
- }\r
- }\r
-\r
- ogeom = lwpoly_construct(curvepoly->SRID, NULL, curvepoly->nrings, ptarray);\r
- return ogeom;\r
-}\r
-\r
-LWMLINE *\r
-lwmcurve_segmentize(LWMCURVE *mcurve, uint32 perQuad)\r
-{\r
- LWMLINE *ogeom;\r
- LWGEOM *tmp;\r
- LWGEOM **lines;\r
- int i;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwmcurve_segmentize called, geoms=%d, dim=%d.", mcurve->ngeoms, TYPE_NDIMS(mcurve->type));\r
-#endif\r
-\r
- lines = lwalloc(sizeof(LWGEOM *)*mcurve->ngeoms);\r
-\r
- for(i = 0; i < mcurve->ngeoms; i++)\r
- {\r
- tmp = mcurve->geoms[i];\r
- if(lwgeom_getType(tmp->type) == CURVETYPE)\r
- {\r
- lines[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad);\r
- }\r
- else if(lwgeom_getType(tmp->type) == LINETYPE)\r
- {\r
- lines[i] = (LWGEOM *)lwline_construct(mcurve->SRID, NULL, ptarray_clone(((LWLINE *)tmp)->points));\r
- }\r
- else\r
- {\r
- lwerror("Unsupported geometry found in MultiCurve.");\r
- return NULL;\r
- }\r
- }\r
-\r
- ogeom = (LWMLINE *)lwcollection_construct(MULTILINETYPE, mcurve->SRID, NULL, mcurve->ngeoms, lines);\r
- return ogeom;\r
-}\r
-\r
-LWMPOLY *\r
-lwmsurface_segmentize(LWMSURFACE *msurface, uint32 perQuad)\r
-{\r
- LWMPOLY *ogeom;\r
- LWGEOM *tmp;\r
- LWPOLY *poly;\r
- LWGEOM **polys;\r
- POINTARRAY **ptarray;\r
- int i, j;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwmsurface_segmentize called.");\r
-#endif\r
-\r
- polys = lwalloc(sizeof(LWGEOM *)*msurface->ngeoms);\r
-\r
- for(i = 0; i < msurface->ngeoms; i++)\r
- {\r
- tmp = msurface->geoms[i];\r
- if(lwgeom_getType(tmp->type) == CURVEPOLYTYPE)\r
- {\r
- polys[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad);\r
- }\r
- else if(lwgeom_getType(tmp->type) == POLYGONTYPE)\r
- {\r
- poly = (LWPOLY *)tmp;\r
- ptarray = lwalloc(sizeof(POINTARRAY *)*poly->nrings);\r
- for(j = 0; j < poly->nrings; j++)\r
- {\r
- ptarray[j] = ptarray_clone(poly->rings[j]);\r
- }\r
- polys[i] = (LWGEOM *)lwpoly_construct(msurface->SRID, NULL, poly->nrings, ptarray);\r
- }\r
- }\r
- ogeom = (LWMPOLY *)lwcollection_construct(MULTIPOLYGONTYPE, msurface->SRID, NULL, msurface->ngeoms, polys);\r
- return ogeom;\r
-}\r
-\r
-LWCOLLECTION *\r
-lwcollection_segmentize(LWCOLLECTION *collection, uint32 perQuad)\r
-{\r
- LWCOLLECTION *ocol;\r
- LWGEOM *tmp;\r
- LWGEOM **geoms;\r
- int i;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwcollection_segmentize called.");\r
-#endif\r
-\r
- if(has_arc((LWGEOM *)collection) == 0)\r
- {\r
- return collection;\r
- }\r
-\r
- geoms = lwalloc(sizeof(LWGEOM *)*collection->ngeoms);\r
-\r
- for(i=0; i<collection->ngeoms; i++)\r
- {\r
- tmp = collection->geoms[i];\r
- switch(lwgeom_getType(tmp->type)) {\r
- case CURVETYPE:\r
- geoms[i] = (LWGEOM *)lwcurve_segmentize((LWCURVE *)tmp, perQuad);\r
- break;\r
- case COMPOUNDTYPE:\r
- geoms[i] = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)tmp, perQuad);\r
- break;\r
- case CURVEPOLYTYPE:\r
- geoms[i] = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)tmp, perQuad);\r
- break;\r
- default:\r
- geoms[i] = lwgeom_clone(tmp);\r
- break;\r
- }\r
- }\r
- ocol = lwcollection_construct(COLLECTIONTYPE, collection->SRID, NULL, collection->ngeoms, geoms);\r
- return ocol;\r
-}\r
-\r
-LWGEOM *\r
-lwgeom_segmentize(LWGEOM *geom, uint32 perQuad)\r
-{\r
- LWGEOM * ogeom = NULL;\r
- switch(lwgeom_getType(geom->type)) {\r
- case CURVETYPE:\r
- ogeom = (LWGEOM *)lwcurve_segmentize((LWCURVE *)geom, perQuad);\r
- break;\r
- case COMPOUNDTYPE:\r
- ogeom = (LWGEOM *)lwcompound_segmentize((LWCOMPOUND *)geom, perQuad);\r
- break;\r
- case CURVEPOLYTYPE:\r
- ogeom = (LWGEOM *)lwcurvepoly_segmentize((LWCURVEPOLY *)geom, perQuad);\r
- break;\r
- case MULTICURVETYPE:\r
- ogeom = (LWGEOM *)lwmcurve_segmentize((LWMCURVE *)geom, perQuad);\r
- break;\r
- case MULTISURFACETYPE:\r
- ogeom = (LWGEOM *)lwmsurface_segmentize((LWMSURFACE *)geom, perQuad);\r
- break;\r
- case COLLECTIONTYPE:\r
- ogeom = (LWGEOM *)lwcollection_segmentize((LWCOLLECTION *)geom, perQuad);\r
- break;\r
- default:\r
- ogeom = lwgeom_clone(geom);\r
- }\r
- return ogeom;\r
-}\r
-\r
-/*******************************************************************************\r
- * End curve segmentize functions\r
- ******************************************************************************/\r
-LWGEOM *\r
-append_segment(LWGEOM *geom, POINTARRAY *pts, int type, int SRID)\r
-{\r
- LWGEOM *result; \r
- int currentType, i;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("append_segment called %p, %p, %d, %d", geom, pts, type, SRID);\r
-#endif\r
-\r
- if(geom == NULL)\r
- {\r
- if(type == LINETYPE)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: line to NULL");\r
-#endif\r
- return (LWGEOM *)lwline_construct(SRID, NULL, pts);\r
- }\r
- else if(type == CURVETYPE)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: curve to NULL %d", pts->npoints);\r
- POINT4D tmp;\r
- for(i=0; i<pts->npoints; i++)\r
- {\r
- getPoint4d_p(pts, i, &tmp);\r
- lwnotice("new point: (%.16f,%.16f)",tmp.x,tmp.y);\r
- }\r
-#endif\r
- return (LWGEOM *)lwcurve_construct(SRID, NULL, pts);\r
- }\r
- else\r
- {\r
- lwerror("Invalid segment type %d.", type);\r
- }\r
- }\r
- \r
- currentType = lwgeom_getType(geom->type);\r
- if(currentType == LINETYPE && type == LINETYPE)\r
- {\r
- POINTARRAY *newPoints;\r
- POINT4D pt;\r
- LWLINE *line = (LWLINE *)geom;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: line to line");\r
-#endif\r
- newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + line->points->npoints - 1);\r
- for(i=0; i<line->points->npoints; i++)\r
- {\r
- getPoint4d_p(pts, i, &pt);\r
- setPoint4d(newPoints, i, &pt);\r
- }\r
- for(i=1; i<pts->npoints; i++)\r
- {\r
- getPoint4d_p(line->points, i, &pt);\r
- setPoint4d(newPoints, i + pts->npoints, &pt);\r
- }\r
- result = (LWGEOM *)lwline_construct(SRID, NULL, newPoints);\r
- lwgeom_release(geom);\r
- return result;\r
- }\r
- else if(currentType == CURVETYPE && type == CURVETYPE)\r
- {\r
- POINTARRAY *newPoints;\r
- POINT4D pt;\r
- LWCURVE *curve = (LWCURVE *)geom;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: curve to curve");\r
-#endif\r
- newPoints = ptarray_construct(TYPE_HASZ(pts->dims), TYPE_HASM(pts->dims), pts->npoints + curve->points->npoints - 1);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("New array length: %d", pts->npoints + curve->points->npoints - 1);\r
-#endif\r
- for(i=0; i<curve->points->npoints; i++)\r
- {\r
- getPoint4d_p(curve->points, i, &pt);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("orig point %d: (%.16f,%.16f)", i, pt.x, pt.y);\r
-#endif\r
- setPoint4d(newPoints, i, &pt);\r
- }\r
- for(i=1; i<pts->npoints;i++)\r
- {\r
- getPoint4d_p(pts, i, &pt);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("new point %d: (%.16f,%.16f)", i + curve->points->npoints - 1, pt.x, pt.y);\r
-#endif\r
- setPoint4d(newPoints, i + curve->points->npoints - 1, &pt);\r
- }\r
- result = (LWGEOM *)lwcurve_construct(SRID, NULL, newPoints);\r
- lwgeom_release(geom);\r
- return result;\r
- }\r
- else if(currentType == CURVETYPE && type == LINETYPE)\r
- {\r
- LWLINE *line;\r
- LWGEOM **geomArray;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: line to curve");\r
-#endif\r
- geomArray = lwalloc(sizeof(LWGEOM *)*2);\r
- geomArray[0] = lwgeom_clone(geom);\r
- \r
- line = lwline_construct(SRID, NULL, pts);\r
- geomArray[1] = lwgeom_clone((LWGEOM *)line);\r
-\r
- result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray);\r
- lwfree((LWGEOM *)line);\r
- lwgeom_release(geom);\r
- return result;\r
- }\r
- else if(currentType == LINETYPE && type == CURVETYPE)\r
- {\r
- LWCURVE *curve;\r
- LWGEOM **geomArray;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: curve to line");\r
-#endif\r
- geomArray = lwalloc(sizeof(LWGEOM *)*2);\r
- geomArray[0] = lwgeom_clone(geom);\r
-\r
- curve = lwcurve_construct(SRID, NULL, pts);\r
- geomArray[1] = lwgeom_clone((LWGEOM *)curve);\r
-\r
- result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, 2, geomArray);\r
- lwfree((LWGEOM *)curve);\r
- lwgeom_release(geom);\r
- return result;\r
- }\r
- else if(currentType == COMPOUNDTYPE)\r
- {\r
- LWGEOM *newGeom;\r
- LWCOMPOUND *compound;\r
- int count;\r
- LWGEOM **geomArray;\r
- \r
- compound = (LWCOMPOUND *)geom;\r
- count = compound->ngeoms + 1;\r
- geomArray = lwalloc(sizeof(LWGEOM *)*count);\r
- for(i=0; i<compound->ngeoms; i++)\r
- {\r
- geomArray[i] = lwgeom_clone(compound->geoms[i]);\r
- }\r
- if(type == LINETYPE)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: line to compound");\r
-#endif\r
- newGeom = (LWGEOM *)lwline_construct(SRID, NULL, pts);\r
- }\r
- else if(type == CURVETYPE)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("append_segment: curve to compound");\r
-#endif\r
- newGeom = (LWGEOM *)lwcurve_construct(SRID, NULL, pts);\r
- }\r
- else\r
- {\r
- lwerror("Invalid segment type %d.", type);\r
- return NULL;\r
- }\r
- geomArray[compound->ngeoms] = lwgeom_clone(newGeom);\r
-\r
- result = (LWGEOM *)lwcollection_construct(COMPOUNDTYPE, SRID, NULL, count, geomArray);\r
- lwfree(newGeom);\r
- lwgeom_release(geom);\r
- return result;\r
- }\r
- lwerror("Invalid state %d-%d", currentType, type);\r
- return NULL;\r
-}\r
-\r
-LWGEOM *\r
-pta_desegmentize(POINTARRAY *points, int type, int SRID)\r
-{\r
- int i, j, commit, isline, count;\r
- double last_angle, last_length;\r
- double dxab, dyab, dxbc, dybc, theta, length;\r
- POINT4D a, b, c, tmp;\r
- POINTARRAY *pts;\r
- LWGEOM *geom = NULL;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("pta_desegmentize called.");\r
-#endif\r
-\r
- getPoint4d_p(points, 0, &a);\r
- getPoint4d_p(points, 1, &b);\r
- getPoint4d_p(points, 2, &c);\r
-\r
- dxab = b.x - a.x;\r
- dyab = b.y - a.y;\r
- dxbc = c.x - b.x;\r
- dybc = c.y - b.y;\r
-\r
- theta = atan2(dyab, dxab);\r
- last_angle = theta - atan2(dybc, dxbc);\r
- last_length = sqrt(dxbc*dxbc+dybc*dybc);\r
- length = sqrt(dxab*dxab+dyab*dyab);\r
- if((last_length - length) < EPSILON_SQLMM) \r
- {\r
- isline = -1;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Starting as unknown.");\r
-#endif\r
- }\r
- else \r
- {\r
- isline = 1; \r
-#ifdef PGIS_DEBUG\r
- lwnotice("Starting as line.");\r
-#endif\r
- }\r
-\r
- commit = 0;\r
- for(i=3; i<points->npoints; i++)\r
- {\r
- getPoint4d_p(points, i-2, &a);\r
- getPoint4d_p(points, i-1, &b);\r
- getPoint4d_p(points, i, &c);\r
- \r
- dxab = b.x - a.x;\r
- dyab = b.y - a.y;\r
- dxbc = c.x - b.x;\r
- dybc = c.y - b.y;\r
-\r
-#ifdef PGIS_DEBUG\r
- lwnotice("(dxab, dyab, dxbc, dybc) (%.16f, %.16f, %.16f, %.16f)", dxab, dyab, dxbc, dybc);\r
-#endif\r
-\r
- theta = atan2(dyab, dxab);\r
- theta = theta - atan2(dybc, dxbc);\r
- length = sqrt(dxbc*dxbc+dybc*dybc);\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Last/current length and angle %.16f/%.16f, %.16f/%.16f", last_angle, theta, last_length, length);\r
-#endif\r
- /* Found a line segment */\r
- if(fabs(length - last_length) > EPSILON_SQLMM || \r
- fabs(theta - last_angle) > EPSILON_SQLMM)\r
- {\r
- last_length = length;\r
- last_angle = theta;\r
- /* We are tracking a line, keep going */\r
- if(isline > 0)\r
- {\r
- }\r
- /* We were tracking a curve, commit it and start line*/\r
- else if(isline == 0)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Building curve, %d - %d", commit, i);\r
-#endif\r
- count = i - commit;\r
- pts = ptarray_construct(\r
- TYPE_HASZ(type),\r
- TYPE_HASM(type),\r
- 3);\r
- getPoint4d_p(points, commit, &tmp);\r
- setPoint4d(pts, 0, &tmp);\r
- getPoint4d_p(points, \r
- commit + count/2, &tmp);\r
- setPoint4d(pts, 1, &tmp);\r
- getPoint4d_p(points, i - 1, &tmp);\r
- setPoint4d(pts, 2, &tmp);\r
-\r
- commit = i-1;\r
- geom = append_segment(geom, pts, CURVETYPE, SRID);\r
- isline = -1;\r
-\r
- /* \r
- * We now need to move ahead one point to \r
- * determine if it's a potential new curve, \r
- * since the last_angle value is corrupt.\r
- */\r
- i++;\r
- getPoint4d_p(points, i-2, &a);\r
- getPoint4d_p(points, i-1, &b);\r
- getPoint4d_p(points, i, &c);\r
-\r
- dxab = b.x - a.x;\r
- dyab = b.y - a.y;\r
- dxbc = c.x - b.x;\r
- dybc = c.y - b.y;\r
-\r
- theta = atan2(dyab, dxab);\r
- last_angle = theta - atan2(dybc, dxbc);\r
- last_length = sqrt(dxbc*dxbc+dybc*dybc);\r
- length = sqrt(dxab*dxab+dyab*dyab);\r
- if((last_length - length) < EPSILON_SQLMM) \r
- {\r
- isline = -1;\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Restarting as unknown.");\r
-#endif\r
- }\r
- else \r
- {\r
- isline = 1; \r
-#ifdef PGIS_DEBUG\r
- lwnotice("Restarting as line.");\r
-#endif\r
- }\r
-\r
-\r
- }\r
- /* We didn't know what we were tracking, now we do. */\r
- else\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("It's a line");\r
-#endif\r
- isline = 1;\r
- }\r
- }\r
- /* Found a curve segment */\r
- else\r
- {\r
- /* We were tracking a curve, commit it and start line */\r
- if(isline > 0)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Building line, %d - %d", commit, i-2);\r
-#endif\r
- count = i - commit - 2;\r
-\r
- pts = ptarray_construct(\r
- TYPE_HASZ(type),\r
- TYPE_HASM(type),\r
- count);\r
- for(j=commit;j<i-2;j++)\r
- {\r
- getPoint4d_p(points, j, &tmp);\r
- setPoint4d(pts, j-commit, &tmp);\r
- }\r
-\r
- commit = i-3;\r
- geom = append_segment(geom, pts, LINETYPE, SRID);\r
- isline = -1;\r
- }\r
- /* We are tracking a curve, keep going */\r
- else if(isline == 0)\r
- {\r
- ;\r
- }\r
- /* We didn't know what we were tracking, now we do */\r
- else\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("It's a curve");\r
-#endif \r
- isline = 0;\r
- }\r
- }\r
- }\r
- count = i - commit;\r
- if(isline == 0 && count > 2)\r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Finishing curve %d,%d.", commit, i);\r
-#endif\r
- pts = ptarray_construct(\r
- TYPE_HASZ(type),\r
- TYPE_HASM(type),\r
- 3);\r
- getPoint4d_p(points, commit, &tmp);\r
- setPoint4d(pts, 0, &tmp);\r
- getPoint4d_p(points, commit + count/2, &tmp);\r
- setPoint4d(pts, 1, &tmp);\r
- getPoint4d_p(points, i - 1, &tmp);\r
- setPoint4d(pts, 2, &tmp);\r
-\r
- geom = append_segment(geom, pts, CURVETYPE, SRID);\r
- }\r
- else \r
- {\r
-#ifdef PGIS_DEBUG\r
- lwnotice("Finishing line %d,%d.", commit, i);\r
-#endif\r
- pts = ptarray_construct(\r
- TYPE_HASZ(type),\r
- TYPE_HASM(type),\r
- count);\r
- for(j=commit;j<i;j++)\r
- {\r
- getPoint4d_p(points, j, &tmp);\r
- setPoint4d(pts, j-commit, &tmp);\r
- }\r
- geom = append_segment(geom, pts, LINETYPE, SRID);\r
- }\r
- return geom;\r
-}\r
-\r
-LWGEOM *\r
-lwline_desegmentize(LWLINE *line)\r
-{\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwline_desegmentize called.");\r
-#endif \r
-\r
- return pta_desegmentize(line->points, line->type, line->SRID);\r
-}\r
-\r
-LWGEOM *\r
-lwpolygon_desegmentize(LWPOLY *poly)\r
-{\r
- LWGEOM **geoms;\r
- int i, hascurve = 0;\r
- \r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwpolygon_desegmentize called.");\r
-#endif\r
-\r
- geoms = lwalloc(sizeof(LWGEOM *)*poly->nrings);\r
- for(i=0; i<poly->nrings; i++)\r
- {\r
- geoms[i] = pta_desegmentize(poly->rings[i], poly->type, poly->SRID);\r
- if(lwgeom_getType(geoms[i]->type) == CURVETYPE ||\r
- lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE)\r
- {\r
- hascurve = 1;\r
- }\r
- }\r
- if(hascurve == 0)\r
- {\r
- for(i=0; i<poly->nrings; i++)\r
- {\r
- lwfree(geoms[i]);\r
- }\r
- return lwgeom_clone((LWGEOM *)poly);\r
- }\r
-\r
- return (LWGEOM *)lwcollection_construct(CURVEPOLYTYPE, poly->SRID, NULL, poly->nrings, geoms);\r
-}\r
-\r
-LWGEOM *\r
-lwmline_desegmentize(LWMLINE *mline)\r
-{\r
- LWGEOM **geoms;\r
- int i, hascurve = 0;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwmline_desegmentize called.");\r
-#endif\r
-\r
- geoms = lwalloc(sizeof(LWGEOM *)*mline->ngeoms);\r
- for(i=0; i<mline->ngeoms; i++)\r
- {\r
- geoms[i] = lwline_desegmentize((LWLINE *)mline->geoms[i]);\r
- if(lwgeom_getType(geoms[i]->type) == CURVETYPE ||\r
- lwgeom_getType(geoms[i]->type) == COMPOUNDTYPE)\r
- {\r
- hascurve = 1;\r
- }\r
- }\r
- if(hascurve == 0)\r
- {\r
- for(i=0; i<mline->ngeoms; i++)\r
- {\r
- lwfree(geoms[i]);\r
- }\r
- return lwgeom_clone((LWGEOM *)mline);\r
- }\r
- return (LWGEOM *)lwcollection_construct(MULTICURVETYPE, mline->SRID, NULL, mline->ngeoms, geoms);\r
-}\r
-\r
-LWGEOM *\r
-lwmpolygon_desegmentize(LWMPOLY *mpoly)\r
-{\r
- LWGEOM **geoms;\r
- int i, hascurve = 0;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwmpoly_desegmentize called.");\r
-#endif\r
-\r
- geoms = lwalloc(sizeof(LWGEOM *)*mpoly->ngeoms);\r
- for(i=0; i<mpoly->ngeoms; i++)\r
- {\r
- geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]);\r
- if(lwgeom_getType(geoms[i]->type) == CURVEPOLYTYPE)\r
- {\r
- hascurve = 1;\r
- }\r
- }\r
- if(hascurve == 0)\r
- {\r
- for(i=0; i<mpoly->ngeoms; i++)\r
- {\r
- lwfree(geoms[i]);\r
- }\r
- return lwgeom_clone((LWGEOM *)mpoly);\r
- }\r
- return (LWGEOM *)lwcollection_construct(MULTISURFACETYPE, mpoly->SRID, NULL, mpoly->ngeoms, geoms);\r
-}\r
-\r
-LWGEOM *\r
-lwgeom_desegmentize(LWGEOM *geom)\r
-{\r
- int type = lwgeom_getType(geom->type);\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("lwgeom_desegmentize called.");\r
-#endif\r
-\r
- switch(type) {\r
- case LINETYPE:\r
- return lwline_desegmentize((LWLINE *)geom);\r
- case POLYGONTYPE:\r
- return lwpolygon_desegmentize((LWPOLY *)geom);\r
- case MULTILINETYPE:\r
- return lwmline_desegmentize((LWMLINE *)geom);\r
- case MULTIPOLYGONTYPE:\r
- return lwmpolygon_desegmentize((LWMPOLY *)geom);\r
- default:\r
- return lwgeom_clone(geom);\r
- }\r
-}\r
-\r
-/*******************************************************************************\r
- * Begin PG_FUNCTIONs\r
- ******************************************************************************/\r
-\r
-PG_FUNCTION_INFO_V1(LWGEOM_has_arc);\r
-Datum LWGEOM_has_arc(PG_FUNCTION_ARGS)\r
-{\r
- PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));\r
- uint32 result = has_arc(lwgeom_deserialize(SERIALIZED_FORM(geom)));\r
- PG_RETURN_BOOL(result == 1);\r
-}\r
-\r
-/*\r
- * Converts any curve segments of the geometry into a linear approximation.\r
- * Curve centers are determined by projecting the defining points into the 2d\r
- * plane. Z and M values are assigned by linear interpolation between \r
- * defining points.\r
- */\r
-PG_FUNCTION_INFO_V1(LWGEOM_curve_segmentize);\r
-Datum LWGEOM_curve_segmentize(PG_FUNCTION_ARGS)\r
-{\r
- PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));\r
- uint32 perQuad = PG_GETARG_INT32(1);\r
- PG_LWGEOM *ret;\r
- LWGEOM *igeom = NULL, *ogeom = NULL;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("LWGEOM_curve_segmentize called.");\r
-#endif\r
-\r
- if(perQuad < 0) \r
- {\r
- elog(ERROR, "2nd argument must be positive.");\r
- PG_RETURN_NULL();\r
- }\r
-#ifdef PGIS_DEBUG\r
- else\r
- {\r
- lwnotice("perQuad = %d", perQuad);\r
- }\r
-#endif\r
- igeom = lwgeom_deserialize(SERIALIZED_FORM(geom));\r
- ogeom = lwgeom_segmentize(igeom, perQuad);\r
- if(ogeom == NULL) PG_RETURN_NULL();\r
- ret = pglwgeom_serialize(ogeom);\r
- lwgeom_release(igeom);\r
- lwgeom_release(ogeom);\r
- PG_FREE_IF_COPY(geom, 0);\r
- PG_RETURN_POINTER(ret);\r
-}\r
-\r
-PG_FUNCTION_INFO_V1(LWGEOM_line_desegmentize);\r
-Datum LWGEOM_line_desegmentize(PG_FUNCTION_ARGS)\r
-{\r
- PG_LWGEOM *geom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));\r
- PG_LWGEOM *ret;\r
- LWGEOM *igeom = NULL, *ogeom = NULL;\r
-\r
-#ifdef PGIS_DEBUG_CALLS\r
- lwnotice("LWGEOM_line_desegmentize.");\r
-#endif\r
-\r
- igeom = lwgeom_deserialize(SERIALIZED_FORM(geom));\r
- ogeom = lwgeom_desegmentize(igeom);\r
- if(ogeom == NULL)\r
- {\r
- lwgeom_release(igeom);\r
- PG_RETURN_NULL();\r
- }\r
- ret = pglwgeom_serialize(ogeom);\r
- lwgeom_release(igeom);\r
- lwgeom_release(ogeom);\r
- PG_FREE_IF_COPY(geom, 0);\r
- PG_RETURN_POINTER(ret);\r
-}\r
-\r
-/*******************************************************************************\r
- * End PG_FUNCTIONs\r
- ******************************************************************************/\r
+/**********************************************************************
+ * $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 <stdio.h>
+#include <stdlib.h>
+#include <stdarg.h>
+#include <string.h>
+#include <math.h>
+
+#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; i<col->ngeoms; 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; i<collection->ngeoms; 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; i<pts->npoints; 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; i<line->points->npoints; i++)
+ {
+ getPoint4d_p(pts, i, &pt);
+ setPoint4d(newPoints, i, &pt);
+ }
+ for(i=1; i<pts->npoints; 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; i<curve->points->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; i<pts->npoints;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; i<compound->ngeoms; 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; i<points->npoints; 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<i-2;j++)
+ {
+ getPoint4d_p(points, j, &tmp);
+ setPoint4d(pts, j-commit, &tmp);
+ }
+
+ commit = i-3;
+ geom = append_segment(geom, pts, LINETYPE, SRID);
+ isline = -1;
+ }
+ /* We are tracking a curve, keep going */
+ else if(isline == 0)
+ {
+ ;
+ }
+ /* We didn't know what we were tracking, now we do */
+ else
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("It's a curve");
+#endif
+ isline = 0;
+ }
+ }
+ }
+ count = i - commit;
+ if(isline == 0 && count > 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;j<i;j++)
+ {
+ getPoint4d_p(points, j, &tmp);
+ setPoint4d(pts, j-commit, &tmp);
+ }
+ geom = append_segment(geom, pts, LINETYPE, SRID);
+ }
+ return geom;
+}
+
+LWGEOM *
+lwline_desegmentize(LWLINE *line)
+{
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("lwline_desegmentize called.");
+#endif
+
+ return pta_desegmentize(line->points, 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; i<poly->nrings; 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; i<poly->nrings; 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; i<mline->ngeoms; 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; i<mline->ngeoms; 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; i<mpoly->ngeoms; i++)
+ {
+ geoms[i] = lwpolygon_desegmentize((LWPOLY *)mpoly->geoms[i]);
+ if(lwgeom_getType(geoms[i]->type) == CURVEPOLYTYPE)
+ {
+ hascurve = 1;
+ }
+ }
+ if(hascurve == 0)
+ {
+ for(i=0; i<mpoly->ngeoms; 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
+ ******************************************************************************/
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
/* Search path for PROJ.4 library */
static bool IsPROJ4LibPathSet = false;
-void SetPROJ4LibPath();
+void SetPROJ4LibPath(void);
/* Memory context cache functions */
static void PROJ4SRSCacheInit(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)
*/
}
{
/*
* 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;
}
{
/*
* 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);
* 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;
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);
projects / postgis / commitdiff