What is most efficient way to find the intersection of a line and a circle in python?

Question

I have a polygon consists of lots of points. I want to find the intersection of the polygon and a circle. Providing the circle center of [x0,y0] and radius of r0, I have wrote a rough function to simply solve the quadratic equation of the circle and a line. But what about the efficiency of find the intersection of every line segment of the polygon one by one? Is there more efficient way?

I know sympy already provide the feature to get the intersections of different geometry. But also what about the efficiency of external library like sympy compared to calculate it by my own function, if I want to deal with lots of polygons?

def LineIntersectCircle(p,lsp,lep):
# p is the circle parameter, lsp and lep is the two end of the line
  x0,y0,r0 = p
  x1,y1 = lsp
  x2,y2 = esp
  if x1 == x2:
    if abs(r0) >= abs(x1 - x0):
        p1 = x1, y0 - sqrt(r0**2 - (x1-x0)**2)
        p2 = x1, y0 + sqrt(r0**2 - (x1-x0)**2)
        inp = [p1,p2]
        # select the points lie on the line segment
        inp = [p for p in inp if p[1]>=min(y1,y2) and p[1]<=max(y1,y2)]
    else:
        inp = []
  else:
    k = (y1 - y2)/(x1 - x2)
    b0 = y1 - k*x1
    a = k**2 + 1
    b = 2*k*(b0 - y0) - 2*x0
    c = (b0 - y0)**2 + x0**2 - r0**2
    delta = b**2 - 4*a*c
    if delta >= 0:
        p1x = (-b - sqrt(delta))/(2*a)
        p2x = (-b + sqrt(delta))/(2*a)
        p1y = k*x1 + b0
        p2y = k*x2 + b0
        inp = [[p1x,p1y],[p2x,p2y]]
        # select the points lie on the line segment
        inp = [p for p in inp if p[0]>=min(x1,x2) and p[0]<=max(x1,x2)]
    else:
        inp = []
  return inp

Answer 1

I guess maybe your question is about how to in theory do this in the fastest manner. But if you want to do this quickly, you should really use something which is written in C/C++.

I am quite used to Shapely, so I will provide an example of how to do this with this library. There are many geometry libraries for python. I will list them at the end of this answer.

from shapely.geometry import LineString
from shapely.geometry import Point

p = Point(5,5)
c = p.buffer(3).boundary
l = LineString([(0,0), (10, 10)])
i = c.intersection(l)

print i.geoms[0].coords[0]
(2.8786796564403576, 2.8786796564403576)

print i.geoms[1].coords[0]
(7.121320343559642, 7.121320343559642)

What is a little bit counter intuitive in Shapely is that circles are the boundries of points with buffer areas. This is why I do p.buffer(3).boundry

Also the intersection i is a list of geometric shapes, both of them points in this case, this is why I have to get both of them from i.geoms[]

There is another Stackoverflow question which goes into details about these libraries for those interested.

• SymPy
• CGAL Python bindings
• PyEuclid
• PythonOCC
• Geometry-simple

EDIT because comments:

Shapely is based on GEOS (trac.osgeo.org/geos) which is built in C++ and considerably faster than anything you write natively in python. SymPy seems to be based on mpmath (mpmath.org) which also seems to be python, but seems to have lots of quite complex math integrated into it. Implementing that yourself may require a lot of work, and will probably not be as fast as GEOS C++ implementations.