From: Sandro Santilli Date: Mon, 27 Oct 2003 10:21:15 +0000 (+0000) Subject: Initial import. X-Git-Tag: pgis_0_8_0~53 X-Git-Url: https://granicus.if.org/sourcecode?a=commitdiff_plain;h=587b27d593a102f6b334efec4b6f3d3928e5b7f2;p=postgis Initial import. git-svn-id: http://svn.osgeo.org/postgis/trunk@329 b70326c6-7e19-0410-871a-916f4a2858ee --- diff --git a/postgis_algo.c b/postgis_algo.c new file mode 100644 index 000000000..322b2b57c --- /dev/null +++ b/postgis_algo.c @@ -0,0 +1,409 @@ +/********************************************************************** + * $Id$ + * + * PostGIS - Spatial Types for PostgreSQL + * http://postgis.refractions.net + * Copyright 2001-2003 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. + * + ********************************************************************** + * $Log$ + * Revision 1.1 2003/10/27 10:21:15 strk + * Initial import. + * + **********************************************************************/ + +#include "postgres.h" +#include "postgis.h" + + +/*********************************************************************** + * Simple Douglas-Peucker line simplification. + * No checks are done to avoid introduction of self-intersections. + * No topology relations are considered. + * + * --strk@keybit.net; + ***********************************************************************/ + +#define SAMEPOINT(a,b) ((a)->x==(b)->x&&(a)->y==(b)->y&&(a)->z==(b)->z) +#define VERBOSE 0 + +#if VERBOSE > 0 +#undef REPORT_POINTS_REDUCTION +#undef REPORT_RINGS_REDUCTION +#define REPORT_RINGS_ADJUSTMENTS +#endif + +#undef CHECK_RING_IS_CLOSE + + +/* + * Search farthest point from segment p1-p2 + * returns squared distance an int pointer + */ +void DP_findsplit(POINT3D *pts, int npts, int p1, int p2, + int *split, double *dist) +{ + int k; + POINT3D *pa, *pb, *pk; + double tmp; + +#if VERBOSE > 4 +elog(NOTICE, "DP_findsplit called"); +#endif + + *dist = -1; + *split = p1; + + if (p1 + 1 < p2) + { + + pa = &(pts[p1]); + pb = &(pts[p2]); + +#if VERBOSE > 4 +elog(NOTICE, "DP_findsplit: P%d(%f,%f) to P%d(%f,%f)", + p1, pa->x, pa->y, p2, pb->x, pb->y); +#endif + + for (k=p1+1; k 4 +elog(NOTICE, "DP_findsplit: P%d(%f,%f)", k, pk->x, pk->y); +#endif + + /* distance computation */ + tmp = distance_pt_seg(pk, pa, pb); + + if (tmp > *dist) + { + *dist = tmp; /* record the maximum */ + *split = k; +#if VERBOSE > 4 +elog(NOTICE, "DP_findsplit: P%d is farthest (%g)", k, *dist); +#endif + } + } + + } /* length---should be redone if can == 0 */ + + else + { +#if VERBOSE > 3 +elog(NOTICE, "DP_findsplit: segment too short, no split/no dist"); +#endif + } + +} + + +void +DP_simplify(POINT3D *inpts, int inptsn, POINT3D **outpts, int *outptsn, double epsilon) +{ + int stack[inptsn]; /* recursion stack */ + int sp=-1; /* recursion stack pointer */ + int p1, split; + double dist_sq; + double epsilon_sq = epsilon*epsilon; + POINT3D *outpoints; + int numoutpoints=0; + + p1 = 0; + stack[++sp] = inptsn-1; + +#if VERBOSE > 4 + elog(NOTICE, "DP_simplify called"); +#endif + + outpoints = (POINT3D *)palloc( sizeof(POINT3D) * inptsn); + memcpy(outpoints, inpts, sizeof(POINT3D)); + numoutpoints++; + +#if VERBOSE > 3 +elog(NOTICE, "DP_simplify: added P0 to simplified point array (size 1)"); +#endif + + + do + { + + DP_findsplit(inpts, inptsn, p1, stack[sp], &split, &dist_sq); +#if VERBOSE > 3 + elog(NOTICE, "DP_simplify: farthest point from P%d-P%d is P%d (dist. %g)", p1, stack[sp], split, sqrt(dist_sq)); +#endif + + if (dist_sq > epsilon_sq) { + stack[++sp] = split; + } else { + /* + outpoints = repalloc( outpoints, sizeof(POINT3D) * numoutpoints+1 ); + if ( outpoints == NULL ) elog(NOTICE, "Out of virtual memory"); + */ + memcpy(outpoints+numoutpoints, &(inpts[stack[sp]]), + sizeof(POINT3D)); + numoutpoints++; +#if VERBOSE > 3 +elog(NOTICE, "DP_simplify: added P%d to simplified point array (size: %d)", + stack[sp], numoutpoints); +#endif + p1 = stack[sp--]; + } +#if VERBOSE > 5 +elog(NOTICE, "stack pointer = %d", sp); +#endif + } + while (! (sp<0) ); + + /* + * If we have reduced the number of points realloc + * outpoints array to free up some memory. + * Might be turned on and of with a SAVE_MEMORY define ... + */ + if ( numoutpoints < inptsn ) + { + outpoints = (POINT3D *)repalloc(outpoints, sizeof(POINT3D)*numoutpoints); + if ( outpoints == NULL ) { + elog(ERROR, "Out of virtual memory"); + } + } + + *outpts = outpoints; + *outptsn = numoutpoints; +} + +char * +simplify_line3d(LINE3D *iline, double dist) +{ + POINT3D *ipts; + POINT3D *opts; + int iptsn, optsn, olinesize; + LINE3D *oline; + + ipts = iline->points; + iptsn = iline->npoints; + + DP_simplify(ipts, iptsn, &opts, &optsn, dist); + + oline = make_line(optsn, opts, &olinesize); + + return (char *)oline; +} + +char * +simplify_polygon3d(POLYGON3D *ipoly, double dist) +{ + POINT3D *ipts; + POINT3D *opts; + int iptsn, optsn, size; + int nrings; + int pts_per_ring[ipoly->nrings]; + POLYGON3D *opoly; + int ri; + POINT3D *allpts = NULL; + int allptsn = 0; + int pt_off = 0; /* point offset for each ring */ + + nrings = 0; + +#ifdef REPORT_RINGS_REDUCTION +elog(NOTICE, "simplify_polygon3d: simplifying polygon with %d rings", + ipoly->nrings); +#endif + + /* Points start here */ + ipts = (POINT3D *) ((char *)&(ipoly->npoints[ipoly->nrings] )); + + pt_off=0; + for (ri=0; rinrings; pt_off += ipoly->npoints[ri++]) + { + iptsn = ipoly->npoints[ri]; + +#ifdef CHECK_RING_IS_CLOSE + if ( ! SAMEPOINT(ipts+pt_off, ipts+pt_off+iptsn-1) ) + elog(NOTICE, "First point != Last point"); +#endif + + DP_simplify(ipts+pt_off, iptsn, &opts, &optsn, dist); + if ( optsn < 2 ) + { + /* There as to be an error in DP_simplify */ + elog(NOTICE, "DP_simplify returned a <2 pts array"); + pfree(opts); + continue; + } + +#ifdef CHECK_RING_IS_CLOSE + if ( ! SAMEPOINT(opts, opts+optsn-1) ) + elog(NOTICE, "First point != Last point"); +#endif + + + if ( optsn < 4 ) + { + pfree(opts); +#ifdef REPORT_RINGS_ADJUSTMENTS + elog(NOTICE, "simplify_polygon3d: ring%d skipped ( <4 pts )", ri); +#endif + if ( ri ) continue; + else break; + } + + +#ifdef REPORT_POINTS_REDUCTION +elog(NOTICE, "simplify_polygon3d: ring%d simplified from %d to %d points", ri, iptsn, optsn); +#endif + + + /* + * Add ring to simplified ring array + * (TODO: dinamic allocation of pts_per_ring) + */ + pts_per_ring[nrings++] = optsn; + if ( ! allptsn ) { + allptsn = optsn; + allpts = palloc(sizeof(POINT3D)*allptsn); + memcpy(allpts, opts, sizeof(POINT3D)*optsn); + } else { + allptsn += optsn; + allpts = repalloc(allpts, sizeof(POINT3D)*allptsn); + memcpy(allpts+(allptsn-optsn), opts, sizeof(POINT3D)*optsn); + } + pfree(opts); + + if ( ! allpts ) { + elog(NOTICE, "Error allocating memory for all ring points"); + return NULL; + } + + } + +#ifdef REPORT_RINGS_REDUCTION +elog(NOTICE, "simplify_polygon3d: simplified polygon with %d rings", nrings); +#endif + + if ( nrings ) + { + opoly = make_polygon(nrings, pts_per_ring, allpts, allptsn, &size); + pfree(allpts); + return (char *)opoly; + } + else + { + return NULL; + } + +} + +char * +simplify_point3d(POINT3D *ipoint, double dist) +{ + return (char *)ipoint; +} + +PG_FUNCTION_INFO_V1(simplify); +Datum simplify(PG_FUNCTION_ARGS) +{ + Datum datum; + BOX3D *bbox; + GEOMETRY *orig_geom; + GEOMETRY *simp_geom = NULL; + char *orig_obj; /* pointer to each object in orig_geom */ + char *simp_obj; /* pointer to each simplified object */ + int simp_obj_size; /* size of simplified object */ + int32 *offsets; + int i; + double dist; + + if ( PG_ARGISNULL(0) ) PG_RETURN_NULL(); + datum = PG_GETARG_DATUM(0); + orig_geom = (GEOMETRY *)PG_DETOAST_DATUM(datum); + + if ( PG_ARGISNULL(1) ) PG_RETURN_NULL(); + dist = PG_GETARG_FLOAT8(1); + + /* + * Three or more points on a straight line will still collapse! + */ + // if ( dist == 0 ) PG_RETURN_POINTER(orig_geom); + + offsets = (int32 *) ( ((char *) &(orig_geom->objType[0] )) + + sizeof(int32) * orig_geom->nobjs ); + + + /* + * Simplify each subobject indipendently. + * No topology relations kept. + */ + for(i=0;inobjs; i++) + { + orig_obj = (char *) orig_geom+offsets[i]; + + if ( orig_geom->objType[i] == LINETYPE ) + { + simp_obj = simplify_line3d((LINE3D *)orig_obj, dist); + } + else if ( orig_geom->objType[i] == POLYGONTYPE ) + { + simp_obj = simplify_polygon3d((POLYGON3D *)orig_obj, dist); + } + else if ( orig_geom->objType[i] == POINTTYPE ) + { + simp_obj = simplify_point3d((POINT3D *)orig_obj, dist); + } + else + { + elog(NOTICE, "Unknown geometry type"); + PG_RETURN_NULL(); + } + + /* Simplified object degenerated to empty set */ + if ( ! simp_obj ) continue; + + /* Get size of simplified object */ + simp_obj_size = size_subobject(simp_obj, orig_geom->objType[i]); + + /* Create one-object geometry (will add objects later) */ + if ( simp_geom == NULL ) + { + simp_geom = make_oneobj_geometry( + simp_obj_size, simp_obj, orig_geom->objType[i], + orig_geom->is3d, orig_geom->SRID, orig_geom->scale, + orig_geom->offsetX, orig_geom->offsetY); + } + + /* Add object to already initialized geometry */ + else + { + simp_geom = add_to_geometry( + simp_geom, simp_obj_size, simp_obj, + orig_geom->objType[i] ); + } + + /* Error in simplified geometry construction */ + if ( ! simp_geom ) + { + elog(ERROR, "geometry construction failed at iteration %d", i); + PG_RETURN_NULL(); + } + + } + + if ( simp_geom == NULL ) PG_RETURN_NULL(); + + /* Calculate the bounding box */ + bbox = bbox_of_geometry(simp_geom); + memcpy(&(simp_geom->bvol), bbox, sizeof(BOX3D)); + pfree(bbox); + + + simp_geom->type = orig_geom->type; + PG_RETURN_POINTER(simp_geom); +} + +/*********************************************************************** + * --strk@keybit.net; + ***********************************************************************/