169 lines
5.8 KiB
Haskell
169 lines
5.8 KiB
Haskell
{-# LANGUAGE GeneralizedNewtypeDeriving #-}
|
|
|
|
import Data.List (sort)
|
|
|
|
data Vec = Vec Double Double
|
|
|
|
vCross, vDot :: Vec -> Vec -> Double
|
|
vCross (Vec a b) (Vec c d) = a*d - b*c
|
|
vDot (Vec a b) (Vec c d) = a*c + b*d
|
|
|
|
vAdd, vSub :: Vec -> Vec -> Vec
|
|
vAdd (Vec a b) (Vec c d) = Vec (a + c) (b + d)
|
|
vSub (Vec a b) (Vec c d) = Vec (a - c) (b - d)
|
|
|
|
vLen :: Vec -> Double
|
|
vLen x = sqrt $ vDot x x
|
|
|
|
vDist :: Vec -> Vec -> Double
|
|
vDist a b = vLen (a `vSub` b)
|
|
|
|
vScale :: Double -> Vec -> Vec
|
|
vScale s (Vec x y) = Vec (x * s) (y * s)
|
|
|
|
vNorm :: Vec -> Vec
|
|
vNorm v@(Vec x y) = Vec (x / l) (y / l) where l = vLen v
|
|
|
|
|
|
newtype Angle = A Double
|
|
deriving (Eq, Ord, Num, Fractional)
|
|
|
|
aPi = A pi
|
|
|
|
vAngle :: Vec -> Angle
|
|
vAngle (Vec x y) = A (atan2 y x)
|
|
|
|
aNorm :: Angle -> Angle
|
|
aNorm a | a > aPi = a - aPi * 2
|
|
| a < -aPi = a + aPi * 2
|
|
| otherwise = a
|
|
|
|
|
|
data Circle = Circle Double Double Double
|
|
deriving (Eq)
|
|
|
|
circleCross :: Circle -> Circle -> [Angle]
|
|
circleCross (Circle x0 y0 r0) (Circle x1 y1 r1)
|
|
| d >= r0 + r1 || d <= abs(r0 - r1) = []
|
|
| otherwise = map aNorm [ang - da, ang + da]
|
|
where
|
|
d = vDist (Vec x0 y0) (Vec x1 y1)
|
|
s = (r0 + r1 + d) / 2
|
|
a = sqrt $ s * (s - d) * (s - r0) * (s - r1)
|
|
h = 2 * a / d
|
|
dr = Vec (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 aPi + vAngle dr
|
|
da = A (asin (h / r0))
|
|
|
|
|
|
-- Angles of the start and end points of the circle arc.
|
|
data Angle2 = Angle2 Angle Angle
|
|
|
|
data Arc = Arc Circle Angle2
|
|
|
|
arcPoint :: Circle -> Angle -> Vec
|
|
arcPoint (Circle x y r) (A a) =
|
|
vAdd (Vec x y) (Vec (r * cos a) (r * sin a))
|
|
|
|
arcStart, arcMid, arcEnd, arcCenter :: Arc -> Vec
|
|
arcStart (Arc c (Angle2 a0 a1)) = arcPoint c a0
|
|
arcMid (Arc c (Angle2 a0 a1)) = arcPoint c ((a0 + a1) / 2)
|
|
arcEnd (Arc c (Angle2 a0 a1)) = arcPoint c a1
|
|
arcCenter (Arc (Circle x y r) _) = Vec x y
|
|
|
|
arcArea :: Arc -> Double
|
|
arcArea (Arc (Circle _ _ r) (Angle2 a0 a1)) = r ^ 2 * aDiff / 2
|
|
where (A aDiff) = a1 - a0
|
|
|
|
|
|
splitCircles :: [Circle] -> [Arc]
|
|
splitCircles cs = filter (not . inAnyCircle) arcs where
|
|
cSplit :: (Circle, [Angle]) -> [Arc]
|
|
cSplit (c, angs) =
|
|
zipWith Arc (repeat c) $ zipWith Angle2 angs $ tail angs
|
|
|
|
-- 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.
|
|
inCircle :: (Vec, Circle) -> Circle -> Bool
|
|
inCircle (Vec x0 y0, c1) c2@(Circle x y r) =
|
|
c1 /= c2 && vDist (Vec x0 y0) (Vec x y) < r
|
|
|
|
f :: Circle -> (Circle, [Angle])
|
|
f c = (c, sort $ [-aPi, aPi] ++ (concatMap (circleCross c) cs))
|
|
cAngs = map f cs
|
|
arcs = concatMap cSplit cAngs
|
|
|
|
inAnyCircle :: Arc -> Bool
|
|
inAnyCircle arc@(Arc c _) = any (inCircle (arcMid arc, c)) cs
|
|
|
|
|
|
{-
|
|
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.
|
|
-}
|
|
makePaths :: [Arc] -> [[Arc]]
|
|
makePaths arcs = joinArcs [] arcs where
|
|
joinArcs :: [Arc] -> [Arc] -> [[Arc]]
|
|
joinArcs a [] = [a]
|
|
joinArcs [] (x:xs) = joinArcs [x] xs
|
|
joinArcs a (x:xs)
|
|
| vDist (arcStart (head a)) (arcEnd (last a)) < 1e-4
|
|
= a : joinArcs [] (x:xs)
|
|
| vDist (arcEnd (last a)) (arcStart x) < 1e-4
|
|
= joinArcs (a ++ [x]) xs
|
|
| otherwise = joinArcs a (xs ++ [x])
|
|
|
|
|
|
pathArea :: [Arc] -> Double
|
|
pathArea arcs = a + polylineArea e where
|
|
(a, e) = foldl f (0, []) arcs
|
|
f (a, e) arc = (a + arcArea arc, e ++ [arcCenter arc, arcEnd arc])
|
|
|
|
|
|
-- Slice N-polygon into N-2 triangles.
|
|
polylineArea :: [Vec] -> Double
|
|
polylineArea (v:vs) = sum $ zipWith (triArea v) vs (tail vs)
|
|
where triArea a b c = ((b `vSub` a) `vCross` (c `vSub` b)) / 2
|
|
|
|
|
|
circlesArea :: [Circle] -> Double
|
|
circlesArea = sum . map pathArea . makePaths . splitCircles
|
|
|
|
|
|
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 = print $ circlesArea circles
|