/* geo-ext.c * 2D Geometric extentions to Postgres developed by TNO, The Netherlands * derived on the standard Postgres "geo-ops.c" implementation. */ #include #include "utils/log.h" #include "utils/geo-decls.h" #include "geo-ext.h" #define LDELIM '(' #define RDELIM ')' #define DELIM ',' #define MDELIM ':' #define POINT2NARGS 2 #define BOX2NARGS 4 #define LSEG2NARGS 4 #define POLY2NBYTES(N) \ ((unsigned) (sizeof(POLYSTRUCT) + \ (((N)-1) > 0 ? ((N)-1) : 0) \ * sizeof(POINT2))) #define PATHALLOCSIZE(N) \ (long) ((unsigned) (sizeof(PATH) + \ (((N)-1) > 0 ? ((N)-1) : 0) \ * sizeof(POINT))) /* * POINT2 conversion functions. */ POINT2 * Point2_in(str) char *str; { extern double atof(); char *coord[POINT2NARGS], *p; int i; POINT2 *result; if (str == NULL) { elog(WARN, "Empty string POINT2"); return(NULL); } for (i = 0, p = str; *p && i < POINT2NARGS && *p != RDELIM; p++) if (*p == DELIM || (*p == LDELIM && !i)) coord[i++] = p + 1; if (i < POINT2NARGS - 1) { elog(WARN, "Ill formatted POINT2: %s", str); return(NULL); } result = PALLOCTYPE(POINT2); result->x = (float) atof(coord[0]); result->y = (float) atof(coord[1]); return(result); } char * Point2_out(pt) POINT2 *pt; { char *result; if (pt == NULL) { elog(WARN, "NULL pointer POINT2"); return(NULL); } result = (char *) palloc(40); (void) sprintf(result, "(%g,%g)", pt->x, pt->y); return(result); } POINT2 * To2Pnt(pt) POINT *pt; { POINT2 *result; if (pt == NULL) { elog(WARN, "NULL pointer POINT"); return(NULL); } result = PALLOCTYPE(POINT2); result->x = (float) pt->x; result->y = (float) pt->y; return(result); } POINT * From2Pnt(pt) POINT2 *pt; { POINT *result; if (pt == NULL) { elog(WARN, "NULL pointer POINT2"); return(NULL); } result = PALLOCTYPE(POINT); result->x = (double) pt->x; result->y = (double) pt->y; return(result); } BOX2 * Box2Pnt(pt) POINT2 *pt; { BOX2 *result; if (pt == NULL) { elog(WARN, "NULL pointer POINT2"); return(NULL); } result = PALLOCTYPE(BOX2); result->xl = pt->x; result->yl = pt->y; result->xh = pt->x; result->yh = pt->y; return(result); } BOX2 * Box2_in(str) char *str; { char *coord[BOX2NARGS], *p; int i; BOX2 *result; extern double atof(); float tmp; if (str == NULL) { elog(WARN, "Empty string BOX2"); return(NULL); } for (i = 0, p = str; *p && i < BOX2NARGS && *p != RDELIM; p++) if (*p == DELIM || (*p == LDELIM && !i)) coord[i++] = p + 1; if (i < BOX2NARGS - 1) { elog(WARN, "Ill formatted BOX2: %s", str); return(NULL); } result = PALLOCTYPE(BOX2); result->xh = (float) atof(coord[0]); result->yh = (float) atof(coord[1]); result->xl = (float) atof(coord[2]); result->yl = (float) atof(coord[3]); if (result->xh < result->xl) { tmp = result->xh; result->xh = result->xl; result->xl = tmp; } if (result->yh < result->yl) { tmp = result->yh; result->yh = result->yl; result->yl = tmp; } return(result); } char * Box2_out(box) BOX2 *box; { char *result; if (box == NULL) { elog(WARN, "NULL pointer for BOX2"); return(NULL); } result = (char *) palloc(80); (void) sprintf(result, "(%g,%g,%g,%g)", box->xl, box->yl, box->xh, box->yh); return(result); } /* Box2 functions*/ /*===============*/ BOX * From2Box(box) BOX2 *box; { BOX *result; if (box == NULL) return(NULL); result = PALLOCTYPE(BOX); result->xl = (double) box->xl; result->yl = (double) box->yl; result->xh = (double) box->xh; result->yh = (double) box->yh; return(result); } BOX2 * To2Box(box) BOX *box; { BOX2 *result; if (box == NULL) return(NULL); result = PALLOCTYPE(BOX2); result->xl = (float) box->xl; result->yl = (float) box->yl; result->xh = (float) box->xh; result->yh = (float) box->yh; return(result); } POINT2 * Center2Box(box) BOX2 *box; { POINT2 *result; result = PALLOCTYPE(POINT2); result->x = (float)((box->xh + box->xl) / 2.0); result->y = (float)((box->yh + box->yl) / 2.0); return(result); } double * Area2Box(box) BOX2 *box; { double *result; result = PALLOCTYPE(double); *result = (double) ((box->xh - box->xl) * (box->yh - box->yl)); return(result); } BOX2 * Inter2BoxBox(box1, box2) BOX2 *box1, *box2; { BOX2 *result; long Overlap2BoxBox(); if (! Overlap2BoxBox(box1,box2)) { /* elog(WARN, "No intersection BOX2 BOX3"); */ return(NULL); } result = PALLOCTYPE(BOX2); result->xh = (double) MIN(box1->xh, box2->xh); result->xl = (double) MAX(box1->xl, box2->xl); result->yh = (double) MIN(box1->yh, box2->yh); result->yl = (double) MAX(box1->yl, box2->yl); return(result); } BOX2 * Union2BoxBox(box1, box2) BOX2 *box1, *box2; { BOX2 *result; result = PALLOCTYPE(BOX2); result->xh = (double) MAX(box1->xh, box2->xh); result->xl = (double) MIN(box1->xl, box2->xl); result->yh = (double) MAX(box1->yh, box2->yh); result->yl = (double) MIN(box1->yl, box2->yl); return(result); } long Size2Box(box) BOX2 *box; { return ((long)((box->xh-box->xl)*(box->yh-box->yl))); } long Equal2BoxBox(box1, box2) BOX2 *box1, *box2; { return(FPeq((double)box1->xh, (double)box2->xh) && FPeq((double)box1->xl, (double)box2->xl) && FPeq((double)box1->yh, (double)box2->yh) && FPeq((double)box1->yl, (double)box2->yl)); } long Overlap2BoxBox(box1, box2) BOX2 *box1, *box2; { return(((box1->xh >= box2->xh && box1->xl <= box2->xh) || (box2->xh >= box1->xh && box2->xl <= box1->xh)) && ((box1->yh >= box2->yh && box1->yl <= box2->yh) || (box2->yh >= box1->yh && box2->yl <= box1->yh))); } long Overlap2BoxWin(box1, box2) BOX2 *box1; BOX *box2; { return(((box1->xh >= (float) box2->xh && box1->xl <= (float) box2->xh) || ((float) box2->xh >= box1->xh && (float) box2->xl <= box1->xh)) && ((box1->yh >= (float) box2->yh && box1->yl <= (float) box2->yh) || ((float) box2->yh >= box1->yh && (float) box2->yl <= box1->yh))); } long Overleft2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->xh < box2->xl); } long Left2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->xh <= box2->xh); } long Overright2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->xl > box2->xh); } long Right2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->xl >= box2->xl); } long Overabove2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->yl > box2->yh); } long Above2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->yl > box2->yl); } long Overbelow2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->yh < box2->yl); } long Below2BoxBox(box1, box2) BOX2 *box1, *box2; { return(box1->yh < box2->yh); } long Contain2BoxBox(box1, box2) BOX2 *box1, *box2; { return(Contained2BoxBox(box2, box1)); } long Contained2BoxBox(box1, box2) BOX2 *box1, *box2; { return((box1->xh <= box2->xh && box1->xl >= box2->xl && box1->yh <= box2->yh && box1->yl >= box2->yl)); } LSEG2 * Lseg2_in(str) char *str; { char *coord[LSEG2NARGS], *p; int i; LSEG2 *result; extern double atof(); if (str == NULL) { elog(WARN, "Empty string LSEG2"); return(NULL); } for (i = 0, p = str; *p && i < LSEG2NARGS && *p != RDELIM; p++) if (*p == DELIM || (*p == LDELIM && !i)) coord[i++] = p + 1; if (i < LSEG2NARGS - 1) { elog(WARN, "Ill formatted LSEG2: %s", str); return(NULL); } result = PALLOCTYPE(LSEG2); result->p[0].x = (float) atof(coord[0]); result->p[0].y = (float) atof(coord[1]); result->p[1].x = (float) atof(coord[2]); result->p[1].y = (float) atof(coord[3]); return(result); } char * Lseg2_out(ls) LSEG2 *ls; { char *result; if (ls == NULL) { elog(WARN, "NULL pointer LSEG2"); return(NULL); } result = (char *) palloc(80); (void) sprintf(result, "(%g,%g,%g,%g)", ls->p[0].x, ls->p[0].y, ls->p[1].x, ls->p[1].y); return(result); } /* Point2 functions*/ /*=================*/ /* are the two points identical? */ long Equal2PntPnt(pt1, pt2) POINT2 *pt1, *pt2; { return(FPeq(pt1->x, pt2->x) && FPeq(pt1->x, pt2->x)); } long NotEqual2PntPnt(pt1, pt2) POINT2 *pt1, *pt2; { return !(Equal2PntPnt(pt1, pt2)); } /* DISTANCE */ /* ======== */ double point2_dt(pt1, pt2) POINT2 *pt1, *pt2; { return( HYPOT( (double) (pt1->x - pt2->x), (double) (pt1->y - pt2->y)) ); } double * Distance2PntPnt(p1, p2) POINT2 *p1, *p2; { double *result, point2_dt(); result = PALLOCTYPE(double); *result = point2_dt(p1, p2); return( result ); } /* what is the position of a certain point towards another point? */ /* the point can be situated above of, beneath of, left of or */ /* rigth of the given point */ long Overright2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->x > pnt2->x); return(result); } long Right2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->x > pnt2->x) || FPeq((double) pnt1->x, (double) pnt2->x); return(result); } long Overleft2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->x < pnt2->x); return(result); } long Left2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->x < pnt2->x) || FPeq((double) pnt1->x, (double) pnt2->x); return(result); } long Overabove2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->y > pnt2->y); return(result); } long Overbelow2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { long result; result = (pnt1->y < pnt2->y); return(result); } BOX2 * Union2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { BOX2 *result; result = PALLOCTYPE(BOX2); result->xl = (double) MIN(pnt1->x, pnt2->x); result->yl = (double) MIN(pnt1->y, pnt2->y); result->xh = (double) MAX(pnt1->x, pnt2->x); result->yh = (double) MAX(pnt1->y, pnt2->y); return result; } BOX2 * Inter2PntPnt(pnt1, pnt2) POINT2 *pnt1, *pnt2; { BOX2 *result; result = PALLOCTYPE(BOX2); if (Equal2PntPnt(pnt1, pnt2)) { result->xl = (double) pnt1->x; result->yl = (double) pnt1->y; result->xh = (double) pnt1->x; result->yh = (double) pnt1->y; } else { result = ((BOX2 *) NULL); } return (result); } int Size2Pnt(pnt) POINT2 *pnt; { return(0); } /* POLY STRUCTURES */ /*==================*/ /*---------------------------------------------------------- * External format: * "(poly.npts: xcoord, ycoord,... )" *---------------------------------------------------------*/ /* * POLY2 conversion functions. */ POLYSTRUCT * path2_in(str) char *str; { float coord; extern double atof(); extern long atol(); long field; char *s; int ct, i; POLYSTRUCT *result; long len; if (str == NULL) { elog(WARN, "Empty string POLYSTRUCT"); return(NULL); } /* read the path2 header information */ for (i = 0, s = str; *s && i < 1 && *s != RDELIM ; ++s) if (*s == DELIM || (*s == LDELIM && !i)) { field = atol(s + 1); ++i; } if (i < 0) { elog(WARN, "Ill formatted POLYSTRUCT2: %s", str); return(NULL); } len = POLY2NBYTES(field); result = (POLYSTRUCT*) PALLOC(len); result->length = len; result->npts = field; /* read the path2 points */ ct = result->npts * 2; /* two coords for every point */ for (i = 0; *s && i < ct && *s != RDELIM; ++s) { if (*s == DELIM || *s == MDELIM) { coord = (float) atof(s + 1); if (i % 2) (result->p[i/2]).y = coord; else (result->p[i/2]).x = coord; ++i; } } if (i % 2 || i < --ct) { PFREE(result); elog(WARN, "Ill formatted POLYSTRUCT2: %s", str); return(NULL); } return(result); } char * path2_out(pth) POLYSTRUCT *pth; { char buf[50000], *result, *s; int i; static char tmp[64]; if (pth == NULL) { elog(WARN,"NULL pointer POLYSTRUCT"); return(NULL); } (void) sprintf(buf,"%c%d\0", LDELIM, pth->npts); s = buf + strlen(buf); (void) sprintf(tmp, ":%g,%g", pth->p[0].x, pth->p[0].y); (void) strcpy(s, tmp); s += strlen(tmp); for (i = 1; i < pth->npts; ++i) { (void) sprintf(tmp, ",%g,%g", pth->p[i].x, pth->p[i].y); (void) strcpy(s, tmp); s += strlen(tmp); } *s++ = RDELIM; *s = '\0'; result = (char *) palloc(strlen(buf) + 1); (void) strcpy(result, buf); return(result); } POLYLINE2 * Polyline2_in(str) char *str; { return((POLYLINE2 *) path2_in (str)); } POLYGON2 * Polygon2_in(str) char *str; { return((POLYGON2 *) path2_in(str)); } char * Polyline2_out(pline) POLYLINE2 *pline; { return(path2_out((POLYSTRUCT *) pline)); } char * Polygon2_out(pgon) POLYGON2 *pgon; { return(path2_out((POLYSTRUCT *) pgon)); } POLYSTRUCT * To2Poly(poly) PATH *poly; { int i; POLYSTRUCT *result; long len; if (poly == NULL) { elog(WARN,"NULL pointer PATH"); return(NULL); } len = POLY2NBYTES(poly->npts); result = (POLYSTRUCT*) PALLOC(len); result->length= len; result->npts = poly->npts; for (i=0;inpts;i++){ result->p[i].x = (float) poly->p[i].x; result->p[i].y = (float) poly->p[i].y; } return(result); } PATH * From2Poly(poly) POLYSTRUCT *poly; { int i; PATH *result; long len; if (poly == NULL) { elog(WARN,"NULL pointer POLYSTRUCT"); return(NULL); } len = PATHALLOCSIZE(poly->npts); result = (PATH*) PALLOC(len); result->length= len; result->npts = poly->npts; for (i=0;inpts;i++){ result->p[i].x = (float) poly->p[i].x; result->p[i].y = (float) poly->p[i].y; } return(result); } void Box2Poly_alloc(poly, result) POLYSTRUCT *poly; BOX2 *result; { int i; result->xl = result->yl = (float)HUGE; result->xh = result->yh = (float)-HUGE; for (i=0;inpts;i++){ result->xl = MIN(result->xl, poly->p[i].x); result->yl = MIN(result->yl, poly->p[i].y); result->xh = MAX(result->xh, poly->p[i].x); result->yh = MAX(result->yh, poly->p[i].y); } } BOX2 * Box2Poly(poly) POLYSTRUCT *poly; { BOX2 *result; if (poly == NULL) { elog(WARN,"NULL pointer POLYSTRUCT"); return(NULL); } result = PALLOCTYPE(BOX2); Box2Poly_alloc(poly,result); return(result); } BOX2 * Box2Pln(poly) POLYLINE2 *poly; { return(Box2Poly((POLYSTRUCT *)poly)); } BOX2 * Box2Pgn(poly) POLYGON2 *poly; { return(Box2Poly((POLYSTRUCT *)poly)); } POLYGON2 * To2Pgn(poly) PATH *poly; { return((POLYGON2 *) To2Poly(poly)); } PATH * From2Pgn(poly) POLYGON2 *poly; { PATH *result; result = From2Poly((POLYSTRUCT *)poly); result->closed = 1; return(result); } POLYLINE2 * To2Pln(poly) PATH *poly; { return((POLYLINE2 *) To2Poly(poly)); } PATH * From2Pln(poly) POLYLINE2 *poly; { PATH *result; result = From2Poly((POLYSTRUCT *)poly); result->closed = 0; return(result); } POINT2 * point2_construct(x, y) double x, y; { POINT2 *result; result = PALLOCTYPE(POINT2); result->x = (float) (x); result->y = (float) (y); return(result); } long On2PntBox(pnt, box) POINT2 *pnt; BOX2 *box; { return((pnt->x >= box->xl && pnt->x <= box->xh && pnt->y >= box->yl && pnt->y <= box->yh)); } /* On2PntPgn * Whether a point lies within a polygon. * We use the old O(n) ray method for point-in-polygon. */ long On2PntPgn(pt, pth) POINT2 *pt; POLYGON2 *pth; { #define NEXT(A) ((A+1) % pth->npts) int i, inter, ascend, prev_ascend; double a,b,c; BOX2 box; if (!pt || !pth) return(0); Box2Poly_alloc(pth, &box); if (! On2PntBox(pt,&box)) return(0); inter = 0; prev_ascend = pth->p[pth->npts-1].y < pth->p[0].y; for (i = 0; i < pth->npts; i++) { if(FPeq(pt->y, pth->p[i].y)&&FPeq(pt->x, pth->p[i].x)) return(1); if(FPeq(pth->p[NEXT(i)].y, pth->p[i].y) && FPeq(pt->y, pth->p[i].y) && (pt->x >= MIN(pth->p[NEXT(i)].x, pth->p[i].x)) && (pt->x <= MAX(pth->p[NEXT(i)].x, pth->p[i].x))) return(1); if(!(((pt->y < pth->p[i].y) && (pt->y <= pth->p[NEXT(i)].y)) || ((pt->y > pth->p[i].y) && (pt->y >= pth->p[NEXT(i)].y)) || FPeq(pth->p[i].y, pth->p[NEXT(i)].y))) { ascend = pth->p[i].y < pth->p[NEXT(i)].y; if(!(FPeq(pt->y, pth->p[i].y) && (prev_ascend != ascend))) { a = ((double)(pth->p[i].x-pth->p[NEXT(i)].x)) / ((double)(pth->p[i].y-pth->p[NEXT(i)].y)); b = pth->p[NEXT(i)].x-a*pth->p[NEXT(i)].y; c = a*pt->y+b; if(FPeq(pt->x, c)) return(1); if(pt->x < c) ++inter; } } if(!FPeq(pth->p[i].y, pth->p[NEXT(i)].y)) prev_ascend = pth->p[i].y < pth->p[NEXT(i)].y; } return(inter % 2); /* odd # of intersections */ } long On2PntPln(pt, pth) POINT2 *pt; POLYLINE2 *pth; { int i, n; double point2_dt(), a, b; n = pth->npts - 1; a = point2_dt(pt, &pth->p[0]); for (i = 0; i < n; i++) { b = point2_dt(pt, &pth->p[i+1]); if (FPeq(a+b, point2_dt(&pth->p[i], &pth->p[i+1]))) return(1); a = b; } return(0); } double area(pgon) POLYGON2 *pgon; { double result; int i; result = 0; for (i = 0; i < pgon->npts; ++i) result += ( (pgon->p[i].x) - pgon->p[(i+1) % (pgon->npts)].x) * ( (pgon->p[i].y) + pgon->p[(i+1) % (pgon->npts)].y); result /= 2.0; return(result); } double chain_length(pth) POLYSTRUCT *pth; { double result; int i; result = 0; for (i = 0; i < (pth->npts-1); ++i) result += point2_dt(&pth->p[i], &pth->p[i+1]); return(result); } double mindist_ppoly (pt, poly) POINT2 *pt; POLYSTRUCT *poly; { int i; double min, dist, mindist_pl(); LSEG2 *ls; min = HUGE; ls = PALLOCTYPE(LSEG2); for (i = 0; i < poly->npts; i++) { ls->p[0].x = poly->p[i].x; ls->p[0].y = poly->p[i].y; ls->p[1].x = poly->p[(i+1) % (poly->npts)].x; ls->p[1].y = poly->p[(i+1) % (poly->npts)].y; dist = mindist_pl(pt, ls); min = MIN(min, dist); } PFREE(ls); return(min); } double mindist_pl(pt, ls) POINT2 *pt; LSEG2 *ls; { double result, dx, dy, c0, c1, c; if (Equal2PntPnt(&ls->p[0],&ls->p[1])) { result = point2_dt(&ls->p[0], pt); } else { dx = (double)(ls->p[1].x - ls->p[0].x); dy = (double)(ls->p[1].y - ls->p[0].y); c0 = dx * ls->p[0].x + dy * ls->p[0].y; c1 = dx * ls->p[1].x + dy * ls->p[1].y; c = dx * pt->x + dy * pt->y; if (c > c1) { result = point2_dt(&ls->p[1], pt); } else if (c < c0) { result = point2_dt(&ls->p[0], pt); } else { c0 = dy * ls->p[0].x - dx * ls->p[0].y; c = dy * pt->x - dx * pt->y; result = fabs(c-c0)/sqrt(dx*dx+dy*dy); } } return(result); } double maxdist_ppoly (pt, poly) POINT2 *pt; POLYSTRUCT *poly; { int i; double dist, max, point2_dt(); max = -HUGE; for (i = 0; i < poly->npts; i++) { dist = point2_dt(&poly->p[i],pt); max = MAX(max, dist); } return(max); } double * MaxDist2PntPln(pt, pl) POINT2 *pt; POLYLINE2 *pl; { double *result, maxdist_ppoly(); result = PALLOCTYPE(double); *result = maxdist_ppoly(pt, ((POLYSTRUCT *) pl)); return(result); } double * MinDist2PntPln(pt, pline) POINT2 *pt; POLYLINE2 *pline; { double *result, mindist_ppoly(); result = PALLOCTYPE(double); *result = mindist_ppoly(pt, ((POLYSTRUCT *) pline)); return(result); } double * Length2Pln(pln) POLYLINE2 *pln; { double *result, chain_length(); result = PALLOCTYPE(double); *result = chain_length((POLYSTRUCT *) pln); return(result); } long Equal2PlnPln(pl1, pl2) POLYLINE2 *pl1, *pl2; { long n, result; long Equal2PntPnt(); int i; if (pl1->npts == pl2->npts) { result = 1; n = pl1->npts; if (Equal2PntPnt(&(pl1->p[0]), &(pl2->p[0]))) { for (i = 0; i < n; i++) result = (result && Equal2PntPnt(&(pl1->p[i]), &(pl2->p[i]))); } else if (Equal2PntPnt(&(pl1->p[0]), &(pl2->p[n-1]))) { for (i = 0; i < n; i++) result = (result && Equal2PntPnt(&(pl1->p[i]), &(pl2->p[n-1-i]))); } else result = 0; } else result = 0; return (result); } POINT2 * PointSelector(pline, nr) POLYLINE2 *pline; long nr; { POINT2 *result; result = PALLOCTYPE(POINT2); if (0 < nr && nr <= pline->npts) nr--; else if (nr < 0 && nr >= -pline->npts) nr = pline->npts + nr; else { elog(WARN, "Not a legal PointSelector: %ld", nr); return(NULL); } result->x = pline->p[nr].x; result->y = pline->p[nr].y; return(result); } double * MaxDist2PntPgn(pt, pgon) POINT2 *pt; POLYGON2 *pgon; { double *result, maxdist_ppoly(); result = PALLOCTYPE(double); *result = maxdist_ppoly(pt, ((POLYSTRUCT *) pgon)); return(result); } double * MinDist2PntPgn(pt, pgon) POINT2 *pt; POLYGON2 *pgon; { double *result, mindist_ppoly(), dist, mindist_pl(); LSEG2 *ls; result = PALLOCTYPE(double); if (On2PntPgn(pt, pgon)) { *result=0.0; } else { ls = PALLOCTYPE(LSEG2); ls->p[0].x = pgon->p[pgon->npts-1].x; ls->p[0].y = pgon->p[pgon->npts-1].y; ls->p[1].x = pgon->p[0].x; ls->p[1].y = pgon->p[0].y; dist = mindist_pl(pt, ls); *result = mindist_ppoly(pt, ((POLYSTRUCT *) pgon)); *result = MIN(*result, dist); PFREE(ls); } return(result); } double * Area2Pgn(pgon) POLYGON2 *pgon; { double area(), *result; result = PALLOCTYPE(double); *result = area(pgon); *result = fabs(*result); return(result); } long Orientation2Pgn(pgon) { double area(); return((area(pgon)>0)); } POINT2 * GravCenter2Pgn(pgon) POLYGON2 *pgon; { POINT2 *result; double area(), resx, resy, pgarea, ar; int n, i; result = PALLOCTYPE(POINT2); pgarea = area(pgon); pgarea *= 6; resx = 0; resy = 0; n = pgon->npts; for (i=0; i< n; ++i) { ar = ((pgon->p[i].x * pgon->p[(i+1) % n].y) - (pgon->p[i].y * pgon->p[(i+1) % n].x)); resx += (pgon->p[i].x + pgon->p[(i+1) % n].x) * ar; resy += (pgon->p[i].y + pgon->p[(i+1) % n].y) * ar; } result->x = resx / pgarea; result->y = resy / pgarea; return(result); } double * Perimeter2Pgn(pgon) POLYGON2 *pgon; { double point2_dt(), *result, chain_length(); result = PALLOCTYPE(double); *result = chain_length((POLYSTRUCT *) pgon); *result += point2_dt(&pgon->p[pgon->npts-1], &pgon->p[0]); return(result); } long Equal2PgnPgn(pl1, pl2) POLYGON2 *pl1, *pl2; { long n, result; long Equal2PntPnt(); int i, j; if (pl1->npts == pl2->npts) { result = 1; n = pl1->npts; i = 0; while ((!Equal2PntPnt(&(pl1->p[0]), &(pl2->p[i]))) && (i < n)) i++; if (i != n) { if (Equal2PntPnt(&(pl1->p[1]), &(pl2->p[(i+1) % n]))) for (j = 0; j < n; j++) result=(result && Equal2PntPnt(&pl1->p[j], &pl2->p[(i+j) % n])); else if (Equal2PntPnt(&pl1->p[1], &pl2->p[(i-1) % n])) for (j = 0; j < n; j++) result=(result && Equal2PntPnt(&pl1->p[j], &pl2->p[(i-j) % n])); } else result = 0; } else result = 0; return (result); } long lseg2_overlap(ln1, ln2) LSEG2 *ln1, *ln2; { double dx1, dx2, dy1, dy2; long result; dx1 = (double) ln1->p[1].x - ln1->p[0].x; dx2 = (double) ln2->p[1].x - ln2->p[0].x; dy1 = (double) ln1->p[1].y - ln1->p[0].y; dy2 = (double) ln2->p[1].y - ln2->p[0].y; if (FPzero(dx1) || FPzero(dx2)) { result = FPeq(dx1,dx2); } else { result = FPeq(dy1/dx1,dy2/dx2); } return(result && (FPzero(mindist_pl(&ln1->p[0], ln2)) || FPzero(mindist_pl(&ln1->p[1], ln2)) || FPzero(mindist_pl(&ln2->p[0], ln1)) || FPzero(mindist_pl(&ln2->p[1], ln1)))); } long Neighbour2PgnPgn(pb1, pb2) POLYGON2 *pb1, *pb2; { LSEG2 border1, border2; BOX2 box1, box2; long result, n1, n2; int i, j; if (!pb1 || !pb2) return(NULL); if (Equal2PgnPgn(pb1, pb2)) return(NULL); Box2Poly_alloc(pb1, &box1); Box2Poly_alloc(pb2, &box2); if (! Overlap2BoxBox(&box1,&box2)) return(NULL); n1 = pb1->npts; n2 = pb2->npts; result = 0; for (i = 0; !result && i < n1; i++) { border1.p[0] = pb1->p[i]; border1.p[1] = pb1->p[(i+1) % n1]; for (j = 0; !result && j < n2; j++) { border2.p[0] = pb2->p[j]; border2.p[1] = pb2->p[(j+1) % n2]; if (result= (lseg2_overlap(&border1, &border2))) break; } } return(result); } long NotNeighb2PgnPgn(pb1, pb2) POLYGON2 *pb1, *pb2; { return ! Neighbour2PgnPgn(pb1, pb2); }