28

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
Improve this question
6
  • 1
    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? Commented Jun 15, 2015 at 12:22
  • 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. Commented Jun 15, 2015 at 13:02
  • @EmanuelePaolini,Thank you.I've modified the script according to your concern. Commented Jun 15, 2015 at 13:52
  • 1
    stackoverflow.com/questions/1073336/… - may be relevant Commented Jun 15, 2015 at 14:00
  • Are there still any open issues not answered by any of the answers posted? Commented Jun 30, 2015 at 12:34

5 Answers 5

Reset to default
22

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.

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.

Improve this answer
Sign up to request clarification or add additional context in comments.

9 Comments

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.
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.
Sometimes this returns just a single Point, and sometimes MultiplePoints. How can I test the return value to know which is which?
@albrnick 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)
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.
|
13

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
Improve this answer

6 Comments

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.
Yeah something definitely looks fishy there, what was he thinking?...
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...
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)`
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.
|
1

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.

Improve this answer

Comments

0

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.

Improve this answer

Comments

-1 
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

Improve this answer

1 Comment

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.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.