[postgis-users] Re: C question, PG_RETURN_POINTER: GPC integration
Nicolas Ribot
nicolas.ribot at scot.cnes.fr
Tue Apr 30 07:40:30 PDT 2002
Oooooopss !!
I made a mistake with attachements !
Here are the correct files
Sorry for the inconvenience
;-(
Nico
-------------- next part --------------
#include "utils/geo_decls.h"
#include "postgis.h"
#include "gpc.h"
// taken from gpc.c to avoid inclusion:
#define MALLOC(p, b, s) {if ((b) > 0) { \
p= malloc(b); if (!(p)) { \
fprintf(stderr, "gpc malloc failure: %s\n", s); \
exit(0);}} else p= NULL;}
/*************************************************
***************EXPOSED TO PSQL********************
**************************************************/
Datum union_pg_pg (PG_FUNCTION_ARGS);
Datum intersection_pg_pg (PG_FUNCTION_ARGS);
Datum difference_pg_pg (PG_FUNCTION_ARGS);
Datum sym_difference_pg_pg(PG_FUNCTION_ARGS);
/*************************************************
**********CONVERSION GPC <-> POSTGIS types********
**************************************************/
GEOMETRY* scot_gpc_clipping (POLYGON3D *poly1, POLYGON3D *poly2, gpc_op operation, int SRID);
GEOMETRY* gpc_to_GEOMETRY (gpc_polygon *p, int SRID);
void POLYGON3D_to_gpc (POLYGON3D *ppoly, gpc_polygon *p);
POLYGON3D* gpc_to_POLYGON3D (gpc_polygon *p, int *rings, int num_rings, int SRID, int *size);
bool is_ccw (double x0, double y0, double x1, double y1, double x2, double y2);
bool contour_in_POLYGON3D (POLYGON3D *ppoly, gpc_polygon *p, int i);
int point_in_gpc_poly ( gpc_vertex *P, gpc_vertex *V, int n );
bool gpc_hole_in_gpc_contour(gpc_polygon *p, int c, int v);
char* gpc_polygon_to_text (gpc_polygon *p);
GEOMETRY* add_geometry(GEOMETRY *geom1, GEOMETRY *geom2);
// the error message generated by one conversion function: not used
char *error_message = "no error";
/*************************************************
***************TEST ZONE**************************
**************************************************/
// geo functions
Datum test_geogpc(PG_FUNCTION_ARGS);
// transformation functions:
void scot_gpc_read_polygon(POLYGON3D *ppoly, gpc_polygon *p);
//void scot_gpc_write_polygon(gpc_polygon *p, POLYGON3D *ppoly, int *ppoly_size);
POLYGON3D * scot_gpc_write_polygon(gpc_polygon *p, int *ppoly_size);
-------------- next part --------------
/*************************************************************
* File: scot_gpc.c
* Author: Nicolas Ribot
* Date: 25 avril 2002
* Version: 0.3
* Description: Providing OpenGIS UDF for polygon manipulation:
* intersection, union, difference, symmetric_difference.
* Some useful method to achieve that:
* is_ccw, gpc<->postigs conversion functions
based on GPC clipping library.
*
* History: 0.1 (25 apr 2002) Original definition with a test function
* 0.2 (28 apr 2002) Added ogis UDF, conversion functions. Testing.
*
* Bugs: gps does not consider that a hole touching a contour at one point is a hole.
* Opengis does... => with such objects, the resulting postgis poly will be
* invalid.
*
* ToDo: Manage objects deallocation
* Provide a dynamic memory allocation for resulting postgis polygon
* Write a debug function for gpc polygons (char*)
-> done
* Merge gpc_to_POLYGON3D_tmp code into gpc_to_POLYGON3D
-> done
* Implements a point_in_poly for gpc polygon, assuming only one contour (no holes)
*
Copyright: ? 2002 Nicolas Ribot.
?*************************************************************/
#include "postgres.h"
#include <math.h>
#include <float.h>
#include <string.h>
#include <stdio.h>
#include <errno.h>
#include "access/gist.h"
#include "access/itup.h"
#include "utils/elog.h"
#include "scot_gpc.h"
/*************************************************
***************EXPOSED TO PSQL********************
**************************************************/
/*
* performs the union of 2 polygons
* supports MULTIPOLYGON with only one polygon in it
* Returns null if input types are not polygons, or
* different SRIDS
* Returns a polygon or multipolygon according to result
*/
PG_FUNCTION_INFO_V1(union_pg_pg);
Datum union_pg_pg(PG_FUNCTION_ARGS){
GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
GEOMETRY *result;
// the conversion of geometry to POLYGON3D
POLYGON3D *poly1;
POLYGON3D *poly2;
int32 *offset1;
int32 *offset2;
char *o1;
char *o2;
elog(LOG,"\n--------------------------------------------------------------");
// looks at the input geometries type
if (geom1->type != POLYGONTYPE || geom2->type != POLYGONTYPE) {
elog(ERROR,"Incompatible types for 2 geometries (expected POLYGON)");
PG_RETURN_NULL();
}
// should add support for multipolygon with 1 object in it
// should provide a function to add 2 input geom into the result in case of
// disjoint geometries.
// creates the 2 polygons
offset1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
offset2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
o1 = (char *) geom1 +offset1[0];
o2 = (char *) geom2 +offset2[0];
poly1 = (POLYGON3D *)o1;
poly2 = (POLYGON3D *)o2;
// are geometries equals ?
if (is_same_polygon(poly1, poly2)) {
// A U B = A or B if A=B
elog(LOG, "union_pg_pg: same input polygons, returning first one");
PG_RETURN_POINTER(geom1);
}
elog(LOG,"union_pg_pg: input polygon 1: %s", geometry_to_text(geom1));
elog(LOG,"union_pg_pg: input polygon 2: %s", geometry_to_text(geom2));
result = scot_gpc_clipping(poly1, poly2, GPC_UNION, geom1->SRID);
if (result == NULL) {
// empty result set, not an error.
PG_RETURN_NULL();
} else {
result->SRID = geom1->SRID;
PG_RETURN_POINTER(result);
}
}
/*
* performs the intersection of 2 polygons
* supports MULTIPOLYGON with only one polygon in it
* Returns null if input types are not polygons, or
* different SRIDS
* Returns a polygon or multipolygon according to result
*/
PG_FUNCTION_INFO_V1(intersection_pg_pg);
Datum intersection_pg_pg(PG_FUNCTION_ARGS){
GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
GEOMETRY *result;
// the conversion of geometry to POLYGON3D
POLYGON3D *poly1;
POLYGON3D *poly2;
int32 *offset1;
int32 *offset2;
char *o1;
char *o2;
elog(LOG,"\n--------------------------------------------------------------");
// looks at the input geometries type
if (geom1->type != POLYGONTYPE || geom2->type != POLYGONTYPE) {
elog(ERROR,"Incompatible types for 2 geometries (expected POLYGON)");
PG_RETURN_NULL();
}
// should add support for multipolygon with 1 object in it
if (geom1->SRID != geom2->SRID) {
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs");
PG_RETURN_NULL();
}
// the basic optimizations:
// are their bbox disjoints ?
if (!box3d_ov(&(geom1->bvol), &(geom2->bvol))) {
elog(LOG,"Geometries bbox disjoint");
PG_RETURN_NULL();
}
// creates the 2 polygons
offset1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
offset2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
o1 = (char *) geom1 +offset1[0];
o2 = (char *) geom2 +offset2[0];
poly1 = (POLYGON3D *)o1;
poly2 = (POLYGON3D *)o2;
// are geometries equals ?
if (is_same_polygon(poly1, poly2)) {
elog(LOG, "same input polygons, returning first one");
PG_RETURN_POINTER(geom1);
}
elog(LOG,"intersection_pg_pg: input polygon 1: %s", geometry_to_text(geom1));
elog(LOG,"intersection_pg_pg: input polygon 2: %s", geometry_to_text(geom2));
result = scot_gpc_clipping(poly1, poly2, GPC_INT, geom1->SRID);
if (result == NULL) {
// empty result set, not an error.
PG_RETURN_NULL();
} else {
result->SRID = geom1->SRID;
PG_RETURN_POINTER(result);
}
}
/*
* performs the difference of polygon1 - polygon2
* (parameters order)
* supports MULTIPOLYGON with only one polygon in it
* Returns null if input types are not polygons, or
* different SRIDS
* Returns a polygon or multipolygon according to result
*/
PG_FUNCTION_INFO_V1(difference_pg_pg);
Datum difference_pg_pg(PG_FUNCTION_ARGS) {
GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
GEOMETRY *result;
// the conversion of geometry to POLYGON3D
POLYGON3D *poly1;
POLYGON3D *poly2;
int32 *offset1;
int32 *offset2;
char *o1;
char *o2;
elog(LOG,"\n--------------------------------------------------------------");
// looks at the input geometries type
if (geom1->type != POLYGONTYPE || geom2->type != POLYGONTYPE) {
elog(ERROR,"Incompatible types for 2 geometries (expected POLYGON)");
PG_RETURN_NULL();
}
// should add support for multipolygon with 1 object in it
if (geom1->SRID != geom2->SRID) {
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs");
PG_RETURN_NULL();
}
// the basic optimizations:
// are their bbox disjoints ?
if (!box3d_ov(&(geom1->bvol), &(geom2->bvol))) {
// result is first geom: A-B = A
elog(LOG,"difference_pg_pg: Geometries bbox disjoint");
PG_RETURN_POINTER(geom1);
}
// creates the 2 polygons
offset1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
offset2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
o1 = (char *) geom1 +offset1[0];
o2 = (char *) geom2 +offset2[0];
poly1 = (POLYGON3D *)o1;
poly2 = (POLYGON3D *)o2;
// are geometries equals ?
if (is_same_polygon(poly1, poly2)) {
elog(LOG, "difference_pg_pg: same input polygons, returning first one");
PG_RETURN_NULL();
}
elog(LOG,"difference_pg_pg: input polygon 1: %s", geometry_to_text(geom1));
elog(LOG,"difference_pg_pg: input polygon 2: %s", geometry_to_text(geom2));
result = scot_gpc_clipping(poly1, poly2, GPC_DIFF, geom1->SRID);
if (result == NULL) {
// empty result set, not an error.
PG_RETURN_NULL();
} else {
result->SRID = geom1->SRID;
PG_RETURN_POINTER(result);
}
}
/**
* performs the symmetric difference of polygon1 - polygon2
* (parameters order)
* supports MULTIPOLYGON with only one polygon in it
* Returns null if input types are not polygons, or
* different SRIDS
* Returns a polygon or multipolygon according to result
*/
PG_FUNCTION_INFO_V1(sym_difference_pg_pg);
Datum sym_difference_pg_pg(PG_FUNCTION_ARGS) {
GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
GEOMETRY *result;
// the conversion of geometry to POLYGON3D
POLYGON3D *poly1;
POLYGON3D *poly2;
int32 *offset1;
int32 *offset2;
char *o1;
char *o2;
elog(LOG,"\n--------------------------------------------------------------");
// looks at the input geometries type
if (geom1->type != POLYGONTYPE || geom2->type != POLYGONTYPE) {
elog(ERROR,"sym_difference_pg_pg: Incompatible types for 2 geometries (expected POLYGON)");
PG_RETURN_NULL();
}
// should add support for multipolygon with 1 object in it
if (geom1->SRID != geom2->SRID) {
elog(ERROR,"sym_difference_pg_pg: Operation on two GEOMETRIES with different SRIDs");
PG_RETURN_NULL();
}
// should provide a function to add 2 input geom into the result in case of
// disjoint geometries.
// creates the 2 polygons
offset1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
offset2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
o1 = (char *) geom1 +offset1[0];
o2 = (char *) geom2 +offset2[0];
poly1 = (POLYGON3D *)o1;
poly2 = (POLYGON3D *)o2;
// are geometries equals ?
if (is_same_polygon(poly1, poly2)) {
// A U B = A or B if A=B
elog(LOG, "sym_difference_pg_pg: same input polygons, returning first one");
PG_RETURN_POINTER(geom1);
}
elog(LOG,"sym_difference_pg_pg: input polygon 1: %s", geometry_to_text(geom1));
elog(LOG,"sym_difference_pg_pg: input polygon 2: %s", geometry_to_text(geom2));
result = scot_gpc_clipping(poly1, poly2, GPC_XOR, geom1->SRID);
if (result == NULL) {
// empty result set, not an error.
PG_RETURN_NULL();
} else {
result->SRID = geom1->SRID;
PG_RETURN_POINTER(result);
}
}
/*************************************************
**********CONVERSION GPC <-> POSTGIS types********
**************************************************/
/**
* Performs the gpc clipping operation on the 2 given polygons
* and returns a pointer to resulting geometry (POLYGON or MULTIPOLYGON)
*/
GEOMETRY * scot_gpc_clipping(POLYGON3D *ppoly1, POLYGON3D *ppoly2, gpc_op operation, int SRID) {
// the gpc object on which real computation will be done
gpc_polygon gpcsubject, gpcclip, gpcresult;
// creates gpc polygon from postgis geometry:
POLYGON3D_to_gpc(ppoly1, &gpcsubject);
POLYGON3D_to_gpc(ppoly2, &gpcclip);
// perform the clipping operation (union geom1-geom2 for test)
gpc_polygon_clip(operation, &gpcsubject, &gpcclip, &gpcresult);
// converts and returns the result
return ( gpc_to_GEOMETRY(&gpcresult, SRID) );
}
/**
* Converts the given POLYGON3D object into a gpc polygon
*/
void POLYGON3D_to_gpc(POLYGON3D *ppoly, gpc_polygon *p) {
int c, v;
// is the current polygon ring a hole ?
int isHole;
// the position of the current point in the array of points for the poly
int point_offset;
// the array of points for each ring
POINT3D *pts;
// First ring is always exterior for OpenGIS polygon
isHole = 0;
// copy the number of rings of the input polygon
p->num_contours = ppoly->nrings;
#ifdef _DEBUG_
//elog(LOG,"POLYGON3D_to_gpc: num rings to convert: %d", ppoly->nrings);
#endif
// allocates rings
MALLOC(p->hole, p->num_contours * sizeof(int), "hole flag array creation");
MALLOC(p->contour, p->num_contours * sizeof(gpc_vertex_list), "contour creation");
// gets the array of points for all the rings forming the polygon:
pts = (POINT3D *) ( (char *)&(ppoly->npoints[ppoly->nrings] ) );
pts = (POINT3D *) MAXALIGN(pts);
point_offset = 0;
// loop for each ring, num_countours is the number of rings in ppoly
for (c= 0; c < p->num_contours; c++) {
// num vertices for this ring: 1 less than postgis polygon
p->contour[c].num_vertices = ppoly->npoints[c] - 1;
// set hole feature of this poly
p->hole[c] = isHole;
if (c == 0) {
// now first ring is treated, all other rings are holes
isHole = 1;
}
// allocates vertex for this ring
MALLOC(p->contour[c].vertex, p->contour[c].num_vertices * sizeof(gpc_vertex), "vertex creation");
// copies each vertex coordinates from opengis poly to gpc poly
for (v= 0; v < p->contour[c].num_vertices; v++) {
p->contour[c].vertex[v].x = pts[v+point_offset].x;
p->contour[c].vertex[v].y = pts[v+point_offset].y;
//elog(LOG,"POLYGON3D_to_gpc: vertex read from ring %d: %f %f", c, p->contour[c].vertex[v].x, p->contour[c].vertex[v].y);
}
//elog(LOG,"POLYGON3D_to_gpc: read %d vertices for ring %d", v, c);
// increments current ring offset: must go forward in the array of
// points to the next ring's first point,
// be careful that postgis polygons contain 1 more vertex than gpc polygon do
// point_offset is here the position inside pts array. The last opengis polygon
// vertex is discarded this way
point_offset += ppoly->npoints[c];
}
elog(LOG,"POLYGON3D_to_gpc: gpc_poly: %s", gpc_polygon_to_text(p));
}
/*
* Converts the given gpc polygon into a GEOMETRY.
* As gpc can represent a multipolygon, we cannot return
* a polygon3D object, to stay symmetric with POLYGON3D_to_gpc function.
* Should provide an efficient way to dynamically allocates memory for
* resulting polygon(s).
* CAUTION: For the moment, max number of resulting polygons is limited to 100,
* max number of holes per polygon is limited to 10. (+ 1 ext)
* MUST change that for a dynamic allocation mechanism.
*/
#define MAX_HOLES 10
#define MAX_RINGS 100
GEOMETRY* gpc_to_GEOMETRY(gpc_polygon *p, int SRID) {
GEOMETRY *result;
// the bbox of the result.
BOX3D *the_bbox;
POLYGON3D *tmp_obj;
// POLYGON3D *ppoly;
int tmp_obj_size=0;
// some loop counters
int c, v;
// number of resulting postgis polygons (exterior contours in the given gpc poly);
int num_polygons = 0;
// the array of polygons contained in the gpc polygon
// each line is an array of rings forming a postgis polygon.
// All values in this bi-d array is an index of p->contour[];
// The indices are sorted so that all holes belong to the right exterior.
// the postgis point_in_polygon function is used to determine that
// First entry in a line is the polygon's exterior. Other are holes
int polygons [MAX_RINGS] [MAX_HOLES+1];
// the array of polygon's size, ie the number of rings forming each polygon.
// this array is // to polygons
int num_rings[MAX_RINGS];
// the number of rings for the current polygon
int tmp_num_rings = 0;
// test input gpc emptiness
if (p->num_contours == 0) {
elog(LOG, "********gpc result polygon is empty****************");
return NULL;
}
//elog(LOG,"gpc_to_GEOMETRY: num countours2: %d", p->num_contours);
elog(LOG,"gpc_to_GEOMETRY: gpc_poly: %s", gpc_polygon_to_text(p));
// construct the polygons array according to p->contours parameters:
// is it a hole or not, if so, which exterior does it belong to ?
for (c = 0; c < p->num_contours; c++) {
if (p->hole[c] == 0) {
//elog(LOG,"gpc_to_GEOMETRY: find ext at %d", c);
// this is a new polygon:
num_polygons++;
// this is also a new ring for this current polygon
tmp_num_rings++;
// stores the index of this ring into our array
polygons[num_polygons-1][tmp_num_rings-1] = c;
// reloop for each interior, to see if it belongs to the polygon
for (v = 0; v < p->num_contours; v++) {
if (p->hole[v] == 1) {
//elog(LOG,"gpc_to_GEOMETRY: find int at %d", v);
// this is a hole, look if it belongs to the current exterior
if (gpc_hole_in_gpc_contour(p, c, v)) {
//elog(LOG,"gpc_to_GEOMETRY: Current contour (%d) is inside Exterior %d",v, c);
// a new ring for the current polygon
tmp_num_rings++;
// stores its index into our array
// contour at index v is inside ppoly: a new ring
polygons[num_polygons-1][tmp_num_rings-1] = v;
}
}
} // end for interiors
// all interiors for this current ring are treaded: store their count:
num_rings[num_polygons-1] = tmp_num_rings;
// can pass to next exterior, with 0 rings
tmp_num_rings = 0;
}
} // end for exteriors
elog(LOG,"gpc_to_GEOMETRY: after ordering: found %d polys", num_polygons);
elog(LOG,"gpc_to_GEOMETRY: final matrix:");
for (c = 0; c < num_polygons; c++) {
elog(LOG,"polygon %d: ", c);
for (v = 0; v < num_rings[c]; v++) {
elog(LOG,"\tring %d", polygons[c][v]);
}
}
//construct a geometry with first polygon,
tmp_obj = gpc_to_POLYGON3D(p, polygons[0], num_rings[0], SRID, &tmp_obj_size);
result = make_oneobj_geometry(tmp_obj_size,
(char *)tmp_obj,
POLYGONTYPE,
FALSE,
SRID,
1.0,
0.0,
0.0
);
//elog(LOG,"gpc_to_GEOMETRY: first resulting polygon built: %s", geometry_to_text(result));
// adds all other polygons into it
for (c = 1; c < num_polygons; c++) {
tmp_obj = gpc_to_POLYGON3D(p, polygons[c], num_rings[c], SRID, &tmp_obj_size);
result = add_to_geometry(result, tmp_obj_size, (char *)tmp_obj, POLYGONTYPE);
}
//computes the bbox
the_bbox = bbox_of_geometry(result);
if (the_bbox != NULL) {
memcpy( &result->bvol, the_bbox, sizeof(BOX3D) );
pfree(the_bbox);
}
elog(LOG,"gpc_to_GEOMETRY: resulting polygon %d built: %s", c, geometry_to_text(result));
return result;
}
/*
Builds a geometry from the given gpc_polygon,
rings is the array of indices into p->contour[] forming the polygon
MUST control the contour orientation: if it is cw or ccw.
Have to add vertices in good order to opengis ccw polygon orientation
This code is more or less a copy of the one in gpc_to_POLYGON3D:
Should change that to merge codes.
*/
POLYGON3D * gpc_to_POLYGON3D(gpc_polygon *p, int *rings, int num_rings, int SRID, int *ppoly_size) {
POLYGON3D *ppoly;
// The array of points: total number of points for all rings
POINT3D *pts;
int c, v, ci;
int num_vertices = 0;
// array of points per ring
int *pts_per_ring;
// the total number of points, used to set pts position correctly
int point_offset = 0;
// first, creates a new polygon from scratch.
pts_per_ring = malloc(num_rings * (sizeof(int)));
// loop through all given rings to get geometry's num points
for (c = 0; c < num_rings; c++) {
// contour's index
ci = rings[c];
num_vertices += p->contour[ci].num_vertices + 1; /* reserve an extra point */
pts_per_ring[c] = p->contour[ci].num_vertices + 1;
}
pts = malloc(num_vertices * sizeof(POINT3D));
//loop through all contours to build the polygon
for (c = 0; c < num_rings; c++) {
// contour's index
ci = rings[c];
//elog(LOG,"gpc_to_POLYGON3D: treating gpc contour %d", ci);
// for each vertex in the current contour:
// a contour has at least 3 vertices
if (is_ccw(p->contour[ci].vertex[0].x,
p->contour[ci].vertex[0].y,
p->contour[ci].vertex[1].x,
p->contour[ci].vertex[1].y,
p->contour[ci].vertex[2].x,
p->contour[ci].vertex[2].y)) {
//elog(LOG,"gpc_to_POLYGON3D: ring: %d is ccw: add coordinates normally", c);
// gpc contour is ccw, can add its vertices normally
for (v = 0; v < p->contour[ci].num_vertices; v++) {
set_point( &pts[v + point_offset],
p->contour[ci].vertex[v].x,
p->contour[ci].vertex[v].y, 0.0);
}
// adds the first point at the end to close the postgis ring
set_point( &pts[v+point_offset],
p->contour[ci].vertex[0].x,
p->contour[ci].vertex[0].y,
0.0);
} else {
//elog(LOG,"gpc_to_POLYGON3D: ring: %d is cw: add coordinates in reverse order", c);
// gpc contour is cw, must loop back to construct a postgis poly
// with a good orientation
for (v = p->contour[ci].num_vertices-1; v >= 0 ; v--) {
set_point( &pts[(p->contour[ci].num_vertices-v-1)+point_offset],
p->contour[ci].vertex[v].x,
p->contour[ci].vertex[v].y, 0.0);
}
// adds the first point at the end to close the postgis ring
set_point( &pts[p->contour[ci].num_vertices+point_offset],
p->contour[ci].vertex[p->contour[ci].num_vertices-1].x,
p->contour[ci].vertex[p->contour[ci].num_vertices-1].y,
0.0);
}
// increments the point offset= total number of vertices since here
point_offset += p->contour[ci].num_vertices + 1;
} // end loop for all contours
//make the polygon
elog(LOG,"gpc_to_POLYGON3D: making final polygon for %d rings, %d vertices",num_rings, num_vertices);
ppoly = make_polygon(num_rings, pts_per_ring, pts, num_vertices, ppoly_size);
return ppoly;
}
/*
Returns geom1 in which geom2 was added.
*/
GEOMETRY* add_geometry(GEOMETRY *geom1, GEOMETRY *geom2) {
BOX3D *the_bbox;
geom1 = make_oneobj_geometry(tmp_obj_size,
(char *)tmp_obj,
POLYGONTYPE,
FALSE,
SRID,
1.0,
0.0,
0.0
);
}
/*
* Returns true if x0y0, x1y1, x2y2 is ccw (trigo orientation)
* computes the cross product of the 3 points
*/
bool is_ccw(double x0, double y0, double x1, double y1, double x2, double y2) {
return (x0*y1 + x2*y0 + x1*y2 - x2*y1 - x0*y2 - x1*y0) >= 0;
}
/*
* Returns true if first point of given gpc polygon contour lies
* in given polygon3d.
* i is the contour index into c.
* as a opengis polygon hole may touch the exterior, performs the test
* on 2 consecutive points in case of negative result.
* (2 consecutive points cannot lie on the exterior)
*/
bool contour_in_POLYGON3D(POLYGON3D *ppoly, gpc_polygon *p, int i) {
POINT3D pt;
pt.x = p->contour[i].vertex[0].x;
pt.y = p->contour[i].vertex[0].y;
pt.z = 0.0;
if (point_in_poly(&pt, ppoly)) {
return TRUE;
}
// test with next point
pt.x = p->contour[i].vertex[1].x;
pt.y = p->contour[i].vertex[1].y;
pt.z = 0.0;
return point_in_poly(&pt, ppoly);
}
/*
* Return true if the given gpc_hole is inside the given gpc_contour:
* p the gpc_polygon containing all holes and contour.
* c, the index of exterior contour to search in, in the p->contour array.
* v, the index of interior contour to search for, in the p->contour array.
* Assuming that 2 consecutive vertices cannot be on the exterior, searches
* for the 2 first interior vertices to see if it is inside.
* this assumption has to be confirmed by examples.
*/
bool gpc_hole_in_gpc_contour(gpc_polygon *p, int c, int v) {
#ifdef _DEBUG_
//elog(ERROR, "gpc_hole_in_gpc_contour");
#endif
if (point_in_gpc_poly(p->contour[v].vertex,
p->contour[c].vertex,
p->contour[c].num_vertices) == 1) {
return true;
}
return (point_in_gpc_poly(p->contour[v].vertex,
p->contour[c].vertex,
p->contour[c].num_vertices) == 1);
}
// copied from postgis_fn.c:
// adapted to work with gpc types:
// cn_PnPoly(): crossing number test for a point in a polygon
// input: P = a vertex,
// V[] = vertex points of a contour V[n] with V[n]!=V[0]
// returns: 0 = outside, 1 = inside
//
// Our polygons DOES NOT have first and last point the same:
// Must take first point when reaching the end of the ring
//
int point_in_gpc_poly( gpc_vertex *P, gpc_vertex *V, int n ) {
int cn = 0; // the crossing number counter
int i;
double vt;
gpc_vertex last_point;
//elog(ERROR, "point_in_gpc_poly: testing point [%g %g]", P->x, P->y);
// loop through all edges of the polygon
for (i=0; i< n; i++) {
// compute last point:
if ( i == n-1 ) {
// last point is first one, last segment is treated
last_point = V[0];
} else {
last_point = V[i+1];
}
if (((V[i].y <= P->y) && (last_point.y > P->y)) // an upward crossing
|| ((V[i].y > P->y) && (last_point.y <= P->y))) { // a downward crossing
vt = (double)(P->y - V[i].y) / (last_point.y - V[i].y);
if (P->x < V[i].x + vt * (last_point.x - V[i].x)) {
// P.x <intersect
++cn; // a valid crossing of y=P.y right of P.x
}
}
}
return (cn&1); // 0 if even (out), and 1 if odd (in)
}
/*
To debug gpc_polygon on screen.
From postgis_inout.c:
*/
#define SHOW_DIGS_DOUBLE 15
#define MAX_DIGS_DOUBLE (SHOW_DIGS_DOUBLE + 4 + 1 + 2 +1)
char* gpc_polygon_to_text(gpc_polygon *p) {
int size = 12;// "gpc_polygon" + null
char *res;
int c = 0, v;
char temp[MAX_DIGS_DOUBLE*2 + 4];
res = (char *) palloc(size);
strcpy(res, "gpc_polygon");
for (c=0; c < p->num_contours; c++) {
// reallocates for this new contour:
size += (p->contour[c].num_vertices * (2*MAX_DIGS_DOUBLE +1+1)) + 3 + 2; //(int,x y,x y)
res = (char *)repalloc(res, size);
if (p->hole[c]) {
// interior
strcat(res, "(int,");
} else {
strcat(res, "(ext,");
}
for (v=0; v < p->contour[c].num_vertices; v++) {
sprintf(temp, "%.15g %.15g", p->contour[c].vertex[v].x, p->contour[c].vertex[v].y);
strcat(res, temp);
strcat(res, ",");
}
// closes this contour
strcat(res, ")");
}
return res;
}
/*************************************************
***************TEST ZONE**************************
**************************************************/
/*
tests
very first gpc operation on 2 geometries
Works only on 2 real POLYGON object,
*/
PG_FUNCTION_INFO_V1(test_geogpc);
Datum test_geogpc(PG_FUNCTION_ARGS) {
GEOMETRY *geom1 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
GEOMETRY *geom2 = (GEOMETRY *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
GEOMETRY *result = (GEOMETRY *) palloc(sizeof(GEOMETRY) );
POLYGON3D *poly1;
POLYGON3D *poly2;
POLYGON3D *poly_res;
int32 *offset1;
int32 *offset2;
int poly_res_size;
char *o1;
char *o2;
// for debug only
/* gpc types */
gpc_polygon gpcsubject, gpcclip, gpcresult;
// looks at the input geometries type
if (geom1->type != POLYGONTYPE || geom2->type != POLYGONTYPE) {
elog(ERROR,"Incompatible types for 2 geometries (expected POLYGON)");
PG_RETURN_NULL();
}
#ifdef _DEBUG_
elog(LOG,"geom type is polygon ");
#endif
if (geom1->SRID != geom2->SRID) {
elog(ERROR,"Operation on two GEOMETRIES with different SRIDs");
PG_RETURN_NULL();
}
#ifdef _DEBUG_
elog(LOG,"SRID's are equals ");
#endif
// creates the 2 polygons and 2 gpc polygons on which operations will be done
offset1 = (int32 *) ( ((char *) &(geom1->objType[0] ))+ sizeof(int32) * geom1->nobjs );
offset2 = (int32 *) ( ((char *) &(geom2->objType[0] ))+ sizeof(int32) * geom2->nobjs );
o1 = (char *) geom1 +offset1[0];
o2 = (char *) geom2 +offset2[0];
poly1 = (POLYGON3D *)o1;
poly2 = (POLYGON3D *)o2;
// creates gpc polygon from postgis polygons:
scot_gpc_read_polygon(poly1, &gpcsubject);
scot_gpc_read_polygon(poly2, &gpcclip);
// perform the clipping operation (union geom1-geom2 for test)
gpc_polygon_clip(GPC_UNION, &gpcsubject, &gpcclip, &gpcresult);
//Convert result to opengis polygon:
poly_res = scot_gpc_write_polygon(&gpcresult, &poly_res_size);
result = make_oneobj_geometry(poly_res_size,
(char *)poly_res,
POLYGONTYPE,
FALSE,
geom1->SRID,
geom1->scale,
geom1->offsetX,
geom1->offsetY
);
#ifdef _DEBUG_
elog(LOG,"test_geogpc: WKT representation of result: %s\n", geometry_to_text(result));
#endif
// free memory for gpc polygons.
gpc_free_polygon(&gpcsubject);
gpc_free_polygon(&gpcclip);
gpc_free_polygon(&gpcresult);
// free memory for postgis polygons ?
PG_FREE_IF_COPY(geom1, 0);
PG_FREE_IF_COPY(geom2, 1);
PG_RETURN_POINTER(result);
}
/*
add-on to GPC to allow the construction of a polygon from an array of coordinates instead of file
method should allow the construction of a GPC polygon from a postgis POLYGON or MULTIPOLYGON
// original method from gpc.c
void gpc_read_polygon(FILE *fp, int read_hole_flags, gpc_polygon *p) {
int c, v;
fscanf(fp, "%d", &(p->num_contours));
MALLOC(p->hole, p->num_contours * sizeof(int),
"hole flag array creation");
MALLOC(p->contour, p->num_contours
* sizeof(gpc_vertex_list), "contour creation");
for (c= 0; c < p->num_contours; c++)
{
fscanf(fp, "%d", &(p->contour[c].num_vertices));
if (read_hole_flags)
fscanf(fp, "%d", &(p->hole[c]));
else
p->hole[c]= FALSE;
MALLOC(p->contour[c].vertex, p->contour[c].num_vertices
* sizeof(gpc_vertex), "vertex creation");
for (v= 0; v < p->contour[c].num_vertices; v++)
fscanf(fp, "%lf %lf", &(p->contour[c].vertex[v].x),
&(p->contour[c].vertex[v].y));
}
}
No control is made on the input polygon3d geometry. Use this function with caution
Convention: when referening to gpc polygons, contour is the ring
for opengis poly, ring is a contour, first one is exterior, all other are holes.
CAUTION: postgis polygons have their last ring point = first one, must deal with this
when building gpc polygons.
*/
void scot_gpc_read_polygon(POLYGON3D *ppoly, gpc_polygon *p) {
int c, v;
// is the current polygon ring a hole ?
int isHole;
// the position of the current point in the array of points for the poly
int point_offset;
// the array of points for each ring
POINT3D *pts;
// First ring is always exterior for OpenGIS polygon
isHole = 0;
// copy the number of rings of the input polygon
p->num_contours = ppoly->nrings;
#ifdef _DEBUG_
elog(LOG,"gpc_read_polygon: num rings to convert: %d", ppoly->nrings);
#endif
// allocates rings
MALLOC(p->hole, p->num_contours * sizeof(int), "hole flag array creation");
MALLOC(p->contour, p->num_contours * sizeof(gpc_vertex_list), "contour creation");
// gets the array of points for all the rings forming the polygon:
pts = (POINT3D *) ( (char *)&(ppoly->npoints[ppoly->nrings] ) );
pts = (POINT3D *) MAXALIGN(pts);
point_offset = 0;
// loop for each ring, num_countours is the number of rings in ppoly
for (c= 0; c < p->num_contours; c++) {
// num vertices for this ring: 1 less thant postgis polygon
p->contour[c].num_vertices = ppoly->npoints[c] - 1;
// set hole feature of this poly
p->hole[c] = isHole;
if (c == 0) {
// now first ring is treated, all other rings are holes
isHole = 1;
}
// allocates vertex for this ring
MALLOC(p->contour[c].vertex, p->contour[c].num_vertices * sizeof(gpc_vertex), "vertex creation");
// copies each vertex coordinates from opengis poly to gpc poly
for (v= 0; v < p->contour[c].num_vertices; v++) {
p->contour[c].vertex[v].x = pts[v].x;
p->contour[c].vertex[v].y = pts[v].y;
#ifdef _DEBUG_
elog(LOG,"gpc_read_polygon: vertex read from ring %d: %f %f", c, p->contour[c].vertex[v].x, p->contour[c].vertex[v].y);
#endif
}
#ifdef _DEBUG_
elog(LOG,"gpc_read_polygon: read %d vertices for ring %d", v, c);
#endif
//increments current ring offset: must go forward in the array of points to the next ring's first point,
// be careful that postgis polygons contain 1 more vertex than gpc polygon do
point_offset += ppoly->npoints[c];
}
}
/*
Write function, from GPC code
// original method from gpc.c
void gpc_write_polygon(FILE *fp, int write_hole_flags, gpc_polygon *p)
{
int c, v;
fprintf(fp, "%d\n", p->num_contours);
for (c= 0; c < p->num_contours; c++)
{
fprintf(fp, "%d\n", p->contour[c].num_vertices);
if (write_hole_flags)
fprintf(fp, "%d\n", p->hole[c]);
for (v= 0; v < p->contour[c].num_vertices; v++)
fprintf(fp, "% .*lf % .*lf\n",
DBL_DIG, p->contour[c].vertex[v].x,
DBL_DIG, p->contour[c].vertex[v].y);
}
}
Writes an existing gpc polygon into a newly created postgis polygon
No control is made on the input gpc polygon geometry. Use this function with caution
Convention: when referening to gpc polygons, contour is the ring
for opengis poly, ring is a contour, first one is exterior, all other are holes.
caution: the number of vertices in each gpc polygon is 1 less than the number in
postgis polygon, as OpenGIS says a ring must have its last point=first point.
Force this point in this code.
*/
POLYGON3D * scot_gpc_write_polygon(gpc_polygon *p, int *ppoly_size) {
int c, v;
// The array of points: total number of points for all rings
POINT3D *pts;
POLYGON3D *ppoly; /*the result*/
int num_vertices = 0;
// array of points per ring
int *pts_per_ring;
#ifdef _DEBUG_
elog(LOG,"gpc_write_polygon: num rings to write: %d", p->num_contours);
#endif
// first, creates a new polygon from scratch.
pts_per_ring = malloc(p->num_contours * sizeof(int));
for (c = 0; c < p->num_contours; c++) {
num_vertices += p->contour[c].num_vertices + 1; /* reserve an extra point */
pts_per_ring[c] = p->contour[c].num_vertices + 1;
}
pts = malloc(num_vertices * sizeof(POINT3D));
//loop for all contours to take their vertices
// no holes to take into account as postgis will treat them correctly
// adds the first point at the end to close the postgis ring
for (c= 0; c < p->num_contours; c++) {
#ifdef _DEBUG_
elog(LOG,"gpc_write_polygon: start of the loop: %d", p->contour[c].num_vertices);
#endif
for (v = p->contour[c].num_vertices-1; v >= 0 ; v--) {
set_point( &pts[c+(p->contour[c].num_vertices-v-1)], p->contour[c].vertex[v].x, p->contour[c].vertex[v].y, 0.0);
#ifdef _DEBUG_
elog(LOG,"gpc_write_polygon: coord (%d) written: %f %f",v, p->contour[c].vertex[v].x, p->contour[c].vertex[v].y);
#endif
}
// adds the first point at the end to close the postgis ring
set_point( &pts[c + p->contour[c].num_vertices], p->contour[c].vertex[p->contour[c].num_vertices-1].x, p->contour[c].vertex[p->contour[c].num_vertices-1].y, 0.0);
#ifdef _DEBUG_
elog(LOG,"gpc_write_polygon: LAST coord written: %f %f", p->contour[c].vertex[p->contour[c].num_vertices-1].x, p->contour[c].vertex[p->contour[c].num_vertices-1].y);
#endif
#ifdef _DEBUG_
elog(LOG,"gpc_write_polygon: written %d vertices for ring: %d", v+1, c);
#endif
}
//make the polygon
ppoly = make_polygon(p->num_contours, pts_per_ring, pts, num_vertices, ppoly_size);
return ppoly;
}
/* working code backup
for (c= 0; c < p->num_contours; c++) {
for (v = 0; v < p->contour[c].num_vertices; v++) {
set_point( &pts[c+v], p->contour[c].vertex[v].x, p->contour[c].vertex[v].y, 0.0);
}
// adds the first point at the end to close the postgis ring
set_point( &pts[c+v+1], p->contour[c].vertex[0].x, p->contour[c].vertex[0].y, 0.0);
}
*/
More information about the postgis-users
mailing list