3

I plan on writing this code in python, but this is more of a general algorithm-design problem than one having to do with any particular language. I'm writing code to simulate a hybrid rocket engine, and long story short, the problem involves finding both the perimeter and the area of a figure composed of many (possibly thousands) of overlapping circles.

I saw this on stackexchange:

Combined area of overlapping circles

except I need to find the perimeter also. Someone in that thread mentioned Monte Carlo (random point guess) methods, but can that be used to find the perimeter as well as the area?

Thanks in advance.

Community
  • 1
  • 1
user2948994
  • 39
  • 2
  • The random sample method for area works well because you only care about the proportion of points in and out of the shape. The same thing wouldn't work for the perimeter as too few points would fall exactly on it and the "area of the perimeter" has no meaning. I think you would have to work out the lengths of perimeter arcs based on points of overlap. – jonrsharpe Nov 24 '13 at 23:00
  • http://math.stackexchange.com is your site for this! – shad0w_wa1k3r Nov 24 '13 at 23:39
  • Have you seen [this](http://stackoverflow.com/questions/13517211/python-how-calculate-a-polygon-perimeter-using-an-osgeo-ogr-geometry-object) question? It might point you in the right direction. – Steinar Lima Nov 24 '13 at 23:49

2 Answers2

1

Let's say your have the circles c1, c2, ..., cn.

Take c1 and intersect it with each other circle. Initialize an empty list of circle sectors for c1.

If the intersection is zero points or one point and c1 is outside the other circle, then add a 360 degree sector to the list.

If the intersection is zero points or one point and c1 is inside the other circle, then the effective outline of c1 is 0, stop processing c1 with an outline of 0 and take the next circle.

If the intersection is two points, add the sector of c1, which is not in the other circle to the list of circle sectors of c1.

The effective outline of c1, is the length of the intersection of all sectors within the list of sectors of c1. This intersection of sectors can consist of various disjoint sectors.

Rinse and repeat for all circles and then sum up the resulting values.

Hyperboreus
  • 31,997
  • 9
  • 47
  • 87
1

I am adding a second answer instead of expanding the first one, because I am quite positive that the other answer is correct, but I am not so sure if this answer is as correct. I have done some simple testing, and it seems to work, but please point out my errors.

This is basically a quick-n-dirty implementation of what I stated before:

from math import atan2, pi

def intersectLine (a, b):
    s0, e0, s1, e1 = a [0], a [1], b [0], b [1]
    s = max (s0, s1)
    e = min (e0, e1)
    if e <= s: return (0, 0)
    return (s, e)

class SectorSet:
    def __init__ (self, sectors):
        self.sectors = sectors [:]

    def __repr__ (self):
        return repr (self.sectors)

    def __iadd__ (self, other):
        acc = []
        for mine in self.sectors:
            for others in other.sectors:
                acc.append (intersectLine (mine, others) )
        self.sectors = [x for x in acc if x [0] != x [1] ]
        return self

class Circle:
    CONTAINS = 0
    CONTAINED = 1
    DISJOINT = 2
    INTERSECT = 3

    def __init__ (self, x, y, r):
        self.x = float (x)
        self.y = float (y)
        self.r = float (r)

    def intersect (self, other):
        a, b, c, d, r0, r1 = self.x, self.y, other.x, other.y, self.r, other.r
        r0sq, r1sq = r0 ** 2, r1 ** 2
        Dsq = (c - a) ** 2 + (d - b) ** 2
        D = Dsq ** .5
        if D >= r0 + r1:
            return Circle.DISJOINT, None, None
        if D <= abs (r0 - r1):
            return Circle.CONTAINED if r0 < r1 else Circle.CONTAINS, None, None
        dd = .25 * ( (D + r0 + r1) * (D + r0 - r1) * (D - r0 + r1) * (-D + r0 + r1) ) ** .5
        x = (a + c) / 2. + (c - a) * (r0sq - r1sq) / 2. / Dsq
        x1 = x + 2 * (b - d) / Dsq * dd
        x2 = x - 2 * (b - d) / Dsq * dd
        y = (b + d) / 2. + (d - b) * (r0sq - r1sq) / 2. / Dsq
        y1 = y - 2 * (a - c) / Dsq * dd
        y2 = y + 2 * (a - c) / Dsq * dd
        return Circle.INTERSECT, (x1, y1), (x2, y2)

    def angle (self, point):
        x0, y0, x, y = self.x, self.y, point [0], point [1]
        a = atan2 (y - y0, x - x0) + 1.5 * pi
        if a >= 2 * pi: a -= 2 * pi
        return a / pi * 180

    def sector (self, other):
        typ, i1, i2 = self.intersect (other)
        if typ == Circle.DISJOINT: return SectorSet ( [ (0, 360) ] )
        if typ == Circle.CONTAINS: return SectorSet ( [ (0, 360) ] )
        if typ == Circle.CONTAINED: return SectorSet ( [] )
        a1 = self.angle (i1)
        a2 = self.angle (i2)
        if a1 > a2: return SectorSet ( [ (0, a2), (a1, 360) ] )
        return SectorSet ( [ (a1, a2) ] )

    def outline (self, others):
        sectors = SectorSet ( [ (0, 360) ] )
        for other in others:
            sectors += self.sector (other)
        u = 2 * self.r * pi
        total = 0
        for start, end in sectors.sectors:
            total += (end - start) / 360. * u
        return total

def outline (circles):
    total = 0
    for circle in circles:
        others = [other for other in circles if circle != other]
        total += circle.outline (others)
    return total

a = Circle (0, 0, 2)
b = Circle (-2, -1, 1)
c = Circle (2, -1, 1)
d = Circle (0, 2, 1)
print (outline ( [a, b, c, d] ) )

Formula for calculating the intersecting points of two circles shamelessly stolen from here.

Hyperboreus
  • 31,997
  • 9
  • 47
  • 87