gistfile1.txt

// MIT license
// Copyright 2022, Even Rouault

#include <cmath>
#include <algorithm>
#include <utility>
#include <vector>

/************************************************************************/
/*                         getOrientation()                             */
/************************************************************************/

typedef std::pair<double, double> XYPair;

// Returns 1 whether (p1,p2,p3) is clockwise oriented,
// -1 if it is counter-clockwise oriented,
// or 0 if it is colinear.
static int getOrientation(const XYPair& p1, const XYPair& p2, const XYPair& p3)
{
    const double p1x = p1.first;
    const double p1y = p1.second;
    const double p2x = p2.first;
    const double p2y = p2.second;
    const double p3x = p3.first;
    const double p3y = p3.second;
    const double val = (p2y - p1y) * (p3x - p2x) - (p2x - p1x) * (p3y - p2y);
    if( std::abs(val) < 1e-20 )
        return 0;
    else if( val > 0 )
        return 1;
    else
        return -1;
}

/************************************************************************/
/*                          isConvex()                                  */
/************************************************************************/

typedef std::vector<XYPair> XYPoly;

// poly must be closed
static bool isConvex(const XYPoly& poly)
{
    const size_t n = poly.size();
    size_t i = 0;
    int last_orientation = getOrientation(poly[i], poly[i + 1], poly[i + 2]);
    ++ i;
    for(; i < n - 2; ++i )
    {
        const int orientation = getOrientation(poly[i], poly[i + 1], poly[i + 2]);
        if( orientation != 0 )
        {
            if( last_orientation == 0 )
                last_orientation = orientation;
            else if( orientation != last_orientation )
                return false;
        }
    }
    return true;
}

/************************************************************************/
/*                     pointIntersectsConvexPoly()                      */
/************************************************************************/

// Returns whether xy intersects poly, that must be closed and convex.
static bool pointIntersectsConvexPoly(const XYPair& xy, const XYPoly& poly)
{
    const size_t n = poly.size();
    double dx1 = xy.first - poly[0].first;
    double dy1 = xy.second - poly[0].second;
    double dx2 = poly[1].first - poly[0].first;
    double dy2 = poly[1].second - poly[0].second;
    double prevCrossProduct = dx1 * dy2 - dx2 * dy1;

    // Check if the point remains on the same side (left/right) of all edges
    for( size_t i = 2; i < n; i++ )
    {
        dx1 = xy.first - poly[i-1].first;
        dy1 = xy.second - poly[i-1].second;

        dx2 = poly[i].first - poly[i-1].first;
        dy2 = poly[i].second - poly[i-1].second;

        double crossProduct = dx1 * dy2 - dx2 * dy1;
        if( std::abs(prevCrossProduct) < 1e-20 )
            prevCrossProduct = crossProduct;
        else if( prevCrossProduct * crossProduct < 0 )
            return false;
    }

    return true;
}

/************************************************************************/
/*                     getIntersection()                                */
/************************************************************************/

/* Returns intersection of [p1,p2] with [p3,p4], if
 * it is a single point, and the 2 segments are not colinear.
 */
static bool getIntersection(const XYPair& p1, const XYPair& p2,
                            const XYPair& p3, const XYPair& p4,
                            XYPair& xy)
{
    const double x1 = p1.first;
    const double y1 = p1.second;
    const double x2 = p2.first;
    const double y2 = p2.second;
    const double x3 = p3.first;
    const double y3 = p3.second;
    const double x4 = p4.first;
    const double y4 = p4.second;
    const double t_num = (x1 - x3) * (y3 - y4) - (y1 - y3) * (x3 - x4);
    const double denom = (x1 - x2) * (y3 - y4) - (y1 - y2) * (x3 - x4);
    if( t_num * denom < 0 || std::abs(t_num) > std::abs(denom) || denom == 0 )
        return false;

    const double u_num = (x1 - x3) * (y1 - y2) - (y1 - y3) * (x1 - x2);
    if( u_num * denom < 0 || std::abs(u_num) > std::abs(denom) )
        return false;

    const double t = t_num / denom;
    xy.first = x1 + t * (x2 - x1);
    xy.second = y1 + t * (y2 - y1);
    return true;
}

/************************************************************************/
/*                     getConvexPolyIntersection()                      */
/************************************************************************/

// poly1 and poly2 must be closed and convex.
// The returned intersection will not necessary be closed.
static void getConvexPolyIntersection(const XYPoly& poly1,
                                      const XYPoly& poly2,
                                      XYPoly& intersection)
{
    intersection.clear();

    // Add all points of poly1 inside poly2
    for( size_t i = 0; i < poly1.size() - 1; ++i)
    {
        if( pointIntersectsConvexPoly(poly1[i], poly2) )
            intersection.push_back(poly1[i]);
    }
    if( intersection.size() == poly1.size() - 1 )
    {
        // poly1 is inside poly2
        return;
    }

    // Add all points of poly2 inside poly1
    for( size_t i = 0; i < poly2.size() - 1; ++i)
    {
        if( pointIntersectsConvexPoly(poly2[i], poly1) )
            intersection.push_back(poly2[i]);
    }

    // Compute the intersection of all edges of both polygons
    XYPair xy;
    for( size_t i1 = 0; i1 < poly1.size() - 1; ++i1 )
    {
        for( size_t i2 = 0; i2 < poly2.size() - 1; ++i2)
        {
            if( getIntersection(poly1[i1],
                                poly1[i1+1],
                                poly2[i2],
                                poly2[i2+1],
                                xy) )
            {
                intersection.push_back(xy);
            }
        }
    }

    if( intersection.empty() )
        return;

    // Find lowest-left point in intersection set
    double lowest_x = std::numeric_limits<double>::max();
    double lowest_y = std::numeric_limits<double>::max();
    for( const auto& pair: intersection )
    {
        const double x = pair.first;
        const double y = pair.second;
        if( y < lowest_y || (y == lowest_y && x < lowest_x) )
        {
            lowest_x = x;
            lowest_y = y;
        }
    }

    const auto sortFunc = [&](const XYPair& p1,
                              const XYPair& p2)
    {
        const double p1x_diff = p1.first - lowest_x;
        const double p1y_diff = p1.second - lowest_y;
        const double p2x_diff = p2.first - lowest_x;
        const double p2y_diff = p2.second - lowest_y;
        if( p2y_diff == 0.0 && p1y_diff == 0.0 )
        {
            if( p1x_diff >= 0 )
            {
                if( p2x_diff >= 0 )
                    return p1.first < p2.first;
                return true;
            }
            else
            {
                if( p2x_diff >= 0)
                    return false;
                return p1.first < p2.first;
            }
        }

        if( p2x_diff == 0.0 && p1x_diff == 0.0 )
            return p1.second < p2.second;

        double tan_p1;
        if( p1x_diff == 0.0 )
            tan_p1 = p1y_diff == 0.0 ? 0.0 : std::numeric_limits<double>::max();
        else
            tan_p1 = p1y_diff / p1x_diff;

        double tan_p2;
        if( p2x_diff == 0.0 )
            tan_p2 = p2y_diff == 0.0 ? 0.0 : std::numeric_limits<double>::max();
        else
            tan_p2 = p2y_diff / p2x_diff;

        if( tan_p1 >= 0 )
        {
            if( tan_p2 >= 0 )
                return tan_p1 < tan_p2;
            else
                return true;
        }
        else
        {
            if( tan_p2 >= 0 )
                return false;
            else
                return tan_p1 < tan_p2;
        }
    };

    // Sort points by increasing atan2(y-lowest_y, x-lowest_x) to form a convex hull
    std::sort(intersection.begin(), intersection.end(), sortFunc);

    // Remove duplicated points
    size_t j = 1;
    for( size_t i = 1; i < intersection.size(); ++i )
    {
        if( intersection[i] != intersection[i-1] )
        {
            if( j < i )
                intersection[j] = intersection[i];
            ++j;
        }
    }
    intersection.resize(j);
}

/************************************************************************/
/*                            getArea()                                 */
/************************************************************************/

// poly may or may not be closed.
static double getArea(const XYPoly& poly)
{
    // CPLAssert(poly.size() >= 2);
    const size_t nPointCount = poly.size();
    double dfAreaSum =
        poly[0].first * (poly[1].second - poly[nPointCount-1].second);

    for( size_t i = 1; i < nPointCount-1; i++ )
    {
        dfAreaSum += poly[i].first * (poly[i+1].second - poly[i-1].second);
    }

    dfAreaSum +=
        poly[nPointCount-1].first *
        (poly[0].second - poly[nPointCount-2].second);

    return 0.5 * std::fabs(dfAreaSum);
}