[postgis-commits] svn - r2781 - trunk/lwgeom
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Thu May 22 13:43:00 PDT 2008
Author: mcayland
Date: 2008-05-22 13:43:00 -0700 (Thu, 22 May 2008)
New Revision: 2781
Modified:
trunk/lwgeom/liblwgeom.h
trunk/lwgeom/long_xact.c
trunk/lwgeom/lwcurve.c
trunk/lwgeom/lwgeom_box.c
trunk/lwgeom/lwgeom_chip.c
trunk/lwgeom/lwgeom_functions_analytic.c
trunk/lwgeom/lwgeom_functions_basic.c
trunk/lwgeom/lwgeom_inout.c
trunk/lwgeom/lwgeom_ogc.c
trunk/lwgeom/lwgeom_rtree.c
trunk/lwgeom/lwgeom_sqlmm.c
trunk/lwgeom/lwgeom_transform.c
trunk/lwgeom/wktunparse.c
Log:
Since PGXS compiles libraries with -Wall, attempt to remove as many warnings as possible. Most of these are missing function prototypes at the top of each file.
Modified: trunk/lwgeom/liblwgeom.h
===================================================================
--- trunk/lwgeom/liblwgeom.h 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/liblwgeom.h 2008-05-22 20:43:00 UTC (rev 2781)
@@ -912,6 +912,7 @@
extern double nextDown_d(float d);
extern double nextUp_d(float d);
+extern float nextafterf_custom(float x, float y);
#define LW_MAX(a,b) ((a) > (b) ? (a) : (b))
@@ -1070,6 +1071,7 @@
extern void ptarray_reverse(POINTARRAY *pa);
extern POINTARRAY *ptarray_substring(POINTARRAY *, double, double);
extern double ptarray_locate_point(POINTARRAY *, POINT2D *);
+extern void closest_point_on_segment(POINT2D *p, POINT2D *A, POINT2D *B, POINT2D *ret);
/*
* Ensure every segment is at most 'dist' long.
@@ -1163,6 +1165,10 @@
* 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);
/*
Modified: trunk/lwgeom/long_xact.c
===================================================================
--- trunk/lwgeom/long_xact.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/long_xact.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -8,6 +8,7 @@
#define ABORT_ON_AUTH_FAILURE 1
Datum check_authorization(PG_FUNCTION_ARGS);
+Datum getTransactionID(PG_FUNCTION_ARGS);
/*
* This trigger will check for authorization before
Modified: trunk/lwgeom/lwcurve.c
===================================================================
--- trunk/lwgeom/lwcurve.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwcurve.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -1,750 +1,762 @@
-/**********************************************************************
- * $Id$
- *
- * PostGIS - Spatial Types for PostgreSQL
- * http://postgis.refractions.net
- * Copyright 2001-2006 Refractions Research Inc.
- *
- * This is free software; you can redistribute and/or modify it under
- * the terms of the GNU General Public Licence. See the COPYING file.
- *
- **********************************************************************/
-
-/* basic LWCURVE functions */
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include "liblwgeom.h"
-
-/*#define PGIS_DEBUG_CALLS 1 */
-/*#define PGIS_DEBUG 1 */
-
-#ifndef MAXFLOAT
- #define MAXFLOAT 3.402823466e+38F
-#endif
-
-/*
- * Construct a new LWCURVE. points will *NOT* be copied
- * use SRID=-1 for unknown SRID (will have 8bit type's S = 0)
- */
-LWCURVE *
-lwcurve_construct(int SRID, BOX2DFLOAT4 *bbox, POINTARRAY *points)
-{
- LWCURVE *result;
-
- /*
- * The first arc requires three points. Each additional
- * arc requires two more points. Thus the minimum point count
- * is three, and the count must be odd.
- */
- if(points->npoints % 2 != 1 || points->npoints < 3)
- {
- lwerror("lwcurve_construct: invalid point count %d", points->npoints);
- return NULL;
- }
-
- result = (LWCURVE*) lwalloc(sizeof(LWCURVE));
-
- result->type = lwgeom_makeType_full(
- TYPE_HASZ(points->dims),
- TYPE_HASM(points->dims),
- (SRID!=-1), CURVETYPE, 0);
- result->SRID = SRID;
- result->points = points;
- result->bbox = bbox;
-
- return result;
-}
-
-/*
- * given the LWGEOM serialized form (or a point into a multi* one)
- * construct a propert LWCURVE.
- * serialized_form should point to the 8bit type format (with type = 8)
- * See serialized form doc
- */
-LWCURVE *
-lwcurve_deserialize(uchar *serialized_form)
-{
- uchar type;
- LWCURVE *result;
- uchar *loc=NULL;
- uint32 npoints;
- POINTARRAY *pa;
-
- type = (uchar)serialized_form[0];
- if(lwgeom_getType(type) != CURVETYPE)
- {
- lwerror("lwcurve_deserialize: attempt to deserialize a curve which is really a %s", lwgeom_typename(type));
- return NULL;
- }
-
- result = (LWCURVE*) lwalloc(sizeof(LWCURVE));
- result->type = type;
-
- loc = serialized_form + 1;
-
- if(lwgeom_hasBBOX(type))
- {
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_deserialize: input has bbox");
-#endif
- result->bbox = lwalloc(sizeof(BOX2DFLOAT4));
- memcpy(result->bbox, loc, sizeof(BOX2DFLOAT4));
- loc += sizeof(BOX2DFLOAT4);
- }
- else
- {
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_deserialize: input lacks bbox");
-#endif
- result->bbox = NULL;
- }
-
- if(lwgeom_hasSRID(type))
- {
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_deserialize: input has srid");
-#endif
- result->SRID = lw_get_int32(loc);
- loc += 4; /* type + SRID */
- }
- else
- {
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_deserialize: input lacks srid");
-#endif
- result->SRID = -1;
- }
-
- /* we've read the type (1 byte) and SRID (4 bytes, if present) */
-
- npoints = lw_get_uint32(loc);
-#ifdef PGIS_DEBUG
- lwnotice("curve npoints = %d", npoints);
-#endif
- loc += 4;
- pa = pointArray_construct(loc, TYPE_HASZ(type), TYPE_HASM(type), npoints);
- result->points = pa;
- return result;
-}
-
-/*
- * convert this curve into its serialized form
- * result's first char will be the 8bit type. See serialized form doc
- */
-uchar *
-lwcurve_serialize(LWCURVE *curve)
-{
- size_t size, retsize;
- uchar * result;
-
- if(curve == NULL) {
- lwerror("lwcurve_serialize:: given null curve");
- return NULL;
- }
-
- size = lwcurve_serialize_size(curve);
- result = lwalloc(size);
- lwcurve_serialize_buf(curve, result, &retsize);
- if(retsize != size)
- lwerror("lwcurve_serialize_size returned %d, ..selialize_buf returned %d", size, retsize);
- return result;
-}
-
-/*
- * convert this curve into its serialized form writing it into
- * the given buffer, and returning number of bytes written into
- * the given int pointer.
- * result's first char will be the 8bit type. See serialized form doc
- */
-void lwcurve_serialize_buf(LWCURVE *curve, uchar *buf, size_t *retsize)
-{
- char hasSRID;
- uchar *loc;
- int ptsize;
- size_t size;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_serialize_buf(%p, %p, %p) called",
- curve, buf, retsize);
-#endif
-
- if(curve == NULL)
- {
- lwerror("lwcurve_serialize:: given null curve");
- return;
- }
-
- if(TYPE_GETZM(curve->type) != TYPE_GETZM(curve->points->dims))
- {
- lwerror("Dimensions mismatch in lwcurve");
- return;
- }
-
- ptsize = pointArray_ptsize(curve->points);
-
- hasSRID = (curve->SRID != -1);
-
- buf[0] = (uchar)lwgeom_makeType_full(
- TYPE_HASZ(curve->type), TYPE_HASM(curve->type),
- hasSRID, CURVETYPE, curve->bbox ? 1 : 0);
- loc = buf+1;
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_serialize_buf added type (%d)", curve->type);
-#endif
-
- if(curve->bbox)
- {
- memcpy(loc, curve->bbox, sizeof(BOX2DFLOAT4));
- loc += sizeof(BOX2DFLOAT4);
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_serialize_buf added BBOX");
-#endif
-
- }
-
- if(hasSRID)
- {
- memcpy(loc, &curve->SRID, sizeof(int32));
- loc += sizeof(int32);
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_serialize_buf added SRID");
-#endif
-
- }
-
- memcpy(loc, &curve->points->npoints, sizeof(uint32));
- loc += sizeof(uint32);
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_serialize_buf added npoints (%d)",
- curve->points->npoints);
-#endif
-
- /* copy in points */
- size = curve->points->npoints * ptsize;
- memcpy(loc, getPoint_internal(curve->points, 0), size);
- loc += size;
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_serialize_buf copied serialized_pointlist (%d bytes)",
- ptsize * curve->points->npoints);
-#endif
-
- if(retsize) *retsize = loc-buf;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_serialize_buf returning (loc: %p, size: %d)",
- loc, loc-buf);
-#endif
-}
-
-/* find length of this deserialized curve */
-size_t
-lwcurve_serialize_size(LWCURVE *curve)
-{
- size_t size = 1; /* type */
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_serialize_size called");
-#endif
-
- if(curve->SRID != -1) size += 4; /* SRID */
- if(curve->bbox) size += sizeof(BOX2DFLOAT4);
-
- size += 4; /* npoints */
- size += pointArray_ptsize(curve->points) * curve->points->npoints;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_serialize_size returning %d", size);
-#endif
-
- return size;
-}
-
-BOX3D *
-lwcircle_compute_box3d(POINT4D *p1, POINT4D *p2, POINT4D *p3)
-{
- double x1, x2, y1, y2, z1, z2;
- double angle, radius, sweep;
- /*
- double top, left;
- */
- double a1, a2, a3;
- double xe = 0.0, ye = 0.0;
- POINT4D *center;
- int i;
- BOX3D *box;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcircle_compute_box3d called.");
-#endif
-
- center = lwalloc(sizeof(POINT4D));
- radius = lwcircle_center(p1, p2, p3, ¢er);
- if(radius < 0.0) return NULL;
-
- /*
- top = center->y + radius;
- left = center->x - radius;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcircle_compute_box3d: top=%.16f, left=%.16f", top, left);
-#endif
- */
-
- x1 = MAXFLOAT;
- x2 = -1 * MAXFLOAT;
- y1 = MAXFLOAT;
- y2 = -1 * MAXFLOAT;
-
- a1 = atan2(p1->y - center->y, p1->x - center->x);
- a2 = atan2(p2->y - center->y, p2->x - center->x);
- a3 = atan2(p3->y - center->y, p3->x - center->x);
-
- /* Determine sweep angle */
- if(a1 > a2 && a2 > a3)
- {
- sweep = a3 - a1;
- }
- /* Counter-clockwise */
- else if(a1 < a2 && a2 < a3)
- {
- sweep = a3 - a1;
- }
- /* Clockwise, wrap */
- else if((a1 < a2 && a1 > a3) || (a2 < a3 && a1 > a3))
- {
- sweep = a3 - a1 + 2*M_PI;
- }
- /* Counter-clockwise, wrap */
- else if((a1 > a2 && a1 < a3) || (a2 > a3 && a1 < a3))
- {
- sweep = a3 - a1 - 2*M_PI;
- }
- else
- {
- sweep = 0.0;
- }
-
-#ifdef PGIS_DEBUG
- lwnotice("a1 %.16f, a2 %.16f, a3 %.16f, sweep %.16f", a1, a2, a3, sweep);
-#endif
-
- angle = 0.0;
- for(i=0; i < 6; i++)
- {
- switch(i) {
- case 0:
- angle = 0.0;
- xe = center->x + radius;
- ye = center->y;
- break;
- case 1:
- angle = M_PI_2;
- xe = center->x;
- ye = center->y + radius;
- break;
- case 2:
- angle = M_PI;
- xe = center->x - radius;
- ye = center->y;
- break;
- case 3:
- angle = -1 * M_PI_2;
- xe = center->x;
- ye = center->y - radius;
- break;
- case 4:
- angle = a1;
- xe = p1->x;
- ye = p1->y;
- break;
- case 5:
- angle = a3;
- xe = p3->x;
- ye = p3->y;
- break;
- }
- if(i < 4)
- {
- if(sweep > 0.0 && (angle > a3 || angle < a1)) continue;
- if(sweep < 0.0 && (angle < a3 || angle > a1)) continue;
- }
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcircle_compute_box3d: potential extreame %d (%.16f, %.16f)", i, xe, ye);
-#endif
- x1 = (x1 < xe) ? x1 : xe;
- y1 = (y1 < ye) ? y1 : ye;
- x2 = (x2 > xe) ? x2 : xe;
- y2 = (y2 > ye) ? y2 : ye;
- }
-#ifdef PGIS_DEBUG
- lwnotice("lwcircle_compute_box3d: extreames found (%.16f %.16f, %.16f %.16f)", x1, y1, x2, y2);
-#endif
-
- /*
- x1 = center->x + x1 * radius;
- x2 = center->x + x2 * radius;
- y1 = center->y + y1 * radius;
- y2 = center->y + y2 * radius;
- */
- z1 = (p1->z < p2->z) ? p1->z : p2->z;
- z1 = (z1 < p3->z) ? z1 : p3->z;
- z2 = (p1->z > p2->z) ? p1->z : p2->z;
- z2 = (z2 > p3->z) ? z2 : p3->z;
-
- box = lwalloc(sizeof(BOX3D));
- box->xmin = x1; box->xmax = x2;
- box->ymin = y1; box->ymax = y2;
- box->zmin = z1; box->zmax = z2;
-
- lwfree(center);
-
- return box;
-}
-
-/*
- * Find bounding box (standard one)
- * zmin=zmax=NO_Z_VALUE if 2d
- * TODO: This ignores curvature, which should be taken into account.
- */
-BOX3D *
-lwcurve_compute_box3d(LWCURVE *curve)
-{
- BOX3D *box, *tmp;
- int i;
- POINT4D *p1 = lwalloc(sizeof(POINT4D));
- POINT4D *p2 = lwalloc(sizeof(POINT4D));
- POINT4D *p3 = lwalloc(sizeof(POINT4D));
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_compute_box3d called.");
-#endif
-
- /* initialize box values */
- box = lwalloc(sizeof(BOX3D));
- box->xmin = MAXFLOAT; box->xmax = -1 * MAXFLOAT;
- box->ymin = MAXFLOAT; box->ymax = -1 * MAXFLOAT;
- box->zmin = MAXFLOAT; box->zmax = -1 * MAXFLOAT;
-
- for(i = 2; i < curve->points->npoints; i+=2)
- {
- getPoint4d_p(curve->points, i-2, p1);
- getPoint4d_p(curve->points, i-1, p2);
- getPoint4d_p(curve->points, i, p3);
- tmp = lwcircle_compute_box3d(p1, p2, p3);
- if(tmp == NULL) continue;
- box->xmin = (box->xmin < tmp->xmin) ? box->xmin : tmp->xmin;
- box->xmax = (box->xmax > tmp->xmax) ? box->xmax : tmp->xmax;
- box->ymin = (box->ymin < tmp->ymin) ? box->ymin : tmp->ymin;
- box->ymax = (box->ymax > tmp->ymax) ? box->ymax : tmp->ymax;
- box->zmin = (box->zmin < tmp->zmin) ? box->zmin : tmp->zmin;
- box->zmax = (box->zmax > tmp->zmax) ? box->zmax : tmp->zmax;
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("curve %d x=(%.16f,%.16f) y=(%.16f,%.16f) z=(%.16f,%.16f)", i/2, box->xmin, box->xmax, box->ymin, box->ymax, box->zmin, box->zmax);
-#endif
- }
-
-
- return box;
-}
-
-int
-lwcurve_compute_box2d_p(LWCURVE *curve, BOX2DFLOAT4 *result)
-{
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwcurve_compute_box2d_p called.");
-#endif
-
- BOX3D *box = lwcurve_compute_box3d(curve);
- if(box == NULL) return 0;
- box3d_to_box2df_p(box, result);
- return 1;
-}
-
-void pfree_curve(LWCURVE *curve)
-{
- lwfree(curve->points);
- lwfree(curve);
-}
-
-/* find length of this serialized curve */
-size_t
-lwgeom_size_curve(const uchar *serialized_curve)
-{
- int type = (uchar)serialized_curve[0];
- uint32 result = 1; /* type */
- const uchar *loc;
- uint32 npoints;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwgeom_size_curve called");
-#endif
- if(lwgeom_getType(type) != CURVETYPE)
- lwerror("lwgeom_size_curve::attempt to find the length of a non-curve");
-
- loc = serialized_curve + 1;
- if(lwgeom_hasBBOX(type))
- {
- loc += sizeof(BOX2DFLOAT4);
- result += sizeof(BOX2DFLOAT4);
- }
-
- if(lwgeom_hasSRID(type))
- {
- loc += 4; /* type + SRID */
- result += 4;
- }
-
- /* we've read the type (1 byte) and SRID (4 bytes, if present) */
- npoints = lw_get_uint32(loc);
- result += sizeof(uint32); /* npoints */
-
- result += TYPE_NDIMS(type) * sizeof(double) * npoints;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("lwgeom_size_curve returning %d", result);
-#endif
-
- return result;
-}
-
-void printLWCURVE(LWCURVE *curve)
-{
- lwnotice("LWCURVE {");
- lwnotice(" ndims = %i", (int)TYPE_NDIMS(curve->type));
- lwnotice(" SRID = %i", (int)curve->SRID);
- printPA(curve->points);
- lwnotice("}");
-}
-
-/* Clone LWCURVE object. POINTARRAY is not copied. */
-LWCURVE *
-lwcurve_clone(const LWCURVE *g)
-{
- LWCURVE *ret = lwalloc(sizeof(LWCURVE));
- memcpy(ret, g, sizeof(LWCURVE));
- if(g->bbox) ret->bbox = box2d_clone(g->bbox);
- return ret;
-}
-
-/*
- * Add 'what' to this curve at position 'where'.
- * where=0 == prepend
- * where=-1 == append
- * Returns a MULTICURVE or a GEOMETRYCOLLECTION
- */
-LWGEOM *
-lwcurve_add(const LWCURVE *to, uint32 where, const LWGEOM *what)
-{
- LWCOLLECTION *col;
- LWGEOM **geoms;
- int newtype;
-
- if(where != -1 && where != 0)
- {
- lwerror("lwcurve_add only supports 0 or -1 as second argument %d", where);
- return NULL;
- }
-
- /* dimensions compatibility are checked by caller */
-
- /* Construct geoms array */
- geoms = lwalloc(sizeof(LWGEOM *)*2);
- if(where == -1) /* append */
- {
- geoms[0] = lwgeom_clone((LWGEOM *)to);
- geoms[1] = lwgeom_clone(what);
- }
- else /* prepend */
- {
- geoms[0] = lwgeom_clone(what);
- geoms[1] = lwgeom_clone((LWGEOM *)to);
- }
-
- /* reset SRID and wantbbox flag from component types */
- geoms[0]->SRID = geoms[1]->SRID = -1;
- TYPE_SETHASSRID(geoms[0]->type, 0);
- TYPE_SETHASSRID(geoms[1]->type, 0);
- TYPE_SETHASBBOX(geoms[0]->type, 0);
- TYPE_SETHASBBOX(geoms[1]->type, 0);
-
- /* Find appropriate geom type */
- if(TYPE_GETTYPE(what->type) == CURVETYPE || TYPE_GETTYPE(what->type) == LINETYPE) newtype = MULTICURVETYPE;
- else newtype = COLLECTIONTYPE;
-
- col = lwcollection_construct(newtype,
- to->SRID, NULL,
- 2, geoms);
-
- return (LWGEOM *)col;
-}
-
-void lwcurve_reverse(LWCURVE *curve)
-{
- ptarray_reverse(curve->points);
-}
-
-/*
- * TODO: Invalid segmentization. This should accomodate the curvature.
- */
-LWCURVE *
-lwcurve_segmentize2d(LWCURVE *curve, double dist)
-{
- return lwcurve_construct(curve->SRID, NULL,
- ptarray_segmentize2d(curve->points, dist));
-}
-
-/* check coordinate equality */
-char
-lwcurve_same(const LWCURVE *me, const LWCURVE *you)
-{
- return ptarray_same(me->points, you->points);
-}
-
-/*
- * Construct a LWCURVE from an array of LWPOINTs
- * LWCURVE dimensions are large enough to host all input dimensions.
- */
-LWCURVE *
-lwcurve_from_lwpointarray(int SRID, unsigned int npoints, LWPOINT **points)
-{
- int zmflag=0;
- unsigned int i;
- POINTARRAY *pa;
- uchar *newpoints, *ptr;
- size_t ptsize, size;
-
- /*
- * Find output dimensions, check integrity
- */
- for(i = 0; i < npoints; i++)
- {
- if(TYPE_GETTYPE(points[i]->type) != POINTTYPE)
- {
- lwerror("lwcurve_from_lwpointarray: invalid input type: %s",
- lwgeom_typename(TYPE_GETTYPE(points[i]->type)));
- return NULL;
- }
- if(TYPE_HASZ(points[i]->type)) zmflag |= 2;
- if(TYPE_HASM(points[i]->type)) zmflag |=1;
- if(zmflag == 3) break;
- }
-
- if(zmflag == 0) ptsize = 2 * sizeof(double);
- else if(zmflag == 3) ptsize = 4 * sizeof(double);
- else ptsize = 3 * sizeof(double);
-
- /*
- * Allocate output points array
- */
- size = ptsize * npoints;
- newpoints = lwalloc(size);
- memset(newpoints, 0, size);
-
- ptr = newpoints;
- for(i = 0; i < npoints; i++)
- {
- size = pointArray_ptsize(points[i]->point);
- memcpy(ptr, getPoint_internal(points[i]->point, 0), size);
- ptr += ptsize;
- }
- pa = pointArray_construct(newpoints, zmflag&2, zmflag&1, npoints);
-
- return lwcurve_construct(SRID, NULL, pa);
-}
-
-/*
- * Construct a LWCURVE from a LWMPOINT
- */
-LWCURVE *
-lwcurve_from_lwmpoint(int SRID, LWMPOINT *mpoint)
-{
- unsigned int i;
- POINTARRAY *pa;
- char zmflag = TYPE_GETZM(mpoint->type);
- size_t ptsize, size;
- uchar *newpoints, *ptr;
-
- if(zmflag == 0) ptsize = 2 * sizeof(double);
- else if(zmflag == 3) ptsize = 4 * sizeof(double);
- else ptsize = 3 * sizeof(double);
-
- /* Allocate space for output points */
- size = ptsize * mpoint->ngeoms;
- newpoints = lwalloc(size);
- memset(newpoints, 0, size);
-
- ptr = newpoints;
- for(i = 0; i < mpoint->ngeoms; i++)
- {
- memcpy(ptr,
- getPoint_internal(mpoint->geoms[i]->point, 0),
- ptsize);
- ptr += ptsize;
- }
-
- pa = pointArray_construct(newpoints, zmflag&2, zmflag&1,
- mpoint->ngeoms);
-
-#ifdef PGIS_DEBUG
- lwnotice("lwcurve_from_lwmpoint: constructed pointarray for %d points, %d zmflag", mpoint->ngeoms, zmflag);
-#endif
-
- return lwcurve_construct(SRID, NULL, pa);
-}
-
-LWCURVE *
-lwcurve_addpoint(LWCURVE *curve, LWPOINT *point, unsigned int where)
-{
- POINTARRAY *newpa;
- LWCURVE *ret;
-
- newpa = ptarray_addPoint(curve->points,
- getPoint_internal(point->point, 0),
- TYPE_NDIMS(point->type), where);
- ret = lwcurve_construct(curve->SRID, NULL, newpa);
-
- return ret;
-}
-
-LWCURVE *
-lwcurve_removepoint(LWCURVE *curve, unsigned int index)
-{
- POINTARRAY *newpa;
- LWCURVE *ret;
-
- newpa = ptarray_removePoint(curve->points, index);
- ret = lwcurve_construct(curve->SRID, NULL, newpa);
-
- return ret;
-}
-
-/*
- * Note: input will be changed, make sure you have permissions for this.
- * */
-void
-lwcurve_setPoint4d(LWCURVE *curve, unsigned int index, POINT4D *newpoint)
-{
- setPoint4d(curve->points, index, newpoint);
-}
-
-
-
-
-
-
-
-
-
-
-
-
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2001-2006 Refractions Research Inc.
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+/* basic LWCURVE functions */
+
+#include <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);
+}
+
+
+
+
+
+
+
+
+
+
+
+
Modified: trunk/lwgeom/lwgeom_box.c
===================================================================
--- trunk/lwgeom/lwgeom_box.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_box.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -5,7 +5,10 @@
#include "liblwgeom.h"
+void box_to_box2df(BOX *box, BOX2DFLOAT4 *out);
+
+
/* convert postgresql BOX to BOX2D */
void
box_to_box2df(BOX *box, BOX2DFLOAT4 *out)
Modified: trunk/lwgeom/lwgeom_chip.c
===================================================================
--- trunk/lwgeom/lwgeom_chip.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_chip.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -18,10 +18,47 @@
/* Define this to debug CHIP ops */
/*#define DEBUG_CHIP 1*/
+typedef unsigned short int UINT16;
+typedef float FLOAT32;
+
+typedef struct PIXEL_T {
+ int type; /* 1=float32, 5=int24, 6=int16 */
+ uchar val[4];
+} PIXEL;
+
+typedef struct RGB_T {
+ uchar red;
+ uchar green;
+ uchar blue;
+} RGB;
+
+
/* Internal funcs */
void swap_char(char *a,char *b);
void flip_endian_double(char *d);
void flip_endian_int32(char *i);
+const char* pixelOpName(int op);
+const char* pixelHEX(PIXEL* p);
+UINT16 pixel_readUINT16(PIXEL *p);
+void pixel_writeUINT16(PIXEL *p, UINT16 i);
+PIXEL pixel_readval(char *buf);
+void pixel_writeval(PIXEL *p, char *buf, size_t maxlen);
+void pixel_add_float32(PIXEL *where, PIXEL *what);
+void pixel_add_int24(PIXEL *where, PIXEL *what);
+void pixel_add_int16(PIXEL *where, PIXEL *what);
+void pixel_add(PIXEL *where, PIXEL *what);
+size_t chip_xy_off(CHIP *c, size_t x, size_t y);
+void chip_setPixel(CHIP *c, int x, int y, PIXEL *p);
+PIXEL chip_getPixel(CHIP *c, int x, int y);
+void chip_draw_pixel(CHIP *chip, int x, int y, PIXEL *pixel, int op);
+void chip_draw_segment(CHIP *chip, int x1, int y1, int x2, int y2, PIXEL *pixel, int op);
+void chip_fill(CHIP *chip, PIXEL *pixel, int op);
+CHIP * pgchip_construct(BOX3D *bvol, int SRID, int width, int height, int datatype, PIXEL *initvalue);
+void chip_draw_ptarray(CHIP *chip, POINTARRAY *pa, PIXEL *pixel, int op);
+void chip_draw_lwpoint(CHIP *chip, LWPOINT *lwpoint, PIXEL* pixel, int op);
+void chip_draw_lwline(CHIP *chip, LWLINE *lwline, PIXEL* pixel, int op);
+void chip_draw_lwgeom(CHIP *chip, LWGEOM *lwgeom, PIXEL *pixel, int op);
+char * text_to_cstring(text *t);
/* Prototypes */
@@ -36,6 +73,13 @@
Datum CHIP_getHeight(PG_FUNCTION_ARGS);
Datum CHIP_getWidth(PG_FUNCTION_ARGS);
Datum CHIP_setSRID(PG_FUNCTION_ARGS);
+Datum CHIP_send(PG_FUNCTION_ARGS);
+Datum CHIP_dump(PG_FUNCTION_ARGS);
+Datum CHIP_construct(PG_FUNCTION_ARGS);
+Datum CHIP_getpixel(PG_FUNCTION_ARGS);
+Datum CHIP_setpixel(PG_FUNCTION_ARGS);
+Datum CHIP_draw(PG_FUNCTION_ARGS);
+Datum CHIP_fill(PG_FUNCTION_ARGS);
/*
@@ -355,21 +399,6 @@
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)
{
Modified: trunk/lwgeom/lwgeom_functions_analytic.c
===================================================================
--- trunk/lwgeom/lwgeom_functions_analytic.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_functions_analytic.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -43,7 +43,16 @@
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
* returns distance in an int pointer
Modified: trunk/lwgeom/lwgeom_functions_basic.c
===================================================================
--- trunk/lwgeom/lwgeom_functions_basic.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_functions_basic.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -73,8 +73,15 @@
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);
-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 */
Modified: trunk/lwgeom/lwgeom_inout.c
===================================================================
--- trunk/lwgeom/lwgeom_inout.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_inout.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -38,6 +38,7 @@
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);
Modified: trunk/lwgeom/lwgeom_ogc.c
===================================================================
--- trunk/lwgeom/lwgeom_ogc.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_ogc.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -81,6 +81,8 @@
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);
/*------------------------------------------------------------------*/
Modified: trunk/lwgeom/lwgeom_rtree.c
===================================================================
--- trunk/lwgeom/lwgeom_rtree.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_rtree.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -1,475 +1,480 @@
-/**********************************************************************
- * $Id$
- *
- * PostGIS - Spatial Types for PostgreSQL
- * http://postgis.refractions.net
- * Copyright 2001-2005 Refractions Research Inc.
- *
- * This is free software; you can redistribute and/or modify it under
- * the terms of the GNU General Public Licence. See the COPYING file.
- *
- **********************************************************************/
-
-#include "lwgeom_pg.h"
-#include "liblwgeom.h"
-#include "lwgeom_rtree.h"
-
-/*
- * Creates an rtree given a pointer to the point array.
- * Must copy the point array.
- */
-RTREE_NODE *createTree(POINTARRAY *pointArray)
-{
- RTREE_NODE *root;
- RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints);
- int i, nodeCount;
- int childNodes, parentNodes;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createTree called with pointarray %p", pointArray);
-#endif
-
- nodeCount = pointArray->npoints - 1;
-
-#ifdef PGIS_DEBUG
- lwnotice("Total leaf nodes: %d", nodeCount);
-#endif
-
- /*
- * Create a leaf node for every line segment.
- */
- for(i = 0; i < nodeCount; i++)
- {
- nodes[i] = createLeafNode(pointArray, i);
- }
-
- /*
- * Next we group nodes by pairs. If there's an odd number of nodes,
- * we bring the last node up a level as is. Continue until we have
- * a single top node.
- */
- childNodes = nodeCount;
- parentNodes = nodeCount / 2;
- while(parentNodes > 0)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Merging %d children into %d parents.", childNodes, parentNodes);
-#endif
- i = 0;
- while(i < parentNodes)
- {
- nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]);
- i++;
- }
- /*
- * Check for an odd numbered final node.
- */
- if(parentNodes * 2 < childNodes)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Shuffling child %d to parent %d", childNodes - 1, i);
-#endif
- nodes[i] = nodes[childNodes - 1];
- parentNodes++;
- }
- childNodes = parentNodes;
- parentNodes = parentNodes / 2;
- }
-
- root = nodes[0];
-
-#ifdef PGIS_DEBUG
- lwnotice("createTree returning %p", root);
-#endif
- return root;
-}
-
-/*
- * Creates an interior node given the children.
- */
-RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){
- RTREE_NODE *parent;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createInteriorNode called for children %p, %p", left, right);
-#endif
- parent = lwalloc(sizeof(RTREE_NODE));
- parent->leftNode = left;
- parent->rightNode = right;
- parent->interval = mergeIntervals(left->interval, right->interval);
- parent->segment = NULL;
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createInteriorNode returning %p", parent);
-#endif
- return parent;
-}
-
-/*
- * Creates a leaf node given the pointer to the start point of the segment.
- */
-RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint)
-{
- RTREE_NODE *parent;
- LWLINE *line;
- double value1;
- double value2;
- POINT4D tmp;
- POINTARRAY *npa;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createLeafNode called for point %d of %p", startPoint, pa);
-#endif
- if(pa->npoints < startPoint + 2)
- {
- lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
- }
-
- /*
- * The given point array will be part of a geometry that will be freed
- * independently of the index. Since we may want to cache the index,
- * we must create independent arrays.
- */
- npa = lwalloc(sizeof(POINTARRAY));
- npa->dims = 0;
- npa->npoints = 2;
- TYPE_SETZM(npa->dims, 0, 0);
- npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2);
-
- getPoint4d_p(pa, startPoint, &tmp);
-
- setPoint4d(npa, 0, &tmp);
- value1 = tmp.y;
-
- getPoint4d_p(pa, startPoint + 1, &tmp);
-
- setPoint4d(npa, 1, &tmp);
- value2 = tmp.y;
-
- line = lwline_construct(-1, NULL, npa);
-
- parent = lwalloc(sizeof(RTREE_NODE));
- parent->interval = createInterval(value1, value2);
- parent->segment = line;
- parent->leftNode = NULL;
- parent->rightNode = NULL;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createLeafNode returning %p", parent);
-#endif
- return parent;
-}
-
-/*
- * Creates an interval with the total extents of the two given intervals.
- */
-INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2)
-{
- INTERVAL *interval;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("mergeIntervals called with %p, %p", inter1, inter2);
-#endif
- interval = lwalloc(sizeof(INTERVAL));
- interval->max = FP_MAX(inter1->max, inter2->max);
- interval->min = FP_MIN(inter1->min, inter2->min);
-#ifdef PGIS_DEBUG
- lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);
-#endif
- return interval;
-}
-
-/*
- * Creates an interval given the min and max values, in arbitrary order.
- */
-INTERVAL *createInterval(double value1, double value2)
-{
- INTERVAL *interval;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createInterval called with %8.3f, %8.3f", value1, value2);
-#endif
- interval = lwalloc(sizeof(INTERVAL));
- interval->max = FP_MAX(value1, value2);
- interval->min = FP_MIN(value1, value2);
-#ifdef PGIS_DEBUG
- lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);
-#endif
- return interval;
-}
-
-/*
- * Recursively frees the child nodes, the interval and the line before
- * freeing the root node.
- */
-void freeTree(RTREE_NODE *root)
-{
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("freeTree called for %p", root);
-#endif
- if(root->leftNode)
- freeTree(root->leftNode);
- if(root->rightNode)
- freeTree(root->rightNode);
- lwfree(root->interval);
- if(root->segment)
- lwfree(root->segment);
- lwfree(root);
-}
-
-/*
- * Retrieves a collection of line segments given the root and crossing value.
- * The collection is a multilinestring consisting of two point lines
- * representing the segments of the ring that may be crossed by the
- * horizontal projection line at the given y value.
- */
-LWMLINE *findLineSegments(RTREE_NODE *root, double value)
-{
- LWMLINE *tmp, *result;
- LWGEOM **lwgeoms;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("findLineSegments called for tree %p and value %8.3f", root, value);
-#endif
- result = NULL;
-
- if(!isContained(root->interval, value))
- {
-#ifdef PGIS_DEBUG
- lwnotice("findLineSegments %p: not contained.", root);
-#endif
- return NULL;
- }
-
- /* If there is a segment defined for this node, include it. */
- if(root->segment)
- {
-#ifdef PGIS_DEBUG
- lwnotice("findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type));
-#endif
- lwgeoms = lwalloc(sizeof(LWGEOM *));
- lwgeoms[0] = (LWGEOM *)root->segment;
-
-#ifdef PGIS_DEBUG
- lwnotice("Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type));
-#endif
- result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms);
- }
-
- /* If there is a left child node, recursively include its results. */
- if(root->leftNode)
- {
-#ifdef PGIS_DEBUG
- lwnotice("findLineSegments %p: recursing left.", root);
-#endif
- tmp = findLineSegments(root->leftNode, value);
- if(tmp)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));
-#endif
- if(result)
- result = mergeMultiLines(result, tmp);
- else
- result = tmp;
- }
- }
-
- /* Same for any right child. */
- if(root->rightNode)
- {
-#ifdef PGIS_DEBUG
- lwnotice("findLineSegments %p: recursing right.", root);
-#endif
- tmp = findLineSegments(root->rightNode, value);
- if(tmp)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));
-#endif
- if(result)
- result = mergeMultiLines(result, tmp);
- else
- result = tmp;
- }
- }
-
- return result;
-}
-
-/* Merges two multilinestrings into a single multilinestring. */
-LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2)
-{
- LWGEOM **geoms;
- LWCOLLECTION *col;
- int i, j, ngeoms;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type));
-#endif
- ngeoms = line1->ngeoms + line2->ngeoms;
- geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);
-
- j = 0;
- for(i = 0; i < line1->ngeoms; i++, j++)
- {
- geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);
- }
- for(i = 0; i < line2->ngeoms; i++, j++)
- {
- geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);
- }
- col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms);
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type));
-#endif
-
- return (LWMLINE *)col;
-}
-
-/*
- * Returns 1 if min < value <= max, 0 otherwise. */
-uint32 isContained(INTERVAL *interval, double value)
-{
- return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
-}
-
-PG_FUNCTION_INFO_V1(LWGEOM_polygon_index);
-Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS)
-{
- PG_LWGEOM *igeom, *result;
- LWGEOM *geom;
- LWPOLY *poly;
- LWMLINE *mline;
- RTREE_NODE *root;
- double yval;
-
-#ifdef PGIS_DEBUG_CALLS
- int i;
- lwnotice("polygon_index called.");
-#endif
- result = NULL;
- igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
- yval = PG_GETARG_FLOAT8(1);
- geom = lwgeom_deserialize(SERIALIZED_FORM(igeom));
- if(TYPE_GETTYPE(geom->type) != POLYGONTYPE)
- {
- lwgeom_release(geom);
- PG_FREE_IF_COPY(igeom, 0);
- PG_RETURN_NULL();
- }
- poly = (LWPOLY *)geom;
- root = createTree(poly->rings[0]);
-
- mline = findLineSegments(root, yval);
-
-#ifdef PGIS_DEBUG
- lwnotice("mline returned %p %d", mline, TYPE_GETTYPE(mline->type));
- for(i = 0; i < mline->ngeoms; i++)
- {
- lwnotice("geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type));
- }
-#endif
-
- if(mline)
- result = pglwgeom_serialize((LWGEOM *)mline);
-
-#ifdef PGIS_DEBUG
- lwnotice("returning result %p", result);
-#endif
- lwfree(root);
-
- PG_FREE_IF_COPY(igeom, 0);
- lwgeom_release((LWGEOM *)poly);
- lwgeom_release((LWGEOM *)mline);
- PG_RETURN_POINTER(result);
-
-}
-
-RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly)
-{
- RTREE_POLY_CACHE *result;
- int i, length;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("createNewCache called with %p", poly);
-#endif
- result = lwalloc(sizeof(RTREE_POLY_CACHE));
- result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
- result->ringCount = poly->nrings;
- length = lwgeom_size_poly(serializedPoly);
- result->poly = lwalloc(length);
- memcpy(result->poly, serializedPoly, length);
- for(i = 0; i < result->ringCount; i++)
- {
- result->ringIndices[i] = createTree(poly->rings[i]);
- }
-#ifdef PGIS_DEBUG
- lwnotice("createNewCache returning %p", result);
-#endif
- return result;
-}
-
-/*
- * Creates a new cachable index if needed, or returns the current cache if
- * it is applicable to the current polygon.
- * The memory context must be changed to function scope before calling this
- * method. The method will allocate memory for the cache it creates,
- * as well as freeing the memory of any cache that is no longer applicable.
- */
-RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly,
- RTREE_POLY_CACHE *currentCache)
-{
- int i, length;
-
-#ifdef PGIS_DEBUG_CALLS
- lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache);
-#endif
- if(!currentCache)
- {
-#ifdef PGIS_DEBUG
- lwnotice("No existing cache, create one.");
-#endif
- return createNewCache(poly, serializedPoly);
- }
- if(!(currentCache->poly))
- {
-#ifdef PGIS_DEBUG
- lwnotice("Cache contains no polygon, creating new cache.");
-#endif
- return createNewCache(poly, serializedPoly);
- }
-
- length = lwgeom_size_poly(serializedPoly);
-
- if(lwgeom_size_poly(currentCache->poly) != length)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Polygon size mismatch, creating new cache.");
-#endif
- lwfree(currentCache->poly);
- lwfree(currentCache);
- return createNewCache(poly, serializedPoly);
- }
- for(i = 0; i < length; i++)
- {
- uchar a = serializedPoly[i];
- uchar b = currentCache->poly[i];
- if(a != b)
- {
-#ifdef PGIS_DEBUG
- lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b);
-#endif
- lwfree(currentCache->poly);
- lwfree(currentCache);
- return createNewCache(poly, serializedPoly);
- }
- }
-
-#ifdef PGIS_DEBUG
- lwnotice("Polygon match, retaining current cache, %p.", currentCache);
-#endif
- return currentCache;
-}
-
+/**********************************************************************
+ * $Id$
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2001-2005 Refractions Research Inc.
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU General Public Licence. See the COPYING file.
+ *
+ **********************************************************************/
+
+#include "lwgeom_pg.h"
+#include "liblwgeom.h"
+#include "lwgeom_rtree.h"
+
+
+Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS);
+RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly);
+
+
+/*
+ * Creates an rtree given a pointer to the point array.
+ * Must copy the point array.
+ */
+RTREE_NODE *createTree(POINTARRAY *pointArray)
+{
+ RTREE_NODE *root;
+ RTREE_NODE** nodes = lwalloc(sizeof(RTREE_NODE*) * pointArray->npoints);
+ int i, nodeCount;
+ int childNodes, parentNodes;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createTree called with pointarray %p", pointArray);
+#endif
+
+ nodeCount = pointArray->npoints - 1;
+
+#ifdef PGIS_DEBUG
+ lwnotice("Total leaf nodes: %d", nodeCount);
+#endif
+
+ /*
+ * Create a leaf node for every line segment.
+ */
+ for(i = 0; i < nodeCount; i++)
+ {
+ nodes[i] = createLeafNode(pointArray, i);
+ }
+
+ /*
+ * Next we group nodes by pairs. If there's an odd number of nodes,
+ * we bring the last node up a level as is. Continue until we have
+ * a single top node.
+ */
+ childNodes = nodeCount;
+ parentNodes = nodeCount / 2;
+ while(parentNodes > 0)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Merging %d children into %d parents.", childNodes, parentNodes);
+#endif
+ i = 0;
+ while(i < parentNodes)
+ {
+ nodes[i] = createInteriorNode(nodes[i*2], nodes[i*2+1]);
+ i++;
+ }
+ /*
+ * Check for an odd numbered final node.
+ */
+ if(parentNodes * 2 < childNodes)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Shuffling child %d to parent %d", childNodes - 1, i);
+#endif
+ nodes[i] = nodes[childNodes - 1];
+ parentNodes++;
+ }
+ childNodes = parentNodes;
+ parentNodes = parentNodes / 2;
+ }
+
+ root = nodes[0];
+
+#ifdef PGIS_DEBUG
+ lwnotice("createTree returning %p", root);
+#endif
+ return root;
+}
+
+/*
+ * Creates an interior node given the children.
+ */
+RTREE_NODE *createInteriorNode(RTREE_NODE *left, RTREE_NODE *right){
+ RTREE_NODE *parent;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createInteriorNode called for children %p, %p", left, right);
+#endif
+ parent = lwalloc(sizeof(RTREE_NODE));
+ parent->leftNode = left;
+ parent->rightNode = right;
+ parent->interval = mergeIntervals(left->interval, right->interval);
+ parent->segment = NULL;
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createInteriorNode returning %p", parent);
+#endif
+ return parent;
+}
+
+/*
+ * Creates a leaf node given the pointer to the start point of the segment.
+ */
+RTREE_NODE *createLeafNode(POINTARRAY *pa, int startPoint)
+{
+ RTREE_NODE *parent;
+ LWLINE *line;
+ double value1;
+ double value2;
+ POINT4D tmp;
+ POINTARRAY *npa;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createLeafNode called for point %d of %p", startPoint, pa);
+#endif
+ if(pa->npoints < startPoint + 2)
+ {
+ lwerror("createLeafNode: npoints = %d, startPoint = %d", pa->npoints, startPoint);
+ }
+
+ /*
+ * The given point array will be part of a geometry that will be freed
+ * independently of the index. Since we may want to cache the index,
+ * we must create independent arrays.
+ */
+ npa = lwalloc(sizeof(POINTARRAY));
+ npa->dims = 0;
+ npa->npoints = 2;
+ TYPE_SETZM(npa->dims, 0, 0);
+ npa->serialized_pointlist = lwalloc(pointArray_ptsize(pa) * 2);
+
+ getPoint4d_p(pa, startPoint, &tmp);
+
+ setPoint4d(npa, 0, &tmp);
+ value1 = tmp.y;
+
+ getPoint4d_p(pa, startPoint + 1, &tmp);
+
+ setPoint4d(npa, 1, &tmp);
+ value2 = tmp.y;
+
+ line = lwline_construct(-1, NULL, npa);
+
+ parent = lwalloc(sizeof(RTREE_NODE));
+ parent->interval = createInterval(value1, value2);
+ parent->segment = line;
+ parent->leftNode = NULL;
+ parent->rightNode = NULL;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createLeafNode returning %p", parent);
+#endif
+ return parent;
+}
+
+/*
+ * Creates an interval with the total extents of the two given intervals.
+ */
+INTERVAL *mergeIntervals(INTERVAL *inter1, INTERVAL *inter2)
+{
+ INTERVAL *interval;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("mergeIntervals called with %p, %p", inter1, inter2);
+#endif
+ interval = lwalloc(sizeof(INTERVAL));
+ interval->max = FP_MAX(inter1->max, inter2->max);
+ interval->min = FP_MIN(inter1->min, inter2->min);
+#ifdef PGIS_DEBUG
+ lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);
+#endif
+ return interval;
+}
+
+/*
+ * Creates an interval given the min and max values, in arbitrary order.
+ */
+INTERVAL *createInterval(double value1, double value2)
+{
+ INTERVAL *interval;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createInterval called with %8.3f, %8.3f", value1, value2);
+#endif
+ interval = lwalloc(sizeof(INTERVAL));
+ interval->max = FP_MAX(value1, value2);
+ interval->min = FP_MIN(value1, value2);
+#ifdef PGIS_DEBUG
+ lwnotice("interval min = %8.3f, max = %8.3f", interval->min, interval->max);
+#endif
+ return interval;
+}
+
+/*
+ * Recursively frees the child nodes, the interval and the line before
+ * freeing the root node.
+ */
+void freeTree(RTREE_NODE *root)
+{
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("freeTree called for %p", root);
+#endif
+ if(root->leftNode)
+ freeTree(root->leftNode);
+ if(root->rightNode)
+ freeTree(root->rightNode);
+ lwfree(root->interval);
+ if(root->segment)
+ lwfree(root->segment);
+ lwfree(root);
+}
+
+/*
+ * Retrieves a collection of line segments given the root and crossing value.
+ * The collection is a multilinestring consisting of two point lines
+ * representing the segments of the ring that may be crossed by the
+ * horizontal projection line at the given y value.
+ */
+LWMLINE *findLineSegments(RTREE_NODE *root, double value)
+{
+ LWMLINE *tmp, *result;
+ LWGEOM **lwgeoms;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("findLineSegments called for tree %p and value %8.3f", root, value);
+#endif
+ result = NULL;
+
+ if(!isContained(root->interval, value))
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("findLineSegments %p: not contained.", root);
+#endif
+ return NULL;
+ }
+
+ /* If there is a segment defined for this node, include it. */
+ if(root->segment)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("findLineSegments %p: adding segment %p %d.", root, root->segment, TYPE_GETTYPE(root->segment->type));
+#endif
+ lwgeoms = lwalloc(sizeof(LWGEOM *));
+ lwgeoms[0] = (LWGEOM *)root->segment;
+
+#ifdef PGIS_DEBUG
+ lwnotice("Found geom %p, type %d, dim %d", root->segment, TYPE_GETTYPE(root->segment->type), TYPE_GETZM(root->segment->type));
+#endif
+ result = (LWMLINE *)lwcollection_construct(lwgeom_makeType_full(0, 0, 0, MULTILINETYPE, 0), -1, NULL, 1, lwgeoms);
+ }
+
+ /* If there is a left child node, recursively include its results. */
+ if(root->leftNode)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("findLineSegments %p: recursing left.", root);
+#endif
+ tmp = findLineSegments(root->leftNode, value);
+ if(tmp)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));
+#endif
+ if(result)
+ result = mergeMultiLines(result, tmp);
+ else
+ result = tmp;
+ }
+ }
+
+ /* Same for any right child. */
+ if(root->rightNode)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("findLineSegments %p: recursing right.", root);
+#endif
+ tmp = findLineSegments(root->rightNode, value);
+ if(tmp)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Found geom %p, type %d, dim %d", tmp, TYPE_GETTYPE(tmp->type), TYPE_GETZM(tmp->type));
+#endif
+ if(result)
+ result = mergeMultiLines(result, tmp);
+ else
+ result = tmp;
+ }
+ }
+
+ return result;
+}
+
+/* Merges two multilinestrings into a single multilinestring. */
+LWMLINE *mergeMultiLines(LWMLINE *line1, LWMLINE *line2)
+{
+ LWGEOM **geoms;
+ LWCOLLECTION *col;
+ int i, j, ngeoms;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("mergeMultiLines called on %p, %d, %d; %p, %d, %d", line1, line1->ngeoms, TYPE_GETTYPE(line1->type), line2, line2->ngeoms, TYPE_GETTYPE(line2->type));
+#endif
+ ngeoms = line1->ngeoms + line2->ngeoms;
+ geoms = lwalloc(sizeof(LWGEOM *) * ngeoms);
+
+ j = 0;
+ for(i = 0; i < line1->ngeoms; i++, j++)
+ {
+ geoms[j] = lwgeom_clone((LWGEOM *)line1->geoms[i]);
+ }
+ for(i = 0; i < line2->ngeoms; i++, j++)
+ {
+ geoms[j] = lwgeom_clone((LWGEOM *)line2->geoms[i]);
+ }
+ col = lwcollection_construct(MULTILINETYPE, -1, NULL, ngeoms, geoms);
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("mergeMultiLines returning %p, %d, %d", col, col->ngeoms, TYPE_GETTYPE(col->type));
+#endif
+
+ return (LWMLINE *)col;
+}
+
+/*
+ * Returns 1 if min < value <= max, 0 otherwise. */
+uint32 isContained(INTERVAL *interval, double value)
+{
+ return FP_CONTAINS_INCL(interval->min, value, interval->max) ? 1 : 0;
+}
+
+PG_FUNCTION_INFO_V1(LWGEOM_polygon_index);
+Datum LWGEOM_polygon_index(PG_FUNCTION_ARGS)
+{
+ PG_LWGEOM *igeom, *result;
+ LWGEOM *geom;
+ LWPOLY *poly;
+ LWMLINE *mline;
+ RTREE_NODE *root;
+ double yval;
+
+#ifdef PGIS_DEBUG_CALLS
+ int i;
+ lwnotice("polygon_index called.");
+#endif
+ result = NULL;
+ igeom = (PG_LWGEOM *)PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
+ yval = PG_GETARG_FLOAT8(1);
+ geom = lwgeom_deserialize(SERIALIZED_FORM(igeom));
+ if(TYPE_GETTYPE(geom->type) != POLYGONTYPE)
+ {
+ lwgeom_release(geom);
+ PG_FREE_IF_COPY(igeom, 0);
+ PG_RETURN_NULL();
+ }
+ poly = (LWPOLY *)geom;
+ root = createTree(poly->rings[0]);
+
+ mline = findLineSegments(root, yval);
+
+#ifdef PGIS_DEBUG
+ lwnotice("mline returned %p %d", mline, TYPE_GETTYPE(mline->type));
+ for(i = 0; i < mline->ngeoms; i++)
+ {
+ lwnotice("geom[%d] %p %d", i, mline->geoms[i], TYPE_GETTYPE(mline->geoms[i]->type));
+ }
+#endif
+
+ if(mline)
+ result = pglwgeom_serialize((LWGEOM *)mline);
+
+#ifdef PGIS_DEBUG
+ lwnotice("returning result %p", result);
+#endif
+ lwfree(root);
+
+ PG_FREE_IF_COPY(igeom, 0);
+ lwgeom_release((LWGEOM *)poly);
+ lwgeom_release((LWGEOM *)mline);
+ PG_RETURN_POINTER(result);
+
+}
+
+RTREE_POLY_CACHE *createNewCache(LWPOLY *poly, uchar *serializedPoly)
+{
+ RTREE_POLY_CACHE *result;
+ int i, length;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("createNewCache called with %p", poly);
+#endif
+ result = lwalloc(sizeof(RTREE_POLY_CACHE));
+ result->ringIndices = lwalloc(sizeof(RTREE_NODE *) * poly->nrings);
+ result->ringCount = poly->nrings;
+ length = lwgeom_size_poly(serializedPoly);
+ result->poly = lwalloc(length);
+ memcpy(result->poly, serializedPoly, length);
+ for(i = 0; i < result->ringCount; i++)
+ {
+ result->ringIndices[i] = createTree(poly->rings[i]);
+ }
+#ifdef PGIS_DEBUG
+ lwnotice("createNewCache returning %p", result);
+#endif
+ return result;
+}
+
+/*
+ * Creates a new cachable index if needed, or returns the current cache if
+ * it is applicable to the current polygon.
+ * The memory context must be changed to function scope before calling this
+ * method. The method will allocate memory for the cache it creates,
+ * as well as freeing the memory of any cache that is no longer applicable.
+ */
+RTREE_POLY_CACHE *retrieveCache(LWPOLY *poly, uchar *serializedPoly,
+ RTREE_POLY_CACHE *currentCache)
+{
+ int i, length;
+
+#ifdef PGIS_DEBUG_CALLS
+ lwnotice("retrieveCache called with %p %p %p", poly, serializedPoly, currentCache);
+#endif
+ if(!currentCache)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("No existing cache, create one.");
+#endif
+ return createNewCache(poly, serializedPoly);
+ }
+ if(!(currentCache->poly))
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Cache contains no polygon, creating new cache.");
+#endif
+ return createNewCache(poly, serializedPoly);
+ }
+
+ length = lwgeom_size_poly(serializedPoly);
+
+ if(lwgeom_size_poly(currentCache->poly) != length)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Polygon size mismatch, creating new cache.");
+#endif
+ lwfree(currentCache->poly);
+ lwfree(currentCache);
+ return createNewCache(poly, serializedPoly);
+ }
+ for(i = 0; i < length; i++)
+ {
+ uchar a = serializedPoly[i];
+ uchar b = currentCache->poly[i];
+ if(a != b)
+ {
+#ifdef PGIS_DEBUG
+ lwnotice("Polygon mismatch, creating new cache. %c, %c", a, b);
+#endif
+ lwfree(currentCache->poly);
+ lwfree(currentCache);
+ return createNewCache(poly, serializedPoly);
+ }
+ }
+
+#ifdef PGIS_DEBUG
+ lwnotice("Polygon match, retaining current cache, %p.", currentCache);
+#endif
+ return currentCache;
+}
+
Modified: trunk/lwgeom/lwgeom_sqlmm.c
===================================================================
--- trunk/lwgeom/lwgeom_sqlmm.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_sqlmm.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -1,1176 +1,1197 @@
-/**********************************************************************
- * $Id$
- *
- * PostGIS - Spatial Types for PostgreSQL
- * http://postgis.refractions.net
- * Copyright 2001-2006 Refractions Research Inc.
- *
- * This is free software; you can redistribute and/or modify it under
- * the terms of the GNU General Public Licence. See the COPYING file.
- *
- **********************************************************************/
-
-#include <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"
-
-/*
- * 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
- ******************************************************************************/
+/**********************************************************************
+ * $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
+ ******************************************************************************/
Modified: trunk/lwgeom/lwgeom_transform.c
===================================================================
--- trunk/lwgeom/lwgeom_transform.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/lwgeom_transform.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -28,6 +28,7 @@
Datum transform_geom(PG_FUNCTION_ARGS);
Datum postgis_proj_version(PG_FUNCTION_ARGS);
+
/* if POSTGIS_PROJ_VERSION undefined, we get a do-nothing transform() function */
#ifndef POSTGIS_PROJ_VERSION
@@ -145,7 +146,7 @@
/* Search path for PROJ.4 library */
static bool IsPROJ4LibPathSet = false;
-void SetPROJ4LibPath();
+void SetPROJ4LibPath(void);
/* Memory context cache functions */
static void PROJ4SRSCacheInit(MemoryContext context);
@@ -244,7 +245,7 @@
{
/*
* Do nothing, but we must supply a function since this call is mandatory according to tgl
- * (see postgis-devel archives July 2007)
+ * (see postgis-devel archives July 2007)
*/
}
@@ -253,7 +254,7 @@
{
/*
* Always return false since this call is mandatory according to tgl
- * (see postgis-devel archives July 2007)
+ * (see postgis-devel archives July 2007)
*/
return FALSE;
}
@@ -263,7 +264,7 @@
{
/*
* Simple stats display function - we must supply a function since this call is mandatory according to tgl
- * (see postgis-devel archives July 2007)
+ * (see postgis-devel archives July 2007)
*/
fprintf(stderr, "%s: PROJ4 context\n", context->name);
@@ -580,7 +581,7 @@
* 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;
Modified: trunk/lwgeom/wktunparse.c
===================================================================
--- trunk/lwgeom/wktunparse.c 2008-05-22 14:34:32 UTC (rev 2780)
+++ trunk/lwgeom/wktunparse.c 2008-05-22 20:43:00 UTC (rev 2781)
@@ -42,6 +42,8 @@
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);
More information about the postgis-commits
mailing list