*/
extern int getPoint2d_p_ro(const POINTARRAY *pa, int n, POINT2D **point);
-/**
-* @brief Check whether or not a lwgeom is big enough to warrant a bounding box.
-*
-* Check whether or not a lwgeom is big enough to warrant a bounding box
-* when stored in the serialized form on disk. Currently only points are
-* considered small enough to not require a bounding box, because the
-* index operations can generate a large number of box-retrieval operations
-* when scanning keys.
-*/
-extern int lwgeom_needs_bbox(LWGEOM *geom);
-/**
-* Count the total number of vertices in any #LWGEOM.
-*/
-extern int lwgeom_count_vertices(LWGEOM *geom);
/**
* Calculate box and add values to gbox. Return #G_SUCCESS on success.
** Floating point comparitors.
*/
#define FP_TOLERANCE 1e-12
-#define FP_ZERO(A) (fabs(A) <= FP_TOLERANCE)
+#define FP_IS_ZERO(A) (fabs(A) <= FP_TOLERANCE)
#define FP_MAX(A, B) (((A) > (B)) ? (A) : (B))
#define FP_MIN(A, B) (((A) < (B)) ? (A) : (B))
#define FP_EQUALS(A, B) (fabs((A)-(B)) <= FP_TOLERANCE)
extern double lwgeom_mindistance2d_recursive(uchar *lw1, uchar *lw2);
extern double lwgeom_mindistance2d_recursive_tolerance(uchar *lw1, uchar *lw2, double tolerance);
extern int lwgeom_pt_inside_circle(POINT2D *p, double cx, double cy, double rad);
-extern int32 lwgeom_npoints(uchar *serialized);
extern char ptarray_isccw(const POINTARRAY *pa);
extern void lwgeom_reverse(LWGEOM *lwgeom);
extern void lwline_reverse(LWLINE *line);
extern int lwgeom_compute_box2d_p(LWGEOM *lwgeom, BOX2DFLOAT4 *box);
void lwgeom_longitude_shift(LWGEOM *lwgeom);
+
+/**
+* @brief Check whether or not a lwgeom is big enough to warrant a bounding box.
+*
+* Check whether or not a lwgeom is big enough to warrant a bounding box
+* when stored in the serialized form on disk. Currently only points are
+* considered small enough to not require a bounding box, because the
+* index operations can generate a large number of box-retrieval operations
+* when scanning keys.
+*/
+extern int lwgeom_needs_bbox(LWGEOM *geom);
+
+/**
+* Count the total number of vertices in any #LWGEOM.
+*/
+extern int lwgeom_count_vertices(LWGEOM *geom);
+extern int32 lwgeom_npoints(uchar *serialized);
+
/* Is lwgeom1 geometrically equal to lwgeom2 ? */
char lwgeom_same(const LWGEOM *lwgeom1, const LWGEOM *lwgeom2);
char ptarray_same(const POINTARRAY *pa1, const POINTARRAY *pa2);
--- /dev/null
+/**********************************************************************
+ * $Id: liblwgeom.h 4497 2009-09-14 18:33:54Z pramsey $
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * http://postgis.refractions.net
+ * Copyright 2001-2006 Refractions Research Inc.
+ * Copyright 2007-2008 Mark Cave-Ayland
+ * Copyright 2008 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * 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 <assert.h>
+
+#include "liblwgeom.h"
+
+
**********************************************************************/
#include <math.h>
-#include "liblwgeom.h"
+#include "liblwgeom_internal.h"
enum CG_SEGMENT_INTERSECTION_TYPE {
SEG_ERROR = -1,
*
**********************************************************************/
-#include <math.h>
-#include "libgeom.h"
+#include "lwgeodetic.h"
+
+
+/**
+* Given two points on a unit sphere, calculate their distance apart.
+*/
+double sphere_distance(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
+{
+ double latS = s.lat;
+ double latE = e.lon;
+ double dlng = e.lon - s.lon;
+ double a1 = pow(cos(latE) * sin(dlng), 2.0);
+ double a2 = pow(cos(latS) * sin(latE) - sin(latS) * cos(latE) * cos(dlng), 2.0);
+ double a = sqrt(a1 + a2);
+ double b = sin(latS) * sin(latE) + cos(latS) * cos(latE) * cos(dlng);
+ return atan2(a, b);
+}
+
+/**
+* Given two points on a unit sphere, calculate the direction from s to e.
+*/
+double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e)
+{
+ double latS = s.lat;
+ double latE = e.lon;
+ double dlng = e.lat - s.lon;
+ double heading = atan2(sin(dlng) * cos(latE),
+ cos(latS) * sin(latE) -
+ sin(latS) * cos(latE) * cos(dlng)) / PI;
+ return heading;
+}
+
+/**
+* Given a starting location r, a distance and an azimuth
+* to the new point, compute the location of the projected point on the unit sphere.
+*/
+void sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n)
+{
+ double d = distance;
+ double lat1 = r.lat;
+ n->lat = asin(sin(lat1) * cos(d) +
+ cos(lat1) * sin(d) * cos(azimuth));
+ double a = cos(lat1) * cos(d) - sin(lat1) * sin(d) * cos(azimuth);
+ double b = signum(d) * sin(azimuth);
+ n->lon = atan(b/a) + r.lon;
+}
+
+void sphere_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox)
+{
+}
+
/*
* This function can only be used on LWGEOM that is built on top of
}
-/*
-** Count points in an LWGEOM.
-*/
-static int lwcollection_count_vertices(LWCOLLECTION *col)
-{
- int i = 0;
- int v = 0; /* vertices */
- assert(col);
- for ( i = 0; i < col->ngeoms; i++ )
- {
- v += lwgeom_count_vertices(col->geoms[i]);
- }
- return v;
-}
-
-static int lwpolygon_count_vertices(LWPOLY *poly)
-{
- int i = 0;
- int v = 0; /* vertices */
- assert(poly);
- for ( i = 0; i < poly->nrings; i ++ )
- {
- v += poly->rings[i]->npoints;
- }
- return v;
-}
-
-static int lwline_count_vertices(LWLINE *line)
-{
- assert(line);
- if ( ! line->points )
- return 0;
- return line->points->npoints;
-}
-
-static int lwpoint_count_vertices(LWPOINT *point)
-{
- assert(point);
- if ( ! point->point )
- return 0;
- return 1;
-}
-
-int lwgeom_count_vertices(LWGEOM *geom)
-{
- int result = 0;
- LWDEBUGF(4, "got type %d", TYPE_GETTYPE(geom->type));
- switch (TYPE_GETTYPE(geom->type))
- {
- case POINTTYPE:
- result = lwpoint_count_vertices((LWPOINT *)geom);;
- break;
- case LINETYPE:
- result = lwline_count_vertices((LWLINE *)geom);
- break;
- case POLYGONTYPE:
- result = lwpolygon_count_vertices((LWPOLY *)geom);
- break;
- case MULTIPOINTTYPE:
- case MULTILINETYPE:
- case MULTIPOLYGONTYPE:
- case COLLECTIONTYPE:
- result = lwcollection_count_vertices((LWCOLLECTION *)geom);
- break;
- default:
- lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type));
- break;
- }
- LWDEBUGF(3, "counted %d vertices", result);
- return result;
-}
-
-int lwgeom_needs_bbox(LWGEOM *geom)
-{
- assert(geom);
- if( TYPE_GETTYPE(geom->type) == POINTTYPE )
- {
- return G_FALSE;
- }
- return G_TRUE;
-}
--- /dev/null
+/**********************************************************************
+ * $Id: lwgeodetic.c 4494 2009-09-14 10:54:33Z mcayland $
+ *
+ * PostGIS - Spatial Types for PostgreSQL
+ * Copyright 2009 Paul Ramsey <pramsey@cleverelephant.ca>
+ *
+ * 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 <math.h>
+#include "libgeom.h"
+
+/**
+* Point in spherical coordinates on the world. Units of radians.
+*/
+typedef struct {
+ double lon;
+ double lat;
+} GEOGRAPHIC_POINT;
+
+/**
+* Two-point great circle segment from a to b.
+*/
+typedef struct {
+ GEOGRAPHIC_POINT a;
+ GEOGRAPHIC_POINT b;
+} GEOGRAPHIC_EDGE;
+
+/**
+* Conversion functions
+*/
+#define deg2rad(d) (PI * (d) / 180.0)
+#define rad2deg(r) (180.0 * (r) / PI)
+
+/**
+* Ape a java function
+*/
+#define signum(a) ((a) < 0 ? -1 : ((a) > 0 ? 1 : (a)))
+
+/*
+** Prototypes for internal functions.
+*/
+double sphere_distance(GEOGRAPHIC_POINT a, GEOGRAPHIC_POINT b);
+double sphere_direction(GEOGRAPHIC_POINT s, GEOGRAPHIC_POINT e);
+void sphere_project(GEOGRAPHIC_POINT r, double distance, double azimuth, GEOGRAPHIC_POINT *n);
+void sphere_gbox(GEOGRAPHIC_EDGE e, GBOX *gbox);
+
+
#include <stdlib.h>
#include <stdarg.h>
-#include "liblwgeom.h"
+#include "liblwgeom_internal.h"
#include "wktparse.h"
return;
};
+
+
+int lwgeom_needs_bbox(LWGEOM *geom)
+{
+ assert(geom);
+ if( TYPE_GETTYPE(geom->type) == POINTTYPE )
+ {
+ return LW_FALSE;
+ }
+ return LW_TRUE;
+}
+
+
+/*
+** Count points in an LWGEOM.
+*/
+
+static int lwcollection_count_vertices(LWCOLLECTION *col)
+{
+ int i = 0;
+ int v = 0; /* vertices */
+ assert(col);
+ for ( i = 0; i < col->ngeoms; i++ )
+ {
+ v += lwgeom_count_vertices(col->geoms[i]);
+ }
+ return v;
+}
+
+static int lwpolygon_count_vertices(LWPOLY *poly)
+{
+ int i = 0;
+ int v = 0; /* vertices */
+ assert(poly);
+ for ( i = 0; i < poly->nrings; i ++ )
+ {
+ v += poly->rings[i]->npoints;
+ }
+ return v;
+}
+
+static int lwline_count_vertices(LWLINE *line)
+{
+ assert(line);
+ if ( ! line->points )
+ return 0;
+ return line->points->npoints;
+}
+
+static int lwpoint_count_vertices(LWPOINT *point)
+{
+ assert(point);
+ if ( ! point->point )
+ return 0;
+ return 1;
+}
+
+int lwgeom_count_vertices(LWGEOM *geom)
+{
+ int result = 0;
+ LWDEBUGF(4, "got type %d", TYPE_GETTYPE(geom->type));
+ switch (TYPE_GETTYPE(geom->type))
+ {
+ case POINTTYPE:
+ result = lwpoint_count_vertices((LWPOINT *)geom);;
+ break;
+ case LINETYPE:
+ result = lwline_count_vertices((LWLINE *)geom);
+ break;
+ case POLYGONTYPE:
+ result = lwpolygon_count_vertices((LWPOLY *)geom);
+ break;
+ case MULTIPOINTTYPE:
+ case MULTILINETYPE:
+ case MULTIPOLYGONTYPE:
+ case COLLECTIONTYPE:
+ result = lwcollection_count_vertices((LWCOLLECTION *)geom);
+ break;
+ default:
+ lwerror("unsupported input geometry type: %d", TYPE_GETTYPE(geom->type));
+ break;
+ }
+ LWDEBUGF(3, "counted %d vertices", result);
+ return result;
+}
*result = size_union - size_orig;
/* All things being equal, we prefer to expand small boxes rather than large boxes. */
- if( FP_ZERO(*result) )
- if( FP_ZERO(size_orig) )
+ if( FP_IS_ZERO(*result) )
+ if( FP_IS_ZERO(size_orig) )
*result = 0.0;
else
*result = 1.0 - (1.0/(1.0 + size_orig));
projects / postgis / commitdiff