import numpy as np


def approximate_polygon(coords, tolerance, target_count):
    """
    Approximate a polygonal chain with the specified tolerance
    or until target count of points is attained.
    It is based on the Douglas-Peucker algorithm.
    Note that the approximated polygon is always within the convex hull of the
    original polygon.
    Parameters
    ----------
    coords : (N, 2) array
        Coordinate array.
    tolerance : float
        Maximum distance from original points of polygon to approximated
        polygonal chain. If tolerance is 0, the original coordinate array
        is returned.
    target_count: int
        Maximum count of polygon points.
    satisfy: {all, any}
        Simplify until all constraints (tolerance and target_count)
        or if any of them are fulfilled.
        Default `all`

    Returns
    -------
    coords : (M, 2) array
        Approximated polygonal chain where M <= N.

    References
    ----------
    .. [1] http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
    """

    validTolerance = tolerance >= 0
    validCount = target_count < len(coords)

    if not validTolerance and not validCount:
        return coords
    elif not validTolerance and validCount:
        # TODO: Calculate based on target_count
        return coords
    elif validTolerance and not validCount:
        # TODO: Calculate based on tolerance
        return coords
    elif validTolerance and validCount:
        # TODO: Calculate which of tolerance and target count converge first
        return coords

    chain = np.zeros(coords.shape[0], 'bool')
    # pre-allocate distance array for all points
    dists = np.zeros(coords.shape[0])
    chain[0] = True
    chain[-1] = True
    pos_stack = [(0, chain.shape[0] - 1)]
    end_of_chain = False

    while not end_of_chain:
        start, end = pos_stack.pop()
        # determine properties of current line segment
        r0, c0 = coords[start, :]
        r1, c1 = coords[end, :]
        dr = r1 - r0
        dc = c1 - c0
        segment_angle = - np.arctan2(dr, dc)
        segment_dist = c0 * np.sin(segment_angle) + r0 * np.cos(segment_angle)

        # select points in-between line segment
        segment_coords = coords[start + 1:end, :]
        segment_dists = dists[start + 1:end]

        # check whether to take perpendicular or euclidean distance with
        # inner product of vectors

        # vectors from points -> start and end
        dr0 = segment_coords[:, 0] - r0
        dc0 = segment_coords[:, 1] - c0
        dr1 = segment_coords[:, 0] - r1
        dc1 = segment_coords[:, 1] - c1
        # vectors points -> start and end projected on start -> end vector
        projected_lengths0 = dr0 * dr + dc0 * dc
        projected_lengths1 = - dr1 * dr - dc1 * dc
        perp = np.logical_and(projected_lengths0 > 0,
                              projected_lengths1 > 0)
        eucl = np.logical_not(perp)
        segment_dists[perp] = np.abs(
            segment_coords[perp, 0] * np.cos(segment_angle)
            + segment_coords[perp, 1] * np.sin(segment_angle)
            - segment_dist
        )
        segment_dists[eucl] = np.minimum(
            # distance to start point
            np.sqrt(dc0[eucl] ** 2 + dr0[eucl] ** 2),
            # distance to end point
            np.sqrt(dc1[eucl] ** 2 + dr1[eucl] ** 2)
        )

        if np.any(segment_dists > tolerance):
            # select point with maximum distance to line
            new_end = start + np.argmax(segment_dists) + 1
            pos_stack.append((new_end, end))
            pos_stack.append((start, new_end))
            chain[new_end] = True

        if len(pos_stack) == 0:
            end_of_chain = True

    return coords[chain, :]
