What is most efficient way to find the intersection of a line and a circle in python?
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 LineStringfrom shapely.geometry import Pointp = Point(5,5)c = p.buffer(3).boundaryl = 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.
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.
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
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 => OutsideAB => OutsideAC => Outside...AJ => IntersectsBJ => 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.