Adds a geometry module with a segment_intersection function

geometry.py

@@ -0,0 +1,208 @@
+#!/usr/bin/env python
+
+from __future__ import division 
+import math
+
+epsilon = 1e-6
+
+def x( point ):
+    return point[0]
+
+def y( point ):
+    return point[1]
+
+
+def mid( xy, pa, pb ):
+    return ( xy(pa) + xy(pb) ) / 2.0
+
+def middle( pa, pb ):
+    return mid(x,pa,pb),mid(y,pa,pb)
+
+
+def euclidian_distance( ci, cj, graph = None):
+    return math.sqrt( float(ci[0] - cj[0])**2 + float(ci[1] - cj[1])**2 )
+
+
+def linear_equation( p0, p1 ):
+    """Return the linear equation coefficients of a line given by two points.
+       Use the general form: c=a*x+b*y """
+    assert( len(p0) == 2 )
+    assert( len(p1) == 2 )
+
+    a = y(p0) - y(p1)
+    b = x(p1) - x(p0)
+    c = x(p0) * y(p1) - x(p1) * y(p0)
+    return a, b, -c
+
+
+def is_null( x, e = epsilon ):
+    return -e <= x <= e
+
+
+def is_vertical( leq ):
+    a,b,c = leq
+    return is_null(b)
+
+
+def is_point( segment ):
+    """Return True if the given segment is degenerated (i.e. is a single point)."""
+    return segment[0] == segment[1]
+
+
+def collinear( p, q, r ):
+    """Returns True if the 3 given points are collinear.
+
+    Note: there is a lot of algorithm to test collinearity, the most known involving linear algebra.
+    This one has been found in Jonathan Shewchuk's "Lecture Notes on Geometric Robustness".
+    It is maybe the most elegant one: just arithmetic on x and y, without ifs, sqrt or risk of divide-by-zero error.
+    """
+    return (x(p)-x(r)) * (y(q)-y(r)) == (x(q)-x(r)) * (y(p)-y(r))
+
+
+def line_intersection( seg0, seg1 ):
+    """Return the coordinates of the intersection point of two lines given by pairs of points, or None."""
+
+    # Degenerated segments
+    def on_line(p,seg):
+        if collinear(p,*seg):
+            return p
+        else:
+            return None
+
+    # Segments degenerated as a single points,
+    if seg0[0] == seg0[1] == seg1[0] == seg1[1]:
+        return seg0[0]
+    # as two different points,
+    elif is_point(seg0) and is_point(seg1) and seg0 != seg1:
+        return None
+    # as a point and a line.
+    elif is_point(seg0) and not is_point(seg1):
+        return on_line(seg0[0],seg1)
+    elif is_point(seg1) and not is_point(seg0):
+        return on_line(seg1[0],seg0)
+
+
+    leq0 = linear_equation(*seg0)
+    leq1 = linear_equation(*seg1)
+
+    # Collinear lines.
+    if leq0 == leq1:
+        return None
+
+    # Vertical line
+    def on_vertical( seg, leq ):
+        a,b,c = leq
+        assert( not is_null(b) )
+        assert( is_null( x(seg[0])-x(seg[1]) ) )
+        px = x(seg[0])
+        # Remember that the leq form is c=a*x+b*y, thus y=(c-ax)/b
+        py = (c-a*px)/b
+        return px,py
+
+    if   is_vertical(leq0) and not is_vertical(leq1):
+        return on_vertical( seg0, leq1 )
+    elif is_vertical(leq1) and not is_vertical(leq0):
+        return on_vertical( seg1, leq0 )
+    elif is_vertical(leq0) and is_vertical(leq1):
+        return None
+
+    # Generic case.
+    a0,b0,c0 = leq0
+    a1,b1,c1 = leq1
+
+    d  = a0 * b1 - b0 * a1
+    dx = c0 * b1 - b0 * c1
+    dy = a0 * c1 - c0 * a1
+    if not is_null(d):
+        px = dx / d
+        py = dy / d
+        return px,py
+    else:
+        # Parallel lines
+        return None
+
+
+def box( points ):
+    """Return the min and max points of the bounding box enclosing the given set of points."""
+    minp = min( [ x(p) for p in points ] ), min( [ y(p) for p in points ] )
+    maxp = max( [ x(p) for p in points ] ), max( [ y(p) for p in points ] )
+    return minp,maxp
+
+
+def in_box( point, box, exclude_edges = False ):
+    """Return True if the given point is located within the given box."""
+    pmin,pmax = box
+    if exclude_edges:
+        return x(pmin)-epsilon <  x(point) <  x(pmax)+epsilon and y(pmin)-epsilon <  y(point) <  y(pmax)+epsilon
+    else:
+        return x(pmin)-epsilon <= x(point) <= x(pmax)+epsilon and y(pmin)-epsilon <= y(point) <= y(pmax)+epsilon
+
+
+def segment_intersection( seg0, seg1 ):
+    """Return the coordinates of the intersection point of two segments, or None.
+    If segments are colinear, returns colinear_value."""
+    assert( len(seg0) == 2 )
+    assert( len(seg1) == 2 )
+
+    p = line_intersection(seg0,seg1)
+
+    if p is None:
+        return None
+    else:
+        if in_box(p,box(seg0)) and in_box(p,box(seg1)):
+            return p
+        else:
+            return None
+
+
+if __name__ == "__main__":
+    import sys
+    import random
+    import utils
+    import uberplot
+    import matplotlib.pyplot as plot
+
+    if len(sys.argv) > 1:
+        scale = 100
+        nb = int(sys.argv[1])
+        points = [ (scale*random.random(),scale*random.random()) for i in range(nb)]
+    else:
+        points = [
+                (10,0),
+                (-190,0),
+                (10,200),
+                (110,100),
+                (110,-100),
+                (-90,100),
+                (-90,-100),
+                ]
+
+    segments = []
+    for p0 in points:
+        for p1 in points:
+            if p0 != p1:
+                segments.append( (p0,p1) )
+
+    seg_inter = []
+    line_inter = []
+    for s0 in segments:
+        for s1 in segments:
+            if s0 != s1:
+                s = segment_intersection( s0, s1 )
+                if s is not None:
+                    seg_inter.append(s)
+                l = line_intersection( s0, s1 )
+                if l is not None:
+                    line_inter.append(l)
+
+
+    fig = plot.figure()
+
+    ax = fig.add_subplot(111)
+    ax.set_aspect('equal')
+    uberplot.plot_segments( ax, segments, linewidth=0.5, edgecolor = "blue" )
+    uberplot.scatter_points( ax, points,     edgecolor="blue", facecolor="blue",  s=120, alpha=1, linewidth=1 )
+    uberplot.scatter_points( ax, line_inter, edgecolor="none", facecolor="green", s=60,  alpha=0.5 )
+    uberplot.scatter_points( ax, seg_inter,  edgecolor="none", facecolor="red",   s=60,  alpha=0.5 )
+    plot.show()
+

hull.py

@@ -1,6 +1,7 @@
 
 import operator
-from utils import x,y,euclidian_distance,LOG,LOGN
+from utils import LOG,LOGN
+from geometry import x,y,euclidian_distance
 
 # Based on the excellent article by Tom Switzer <thomas.switzer@gmail.com>
 # http://tomswitzer.net/2010/12/2d-convex-hulls-chans-algorithm/

shortpath.py

@@ -1,13 +1,6 @@
 
 import math
-
-def euclidian_distance(node, goal):
-    def x(node):
-        return node[0]
-    def y(node):
-        return node[1]
-
-    return math.sqrt( (x(node) - x(goal))**2 + (y(node) - y(goal))**2)
+from geometry import x,y,euclidian_distance
 
 
 def astar(graph, start, goal, cost = euclidian_distance, heuristic = euclidian_distance):

triangulation.py

@@ -1,9 +1,11 @@
 
 import sys
 import math
-from utils import tour,LOG,LOGN,x,y
 from itertools import ifilterfalse as filter_if_not
 
+from utils import tour,LOG,LOGN,x,y
+from geometry import mid, middle
+
 # Based on http://paulbourke.net/papers/triangulate/
 # Efficient Triangulation Algorithm Suitable for Terrain Modelling
 #     An Algorithm for Interpolating Irregularly-Spaced Data
@@ -12,12 +14,6 @@ from itertools import ifilterfalse as filter_if_not
 # Presented at Pan Pacific Computer Conference, Beijing, China.
 # January 1989
 
-def mid( xy, pa, pb ):
-    return ( xy(pa) + xy(pb) ) / 2.0
-
-def middle( pa, pb ):
-    return mid(x,pa,pb),mid(y,pa,pb)
-
 def mtan( pa, pb ):
     return -1 * ( x(pa) - x(pb) ) / ( y(pa) - y(pb) )
 
@@ -97,7 +93,7 @@ def in_circumcircle( p, triangle, epsilon  = sys.float_info.epsilon ):
     return in_circle( p, (cx,cy), r, epsilon )
 
 
-def in_triangle( p0, triangle, exclude_edges = True ):
+def in_triangle( p0, triangle, exclude_edges = False ):
     """Return True if the given point lies inside the given triangle"""
 
     p1,p2,p3 = triangle

utils.py

@@ -1,13 +1,7 @@
 
 import sys
 import math
-
-def x( point ):
-    return point[0]
-
-def y( point ):
-    return point[1]
-
+from geometry import x,y
 
 def LOG( *args ):
     """Print something on stderr and flush"""
@@ -84,7 +78,7 @@ def adjacency_from_set( segments ):
     return graph
 
 
-def vertices_from_set( segments ):
+def vertices_of( segments ):
     vertices = set()
     for start,end in segments:
         vertices.add(start)
@@ -98,6 +92,3 @@ def tour(lst):
         yield (a,b)
 
 
-def euclidian_distance( ci, cj, graph = None):
-    return math.sqrt( float(ci[0] - cj[0])**2 + float(ci[1] - cj[1])**2 )
-