# pointinpoly.py -- V.1.2, Feb 12 2007, by leonardo maffi


def triangleContainsQ( ((x1,y1), (x2,y2), (x3,y3)), (px,py) ):
    """triangleContainsQ( ((x1,y1), (x2,y2), (x3,y3)), (px,py) ): True if point
    (px,py) is contained in the triangle of given vertices.
    Sometimes the edges aren't part of the triangle.

    >>> triangle1 = [(0,0), (3,0), (3,4)]
    >>> triangle1_test = [(-1,1), (2,-1), (0,0), (3,0), (2,1), (3,1), (2,2), (1,3), (3.5,3), (4,3)]
    >>> [triangleContainsQ(triangle1, p) for p in triangle1_test]
    [False, False, False, False, True, False, True, False, False, False]
    >>> triangle2 = [(0,0), (3,0), (5,0)]
    >>> triangle2_test = [(0,0), (2,0), (3,0), (2,2)]
    >>> [triangleContainsQ(triangle2, p) for p in triangle2_test]
    [False, False, False, False]
    """
    # Adapted from Patrick Nichols (pnichols@mit.edu) www.muxalicious.com/programming.php
    fAB = (py-y1)*(x2-x1) - (px-x1)*(y2-y1)
    fCA = (py-y3)*(x1-x3) - (px-x3)*(y1-y3)
    fBC = (py-y2)*(x3-x2) - (px-x2)*(y3-y2)
    return fAB * fBC > 0 and fBC * fCA > 0


def polyContainsQ(vertices, (x, y)):
    """polyContainsQ(vertices, (x, y)): return True if the given (x,y) point
    is inside the given polygon. Polygon is a list of (x,y) pairs.
    Points that are on the boundary are classified as either inside or outside.
    Any particular point is always classified consistently the same way.
    It considers each polygon to be topologically a semi-open set.

    >>> triangle1 = [(0,0), (3,0), (3,4)]
    >>> triangle2 = [(0,0), (3,0), (5,0)]
    >>> rectangle = [(0,0), (4,0), (4,3), (0,3)]
    >>> octangle = [(1,0), (4,0), (5,2), (4,3), (3,3), (2,3), (1,3), (0,2)]

    >>> triangle1_test = [(-1,1), (2,-1), (0,0), (3,0), (2,1), (3,1), (2,2), (1,3), (3.5,3), (4,3)]
    >>> triangle2_test = [(0,0), (2,0), (3,0), (2,2)]
    >>> rectangle_test = [(1,2), (4,2), (5,2), (4,3)]
    >>> octangle_test = [(3,2), (4,2), (5,2), (5.5,2), (6,2)]

    # Note that on this triangle on (0,0) polyContainsQ and triangleContainsQ give different results
    >>> [polyContainsQ(triangle1, p) for p in triangle1_test]
    [False, False, True, False, True, False, True, False, False, False]
    >>> [polyContainsQ(triangle2, p) for p in triangle2_test]
    [False, False, False, False]
    >>> [polyContainsQ(rectangle, p) for p in rectangle_test]
    [True, False, False, False]
    >>> [polyContainsQ(octangle, p) for p in octangle_test]
    [True, True, False, False, False]
    """
    # From Comp.Graphics.Algorithms, Frequently Asked Questions,
    # Section 2. 2D Polygon Computations
    # The essence of the ray-crossing method is as follows.
    # Think of standing inside a field with a fence representing the polygon.
    # Then walk north. If you have to jump the fence you know you are now
    # outside the poly. If you have to cross again you know you are now
    # inside again; i.e., if you were inside the field to start with, the total
    # number of fence jumps you would make will be odd, whereas if you were
    # ouside the jumps will be even.
    #
    # The code below is from Wm. Randolph Franklin <wrf@ecse.rpi.edu>
    # (see URL below) with some minor modifications for speed.  It returns
    # 1 for strictly interior points, 0 for strictly exterior, and 0 or 1
    # for points on the boundary.  The boundary behavior is complex but
    # determined; in particular, for a partition of a region into polygons,
    # each point is "in" exactly one polygon.
    # (See p.243 of [O'Rourke (C)] for a discussion of boundary behavior.)
    #
    # int pnpoly(int npol, float *xp, float *yp, float x, float y) {
    #   int i, j, c = 0;
    #   for (i = 0, j = npol-1; i < npol; j = i++) {
    #     if ((((yp[i]<=y) && (y<yp[j])) ||
    #          ((yp[j]<=y) && (y<yp[i]))) &&
    #         (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i]))
    #       c = !c;
    #   }
    #   return c;
    # }
    #
    # The code may be further accelerated, at some loss in clarity, by
    # avoiding the central computation when the inequality can be deduced,
    # and by replacing the division by a multiplication for those processors
    # with slow divides.  For code that distinguishes strictly interior
    # points from those on the boundary, see [O'Rourke (C)] pp. 239-245.
    # For a method based on winding number, see Dan Sunday,
    # "Fast Winding Number Test for Point Inclusion in a Polygon,"
    # http://softsurfer.com/algorithms.htm, March 2001.
    #
    # References:
    # Franklin's code:
    #    http://www.ecse.rpi.edu/Homepages/wrf/research/geom/pnpoly.html
    # [Gems IV]  pp. 24-46
    inside = False
    vxj = vertices[-1][0]
    vyj = vertices[-1][1]
    for vxi,vyi in vertices:
        if (((vyi <= y) and (y < vyj)) or ((vyj <= y) and (y < vyi))) and \
            (x < (vxj - vxi) * (y - vyi) / (vyj - vyi) + vxi):
            inside = not inside
        vxj = vxi
        vyj = vyi
    return inside


if __name__ == "__main__":
    import doctest
    doctest.testmod()
    print "Doctests finished.\n"