Total circles area

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

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