Total circles area
![Task](http://static.miraheze.org/rosettacodewiki/thumb/b/ba/Rcode-button-task-crushed.png/64px-Rcode-button-task-crushed.png)
You are encouraged to solve this task according to the task description, using any language you may know.
Given some partially overlapping circles on the plane, compute and show the total area covered by them, with four or six (or a little more) decimal digits of precision. The area covered by two or more disks needs to be counted only once.
One point of this Task is also to compare and discuss the relative merits of various solution strategies, their performance, precision and simplicity. This means keeping both slower and faster solutions for a language (like C) is welcome.
To allow a better comparison of the different implementations, solve the problem with this standard dataset, each line contains the x and y coordinates of the centers of the disks and their radii (11 disks are fully contained inside other disks):
![](http://static.miraheze.org/rosettacodewiki/thumb/f/f9/Total_circles_area_full.png/200px-Total_circles_area_full.png)
![](http://static.miraheze.org/rosettacodewiki/thumb/6/6a/Total_circles_area_filtered.png/200px-Total_circles_area_filtered.png)
xc yc radius 1.6417233788 1.6121789534 0.0848270516 -1.4944608174 1.2077959613 1.1039549836 0.6110294452 -0.6907087527 0.9089162485 0.3844862411 0.2923344616 0.2375743054 -0.2495892950 -0.3832854473 1.0845181219 1.7813504266 1.6178237031 0.8162655711 -0.1985249206 -0.8343333301 0.0538864941 -1.7011985145 -0.1263820964 0.4776976918 -0.4319462812 1.4104420482 0.7886291537 0.2178372997 -0.9499557344 0.0357871187 -0.6294854565 -1.3078893852 0.7653357688 1.7952608455 0.6281269104 0.2727652452 1.4168575317 1.0683357171 1.1016025378 1.4637371396 0.9463877418 1.1846214562 -0.5263668798 1.7315156631 1.4428514068 -1.2197352481 0.9144146579 1.0727263474 -0.1389358881 0.1092805780 0.7350208828 1.5293954595 0.0030278255 1.2472867347 -0.5258728625 1.3782633069 1.3495508831 -0.1403562064 0.2437382535 1.3804956588 0.8055826339 -0.0482092025 0.3327165165 -0.6311979224 0.7184578971 0.2491045282 1.4685857879 -0.8347049536 1.3670667538 -0.6855727502 1.6465021616 1.0593087096 0.0152957411 0.0638919221 0.9771215985
The result is 21.5650366038563989590842249288781480183977... .
- See also
C
Montecarlo Sampling
This program uses a Montecarlo sampling. For this problem this is less efficient (converges more slowly) than a regular grid sampling, like in the Python entry. <lang c>#include <stdio.h>
- include <stdlib.h>
- include <math.h>
- include <time.h>
- include <stdbool.h>
typedef double Fp; typedef struct { Fp x, y, r; } Circle;
Circle circles[] = {
{ 1.6417233788, 1.6121789534, 0.0848270516}, {-1.4944608174, 1.2077959613, 1.1039549836}, { 0.6110294452, -0.6907087527, 0.9089162485}, { 0.3844862411, 0.2923344616, 0.2375743054}, {-0.2495892950, -0.3832854473, 1.0845181219}, { 1.7813504266, 1.6178237031, 0.8162655711}, {-0.1985249206, -0.8343333301, 0.0538864941}, {-1.7011985145, -0.1263820964, 0.4776976918}, {-0.4319462812, 1.4104420482, 0.7886291537}, { 0.2178372997, -0.9499557344, 0.0357871187}, {-0.6294854565, -1.3078893852, 0.7653357688}, { 1.7952608455, 0.6281269104, 0.2727652452}, { 1.4168575317, 1.0683357171, 1.1016025378}, { 1.4637371396, 0.9463877418, 1.1846214562}, {-0.5263668798, 1.7315156631, 1.4428514068}, {-1.2197352481, 0.9144146579, 1.0727263474}, {-0.1389358881, 0.1092805780, 0.7350208828}, { 1.5293954595, 0.0030278255, 1.2472867347}, {-0.5258728625, 1.3782633069, 1.3495508831}, {-0.1403562064, 0.2437382535, 1.3804956588}, { 0.8055826339, -0.0482092025, 0.3327165165}, {-0.6311979224, 0.7184578971, 0.2491045282}, { 1.4685857879, -0.8347049536, 1.3670667538}, {-0.6855727502, 1.6465021616, 1.0593087096}, { 0.0152957411, 0.0638919221, 0.9771215985}};
const size_t n_circles = sizeof(circles) / sizeof(Circle);
static inline Fp min(const Fp a, const Fp b) { return a <= b ? a : b; }
static inline Fp max(const Fp a, const Fp b) { return a >= b ? a : b; }
static inline Fp sq(const Fp a) { return a * a; }
// Return an uniform random value in [a, b). static inline double uniform(const double a, const double b) {
const double r01 = rand() / (double)RAND_MAX; return a + (b - a) * r01;
}
static inline bool is_inside_circles(const Fp x, const Fp y) {
for (size_t i = 0; i < n_circles; i++) if (sq(x - circles[i].x) + sq(y - circles[i].y) < circles[i].r) return true; return false;
}
int main() {
// Initialize the bounding box (bbox) of the circles. Fp x_min = INFINITY, x_max = -INFINITY; Fp y_min = x_min, y_max = x_max;
// Compute the bounding box of the circles. for (size_t i = 0; i < n_circles; i++) { Circle *c = &circles[i]; x_min = min(x_min, c->x - c->r); x_max = max(x_max, c->x + c->r); y_min = min(y_min, c->y - c->r); y_max = max(y_max, c->y + c->r);
c->r *= c->r; // Square the radii to speed up testing. }
const Fp bbox_area = (x_max - x_min) * (y_max - y_min);
// Montecarlo sampling. srand(time(0)); size_t to_try = 1U << 16; size_t n_tries = 0; size_t n_hits = 0;
while (true) { n_hits += is_inside_circles(uniform(x_min, x_max), uniform(y_min, y_max)); n_tries++;
if (n_tries == to_try) { const Fp area = bbox_area * n_hits / n_tries; const Fp r = (Fp)n_hits / n_tries; const Fp s = area * sqrt(r * (1 - r) / n_tries); printf("%.4f +/- %.4f (%zd samples)\n", area, s, n_tries); if (s * 3 <= 1e-3) // Stop at 3 sigmas. break; to_try *= 2; } }
return 0;
}</lang>
- Output:
21.4498 +/- 0.0370 (65536 samples) 21.5031 +/- 0.0262 (131072 samples) 21.5170 +/- 0.0185 (262144 samples) 21.5442 +/- 0.0131 (524288 samples) 21.5477 +/- 0.0093 (1048576 samples) 21.5531 +/- 0.0065 (2097152 samples) 21.5624 +/- 0.0046 (4194304 samples) 21.5631 +/- 0.0033 (8388608 samples) 21.5602 +/- 0.0023 (16777216 samples) 21.5632 +/- 0.0016 (33554432 samples) 21.5617 +/- 0.0012 (67108864 samples) 21.5628 +/- 0.0008 (134217728 samples) 21.5639 +/- 0.0006 (268435456 samples) 21.5637 +/- 0.0004 (536870912 samples) 21.5637 +/- 0.0003 (1073741824 samples)
Scanline Method
This version performs about 5 million scanlines in about a second, result should be accurate to maybe 10 decimal points. <lang c>#include <stdio.h>
- include <string.h>
- include <stdlib.h>
- include <math.h>
typedef double flt; typedef struct { flt x, y, r, r2; flt y0, y1; // extent of circle y+r and y-r flt x0, x1; // where scanline intersects circle } circle_t;
- define SZ sizeof(circle_t)
circle_t circles[] = { { 1.6417233788, 1.6121789534, 0.0848270516},
- error data snipped for space; copy from previous C example
};
flt max(flt x, flt y) { return x < y ? y : x; } flt min(flt x, flt y) { return x > y ? y : x; } flt sq(flt x) { return x * x; } flt cdist(circle_t *c1, circle_t *c2) { return sqrt(sq(c1->x - c2->x) + sq(c1->y - c2->y)); }
inline void swap_c(circle_t *c) { circle_t tmp = c[0]; c[0] = c[1], c[1] = tmp; }
flt area(circle_t *circs, int n_circ, flt ymin, flt ymax, flt step) { int i, n = n_circ;
circle_t *c = malloc(SZ * n); memcpy(c, circs, SZ * n);
while (n--) for (i = 0; i < n; i++) if (c[i].y1 < c[i+1].y1) swap_c(c + i);
flt total = 0;
int row = 1 + ceil((ymax - ymin) / step); while (row--) { flt y = ymin + step * row; for (n = 0; n < n_circ; n++) if (y >= c[n].y1) // rest of circles below scanline, ignore break; else if (y > c[n].y0) { flt dx = sqrt(c[n].r2 - sq(y - c[n].y)); c[n].x0 = c[n].x - dx; c[n].x1 = c[n].x + dx;
// keep circles sorted by left intersection for (i = n; i-- && c[i].x0 > c[i+1].x0; swap_c(c + i));
} else {// remove a circle when scanline has passed it memmove(c + n, c + n + 1, SZ * (--n_circ - n)); n--; }
if (!n) continue;
flt right = c->x1; total += c->x1 - c->x0;
for (i = 1; i < n; i++) { if (c[i].x1 <= right) continue; total += c[i].x1 - max(c[i].x0, right); right = c[i].x1; } }
free(c); return total * step; }
int main(void) { int n_circ = sizeof(circles) / SZ; flt ymin = INFINITY, ymax = -INFINITY;
circle_t *c1, *c2; for (c1 = circles + n_circ; c1-- > circles; ) { for (c2 = circles + n_circ; c2-- > circles; ) // throw out circles inside another circle if (c1 != c2 && cdist(c1, c2) + c1->r <= c2->r) { *c1 = circles[--n_circ]; break; } ymin = min(ymin, c1->y0 = c1->y - c1->r); ymax = max(ymax, c1->y1 = c1->y + c1->r); c1->r2 = sq(c1->r); }
flt s = 1. / (1 << 20); flt y0 = floor(ymin / s) * s; flt a = area(circles, n_circ, y0, ymax, s); int nlines = (ymax - y0) / s;
// roughly, it cease to make sense if sqrt(nlines) * 1e-14 >> s * s printf("area = %.10f\tat %d scanlines\n", a, nlines);
return 0; }</lang>
- Output:
area = 21.5650366037 at 5637290 scanlines
D
This version converges much faster than both the ordered grid and Montecarlo sampling solutions.
Scanline Method
<lang d>import std.stdio, std.math, std.algorithm, std.typecons, std.range;
alias Fp = real; struct Circle { Fp x, y, r; }
void removeInternalDisks(ref Circle[] circles) pure nothrow @safe {
static bool isFullyInternal(in Circle c1, in Circle c2) pure nothrow @safe @nogc { if (c1.r > c2.r) // Quick exit. return false; return (c1.x - c2.x) ^^ 2 + (c1.y - c2.y) ^^ 2 < (c2.r - c1.r) ^^ 2; }
// Heuristics for performance: large radii first. circles.sort!q{ a.r > b.r };
// Remove circles inside another circle. for (auto i = circles.length; i-- > 0; ) for (auto j = circles.length; j-- > 0; ) if (i != j && isFullyInternal(circles[i], circles[j])) { circles[i] = circles[$ - 1]; circles.length--; break; }
}
void main() {
Circle[] circles = [ { 1.6417233788, 1.6121789534, 0.0848270516}, {-1.4944608174, 1.2077959613, 1.1039549836}, { 0.6110294452, -0.6907087527, 0.9089162485}, { 0.3844862411, 0.2923344616, 0.2375743054}, {-0.2495892950, -0.3832854473, 1.0845181219}, { 1.7813504266, 1.6178237031, 0.8162655711}, {-0.1985249206, -0.8343333301, 0.0538864941}, {-1.7011985145, -0.1263820964, 0.4776976918}, {-0.4319462812, 1.4104420482, 0.7886291537}, { 0.2178372997, -0.9499557344, 0.0357871187}, {-0.6294854565, -1.3078893852, 0.7653357688}, { 1.7952608455, 0.6281269104, 0.2727652452}, { 1.4168575317, 1.0683357171, 1.1016025378}, { 1.4637371396, 0.9463877418, 1.1846214562}, {-0.5263668798, 1.7315156631, 1.4428514068}, {-1.2197352481, 0.9144146579, 1.0727263474}, {-0.1389358881, 0.1092805780, 0.7350208828}, { 1.5293954595, 0.0030278255, 1.2472867347}, {-0.5258728625, 1.3782633069, 1.3495508831}, {-0.1403562064, 0.2437382535, 1.3804956588}, { 0.8055826339, -0.0482092025, 0.3327165165}, {-0.6311979224, 0.7184578971, 0.2491045282}, { 1.4685857879, -0.8347049536, 1.3670667538}, {-0.6855727502, 1.6465021616, 1.0593087096}, { 0.0152957411, 0.0638919221, 0.9771215985}];
writeln("Input Circles: ", circles.length); removeInternalDisks(circles); writeln("Circles left: ", circles.length);
immutable Fp xMin = reduce!((acc, c) => min(acc, c.x - c.r)) (Fp.max, circles[]); immutable Fp xMax = reduce!((acc, c) => max(acc, c.x + c.r)) (Fp(0), circles[]);
alias YRange = Tuple!(Fp,"y0", Fp,"y1"); auto yRanges = new YRange[circles.length];
Fp computeTotalArea(in Fp nSlicesX) nothrow @safe { Fp total = 0;
// Adapted from an idea by Cosmologicon. foreach (immutable p; cast(int)(xMin * nSlicesX) .. cast(int)(xMax * nSlicesX) + 1) { immutable Fp x = p / nSlicesX; size_t nPairs = 0;
// Look for the circles intersecting the current // vertical secant: foreach (const ref c; circles) { immutable Fp d = c.r ^^ 2 - (c.x - x) ^^ 2; immutable Fp sd = d.sqrt; if (d > 0) // And keep only the intersection chords. yRanges[nPairs++] = YRange(c.y - sd, c.y + sd); }
// Merge the ranges, counting the overlaps only once. yRanges[0 .. nPairs].sort(); Fp y = -Fp.max; foreach (immutable r; yRanges[0 .. nPairs]) if (y < r.y1) { total += r.y1 - max(y, r.y0); y = r.y1; } }
return total / nSlicesX; }
// Iterate to reach some precision. enum Fp epsilon = 1e-9; Fp nSlicesX = 1_000; Fp oldArea = -1; while (true) { immutable Fp newArea = computeTotalArea(nSlicesX); if (abs(oldArea - newArea) < epsilon) { writeln("N. vertical slices: ", nSlicesX); writefln("Approximate area: %.17f", newArea); return; } oldArea = newArea; nSlicesX *= 2; }
}</lang>
- Output:
Input Circles: 25 Circles left: 14 N. vertical slices: 256000 Approximate area: 21.56503660593628004
Run-time is about 4.13 seconds with ldc2 compiler.
Analytical Solution
This version is not fully idiomatic D, it retains some of the style of the Haskell version. <lang d>import std.stdio, std.typecons, std.math, std.algorithm, std.range;
struct Vec { double x, y; }
alias VF = double function(in Vec, in Vec) pure nothrow @safe; enum VF vCross = (v1, v2) => v1.x * v2.y - v1.y * v2.x; enum VF vDot = (v1, v2) => v1.x * v2.x + v1.y * v2.y;
alias VV = Vec function(in Vec, in Vec) pure nothrow @safe; enum VV vAdd = (v1, v2) => Vec(v1.x + v2.x, v1.y + v2.y); enum VV vSub = (v1, v2) => Vec(v1.x - v2.x, v1.y - v2.y);
enum vLen = (in Vec v) pure nothrow @safe => vDot(v, v).sqrt; enum VF vDist = (a, b) => vSub(a, b).vLen; enum vScale = (in double s, in Vec v) pure nothrow @safe =>
Vec(v.x * s, v.y * s);
enum vNorm = (in Vec v) pure nothrow @safe =>
Vec(v.x / v.vLen, v.y / v.vLen);
alias A = Typedef!double;
enum vAngle = (in Vec v) pure nothrow @safe => atan2(v.y, v.x).A;
A aNorm(in A a) pure nothrow @safe {
if (a > PI) return A(a - PI * 2.0); if (a < -PI) return A(a + PI * 2.0); return a;
}
struct Circle { double x, y, r; }
A[] circleCross(in Circle c0, in Circle c1) pure nothrow {
immutable d = vDist(Vec(c0.x, c0.y), Vec(c1.x, c1.y)); if (d >= c0.r + c1.r || d <= abs(c0.r - c1.r)) return [];
immutable s = (c0.r + c1.r + d) / 2.0; immutable a = sqrt(s * (s - d) * (s - c0.r) * (s - c1.r)); immutable h = 2.0 * a / d; immutable dr = Vec(c1.x - c0.x, c1.y - c0.y); immutable dx = vScale(sqrt(c0.r ^^ 2 - h ^^ 2), dr.vNorm); immutable ang = (c0.r ^^ 2 + d ^^ 2 > c1.r ^^ 2) ? dr.vAngle : A(PI + dr.vAngle); immutable da = asin(h / c0.r).A; return [A(ang - da), A(ang + da)].map!aNorm.array;
}
// Angles of the start and end points of the circle arc. alias Angle2 = Tuple!(A,"a0", A,"a1");
alias Arc = Tuple!(Circle,"c", Angle2,"aa");
enum arcPoint = (in Circle c, in A a) pure nothrow @safe =>
vAdd(Vec(c.x, c.y), Vec(c.r * cos(cast(double)a), c.r * sin(cast(double)a)));
alias ArcF = Vec function(in Arc) pure nothrow @safe; enum ArcF arcStart = ar => arcPoint(ar.c, ar.aa.a0); enum ArcF arcMid = ar => arcPoint(ar.c, A((ar.aa.a0 + ar.aa.a1) / 2)); enum ArcF arcEnd = ar => arcPoint(ar.c, ar.aa.a1); enum ArcF arcCenter = ar => Vec(ar.c.x, ar.c.y);
enum arcArea = (in Arc ar) pure nothrow @safe =>
ar.c.r ^^ 2 * (ar.aa.a1 - ar.aa.a0) / 2.0;
Arc[] splitCircles(in Circle[] cs) /*pure nothrow*/ {
static enum cSplit = (in Circle c, in A[] angs) pure nothrow => c.repeat.zip(angs.zip(angs.dropOne).map!Angle2).map!Arc;
// If an arc that was part of one circle is inside *another* circle, // it will not be part of the zero-winding path, so reject it. static bool inCircle(VC)(in VC vc, in Circle c) pure nothrow { return vc[1] != c && vDist(Vec(vc[0].x, vc[0].y), Vec(c.x, c.y)) < c.r; }
enum inAnyCircle = (in Arc arc) nothrow @safe => cs.map!(c => inCircle(tuple(arc.arcMid, arc.c), c)).any;
auto f(in Circle c) /*pure nothrow*/ { auto angs = cs.map!(c1 => circleCross(c, c1)).join; return tuple(c, ([A(-PI), A(PI)] ~ angs) .sort().release); }
return cs.map!f.map!(ca => cSplit(ca[])).join .filter!(ar => !inAnyCircle(ar)).array;
}
/** Given a list of arcs, build sets of closed paths from them. If
one arc's end point is no more than 1e-4 from another's start point,
they are considered connected. Since these start/end points resulted
from intersecting circles earlier, they *should* be exactly the same,
but floating point precision may cause small differences, hence the
1e-4 error margin. When there are genuinely different intersections
closer than this margin, the method will backfire, badly. */
const(Arc[])[] makePaths(in Arc[] arcs) pure nothrow @safe {
static const(Arc[])[] joinArcs(in Arc[] a, in Arc[] xxs) pure nothrow { static enum eps = 1e-4; if (xxs.empty) return [a]; immutable x = xxs[0]; const xs = xxs.dropOne; if (a.empty) return joinArcs([x], xs); if (vDist(a[0].arcStart, a.back.arcEnd) < eps) return [a] ~ joinArcs([], xxs); if (vDist(a.back.arcEnd, x.arcStart) < eps) return joinArcs(a ~ [x], xs); return joinArcs(a, xs ~ [x]); } return joinArcs([], arcs);
}
// Slice N-polygon into N-2 triangles. double polylineArea(in Vec[] vvs) /*pure*/ nothrow {
static enum triArea= (in Vec a, in Vec b, in Vec c) pure nothrow => vCross(vSub(b, a), vSub(c, b)) / 2.0; const vs = vvs.dropOne; // sum(double[]) isn't pure yet. return zip(vs, vs.dropOne).map!(vv => triArea(vvs[0], vv[])).sum;
}
double pathArea(in Arc[] arcs) /*pure*/ nothrow {
static f(/*in*/ Tuple!(double, Vec[]) ae, in Arc arc) pure nothrow { return tuple(ae[0] + arc.arcArea, ae[1] ~ [arc.arcCenter, arc.arcEnd]); } const ae = reduce!f(tuple(0.0, (Vec[]).init), arcs); return ae[0] + ae[1].polylineArea;
}
enum circlesArea = (in Circle[] cs) /*pure nothrow*/ =>
cs.splitCircles.makePaths.map!pathArea.sum;
void main() {
immutable circles = [ Circle( 1.6417233788, 1.6121789534, 0.0848270516), Circle(-1.4944608174, 1.2077959613, 1.1039549836), Circle( 0.6110294452, -0.6907087527, 0.9089162485), Circle( 0.3844862411, 0.2923344616, 0.2375743054), Circle(-0.2495892950, -0.3832854473, 1.0845181219), Circle( 1.7813504266, 1.6178237031, 0.8162655711), Circle(-0.1985249206, -0.8343333301, 0.0538864941), Circle(-1.7011985145, -0.1263820964, 0.4776976918), Circle(-0.4319462812, 1.4104420482, 0.7886291537), Circle( 0.2178372997, -0.9499557344, 0.0357871187), Circle(-0.6294854565, -1.3078893852, 0.7653357688), Circle( 1.7952608455, 0.6281269104, 0.2727652452), Circle( 1.4168575317, 1.0683357171, 1.1016025378), Circle( 1.4637371396, 0.9463877418, 1.1846214562), Circle(-0.5263668798, 1.7315156631, 1.4428514068), Circle(-1.2197352481, 0.9144146579, 1.0727263474), Circle(-0.1389358881, 0.1092805780, 0.7350208828), Circle( 1.5293954595, 0.0030278255, 1.2472867347), Circle(-0.5258728625, 1.3782633069, 1.3495508831), Circle(-0.1403562064, 0.2437382535, 1.3804956588), Circle( 0.8055826339, -0.0482092025, 0.3327165165), Circle(-0.6311979224, 0.7184578971, 0.2491045282), Circle( 1.4685857879, -0.8347049536, 1.3670667538), Circle(-0.6855727502, 1.6465021616, 1.0593087096), Circle( 0.0152957411, 0.0638919221, 0.9771215985)];
writefln("Area: %1.13f", circles.circlesArea);
}</lang>
- Output:
Area: 21.5650366038564
The run-time is minimal (0.03 seconds or less).
Haskell
Grid Sampling Version
<lang haskell>data Circle = Circle { cx :: Double, cy :: Double, cr :: Double }
isInside :: Double -> Double -> Circle -> Bool isInside x y c = (x - cx c) ^ 2 + (y - cy c) ^ 2 <= (cr c ^ 2)
isInsideAny :: Double -> Double -> [Circle] -> Bool isInsideAny x y = any (isInside x y)
approximatedArea :: [Circle] -> Int -> Double approximatedArea cs box_side = (fromIntegral count) * dx * dy
where -- compute the bounding box of the circles x_min = minimum [cx c - cr c | c <- circles] x_max = maximum [cx c + cr c | c <- circles] y_min = minimum [cy c - cr c | c <- circles] y_max = maximum [cy c + cr c | c <- circles] dx = (x_max - x_min) / (fromIntegral box_side) dy = (y_max - y_min) / (fromIntegral box_side) count = length [0 | r <- [0 .. box_side - 1], c <- [0 .. box_side - 1], isInsideAny (posx c) (posy r) circles] posy r = y_min + (fromIntegral r) * dy posx c = x_min + (fromIntegral c) * dx
circles :: [Circle] circles = [Circle ( 1.6417233788) ( 1.6121789534) 0.0848270516,
Circle (-1.4944608174) ( 1.2077959613) 1.1039549836, Circle ( 0.6110294452) (-0.6907087527) 0.9089162485, Circle ( 0.3844862411) ( 0.2923344616) 0.2375743054, Circle (-0.2495892950) (-0.3832854473) 1.0845181219, Circle ( 1.7813504266) ( 1.6178237031) 0.8162655711, Circle (-0.1985249206) (-0.8343333301) 0.0538864941, Circle (-1.7011985145) (-0.1263820964) 0.4776976918, Circle (-0.4319462812) ( 1.4104420482) 0.7886291537, Circle ( 0.2178372997) (-0.9499557344) 0.0357871187, Circle (-0.6294854565) (-1.3078893852) 0.7653357688, Circle ( 1.7952608455) ( 0.6281269104) 0.2727652452, Circle ( 1.4168575317) ( 1.0683357171) 1.1016025378, Circle ( 1.4637371396) ( 0.9463877418) 1.1846214562, Circle (-0.5263668798) ( 1.7315156631) 1.4428514068, Circle (-1.2197352481) ( 0.9144146579) 1.0727263474, Circle (-0.1389358881) ( 0.1092805780) 0.7350208828, Circle ( 1.5293954595) ( 0.0030278255) 1.2472867347, Circle (-0.5258728625) ( 1.3782633069) 1.3495508831, Circle (-0.1403562064) ( 0.2437382535) 1.3804956588, Circle ( 0.8055826339) (-0.0482092025) 0.3327165165, Circle (-0.6311979224) ( 0.7184578971) 0.2491045282, Circle ( 1.4685857879) (-0.8347049536) 1.3670667538, Circle (-0.6855727502) ( 1.6465021616) 1.0593087096, Circle ( 0.0152957411) ( 0.0638919221) 0.9771215985]
main = putStrLn $ "Approximated area: " ++
(show $ approximatedArea circles 5000)</lang>
- Output:
Approximated area: 21.564955642878786
Analytical Solution
Breaking down circles to non-intersecting arcs and assemble zero winding paths, then calculate their areas. Pro: precision doesn't depend on a step size, so no need to wait longer for a more precise result; Con: probably not numerically stable in marginal situations, which can be catastrophic. <lang haskell>{-# LANGUAGE GeneralizedNewtypeDeriving #-}
import Data.List (sort)
data Vec = Vec Double Double
vCross, vDot :: Vec -> Vec -> Double vCross (Vec a b) (Vec c d) = a*d - b*c vDot (Vec a b) (Vec c d) = a*c + b*d
vAdd, vSub :: Vec -> Vec -> Vec vAdd (Vec a b) (Vec c d) = Vec (a + c) (b + d) vSub (Vec a b) (Vec c d) = Vec (a - c) (b - d)
vLen :: Vec -> Double vLen x = sqrt $ vDot x x
vDist :: Vec -> Vec -> Double vDist a b = vLen (a `vSub` b)
vScale :: Double -> Vec -> Vec vScale s (Vec x y) = Vec (x * s) (y * s)
vNorm :: Vec -> Vec vNorm v@(Vec x y) = Vec (x / l) (y / l) where l = vLen v
newtype Angle = A Double
deriving (Eq, Ord, Num, Fractional)
aPi = A pi
vAngle :: Vec -> Angle vAngle (Vec x y) = A (atan2 y x)
aNorm :: Angle -> Angle aNorm a | a > aPi = a - aPi * 2
| a < -aPi = a + aPi * 2 | otherwise = a
data Circle = Circle Double Double Double
deriving (Eq)
circleCross :: Circle -> Circle -> [Angle] circleCross (Circle x0 y0 r0) (Circle x1 y1 r1)
| d >= r0 + r1 || d <= abs(r0 - r1) = [] | otherwise = map aNorm [ang - da, ang + da] where d = vDist (Vec x0 y0) (Vec x1 y1) s = (r0 + r1 + d) / 2 a = sqrt $ s * (s - d) * (s - r0) * (s - r1) h = 2 * a / d dr = Vec (x1 - x0) (y1 - y0) dx = vScale (sqrt $ r0 ^ 2 - h ^ 2) $ vNorm dr ang = if r0 ^ 2 + d ^ 2 > r1 ^ 2 then vAngle dr else aPi + vAngle dr da = A (asin (h / r0))
-- Angles of the start and end points of the circle arc.
data Angle2 = Angle2 Angle Angle
data Arc = Arc Circle Angle2
arcPoint :: Circle -> Angle -> Vec arcPoint (Circle x y r) (A a) =
vAdd (Vec x y) (Vec (r * cos a) (r * sin a))
arcStart, arcMid, arcEnd, arcCenter :: Arc -> Vec arcStart (Arc c (Angle2 a0 a1)) = arcPoint c a0 arcMid (Arc c (Angle2 a0 a1)) = arcPoint c ((a0 + a1) / 2) arcEnd (Arc c (Angle2 a0 a1)) = arcPoint c a1 arcCenter (Arc (Circle x y r) _) = Vec x y
arcArea :: Arc -> Double arcArea (Arc (Circle _ _ r) (Angle2 a0 a1)) = r ^ 2 * aDiff / 2
where (A aDiff) = a1 - a0
splitCircles :: [Circle] -> [Arc]
splitCircles cs = filter (not . inAnyCircle) arcs where
cSplit :: (Circle, [Angle]) -> [Arc] cSplit (c, angs) = zipWith Arc (repeat c) $ zipWith Angle2 angs $ tail angs
-- If an arc that was part of one circle is inside *another* circle, -- it will not be part of the zero-winding path, so reject it. inCircle :: (Vec, Circle) -> Circle -> Bool inCircle (Vec x0 y0, c1) c2@(Circle x y r) = c1 /= c2 && vDist (Vec x0 y0) (Vec x y) < r
f :: Circle -> (Circle, [Angle]) f c = (c, sort $ [-aPi, aPi] ++ (concatMap (circleCross c) cs)) cAngs = map f cs arcs = concatMap cSplit cAngs
inAnyCircle :: Arc -> Bool inAnyCircle arc@(Arc c _) = any (inCircle (arcMid arc, c)) cs
{-
Given a list of arcs, build sets of closed paths from them.
If one arc's end point is no more than 1e-4 from another's
start point, they are considered connected. Since these
start/end points resulted from intersecting circles earlier,
they *should* be exactly the same, but floating point
precision may cause small differences, hence the 1e-4 error
margin. When there are genuinely different intersections
closer than this margin, the method will backfire, badly.
-}
makePaths :: [Arc] -> Arc
makePaths arcs = joinArcs [] arcs where
joinArcs :: [Arc] -> [Arc] -> Arc joinArcs a [] = [a] joinArcs [] (x:xs) = joinArcs [x] xs joinArcs a (x:xs) | vDist (arcStart (head a)) (arcEnd (last a)) < 1e-4 = a : joinArcs [] (x:xs) | vDist (arcEnd (last a)) (arcStart x) < 1e-4 = joinArcs (a ++ [x]) xs | otherwise = joinArcs a (xs ++ [x])
pathArea :: [Arc] -> Double
pathArea arcs = a + polylineArea e where
(a, e) = foldl f (0, []) arcs f (a, e) arc = (a + arcArea arc, e ++ [arcCenter arc, arcEnd arc])
-- Slice N-polygon into N-2 triangles.
polylineArea :: [Vec] -> Double
polylineArea (v:vs) = sum $ zipWith (triArea v) vs (tail vs)
where triArea a b c = ((b `vSub` a) `vCross` (c `vSub` b)) / 2
circlesArea :: [Circle] -> Double
circlesArea = sum . map pathArea . makePaths . splitCircles
circles :: [Circle]
circles = [Circle ( 1.6417233788) ( 1.6121789534) 0.0848270516,
Circle (-1.4944608174) ( 1.2077959613) 1.1039549836, Circle ( 0.6110294452) (-0.6907087527) 0.9089162485, Circle ( 0.3844862411) ( 0.2923344616) 0.2375743054, Circle (-0.2495892950) (-0.3832854473) 1.0845181219, Circle ( 1.7813504266) ( 1.6178237031) 0.8162655711, Circle (-0.1985249206) (-0.8343333301) 0.0538864941, Circle (-1.7011985145) (-0.1263820964) 0.4776976918, Circle (-0.4319462812) ( 1.4104420482) 0.7886291537, Circle ( 0.2178372997) (-0.9499557344) 0.0357871187, Circle (-0.6294854565) (-1.3078893852) 0.7653357688, Circle ( 1.7952608455) ( 0.6281269104) 0.2727652452, Circle ( 1.4168575317) ( 1.0683357171) 1.1016025378, Circle ( 1.4637371396) ( 0.9463877418) 1.1846214562, Circle (-0.5263668798) ( 1.7315156631) 1.4428514068, Circle (-1.2197352481) ( 0.9144146579) 1.0727263474, Circle (-0.1389358881) ( 0.1092805780) 0.7350208828, Circle ( 1.5293954595) ( 0.0030278255) 1.2472867347, Circle (-0.5258728625) ( 1.3782633069) 1.3495508831, Circle (-0.1403562064) ( 0.2437382535) 1.3804956588, Circle ( 0.8055826339) (-0.0482092025) 0.3327165165, Circle (-0.6311979224) ( 0.7184578971) 0.2491045282, Circle ( 1.4685857879) (-0.8347049536) 1.3670667538, Circle (-0.6855727502) ( 1.6465021616) 1.0593087096, Circle ( 0.0152957411) ( 0.0638919221) 0.9771215985]
main = print $ circlesArea circles</lang>
- Output:
21.5650366038564
This is how this solution works:
- Given a list of circles, give all of them the same winding. Say, imagine every circle turns clockwise.
- Find all the intersection points among all circles. Between each pair, there may be 0, 1 or 2 intersections. Ignore 0 and 1, we need only the 2 case.
- For each circle, sort its intersection points in clockwise order (keep in mind it's a cyclic list), and split it into arc segments between neighboring point pairs. Circles that don't cross other circles are treated as a single 360 degree arc.
- Imagine all circles are on a white piece of paper, and have their interiors inked black. We are only interested in the arcs separating black and white areas, so we get rid of the arcs that aren't so. These are the arcs that lie entirely within another circle. Because we've taken care of all intersections eariler, we only need to check if any point on an arc segment is in any other circle; if so, remove this arc. (The haskell code uses the middle point of an arc for this, but anything other than the two end points would do).
- The remaining arcs form one or more closed paths. Each path's area can be calculated, and the sum of them all is the area needed. Each path is done like the way mentioned on that stackexchange page cited somewhere. This works for holes, too. Suppose the circles form an area with a hole in it; you'd end up with two paths, where the outer one winds clockwise, and the inner one ccw. Use the same method to calculate the areas for both, and the outer one would have a positive area, the inner one negative. Just add them up.
There are a few concerns about this algorithm: firstly, it's fast only if there are few circles. Its complexity is maybe O(N^3) with N = number of circles, while normal scanline method is probably O(N * n) or less, with n = number of scanlines. Secondly, step 4 needs to be accurate; a small precision error there may cause an arc to remain or be removed by mistake, with disastrous consequences. Also, it's difficult to estimate the error in the final result. The scanline or Monte Carlo methods have errors mostly due to statistics, while this method's error is due to floating point precision loss, which is a very different can of worms.
J
Uniform Grid
We're missing an error estimate. Because it happened to be fairly accurate isn't proof that it's good. Runtime is 16 seconds on a Lenovo T500 with plenty of memory and connected to the power grid. <lang J>NB. check points on a regular grid within the bounding box
N=: 400 NB. grids in each dimension. Controls accuracy.
'X Y R'=: |: XYR=: (_&".;._2~ LF&=)0 :0
1.6417233788 1.6121789534 0.0848270516
-1.4944608174 1.2077959613 1.1039549836
0.6110294452 -0.6907087527 0.9089162485 0.3844862411 0.2923344616 0.2375743054
-0.2495892950 -0.3832854473 1.0845181219
1.7813504266 1.6178237031 0.8162655711
-0.1985249206 -0.8343333301 0.0538864941 -1.7011985145 -0.1263820964 0.4776976918 -0.4319462812 1.4104420482 0.7886291537
0.2178372997 -0.9499557344 0.0357871187
-0.6294854565 -1.3078893852 0.7653357688
1.7952608455 0.6281269104 0.2727652452 1.4168575317 1.0683357171 1.1016025378 1.4637371396 0.9463877418 1.1846214562
-0.5263668798 1.7315156631 1.4428514068 -1.2197352481 0.9144146579 1.0727263474 -0.1389358881 0.1092805780 0.7350208828
1.5293954595 0.0030278255 1.2472867347
-0.5258728625 1.3782633069 1.3495508831 -0.1403562064 0.2437382535 1.3804956588
0.8055826339 -0.0482092025 0.3327165165
-0.6311979224 0.7184578971 0.2491045282
1.4685857879 -0.8347049536 1.3670667538
-0.6855727502 1.6465021616 1.0593087096
0.0152957411 0.0638919221 0.9771215985
)
bbox=: (<./@:- , >./@:+)&R BBOXX=: bbox X BBOXY=: bbox Y
grid=: 3 : 0 'MN MX N'=. y D=. MX-MN EDGE=. D%N (MN(+ -:)EDGE)+(D-EDGE)*(i. % <:)N )
assert 2.2 2.6 3 3.4 3.8 -: grid 2 4 5
GRIDDED_SAMPLES=: BBOXX {@:;&(grid@:(,&N)) BBOXY
Note '4 4{.GRIDDED_SAMPLES' NB. example ┌─────────────────┬─────────────────┬─────────────────┬─────────────────┐ │_2.59706 _2.20043│_2.59706 _2.19774│_2.59706 _2.19505│_2.59706 _2.19236│ ├─────────────────┼─────────────────┼─────────────────┼─────────────────┤ │_2.59434 _2.20043│_2.59434 _2.19774│_2.59434 _2.19505│_2.59434 _2.19236│ ├─────────────────┼─────────────────┼─────────────────┼─────────────────┤ │_2.59162 _2.20043│_2.59162 _2.19774│_2.59162 _2.19505│_2.59162 _2.19236│ ├─────────────────┼─────────────────┼─────────────────┼─────────────────┤ │_2.58891 _2.20043│_2.58891 _2.19774│_2.58891 _2.19505│_2.58891 _2.19236│ └─────────────────┴─────────────────┴─────────────────┴─────────────────┘ ) XY=: >,GRIDDED_SAMPLES NB. convert to an usual array of floats.
mp=: $:~ :(+/ .*) NB. matrix product assert (*: 5 13) -: (mp"1) 3 4,:5 12
in=: *:@:{:@:] >: [: mp (- }:) NB. logical function assert 0 0 in 1 0 2 NB. X Y in X Y R assert 0 0 (-.@:in) 44 2 3
CONTAINED=: XY in"1/XYR NB. logical table of circles containing each grid FRACTION=: CONTAINED (+/@:(+./"1)@:[ % *:@:]) N AREA=: BBOXX*&(-/)BBOXY NB. area of the bounding box. FRACTION*AREA
NB. result is 21.5645</lang>
Perl 6
This subdivides the outer rectangle repeatedly into subrectangles, and classifies them into wet, dry, or unknown. The knowns are summed to provide an inner bound and an outer bound, while the unknowns are further subdivided. The estimate is the average of the outer bound and the inner bound. Not the simplest algorithm, but converges fairly rapidly because it can treat large areas sparsely, saving the fine subdivisions for the circle boundaries. The number of unknown rectangles roughly doubles each pass, but the area of those unknowns is about half. <lang perl6>class Point {
has Real $.x; has Real $.y; has Int $!cbits; # bitmap of circle membership
method cbits { $!cbits //= set_cbits(self) } method gist { $!x ~ "\t" ~ $!y }
}
multi infix:<to>(Point $p1, Point $p2) {
sqrt ($p1.x - $p2.x) ** 2 + ($p1.y - $p2.y) ** 2;
}
multi infix:<mid>(Point $p1, Point $p2) {
Point.new(x => ($p1.x + $p2.x) / 2, y => ($p1.y + $p2.y) / 2);
}
class Circle {
has Point $.center; has Real $.radius;
has Point $.north = Point.new(x => $!center.x, y => $!center.y + $!radius); has Point $.west = Point.new(x => $!center.x - $!radius, y => $!center.y); has Point $.south = Point.new(x => $!center.x, y => $!center.y - $!radius); has Point $.east = Point.new(x => $!center.x + $!radius, y => $!center.y);
multi method contains(Circle $c) { $!center to $c.center <= $!radius - $c.radius } multi method contains(Point $p) { $!center to $p <= $!radius } method gist { $!center.gist ~ "\t" ~ $.radius }
}
class Rect {
has Point $.nw; has Point $.ne; has Point $.sw; has Point $.se;
method diag { $!ne to $!se } method area { ($!ne.x - $!nw.x) * ($!nw.y - $!sw.y) } method contains(Point $p) {
$!nw.x < $p.x < $!ne.x and $!sw.y < $p.y < $!nw.y;
}
}
my @rawcircles = sort -*.radius,
map -> $x, $y, $radius { Circle.new(:center(Point.new(:$x, :$y)), :$radius) }, <
1.6417233788 1.6121789534 0.0848270516 -1.4944608174 1.2077959613 1.1039549836 0.6110294452 -0.6907087527 0.9089162485 0.3844862411 0.2923344616 0.2375743054 -0.2495892950 -0.3832854473 1.0845181219 1.7813504266 1.6178237031 0.8162655711 -0.1985249206 -0.8343333301 0.0538864941 -1.7011985145 -0.1263820964 0.4776976918 -0.4319462812 1.4104420482 0.7886291537 0.2178372997 -0.9499557344 0.0357871187 -0.6294854565 -1.3078893852 0.7653357688 1.7952608455 0.6281269104 0.2727652452 1.4168575317 1.0683357171 1.1016025378 1.4637371396 0.9463877418 1.1846214562 -0.5263668798 1.7315156631 1.4428514068 -1.2197352481 0.9144146579 1.0727263474 -0.1389358881 0.1092805780 0.7350208828 1.5293954595 0.0030278255 1.2472867347 -0.5258728625 1.3782633069 1.3495508831 -0.1403562064 0.2437382535 1.3804956588 0.8055826339 -0.0482092025 0.3327165165 -0.6311979224 0.7184578971 0.2491045282 1.4685857879 -0.8347049536 1.3670667538 -0.6855727502 1.6465021616 1.0593087096 0.0152957411 0.0638919221 0.9771215985
>».Num;
- remove redundant circles
my @circles; while @rawcircles {
my $c = @rawcircles.shift; next if @circles.any.contains($c); push @circles, $c;
}
sub set_cbits(Point $p) {
my $cbits = 0; for @circles Z (1,2,4...*) -> $c, $b {
$cbits += $b if $c.contains($p);
} $cbits;
}
my $xmin = [min] @circles.map: { .center.x - .radius } my $xmax = [max] @circles.map: { .center.x + .radius } my $ymin = [min] @circles.map: { .center.y - .radius } my $ymax = [max] @circles.map: { .center.y + .radius }
my $min-radius = @circles[*-1].radius;
my $outer-rect = Rect.new:
nw => Point.new(x => $xmin, y => $ymax), ne => Point.new(x => $xmax, y => $ymax), sw => Point.new(x => $xmin, y => $ymin), se => Point.new(x => $xmax, y => $ymin);
my $outer-area = $outer-rect.area;
my @unknowns = $outer-rect; my $known-dry = 0e0; my $known-wet = 0e0; my $div = 1;
- divide current rects each into four rects, analyze each
sub divide(@old) {
$div *= 2;
# rects too small to hold circle? my $smallish = @old[0].diag < $min-radius;
my @unk; for @old {
my $center = .nw mid .se; my $north = .nw mid .ne; my $south = .sw mid .se; my $west = .nw mid .sw; my $east = .ne mid .se;
for Rect.new(nw => .nw, ne => $north, sw => $west, se => $center), Rect.new(nw => $north, ne => .ne, sw => $center, se => $east), Rect.new(nw => $west, ne => $center, sw => .sw, se => $south), Rect.new(nw => $center, ne => $east, sw => $south, se => .se) { my @bits = .nw.cbits, .ne.cbits, .sw.cbits, .se.cbits;
# if all 4 points wet by same circle, guaranteed wet if [+&] @bits { $known-wet += .area; next; }
# if all 4 corners are dry, must check further if not [+|] @bits and $smallish {
# check that no circle bulges into this rect my $ok = True; for @circles -> $c { if .contains($c.east) or .contains($c.west) or .contains($c.north) or .contains($c.south) { $ok = False; last; } } if $ok { $known-dry += .area; next; } } push @unk, $_; # dunno yet }
} @unk;
}
my $delta = 0.001; repeat until my $diff < $delta {
@unknowns = divide(@unknowns);
$diff = $outer-area - $known-dry - $known-wet; say 'div: ', $div.fmt('%-5d'),
' unk: ', (+@unknowns).fmt('%-6d'), ' est: ', ($known-wet + $diff/2).fmt('%9.6f'), ' wet: ', $known-wet.fmt('%9.6f'), ' dry: ', ($outer-area - $known-dry).fmt('%9.6f'), ' diff: ', $diff.fmt('%9.6f'), ' error: ', ($diff - @unknowns * @unknowns[0].area).fmt('%e'); }</lang>
- Output:
div: 2 unk: 4 est: 14.607153 wet: 0.000000 dry: 29.214306 diff: 29.214306 error: 0.000000e+000 div: 4 unk: 15 est: 15.520100 wet: 1.825894 dry: 29.214306 diff: 27.388411 error: 0.000000e+000 div: 8 unk: 39 est: 20.313072 wet: 11.411838 dry: 29.214306 diff: 17.802467 error: 7.105427e-015 div: 16 unk: 107 est: 23.108972 wet: 17.003639 dry: 29.214306 diff: 12.210667 error: -1.065814e-014 div: 32 unk: 142 est: 21.368667 wet: 19.343066 dry: 23.394268 diff: 4.051203 error: 8.881784e-014 div: 64 unk: 290 est: 21.504182 wet: 20.469985 dry: 22.538380 diff: 2.068396 error: 1.771916e-013 div: 128 unk: 582 est: 21.534495 wet: 21.015613 dry: 22.053377 diff: 1.037764 error: -3.175238e-013 div: 256 unk: 1169 est: 21.557898 wet: 21.297343 dry: 21.818454 diff: 0.521111 error: -2.501332e-013 div: 512 unk: 2347 est: 21.563415 wet: 21.432636 dry: 21.694194 diff: 0.261558 error: -1.046996e-012 div: 1024 unk: 4700 est: 21.564111 wet: 21.498638 dry: 21.629584 diff: 0.130946 error: 1.481315e-013 div: 2048 unk: 9407 est: 21.564804 wet: 21.532043 dry: 21.597565 diff: 0.065522 error: 1.781700e-012 div: 4096 unk: 18818 est: 21.564876 wet: 21.548492 dry: 21.581260 diff: 0.032768 error: 1.098372e-011 div: 8192 unk: 37648 est: 21.564992 wet: 21.556797 dry: 21.573187 diff: 0.016389 error: -1.413968e-011 div: 16384 unk: 75301 est: 21.565017 wet: 21.560920 dry: 21.569115 diff: 0.008195 error: -7.683898e-011 div: 32768 unk: 150599 est: 21.565031 wet: 21.562982 dry: 21.567080 diff: 0.004097 error: -1.247991e-010 div: 65536 unk: 301203 est: 21.565035 wet: 21.564010 dry: 21.566059 diff: 0.002049 error: -2.830591e-010 div: 131072 unk: 602411 est: 21.565036 wet: 21.564524 dry: 21.565548 diff: 0.001024 error: -1.607121e-010
Here the "diff" is calculated by subtracting the known wet and dry areas from the total area, and the "error" is the difference between that and the sum of the areas of the unknown blocks, to give a rough idea of how much floating point roundoff error we've accumulated.
Python
Grid Sampling Version
This implements a regular grid sampling. For this problems this is more efficient than a Montecarlo sampling. <lang python>from collections import namedtuple
Circle = namedtuple("Circle", "x y r")
circles = [
Circle( 1.6417233788, 1.6121789534, 0.0848270516), Circle(-1.4944608174, 1.2077959613, 1.1039549836), Circle( 0.6110294452, -0.6907087527, 0.9089162485), Circle( 0.3844862411, 0.2923344616, 0.2375743054), Circle(-0.2495892950, -0.3832854473, 1.0845181219), Circle( 1.7813504266, 1.6178237031, 0.8162655711), Circle(-0.1985249206, -0.8343333301, 0.0538864941), Circle(-1.7011985145, -0.1263820964, 0.4776976918), Circle(-0.4319462812, 1.4104420482, 0.7886291537), Circle( 0.2178372997, -0.9499557344, 0.0357871187), Circle(-0.6294854565, -1.3078893852, 0.7653357688), Circle( 1.7952608455, 0.6281269104, 0.2727652452), Circle( 1.4168575317, 1.0683357171, 1.1016025378), Circle( 1.4637371396, 0.9463877418, 1.1846214562), Circle(-0.5263668798, 1.7315156631, 1.4428514068), Circle(-1.2197352481, 0.9144146579, 1.0727263474), Circle(-0.1389358881, 0.1092805780, 0.7350208828), Circle( 1.5293954595, 0.0030278255, 1.2472867347), Circle(-0.5258728625, 1.3782633069, 1.3495508831), Circle(-0.1403562064, 0.2437382535, 1.3804956588), Circle( 0.8055826339, -0.0482092025, 0.3327165165), Circle(-0.6311979224, 0.7184578971, 0.2491045282), Circle( 1.4685857879, -0.8347049536, 1.3670667538), Circle(-0.6855727502, 1.6465021616, 1.0593087096), Circle( 0.0152957411, 0.0638919221, 0.9771215985)]
def main():
# compute the bounding box of the circles x_min = min(c.x - c.r for c in circles) x_max = max(c.x + c.r for c in circles) y_min = min(c.y - c.r for c in circles) y_max = max(c.y + c.r for c in circles)
box_side = 500
dx = (x_max - x_min) / box_side dy = (y_max - y_min) / box_side
count = 0
for r in xrange(box_side): y = y_min + r * dy for c in xrange(box_side): x = x_min + c * dx if any((x-circle.x)**2 + (y-circle.y)**2 <= (circle.r ** 2) for circle in circles): count += 1
print "Approximated area:", count * dx * dy
main()</lang>
- Output:
Approximated area: 21.561559772
Scanline Conversion
<lang python>from math import floor, ceil, sqrt
def area_scan(prec, circs):
def sect((cx, cy, r), y): dr = sqrt(r ** 2 - (y - cy) ** 2) return (cx - dr, cx + dr)
ys = [a[1] + a[2] for a in circs] + [a[1] - a[2] for a in circs] mins = int(floor(min(ys) / prec)) maxs = int(ceil(max(ys) / prec))
total = 0 for y in (prec * x for x in xrange(mins, maxs + 1)): right = -float("inf")
for (x0, x1) in sorted(sect((cx, cy, r), y) for (cx, cy, r) in circs if abs(y - cr) < r): if x1 <= right: continue total += x1 - max(x0, right) right = x1
return total * prec
def main():
circles = [ ( 1.6417233788, 1.6121789534, 0.0848270516), (-1.4944608174, 1.2077959613, 1.1039549836), ( 0.6110294452, -0.6907087527, 0.9089162485), ( 0.3844862411, 0.2923344616, 0.2375743054), (-0.2495892950, -0.3832854473, 1.0845181219), ( 1.7813504266, 1.6178237031, 0.8162655711), (-0.1985249206, -0.8343333301, 0.0538864941), (-1.7011985145, -0.1263820964, 0.4776976918), (-0.4319462812, 1.4104420482, 0.7886291537), ( 0.2178372997, -0.9499557344, 0.0357871187), (-0.6294854565, -1.3078893852, 0.7653357688), ( 1.7952608455, 0.6281269104, 0.2727652452), ( 1.4168575317, 1.0683357171, 1.1016025378), ( 1.4637371396, 0.9463877418, 1.1846214562), (-0.5263668798, 1.7315156631, 1.4428514068), (-1.2197352481, 0.9144146579, 1.0727263474), (-0.1389358881, 0.1092805780, 0.7350208828), ( 1.5293954595, 0.0030278255, 1.2472867347), (-0.5258728625, 1.3782633069, 1.3495508831), (-0.1403562064, 0.2437382535, 1.3804956588), ( 0.8055826339, -0.0482092025, 0.3327165165), (-0.6311979224, 0.7184578971, 0.2491045282), ( 1.4685857879, -0.8347049536, 1.3670667538), (-0.6855727502, 1.6465021616, 1.0593087096), ( 0.0152957411, 0.0638919221, 0.9771215985)]
p = 1e-3 print "@stepsize", p, "area = %.4f" % area_scan(p, circles)
main()</lang>
2D Van der Corput sequence
![](http://static.miraheze.org/rosettacodewiki/thumb/7/71/Van_der_Corput_2D.png/200px-Van_der_Corput_2D.png)
Remembering that the Van der Corput sequence is used for Monte Carlo-like simulations. This example uses a Van der Corput sequence generator of base 2 for the first dimension, and base 3 for the second dimension of the 2D space which seems to cover evenly.
To aid in efficiency:
- Circles are uniquified,
- Sorted in descending order of size,
- Wholly obscured circles removed,
- And the square of the radius computed outside the main loops.
<lang python>from __future__ import division from math import sqrt from itertools import count from pprint import pprint as pp try:
from itertools import izip as zip
except:
pass
- Remove duplicates and sort, largest first.
circles = sorted(set([
# xcenter ycenter radius (1.6417233788, 1.6121789534, 0.0848270516), (-1.4944608174, 1.2077959613, 1.1039549836), (0.6110294452, -0.6907087527, 0.9089162485), (0.3844862411, 0.2923344616, 0.2375743054), (-0.2495892950, -0.3832854473, 1.0845181219), (1.7813504266, 1.6178237031, 0.8162655711), (-0.1985249206, -0.8343333301, 0.0538864941), (-1.7011985145, -0.1263820964, 0.4776976918), (-0.4319462812, 1.4104420482, 0.7886291537), (0.2178372997, -0.9499557344, 0.0357871187), (-0.6294854565, -1.3078893852, 0.7653357688), (1.7952608455, 0.6281269104, 0.2727652452), (1.4168575317, 1.0683357171, 1.1016025378), (1.4637371396, 0.9463877418, 1.1846214562), (-0.5263668798, 1.7315156631, 1.4428514068), (-1.2197352481, 0.9144146579, 1.0727263474), (-0.1389358881, 0.1092805780, 0.7350208828), (1.5293954595, 0.0030278255, 1.2472867347), (-0.5258728625, 1.3782633069, 1.3495508831), (-0.1403562064, 0.2437382535, 1.3804956588), (0.8055826339, -0.0482092025, 0.3327165165), (-0.6311979224, 0.7184578971, 0.2491045282), (1.4685857879, -0.8347049536, 1.3670667538), (-0.6855727502, 1.6465021616, 1.0593087096), (0.0152957411, 0.0638919221, 0.9771215985), ]), key=lambda x: -x[-1])
def vdcgen(base=2):
'Van der Corput sequence generator' for n in count(): vdc, denom = 0,1 while n: denom *= base n, remainder = divmod(n, base) vdc += remainder / denom yield vdc
def vdc_2d():
'Two dimensional Van der Corput sequence generator' for x, y in zip(vdcgen(base=2), vdcgen(base=3)): yield x, y
def bounding_box(circles):
'Return minx, maxx, miny, maxy' return (min(x - r for x,y,r in circles), max(x + r for x,y,r in circles), min(y - r for x,y,r in circles), max(y + r for x,y,r in circles) )
def circle_is_in_circle(c1, c2):
x1, y1, r1 = c1 x2, y2, r2 = c2 return sqrt((x2 - x1)**2 + (y2 - y1)**2) <= r1 - r2
def remove_covered_circles(circles):
'Takes circles in decreasing radius order. Removes those covered by others' covered = [] for i, c1 in enumerate(circles): eliminate = [c2 for c2 in circles[i+1:] if circle_is_in_circle(c1, c2)] if eliminate: covered += [c1, eliminate] for c in eliminate: circles.remove(c) #pp(covered)
def main(circles):
print('Originally %i circles' % len(circles)) print('Bounding box: %r' % (bounding_box(circles),)) remove_covered_circles(circles) print(' down to %i due to some being wholly covered by others' % len(circles)) minx, maxx, miny, maxy = bounding_box(circles) # Shift to 0,0 and compute r**2 once circles2 = [(x - minx, y - miny, r*r) for x, y, r in circles] scalex, scaley = abs(maxx - minx), abs(maxy - miny) pcount, inside, last = 0, 0, for px, py in vdc_2d(): pcount += 1 px *= scalex; py *= scaley if any((px-cx)**2 + (py-cy)**2 <= cr2 for cx, cy, cr2 in circles2): inside += 1 if not pcount % 100000: area = (inside/pcount) * scalex * scaley print('Points: %8i, Area estimate: %r' % (pcount, area)) # Hack to check if precision OK this = '%.4f' % area if this == last: break else: last = this print('The value has settled to %s' % this)
if __name__ == '__main__':
main(circles)</lang>
The above is tested to work with Python v.2.7, Python3 and PyPy.
- Output:
python3 total_circle_area.py Originally 25 circles Bounding box: (-2.598415801, 2.8356525417, -2.2017717074, 3.1743670698999997) down to 14 due to some being wholly covered by others Points: 100000, Area estimate: 21.57125892144117 Points: 200000, Area estimate: 21.565708203389384 Points: 300000, Area estimate: 21.56668201357391 Points: 400000, Area estimate: 21.566949811374652 Points: 500000, Area estimate: 21.567694776165812 Points: 600000, Area estimate: 21.566097727463195 Points: 700000, Area estimate: 21.565374325611838 Points: 800000, Area estimate: 21.565963828562822 Points: 900000, Area estimate: 21.56596788610526 The value has settled to 21.5660
- Comparison
As mentioned on the talk page, this version does more with its 200K points than the rectangular grid implementation does with its 250K points, (and the C coded Monte Carlo method does with 100 million points).
Analytical Solution
This requires the dmath module. Because of the usage of Decimal and trigonometric functions implemented in Python, this program takes few minutes to run.
<lang python>from collections import namedtuple from functools import partial from itertools import repeat, imap, izip from decimal import Decimal, getcontext
- Requires the egg: https://pypi.python.org/pypi/dmath/
from dmath import atan2, asin, sin, cos, pi as piCompute
getcontext().prec = 40 # Set FP precision. sqrt = Decimal.sqrt pi = piCompute() D2 = Decimal(2)
Vec = namedtuple("Vec", "x y") vcross = lambda (a, b), (c, d): a*d - b*c vdot = lambda (a, b), (c, d): a*c + b*d vadd = lambda (a, b), (c, d): Vec(a + c, b + d) vsub = lambda (a, b), (c, d): Vec(a - c, b - d) vlen = lambda x: sqrt(vdot(x, x)) vdist = lambda a, b: vlen(vsub(a, b)) vscale = lambda s, (x, y): Vec(x * s, y * s)
def vnorm(v):
l = vlen(v) return Vec(v.x / l, v.y / l)
vangle = lambda (x, y): atan2(y, x)
def anorm(a):
if a > pi: return a - pi * D2 if a < -pi: return a + pi * D2 return a
Circle = namedtuple("Circle", "x y r")
def circle_cross((x0, y0, r0), (x1, y1, r1)):
d = vdist(Vec(x0, y0), Vec(x1, y1)) if d >= r0 + r1 or d <= abs(r0 - r1): return []
s = (r0 + r1 + d) / D2 a = sqrt(s * (s - d) * (s - r0) * (s - r1)) h = D2 * a / d dr = Vec(x1 - x0, y1 - y0) dx = vscale(sqrt(r0 ** 2 - h ** 2), vnorm(dr)) ang = vangle(dr) if \ r0 ** 2 + d ** 2 > r1 ** 2 \ else pi + vangle(dr) da = asin(h / r0) return map(anorm, [ang - da, ang + da])
- Angles of the start and end points of the circle arc.
Angle2 = namedtuple("Angle2", "a1 a2")
Arc = namedtuple("Arc", "c aa")
arcPoint = lambda (x, y, r), a: \
vadd(Vec(x, y), Vec(r * cos(a), r * sin(a)))
arc_start = lambda (c, (a0, a1)): arcPoint(c, a0) arc_mid = lambda (c, (a0, a1)): arcPoint(c, (a0 + a1) / D2) arc_end = lambda (c, (a0, a1)): arcPoint(c, a1) arc_center = lambda ((x, y, r), _): Vec(x, y)
arc_area = lambda ((_0, _1, r), (a0, a1)): r ** 2 * (a1 - a0) / D2
def split_circles(cs):
cSplit = lambda (c, angs): \ imap(Arc, repeat(c), imap(Angle2, angs, angs[1:]))
# If an arc that was part of one circle is inside *another* circle, # it will not be part of the zero-winding path, so reject it. in_circle = lambda ((x0, y0), c), (x, y, r): \ c != Circle(x, y, r) and vdist(Vec(x0, y0), Vec(x, y)) < r
def in_any_circle(arc): return any(in_circle((arc_mid(arc), arc.c), c) for c in cs)
concat_map = lambda f, xs: [y for x in xs for y in f(x)]
f = lambda c: \ (c, sorted([-pi, pi] + concat_map(partial(circle_cross, c), cs))) cAngs = map(f, cs) arcs = concat_map(cSplit, cAngs) return filter(lambda ar: not in_any_circle(ar), arcs)
- Given a list of arcs, build sets of closed paths from them.
- If one arc's end point is no more than 1e-4 from another's
- start point, they are considered connected. Since these
- start/end points resulted from intersecting circles earlier,
- they *should* be exactly the same, but floating point
- precision may cause small differences, hence the 1e-4 error
- margin. When there are genuinely different intersections
- closer than this margin, the method will backfire, badly.
def make_paths(arcs):
eps = Decimal("0.0001") def join_arcs(a, xxs): if not xxs: return [a] x, xs = xxs[0], xxs[1:] if not a: return join_arcs([x], xs) if vdist(arc_start(a[0]), arc_end(a[-1])) < eps: return [a] + join_arcs([], xxs) if vdist(arc_end(a[-1]), arc_start(x)) < eps: return join_arcs(a + [x], xs) return join_arcs(a, xs + [x]) return join_arcs([], arcs)
- Slice N-polygon into N-2 triangles.
def polyline_area(vvs):
tri_area = lambda a, b, c: vcross(vsub(b, a), vsub(c, b)) / D2 v, vs = vvs[0], vvs[1:] return sum(tri_area(v, v1, v2) for v1, v2 in izip(vs, vs[1:]))
def path_area(arcs):
f = lambda (a, e), arc: \ (a + arc_area(arc), e + [arc_center(arc), arc_end(arc)]) (a, e) = reduce(f, arcs, (0, [])) return a + polyline_area(e)
circles_area = lambda cs: \
sum(imap(path_area, make_paths(split_circles(cs))))
def main():
raw_circles = """\ 1.6417233788 1.6121789534 0.0848270516 -1.4944608174 1.2077959613 1.1039549836 0.6110294452 -0.6907087527 0.9089162485 0.3844862411 0.2923344616 0.2375743054 -0.2495892950 -0.3832854473 1.0845181219 1.7813504266 1.6178237031 0.8162655711 -0.1985249206 -0.8343333301 0.0538864941 -1.7011985145 -0.1263820964 0.4776976918 -0.4319462812 1.4104420482 0.7886291537 0.2178372997 -0.9499557344 0.0357871187 -0.6294854565 -1.3078893852 0.7653357688 1.7952608455 0.6281269104 0.2727652452 1.4168575317 1.0683357171 1.1016025378 1.4637371396 0.9463877418 1.1846214562 -0.5263668798 1.7315156631 1.4428514068 -1.2197352481 0.9144146579 1.0727263474 -0.1389358881 0.1092805780 0.7350208828 1.5293954595 0.0030278255 1.2472867347 -0.5258728625 1.3782633069 1.3495508831 -0.1403562064 0.2437382535 1.3804956588 0.8055826339 -0.0482092025 0.3327165165 -0.6311979224 0.7184578971 0.2491045282 1.4685857879 -0.8347049536 1.3670667538 -0.6855727502 1.6465021616 1.0593087096 0.0152957411 0.0638919221 0.9771215985""".splitlines()
circles = [Circle(*imap(Decimal, row.split())) for row in raw_circles] print "Total Area:", circles_area(circles)
main()</lang>
- Output:
Total Area: 21.56503660385639895908422492887814801839
REXX
These REXX programs use the grid sampling method.
using all circles
<lang rexx>/*REXX pgm calculates the total area of (possibly overlapping) circles. */ parse arg box . /*obtain possible #boxes from CL.*/ if box== then box=500 /*allow specification of # boxes,*/ numeric digits 20 /*ensure enough digits for points*/
data = ' 1.6417233788 1.6121789534 0.0848270516', '-1.4944608174 1.2077959613 1.1039549836', ' 0.6110294452 -0.6907087527 0.9089162485', ' 0.3844862411 0.2923344616 0.2375743054', '-0.2495892950 -0.3832854473 1.0845181219', ' 1.7813504266 1.6178237031 0.8162655711', '-0.1985249206 -0.8343333301 0.0538864941', '-1.7011985145 -0.1263820964 0.4776976918', '-0.4319462812 1.4104420482 0.7886291537', ' 0.2178372997 -0.9499557344 0.0357871187', '-0.6294854565 -1.3078893852 0.7653357688', ' 1.7952608455 0.6281269104 0.2727652452', ' 1.4168575317 1.0683357171 1.1016025378', ' 1.4637371396 0.9463877418 1.1846214562', '-0.5263668798 1.7315156631 1.4428514068', '-1.2197352481 0.9144146579 1.0727263474', '-0.1389358881 0.1092805780 0.7350208828', ' 1.5293954595 0.0030278255 1.2472867347', '-0.5258728625 1.3782633069 1.3495508831', '-0.1403562064 0.2437382535 1.3804956588', ' 0.8055826339 -0.0482092025 0.3327165165', '-0.6311979224 0.7184578971 0.2491045282', ' 1.4685857879 -0.8347049536 1.3670667538', '-0.6855727502 1.6465021616 1.0593087096', ' 0.0152957411 0.0638919221 0.9771215985'
circles=words(data)%3 /* ══x══ ══y══ ══radius══ */ parse var data minX minY . 1 maxX maxY . /*assign some min & max vals.*/
do j=1 for circles; _=j*3-2 /*assign circles with datam. */ @x.j=word(data,_); @y.j=word(data,_+1) @r.j=word(data,_+2); @rr.j=@r.j**2 minX=min(minX, @x.j-@r.j); maxX=max(maxX, @x.j+@r.j) minY=min(minY, @y.j-@r.j); maxY=max(maxY, @y.j+@r.j) end /*j*/
dx=(maxX-minX) / box dy=(maxY-minY) / box
- =0 /*count of sample points. */
do row=0 to box; y=minY+row*dy /*process the grid rows. */ do col=0 to box; x=minX+col*dx /*process the grid cols. */ do k=1 for circles /*now process each circle. */ if (x-@x.k)**2+(y-@y.k)**2 <= @rr.k then do; #=#+1; leave; end end /*k*/ end /*col*/ end /*row*/ /*stick a fork in it, we're done.*/
say 'The approximate area is: ' #*dx*dy " using" box "boxes."</lang> outputs when using the default input:
The approximate area is: 21.561559772003320766 using 500 boxes.
eliding contained circles
This version elides any circle that is completely contained in another circle.
This version also has additional information displayed. <lang rexx>/*REXX pgm calculates the total area of (possibly overlapping) circles. */ parse arg box . /*obtain possible #boxes from CL.*/ if box== then box=-500 /*allow specification of # boxes,*/ verbose= box<0 /*set a flag if in verbose mode. */ box=abs(box); boxen=box+1 /*use |box| value from here on.*/ numeric digits 20 /*ensure enough digits for points*/
data = ' 1.6417233788 1.6121789534 0.0848270516', '-1.4944608174 1.2077959613 1.1039549836', ' 0.6110294452 -0.6907087527 0.9089162485', ' 0.3844862411 0.2923344616 0.2375743054', '-0.2495892950 -0.3832854473 1.0845181219', ' 1.7813504266 1.6178237031 0.8162655711', '-0.1985249206 -0.8343333301 0.0538864941', '-1.7011985145 -0.1263820964 0.4776976918', '-0.4319462812 1.4104420482 0.7886291537', ' 0.2178372997 -0.9499557344 0.0357871187', '-0.6294854565 -1.3078893852 0.7653357688', ' 1.7952608455 0.6281269104 0.2727652452', ' 1.4168575317 1.0683357171 1.1016025378', ' 1.4637371396 0.9463877418 1.1846214562', '-0.5263668798 1.7315156631 1.4428514068', '-1.2197352481 0.9144146579 1.0727263474', '-0.1389358881 0.1092805780 0.7350208828', ' 1.5293954595 0.0030278255 1.2472867347', '-0.5258728625 1.3782633069 1.3495508831', '-0.1403562064 0.2437382535 1.3804956588', ' 0.8055826339 -0.0482092025 0.3327165165', '-0.6311979224 0.7184578971 0.2491045282', ' 1.4685857879 -0.8347049536 1.3670667538', '-0.6855727502 1.6465021616 1.0593087096', ' 0.0152957411 0.0638919221 0.9771215985'
circles=words(data)%3 /* ══x══ ══y══ ══radius══ */ if verbose then say 'There are' circles "circles." parse var data minX minY . 1 maxX maxY . /*assign some min & max vals.*/
do j=1 for circles; _=j*3-2 /*assign circles with datam. */ @x.j=word(data,_); @y.j=word(data,_+1) @r.j=word(data,_+2); @rr.j=@r.j**2 minX=min(minX, @x.j-@r.j); maxX=max(maxX, @x.j+@r.j) minY=min(minY, @y.j-@r.j); maxY=max(maxY, @y.j+@r.j) end /*j*/
dx=(maxX-minX) / box dy=(maxY-minY) / box
- in=0 /*# fully contained circles. */
/* [↓] find contained circles*/ do j=1 for circles /*traipse through J circles*/ do k=1 for circles /* " " K " */ if k==j | @r.j==0 then iterate /*ignore self or zero radius.*/ if @x.j+@r.j > @x.k+@r.k |, /*is circle J outside cir K ?*/ @x.j-@r.j < @x.k-@r.k |, @y.j+@r.j > @y.k+@r.k |, @y.j-@r.j < @y.k-@r.k then iterate /*circle ¬ contained.*/ if verbose then say 'Circle ' j " is contained in circle " k @r.j=0; #in=#in+1 /*elide this circle; bump #in*/ end /*k*/ end /*j*/
if #in==0 then #in='no' /*use gooder English. (joke)*/ if verbose then say #in "circles that are totally within another circle." nC=0 /*number of "new" circles. */
do n=1 for circles; if @r.n==0 then iterate /*skip if 0.*/ nC=nC+1; @x.nC=@x.n; @y.nC=@y.n; @r.nC=@r.n; @rr.nC=@r.n**2 end /*n*/ /* [↑] elide overlapping cir*/
- =0 /*count of sample points. */
do row=0 for boxen; y=minY+row*dy /*process the grid rows. */ do col=0 for boxen; x=minX+col*dx /*process the grid cols. */ do k=1 for nC /*now process each new circle*/ if (x-@x.k)**2+(y-@y.k)**2 <= @rr.k then do; #=#+1; leave; end end /*k*/ end /*col*/ end /*row*/
if verbose then say /*stick a fork in it, we're done.*/ say 'The approximate area is: ' #*dx*dy " using" box "boxes."</lang> output using various number of boxes:
There are 25 circles. Circle 1 is contained in circle 6 Circle 4 is contained in circle 5 Circle 7 is contained in circle 3 Circle 9 is contained in circle 15 Circle 10 is contained in circle 3 Circle 12 is contained in circle 13 Circle 17 is contained in circle 20 Circle 21 is contained in circle 18 Circle 22 is contained in circle 15 Circle 24 is contained in circle 15 Circle 25 is contained in circle 20 11 circles that are totally within another circle. The approximate area is: 21.561559772003320766 using 500 boxes.
The approximate area is: 21.563838487835104416 using 1000 boxes.
The approximate area is: 21.564656488390103675 using 2000 boxes.
The approximate area is: 21.564912113563540943 using 4000 boxes.
The approximate area is: 21.565000669427195997 using 8000 boxes.
Tcl
This presents three techniques, a regular grid sampling, a pure Monte Carlo sampling, and a technique where there is a regular grid but then a vote of uniformly-distributed samples within each grid cell is taken to see if the cell is “in” or “out”. <lang tcl>set circles {
1.6417233788 1.6121789534 0.0848270516 -1.4944608174 1.2077959613 1.1039549836 0.6110294452 -0.6907087527 0.9089162485 0.3844862411 0.2923344616 0.2375743054 -0.2495892950 -0.3832854473 1.0845181219 1.7813504266 1.6178237031 0.8162655711 -0.1985249206 -0.8343333301 0.0538864941 -1.7011985145 -0.1263820964 0.4776976918 -0.4319462812 1.4104420482 0.7886291537 0.2178372997 -0.9499557344 0.0357871187 -0.6294854565 -1.3078893852 0.7653357688 1.7952608455 0.6281269104 0.2727652452 1.4168575317 1.0683357171 1.1016025378 1.4637371396 0.9463877418 1.1846214562 -0.5263668798 1.7315156631 1.4428514068 -1.2197352481 0.9144146579 1.0727263474 -0.1389358881 0.1092805780 0.7350208828 1.5293954595 0.0030278255 1.2472867347 -0.5258728625 1.3782633069 1.3495508831 -0.1403562064 0.2437382535 1.3804956588 0.8055826339 -0.0482092025 0.3327165165 -0.6311979224 0.7184578971 0.2491045282 1.4685857879 -0.8347049536 1.3670667538 -0.6855727502 1.6465021616 1.0593087096 0.0152957411 0.0638919221 0.9771215985
}
proc init {} {
upvar 1 circles circles set xMin [set yMin inf] set xMax [set yMax -inf] set i -1 foreach {xc yc rad} $circles {
set xMin [expr {min($xMin, $xc-$rad)}] set xMax [expr {max($xMax, $xc+$rad)}] set yMin [expr {min($yMin, $yc-$rad)}] set yMax [expr {max($yMax, $yc+$rad)}] lset circles [incr i 3] [expr {$rad**2}]
} return [list $xMin $xMax $yMin $yMax]
}
proc sampleGrid {circles steps} {
lassign [init] xMin xMax yMin yMax set dx [expr {($xMax-$xMin)/$steps}] set dy [expr {($yMax-$yMin)/$steps}] set n 0 for {set i 0} {$i < $steps} {incr i} {
set x [expr {$xMin + $i * $dx}] for {set j 0} {$j < $steps} {incr j} { set y [expr {$yMin + $j * $dy}] foreach {xc yc rad2} $circles { if {($x-$xc)**2 + ($y-$yc)**2 <= $rad2} { incr n break } } }
} return [expr {$dx * $dy * $n}]
}
proc sampleMC {circles samples} {
lassign [init] xMin xMax yMin yMax set n 0 for {set i 0} {$i < $samples} {incr i} {
set x [expr {$xMin+rand()*($xMax-$xMin)}] set y [expr {$yMin+rand()*($yMax-$yMin)}] foreach {xc yc rad2} $circles { if {($x-$xc)**2 + ($y-$yc)**2 <= $rad2} { incr n break } }
} return [expr {($xMax-$xMin) * ($yMax-$yMin) * $n / $samples}]
}
proc samplePerturb {circles steps votes} {
lassign [init] xMin xMax yMin yMax set dx [expr {($xMax-$xMin)/$steps}] set dy [expr {($yMax-$yMin)/$steps}] set n 0 for {set i 0} {$i < $steps} {incr i} {
set x [expr {$xMin + $i * $dx}] for {set j 0} {$j < $steps} {incr j} { set y [expr {$yMin + $j * $dy}] foreach {xc yc rad2} $circles { set in 0 for {set v 0} {$v < $votes} {incr v} { set xr [expr {$x + (rand()-0.5)*$dx}] set yr [expr {$y + (rand()-0.5)*$dy}] if {($xr-$xc)**2 + ($yr-$yc)**2 <= $rad2} { incr in } } if {$in*2 >= $votes} { incr n break } } }
} return [expr {$dx * $dy * $n}]
}
puts [format "estimated area (grid): %.4f" [sampleGrid $circles 500]] puts [format "estimated area (monte carlo): %.2f" [sampleMC $circles 1000000]] puts [format "estimated area (perturbed sample): %.4f" [samplePerturb $circles 500 5]]</lang>
- Output:
estimated area (grid): 21.5616 estimated area (monte carlo): 21.56 estimated area (perturbed sample): 21.5645
Note that the error on the Monte Carlo sampling is actually very high; the above run happened to deliver a figure closer to the real value than usual.