Source

geoalg / geoalg / convexhull.py

Full commit
#!/usr/bin/env python
#coding:utf-8
# Author:  mozman
# Purpose: convex hull algorithm
# Created: 28.02.2010

from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import

__all__ = ['ConvexHull']

class ConvexHull(object):
    def __init__(self, points):
        self._points = ConvexHull._construct(points)

    def __iter__(self):
        return iter(self._points)

    def values(self):
        return self._points[:] # list shallow copy

    @staticmethod
    def _construct(points):
        def convex_hull(hull):
            while len(hull) > 2:
                start_point, check_point, destination_point = hull[-3:] # the last three points
                if not left_of_line(check_point, start_point, destination_point): # curve not turns right
                    del hull[-2] # remove the penultimate point
                else:
                    break
            return hull

        points = sorted(set(points)) #remove duplicate points

        if len(points) < 3:
            raise ValueError("ConvexHull(): Less than 3 unique points given!")

        upper_hull = points[:2] # first two points
        for next_point in points[2:]:
            upper_hull.append(next_point)
            upper_hull = convex_hull(upper_hull)
        lower_hull = [points[-1], points[-2]] # last two points

        for next_point in reversed(points[:-2]):
            lower_hull.append(next_point)
            lower_hull = convex_hull(lower_hull)
        upper_hull.extend(lower_hull[1:-1])
        return upper_hull

def left_of_line(point, p1, p2):
    """ True if the point self is left of the line p1 -> p2
    """
    # check if a and b are on the same vertical line
    if p1[0] == p2[0]:
        # compute # on which site of the line self should be
        should_be_left = p1[1] < p2[1]
        if should_be_left:
            return point[0] < p1[0]
        else:
            return point[0] > p1[0]
    else:
        # get pitch of line
        pitch = (p2[1] - p1[1]) / (p2[0] - p1[0])

        # get y-value at c's x-position
        y = pitch * (point[0] - p1[0]) + p1[1]

        # compute if point should be above or below the line
        should_be_above = p1[0] < p2[0]
        if should_be_above :
            return point[1] > y
        else:
            return point[1] < y