sympy/geometry/util.py

   1
   2 def intersection(*entities):
   3     """
   4     Finds the intersection between a list GeometryEntity instances. Returns a
   5     list of all the intersections, Will raise a NotImplementedError exception
   6     if unable to calculate the intersection.
   7
   8     Examples:
   9     =========
  10         >>> from sympy.geometry import *
  11         >>> p1,p2,p3 = Point(0,0), Point(1,1), Point(-1, 5)
  12         >>> l1, l2 = Line(p1, p2), Line(p3, p2)
  13         >>> c = Circle(p2, 1)
  14         >>> intersection(l1, p2)
  15         [Point(1, 1)]
  16         >>> intersection(l1, l2)
  17         [Point(1, 1)]
  18         >>> intersection(c, p2)
  19         []
  20         >>> intersection(c, Point(1, 0))
  21         [Point(1, 0)]
  22         >>> intersection(c, l2)
  23         [Point(1 - 1/5*5**(1/2), 1 + 2*5**(1/2)/5), Point(1 + 1/5*5**(1/2), 1 - 2*5**(1/2)/5)]
  24
  25     Notes:
  26     ======
  27         - The intersection of any geometrical entity with itself should return
  28           a list with one item: the entity in question.
  29         - An intersection requires two or more entities. If only a single
  30           entity is given then one will receive an empty intersection list.
  31         - It is possible for intersection() to miss intersections that one
  32           knows exists because the required quantities were not fully
  33           simplified internally.
  34     """
  35     from entity import GeometryEntity
  36
  37     entities = GeometryEntity.extract_entities(entities, False)
  38     if len(entities) <= 1: return []
  39
  40     res = GeometryEntity.do_intersection(entities[0], entities[1])
  41     for entity in entities[2:]:
  42         newres = []
  43         for x in res:
  44             newres.extend( GeometryEntity.do_intersection(x, entity) )
  45         res = newres
  46     return res
  47
  48
  49 def convex_hull(*args):
  50     """
  51     Returns a Polygon representing the convex hull of a set of 2D points.
  52
  53     Notes:
  54     ======
  55         This can only be performed on a set of non-symbolic points.
  56
  57     Example:
  58     ========
  59         >>> from sympy.geometry import Point
  60         >>> points = [ Point(x) for x in [(1,1), (1,2), (3,1), (-5,2), (15,4)] ]
  61         >>> convex_hull(points)
  62         Polygon(Point(3, 1), Point(15, 4), Point(-5, 2), Point(1, 1))
  63
  64     Description of method used:
  65     ===========================
  66         See http://en.wikipedia.org/wiki/Graham_scan.
  67     """
  68     from point import Point
  69     from line import Segment
  70     from polygon import Polygon
  71
  72     p = args[0]
  73     if isinstance(p, Point):
  74         p = args
  75
  76     # Basic checks
  77     if len(p) == 1:
  78         return p[0]
  79     elif len(p) == 2:
  80         return Segment(p[0], p[1])
  81
  82     # Find lowest+rightmost point
  83     m = 0
  84     for i in xrange(1, len(p)):
  85         if (p[i][1] < p[m][1]) or ((p[i][1] == p[m][1]) and (p[i][0] > p[m][0])):
  86             m = i
  87     p[0], p[m] = p[m], p[0]
  88
  89     def tarea(a, b, c):
  90         return (b[0] - a[0])*(c[1] - a[1]) - (c[0] - a[0])*(b[1] - a[1])
  91
  92     # Radial sort of points with respect to p[0] (our pivot)
  93     destroy = {}
  94     p0 = p[0]
  95     def pcompare(p1, p2):
  96         a = tarea(p0, p1, p2)
  97         if a > 0:
  98             return -1
  99         elif a < 0:
 100             return 1
 101         else:
 102             x = abs(p1[0] - p0[0]) - abs(p2[0] - p0[0])
 103             y = abs(p1[1] - p0[1]) - abs(p2[1] - p0[1])
 104             if (x < 0) or (y < 0):
 105                 destroy[p1] = True
 106                 return -1
 107             elif (x > 0) or (y > 0):
 108                 destroy[p2] = True
 109                 return 1
 110             else:
 111                 destroy[p1] = True
 112                 return 0
 113     p = p[1:]
 114     p.sort(pcompare)
 115     p.insert(0, p0)
 116
 117     # Destroy points as found by sorting
 118     for i in xrange(len(p)-1, -1, -1):
 119         if p[i] in destroy:
 120             del p[i]
 121
 122     # Graham scan
 123     def isleft(a, b, c):
 124         return (tarea(a, b, c) > 0)
 125
 126     top = [p[0], p[1]]
 127     i = 2
 128     while i < len(p):
 129         p1 = top[-2]
 130         p2 = top[-1]
 131         if isleft(p1, p2, p[i]):
 132             top.append(p[i])
 133             i += 1
 134         else:
 135             top.pop()
 136     return Polygon(top)
 137
 138 def are_similar(e1, e2):
 139     """
 140     Returns True if e1 and e2 are similar (one can be uniformly scaled to
 141     the other) or False otherwise.
 142
 143     Notes:
 144     ======
 145         - If the two objects are equal then they are always similar.
 146     """
 147     if e1 == e2: return True
 148     try:
 149         return e1.is_similar(e2)
 150     except AttributeError:
 151         try:
 152             return e2.is_similar(e1)
 153         except AttributeError:
 154             n1 = e1.__class__.__name__
 155             n2 = e2.__class__.__name__
 156             raise GeometryError("Cannot test similarity between %s and %s" % (n1, n2))