/**********************************************************************
 * $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 5

int 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;
}