How is this at all specific to Python?
This also sounds suspiciously like homework, which you should know
this list is unlikely to give direct answers to, though you might be
able to get a few pointers or some general suggestions.
Cheers,
Chris
--
Back to toiling away on CSE 105 HW#3
http://blog.rebertia.com
I can see two possible solutions, both approximate. In both solutions
you first look if there are a pair of circles that don't intersect, in
this case the intersect area is zero. You also remove all circles
fully contained in another circle.
The first strategy is easy to do, but it probably leads to a lower
precision. Then you can sample randomly many times the rectangular
area that surely contains all circles. You count how many of those
random points are inside all circles. This gives you an approximate
area of their intersection. Increasing the points numbers slowly
increases the result precision.
The second strategy is more complex. You convert your circles into
polygons with a fixed number of vertices (like 50 or 100 or 1000 or
more. This number is constant even if the circles don't have all the
same radius). So you can turn this circle into a simple mesh of
triangles (all vertices are on the circumference). Then you "subtract"
the second polygonalized circle, this can create a hole and split
triangles in pieces, and so on with successive circles. At the end you
can compute the total area of the triangles left. This is doable, but
you need time to do implement this. The advantage is that the
numerical precision of the result is probably higher. If you implement
this second solution you can implement the first one too, and use it
as a test to avoid bugs. Visualizing the triangles with Pygame or
MatPlotLib can be useful to look for bugs.
Bye,
bearophile
There is a method for calculating the area of a polygon by traversing
its boundary. I forget the formula, but you should be able to find it.
So *if* you can find the arcs that bound the area, you can approximate
each by a series of short lines and apply the polygon method.
Terry Jan Reedy
A brute force approach - create a grid of small squares and calculate
which squares are in all circles. I don't know whether it is any better
than monte-carlo: for two circles, delta=0.001 takes about a minute and
delta=0.0001 is still running after 30 minutes :-)
1.22
1.2281
1.22834799999
------------------------------------------------------------------
import math
class Circle:
def __init__(self, x, y, r=1):
self.x = float(x)
self.y = float(y)
self.r = float(r)
def contains(self, a, b):
return math.sqrt((self.x-a)**2 + (self.y-b)**2) <= self.r
class Grid:
def __init__(self, circles):
self.X1 = min(c.x-c.r for c in circles)
self.Y1 = max(c.y+c.r for c in circles)
self.X2 = max(c.x+c.r for c in circles)
self.Y2 = min(c.y-c.r for c in circles)
self.circles = circles
def iter_boxes(self, delta):
X2 = self.X2
Y2 = self.Y2
a = self.X1
while a < X2:
a += delta
b = self.Y1
while b > Y2:
b -= delta
#print a, b
yield a, b
def intersection(self, delta=0.1):
S = 0
s = delta**2 #box area
circles = self.circles
for a, b in self.iter_boxes(delta):
if all(c.contains(a, b) for c in circles):
S += s
return S
c = Circle(0, 1)
assert c.contains(0, 1)
assert c.contains(0, 1.5)
assert c.contains(math.sqrt(2)/2, math.sqrt(2)/2)
assert c.contains(0, 2)
assert not c.contains(0, 2.01)
assert not c.contains(0, 2.1)
assert not c.contains(0, 3)
circles = [
Circle(0, 0),
Circle(0, 1),
]
g = Grid(circles)
print '-'*30
print g.intersection()
print g.intersection(0.01)
print g.intersection(0.001)
print g.intersection(0.0001)
------------------------------------------------------------------
That's just what the monte-carlo method is -- except the full family of
monte-carlo methods can be quite sophisticated, including being more
general (pseudo or quasi random points instead of a regular grid) and
can provide a theoretical basis for calculating the rate of convergence,
and so on.
Gary Herron
I would have said a Monte Carlo method necessarily involves random
sampling in some shape or form, no? (Casinos, Chance etc..) And as I've
previously understood Monte Carlo Integration, the technique involves
rather the ratio "hits/total sample" rather than the sum of small areas.
That's how I remember it, at least - it's been a while though.
Regards
G.F
This is untested, and probably contains typos, at least.
I would do it with successively smaller squares. A square can be
described using center coordinates (x,y) , and length (s) of a side.
Each square can have three states: not included, partially included,
and fully included.
You build a list of squares, starting with one, the bounding square for
the first circle.
For each square in the list, you go through all the circles (each with
center a,b, and radius r), and for each circle decide if the square is
entirely within the circle, or partially. If it's fully included in all
of them, you add its area to the the accumulator, and drop it from your
list. if it's not included at all in one of them, you just drop it
from the list. (Note, be careful about modifying a list you're
iterating through -- make a copy)
After a pass through the list, the accumulator is a lower-bound for the
area, and the accumulator plus the areas of all the squares is an upper
bound. If these are close enough together, you terminate.
Otherwise you take the remaining squares in the list, split each into
four pieces, and end up with a list 4 times as long. And repeat the
loop, checking each square against each circle.
The only remaining question is how you decide for a given circle and a
given square which of the three states it's in. And it turns out you
can be pretty approximate in this, and it'll still converge. So you can
take a distance between (a,b) and (x,y), and compare it to r-s/2 and to
r+(s*0.707106781187). If it's greater than r+s*0.707106781187, no
intersection. If it's less than r-s/2, then the square's entirely inside.
Note that I used 0.707106781187 above, but you'd actually use sqr(0.5)
to whatever precision you needed. And if it's convenient to round it,
round it upwards; again it won't affect the result. I doubt if it
would matter in Python code, but you could even do the whole thing in
integers (avoiding the sqrt for distance calc), if you were careful
about the bounding formulas. That's premature optimization, however.
HTH,
DaveA
How much accuracy do you need? How fast does the algorithm need to
be? How many circles are typically involved? Is the intersection
necessarily non-empty for the configurations of circles that you use?
How much code are you prepared to write? How much mathematics do you
want to use?
Besides the Monte-Carlo and grid-based approaches already mentioned, I
see a couple of other possibilities:
(1) It should be entirely possible to do this analytically. The hard
part would be identifying the intersection points that form the
vertices of the intersection (and especially, doing this robustly in
the face of numerical errors and triple intersections or near-triple
intersections). Once you've got the vertices, it should be
straightforward to compute the area: compute the area of the polygon
given by the convex hull of the vertices, then add on the extra areas
for the segments bounded by the (straight) sides of the polygon and
the corresponding outer arcs.
(2) For a numerical approach that might be slightly better than the 2d
brute-force approaches described by others: make use of the fact that
the intersection is convex. Pick a point P in the intersection (if it
exists: finding such a point is an issue in itself). For each angle
t, 0 <= t <= 2*pi, define f(t) to be the distance from P to the
boundary of the region, in the direction of the angle t. We always
have 0 <= f(t) <= 1, so f(t) can be quickly and accurately computed
with a bisection search. Now find the definite integral of 0.5 *
(f(t))**2, as t goes from 0 to 2*pi, using your favourite quadrature
method (but avoid methods that would be very sensitive to continuity
of derivatives: f(t) will be continuous, but f'(t) will not).
Simpson's rule is probably good enough.
Good luck!
--
Mark
Good point.
This is actually a problem in what's called "constructive geometry".
Usually, this is a 3D problem, for which the term "constructive solid
geometry", or CSG, is used. That's where to look for algorithms.
Yes, there are efficient algorithms and near-exact solutions.
The basic idea is that you find all the points at which circles
intersect, sort those by one coordinate, and advance across that
coordinate, from one value of interest to the next. Between any
two values of interest you're dealing with areas bounded by arcs
and lines, so the areas can be calculated analytically. It's a
lot like a rasterized polygon fill algorithm.
(I used to do stuff like this back when I did physics engines
with efficient 3D collision detection.)
John Nagle