Demonstrates algorithms for solving several well-known problems in computational geometry including:1. Determining whether a polygon is convex2. Determining whether an ordered list of points represent a simple polygon3. Finding intersection points of line segments4. Finding the convex hull of a set of points5. Finding the minimum enclosing disc …

Comp_Geometry.py

__author__ = 'Meir'

import random
import numpy as np
import pylab as pl
from matplotlib import collections as mc
from itertools import combinations, permutations

N = 50 # number of points

def r(): # returns a random number
    return random.randint(0,400)

def p(): # returns a random point
    return (r(),r())

def s(): # returns a random line segment
    return (p(), p())

def normalize(vec):
    """converts a vector into a unit vector"""
    return [i/(sum([j ** 2 for j in vec]) ** .5) for i in vec]

def dist(a, b):
    """Computes the distance between two points"""
    return (sum([(a[i] - b[i]) ** 2 for i in range(len(a))]) ** .5)

def duplicates_removed(lis):
    """is like set(), but preserves order. keeps the first of each element"""
    new_list = []
    for i in lis:
        if i not in new_list:
            new_list.append(i)
    return new_list

def determinant(a,b,c,d,e,f,g,h,i):
    """
    returns the determinant

    | a b c |
    | d e f |
    | g h i |
    """
    return a*e*i + b*f*g + c*d*h - b*d*i - a*f*h - c*e*g

def ccw(p1, p2, p3):
    """Checks whether a list of 3 points in order have counter-clockwise orientation"""
    x1, y1, x2, y2, x3, y3 = p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]
    return ((x2-x1)*(y3-y2) > (x3-x2)*(y2-y1))

def ccw_or_collinear(p1,p2,p3):
    x1, y1, x2, y2, x3, y3 = p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]
    return ((x2-x1)*(y3-y2) >= (x3-x2)*(y2-y1))

def collinear(p1,p2,p3):
    x1, y1, x2, y2, x3, y3 = p1[0], p1[1], p2[0], p2[1], p3[0], p3[1]
    return ((x2-x1)*(y3-y2) == (x3-x2)*(y2-y1))

def opp(point): #returns point reflected over the line y = x
    return (point[1], point[0])

def midpoint(a, b):
    return ((a[0]+b[0])/2, (a[1]+b[1])/2)

def intersects(seg1, seg2):
    p1, p2, p3, p4 = seg1[0], seg1[1], seg2[0], seg2[1]
    if len(list(set([p1,p2,p3,p4]))) != 4: # if the segments share a common endpoint
        return False
    return (ccw(p1, p2, p3) != ccw(p1, p2, p4)) and (ccw(p3,p4,p1) != ccw(p3,p4,p2))

def circumcenter(a,b,c):
    """returns the center of the circle which contains a, b, and c"""
    if collinear(a,b,c):
        print("The points are collinear. No circumcircle exists.")
        return None
    if a[1] == b[1]: #to avoid division by zero
        b,c = c,b
    v0x, v0y = (b[0]-a[0], b[1]-a[1]) # vector from a to b
    v1x, v1y = (c[0]-b[0], c[1]-b[1]) # vector from b to c
    x0, y0 = midpoint(a,b)
    x1, y1 = midpoint(b,c)
    s = (y0 + (v0x*x0/v0y) - (v0x*x1/v0y) - y1) / (v1x - (v0x*v1y/v0y))
    center = (x1 - v1y*s, y1 + v1x*s)
    return center

def circumcircle(a,b,c):
    center = circumcenter(a,b,c)
    radius = dist(a,center)
    return (center, radius)

def convex(points):
    # points is a list of tuples
    full_cycle = points[:] + points[:2]
    # make sure that either all consecutive triples of points are either ccw, or all are cw
    return len(list(set([ccw_or_collinear(full_cycle[i], full_cycle[i+1], full_cycle[i+2]) for i in range(len(full_cycle)-2)]))) == 1

def simple(points):
    """Checks whether a list of points represents a simple polygon"""
    full_cycle = points[:] + [points[0]]
    segs = [(full_cycle[i], full_cycle[i+1]) for i in range(len(full_cycle)-1)]
    combs = []
    for comb in combinations(segs,2):
        combs.append(comb)
    for comb in combs:
        if len(list(set((comb[0][0], comb[0][1], comb[1][0], comb[1][1])))) == 4 and intersects(comb[0],comb[1]):
            #print(comb)
            return False
    return True

def dot_product(v1,v2):
    """returns the dot product of two vectors"""
    return sum([v1[i]*(v2[1]) for i in range(len(v1))])

def cross_product(v1,v2):
    """returns the magnitude of the cross product of two vectors.
    sign denotes whether positive or negative z-direction"""
    return v1[0]*v2[1] - v2[0]*v1[1]

def polygon_area(polygon): #polygon is an ordered list of points
    assert simple(polygon)
    area = 0
    n = len(polygon)
    for i in range(1,n-1):
        x1 = polygon[i][0] - polygon[0][0]
        y1 = polygon[i][1] - polygon[0][1]
        x2 = polygon[i+1][0] - polygon[0][0]
        y2 = polygon[i+1][1] - polygon[0][1]
        cross = x1*y2 - x2*y1
        area += cross / 2
        print(cross / 2)
    print()
    return abs(area)

def my_dot_product(p, a):
    """takes the dot product of the unit vector p -> a and (0,0) -> (1,0)"""
    v = normalize((a[0]-p[0], a[1] - p[1]))
    return v[0]

def convex_hull(points):
    """Graham's scan"""
    points = list(set(points)) #is now unique; order not preserved but that doesn't matter
    p = min(points, key=(lambda point: (point[1], point[0])))
    high_points = points
    high_points.remove(p)
    #now sort them according to the angle they and p make with the positive x-axis and their distance from p
    #priority is in that order
    high_points = list(sorted(high_points, key=lambda point:
    (-round(my_dot_product(p, point),7), dist(p, point))))
    # for hp in high_points:
    #     print(hp, "dot product: "+ str(my_dot_product(p, hp)))
    hull = [p, high_points[0]]
    counter = 1 # start at second high point
    while counter < len(high_points): # allows more flexibility than a for loop (might use for animated graphing)
        hull.append(high_points[counter])
        while len(hull) >= 3 and (not ccw(hull[-3], hull[-2], hull[-1])): #must be in this order! otherwise ccw might be called when len(hull) < 3
            hull = hull[:-2] + [hull[-1]]
        counter += 1
    if not ccw(hull[-2], hull[-1], p):
        hull = hull[:-1]
    assert convex(hull)
    assert simple(hull)
    return hull

def slow_triangulation_with_skinnies(points): #naive, adds segments as quickly and easily as it can to the triangulation
    points = list(set(points)) # order not preserved, doesn't matter
    print("points:", points)
    triangulation = []
    print(convex_hull(points))
    h = len(convex_hull(points)) # number of points in hull
    print('h:', h)
    v = len(points) # number of points
    print('v:', v)
    e = 3*v - h - 3 # number of edges in triangulation. I have not proven it, just experimented many times
    print("e:", e)
    # while len(triangulation) < e:
    #     print('starting...')
    # here we don't need e
    for p1 in points:
        for p2 in points:
            if p1 != p2:
                if not any([intersects((p1,p2),edge) for edge in triangulation]):
                    if not ((p1,p2) in triangulation) and not ((p2,p1) in triangulation):
                        triangulation.append((p1,p2))
    return triangulation

def slow_triangulation(points): #smarter, tries to minimize edge length
    """Until the triangulation is complete, go through the points sequentially, and for a point p,
    and add an edge to the triangulation consisting of p and the viable point closest to it. A point is viable if
    it does not create an intersection or an edge already in the triangulation.
    """
    points = duplicates_removed(points)
    triangulation = []
    h = len(convex_hull(points)) # number of points in hull
    v = len(points) # number of points
    e = 3*v - h - 3 # number of edges in triangulation. found by experimentation, verified mathematically
    index = 0
    while len(triangulation) < e:
        p = points[index]
        other_points = points[:]
        other_points.remove(p)
        other_points = sorted(other_points, key=lambda other_point: dist(other_point, p))
        other_points_index = 0
        while other_points_index < len(other_points) and \
        (any([intersects((p,other_points[other_points_index]),edge) for edge in triangulation]) or \
        ((p,other_points[other_points_index]) in triangulation) or \
        ((other_points[other_points_index],p) in triangulation)): # 'and' expression must be ordered to avoid IndexError; god i hate doing that
            other_points_index += 1
        if other_points_index < len(other_points):
            triangulation.append((p,other_points[other_points_index]))
        index = (index+1) % v
    return triangulation

def almostEquals(a,b, epsilon=.0001):
    return abs(a-b) < epsilon

def contains(circle, point):
    """Checks whether a circle in (center, radius) form contains a point"""
    return dist(circle[0], point) <=  circle[1] or almostEquals(dist(circle[0], point),circle[1])

def random_permutation(elements):
    """like random.shuffle() but functional instead of in-place"""
    new_elements = elements[:]
    random.shuffle(new_elements)
    return new_elements

def mini_disc(points):
    """Computes the minimum enclosing disc of a set of points. Returns (center, radius)"""
    points = list(set(points))
    points = random_permutation(points)
    n = len(points)
    D = (midpoint(points[0], points[1]), dist(points[0], points[1])/2)
    for i in range(2,n):
        if not contains(D, points[i]):
            D = mini_disc_with_point(points[:i], points[i])
    return D

def mini_disc_with_point(points, q):
    """Computes the smallest enclosing disc for a set of points with q on its boundary"""
    points = random_permutation(points)
    n = len(points)
    D = (midpoint(points[0], q), dist(points[0], q)/2)
    for j in range(1,n):
        if not contains(D, points[j]):
            D = mini_disc_with_2_points(points[:j], points[j], q)
    return D

def mini_disc_with_2_points(points, q1, q2):
    """Computes the smallest enclosing disc for a set of points with q1 and q2 on its boundary."""
    D = (midpoint(q1, q2), dist(q1, q2)/2)
    n = len(points)
    for k in range(n):
        if not contains(D, points[k]):
            D = circumcircle(q1,q2,points[k])
    return D

###################################################
###################################################
###################################################
#                                                 #
#                                                 #
#                                                 #
#                      TESTS                      #
#                                                 #
#                                                 #
#                                                 #
###################################################
###################################################
###################################################


def convexity_test():
    """Tests whether 'points' represents a convex polygon. Does not check simplicity."""
    points = [p() for i in range(3)]
    while len(points) < N:
        new_point = p()
        while not simple(points + [new_point]):
            new_point = p()
        points.append(new_point)
    points += [points[0]] # complete the cycle
    lines = []
    for i in range(len(points)-1):
        lines.append([points[i], points[i+1]])
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=2)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(0.1)
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.axis([min(i[0] for i in points) - 1.5, max(i[0] for i in points) + 1.5,
             min(i[1] for i in points) - 1.5, max(i[1] for i in points) + 1.5])
    ax.text(min([i[0] for i in points]) - 1, min(i[1] for i in points) - 1,
            "CONVEXITY: " + str(convex(points[:-1])), fontsize=15)  # remember to omit the cycle completer at the end
    pl.show()

def simplicity_test():
    points = [p() for i in range(N)]
    points += [points[0]] # complete the cycle
    lines = []
    for i in range(len(points)-1):
        lines.append([points[i], points[i+1]])
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=2)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(0.1)
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.axis([min(i[0] for i in points) - 1.5, max(i[0] for i in points) + 1.5,
             min(i[1] for i in points) - 1.5, max(i[1] for i in points) + 1.5])
    ax.text(min([i[0] for i in points]) - 1, min(i[1] for i in points) - 1,
            "SIMPLE: " + str(simple(points[:-1])), fontsize=15) # remember to omit the cycle completer at the end
    pl.show()

def intersection_test():
    segs = (s(),s())
    p1, p2, p3, p4 = segs[0][0], segs[0][1], segs[1][0], segs[1][1]
    segs = ((p1,p2),(p3,p4))
    points = (p1,p2,p3,p4)
    lines = []
    lines.append((p1,p2))
    lines.append((p3,p4))
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=2)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(0.1)
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.axis([min(i[0] for i in points) - 1.5, max(i[0] for i in points) + 1.5,
             min(i[1] for i in points) - 1.5, max(i[1] for i in points) + 1.5])
    ax.text(min([i[0] for i in points]) - 1, min(i[1] for i in points) - 1,
            "Intersection: " + str(intersects((p1,p2),(p3,p4))), fontsize=15)
    pl.show()

def hull_test():
    #random points within a square:
    #points = [p() for i in range(N)]
    #random points within a circle:
    points = []
    for i in range(N):
        x,y = (0,0)
        while (x-.5)**2 + (y-.5)**2 > .25:
            x, y = random.random(), random.random()
        points.append((x,y))
    #points = [(0,0),(-1,1),(0,1),(1,1),(-2,2),(-1,2),(0,2),(1,2),(2,2),(-3,3),(-2,3),(-1,3),(0,3),(1,3),(2,3),(3,3)] # diagonal
    #points = [(3,4),(7,4),(5,4),(8.5,4),(5,9),(10,10),(-10,10)] #horizontal
    print("points:",points)
    hull = convex_hull(points)
    print('hull:', hull)
    hull += [hull[0]] # complete the cycle
    lines = []
    for i in range(len(hull)-1):
        lines.append([hull[i], hull[i+1]])
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=2)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(.1)
    xrange = max(points, key=lambda point: point[0])[0] - min(points, key=lambda point: point[0])[0]
    yrange = max(points, key=lambda point: point[1])[1] - min(points, key=lambda point: point[1])[1]
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.plot([i[0] for i in hull], [i[1] for i in hull], 'ro')
    pl.axis([min(i[0] for i in points) - xrange/10, max(i[0] for i in points) + xrange/10,
             min(i[1] for i in points) - yrange/10, max(i[1] for i in points) + yrange/10])
    pl.show()

def circumcenter_test():
    #a, b, c, = p(), p(), p()
    a, b, c = (193,-48),(0,19),(58,-77)
    points = (a, b, c)
    if collinear(*points):
        print("The points are collinear. No circumcircle exists.")
        return
    print("points:", points)
    center = circumcenter(a, b, c)
    print("center:", center)
    radius = dist(center, a)
    print("radius:", radius)
    lines = [(point, center) for point in points]
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=2)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    circle = pl.Circle(center, radius=radius, alpha = .1)
    ax.add_patch(circle)
    ax.autoscale()
    ax.margins(0.1)
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.plot([center[0]], [center[1]], 'ro')
    margin = radius / 10
    pl.axis([center[0]-radius-margin, center[0]+radius+margin,
             center[1]-radius-margin, center[1]+radius+margin])
    pl.show()

def triangulation_test():
    #points = [p() for i in range(N)]
    #random points within a circle:
    points = []
    for i in range(N):
        x,y = (0,0)
        while (x-.5)**2 + (y-.5)**2 > .25:
            x, y = random.random(), random.random()
        points.append((x,y))

    #slightly staggered grid
    # for x in range(8):
    #     for y in range(8):
    #         delta = .5
    #         points.append((x+random.uniform(-delta,delta),y+random.uniform(-delta,delta)))

    triangulation_function = slow_triangulation#_with_skinnies
    lines = triangulation_function(points)
    col = np.array([(0, 0, 0, 1)])
    lc = mc.LineCollection(lines, colors=col, linewidths=1.5)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(.1)
    xrange = max(points, key=lambda point: point[0])[0] - min(points, key=lambda point: point[0])[0]
    yrange = max(points, key=lambda point: point[1])[1] - min(points, key=lambda point: point[1])[1]
    pl.plot([i[0] for i in points], [i[1] for i in points], 'ro')
    pl.axis([min(i[0] for i in points) - xrange/10, max(i[0] for i in points) + xrange/10,
             min(i[1] for i in points) - yrange/10, max(i[1] for i in points) + yrange/10])
    pl.show()

def mini_disc_test():
    points = [p() for i in range(N)]
    print("points:", points)
    center, radius = zach_mini_disc(points)
    print("center:", center)
    print("radius:", radius)
    col = np.array([(0, 0, 0, 1)])
    fig, ax = pl.subplots()
    circle = pl.Circle(center, radius=radius, alpha = .1)
    ax.add_patch(circle)
    ax.autoscale()
    ax.margins(0.1)
    pl.plot([i[0] for i in points], [i[1] for i in points], 'bo')
    pl.plot([center[0]], [center[1]], 'ro')
    margin = radius / 10
    scale = 1
    # xrange = max(points, key=lambda point: point[0])[0] - min(points, key=lambda point: point[0])[0]
    # yrange = max(points, key=lambda point: point[1])[1] - min(points, key=lambda point: point[1])[1]
    # pl.axis([min(i[0] for i in points) - xrange/10, max(i[0] for i in points) + xrange/10,
    #          min(i[1] for i in points) - yrange/10, max(i[1] for i in points) + yrange/10])
    pl.axis([center[0]-radius-margin*scale, center[0]+radius+margin*scale,
             center[1]-radius-margin, center[1]+radius+margin])
    pl.show()


#convexity_test()
#simplicity_test()
#intersection_test()
#hull_test()
#circumcenter_test()
#triangulation_test()
#mini_disc_test()
#print(intersects(((0,0),(0,9)), ((0,0),(0,7))))