Total circles area: Difference between revisions

From Rosetta Code
Content added Content deleted
(→‎{{header|Haskell}}: analytical solution)
Line 374: Line 374:


in_circle ((x0,y0),c) (x,y,r) = c /= (x,y,r) && vdist (x0,y0) (x,y) < r
in_circle ((x0,y0),c) (x,y,r) = c /= (x,y,r) && vdist (x0,y0) (x,y) < r
in_any_circle arc = or $ map (in_circle (arc_mid arc, c)) cs
in_any_circle arc = any (in_circle (arc_mid arc, c)) cs
where (c,_) = arc
where (c,_) = arc


Revision as of 07:04, 18 September 2012

Total circles area is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

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): 

      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

According to one algorithm, the approximate solution is 21.56503660.

Beside solving the standard dataset, optionally solve a larger random dataset.

See also (idea originally from Steve132): http://www.reddit.com/r/dailyprogrammer/comments/zff9o/9062012_challenge_96_difficult_water_droplets/

C

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>

  1. include <stdlib.h>

typedef struct { double x, y, r; } Circle;

const 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}};

double min(const double a, const double b) { return a <= b ? a : b; }

double max(const double a, const double b) { return a >= b ? a : b; }

double uniform(const double a, const double b) { 

   const double r01 = rand() / (double)RAND_MAX;
   return a + (b - a) * r01;

}

int main() { 

   const size_t n_circles = sizeof(circles) / sizeof(Circle);

   // Initialize the bounding box of the circles.
   double x_min =  1e100;
   double x_max = -1e100;
   double y_min =  1e100;
   double y_max = -1e100;

   // Compute the bounding box of the circles.
   for (size_t i = 0; i < n_circles; i++) {
       const 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);
   }

   // Montecarlo sampling.
   srand(1);
   const size_t n_samples = 100 * 1000 * 1000;

   size_t hits = 0;

   for (size_t i = 0; i < n_samples; i++) {
       const double x = uniform(x_min, x_max);
       const double y = uniform(y_min, y_max);
       for (size_t j = 0; j < n_circles; j++) {
           const double dx = x - circles[j].x;
           const double dy = y - circles[j].y;
           if ((dx * dx + dy * dy) <= (circles[j].r * circles[j].r)) {
               hits++;
               break;
           }
       }
   }

   printf("Approximated area: %.8f\n",
          (double)(x_max - x_min) *
          (double)(y_max - y_min) *
          ((double)hits / n_samples));

   return 0;

}</lang>

Output:
Approximated area: 21.56262288

D

This version converges much faster than both the ordered grid and Montecarlo sampling solutions. <lang d>import std.stdio, std.math, std.algorithm, std.typecons, std.range;

alias real Fp; struct Circle { Fp x, y, r; }

void removeInternalDisks(ref Circle[] circles) { 

   static bool isFullyInternal(in Circle c1, in Circle c2)
   pure nothrow {
       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}();

   // What circles will be kept.
   auto keep = new bool[circles.length];
   keep[] = true;

   foreach (i1, c1; circles)
       if (keep[i1])
           foreach (i2, c2; circles)
               if (keep[i2])
                   if (i1 != i2 && isFullyInternal(c2, c1))
                       keep[i2] = false;

   // Pack circles array, removing fully internal circles.
   size_t pos = 0;
   foreach (i, k; keep)
       if (k)
           circles[pos++] = circles[i];
   circles.length = pos;
   // Alternative implementation of the packing:
   // circles = zip(circles, keep)
   //           .filter!(ck => ck[1])()
   //           .map!(ck => ck[0])()
   //           .array();

}


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))
                              (cast(Fp)0, circles[]);

   alias Tuple!(Fp,"y0", Fp,"y1") YRange;
   auto yRanges = new YRange[circles.length];

   Fp computeTotalArea(in Fp nSlicesX) {
       Fp total = 0;

       // Adapted from an idea by Cosmologicon.
       foreach (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 (ref const Circle c; circles) {
               immutable Fp d = c.r ^^ 2 - (c.x - x) ^^ 2;
               immutable Fp sd = sqrt(d);
               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 (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

Runtime is about 7.6 seconds. DMD 2.061, -O -release -inline -noboundscheck.

Haskell

Translation of: Python

<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 cs = any (isInside x y) cs

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]
     where
       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>import Data.List (sort)

vcross (a,b) (c,d) = a*d - b*c vdot (a,b) (c,d) = a*c + b*d vadd (a,b) (c,d) = (a+c, b+d) vsub (a,b) (c,d) = (a-c, b-d) vminus (a,b) (c,d) = (a-c, b-d) vlen2 x = vdot x x vlen x = sqrt $ vlen2 x vdist a b = vlen (a `vminus` b) vscale s (x,y) = (x*s, y*s) vnorm (x,y) = (x/l, y/l) where l = vlen (x,y) vangle (x,y) = atan2 y x anorm a | a > pi = a - pi * 2 | a < (-pi) = a + pi * 2 | otherwise = a

circle_cross (x0,y0,r0) (x1,y1,r1) | d >= r0 + r1 || d <= abs(r0 - r1) = [] | otherwise = map anorm [ang - da, ang + da] where d = vdist (x0,y0) (x1,y1) s = (r0 + r1 + d)/2 a = sqrt $ s * (s - d) * (s - r0) * (s - r1) h = 2 * a / d dr = (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 pi + (vangle dr) da = asin (h / r0)

arc_point (x,y,r) a = vadd (x,y) (r * cos a, r * sin a) arc_start (c,(a0,a1)) = arc_point c a0 arc_mid (c,(a0,a1)) = arc_point c ((a0+a1)/2) arc_end (c,(a0,a1)) = arc_point c a1 arc_center((x,y,r),_) = (x,y) arc_area ((_,_,r),(a0,a1)) = r * r * (a1 - a0) / 2

split_circles cs = filter (not.in_any_circle) arcs where arcs = concatMap csplit circs circs = map f cs where f c = (c, sort $ [-pi,pi] ++ (concatMap (circle_cross c) cs)) csplit (c, angs) = zip (repeat c) $ zip angs $ tail angs

in_circle ((x0,y0),c) (x,y,r) = c /= (x,y,r) && vdist (x0,y0) (x,y) < r in_any_circle arc = any (in_circle (arc_mid arc, c)) cs where (c,_) = arc

make_paths arcs = join_arcs [] arcs where join_arcs a [] = [a] join_arcs [] (x:xs) = join_arcs [x] xs join_arcs a (x:xs) | vdist (arc_start (head a)) (arc_end (last a)) < 1e-3 = a : join_arcs [] (x:xs) | vdist (arc_end (last a)) (arc_start x) < 1e-3 = join_arcs (a++[x]) xs | otherwise = join_arcs a (xs++[x])

path_area arcs = area_ 0 [] arcs where area_ a e [] = a + polyline_area e area_ a e (arc:as) = area_ (a + arc_area arc) (e++[arc_center arc, arc_end arc]) as

polyline_area [a,b] = 0 polyline_area (a:b:c:rest) = tri_area a b c + polyline_area (a:c:rest) where tri_area a b c = (b `vsub` a) `vcross` (c `vsub` b) / 2

circles_area = sum . map path_area . make_paths . split_circles

main = print $ circles_area [ ( 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)]</lang>

Output:
21.565036603856395

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

  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.01; 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

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

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
           for circle in circles:
               if (x-circle.x)**2 + (y-circle.y)**2 <= (circle.r ** 2):
                   count += 1
                   break

   print "Approximated area:", count * dx * dy

main()</lang>

Output:
Approximated area: 21.561559772

Scanline conversion: <lang python>from math import floor, ceil, sqrt

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)]

def area_scan(prec, circs): 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))

def sect((cx, cy, r), y): dr = sqrt(r ** 2 - (y - cy) ** 2) return (cx - dr, cx + dr)

total = 0 for y in [prec * x for x in range(mins, maxs + 1)]: right = -float("inf")

for (x0, x1) in sorted([sect((cx, cy, r), y) for (cx, cy, r) in circs if cy - r < y and y < cy + r ]):

if x1 <= right: continue total += x1 - max(x0, right) right = x1

return total * prec

p = 1e-3 print "@stepsize", p, "area = %.4f" % area_scan(p, circles)</lang>

Retrieved from "https://rosettacode.org/wiki/Total_circles_area?oldid=219405"