/********************************************************************** * $Id: mitab_geometry.cpp,v 1.5 2004/06/30 20:29:04 dmorissette Exp $ * * Name:     mitab_geometry.cpp * Project:  MapInfo TAB Read/Write library * Language: C++ * Purpose:  Geometry manipulation functions. * Author:   Daniel Morissette, dmorissette@dmsolutions.ca *           Based on functions from mapprimitive.c/mapsearch.c in the source *           of UMN MapServer by Stephen Lime (http://mapserver.gis.umn.edu/) ********************************************************************** * Copyright (c) 1999-2001, Daniel Morissette * * Permission is hereby granted, free of charge, to any person obtaining a * copy of this software and associated documentation files (the "Software"), * to deal in the Software without restriction, including without limitation * the rights to use, copy, modify, merge, publish, distribute, sublicense, * and/or sell copies of the Software, and to permit persons to whom the * Software is furnished to do so, subject to the following conditions: *  * The above copyright notice and this permission notice shall be included * in all copies or substantial portions of the Software. *  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER  * DEALINGS IN THE SOFTWARE. ********************************************************************** * * $Log: mitab_geometry.cpp,v $ * Revision 1.5  2004/06/30 20:29:04  dmorissette * Fixed refs to old address danmo@videotron.ca * * Revision 1.4  2001/12/18 23:42:28  daniel * Added a test in OGRPolygonLabelPoint() to prevent returning a point * outside of the polygon MBR (bug 673). * * Revision 1.3  2001/01/22 16:03:58  warmerda * expanded tabs * * Revision 1.2  2000/09/28 16:39:44  warmerda * avoid warnings for unused, and unitialized variables * * Revision 1.1  2000/09/19 17:19:40  daniel * Initial Revision * **********************************************************************/#include "ogr_geometry.h"#define OGR_NUM_RINGS(poly)   (poly->getNumInteriorRings()+1)#define OGR_GET_RING(poly, i) (i==0?poly->getExteriorRing():poly->getInteriorRing(i-1))/********************************************************************** *                   OGRPointInRing() * * Returns TRUE is point is inside ring, FALSE otherwise * * Adapted version of msPointInPolygon() from MapServer's mapsearch.c **********************************************************************/GBool OGRPointInRing(OGRPoint *poPoint, OGRLineString *poRing){    int i, j, numpoints;    GBool status = FALSE;    double x, y;    numpoints = poRing->getNumPoints();    x = poPoint->getX();    y = poPoint->getY();    for (i = 0, j = numpoints-1; i < numpoints; j = i++)     {        if ((((poRing->getY(i)<=y) && (y<poRing->getY(j))) ||              ((poRing->getY(j)<=y) && (y<poRing->getY(i)))) &&             (x < (poRing->getX(j) - poRing->getX(i)) * (y - poRing->getY(i)) /                  (poRing->getY(j) - poRing->getY(i)) + poRing->getX(i)))            status = !status;    }    return status;}/********************************************************************** *                   OGRIntersectPointPolygon() * * Instead of using ring orientation we count the number of parts the * point falls in. If odd the point is in the polygon, if 0 or even * then the point is in a hole or completely outside. * * Returns TRUE is point is inside polygon, FALSE otherwise * * Adapted version of msIntersectPointPolygon() from MapServer's mapsearch.c **********************************************************************/GBool OGRIntersectPointPolygon(OGRPoint *poPoint, OGRPolygon *poPoly) {    int i;    GBool status = FALSE;    for(i=0; i<OGR_NUM_RINGS(poPoly); i++)     {        if (OGRPointInRing(poPoint, OGR_GET_RING(poPoly, i)))         {            /* ok, the point is in a line */            status = !status;        }    }    return status;}/********************************************************************** *                   OGRPolygonLabelPoint() * * Generate a label point on the surface of a polygon. *  * The function is based on a scanline conversion routine used for polygon  * fills.  Instead of processing each line the as with drawing, the * polygon is sampled. The center of the longest sample is chosen for the  * label point. The label point is guaranteed to be in the polygon even if  * it has holes assuming the polygon is properly formed. *  * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise. * * Adapted version of msPolygonLabelPoint() from MapServer's mapprimitive.c **********************************************************************/typedef enum {CLIP_LEFT, CLIP_MIDDLE, CLIP_RIGHT} CLIP_STATE;#define CLIP_CHECK(min, a, max) ((a) < (min) ? CLIP_LEFT : ((a) > (max) ? CLIP_RIGHT : CLIP_MIDDLE));#define ROUND(a)       ( (a) + 0.5 )#define SWAP( a, b, t) ( (t) = (a), (a) = (b), (b) = (t) )#define EDGE_CHECK( x0, x, x1) ((x) < MIN( (x0), (x1)) ? CLIP_LEFT : ((x) > MAX( (x0), (x1)) ? CLIP_RIGHT : CLIP_MIDDLE ))#define NUM_SCANLINES 5int OGRPolygonLabelPoint(OGRPolygon *poPoly, OGRPoint *poLabelPoint){    double        slope;    OGRRawPoint   point1, point2;    int           i, j, k, nfound;    double        x, y, *xintersect, temp;    double        hi_y, lo_y;    int           wrong_order, n;    double        len, max_len=0;    double        skip;    OGREnvelope   oEnv;    if (poPoly == NULL)        return OGRERR_FAILURE;    poPoly->getEnvelope(&oEnv);    poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);    poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);    //if(get_centroid(p, lp, &miny, &maxy) == -1) return(-1);    if(OGRIntersectPointPolygon(poLabelPoint, poPoly) == TRUE) /* cool, done */        return OGRERR_NONE;    /* do it the hard way - scanline */    skip = (oEnv.MaxY - oEnv.MinY)/NUM_SCANLINES;    n=0;    for(j=0; j<OGR_NUM_RINGS(poPoly); j++)     {        /* count total number of points */        n += OGR_GET_RING(poPoly, j)->getNumPoints();    }    xintersect = (double *)calloc(n, sizeof(double));    if (xintersect == NULL)        return OGRERR_FAILURE;    for(k=1; k<=NUM_SCANLINES; k++)     {         /* sample the shape in the y direction */            y = oEnv.MaxY - k*skip;         /* need to find a y that won't intersect any vertices exactly */          hi_y = y - 1; /* first initializing lo_y, hi_y to be any 2 pnts on either side of lp->y */        lo_y = y + 1;        for(j=0; j<OGR_NUM_RINGS(poPoly); j++)         {            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);            if((lo_y < y) && (hi_y >= y))                 break; /* already initialized */            for(i=0; i < poRing->getNumPoints(); i++)             {                   if((lo_y < y) && (hi_y >= y))                     break; /* already initialized */                if(poRing->getY(i) < y)                    lo_y = poRing->getY(i);                if(poRing->getY(i) >= y)                    hi_y = poRing->getY(i);            }        }        n=0;        for(j=0; j<OGR_NUM_RINGS(poPoly); j++)         {            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);            for(i=0; i < poRing->getNumPoints(); i++)             {                if((poRing->getY(i) < y) &&                    ((y - poRing->getY(i)) < (y - lo_y)))                    lo_y = poRing->getY(i);                if((poRing->getY(i) >= y) &&                    ((poRing->getY(i) - y) < (hi_y - y)))                    hi_y = poRing->getY(i);            }              }        if(lo_y == hi_y)             return OGRERR_FAILURE;        else              y = (hi_y + lo_y)/2.0;                nfound = 0;        for(j=0; j<OGR_NUM_RINGS(poPoly); j++)   /* for each line */        {            OGRLinearRing *poRing = OGR_GET_RING(poPoly,j);            point1.x = poRing->getX(poRing->getNumPoints()-1);            point1.y = poRing->getY(poRing->getNumPoints()-1);            for(i=0; i < poRing->getNumPoints(); i++)             {                point2.x = poRing->getX(i);                point2.y = poRing->getY(i);                        if(EDGE_CHECK(point1.y, y, point2.y) == CLIP_MIDDLE)                 {                    if(point1.y == point2.y)                        continue; /* ignore horizontal edges */                    else                        slope = (point2.x - point1.x) / (point2.y - point1.y);                              x = point1.x + (y - point1.y)*slope;                    xintersect[nfound++] = x;                } /* End of checking this edge */                        point1 = point2;  /* Go on to next edge */            }        } /* Finished the scanline */            /* First, sort the intersections */        do         {            wrong_order = 0;            for(i=0; i < nfound-1; i++)             {                if(xintersect[i] > xintersect[i+1])                 {                    wrong_order = 1;                    SWAP(xintersect[i], xintersect[i+1], temp);                }            }        } while(wrong_order);            /* Great, now find longest span */        point1.y = point2.y = y;        for(i=0; i < nfound; i += 2)         {            point1.x = xintersect[i];            point2.x = xintersect[i+1];            /* len = length(point1, point2); */            len = ABS((point2.x - point1.x));            if(len > max_len)             {                max_len = len;                poLabelPoint->setX( (point1.x + point2.x)/2 );                poLabelPoint->setY( y );            }        }    }    free(xintersect);    /* __TODO__ Bug 673     * There seem to be some polygons for which the label is returned     * completely outside of the polygon's MBR and this messes the      * file bounds, etc.     * Until we find the source of the problem, we'll at least validate     * the label point to make sure that it overlaps the polygon MBR.     */    if( poLabelPoint->getX() < oEnv.MinX        || poLabelPoint->getY() < oEnv.MinY        || poLabelPoint->getX() > oEnv.MaxX        || poLabelPoint->getY() > oEnv.MaxY )    {        // Reset label coordinates to center of MBR, just in case        poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);        poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);        // And return an error        return OGRERR_FAILURE;    }    if(max_len > 0)        return OGRERR_NONE;    else        return OGRERR_FAILURE;}/********************************************************************** *                   OGRGetCentroid() * * Calculate polygon gravity center. *  * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise. * * Adapted version of get_centroid() from MapServer's mapprimitive.c **********************************************************************/int OGRGetCentroid(OGRPolygon *poPoly, OGRPoint *poCentroid){    int i,j;    double cent_weight_x=0.0, cent_weight_y=0.0;    double len, total_len=0;    for(i=0; i<OGR_NUM_RINGS(poPoly); i++)     {        double x1, y1, x2, y2;        OGRLinearRing *poRing = OGR_GET_RING(poPoly, i);        x2 = poRing->getX(0);        y2 = poRing->getY(0);        for(j=1; j<poRing->getNumPoints(); j++)         {            x1 = x2;            y1 = y2;            x2 = poRing->getX(j);            y2 = poRing->getY(j);            len = sqrt( pow((x2-x1),2) + pow((y2-y1),2) );            cent_weight_x += len * ((x1 + x2)/2.0);            cent_weight_y += len * ((y1 + y2)/2.0);            total_len += len;        }    }    if(total_len == 0)        return(OGRERR_FAILURE);    poCentroid->setX( cent_weight_x / total_len );    poCentroid->setY( cent_weight_y / total_len );      return OGRERR_NONE;}/********************************************************************** *                   OGRPolylineCenterPoint() * * Return the center point of a polyline. * * In MapInfo, for a simple or multiple polyline (pline), the center point  * in the object definition is supposed to be either the center point of  * the pline or the first section of a multiple pline (if an odd number of  * points in the pline or first section), or the midway point between the  * two central points (if an even number of points involved).  * * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise. **********************************************************************/int OGRPolylineCenterPoint(OGRLineString *poLine, OGRPoint *poLabelPoint){    if (poLine == NULL || poLine->getNumPoints() < 2)        return OGRERR_FAILURE;    if (poLine->getNumPoints() % 2 == 0)    {        // Return the midway between the 2 center points        int i = poLine->getNumPoints()/2;        poLabelPoint->setX( (poLine->getX(i-1) + poLine->getX(i))/2.0 );        poLabelPoint->setY( (poLine->getY(i-1) + poLine->getY(i))/2.0 );    }    else    {        // Return the center point        poLine->getPoint(poLine->getNumPoints()/2, poLabelPoint);    }        return OGRERR_NONE;}/********************************************************************** *                   OGRPolylineLabelPoint() * * Generate a label point on a polyline: The center of the longest segment. * * Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise. **********************************************************************/int OGRPolylineLabelPoint(OGRLineString *poLine, OGRPoint *poLabelPoint){    double      segment_length, max_segment_length = 0.0;    double      x1, y1, x2, y2;    if (poLine == NULL || poLine->getNumPoints() < 2)        return OGRERR_FAILURE;    max_segment_length = -1.0;    x2 = poLine->getX(0);    y2 = poLine->getY(0);    for(int i=1; i<poLine->getNumPoints(); i++)     {        x1 = x2;        y1 = y2;        x2 = poLine->getX(i);        y2 = poLine->getY(i);        segment_length = pow((x2-x1),2) + pow((y2-y1),2);        if (segment_length > max_segment_length)        {            max_segment_length = segment_length;            poLabelPoint->setX( (x1 + x2)/2.0 );            poLabelPoint->setY( (y1 + y2)/2.0 );        }    }        return OGRERR_NONE;}