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

Familial

seems to me that you consider the intersections of the circle with the whole line, not only the line segment between the two given points. Is this what you want? – Fleabitten

The only possible optimisation I can think of involves the use of some kind of spatial partition. E.g. a quad-tree. But there is non-trivial cost associated in computing these, so it depends on your larger problem if that's useful or not. – Hebephrenia

@EmanuelePaolini,Thank you.I've modified the script according to your concern. – Familial

#1073836 - may be relevant – Asseverate

Are there still any open issues not answered by any of the answers posted? – Pellucid

You need to fix the variable "esp" in the line 5. I assume it is "lep" as the input function. – Candicecandid

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.

Pellucid

Yeah...What about the efficiency of calculate the intersection by the library like Shapely or SymPy, compared to define function myself? That's one of what I want to know. – Familial

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. – Pellucid

Sometimes this returns just a single Point, and sometimes MultiplePoints. How can I test the return value to know which is which? – Labarum

@Labarum a straight line can intersect a circle 0-2 times. This is why you get sometimes a single point or a multipoint. It would be pythonic to assume single and if you get an exception you fallback to multi. But I don't know what you want to test. – Pellucid

I was looking to see how to detect which I get. Turns out I can do an isinstance() to check the class. (I was hoping for a method on the object proper) – Labarum

@Labarum I think you can always instantiate it as a multipoint and treat it as a multipoint. just take the result and hand it to MultiPoint(result) – Pellucid

If I change the first point of the LineString to (6,5) or (5,5), the output gives an error stating that the" 'Point' object has no attribute 'geoms' " . Does it mean that the line must start and end outside the circle? Could you please elaborate? – Charil

Since somebody in the comments asked: I benchmarked the python function proposed by Peter below with the shapely library as proposed in this answer python 0.000019 seconds. shapely 0.000101 seconds. The native python implementation is 5 times faster on my computer. This is primarily due to the overhead of creating the shapely objects (~70 mu) but the raw intersection computation is still ~15 mu faster in python. – Shrew

@BernhardJaeger Thanks for checking and providing us with your findings. I think it should be noted that creating the shapely objects can be done once, then you can reuse them to check intersection many times, at which point, using shapely will be faster. But for a one time check, native python can be faster. There are usually many edge cases in geometric math that needs to be accounted for. – Pellucid

Answer 2

Here's a solution that computes the intersection of a circle with either a line or a line segment defined by two (x, y) points:

def circle_line_segment_intersection(circle_center, circle_radius, pt1, pt2, full_line=True, tangent_tol=1e-9):
    """ Find the points at which a circle intersects a line-segment.  This can happen at 0, 1, or 2 points.

    :param circle_center: The (x, y) location of the circle center
    :param circle_radius: The radius of the circle
    :param pt1: The (x, y) location of the first point of the segment
    :param pt2: The (x, y) location of the second point of the segment
    :param full_line: True to find intersections along full line - not just in the segment.  False will just return intersections within the segment.
    :param tangent_tol: Numerical tolerance at which we decide the intersections are close enough to consider it a tangent
    :return Sequence[Tuple[float, float]]: A list of length 0, 1, or 2, where each element is a point at which the circle intercepts a line segment.

    Note: We follow: http://mathworld.wolfram.com/Circle-LineIntersection.html
    """

    (p1x, p1y), (p2x, p2y), (cx, cy) = pt1, pt2, circle_center
    (x1, y1), (x2, y2) = (p1x - cx, p1y - cy), (p2x - cx, p2y - cy)
    dx, dy = (x2 - x1), (y2 - y1)
    dr = (dx ** 2 + dy ** 2)**.5
    big_d = x1 * y2 - x2 * y1
    discriminant = circle_radius ** 2 * dr ** 2 - big_d ** 2

    if discriminant < 0:  # No intersection between circle and line
        return []
    else:  # There may be 0, 1, or 2 intersections with the segment
        intersections = [
            (cx + (big_d * dy + sign * (-1 if dy < 0 else 1) * dx * discriminant**.5) / dr ** 2,
             cy + (-big_d * dx + sign * abs(dy) * discriminant**.5) / dr ** 2)
            for sign in ((1, -1) if dy < 0 else (-1, 1))]  # This makes sure the order along the segment is correct
        if not full_line:  # If only considering the segment, filter out intersections that do not fall within the segment
            fraction_along_segment = [(xi - p1x) / dx if abs(dx) > abs(dy) else (yi - p1y) / dy for xi, yi in intersections]
            intersections = [pt for pt, frac in zip(intersections, fraction_along_segment) if 0 <= frac <= 1]
        if len(intersections) == 2 and abs(discriminant) <= tangent_tol:  # If line is tangent to circle, return just one point (as both intersections have same location)
            return [intersections[0]]
        else:
            return intersections

Outstand

I haven't tried running this code, but in the first line where it's actually computing intersections, am I nuts or are you multiplying "sign" (defined here as the negative of the actual sign of dy) by the actual sign of dy? In which case, these always reduce to -1, and the value of dy is never included as on the Wolfram page giving the equation. There's a bit too much logical flipping here for me to read this clearly, would be nice if it were expanded for clarity rather than using a clever list comprehension. Sorting the returned points could be left to the reader to do. – Brickwork

Yeah something definitely looks fishy there, what was he thinking?... – Outstand

Well, this code works. I ported it to Javascript, and it works. My problem is, I can't reason geometrically about it as I don't understand what all the variables represent, so for instance I can't understand the relation between sign of dy and the solutions to find the intersection points. I hate having code that works that I can't comprehend... – Brickwork

I think the author was attempting to make it so the intersections would be ordered according to their position along the line segment (so that the first would be closer to pt1), but they made the change in two places instead of one, so it will have zero effect. I suspect you can just remove the ` * (-1 if dy < 0 else 1)` – Outstand

As I grok it, the (-1 if dy < 0 else 1) is the actual sign(x) function. What he names "sign", I renamed plusMinus. His code works because the plusMinus does two passes on the equations, to produce the +/- of the quadratic formula it's derived from. What I don't get is how one figures out the relation between the +/- in the x coord to the +/- in the y coord, since we want 2 solutions, not 2*2 solutions. The author gets this but doesn't explain the reasoning. He's checking the sign of dy in the list comprehension, both to sort the result points but also to ensure this relation. – Brickwork

Very useful, thanks. I had the link with the formula, but as a mere computer programmer, i found it hard to put into code. Your code is very useful. – Inexistent

Answer 3

A low cost spacial partition might be to divide the plane into 9 pieces

Here is a crappy diagram. Imagine, if you will, that the lines are just touching the circle.

  | |
__|_|__
__|O|__
  | |
  | |

8 of the areas we are interested in are surrounding the circle. The square in the centre isn't much use for a cheap test, but you can place a square of r/sqrt(2) inside the circle, so it's corners just touch the circle.

Lets label the areas

A |B| C
__|_|__
D_|O|_E
  | |
F |G| H

And the square of r/sqrt(2) in the centre we'll call J

We will call the set of points in the centre square shown in the diagram that aren't in J, Z

Label each vertex of the polygon with it's letter code.

Now we can quickly see

AA => Outside
AB => Outside
AC => Outside
...
AJ => Intersects
BJ => Intersects
...
JJ => Inside

This can turned into a lookup table

So depending on your dataset, you may have saved yourself a load of work. Anything with an endpoint in Z will need to be tested however.

Sewellyn

Answer 4

I think that the formula you use to find the coordinates of the two intersections cannot be optimized further. The only improvement (which is numerically important) is to distinguish the two cases: |x_2-x_1| >= |y_2-y_1| and |x_2-x1| < |y_2-y1| so that the quantity k is always between -1 and 1 (in your computation you can get very high numerical errors if |x_2-x_1| is very small). You can swap x-s and y-s to reduce one case to the other.

You could also implement a preliminary check: if both endpoints are internal to the circle there are no intersection. By computing the squared distance from the points to the center of the circle this becomes a simple formula which does not use the square root function. The other check: "whether the line is outside the circle" is already implemented in your code and corresponds to delta < 0. If you have a lot of small segments these two check should give a shortcut answer (no intersection) in most cases.

Fleabitten

Answer 5

def sg_check(a,b):
    if b != 0:
        return 1 if a/b >= 0 and a/b <= 1 else -1
    else:
        return -1
def circle_line_seg_intersect(l1x,l1y,l2x,l2y,cx,cy,cr):
    m=(l2y-l1y)/(l2x-l1x+1e-4)
    a=(1+m**2)
    b=2*(m*l1y-cy*m-cx-m**2*l1x)
    c=cx**2+m**2*l1x**2+l1y**2+2*cy*m*l1x-2*l1y*m*l1x-2*cy*l1y+cy**2-cr**2
    root=b**2-4*a*c
    if root >= 0:
        x1=(-b+math.sqrt(root))/(2*a)
        y1=m*(x1-l1x)+l1y
        x2=(-b-math.sqrt(root))/(2*a)
        y2=m*(x2-l1x)+l1y
        sx1=sg_check(x1-l1x,l2x-l1x)
        sy1=sg_check(y1-l1y,l2y-l1y)
        sx2=sg_check(x2-l1x,l2x-l1x)
        sy2=sg_check(y2-l1y,l2y-l1y)
        if sx1+sy1 >= 0 and sx2+sy2 >= 0:
            return ([x1,y1],[x2,y2])
        elif sx1+sy1 >= 0:
            return ([([x1,y1])])
        elif sx2+sy2 >= 0:
            return ([([x2,y2])])
    return []

idk if this is efficient but this is how I did it using the quadratic formula

I implemented this from my desmos graph displaying the intersection of a line and circle Desmos

Crankcase

Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. – Zyrian