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