159 lines
5.4 KiB
Python
159 lines
5.4 KiB
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()
|