# -*- coding: utf-8 -*-

"""
Functions dealing with shapes, especially points and polygons.
"""


def _is_a_vertex_of_polygon(x, y, polygon):
    """
    Check if the `x`/`y` coordinate is a vertex of the `polygon`.

        >>> polygon = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))

        >>> _is_a_vertex_of_polygon(0.0, 0.0, polygon)
        True

        >>> _is_a_vertex_of_polygon(0.5, 0.5, polygon)
        False

    Beware:
        - latitude is `y` and longitude is `x`
        - coordinate and `polygon` must use the same geodesic projection system
    """
    return (x, y) in polygon


def _is_within_bounding_box(x, y, minx, miny, maxx, maxy):
    """
    Check if the `x`/`y` coordinate is within the bounding box, defined
    by the `minx`/`miny` coordinate and the `maxx`/`maxy` coordinate.

        >>> _is_within_bounding_box(0.5, 0.5, 0.0, 0.0, 1.0, 1.0)
        True

        >>> _is_within_bounding_box(0.5, 0.5, 1.0, 1.0, 0.0, 0.0)
        True

        >>> _is_within_bounding_box(0.5, 0.5, 1.0, 0.0, 0.0, 1.0)
        True

        >>> _is_within_bounding_box(0.5, 0.5, 0.0, 1.0, 1.0, 0.0)
        True

        >>> _is_within_bounding_box(2.0, 2.0, 0.0, 0.0, 1.0, 1.0)
        False

        >>> _is_within_bounding_box(2.0, 2.0, 1.0, 1.0, 0.0, 0.0)
        False

        >>> _is_within_bounding_box(2.0, 2.0, 1.0, 0.0, 0.0, 1.0)
        False

        >>> _is_within_bounding_box(2.0, 2.0, 0.0, 1.0, 1.0, 0.0)
        False

    Beware:
        - latitude is `y` and longitude is `x`
        - All parameters must use the same geodesic projection system
    """
    return (
        min(minx, maxx) <= x <= max(minx, maxx) and
        min(miny, maxy) <= y <= max(miny, maxy)
    )


def _is_on_line(x, y, ax, ay, bx, by):
    """
    Check if point `x`/`y` is on the line segment from `ax`/`ay` to `bx`/`by`
    or the degenerate case that all 3 points are coincident.

        >>> _is_on_line(0.5, 0.5, 0.0, 0.0, 1.0, 1.0)
        True

        >>> _is_on_line(0.0, 0.0, 0.0, 0.0, 1.0, 1.0)
        True

        >>> _is_on_line(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
        True

        >>> _is_on_line(0.0, 0.5, 0.0, 0.0, 1.0, 1.0)
        False

        >>> _is_on_line(0.5, 0.0, 0.0, 0.0, 1.0, 1.0)
        False

    """
    return ((bx - ax) * (y - ay) == (x - ax) * (by - ay) and
            ((ax <= x <= bx or bx <= x <= ax) if ax != bx else
             (ay <= y <= by or by <= y <= ay)))


def _is_inside_polygon(x, y, polygon):
    """
    Check if the given `x`/`y` coordinate is inside the given `polygon`.

        >>> polygon = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))

        >>> _is_inside_polygon(0.0, 0.0, polygon)
        True

        >>> _is_inside_polygon(0.0, 0.5, polygon)
        ... # Attention: edge case, this should be True
        False

        >>> _is_inside_polygon(0.5, 0.0, polygon)
        ... # Attention: edge case, this should be True
        False

        >>> _is_inside_polygon(0.5, 0.5, polygon)
        True

        >>> _is_inside_polygon(1.0, 0.5, polygon)
        True

        >>> _is_inside_polygon(0.5, 1.0, polygon)
        True

        >>> _is_inside_polygon(1.0, 1.0, polygon)
        ... # Attention: edge case, this should be True
        False

        >>> _is_inside_polygon(2.0, 0.0, polygon)
        False

        >>> _is_inside_polygon(0.0, 2.0, polygon)
        False

        >>> _is_inside_polygon(-1.0, -1.0, polygon)
        False

    Beware:
        - latitude is `y` and longitude is `x`
        - coordinate and `polygon` must use the same geodesic projection system
        - this implementation (delibrately) fails for certain edge cases

    This function implements the ray casting algorithm.
    """

    # _infinity is used to act as infinity if we divide by zero
    _infinity = float('+inf')

    # _offset is used to make sure points are not on the same line as vertexes
    _offset = 0.00000001

    # used a few lines below
    _length = len(polygon)

    # we start on the outside of the polygon
    inside = False
    for index, point in enumerate(polygon):
        A = point
        B = polygon[(index + 1) % _length]

        # Make sure A is the lower point of the edge
        if A[1] > B[1]:
            A, B = B, A

        # Make sure point is not at same height as vertex
        if y == A[1]:
            A = (A[0], A[1] - _offset)
        elif y == B[1]:
            B = (B[0], B[1] + _offset)

        if (y > B[1] or y < A[1] or x > max(A[0], B[0])):
            # The horizontal ray does not intersect with the edge
            continue

        if x < min(A[0], B[0]):
            # The ray intersects with the edge
            inside = not inside
            continue

        try:
            edge = (B[1] - A[1]) / (B[0] - A[0])
        except ZeroDivisionError:
            edge = _infinity

        try:
            point = (y - A[1]) / (x - A[0])
        except ZeroDivisionError:
            point = _infinity

        if point >= edge:
            # The ray intersects with the edge
            inside = not inside
            continue

    return inside


def is_on_polygon(x, y, polygon, bounds=None):
    """
    Check if the given `x`/`y` coordinate is on the given `polygon`.
    The `bounds` can given be given as `(minx, miny, maxx, maxy)` to
    speed up the algorithm.

        >>> polygon = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))

        >>> is_on_polygon(0.0, 0.0, polygon)
        True

        >>> is_on_polygon(0.0, 0.5, polygon)
        True

        >>> is_on_polygon(0.5, 0.0, polygon)
        True

        >>> is_on_polygon(0.5, 0.5, polygon)
        True

        >>> is_on_polygon(1.0, 0.5, polygon)
        True

        >>> is_on_polygon(0.5, 1.0, polygon)
        True

        >>> is_on_polygon(1.0, 1.0, polygon)
        True

        >>> is_on_polygon(2.0, 0.0, polygon)
        False

        >>> is_on_polygon(0.0, 2.0, polygon)
        False

        >>> is_on_polygon(-1.0, -1.0, polygon)
        False

    Beware:
        - latitude is `y` and longitude is `x`
        - coordinate and `polygon` must use the same geodesic projection system
    """
    if _is_a_vertex_of_polygon(x, y, polygon):
        return True

    if bounds is None:
        minx, maxx = float('+inf'), float('-inf')
        miny, maxy = minx, maxx
        for point in polygon:
            minx = min(minx, point[0])
            miny = min(miny, point[1])
            maxx = max(maxx, point[0])
            maxy = max(maxy, point[1])
    else:
        minx, miny, maxx, maxy = bounds

    if not _is_within_bounding_box(x, y, minx, miny, maxx, maxy):
        return False

    if _is_inside_polygon(x, y, polygon):
        return True

    # used a few lines below
    _length = len(polygon)

    # make sure point is not an egde case from _is_inside_polygon
    for index, point in enumerate(polygon):
        A = point
        B = polygon[(index + 1) % _length]

        if _is_on_line(x, y, A[0], A[1], B[0], B[1]):
            return True

    return False
