From 2cbd644c59a03261fe83c4b93975570fc398443c Mon Sep 17 00:00:00 2001 From: Lemmih Date: Mon, 31 Aug 2020 13:59:28 +0000 Subject: [PATCH] =?UTF-8?q?Deploying=20to=20gh-pages=20from=20=20@=201f7c0?= =?UTF-8?q?e705cf581b29e99920bb434c12086f08d0c=20=F0=9F=9A=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- hpc_index.html | 6 +- hpc_index_alt.html | 6 +- hpc_index_exp.html | 6 +- hpc_index_fun.html | 6 +- playground/snippets.js | 2 +- .../Reanimate.Math.Polygon.hs.html | 1491 ++++++++--------- .../Reanimate.Math.SSSP.hs.html | 700 ++++---- .../Reanimate.PolyShape.hs.html | 614 ++++--- 8 files changed, 1415 insertions(+), 1416 deletions(-) diff --git a/hpc_index.html b/hpc_index.html index 07e492d..39ffa3e 100644 --- a/hpc_index.html +++ b/hpc_index.html @@ -68,10 +68,10 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/26
0%0/12
0%0/407
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon -8%7/82
5%4/79
2%58/2494
+8%7/81
5%4/79
2%58/2488
  module reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP -0%0/16
0%0/43
0%0/854
+0%0/19
0%0/53
0%0/911
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Triangulate 0%0/5
- 0/0 0%0/139
@@ -125,5 +125,5 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 100%6/6
50%1/2
95%43/45
  Program Coverage Total -31%248/798
16%131/818
30%4682/15556
+31%248/800
15%131/828
29%4682/15607
diff --git a/hpc_index_alt.html b/hpc_index_alt.html index 0341f45..2e6ed3e 100644 --- a/hpc_index_alt.html +++ b/hpc_index_alt.html @@ -59,7 +59,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 20%20/99
5%2/36
28%102/356
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon -8%7/82
5%4/79
2%58/2494
+8%7/81
5%4/79
2%58/2488
  module reanimate-0.4.3.0-inplace/Reanimate.Morph.Common 33%5/15
5%1/20
40%133/325
@@ -86,7 +86,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/26
0%0/12
0%0/407
  module reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP -0%0/16
0%0/43
0%0/854
+0%0/19
0%0/53
0%0/911
  module reanimate-0.4.3.0-inplace/Reanimate.Misc 0%0/7
0%0/14
0%0/162
@@ -125,5 +125,5 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 15%3/20
- 0/0 14%7/49
  Program Coverage Total -31%248/798
16%131/818
30%4682/15556
+31%248/800
15%131/828
29%4682/15607
diff --git a/hpc_index_exp.html b/hpc_index_exp.html index 7c8bac2..47f9905 100644 --- a/hpc_index_exp.html +++ b/hpc_index_exp.html @@ -92,7 +92,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 8%1/12
1%1/55
3%9/246
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon -8%7/82
5%4/79
2%58/2494
+8%7/81
5%4/79
2%58/2488
  module reanimate-0.4.3.0-inplace/Reanimate.Cache 0%0/8
0%0/12
0%0/160
@@ -113,7 +113,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/26
0%0/12
0%0/407
  module reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP -0%0/16
0%0/43
0%0/854
+0%0/19
0%0/53
0%0/911
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Triangulate 0%0/5
- 0/0 0%0/139
@@ -125,5 +125,5 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/1
0%0/4
0%0/61
  Program Coverage Total -31%248/798
16%131/818
30%4682/15556
+31%248/800
15%131/828
29%4682/15607
diff --git a/hpc_index_fun.html b/hpc_index_fun.html index e03d55e..0beadcc 100644 --- a/hpc_index_fun.html +++ b/hpc_index_fun.html @@ -89,7 +89,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 8%1/12
1%1/55
3%9/246
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon -8%7/82
5%4/79
2%58/2494
+8%7/81
5%4/79
2%58/2488
  module reanimate-0.4.3.0-inplace/Reanimate.Render 5%1/19
0%0/53
4%39/852
@@ -113,7 +113,7 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/26
0%0/12
0%0/407
  module reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP -0%0/16
0%0/43
0%0/854
+0%0/19
0%0/53
0%0/911
  module reanimate-0.4.3.0-inplace/Reanimate.Math.Triangulate 0%0/5
- 0/0 0%0/139
@@ -125,5 +125,5 @@ table.dashboard { border-collapse: collapse ; border: solid 1px black } 0%0/1
0%0/4
0%0/61
  Program Coverage Total -31%248/798
16%131/818
30%4682/15556
+31%248/800
15%131/828
29%4682/15607
diff --git a/playground/snippets.js b/playground/snippets.js index 71b43a1..a4ab53f 100644 --- a/playground/snippets.js +++ b/playground/snippets.js @@ -9,4 +9,4 @@ const snippets = [{"title": "Hello World","url": "https://reanimate.clozecards.c ,{"title": "Object Positions","url": "https://reanimate.clozecards.com/OdYISph5tKE/195.svg","code": "env =\n addStatic (mkBackground \"white\") .\n mapA (withStrokeColor \"black\")\n\nanimation :: Animation\nanimation = env $\n sceneAnimation $ do\n -- Configure objects\n txt <- newText \"Center\"\n top <- newText \"Top\"\n oModifyS top $ \n oTopY .= screenTop\n topR <- newText \"Top right\"\n oModifyS topR $ do\n oTopY .= screenTop\n oRightX .= screenRight\n botR <- newText \"Bottom right\"\n oModifyS botR $ do\n oTranslate .= (0, screenBottom+0.5)\n oRightX .= screenRight\n botL <- newText \"Bottom left\"\n oModifyS botL $ do\n oTranslate .= (0, screenBottom+0.5)\n oLeftX .= screenLeft\n topL <- newText \"Top left\"\n oModifyS topL $ do\n oTopY .= screenTop\n oLeftX .= screenLeft\n -- Show objects\n oShow txt\n wait 1\n switchTo txt top\n switchTo top topR\n switchTo topR botR\n switchTo botR botL\n switchTo botL topL\n switchTo topL txt\n\nswitchTo src dst = do\n fork $ oFadeOut src 1\n oModify dst $ oOpacity .~ 1\n oFadeIn dst 1\n wait 1\n\nnewText txt =\n newObject $ scale 1.5 $ centerX $ latex txt\n"} ,{"title": "Camera","url": "https://reanimate.clozecards.com/Hcx00P+aeph/150.svg","code": "animation :: Animation\nanimation = docEnv $ mapA (withFillOpacity 1) $ sceneAnimation $ do\n cam <- newObject Camera\n\n txt <- newObject $ center $ latex \"Fixed (non-cam)\"\n oModifyS txt $ do\n oTopY .= screenTop \n oZIndex .= 2\n\n circle <- newObject $ withFillColor \"blue\" $ mkCircle 1\n cameraAttach cam circle\n circleRight <- oRead circle oRightX\n\n box <- newObject $ withFillColor \"green\" $ mkRect 2 2\n cameraAttach cam box\n oModify box $ oLeftX .~ circleRight\n boxCenter <- oRead box oCenterXY\n\n small <- newObject $ center $ latex \"This text is very small\"\n cameraAttach cam small\n oModifyS small $ do\n oCenterXY .= boxCenter\n oScale .= 0.1\n \n oShow txt\n oShow small\n oShow circle\n oShow box\n\n wait 1\n\n cameraFocus cam boxCenter\n waitOn $ do\n fork $ cameraPan cam 3 boxCenter\n fork $ cameraZoom cam 3 15\n \n wait 2\n cameraZoom cam 3 1\n cameraPan cam 1 (0,0)\n"} ]; -const playgroundVersion = "2020-08-29 (1ec6d)"; +const playgroundVersion = "2020-08-31 (1f7c0)"; diff --git a/reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon.hs.html b/reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon.hs.html index f3b5449..fdac72b 100644 --- a/reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon.hs.html +++ b/reanimate-0.4.3.0-inplace/Reanimate.Math.Polygon.hs.html @@ -68,756 +68,751 @@ span.spaces { background: white } 49 -- * Single-Source-Shortest-Path 50 , ssspVisibility -- :: Polygon -> Polygon 51 , ssspWindows -- :: Polygon -> [(V2 Rational, V2 Rational)] - 52 -- * Duals - 53 , pdualPolygons -- :: Polygon -> PDual -> [Polygon] - 54 -- * Built-in shapes for testing - 55 , triangle -- :: Polygon - 56 , triangle' -- :: [P] - 57 , shape1 -- :: Polygon - 58 , shape2 -- :: Polygon - 59 , shape3 -- :: Polygon - 60 , shape4 -- :: Polygon - 61 , shape5 -- :: Polygon - 62 , shape6 -- :: Polygon - 63 , shape7 -- :: Polygon - 64 , shape8 -- :: Polygon - 65 , shape9 -- :: Polygon - 66 , shape10 -- :: Polygon - 67 , shape11 -- :: Polygon - 68 , shape12 -- :: Polygon - 69 , shape13 -- :: Polygon - 70 , shape14 -- :: Polygon - 71 , shape15 -- :: Polygon - 72 , shape16 -- :: Polygon - 73 , shape17 -- :: Polygon - 74 , shape18 -- :: Polygon - 75 , shape19 -- :: Polygon - 76 , shape20 -- :: Polygon - 77 , shape21 -- :: Polygon - 78 , shape22 -- :: Polygon - 79 , shape23 -- :: Polygon - 80 , concave -- :: Polygon - 81 -- * Internals - 82 , pRing -- :: APolygon a -> Ring a - 83 , pUnsafeMap -- :: (Ring a -> Ring a) -> APolygon a -> APolygon a - 84 , pCopy -- :: Polygon -> Polygon - 85 , pGenerate -- :: [(Double, Double)] -> Polygon - 86 , pUnGenerate -- :: Polygon -> [(Double, Double)] - 87 , Epsilon - 88 ) where - 89 - 90 -- import Control.Exception - 91 import Data.Hashable - 92 import Data.List (intersect, maximumBy, sort, sortOn, - 93 tails) - 94 import Data.Maybe - 95 import Data.Ratio - 96 import Data.Serialize - 97 import Data.Vector (Vector) - 98 import qualified Data.Vector as V - 99 import Linear.V2 - 100 import Linear.Vector - 101 import Reanimate.Math.Common - 102 -- import Reanimate.Math.EarClip - 103 import Reanimate.Math.SSSP - 104 import Reanimate.Math.Triangulate + 52 -- * Built-in shapes for testing + 53 , triangle -- :: Polygon + 54 , triangle' -- :: [P] + 55 , shape1 -- :: Polygon + 56 , shape2 -- :: Polygon + 57 , shape3 -- :: Polygon + 58 , shape4 -- :: Polygon + 59 , shape5 -- :: Polygon + 60 , shape6 -- :: Polygon + 61 , shape7 -- :: Polygon + 62 , shape8 -- :: Polygon + 63 , shape9 -- :: Polygon + 64 , shape10 -- :: Polygon + 65 , shape11 -- :: Polygon + 66 , shape12 -- :: Polygon + 67 , shape13 -- :: Polygon + 68 , shape14 -- :: Polygon + 69 , shape15 -- :: Polygon + 70 , shape16 -- :: Polygon + 71 , shape17 -- :: Polygon + 72 , shape18 -- :: Polygon + 73 , shape19 -- :: Polygon + 74 , shape20 -- :: Polygon + 75 , shape21 -- :: Polygon + 76 , shape22 -- :: Polygon + 77 , shape23 -- :: Polygon + 78 , concave -- :: Polygon + 79 -- * Internals + 80 , pRing -- :: APolygon a -> Ring a + 81 , pUnsafeMap -- :: (Ring a -> Ring a) -> APolygon a -> APolygon a + 82 , pCopy -- :: Polygon -> Polygon + 83 , pGenerate -- :: [(Double, Double)] -> Polygon + 84 , pUnGenerate -- :: Polygon -> [(Double, Double)] + 85 , Epsilon + 86 ) where + 87 + 88 -- import Control.Exception + 89 import Data.Hashable + 90 import Data.List (intersect, maximumBy, sort, sortOn, + 91 tails) + 92 import Data.Maybe + 93 import Data.Ratio + 94 import Data.Serialize + 95 import Data.Vector (Vector) + 96 import qualified Data.Vector as V + 97 import Linear.V2 + 98 import Linear.Vector + 99 import Reanimate.Math.Common + 100 -- import Reanimate.Math.EarClip + 101 import Reanimate.Math.SSSP + 102 import Reanimate.Math.Triangulate + 103 + 104 -- import Debug.Trace 105 - 106 -- import Debug.Trace - 107 - 108 -- Generate random polygons, options: - 109 -- 1. put corners around a circle. Vary the radius. - 110 -- 2. close a hilbert curve - 111 type FPolygon = APolygon Double - 112 -- Optimize representation? - 113 -- Polygon = (Vector XNumerator, Vector XDenominator - 114 -- ,Vector YNumerator, Vector YDenominator) - 115 data APolygon a = Polygon - 116 { polygonPoints :: Vector (V2 a) - 117 , polygonOffset :: Int - 118 , polygonTriangulation :: Triangulation - 119 , polygonSSSP :: Vector SSSP - 120 } - 121 type Polygon = APolygon Rational - 122 type P = V2 Double - 123 - 124 instance Show a => Show (APolygon a) where - 125 show = show . V.toList . polygonPoints - 126 - 127 instance Hashable a => Hashable (APolygon a) where - 128 hashWithSalt s p = V.foldl' hashWithSalt s (polygonPoints p) - 129 - 130 instance (PolyCtx a, Serialize a) => Serialize (APolygon a) where - 131 put = put . V.toList . polygonPoints - 132 get = mkPolygon . V.fromList <$> get - 133 - 134 pRing :: APolygon a -> Ring a - 135 pRing = ringPack . polygonPoints + 106 -- Generate random polygons, options: + 107 -- 1. put corners around a circle. Vary the radius. + 108 -- 2. close a hilbert curve + 109 type FPolygon = APolygon Double + 110 -- Optimize representation? + 111 -- Polygon = (Vector XNumerator, Vector XDenominator + 112 -- ,Vector YNumerator, Vector YDenominator) + 113 data APolygon a = Polygon + 114 { polygonPoints :: Vector (V2 a) + 115 , polygonOffset :: Int + 116 , polygonTriangulation :: Triangulation + 117 , polygonSSSP :: Vector SSSP + 118 } + 119 type Polygon = APolygon Rational + 120 type P = V2 Double + 121 + 122 instance Show a => Show (APolygon a) where + 123 show = show . V.toList . polygonPoints + 124 + 125 instance Hashable a => Hashable (APolygon a) where + 126 hashWithSalt s p = V.foldl' hashWithSalt s (polygonPoints p) + 127 + 128 instance (PolyCtx a, Serialize a) => Serialize (APolygon a) where + 129 put = put . V.toList . polygonPoints + 130 get = mkPolygon . V.fromList <$> get + 131 + 132 pRing :: APolygon a -> Ring a + 133 pRing = ringPack . polygonPoints + 134 + 135 type PolyCtx a = (Real a, Fractional a, Epsilon a) 136 - 137 type PolyCtx a = (Real a, Fractional a, Epsilon a) - 138 - 139 mkPolygon :: PolyCtx a => V.Vector (V2 a) -> APolygon a - 140 mkPolygon points = Polygon - 141 { polygonPoints = points - 142 , polygonOffset = 0 - 143 , polygonTriangulation = trig - 144 , polygonSSSP = V.generate n $ \i -> sssp ring (dual i trig) - 145 } - 146 where - 147 n = length points - 148 ring = ringPack points - 149 trig = triangulate ring - 150 -- earClip ring - 151 - 152 castPolygon :: (PolyCtx a, PolyCtx b) => APolygon a -> APolygon b - 153 castPolygon = mkPolygon . V.map (fmap realToFrac) . polygonPoints - 154 - 155 mkPolygonFromRing :: PolyCtx a => Ring a -> APolygon a - 156 mkPolygonFromRing = mkPolygon . ringUnpack - 157 - 158 pUnsafeMap :: (Ring a -> Ring a) -> APolygon a -> APolygon a - 159 pUnsafeMap fn p = p{ polygonPoints = ringUnpack (fn (pRing p)) } - 160 - 161 -- pParent p i j = shortest-path parent from j to i - 162 pParent :: APolygon a -> Int -> Int -> Int - 163 pParent p i j = - 164 (sTree V.! mod (j + polygonOffset p) n - polygonOffset p) `mod` n - 165 where - 166 sTree = polygonSSSP p V.! mod (i + polygonOffset p) n - 167 n = pSize p - 168 - 169 pCopy :: Polygon -> Polygon - 170 pCopy p = mkPolygon $ V.generate (pSize p) $ pAccess p - 171 - 172 pSetOffset :: APolygon a -> Int -> APolygon a - 173 pSetOffset p offset = - 174 p { polygonOffset = offset `mod` pSize p } - 175 - 176 pAdjustOffset :: APolygon a -> Int -> APolygon a - 177 pAdjustOffset p offset = - 178 p { polygonOffset = (polygonOffset p + offset) `mod` pSize p } - 179 - 180 {-# INLINE pSize #-} - 181 pSize :: APolygon a -> Int - 182 pSize = length . polygonPoints - 183 - 184 pNull :: APolygon a -> Bool - 185 pNull = V.null . polygonPoints - 186 - 187 pNext :: APolygon a -> Int -> Int - 188 pNext p i = (i+1) `mod` pSize p - 189 - 190 pPrev :: APolygon a -> Int -> Int - 191 pPrev p i = (i-1) `mod` pSize p - 192 - 193 -- When is a polygon valid/simple? - 194 -- It is counter-clockwise. - 195 -- No edges intersect. - 196 -- O(n^2) - 197 -- 'checkEdge' takes 90% of the time. - 198 pIsSimple :: Polygon -> Bool - 199 pIsSimple p | pSize p < 3 = False - 200 pIsSimple p = pIsCCW p && noDups && checkEdge 0 2 - 201 where - 202 noDups = checkForDups (sort (V.toList (polygonPoints p))) - 203 checkForDups (x:y:xs) - 204 = x /= y && checkForDups (y:xs) - 205 checkForDups _ = True - 206 len = pSize p - 207 -- check i,i+1 against j,j+1 - 208 -- j > i+1 - 209 checkEdge i j - 210 | j >= len = (i > len-3) || checkEdge (i+1) (i+3) - 211 | otherwise = - 212 case lineIntersect (pAccess p i, pAccess p $ i+1) - 213 (pAccess p j, pAccess p $ j+1) of - 214 Just u | u /= pAccess p i -> False - 215 _nothing -> checkEdge i (j+1) - 216 - 217 pScale :: Rational -> Polygon -> Polygon - 218 pScale s = pUnsafeMap (ringMap (^* s)) - 219 - 220 pAtCentroid :: Polygon -> Polygon - 221 pAtCentroid p = pTranslate (negate c) p - 222 where c = pCentroid p ^/ 2 - 223 - 224 pAtCenter :: Polygon -> Polygon - 225 pAtCenter p = pTranslate (negate $ pCenter p) p - 226 - 227 pTranslate :: V2 Rational -> Polygon -> Polygon - 228 pTranslate v = pUnsafeMap (ringMap (+v)) - 229 - 230 pCenter :: Polygon -> V2 Rational - 231 pCenter p = V2 (x+w/2) (y+h/2) - 232 where - 233 (x,y,w,h) = pBoundingBox p - 234 - 235 -- Returns (min-x, min-y, width, height) - 236 pBoundingBox :: Polygon -> (Rational, Rational, Rational, Rational) - 237 pBoundingBox = \p -> - 238 let V2 x y = pAccess p 0 in - 239 case V.foldl' worker (x, y, 0, 0) (polygonPoints p) of - 240 (xMin, yMin, xMax, yMax) -> - 241 (xMin, yMin, xMax-xMin, yMax-yMin) - 242 where - 243 worker (xMin,yMin,xMax,yMax) (V2 thisX thisY) = - 244 (min xMin thisX, min yMin thisY - 245 ,max xMax thisX, max yMax thisY) - 246 - 247 -- Place n points on a circle, use one parameter to slide the points back and forth. - 248 -- Use second parameter to move points closer to center circle. - 249 pGenerate :: [(Double, Double)] -> Polygon - 250 pGenerate points - 251 | len < 4 = error "pGenerate: require at least four points" - 252 | otherwise = mkPolygon $ V.fromList - 253 [ V2 (realToFrac $ cos ang * rMod) - 254 (realToFrac $ sin ang * rMod) - 255 | (i,(angMod,rMod)) <- zip [0..] points - 256 , let minAngle = tau / len * i - pi - 257 maxAngle = tau / len * (i+1) - pi - 258 ang = minAngle + (maxAngle-minAngle)*angMod - 259 ] - 260 where - 261 tau = 2*pi - 262 len = fromIntegral (length points) - 263 - 264 pUnGenerate :: Polygon -> [(Double, Double)] - 265 pUnGenerate p = - 266 [ worker i (fmap realToFrac e) - 267 | (i,e) <- zip [0..] (V.toList $ polygonPoints p) ] - 268 where - 269 len = fromIntegral (pSize p) - 270 worker i (V2 x y) = - 271 let ang = atan2 y x - 272 minAngle = tau / len * i - pi - 273 maxAngle = tau / len * (i+1) - pi - 274 in ((ang-minAngle)/(maxAngle-minAngle), sqrt (x*x+y*y)) - 275 tau = 2*pi - 276 - 277 -- When is a triangulation valid? - 278 -- Intersection: No internal edges intersect. - 279 -- Completeness: All edge neighbours share a single internal edge. - 280 isValidTriangulation :: Polygon -> Triangulation -> Bool - 281 isValidTriangulation p t = isComplete && intersectionFree - 282 where - 283 o = polygonOffset p - 284 isComplete = all isProper [0 .. pSize p-1] - 285 isProper i = - 286 let j = pNext p i in - 287 length ((pPrev p i : (t V.! i)) `intersect` (pNext p j : t V.! j)) == 1 - 288 intersectionFree = and - 289 [ case lineIntersect (pAccess p (a-o), pAccess p (b-o)) (pAccess p (c-o), pAccess p (d-o)) of - 290 Nothing -> True - 291 Just u -> u == pAccess p (a-o) || u == pAccess p (b-o) || - 292 u == pAccess p (c-o) || u == pAccess p (d-o) - 293 | ((a,b),(c,d)) <- edgePairs ] - 294 edgePairs = [ (e1, e2) | (e1, rest) <- zip edges (drop 1 $ tails edges), e2 <- rest] - 295 edges = - 296 [ (n, i) - 297 | (n, lst) <- zip [0..] (V.toList t) - 298 , i <- lst - 299 , n < i - 300 ] - 301 - 302 triangulationsToPolygons :: Polygon -> Triangulation -> [Polygon] - 303 triangulationsToPolygons p t = - 304 [ mkPolygon $ V.fromList - 305 [ pAccess p g, pAccess p i, pAccess p j ] - 306 | i <- [0 .. pSize p-1] - 307 , let js = filter (i<) $ t V.! i - 308 , (g, j) <- zip (i-1:js) js - 309 ] - 310 - 311 pIsInside :: Polygon -> V2 Rational -> Bool - 312 pIsInside p point = or - 313 [ isInside (rawAccess g) (rawAccess i) (rawAccess j) point - 314 | i <- [0 .. pSize p-1] - 315 , let js = filter (i<) $ polygonTriangulation p V.! i - 316 , (g, j) <- zip (i-1:js) js - 317 ] - 318 where - 319 rawAccess x = polygonPoints p V.! x - 320 - 321 -- reducePolygons :: Int -> [Polygon] -> [Polygon] - 322 -- reducePolygons n ps - 323 -- | length ps <= n = ps - 324 -- | otherwise = - 325 -- let p = findSmallest ps - 326 -- es = edges p - 327 -- e = findSmallest es - 328 -- in reducePolygons n (merge p e : delete p (delete e ps)) - 329 -- where - 330 -- findSmallest = minimumBy (comparing area2X) - 331 -- shareEdge p1 p2 = - 332 - 333 {-# INLINE pAccess #-} - 334 pAccess :: APolygon a -> Int -> V2 a - 335 pAccess p i = -- polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p) - 336 polygonPoints p `V.unsafeIndex` ((polygonOffset p + i) `mod` pSize p) - 337 - 338 triangle :: Polygon - 339 triangle = mkPolygon $ V.fromList [V2 1 1, V2 0 0, V2 2 0] - 340 - 341 triangle' :: [P] - 342 triangle' = reverse [V2 1 1, V2 0 0, V2 2 0] - 343 - 344 shape1 :: Polygon - 345 shape1 = mkPolygon $ V.fromList - 346 [ V2 0 0, V2 2 0 - 347 , V2 2 1, V2 2 2, V2 2 3, V2 2 4, V2 2 5, V2 2 6 - 348 , V2 1 1, V2 0 1 ] - 349 - 350 shape2 :: Polygon - 351 shape2 = mkPolygon $ V.fromList - 352 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 0 (-1), V2 0 (-2) - 353 , V2 3 (-2), V2 3 2, V2 0 2] - 354 - 355 shape3 :: Polygon - 356 shape3 = mkPolygon $ V.fromList - 357 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 2, V2 0 2] - 358 - 359 shape4 :: Polygon - 360 shape4 = mkPolygon $ V.fromList - 361 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 3 (-1),V2 3 2, V2 0 2] - 362 - 363 shape5 :: Polygon - 364 shape5 = pCycles shape4 !! 2 - 365 - 366 -- square - 367 shape6 :: Polygon - 368 shape6 = mkPolygon $ V.fromList [ V2 0 0, V2 1 0, V2 1 1, V2 0 1 ] - 369 - 370 shape7 :: Polygon - 371 shape7 = pScale 6 $ mkPolygon $ V.fromList - 372 [V2 ((-1567171105775771) % 144115188075855872) ((-7758063241391039) % 1152921504606846976) - 373 ,V2 ((-2711114907999263) % 18014398509481984) ((-3561889280168807) % 18014398509481984) - 374 ,V2 ((-6897139157863177) % 72057594037927936) ((-1632144794297397) % 4503599627370496) - 375 ,V2 (5592137945106423 % 36028797018963968) ((-71351641856107) % 281474976710656) - 376 ,V2 (2568147525079071 % 4503599627370496) ((-4312925637247687) % 18014398509481984) - 377 ,V2 (1291079014395023 % 2251799813685248) (321513444515769 % 2251799813685248) - 378 ,V2 (2071709221627247 % 4503599627370496) (4019115966736491 % 9007199254740992) - 379 ,V2 ((-1589087869859839) % 144115188075855872) (4904023654354179 % 9007199254740992) - 380 ,V2 ((-2328090886101149) % 36028797018963968) (2587887893460759 % 36028797018963968) - 381 ,V2 ((-7990199074159871) % 18014398509481984) (1301850651537745 % 4503599627370496)] - 382 - 383 shape8 :: Polygon - 384 shape8 = pScale 10 $ pGenerate - 385 [(0.36,0.4),(0.7,1.8e-2),(0.7,0.2),(0.1,0.4),(0.2,0.2),(0.7,0.1),(0.4,8.0e-2)] - 386 - 387 shape9 :: Polygon - 388 shape9 = pScale 5 $ pGenerate - 389 [(0.5,0.2),(0.7,0.6),(0.4,0.3),(0.1,0.7),(0.3,1.0e-2),(0.5,0.3),(0.2,0.8),(0.1,0.8),(0.7,6.0e-2),(0.1,0.6)] - 390 - 391 shape10 :: Polygon - 392 shape10 = pGenerate - 393 [(0.4,0.7),(0.2,0.2),(0.3,0.9),(5.0e-2,0.1),(0.7,1.0e-2),(0.7,0.9),(0.2,0.1),(0.5,6.0e-2),(0.6,9.0e-2)] - 394 - 395 shape11 :: Polygon - 396 shape11 = pGenerate - 397 [(0.1,0.8),(0.7,0.6),(0.7,0.4),(0.3,0.5),(0.8,0.9),(0.8,6.0e-2),(1.0e-2,4.0e-2),(0.8,0.1)] - 398 - 399 shape12 :: Polygon - 400 shape12 = mkPolygon $ V.fromList - 401 [ V2 0 0, V2 0.5 1.5, V2 2 2, V2 (-2) 2, V2 (-0.5) 1.5 ] - 402 - 403 -- F shape - 404 shape13 :: Polygon - 405 shape13 = pCycles (mkPolygon $ V.reverse (V.fromList - 406 [ V2 0 0, V2 0 2 - 407 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 - 408 , V2 1 1, V2 1 0.7 - 409 , V2 0.3 0.7, V2 0.3 0 ])) !! 7 - 410 - 411 -- E shape - 412 shape14 :: Polygon - 413 shape14 = pCycles (mkPolygon $ V.reverse $ V.fromList - 414 [ V2 0 0, V2 0 2 -- up - 415 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 -- first prong - 416 , V2 1 1, V2 1 0.7, V2 0.3 0.7, V2 0.3 0.3 -- second prong - 417 , V2 1 0.3, V2 1 0 -- last prong - 418 ]) !! 9 - 419 - 420 -- - 421 shape15 :: Polygon - 422 shape15 = mkPolygon $ V.fromList - 423 [ V2 0 0, V2 2 0 - 424 , V2 2 2, V2 1 2 - 425 , V2 1 1, V2 0 1] - 426 - 427 shape16 :: Polygon - 428 shape16 = mkPolygon $ V.fromList - 429 [ V2 0 0, V2 2 0 - 430 , V2 2 1, V2 1 1 - 431 , V2 1 2, V2 0 2] - 432 - 433 shape17 :: Polygon - 434 shape17 = mkPolygon $ V.fromList - 435 [ V2 2 0, V2 2 1 - 436 , V2 1 1, V2 1 2 - 437 , V2 0 2, V2 0 1, V2 0 0 ] - 438 - 439 shape18 :: Polygon - 440 shape18 = mkPolygon $ V.fromList - 441 [ V2 2 0, V2 2 1, V2 2 2 - 442 , V2 1 2, V2 1 1 - 443 , V2 0 1, V2 0 0 ] - 444 - 445 shape19 :: Polygon - 446 shape19 = mkPolygon $ V.fromList - 447 [ V2 (-3) (-3), V2 0 (-1) - 448 , V2 3 (-3), V2 1 0 - 449 , V2 3 3, V2 0 1 - 450 , V2 (-3) 3, V2 (-1) 0 ] - 451 - 452 shape20 :: Polygon - 453 shape20 = mkPolygon $ V.fromList - 454 [ V2 (-3) (-3) - 455 , V2 0 (-1) - 456 , V2 3 (-3) - 457 , V2 5 0 - 458 , V2 2.5 (-2) - 459 , V2 1 0 - 460 , V2 3 3 - 461 , V2 0 1 - 462 , V2 (-3) 3 - 463 , V2 (-1) 0 ] - 464 - 465 shape21 :: Polygon - 466 shape21 = mkPolygon $ V.fromList - 467 [V2 0.0 0.0,V2 1.0 0.0,V2 1.0 1.0,V2 2.0 1.0,V2 2.0 (-1.0),V2 3.0 (-1.0) - 468 ,V2 3.0 2.0,V2 0.0 2.0] - 469 - 470 shape22 :: Polygon - 471 shape22 = pScale 2 $ mkPolygon $ V.fromList - 472 [V2 (-0.17) (-0.08) - 473 ,V2 (-0.34) (-0.21) - 474 ,V2 0.0 0.0 - 475 ,V2 (-0.10) 0.60 - 476 ,V2 (-0.14) 0.19 - 477 ,V2 (-0.05) 0.03 - 478 ] - 479 - 480 shape23 :: Polygon - 481 shape23 = mkPolygon $ V.fromList - 482 [ V2 0 0, V2 4 0 - 483 , V2 4 3, V2 2 3 - 484 , V2 2 2, V2 3 2 - 485 , V2 3 1, V2 1 1 - 486 , V2 1 2, V2 2 2 - 487 , V2 2 3, V2 0 3 ] - 488 - 489 concave :: Polygon - 490 concave = mkPolygon $ - 491 V.fromList [V2 0 0, V2 2 0, V2 2 2, V2 1 1, V2 0 2] - 492 - 493 pMkWinding :: Int -> Polygon - 494 pMkWinding n | n < 1 = error "Polygon must have at least one winding." - 495 pMkWinding n = mkPolygon $ - 496 V.fromList $ p0 : p1 : walkTo p1 1 n (V2 1 0) ++ reverse (walkTo p0 1 (n+2) (V2 (-1) 0)) - 497 where - 498 p0 = V2 0 0 - 499 p1 = V2 0 1 - 500 walkTo at a b dir - 501 | a == b = [] - 502 | otherwise = - 503 let newAt = at + (dir ^* toRational a) - 504 in newAt : walkTo newAt (a+1) b (rot dir) - 505 rot (V2 x y) = - 506 V2 y (-x) - 507 - 508 pDeoverlap :: Polygon -> Polygon - 509 pDeoverlap p = mkPolygon arr - 510 where - 511 arr = V.generate (pSize p) worker - 512 worker 0 = pAccess p 0 - 513 worker n = - 514 if length (V.elemIndices (pAccess p n) (polygonPoints p)) /= 1 - 515 then - 516 let prev = arr V.! (n-1) - 517 this = pAccess p n - 518 in lerp 0.99999 this prev - 519 else pAccess p n - 520 - 521 pCycles :: APolygon a -> [APolygon a] - 522 pCycles p = map (pAdjustOffset p) [0 .. pSize p-1] - 523 - 524 pCycle :: PolyCtx a => APolygon a -> Double -> APolygon a - 525 pCycle p 0 = p - 526 pCycle p t = mkPolygon $ worker 0 0 - 527 where - 528 worker acc i - 529 | segment + acc > limit = - 530 V.singleton (lerp (realToFrac $ (segment + acc - limit)/segment) x y) <> - 531 -- V.drop (i+1) (polygonPoints p) <> - 532 V.fromList (map (pAccess p) [i+1..pSize p-1]) <> - 533 V.fromList (map (pAccess p) [0 .. i]) - 534 -- V.take (i+1) (polygonPoints p) - 535 | i == pSize p-1 = V.fromList (map (pAccess p) [0 .. pSize p-1]) - 536 | otherwise = worker (acc+segment) (i+1) - 537 where - 538 x = pAccess p i - 539 y = pAccess p $ i+1 - 540 segment = distance' x y - 541 len = pCircumference' p - 542 limit = t * len - 543 - 544 pCentroid :: Fractional a => APolygon a -> V2 a - 545 pCentroid p = V2 cx cy - 546 where - 547 a = pArea p - 548 cx = recip (6*a) * V.sum (pMapEdges fnX p) - 549 cy = recip (6*a) * V.sum (pMapEdges fnY p) - 550 fnX (V2 x y) (V2 x' y') = (x+x')*(x*y' - x'*y) - 551 fnY (V2 x y) (V2 x' y') = (y+y')*(x*y' - x'*y) - 552 - 553 {-# INLINE pMapEdges #-} - 554 pMapEdges :: (V2 a -> V2 a -> b) -> APolygon a -> V.Vector b - 555 pMapEdges fn p = V.generate n $ \i -> - 556 if i == n-1 - 557 then fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) - 558 else fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1)) - 559 where - 560 n = pSize p - 561 arr = polygonPoints p - 562 - 563 {-# SPECIALIZE pArea :: APolygon Double -> Double #-} - 564 {-# SPECIALIZE pArea :: APolygon Rational -> Rational #-} - 565 pArea :: (Fractional a) => APolygon a -> a - 566 pArea p = - 567 -- 0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p) - 568 0.5 * worker 0 0 - 569 where - 570 fn (V2 x y) (V2 x' y') = x*y' - x'*y - 571 arr = polygonPoints p - 572 worker !acc i - 573 | i == pSize p - 1 = acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) - 574 | otherwise = - 575 worker (acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1))) (i+1) - 576 - 577 pCircumference :: (Real a, Fractional a) => APolygon a -> a - 578 pCircumference p = sum - 579 [ approxDist (pAccess p i) (pAccess p $ i+1) - 580 | i <- [0 .. pSize p-1]] - 581 - 582 pCircumference' :: (Real a, Fractional a) => APolygon a -> Double - 583 pCircumference' p = sum - 584 [ distance' (pAccess p i) (pAccess p $ i+1) - 585 | i <- [0 .. pSize p-1]] - 586 - 587 - 588 -- Add points by splitting the longest lines in half repeatedly. - 589 pAddPoints :: PolyCtx a => Int -> APolygon a -> APolygon a - 590 pAddPoints n p | n <= 0 = p - 591 pAddPoints n p = pAddPoints (n-1) $ - 592 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] - 593 where - 594 worker idx - 595 | idx == longestEdge = - 596 let start = pAccess p idx - 597 end = pAccess p $ idx+1 - 598 middle = lerp 0.5 end start - 599 in [start, middle] - 600 | otherwise = [pAccess p idx] - 601 longestEdge = maximumBy cmpLength [0 .. pSize p-1] - 602 cmpLength a b = - 603 distSquared (pAccess p a) (pAccess p $ a+1) `compare` - 604 distSquared (pAccess p b) (pAccess p $ b+1) - 605 - 606 pAddPointsRestricted :: PolyCtx a => [(V2 a, V2 a)] -> Int -> APolygon a -> APolygon a - 607 pAddPointsRestricted _immutableEdges n p | n <= 0 = p - 608 pAddPointsRestricted immutableEdges n p = pAddPointsRestricted immutableEdges (n-1) $ - 609 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] - 610 where - 611 isImmutable idx = - 612 (pAccess p idx, pAccess p $ idx+1) `elem` immutableEdges || - 613 (pAccess p $ idx+1, pAccess p idx) `elem` immutableEdges - 614 worker idx - 615 | idx == longestEdge && not (isImmutable idx) = - 616 let start = pAccess p idx - 617 end = pAccess p $ idx+1 - 618 middle = lerp 0.5 end start - 619 in [start, middle] - 620 | otherwise = [pAccess p idx] - 621 longestEdge = maximumBy cmpLength [0 .. pSize p-1] - 622 cmpLength a _ | isImmutable a = LT - 623 cmpLength _ b | isImmutable b = GT - 624 cmpLength a b = - 625 distSquared (pAccess p a) (pAccess p $ a+1) `compare` - 626 distSquared (pAccess p b) (pAccess p $ b+1) - 627 - 628 pAddPointsBetween :: PolyCtx a => (Int, Int) -> Int -> APolygon a -> APolygon a - 629 pAddPointsBetween _ n p | n <= 0 = p - 630 pAddPointsBetween (i,l) n p = pAddPointsBetween (i,l+1) (n-1) $ - 631 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] - 632 where - 633 worker idx - 634 | idx == longestEdge = - 635 let start = pAccess p idx - 636 end = pAccess p $ idx+1 - 637 middle = lerp 0.5 end start - 638 in [start, middle] - 639 | otherwise = [pAccess p idx] - 640 longestEdge = maximumBy cmpLength [i .. i+l-1] - 641 cmpLength a b = - 642 distSquared (pAccess p a) (pAccess p $ a+1) `compare` - 643 distSquared (pAccess p b) (pAccess p $ b+1) - 644 - 645 -- addPoints :: Int -> Polygon -> Polygon - 646 -- addPoints n p = mkPolygon $ V.fromList $ worker n 0 (map (pAccess p) [0..s]) - 647 -- where - 648 -- worker 0 _ rest = init rest - 649 -- worker i acc (x:y:xs) = - 650 -- let xy = approxDist x y in - 651 -- if acc + xy > limit - 652 -- then x : worker (i-1) 0 (lerp ((limit-acc)/xy) y x : y:xs) - 653 -- else x : worker i (acc+xy) (y:xs) - 654 -- worker _ _ [_] = [] - 655 -- worker _ _ _ = error "addPoints: invalid polygon" - 656 -- s = pSize p - 657 -- len = polygonLength p - 658 -- limit = len / fromIntegral (n+1) - 659 - 660 pIsConvex :: Polygon -> Bool - 661 pIsConvex p = and - 662 [ area2X (pAccess p i) (pAccess p j) (pAccess p k) > 0 - 663 | i <- [0..n-1] - 664 , j <- [i+1..n-1] - 665 , k <- [j+1..n-1] - 666 ] - 667 where n = pSize p - 668 - 669 pIsCCW :: Polygon -> Bool - 670 pIsCCW p | pNull p = False - 671 pIsCCW p = V.sum (pMapEdges fn p) < 0 - 672 where - 673 fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1) - 674 - 675 {-# INLINE pRayIntersect #-} - 676 pRayIntersect :: PolyCtx a => APolygon a -> (Int, Int) -> (Int,Int) -> Maybe (V2 a) - 677 pRayIntersect p (a,b) (c,d) = - 678 rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d) - 679 - 680 pCuts :: (Real a, Fractional a, Epsilon a) => APolygon a -> [(APolygon a,APolygon a)] - 681 pCuts p = - 682 [ pCutAt (pAdjustOffset p i) (j-i) - 683 | i <- [0 .. pSize p-1 ] - 684 , j <- [i+2 .. pSize p-1 ] - 685 , (j+1) `mod` pSize p /= i - 686 , pParent p i j == i ] - 687 - 688 pCutEqual :: PolyCtx a => APolygon a -> (APolygon a, APolygon a) - 689 pCutEqual p = - 690 fromMaybe (p,p) $ listToMaybe $ sortOn f $ pCuts p - 691 where - 692 f (a,b) = abs (pArea a - pArea b) - 693 - 694 -- FIXME: This should be more efficient - 695 pCutAt :: PolyCtx a => APolygon a -> Int -> (APolygon a, APolygon a) - 696 pCutAt p i = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right) - 697 where - 698 n = pSize p - 699 left = map (pAccess p) [0 .. i] - 700 right = map (pAccess p) (0:[i..n-1]) - 701 - 702 pOverlap :: PolyCtx a => APolygon a -> APolygon a -> APolygon a - 703 pOverlap a b = mkPolygon $ V.fromList $ clearDups $ concatMap edgeIntersect [0 .. pSize a-1] - 704 where - 705 clearDups (x:y:xs) - 706 | x == y = clearDups (y:xs) - 707 | otherwise = x : clearDups (y:xs) - 708 clearDups xs = xs - 709 edgeIntersect edge = - 710 sortOn (distSquared (pAccess a edge)) $ catMaybes - 711 [ lineIntersect (aP, aP') (bP, bP') - 712 | i <- [0 .. pSize b-1] - 713 , let aP = pAccess a edge - 714 aP' = pAccess a (edge+1) - 715 bP = pAccess b i - 716 bP' = pAccess b (i+1) - 717 ] - 718 - 719 --------------------------------------------------------- - 720 -- SSSP visibility and SSSP windows - 721 - 722 ssspVisibility :: PolyCtx a => APolygon a -> APolygon a - 723 ssspVisibility p = mkPolygon $ - 724 V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1]) - 725 where - 726 clearDups (x:y:xs) - 727 | x == y = clearDups (y:xs) - 728 | otherwise = x : clearDups (y:xs) - 729 clearDups xs = xs - 730 obstructedBy n = - 731 case pParent p 0 n of - 732 0 -> n - 733 i -> obstructedBy i - 734 go [] = [] - 735 go [x] = [pAccess p x] - 736 go (x:y:xs) = - 737 let xO = obstructedBy x - 738 yO = obstructedBy y - 739 in case () of - 740 () - 741 -- Both ends are visible. - 742 | xO == x && yO == y -> pAccess p x : go (y:xs) - 743 -- X is visible, x to intersect (0,yO) (x,y) - 744 | xO == x -> - 745 pAccess p x : fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y)) : go (y:xs) - 746 -- Y is visible - 747 | yO == y -> fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) : pAccess p y : go (y:xs) - 748 -- Neither is visible and they've obstructed by the same point - 749 -- so the entire edge is hidden. - 750 | xO == yO -> go (y:xs) - 751 -- Neither is visible. Cast shadow from obstruction points to - 752 -- find if a subsection of the edge is visible. - 753 | otherwise -> - 754 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y)) - 755 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y)) - 756 in if a /= b - 757 then a : b : go (y:xs) - 758 else go (y:xs) - 759 - 760 ssspWindows :: Polygon -> [(V2 Rational, V2 Rational)] - 761 ssspWindows p = clearDups $ go (pAccess p 0) [0..pSize p-1] - 762 where - 763 clearDups (x:y:xs) - 764 | x == y = clearDups (y:xs) - 765 | otherwise = x : clearDups (y:xs) - 766 clearDups xs = xs - 767 obstructedBy n = - 768 case pParent p 0 n of - 769 0 -> n - 770 i -> obstructedBy i - 771 go _ [] = [] - 772 go _ [_] = [] - 773 go l (x:y:xs) = - 774 let xO = obstructedBy x - 775 yO = obstructedBy y - 776 in case () of - 777 () - 778 -- Both ends are visible. - 779 | xO == x && yO == y -> go (pAccess p x) (y:xs) - 780 -- X is visible, x to intersect (0,yO) (x,y) - 781 | xO == x -> - 782 go (fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y))) (y:xs) - 783 -- Y is visible - 784 | yO == y -> - 785 let newL = fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) in - 786 (l, newL) : - 787 go newL (y:xs) - 788 -- Neither is visible and they've obstructed by the same point - 789 -- so the entire edge is hidden. - 790 | xO == yO -> go l (y:xs) - 791 -- Neither is visible. Cast shadow from obstruction points to - 792 -- find if a subsection of the edge is visible. - 793 | otherwise -> - 794 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y)) - 795 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y)) - 796 in if a /= b - 797 then (l, a) : (b, pAccess p yO) : go (pAccess p yO) (y:xs) - 798 else go l (y:xs) - 799 - 800 pdualPolygons :: Polygon -> PDual -> [Polygon] - 801 pdualPolygons p pdual = map mkPolygonFromRing (pdualRings (pRing p) pdual) + 137 mkPolygon :: PolyCtx a => V.Vector (V2 a) -> APolygon a + 138 mkPolygon points = Polygon + 139 { polygonPoints = points + 140 , polygonOffset = 0 + 141 , polygonTriangulation = trig + 142 , polygonSSSP = V.generate n $ \i -> sssp ring (dual i trig) + 143 } + 144 where + 145 n = length points + 146 ring = ringPack points + 147 trig = triangulate ring + 148 -- earClip ring + 149 + 150 castPolygon :: (PolyCtx a, PolyCtx b) => APolygon a -> APolygon b + 151 castPolygon = mkPolygon . V.map (fmap realToFrac) . polygonPoints + 152 + 153 mkPolygonFromRing :: PolyCtx a => Ring a -> APolygon a + 154 mkPolygonFromRing = mkPolygon . ringUnpack + 155 + 156 pUnsafeMap :: (Ring a -> Ring a) -> APolygon a -> APolygon a + 157 pUnsafeMap fn p = p{ polygonPoints = ringUnpack (fn (pRing p)) } + 158 + 159 -- pParent p i j = shortest-path parent from j to i + 160 pParent :: APolygon a -> Int -> Int -> Int + 161 pParent p i j = + 162 (sTree V.! mod (j + polygonOffset p) n - polygonOffset p) `mod` n + 163 where + 164 sTree = polygonSSSP p V.! mod (i + polygonOffset p) n + 165 n = pSize p + 166 + 167 pCopy :: Polygon -> Polygon + 168 pCopy p = mkPolygon $ V.generate (pSize p) $ pAccess p + 169 + 170 pSetOffset :: APolygon a -> Int -> APolygon a + 171 pSetOffset p offset = + 172 p { polygonOffset = offset `mod` pSize p } + 173 + 174 pAdjustOffset :: APolygon a -> Int -> APolygon a + 175 pAdjustOffset p offset = + 176 p { polygonOffset = (polygonOffset p + offset) `mod` pSize p } + 177 + 178 {-# INLINE pSize #-} + 179 pSize :: APolygon a -> Int + 180 pSize = length . polygonPoints + 181 + 182 pNull :: APolygon a -> Bool + 183 pNull = V.null . polygonPoints + 184 + 185 pNext :: APolygon a -> Int -> Int + 186 pNext p i = (i+1) `mod` pSize p + 187 + 188 pPrev :: APolygon a -> Int -> Int + 189 pPrev p i = (i-1) `mod` pSize p + 190 + 191 -- When is a polygon valid/simple? + 192 -- It is counter-clockwise. + 193 -- No edges intersect. + 194 -- O(n^2) + 195 -- 'checkEdge' takes 90% of the time. + 196 pIsSimple :: Polygon -> Bool + 197 pIsSimple p | pSize p < 3 = False + 198 pIsSimple p = pIsCCW p && noDups && checkEdge 0 2 + 199 where + 200 noDups = checkForDups (sort (V.toList (polygonPoints p))) + 201 checkForDups (x:y:xs) + 202 = x /= y && checkForDups (y:xs) + 203 checkForDups _ = True + 204 len = pSize p + 205 -- check i,i+1 against j,j+1 + 206 -- j > i+1 + 207 checkEdge i j + 208 | j >= len = (i > len-3) || checkEdge (i+1) (i+3) + 209 | otherwise = + 210 case lineIntersect (pAccess p i, pAccess p $ i+1) + 211 (pAccess p j, pAccess p $ j+1) of + 212 Just u | u /= pAccess p i -> False + 213 _nothing -> checkEdge i (j+1) + 214 + 215 pScale :: Rational -> Polygon -> Polygon + 216 pScale s = pUnsafeMap (ringMap (^* s)) + 217 + 218 pAtCentroid :: Polygon -> Polygon + 219 pAtCentroid p = pTranslate (negate c) p + 220 where c = pCentroid p ^/ 2 + 221 + 222 pAtCenter :: Polygon -> Polygon + 223 pAtCenter p = pTranslate (negate $ pCenter p) p + 224 + 225 pTranslate :: V2 Rational -> Polygon -> Polygon + 226 pTranslate v = pUnsafeMap (ringMap (+v)) + 227 + 228 pCenter :: Polygon -> V2 Rational + 229 pCenter p = V2 (x+w/2) (y+h/2) + 230 where + 231 (x,y,w,h) = pBoundingBox p + 232 + 233 -- Returns (min-x, min-y, width, height) + 234 pBoundingBox :: Polygon -> (Rational, Rational, Rational, Rational) + 235 pBoundingBox = \p -> + 236 let V2 x y = pAccess p 0 in + 237 case V.foldl' worker (x, y, 0, 0) (polygonPoints p) of + 238 (xMin, yMin, xMax, yMax) -> + 239 (xMin, yMin, xMax-xMin, yMax-yMin) + 240 where + 241 worker (xMin,yMin,xMax,yMax) (V2 thisX thisY) = + 242 (min xMin thisX, min yMin thisY + 243 ,max xMax thisX, max yMax thisY) + 244 + 245 -- Place n points on a circle, use one parameter to slide the points back and forth. + 246 -- Use second parameter to move points closer to center circle. + 247 pGenerate :: [(Double, Double)] -> Polygon + 248 pGenerate points + 249 | len < 4 = error "pGenerate: require at least four points" + 250 | otherwise = mkPolygon $ V.fromList + 251 [ V2 (realToFrac $ cos ang * rMod) + 252 (realToFrac $ sin ang * rMod) + 253 | (i,(angMod,rMod)) <- zip [0..] points + 254 , let minAngle = tau / len * i - pi + 255 maxAngle = tau / len * (i+1) - pi + 256 ang = minAngle + (maxAngle-minAngle)*angMod + 257 ] + 258 where + 259 tau = 2*pi + 260 len = fromIntegral (length points) + 261 + 262 pUnGenerate :: Polygon -> [(Double, Double)] + 263 pUnGenerate p = + 264 [ worker i (fmap realToFrac e) + 265 | (i,e) <- zip [0..] (V.toList $ polygonPoints p) ] + 266 where + 267 len = fromIntegral (pSize p) + 268 worker i (V2 x y) = + 269 let ang = atan2 y x + 270 minAngle = tau / len * i - pi + 271 maxAngle = tau / len * (i+1) - pi + 272 in ((ang-minAngle)/(maxAngle-minAngle), sqrt (x*x+y*y)) + 273 tau = 2*pi + 274 + 275 -- When is a triangulation valid? + 276 -- Intersection: No internal edges intersect. + 277 -- Completeness: All edge neighbours share a single internal edge. + 278 isValidTriangulation :: Polygon -> Triangulation -> Bool + 279 isValidTriangulation p t = isComplete && intersectionFree + 280 where + 281 o = polygonOffset p + 282 isComplete = all isProper [0 .. pSize p-1] + 283 isProper i = + 284 let j = pNext p i in + 285 length ((pPrev p i : (t V.! i)) `intersect` (pNext p j : t V.! j)) == 1 + 286 intersectionFree = and + 287 [ case lineIntersect (pAccess p (a-o), pAccess p (b-o)) (pAccess p (c-o), pAccess p (d-o)) of + 288 Nothing -> True + 289 Just u -> u == pAccess p (a-o) || u == pAccess p (b-o) || + 290 u == pAccess p (c-o) || u == pAccess p (d-o) + 291 | ((a,b),(c,d)) <- edgePairs ] + 292 edgePairs = [ (e1, e2) | (e1, rest) <- zip edges (drop 1 $ tails edges), e2 <- rest] + 293 edges = + 294 [ (n, i) + 295 | (n, lst) <- zip [0..] (V.toList t) + 296 , i <- lst + 297 , n < i + 298 ] + 299 + 300 triangulationsToPolygons :: Polygon -> Triangulation -> [Polygon] + 301 triangulationsToPolygons p t = + 302 [ mkPolygon $ V.fromList + 303 [ pAccess p g, pAccess p i, pAccess p j ] + 304 | i <- [0 .. pSize p-1] + 305 , let js = filter (i<) $ t V.! i + 306 , (g, j) <- zip (i-1:js) js + 307 ] + 308 + 309 pIsInside :: Polygon -> V2 Rational -> Bool + 310 pIsInside p point = or + 311 [ isInside (rawAccess g) (rawAccess i) (rawAccess j) point + 312 | i <- [0 .. pSize p-1] + 313 , let js = filter (i<) $ polygonTriangulation p V.! i + 314 , (g, j) <- zip (i-1:js) js + 315 ] + 316 where + 317 rawAccess x = polygonPoints p V.! x + 318 + 319 -- reducePolygons :: Int -> [Polygon] -> [Polygon] + 320 -- reducePolygons n ps + 321 -- | length ps <= n = ps + 322 -- | otherwise = + 323 -- let p = findSmallest ps + 324 -- es = edges p + 325 -- e = findSmallest es + 326 -- in reducePolygons n (merge p e : delete p (delete e ps)) + 327 -- where + 328 -- findSmallest = minimumBy (comparing area2X) + 329 -- shareEdge p1 p2 = + 330 + 331 {-# INLINE pAccess #-} + 332 pAccess :: APolygon a -> Int -> V2 a + 333 pAccess p i = -- polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p) + 334 polygonPoints p `V.unsafeIndex` ((polygonOffset p + i) `mod` pSize p) + 335 + 336 triangle :: Polygon + 337 triangle = mkPolygon $ V.fromList [V2 1 1, V2 0 0, V2 2 0] + 338 + 339 triangle' :: [P] + 340 triangle' = reverse [V2 1 1, V2 0 0, V2 2 0] + 341 + 342 shape1 :: Polygon + 343 shape1 = mkPolygon $ V.fromList + 344 [ V2 0 0, V2 2 0 + 345 , V2 2 1, V2 2 2, V2 2 3, V2 2 4, V2 2 5, V2 2 6 + 346 , V2 1 1, V2 0 1 ] + 347 + 348 shape2 :: Polygon + 349 shape2 = mkPolygon $ V.fromList + 350 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 0 (-1), V2 0 (-2) + 351 , V2 3 (-2), V2 3 2, V2 0 2] + 352 + 353 shape3 :: Polygon + 354 shape3 = mkPolygon $ V.fromList + 355 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 2, V2 0 2] + 356 + 357 shape4 :: Polygon + 358 shape4 = mkPolygon $ V.fromList + 359 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 3 (-1),V2 3 2, V2 0 2] + 360 + 361 shape5 :: Polygon + 362 shape5 = pCycles shape4 !! 2 + 363 + 364 -- square + 365 shape6 :: Polygon + 366 shape6 = mkPolygon $ V.fromList [ V2 0 0, V2 1 0, V2 1 1, V2 0 1 ] + 367 + 368 shape7 :: Polygon + 369 shape7 = pScale 6 $ mkPolygon $ V.fromList + 370 [V2 ((-1567171105775771) % 144115188075855872) ((-7758063241391039) % 1152921504606846976) + 371 ,V2 ((-2711114907999263) % 18014398509481984) ((-3561889280168807) % 18014398509481984) + 372 ,V2 ((-6897139157863177) % 72057594037927936) ((-1632144794297397) % 4503599627370496) + 373 ,V2 (5592137945106423 % 36028797018963968) ((-71351641856107) % 281474976710656) + 374 ,V2 (2568147525079071 % 4503599627370496) ((-4312925637247687) % 18014398509481984) + 375 ,V2 (1291079014395023 % 2251799813685248) (321513444515769 % 2251799813685248) + 376 ,V2 (2071709221627247 % 4503599627370496) (4019115966736491 % 9007199254740992) + 377 ,V2 ((-1589087869859839) % 144115188075855872) (4904023654354179 % 9007199254740992) + 378 ,V2 ((-2328090886101149) % 36028797018963968) (2587887893460759 % 36028797018963968) + 379 ,V2 ((-7990199074159871) % 18014398509481984) (1301850651537745 % 4503599627370496)] + 380 + 381 shape8 :: Polygon + 382 shape8 = pScale 10 $ pGenerate + 383 [(0.36,0.4),(0.7,1.8e-2),(0.7,0.2),(0.1,0.4),(0.2,0.2),(0.7,0.1),(0.4,8.0e-2)] + 384 + 385 shape9 :: Polygon + 386 shape9 = pScale 5 $ pGenerate + 387 [(0.5,0.2),(0.7,0.6),(0.4,0.3),(0.1,0.7),(0.3,1.0e-2),(0.5,0.3),(0.2,0.8),(0.1,0.8),(0.7,6.0e-2),(0.1,0.6)] + 388 + 389 shape10 :: Polygon + 390 shape10 = pGenerate + 391 [(0.4,0.7),(0.2,0.2),(0.3,0.9),(5.0e-2,0.1),(0.7,1.0e-2),(0.7,0.9),(0.2,0.1),(0.5,6.0e-2),(0.6,9.0e-2)] + 392 + 393 shape11 :: Polygon + 394 shape11 = pGenerate + 395 [(0.1,0.8),(0.7,0.6),(0.7,0.4),(0.3,0.5),(0.8,0.9),(0.8,6.0e-2),(1.0e-2,4.0e-2),(0.8,0.1)] + 396 + 397 shape12 :: Polygon + 398 shape12 = mkPolygon $ V.fromList + 399 [ V2 0 0, V2 0.5 1.5, V2 2 2, V2 (-2) 2, V2 (-0.5) 1.5 ] + 400 + 401 -- F shape + 402 shape13 :: Polygon + 403 shape13 = pCycles (mkPolygon $ V.reverse (V.fromList + 404 [ V2 0 0, V2 0 2 + 405 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 + 406 , V2 1 1, V2 1 0.7 + 407 , V2 0.3 0.7, V2 0.3 0 ])) !! 7 + 408 + 409 -- E shape + 410 shape14 :: Polygon + 411 shape14 = pCycles (mkPolygon $ V.reverse $ V.fromList + 412 [ V2 0 0, V2 0 2 -- up + 413 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 -- first prong + 414 , V2 1 1, V2 1 0.7, V2 0.3 0.7, V2 0.3 0.3 -- second prong + 415 , V2 1 0.3, V2 1 0 -- last prong + 416 ]) !! 9 + 417 + 418 -- + 419 shape15 :: Polygon + 420 shape15 = mkPolygon $ V.fromList + 421 [ V2 0 0, V2 2 0 + 422 , V2 2 2, V2 1 2 + 423 , V2 1 1, V2 0 1] + 424 + 425 shape16 :: Polygon + 426 shape16 = mkPolygon $ V.fromList + 427 [ V2 0 0, V2 2 0 + 428 , V2 2 1, V2 1 1 + 429 , V2 1 2, V2 0 2] + 430 + 431 shape17 :: Polygon + 432 shape17 = mkPolygon $ V.fromList + 433 [ V2 2 0, V2 2 1 + 434 , V2 1 1, V2 1 2 + 435 , V2 0 2, V2 0 1, V2 0 0 ] + 436 + 437 shape18 :: Polygon + 438 shape18 = mkPolygon $ V.fromList + 439 [ V2 2 0, V2 2 1, V2 2 2 + 440 , V2 1 2, V2 1 1 + 441 , V2 0 1, V2 0 0 ] + 442 + 443 shape19 :: Polygon + 444 shape19 = mkPolygon $ V.fromList + 445 [ V2 (-3) (-3), V2 0 (-1) + 446 , V2 3 (-3), V2 1 0 + 447 , V2 3 3, V2 0 1 + 448 , V2 (-3) 3, V2 (-1) 0 ] + 449 + 450 shape20 :: Polygon + 451 shape20 = mkPolygon $ V.fromList + 452 [ V2 (-3) (-3) + 453 , V2 0 (-1) + 454 , V2 3 (-3) + 455 , V2 5 0 + 456 , V2 2.5 (-2) + 457 , V2 1 0 + 458 , V2 3 3 + 459 , V2 0 1 + 460 , V2 (-3) 3 + 461 , V2 (-1) 0 ] + 462 + 463 shape21 :: Polygon + 464 shape21 = mkPolygon $ V.fromList + 465 [V2 0.0 0.0,V2 1.0 0.0,V2 1.0 1.0,V2 2.0 1.0,V2 2.0 (-1.0),V2 3.0 (-1.0) + 466 ,V2 3.0 2.0,V2 0.0 2.0] + 467 + 468 shape22 :: Polygon + 469 shape22 = pScale 2 $ mkPolygon $ V.fromList + 470 [V2 (-0.17) (-0.08) + 471 ,V2 (-0.34) (-0.21) + 472 ,V2 0.0 0.0 + 473 ,V2 (-0.10) 0.60 + 474 ,V2 (-0.14) 0.19 + 475 ,V2 (-0.05) 0.03 + 476 ] + 477 + 478 shape23 :: Polygon + 479 shape23 = mkPolygon $ V.fromList + 480 [ V2 0 0, V2 4 0 + 481 , V2 4 3, V2 2 3 + 482 , V2 2 2, V2 3 2 + 483 , V2 3 1, V2 1 1 + 484 , V2 1 2, V2 2 2 + 485 , V2 2 3, V2 0 3 ] + 486 + 487 concave :: Polygon + 488 concave = mkPolygon $ + 489 V.fromList [V2 0 0, V2 2 0, V2 2 2, V2 1 1, V2 0 2] + 490 + 491 pMkWinding :: Int -> Polygon + 492 pMkWinding n | n < 1 = error "Polygon must have at least one winding." + 493 pMkWinding n = mkPolygon $ + 494 V.fromList $ p0 : p1 : walkTo p1 1 n (V2 1 0) ++ reverse (walkTo p0 1 (n+2) (V2 (-1) 0)) + 495 where + 496 p0 = V2 0 0 + 497 p1 = V2 0 1 + 498 walkTo at a b dir + 499 | a == b = [] + 500 | otherwise = + 501 let newAt = at + (dir ^* toRational a) + 502 in newAt : walkTo newAt (a+1) b (rot dir) + 503 rot (V2 x y) = + 504 V2 y (-x) + 505 + 506 pDeoverlap :: Polygon -> Polygon + 507 pDeoverlap p = mkPolygon arr + 508 where + 509 arr = V.generate (pSize p) worker + 510 worker 0 = pAccess p 0 + 511 worker n = + 512 if length (V.elemIndices (pAccess p n) (polygonPoints p)) /= 1 + 513 then + 514 let prev = arr V.! (n-1) + 515 this = pAccess p n + 516 in lerp 0.99999 this prev + 517 else pAccess p n + 518 + 519 pCycles :: APolygon a -> [APolygon a] + 520 pCycles p = map (pAdjustOffset p) [0 .. pSize p-1] + 521 + 522 pCycle :: PolyCtx a => APolygon a -> Double -> APolygon a + 523 pCycle p 0 = p + 524 pCycle p t = mkPolygon $ worker 0 0 + 525 where + 526 worker acc i + 527 | segment + acc > limit = + 528 V.singleton (lerp (realToFrac $ (segment + acc - limit)/segment) x y) <> + 529 -- V.drop (i+1) (polygonPoints p) <> + 530 V.fromList (map (pAccess p) [i+1..pSize p-1]) <> + 531 V.fromList (map (pAccess p) [0 .. i]) + 532 -- V.take (i+1) (polygonPoints p) + 533 | i == pSize p-1 = V.fromList (map (pAccess p) [0 .. pSize p-1]) + 534 | otherwise = worker (acc+segment) (i+1) + 535 where + 536 x = pAccess p i + 537 y = pAccess p $ i+1 + 538 segment = distance' x y + 539 len = pCircumference' p + 540 limit = t * len + 541 + 542 pCentroid :: Fractional a => APolygon a -> V2 a + 543 pCentroid p = V2 cx cy + 544 where + 545 a = pArea p + 546 cx = recip (6*a) * V.sum (pMapEdges fnX p) + 547 cy = recip (6*a) * V.sum (pMapEdges fnY p) + 548 fnX (V2 x y) (V2 x' y') = (x+x')*(x*y' - x'*y) + 549 fnY (V2 x y) (V2 x' y') = (y+y')*(x*y' - x'*y) + 550 + 551 {-# INLINE pMapEdges #-} + 552 pMapEdges :: (V2 a -> V2 a -> b) -> APolygon a -> V.Vector b + 553 pMapEdges fn p = V.generate n $ \i -> + 554 if i == n-1 + 555 then fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) + 556 else fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1)) + 557 where + 558 n = pSize p + 559 arr = polygonPoints p + 560 + 561 {-# SPECIALIZE pArea :: APolygon Double -> Double #-} + 562 {-# SPECIALIZE pArea :: APolygon Rational -> Rational #-} + 563 pArea :: (Fractional a) => APolygon a -> a + 564 pArea p = + 565 -- 0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p) + 566 0.5 * worker 0 0 + 567 where + 568 fn (V2 x y) (V2 x' y') = x*y' - x'*y + 569 arr = polygonPoints p + 570 worker !acc i + 571 | i == pSize p - 1 = acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) + 572 | otherwise = + 573 worker (acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1))) (i+1) + 574 + 575 pCircumference :: (Real a, Fractional a) => APolygon a -> a + 576 pCircumference p = sum + 577 [ approxDist (pAccess p i) (pAccess p $ i+1) + 578 | i <- [0 .. pSize p-1]] + 579 + 580 pCircumference' :: (Real a, Fractional a) => APolygon a -> Double + 581 pCircumference' p = sum + 582 [ distance' (pAccess p i) (pAccess p $ i+1) + 583 | i <- [0 .. pSize p-1]] + 584 + 585 + 586 -- Add points by splitting the longest lines in half repeatedly. + 587 pAddPoints :: PolyCtx a => Int -> APolygon a -> APolygon a + 588 pAddPoints n p | n <= 0 = p + 589 pAddPoints n p = pAddPoints (n-1) $ + 590 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] + 591 where + 592 worker idx + 593 | idx == longestEdge = + 594 let start = pAccess p idx + 595 end = pAccess p $ idx+1 + 596 middle = lerp 0.5 end start + 597 in [start, middle] + 598 | otherwise = [pAccess p idx] + 599 longestEdge = maximumBy cmpLength [0 .. pSize p-1] + 600 cmpLength a b = + 601 distSquared (pAccess p a) (pAccess p $ a+1) `compare` + 602 distSquared (pAccess p b) (pAccess p $ b+1) + 603 + 604 pAddPointsRestricted :: PolyCtx a => [(V2 a, V2 a)] -> Int -> APolygon a -> APolygon a + 605 pAddPointsRestricted _immutableEdges n p | n <= 0 = p + 606 pAddPointsRestricted immutableEdges n p = pAddPointsRestricted immutableEdges (n-1) $ + 607 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] + 608 where + 609 isImmutable idx = + 610 (pAccess p idx, pAccess p $ idx+1) `elem` immutableEdges || + 611 (pAccess p $ idx+1, pAccess p idx) `elem` immutableEdges + 612 worker idx + 613 | idx == longestEdge && not (isImmutable idx) = + 614 let start = pAccess p idx + 615 end = pAccess p $ idx+1 + 616 middle = lerp 0.5 end start + 617 in [start, middle] + 618 | otherwise = [pAccess p idx] + 619 longestEdge = maximumBy cmpLength [0 .. pSize p-1] + 620 cmpLength a _ | isImmutable a = LT + 621 cmpLength _ b | isImmutable b = GT + 622 cmpLength a b = + 623 distSquared (pAccess p a) (pAccess p $ a+1) `compare` + 624 distSquared (pAccess p b) (pAccess p $ b+1) + 625 + 626 pAddPointsBetween :: PolyCtx a => (Int, Int) -> Int -> APolygon a -> APolygon a + 627 pAddPointsBetween _ n p | n <= 0 = p + 628 pAddPointsBetween (i,l) n p = pAddPointsBetween (i,l+1) (n-1) $ + 629 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] + 630 where + 631 worker idx + 632 | idx == longestEdge = + 633 let start = pAccess p idx + 634 end = pAccess p $ idx+1 + 635 middle = lerp 0.5 end start + 636 in [start, middle] + 637 | otherwise = [pAccess p idx] + 638 longestEdge = maximumBy cmpLength [i .. i+l-1] + 639 cmpLength a b = + 640 distSquared (pAccess p a) (pAccess p $ a+1) `compare` + 641 distSquared (pAccess p b) (pAccess p $ b+1) + 642 + 643 -- addPoints :: Int -> Polygon -> Polygon + 644 -- addPoints n p = mkPolygon $ V.fromList $ worker n 0 (map (pAccess p) [0..s]) + 645 -- where + 646 -- worker 0 _ rest = init rest + 647 -- worker i acc (x:y:xs) = + 648 -- let xy = approxDist x y in + 649 -- if acc + xy > limit + 650 -- then x : worker (i-1) 0 (lerp ((limit-acc)/xy) y x : y:xs) + 651 -- else x : worker i (acc+xy) (y:xs) + 652 -- worker _ _ [_] = [] + 653 -- worker _ _ _ = error "addPoints: invalid polygon" + 654 -- s = pSize p + 655 -- len = polygonLength p + 656 -- limit = len / fromIntegral (n+1) + 657 + 658 pIsConvex :: Polygon -> Bool + 659 pIsConvex p = and + 660 [ area2X (pAccess p i) (pAccess p j) (pAccess p k) > 0 + 661 | i <- [0..n-1] + 662 , j <- [i+1..n-1] + 663 , k <- [j+1..n-1] + 664 ] + 665 where n = pSize p + 666 + 667 pIsCCW :: Polygon -> Bool + 668 pIsCCW p | pNull p = False + 669 pIsCCW p = V.sum (pMapEdges fn p) < 0 + 670 where + 671 fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1) + 672 + 673 {-# INLINE pRayIntersect #-} + 674 pRayIntersect :: PolyCtx a => APolygon a -> (Int, Int) -> (Int,Int) -> Maybe (V2 a) + 675 pRayIntersect p (a,b) (c,d) = + 676 rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d) + 677 + 678 pCuts :: (Real a, Fractional a, Epsilon a) => APolygon a -> [(APolygon a,APolygon a)] + 679 pCuts p = + 680 [ pCutAt (pAdjustOffset p i) (j-i) + 681 | i <- [0 .. pSize p-1 ] + 682 , j <- [i+2 .. pSize p-1 ] + 683 , (j+1) `mod` pSize p /= i + 684 , pParent p i j == i ] + 685 + 686 pCutEqual :: PolyCtx a => APolygon a -> (APolygon a, APolygon a) + 687 pCutEqual p = + 688 fromMaybe (p,p) $ listToMaybe $ sortOn f $ pCuts p + 689 where + 690 f (a,b) = abs (pArea a - pArea b) + 691 + 692 -- FIXME: This should be more efficient + 693 pCutAt :: PolyCtx a => APolygon a -> Int -> (APolygon a, APolygon a) + 694 pCutAt p i = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right) + 695 where + 696 n = pSize p + 697 left = map (pAccess p) [0 .. i] + 698 right = map (pAccess p) (0:[i..n-1]) + 699 + 700 pOverlap :: PolyCtx a => APolygon a -> APolygon a -> APolygon a + 701 pOverlap a b = mkPolygon $ V.fromList $ clearDups $ concatMap edgeIntersect [0 .. pSize a-1] + 702 where + 703 clearDups (x:y:xs) + 704 | x == y = clearDups (y:xs) + 705 | otherwise = x : clearDups (y:xs) + 706 clearDups xs = xs + 707 edgeIntersect edge = + 708 sortOn (distSquared (pAccess a edge)) $ catMaybes + 709 [ lineIntersect (aP, aP') (bP, bP') + 710 | i <- [0 .. pSize b-1] + 711 , let aP = pAccess a edge + 712 aP' = pAccess a (edge+1) + 713 bP = pAccess b i + 714 bP' = pAccess b (i+1) + 715 ] + 716 + 717 --------------------------------------------------------- + 718 -- SSSP visibility and SSSP windows + 719 + 720 ssspVisibility :: PolyCtx a => APolygon a -> APolygon a + 721 ssspVisibility p = mkPolygon $ + 722 V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1]) + 723 where + 724 clearDups (x:y:xs) + 725 | x == y = clearDups (y:xs) + 726 | otherwise = x : clearDups (y:xs) + 727 clearDups xs = xs + 728 obstructedBy n = + 729 case pParent p 0 n of + 730 0 -> n + 731 i -> obstructedBy i + 732 go [] = [] + 733 go [x] = [pAccess p x] + 734 go (x:y:xs) = + 735 let xO = obstructedBy x + 736 yO = obstructedBy y + 737 in case () of + 738 () + 739 -- Both ends are visible. + 740 | xO == x && yO == y -> pAccess p x : go (y:xs) + 741 -- X is visible, x to intersect (0,yO) (x,y) + 742 | xO == x -> + 743 pAccess p x : fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y)) : go (y:xs) + 744 -- Y is visible + 745 | yO == y -> fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) : pAccess p y : go (y:xs) + 746 -- Neither is visible and they've obstructed by the same point + 747 -- so the entire edge is hidden. + 748 | xO == yO -> go (y:xs) + 749 -- Neither is visible. Cast shadow from obstruction points to + 750 -- find if a subsection of the edge is visible. + 751 | otherwise -> + 752 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y)) + 753 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y)) + 754 in if a /= b + 755 then a : b : go (y:xs) + 756 else go (y:xs) + 757 + 758 ssspWindows :: Polygon -> [(V2 Rational, V2 Rational)] + 759 ssspWindows p = clearDups $ go (pAccess p 0) [0..pSize p-1] + 760 where + 761 clearDups (x:y:xs) + 762 | x == y = clearDups (y:xs) + 763 | otherwise = x : clearDups (y:xs) + 764 clearDups xs = xs + 765 obstructedBy n = + 766 case pParent p 0 n of + 767 0 -> n + 768 i -> obstructedBy i + 769 go _ [] = [] + 770 go _ [_] = [] + 771 go l (x:y:xs) = + 772 let xO = obstructedBy x + 773 yO = obstructedBy y + 774 in case () of + 775 () + 776 -- Both ends are visible. + 777 | xO == x && yO == y -> go (pAccess p x) (y:xs) + 778 -- X is visible, x to intersect (0,yO) (x,y) + 779 | xO == x -> + 780 go (fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y))) (y:xs) + 781 -- Y is visible + 782 | yO == y -> + 783 let newL = fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) in + 784 (l, newL) : + 785 go newL (y:xs) + 786 -- Neither is visible and they've obstructed by the same point + 787 -- so the entire edge is hidden. + 788 | xO == yO -> go l (y:xs) + 789 -- Neither is visible. Cast shadow from obstruction points to + 790 -- find if a subsection of the edge is visible. + 791 | otherwise -> + 792 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y)) + 793 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y)) + 794 in if a /= b + 795 then (l, a) : (b, pAccess p yO) : go (pAccess p yO) (y:xs) + 796 else go l (y:xs) diff --git a/reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP.hs.html b/reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP.hs.html index 4ff5029..7555e6b 100644 --- a/reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP.hs.html +++ b/reanimate-0.4.3.0-inplace/Reanimate.Math.SSSP.hs.html @@ -19,347 +19,371 @@ span.spaces { background: white } 
     1 {-# LANGUAGE FlexibleInstances     #-}
     2 {-# LANGUAGE MultiParamTypeClasses #-}
-    3 {-# OPTIONS_GHC -fno-warn-orphans #-}
-    4 {-# OPTIONS_HADDOCK hide #-}
-    5 module Reanimate.Math.SSSP
-    6   ( -- * Single-Source-Shortest-Path
-    7     SSSP
-    8   , sssp                -- :: (Fractional a, Ord a) => Ring a -> Dual -> SSSP
-    9   , dual                -- :: Int -> Triangulation -> Dual
-   10   , Dual(..)
-   11   , DualTree(..)
-   12   , PDual
-   13   , toPDual             -- :: Ring Rational -> Dual -> PDual
-   14   , pdualRings          -- :: Ring Rational -> PDual -> [Ring Rational]
-   15     -- * Misc
-   16   , dualToTriangulation -- :: Ring Rational -> Dual -> Triangulation
-   17   , pdualReduce         -- :: Ring Rational -> PDual -> Int -> PDual
-   18   , visibilityArray     -- :: Ring Rational -> V.Vector [Int]
-   19   , naive               -- :: Ring Rational -> SSSP
-   20   , naive2              -- :: Ring Rational -> SSSP
-   21   , drawDual            -- :: Dual -> String
-   22   ) where
-   23 
-   24 import           Control.Monad
-   25 -- import           Control.Exception
-   26 import           Control.Monad.ST
-   27 -- import           Data.FingerTree            (SearchResult (..), (|>))
-   28 -- import qualified Data.FingerTree            as F
-   29 import           Data.Foldable
-   30 import           Data.List
-   31 import qualified Data.Map                   as Map
-   32 import           Data.Maybe
-   33 import           Data.Ord
-   34 import           Data.STRef
-   35 import           Data.Tree
-   36 import qualified Data.Vector                as V
-   37 import qualified Data.Vector.Mutable        as MV
-   38 import           Reanimate.Math.Common
-   39 import           Reanimate.Math.Triangulate
+    3 {-# LANGUAGE RecordWildCards       #-}
+    4 {-# OPTIONS_GHC -fno-warn-orphans #-}
+    5 {-# OPTIONS_HADDOCK hide #-}
+    6 module Reanimate.Math.SSSP
+    7   ( -- * Single-Source-Shortest-Path
+    8     SSSP
+    9   , sssp                -- :: (Fractional a, Ord a) => Ring a -> Dual -> SSSP
+   10   , ssspFinger
+   11   , dual                -- :: Int -> Triangulation -> Dual
+   12   , Dual(..)
+   13   , DualTree(..)
+   14     -- * Misc
+   15   , dualToTriangulation -- :: Ring Rational -> Dual -> Triangulation
+   16   , visibilityArray     -- :: Ring Rational -> V.Vector [Int]
+   17   , naive               -- :: Ring Rational -> SSSP
+   18   , naive2              -- :: Ring Rational -> SSSP
+   19   , drawDual            -- :: Dual -> String
+   20   ) where
+   21 
+   22 import           Control.Monad
+   23 import           Control.Monad.ST
+   24 import qualified Data.FingerTree            as F
+   25 import           Data.Foldable
+   26 import           Data.List
+   27 import qualified Data.Map                   as Map
+   28 import           Data.Maybe
+   29 import           Data.STRef
+   30 import           Data.Tree
+   31 import qualified Data.Vector                as V
+   32 import qualified Data.Vector.Mutable        as MV
+   33 import           Reanimate.Math.Common
+   34 import           Reanimate.Math.Triangulate
+   35 
+   36 -- import           Debug.Trace
+   37 
+   38 type SSSP = V.Vector Int
+   39 
    40 
-   41 -- import           Debug.Trace
-   42 
-   43 type SSSP = V.Vector Int
-   44 
-   45 
-   46 -- ssspParent :: Polygon -> SSSP -> Int -> Int
-   47 -- ssspParent p sTree x =
-   48 --     (sTree V.! ((x - polygonOffset p) `mod` n) + polygonOffset p) `mod` n
-   49 --   where
-   50 --     n = polygonSize p
-   51 
-   52 visibilityArray :: Ring Rational -> V.Vector [Int]
-   53 visibilityArray p = arr
-   54   where
-   55     n = ringSize p
-   56     arr = V.fromList
-   57         [ visibility y
-   58         | y <- [0..n-1]
-   59         ]
-   60     visibility y =
-   61       [ i
-   62       | i <- [0..y-1]
-   63       , y `elem` arr V.! i ] ++
-   64       [ i
-   65       | i <- [y+1 .. n-1]
-   66       , let pI = ringAccess p i
-   67             isOpen = isRightTurn pYp pY pYn
-   68       , ringClamp p (y+1) == i || ringClamp p (y-1) == i || if isOpen
-   69         then isLeftTurnOrLinear pY pYn pI ||
-   70              isLeftTurnOrLinear pYp pY pI
-   71         else not $ isRightTurn pY pYn pI ||
-   72                    isRightTurn pYp pY pI
-   73       , let myEdges = [(e1,e2) | (e1,e2) <- edges, e1/=y, e1/=i, e2/=y,e2/=i]
-   74       , all (isNothing . lineIntersect (pY,pI))
-   75               [ (ringAccess p e1, ringAccess p e2) | (e1,e2) <- myEdges ]]
-   76       where
-   77         pY = ringAccess p y
-   78         pYn = ringAccess p $ y+1
-   79         pYp = ringAccess p $ y-1
-   80         edges = zip [0..n-1] (tail [0..n-1] ++ [0])
-   81 
-   82 
-   83 
-   84 -- Iterative Single Source Shortest Path solver. Quite slow.
-   85 naive :: Ring Rational -> SSSP
-   86 naive p =
-   87     V.fromList $ Map.elems $
-   88     Map.map snd $
-   89     worker initial
-   90   where
-   91     initial = Map.singleton 0 (0,0)
-   92     visibility = visibilityArray p
-   93     worker :: Map.Map Int (Rational, Int) -> Map.Map Int (Rational, Int)
-   94     worker m
-   95         | m==newM   = newM
-   96         | otherwise = worker newM
-   97       where
-   98         ms' = [ Map.fromList
-   99                     [ case Map.lookup v m of
-  100                         Nothing -> (v, (distThroughI, i))
-  101                         Just (otherDist,parent)
-  102                           | otherDist > distThroughI -> (v, (distThroughI, i))
-  103                           | otherwise -> (v, (otherDist, parent))
-  104                     | v <- visibility V.! i
-  105                     , let distThroughI = dist + approxDist (ringAccess p i) (ringAccess p v) ]
-  106               | (i,(dist,_)) <- Map.toList m
-  107               ]
-  108         newM = Map.unionsWith g (m:ms') :: Map.Map Int (Rational,Int)
-  109     g a b = if fst a < fst b then a else b
-  110 
-  111 naive2 :: Ring Rational -> SSSP
-  112 naive2 p = runST $ do
-  113     parents <- MV.replicate (ringSize p) (-1)
-  114     costs <- MV.replicate (ringSize p) (-1)
-  115     MV.write parents 0 0
-  116     MV.write costs 0 0
-  117     changedRef <- newSTRef False
-  118     let loop i
-  119           | i == ringSize p = do
-  120             changed <- readSTRef changedRef
-  121             when changed $ do
-  122               writeSTRef changedRef False
-  123               loop 0
-  124           | otherwise = do
-  125             myCost <- MV.read costs i
-  126             unless (myCost < 0) $
-  127               forM_ (visibility V.! i) $ \n -> do
-  128                 -- n is visible from i.
-  129                 theirCost <- MV.read costs n
-  130                 let throughCost = myCost + approxDist (ringAccess p i) (ringAccess p n)
-  131                 when (throughCost < theirCost || theirCost < 0) $ do
-  132                     MV.write parents n i
-  133                     MV.write costs n throughCost
-  134                     writeSTRef changedRef True
-  135             loop (i+1)
-  136     loop 0
-  137     V.unsafeFreeze parents
-  138   where
-  139     visibility = visibilityArray p
-  140 
-  141 data PDual = PDual (V.Vector Int) Rational [PDual]
-  142   deriving (Show)
-  143 
-  144 toPDual :: Ring Rational -> Dual -> PDual
-  145 toPDual p d =
-  146   case d of
-  147     Dual (a,b,c) l r ->
-  148       PDual (V.fromList [a,b,c])
-  149         (area2X (ringAccess p a) (ringAccess p b) (ringAccess p c))
-  150         (catMaybes [ worker c a l, worker b c r])
-  151   where
-  152     worker _ _ EmptyDual = Nothing
-  153     worker a b (NodeDual x l r) = Just $
-  154       PDual (V.fromList [a,x,b])
-  155         (area2X (ringAccess p a) (ringAccess p x) (ringAccess p b))
-  156         (catMaybes [ worker x b l, worker a x r])
+   41 -- ssspParent :: Polygon -> SSSP -> Int -> Int
+   42 -- ssspParent p sTree x =
+   43 --     (sTree V.! ((x - polygonOffset p) `mod` n) + polygonOffset p) `mod` n
+   44 --   where
+   45 --     n = polygonSize p
+   46 
+   47 visibilityArray :: Ring Rational -> V.Vector [Int]
+   48 visibilityArray p = arr
+   49   where
+   50     n = ringSize p
+   51     arr = V.fromList
+   52         [ visibility y
+   53         | y <- [0..n-1]
+   54         ]
+   55     visibility y =
+   56       [ i
+   57       | i <- [0..y-1]
+   58       , y `elem` arr V.! i ] ++
+   59       [ i
+   60       | i <- [y+1 .. n-1]
+   61       , let pI = ringAccess p i
+   62             isOpen = isRightTurn pYp pY pYn
+   63       , ringClamp p (y+1) == i || ringClamp p (y-1) == i || if isOpen
+   64         then isLeftTurnOrLinear pY pYn pI ||
+   65              isLeftTurnOrLinear pYp pY pI
+   66         else not $ isRightTurn pY pYn pI ||
+   67                    isRightTurn pYp pY pI
+   68       , let myEdges = [(e1,e2) | (e1,e2) <- edges, e1/=y, e1/=i, e2/=y,e2/=i]
+   69       , all (isNothing . lineIntersect (pY,pI))
+   70               [ (ringAccess p e1, ringAccess p e2) | (e1,e2) <- myEdges ]]
+   71       where
+   72         pY = ringAccess p y
+   73         pYn = ringAccess p $ y+1
+   74         pYp = ringAccess p $ y-1
+   75         edges = zip [0..n-1] (tail [0..n-1] ++ [0])
+   76 
+   77 
+   78 
+   79 -- Iterative Single Source Shortest Path solver. Quite slow.
+   80 naive :: Ring Rational -> SSSP
+   81 naive p =
+   82     V.fromList $ Map.elems $
+   83     Map.map snd $
+   84     worker initial
+   85   where
+   86     initial = Map.singleton 0 (0,0)
+   87     visibility = visibilityArray p
+   88     worker :: Map.Map Int (Rational, Int) -> Map.Map Int (Rational, Int)
+   89     worker m
+   90         | m==newM   = newM
+   91         | otherwise = worker newM
+   92       where
+   93         ms' = [ Map.fromList
+   94                     [ case Map.lookup v m of
+   95                         Nothing -> (v, (distThroughI, i))
+   96                         Just (otherDist,parent)
+   97                           | otherDist > distThroughI -> (v, (distThroughI, i))
+   98                           | otherwise -> (v, (otherDist, parent))
+   99                     | v <- visibility V.! i
+  100                     , let distThroughI = dist + approxDist (ringAccess p i) (ringAccess p v) ]
+  101               | (i,(dist,_)) <- Map.toList m
+  102               ]
+  103         newM = Map.unionsWith g (m:ms') :: Map.Map Int (Rational,Int)
+  104     g a b = if fst a < fst b then a else b
+  105 
+  106 naive2 :: Ring Rational -> SSSP
+  107 naive2 p = runST $ do
+  108     parents <- MV.replicate (ringSize p) (-1)
+  109     costs <- MV.replicate (ringSize p) (-1)
+  110     MV.write parents 0 0
+  111     MV.write costs 0 0
+  112     changedRef <- newSTRef False
+  113     let loop i
+  114           | i == ringSize p = do
+  115             changed <- readSTRef changedRef
+  116             when changed $ do
+  117               writeSTRef changedRef False
+  118               loop 0
+  119           | otherwise = do
+  120             myCost <- MV.read costs i
+  121             unless (myCost < 0) $
+  122               forM_ (visibility V.! i) $ \n -> do
+  123                 -- n is visible from i.
+  124                 theirCost <- MV.read costs n
+  125                 let throughCost = myCost + approxDist (ringAccess p i) (ringAccess p n)
+  126                 when (throughCost < theirCost || theirCost < 0) $ do
+  127                     MV.write parents n i
+  128                     MV.write costs n throughCost
+  129                     writeSTRef changedRef True
+  130             loop (i+1)
+  131     loop 0
+  132     V.unsafeFreeze parents
+  133   where
+  134     visibility = visibilityArray p
+  135 
+  136 -- Dual of triangulated polygon
+  137 data Dual = Dual (Int,Int,Int) -- (a,b,c)
+  138                   DualTree -- borders ca
+  139                   DualTree -- borders bc
+  140   deriving (Show)
+  141 
+  142 data DualTree
+  143   = EmptyDual
+  144   | NodeDual Int -- axb triangle, a and b are from parent.
+  145       DualTree -- borders xb
+  146       DualTree -- borders ax
+  147   deriving (Show)
+  148 
+  149 drawDual :: Dual -> String
+  150 drawDual d = drawTree $
+  151   case d of
+  152     Dual (a,b,c) l r -> Node (show (a,b,c)) [worker c a l, worker b c r]
+  153   where
+  154     worker _a _b EmptyDual = Node "Leaf" []
+  155     worker a b (NodeDual x l r) =
+  156       Node (show (b,a,x)) [worker x b l, worker a x r]
   157 
-  158 pdualSize :: PDual -> Int
-  159 pdualSize (PDual _ _ children) = 1 + sum (map pdualSize children)
-  160 
-  161 pdualArea :: PDual -> Rational
-  162 pdualArea (PDual _ faceArea _) = faceArea
-  163 
-  164 -- FIXME: 'origin' isn't used. Remove.
-  165 pdualReduce :: Ring Rational -> PDual -> Int -> PDual
-  166 pdualReduce origin pdual n
-  167   | pdualSize pdual <= n = pdual
-  168   | otherwise =
-  169     let smallest = minimum $ pAreas pdual
-  170     in pdualReduce origin (merge smallest pdual) n
-  171   where
-  172     merge _s (PDual p faceArea []) = PDual p faceArea []
-  173     merge s (PDual p faceArea children)
-  174       | faceArea == s =
-  175         let (PDual p2 area2 children2:xs) = sortBy (comparing pdualArea) children
-  176         in PDual (joinP p p2) (faceArea+area2) (children2++xs)
-  177       | otherwise =
-  178         let (PDual p2 area2 children2:xs) = sortBy (comparing pdualArea) children
-  179         in if area2 == s
-  180             then PDual (joinP p p2) (faceArea+area2) (children2++xs)
-  181             else PDual p faceArea (map (merge s) children)
-  182     pAreas (PDual _ faceArea children) = faceArea : concatMap pAreas children
-  183     joinP a b = V.fromList (sort (V.toList a ++ V.toList b))
-  184 
-  185 pdualRings :: Ring Rational -> PDual -> [Ring Rational]
-  186 pdualRings p (PDual pts _area children) =
-  187   ringPack (V.map (ringAccess p) pts) : concatMap (pdualRings p) children
-  188 
-  189 -- Dual of triangulated polygon
-  190 data Dual = Dual (Int,Int,Int) -- (a,b,c)
-  191                   DualTree -- borders ca
-  192                   DualTree -- borders bc
-  193   deriving (Show)
-  194 
-  195 data DualTree
-  196   = EmptyDual
-  197   | NodeDual Int -- axb triangle, a and b are from parent.
-  198       DualTree -- borders xb
-  199       DualTree -- borders ax
-  200   deriving (Show)
-  201 
-  202 drawDual :: Dual -> String
-  203 drawDual d = drawTree $
-  204   case d of
-  205     Dual (a,b,c) l r -> Node (show (a,b,c)) [worker c a l, worker b c r]
-  206   where
-  207     worker _a _b EmptyDual = Node "Leaf" []
-  208     worker a b (NodeDual x l r) =
-  209       Node (show (b,a,x)) [worker x b l, worker a x r]
+  158 dualToTriangulation :: Ring Rational -> Dual -> Triangulation
+  159 dualToTriangulation p d = edgesToTriangulation (ringSize p) $ filter goodEdge $
+  160     case d of
+  161       Dual (a,b,c) l r ->
+  162         (a,b):(a,c):(b,c):worker c a l ++ worker b c r
+  163   where
+  164     goodEdge (a,b)
+  165       = a /= ringClamp p (b+1) && a /= ringClamp p (b-1)
+  166     worker _a _b EmptyDual = []
+  167     worker a b (NodeDual x l r) =
+  168       (a,x) : (x, b) : worker x b l ++ worker a x r
+  169 
+  170 -- Dual path:
+  171 -- (Int,Int,Int) + V.Vector Int + V.Vector LeftOrRight
+  172 
+  173 -- simplifyDual :: DualTree -> DualTree
+  174 -- -- simplifyDual (NodeDual x EmptyDual EmptyDual) = NodeLeaf x
+  175 -- -- simplifyDual (NodeDual x l EmptyDual) = NodeDualL x l
+  176 -- -- simplifyDual (NodeDual x EmptyDual r) = NodeDualR x r
+  177 -- simplifyDual d = d
+  178 
+  179 dual :: Int -> Triangulation -> Dual
+  180 dual root t =
+  181   case hasTriangle of
+  182     []    -> error "weird triangulation"
+  183     -- [] -> Dual (0,1,V.length t-1) EmptyDual (dualTree t (1, (V.length t-1)) 0)
+  184     (x:_) -> Dual (root,rootNext,x) (dualTree t (x,root) rootNext) (dualTree t (rootNext,x) root)
+  185   where
+  186     rootNext = idx (root+1)
+  187     rootPrev = idx (root-1)
+  188     rootNNext = idx (root+2)
+  189     idx i = i `mod` n
+  190     hasTriangle = (rootPrev : t V.! root) `intersect` (rootNNext : t V.! rootNext)
+  191     n = V.length t
+  192 
+  193 -- a=6, b=0, e=1
+  194 dualTree :: Triangulation -> (Int,Int) -> Int -> DualTree
+  195 dualTree t (a,b) e = -- simplifyDual $
+  196     case hasTriangle of
+  197       [] -> EmptyDual
+  198       [(ab)] ->
+  199         NodeDual ab
+  200           (dualTree t (ab,b) a)
+  201           (dualTree t (a,ab) b)
+  202       _ -> error $ "Invalid triangulation: " ++ show (a,b,e,hasTriangle)
+  203   where
+  204     hasTriangle = (prev a : next a : t V.! a) `intersect` (prev b : next b : t V.! b)
+  205       \\ [e]
+  206     n = V.length t
+  207     next x = (x+1) `mod` n
+  208     prev x = (x-1) `mod` n
+  209 
   210 
-  211 dualToTriangulation :: Ring Rational -> Dual -> Triangulation
-  212 dualToTriangulation p d = edgesToTriangulation (ringSize p) $ filter goodEdge $
-  213     case d of
-  214       Dual (a,b,c) l r ->
-  215         (a,b):(a,c):(b,c):worker c a l ++ worker b c r
-  216   where
-  217     goodEdge (a,b)
-  218       = a /= ringClamp p (b+1) && a /= ringClamp p (b-1)
-  219     worker _a _b EmptyDual = []
-  220     worker a b (NodeDual x l r) =
-  221       (a,x) : (x, b) : worker x b l ++ worker a x r
-  222 
-  223 -- Dual path:
-  224 -- (Int,Int,Int) + V.Vector Int + V.Vector LeftOrRight
-  225 
-  226 -- simplifyDual :: DualTree -> DualTree
-  227 -- -- simplifyDual (NodeDual x EmptyDual EmptyDual) = NodeLeaf x
-  228 -- -- simplifyDual (NodeDual x l EmptyDual) = NodeDualL x l
-  229 -- -- simplifyDual (NodeDual x EmptyDual r) = NodeDualR x r
-  230 -- simplifyDual d = d
-  231 
-  232 dual :: Int -> Triangulation -> Dual
-  233 dual root t =
-  234   case hasTriangle of
-  235     []    -> error "weird triangulation"
-  236     -- [] -> Dual (0,1,V.length t-1) EmptyDual (dualTree t (1, (V.length t-1)) 0)
-  237     (x:_) -> Dual (root,rootNext,x) (dualTree t (x,root) rootNext) (dualTree t (rootNext,x) root)
-  238   where
-  239     rootNext = idx (root+1)
-  240     rootPrev = idx (root-1)
-  241     rootNNext = idx (root+2)
-  242     idx i = i `mod` n
-  243     hasTriangle = (rootPrev : t V.! root) `intersect` (rootNNext : t V.! rootNext)
-  244     n = V.length t
-  245 
-  246 -- a=6, b=0, e=1
-  247 dualTree :: Triangulation -> (Int,Int) -> Int -> DualTree
-  248 dualTree t (a,b) e = -- simplifyDual $
-  249     case hasTriangle of
-  250       [] -> EmptyDual
-  251       [(ab)] ->
-  252         NodeDual ab
-  253           (dualTree t (ab,b) a)
-  254           (dualTree t (a,ab) b)
-  255       _ -> error $ "Invalid triangulation: " ++ show (a,b,e,hasTriangle)
-  256   where
-  257     hasTriangle = (prev a : next a : t V.! a) `intersect` (prev b : next b : t V.! b)
-  258       \\ [e]
-  259     n = V.length t
-  260     next x = (x+1) `mod` n
-  261     prev x = (x-1) `mod` n
-  262 
-  263 -- data MinMax = MinMax Int Int | MinMaxEmpty deriving (Show)
-  264 -- instance Semigroup MinMax where
-  265 --   MinMaxEmpty <> b = b
-  266 --   a <> MinMaxEmpty = a
-  267 --   MinMax a b <> MinMax c d
-  268 --     = MinMax (min a c) (max b d)
-  269 --     -- = MinMax c b
-  270 -- instance Monoid MinMax where
-  271 --   mempty = MinMaxEmpty
-  272 --
-  273 -- instance F.Measured MinMax Int where
-  274 --   measure i = MinMax i i
-  275 
-  276 -- dualRoot :: Dual -> Int
-  277 -- dualRoot (Dual (a,_,_) _ _) = a
-  278 
-  279 -- O(n*ln n), could be O(n) if I could figure out how to use fingertrees...
-  280 sssp :: (Fractional a, Ord a, Epsilon a) => Ring a -> Dual -> SSSP
-  281 sssp p d = toSSSP $
-  282     case d of
-  283       Dual (a,b,c) l r ->
-  284         (a, a) :
-  285         (b, a) :
-  286         (c, a) :
-  287         worker [c] [b] a r ++
-  288         loopLeft a c l
-  289   where
-  290     toSSSP edges =
-  291       (V.fromList . map snd . sortOn fst) edges
-  292     loopLeft a outer l =
-  293       case l of
-  294         EmptyDual -> []
-  295         NodeDual x l' r' ->
-  296           (x,a) :
-  297           worker [x] [outer] a r' ++
-  298           loopLeft a x l'
-  299     searchFn _checkStep _cusp _x [] = Nothing
-  300     searchFn checkStep cusp x (y:ys)
-  301       | not (checkStep (ringAccess p cusp) (ringAccess p y) (ringAccess p x))
-  302         = Just $ helper [] y ys
-  303       | otherwise = Nothing
-  304       where
-  305         helper acc v [] = (v, [], reverse acc)
-  306         helper acc v1 (v2:vs)
-  307           | checkStep (ringAccess p v1) (ringAccess p v2) (ringAccess p x) =
-  308             (v1, v2:vs, reverse acc)
-  309           | otherwise = helper (v1:acc) v2 vs
-  310     searchRight = searchFn isLeftTurn
-  311     searchLeft = searchFn isRightTurn
-  312     -- adj x = x -- ringClamp p (x-dualRoot d)
-  313     -- optTrace msg =
-  314     --   if False -- dualRoot d == 1 || dualRoot d == 0
-  315     --     then trace msg
-  316     --     else id
-  317     worker _ _ _ EmptyDual = []
-  318     worker f1 f2 cusp (NodeDual x l r) =
-  319         -- (optTrace ("Funnel: " ++ show
-  320         --       (map adj $ toList f1
-  321         --       ,adj cusp
-  322         --       ,map adj $ toList f2
-  323         --       ,adj x
-  324         --       , dualRoot d))
-  325         --   ) $
-  326         case searchLeft cusp x (toList f1) of
-  327           Just (v, f1Hi, f1Lo) ->
-  328                 -- optTrace ("  Visble from left: " ++ show (adj x,adj v)) $
-  329                 (x, v::Int) :
-  330                 worker f1Hi [x] v l ++
-  331                 worker (f1Lo ++ [v, x]) f2 cusp r
-  332           Nothing ->
-  333             case searchRight cusp x (toList f2) of
-  334               Just (v, f2Hi, f2Lo) ->
-  335                 -- optTrace ("  Visble from right: " ++ show (adj x,adj v)) $
-  336                 (x, v::Int) :
-  337                 worker f1 (f2Lo ++ [v, x]) cusp l ++
-  338                 worker [x] f2Hi v r
-  339               Nothing ->
-  340                 -- optTrace ("  Visble from cusp: " ++ show (adj x,adj cusp)) $
-  341                 (x, cusp::Int) :
-  342                 worker f1 [x] cusp l ++
-  343                 worker [x] f2 cusp r
+  211 -- dualRoot :: Dual -> Int
+  212 -- dualRoot (Dual (a,_,_) _ _) = a
+  213 
+  214 -- O(n*ln n), could be O(n) if I could figure out how to use fingertrees...
+  215 sssp :: (Fractional a, Ord a, Epsilon a) => Ring a -> Dual -> SSSP
+  216 sssp p d = toSSSP $
+  217     case d of
+  218       Dual (a,b,c) l r ->
+  219         (a, a) :
+  220         (b, a) :
+  221         (c, a) :
+  222         worker [c] [b] a r ++
+  223         loopLeft a c l
+  224   where
+  225     toSSSP edges =
+  226       (V.fromList . map snd . sortOn fst) edges
+  227     loopLeft a outer l =
+  228       case l of
+  229         EmptyDual -> []
+  230         NodeDual x l' r' ->
+  231           (x,a) :
+  232           worker [x] [outer] a r' ++
+  233           loopLeft a x l'
+  234     searchFn _checkStep _cusp _x [] = Nothing
+  235     searchFn checkStep cusp x (y:ys)
+  236       | not (checkStep (ringAccess p cusp) (ringAccess p y) (ringAccess p x))
+  237         = Just $ helper [] y ys
+  238       | otherwise = Nothing
+  239       where
+  240         helper acc v [] = (v, [], reverse acc)
+  241         helper acc v1 (v2:vs)
+  242           | checkStep (ringAccess p v1) (ringAccess p v2) (ringAccess p x) =
+  243             (v1, v2:vs, reverse acc)
+  244           | otherwise = helper (v1:acc) v2 vs
+  245     searchRight = searchFn isLeftTurn
+  246     searchLeft = searchFn isRightTurn
+  247     -- adj x = x -- ringClamp p (x-dualRoot d)
+  248     -- optTrace msg =
+  249     --   if False -- dualRoot d == 1 || dualRoot d == 0
+  250     --     then trace msg
+  251     --     else id
+  252     worker _ _ _ EmptyDual = []
+  253     worker f1 f2 cusp (NodeDual x l r) =
+  254         -- (optTrace ("Funnel: " ++ show
+  255         --       (map adj $ toList f1
+  256         --       ,adj cusp
+  257         --       ,map adj $ toList f2
+  258         --       ,adj x
+  259         --       , dualRoot d))
+  260         --   ) $
+  261         case searchLeft cusp x (toList f1) of
+  262           Just (v, f1Hi, f1Lo) ->
+  263                 -- optTrace ("  Visble from left: " ++ show (adj x,adj v)) $
+  264                 (x, v::Int) :
+  265                 worker f1Hi [x] v l ++
+  266                 worker (f1Lo ++ [v, x]) f2 cusp r
+  267           Nothing ->
+  268             case searchRight cusp x (toList f2) of
+  269               Just (v, f2Hi, f2Lo) ->
+  270                 -- optTrace ("  Visble from right: " ++ show (adj x,adj v)) $
+  271                 (x, v::Int) :
+  272                 worker f1 (f2Lo ++ [v, x]) cusp l ++
+  273                 worker [x] f2Hi v r
+  274               Nothing ->
+  275                 -- optTrace ("  Visble from cusp: " ++ show (adj x,adj cusp)) $
+  276                 (x, cusp::Int) :
+  277                 worker f1 [x] cusp l ++
+  278                 worker [x] f2 cusp r
+  279 
+  280 data MinMax = MinMax Int Int | MinMaxEmpty deriving (Show)
+  281 instance Semigroup MinMax where
+  282   MinMaxEmpty <> b = b
+  283   a <> MinMaxEmpty = a
+  284   MinMax a _b <> MinMax _c d
+  285     = MinMax a d
+  286 instance Monoid MinMax where
+  287   mempty = MinMaxEmpty
+  288 
+  289 type Chain = F.FingerTree MinMax Int
+  290 data Funnel = Funnel
+  291   { funnelLeft  :: Chain
+  292   , funnelCusp  :: Int
+  293   , funnelRight :: Chain
+  294   }
+  295 
+  296 instance F.Measured MinMax Int where
+  297   measure i = MinMax i i
+  298 
+  299 splitFunnel :: (Epsilon a, Fractional a, Ord a) => Ring a -> Int -> Funnel -> (Int, Funnel, Funnel)
+  300 splitFunnel p x Funnel{..}
+  301     | isOnLeftChain =
+  302       case doSearch isRightTurn funnelLeft of
+  303         (lower, t, upper) ->
+  304           ( t
+  305           , Funnel upper t (F.singleton x)
+  306           , Funnel (lower F.|> t F.|> x) funnelCusp funnelRight)
+  307     | isOnRightChain =
+  308       case doSearch isLeftTurn funnelRight of
+  309         (lower, t, upper) ->
+  310           ( t
+  311           , Funnel funnelLeft funnelCusp (lower F.|> t F.|> x)
+  312           , Funnel (F.singleton x) t upper)
+  313     | otherwise =
+  314       ( funnelCusp
+  315       , Funnel funnelLeft funnelCusp (F.singleton x)
+  316       , Funnel (F.singleton x) funnelCusp funnelRight)
+  317   where
+  318     isOnLeftChain  = fromMaybe False $
+  319       isLeftTurnOrLinear cuspElt <$> leftElt <*> pure targetElt
+  320     isOnRightChain = fromMaybe False $
+  321       isRightTurnOrLinear cuspElt <$> rightElt <*> pure targetElt
+  322     doSearch fn chain =
+  323       case F.search (searchChain fn) (chain::Chain) of
+  324         F.Position lower t upper -> (lower, t, upper)
+  325         F.OnLeft                 -> error "cannot happen"
+  326         F.OnRight                -> error "cannot happen"
+  327         F.Nowhere                -> error "cannot happen"
+  328     searchChain _ MinMaxEmpty _             = False
+  329     searchChain _ _ MinMaxEmpty             = True
+  330     searchChain check (MinMax _ l) (MinMax r _) =
+  331       check (ringAccess p l) (ringAccess p r) targetElt
+  332     cuspElt   = ringAccess p funnelCusp
+  333     targetElt = ringAccess p x
+  334     leftElt   = ringAccess p <$> chainLeft funnelLeft
+  335     rightElt  = ringAccess p <$> chainLeft funnelRight
+  336     chainLeft chain =
+  337       case F.viewl chain of
+  338         F.EmptyL   -> Nothing
+  339         elt F.:< _ -> Just elt
+  340 
+  341 -- O(n)
+  342 ssspFinger :: (Epsilon a, Fractional a, Ord a) => Ring a -> Dual -> SSSP
+  343 ssspFinger p d = toSSSP $
+  344     case d of
+  345       Dual (a,b,c) l r ->
+  346         (a, a) :
+  347         (b, a) :
+  348         (c, a) :
+  349         worker (Funnel (F.singleton c) a (F.singleton b)) r ++
+  350         loopLeft a c l
+  351   where
+  352     toSSSP edges =
+  353       (V.fromList . map snd . sortOn fst) edges
+  354     loopLeft a outer l =
+  355       case l of
+  356         EmptyDual -> []
+  357         NodeDual x l' r' ->
+  358           (x,a) :
+  359           worker (Funnel (F.singleton x) a (F.singleton outer)) r' ++
+  360           loopLeft a x l'
+  361     worker _ EmptyDual = []
+  362     worker f (NodeDual x l r) =
+  363       case splitFunnel p x f of
+  364         (v, fL, fR) ->
+  365           (x, v) :
+  366           worker fL l ++
+  367           worker fR r
diff --git a/reanimate-0.4.3.0-inplace/Reanimate.PolyShape.hs.html b/reanimate-0.4.3.0-inplace/Reanimate.PolyShape.hs.html index 8eb0b0d..47fc3fd 100644 --- a/reanimate-0.4.3.0-inplace/Reanimate.PolyShape.hs.html +++ b/reanimate-0.4.3.0-inplace/Reanimate.PolyShape.hs.html @@ -157,326 +157,306 @@ span.spaces { background: white } 138 then [bezierSubsegment c 0 (arcLengthParam c l polyShapeTolerance)] 139 else c : takeLen (l-cLen) cs 140 - 141 -- earClip :: Polygon -> Triangulation - 142 -- dual :: Triangulation -> Dual - 143 -- toPDual :: Polygon -> Dual -> PDual - 144 -- pdualReduce :: Polygon -> PDual -> Int -> PDual - 145 -- pdualPolygons :: Polygon -> PDual -> [Polygon] - 146 -- splitPolyShape :: Double -> Int -> PolyShape -> [PolyShape] - 147 -- splitPolyShape tol n poly = - 148 -- let polygon = toPolygon (plPolygonify tol poly) - 149 -- trig = triangulate $ pRing polygon - 150 -- d = dual 0 trig - 151 -- pd = toPDual (pRing polygon) d - 152 -- reduced = pdualReduce (pRing polygon) pd n - 153 -- polygons = pdualPolygons polygon reduced - 154 -- in map toPolyShape polygons - 155 -- where - 156 -- toPolygon :: [RPoint] -> Polygon - 157 -- toPolygon = mkPolygon . V.fromList . nub . map (fmap realToFrac) - 158 -- toPolyShape :: Polygon -> PolyShape - 159 -- toPolyShape = plFromPolygon . map (fmap realToFrac) . V.toList . polygonPoints - 160 - 161 -- plPartial' :: Double -> ([RPoint], PolyShape) -> PolyShape - 162 -- plPartial' delta (seen', PolyShape (ClosedPath lst)) = - 163 -- case lst of - 164 -- [] -> PolyShape (ClosedPath []) - 165 -- (startP, startJoin) : rest -> PolyShape $ ClosedPath $ - 166 -- (startP, startJoin) : worker startP rest - 167 -- where - 168 -- seen = filter (`elem` plPoints) seen' - 169 -- closestSeen pt = minimumBy (comparing (vectorDistance pt)) seen - 170 -- worker _ [] = [] - 171 -- worker _ ((newP, newJoin) : rest) - 172 -- | newP `elem` seen = (newP, newJoin) : worker newP rest - 173 -- | otherwise = - 174 -- let newAt = interpolateVector (closestSeen newP) newP delta - 175 -- in (newAt, newJoin) : worker newAt rest - 176 -- plPoints = - 177 -- [ p | (p,_) <- lst ] - 178 - 179 -- | Find intersection points. - 180 plGroupTouching :: [PolyShape] -> [[([RPoint],PolyShape)]] - 181 plGroupTouching [] = [] - 182 plGroupTouching pls = worker [polyShapeOrigin (head pls)] pls - 183 where - 184 worker _ [] = [] - 185 worker seen shapes = - 186 let (touching, notTouching) = partition (isTouching seen) shapes - 187 in if null touching - 188 then plGroupTouching notTouching - 189 else map ((,) seen . changeOrigin seen) touching : - 190 worker (seen ++ concatMap plPoints touching) notTouching - 191 isTouching pts = any (`elem` pts) . plPoints - 192 changeOrigin seen (PolyShape (ClosedPath segments)) = PolyShape $ ClosedPath $ helper [] segments - 193 where - 194 helper acc [] = reverse acc - 195 helper acc lst@((startP,startJ):rest) - 196 | startP `elem` seen = lst ++ reverse acc - 197 | otherwise = helper ((startP, startJ):acc) rest - 198 plPoints :: PolyShape -> [RPoint] - 199 plPoints (PolyShape (ClosedPath lst)) = - 200 [ p | (p,_) <- lst ] - 201 - 202 -- | Deconstruct a polyshape into non-intersecting, convex polygons. - 203 plDecompose :: [PolyShape] -> [[RPoint]] - 204 plDecompose = plDecompose' 0.001 - 205 - 206 -- | Deconstruct a polyshape into non-intersecting, convex polygons. - 207 plDecompose' :: Double -> [PolyShape] -> [[RPoint]] - 208 plDecompose' tol = - 209 concatMap (decomposePolygon . plPolygonify tol . mergePolyShapeHoles) . - 210 plGroupShapes . - 211 unionPolyShapes - 212 - 213 -- | Split polygon into smaller, convex polygons. - 214 decomposePolygon :: [RPoint] -> [[RPoint]] - 215 decomposePolygon poly = - 216 [ [ V2 x y - 217 | v <- V.toList (Geo.boundaryVertices f pg) - 218 , let Geo.Point2 x y =(pg^.Geo.vertexDataOf v) ^. Geo.location ] - 219 | (f, Inside) <- V.toList (Geo.internalFaces pg) ] - 220 - 221 where - 222 pg = triangulate' Proxy p - 223 p = Geo.fromPoints $ - 224 [ Geo.Point2 x y :+ () - 225 | V2 x y <- poly ] + 141 -- plPartial' :: Double -> ([RPoint], PolyShape) -> PolyShape + 142 -- plPartial' delta (seen', PolyShape (ClosedPath lst)) = + 143 -- case lst of + 144 -- [] -> PolyShape (ClosedPath []) + 145 -- (startP, startJoin) : rest -> PolyShape $ ClosedPath $ + 146 -- (startP, startJoin) : worker startP rest + 147 -- where + 148 -- seen = filter (`elem` plPoints) seen' + 149 -- closestSeen pt = minimumBy (comparing (vectorDistance pt)) seen + 150 -- worker _ [] = [] + 151 -- worker _ ((newP, newJoin) : rest) + 152 -- | newP `elem` seen = (newP, newJoin) : worker newP rest + 153 -- | otherwise = + 154 -- let newAt = interpolateVector (closestSeen newP) newP delta + 155 -- in (newAt, newJoin) : worker newAt rest + 156 -- plPoints = + 157 -- [ p | (p,_) <- lst ] + 158 + 159 -- | Find intersection points. + 160 plGroupTouching :: [PolyShape] -> [[([RPoint],PolyShape)]] + 161 plGroupTouching [] = [] + 162 plGroupTouching pls = worker [polyShapeOrigin (head pls)] pls + 163 where + 164 worker _ [] = [] + 165 worker seen shapes = + 166 let (touching, notTouching) = partition (isTouching seen) shapes + 167 in if null touching + 168 then plGroupTouching notTouching + 169 else map ((,) seen . changeOrigin seen) touching : + 170 worker (seen ++ concatMap plPoints touching) notTouching + 171 isTouching pts = any (`elem` pts) . plPoints + 172 changeOrigin seen (PolyShape (ClosedPath segments)) = PolyShape $ ClosedPath $ helper [] segments + 173 where + 174 helper acc [] = reverse acc + 175 helper acc lst@((startP,startJ):rest) + 176 | startP `elem` seen = lst ++ reverse acc + 177 | otherwise = helper ((startP, startJ):acc) rest + 178 plPoints :: PolyShape -> [RPoint] + 179 plPoints (PolyShape (ClosedPath lst)) = + 180 [ p | (p,_) <- lst ] + 181 + 182 -- | Deconstruct a polyshape into non-intersecting, convex polygons. + 183 plDecompose :: [PolyShape] -> [[RPoint]] + 184 plDecompose = plDecompose' 0.001 + 185 + 186 -- | Deconstruct a polyshape into non-intersecting, convex polygons. + 187 plDecompose' :: Double -> [PolyShape] -> [[RPoint]] + 188 plDecompose' tol = + 189 concatMap (decomposePolygon . plPolygonify tol . mergePolyShapeHoles) . + 190 plGroupShapes . + 191 unionPolyShapes + 192 + 193 -- | Split polygon into smaller, convex polygons. + 194 decomposePolygon :: [RPoint] -> [[RPoint]] + 195 decomposePolygon poly = + 196 [ [ V2 x y + 197 | v <- V.toList (Geo.boundaryVertices f pg) + 198 , let Geo.Point2 x y =(pg^.Geo.vertexDataOf v) ^. Geo.location ] + 199 | (f, Inside) <- V.toList (Geo.internalFaces pg) ] + 200 + 201 where + 202 pg = triangulate' Proxy p + 203 p = Geo.fromPoints $ + 204 [ Geo.Point2 x y :+ () + 205 | V2 x y <- poly ] + 206 + 207 plPolygonify :: Double -> PolyShape -> [RPoint] + 208 plPolygonify tol shape = + 209 startPoint (head curves) : concatMap worker curves + 210 where + 211 curves = plCurves shape + 212 worker c | endPoint c == startPoint c = + 213 [] -- error $ "Bad bezier: " ++ show c + 214 worker c = + 215 if colinear c tol -- && arcLength c 1 tol < 1 + 216 then [endPoint c] + 217 else + 218 let (lhs,rhs) = splitBezier c 0.5 + 219 in worker lhs ++ worker rhs + 220 endPoint (CubicBezier _ _ _ d) = d + 221 startPoint (CubicBezier a _ _ _) = a + 222 + 223 -- | Convert a polyshape to a list of SVG path commands. + 224 plPathCommands :: PolyShape -> [PathCommand] + 225 plPathCommands = lineToPath . plLineCommands 226 - 227 plPolygonify :: Double -> PolyShape -> [RPoint] - 228 plPolygonify tol shape = - 229 startPoint (head curves) : concatMap worker curves - 230 where - 231 curves = plCurves shape - 232 worker c | endPoint c == startPoint c = - 233 [] -- error $ "Bad bezier: " ++ show c - 234 worker c = - 235 if colinear c tol -- && arcLength c 1 tol < 1 - 236 then [endPoint c] - 237 else - 238 let (lhs,rhs) = splitBezier c 0.5 - 239 in worker lhs ++ worker rhs - 240 endPoint (CubicBezier _ _ _ d) = d - 241 startPoint (CubicBezier a _ _ _) = a - 242 - 243 -- | Convert a polyshape to a list of SVG path commands. - 244 plPathCommands :: PolyShape -> [PathCommand] - 245 plPathCommands = lineToPath . plLineCommands - 246 - 247 -- | Convert a polyshape to a list of line commands. - 248 plLineCommands :: PolyShape -> [LineCommand] - 249 plLineCommands pl = - 250 case curves of - 251 [] -> [] - 252 (CubicBezier start _ _ _:_) -> - 253 LineMove start : - 254 zipWith worker (drop 1 dstList ++ [start]) joinList ++ - 255 [LineEnd start] - 256 where - 257 ClosedPath closedPath = unPolyShape pl - 258 (dstList, joinList) = unzip closedPath - 259 curves = plCurves pl - 260 worker dst JoinLine = - 261 LineBezier [dst] - 262 worker dst (JoinCurve a b) = - 263 LineBezier [a,b,dst] - 264 - 265 -- | Extract all shapes from SVG nodes. Drawing attributes such - 266 -- as stroke and fill color are discarded. - 267 svgToPolyShapes :: Tree -> [PolyShape] - 268 svgToPolyShapes = cmdsToPolyShapes . toLineCommands . extractPath - 269 - 270 -- | Extract all polygons from SVG nodes. Curves are approximated to - 271 -- within the given tolerance. - 272 svgToPolygons :: Double -> SVG -> [Polygon] - 273 svgToPolygons tol = map (toPolygon . plPolygonify tol) . svgToPolyShapes - 274 where - 275 toPolygon :: [RPoint] -> Polygon - 276 toPolygon = mkPolygon . - 277 V.fromList . nub . map (fmap realToFrac) - 278 - 279 cmdsToPolyShapes :: [LineCommand] -> [PolyShape] - 280 cmdsToPolyShapes [] = [] - 281 cmdsToPolyShapes cmds = - 282 case cmds of - 283 (LineMove dst:cont) -> map PolyShape $ worker dst [] cont - 284 _ -> bad - 285 where - 286 bad = error $ "Reanimate.PolyShape: Invalid commands: " ++ show cmds - 287 finalize [] rest = rest - 288 finalize acc rest = ClosedPath (reverse acc) : rest - 289 worker _from acc [] = finalize acc [] - 290 worker _from acc (LineMove newStart : xs) = - 291 finalize acc $ - 292 worker newStart [] xs - 293 worker from acc (LineEnd orig:LineMove dst:xs) | from /= orig = - 294 finalize ((from, JoinLine):acc) $ - 295 worker dst [] xs - 296 worker _from acc (LineEnd{}:LineMove dst:xs) = - 297 finalize acc $ - 298 worker dst [] xs - 299 worker from acc [LineEnd orig] | from /= orig = - 300 finalize ((from, JoinLine):acc) [] - 301 worker _from acc [LineEnd{}] = - 302 finalize acc [] - 303 worker from acc (LineBezier [x]:xs) = - 304 worker x ((from, JoinLine) : acc) xs - 305 worker from acc (LineBezier [a,b]:xs) = - 306 let quad = QuadBezier from a b - 307 CubicBezier _ a' b' c' = quadToCubic quad - 308 in worker from acc (LineBezier [a',b',c']:xs) - 309 worker from acc (LineBezier [a,b,c]:xs) = - 310 worker c ((from, JoinCurve a b) : acc) xs - 311 worker _ _ _ = bad - 312 - 313 -- | Merge overlapping shapes. - 314 unionPolyShapes :: [PolyShape] -> [PolyShape] - 315 unionPolyShapes shapes = - 316 map PolyShape $ - 317 union (map unPolyShape shapes) FillNonZero (polyShapeTolerance/10000) - 318 - 319 -- | Merge overlapping shapes to within given tolerance. - 320 unionPolyShapes' :: Double -> [PolyShape] -> [PolyShape] - 321 unionPolyShapes' tol shapes = - 322 map PolyShape $ - 323 union (map unPolyShape shapes) FillNonZero tol - 324 - 325 -- | True iff lhs is inside of rhs. - 326 -- lhs and rhs may not overlap. - 327 -- Implementation: Trace a vertical line through the origin of A and check - 328 -- of this line intersects and odd number of times on both sides of A. - 329 isInsideOf :: PolyShape -> PolyShape -> Bool - 330 lhs `isInsideOf` rhs = - 331 odd (length upHits) && odd (length downHits) - 332 where - 333 (upHits, downHits) = polyIntersections origin rhs - 334 origin = polyShapeOrigin lhs - 335 - 336 polyIntersections :: RPoint -> PolyShape -> ([RPoint],[RPoint]) - 337 polyIntersections origin rhs = - 338 (nub $ concatMap (intersections rayUp) curves - 339 ,nub $ concatMap (intersections rayDown) curves) - 340 where - 341 curves = plCurves rhs - 342 - 343 intersections line bs = - 344 map (evalBezier bs . fst) (bezierIntersection bs line polyShapeTolerance) - 345 limit = 1000 - 346 rayUp = CubicBezier origin origin origin (V2 limit limit) - 347 rayDown = CubicBezier origin origin origin (V2 (-limit) (-limit)) - 348 - 349 polyShapeOrigin :: PolyShape -> V2 Double - 350 polyShapeOrigin (PolyShape closedPath) = - 351 case closedPath of - 352 ClosedPath [] -> V2 0 0 - 353 ClosedPath ((start,_):_) -> start - 354 - 355 -- | Find holes and group them with their parent. - 356 plGroupShapes :: [PolyShape] -> [PolyShapeWithHoles] - 357 plGroupShapes = worker - 358 where - 359 worker (s:rest) - 360 | null (parents s rest) = - 361 let isOnlyChild x = parents x (s:rest) == [s] - 362 (holes, nonHoles) = partition isOnlyChild rest - 363 prime = PolyShapeWithHoles - 364 { polyShapeParent = s - 365 , polyShapeHoles = holes } - 366 in prime : worker nonHoles - 367 | otherwise = worker (rest ++ [s]) - 368 worker [] = [] - 369 - 370 parents :: PolyShape -> [PolyShape] -> [PolyShape] - 371 parents self = filter (self `isInsideOf`) . filter (/=self) - 372 - 373 instance Eq PolyShape where - 374 a == b = plCurves a == plCurves b - 375 - 376 -- | Cut out holes. - 377 mergePolyShapeHoles :: PolyShapeWithHoles -> PolyShape - 378 mergePolyShapeHoles (PolyShapeWithHoles parent []) = parent - 379 mergePolyShapeHoles (PolyShapeWithHoles parent (child:children)) = - 380 mergePolyShapeHoles $ - 381 PolyShapeWithHoles (mergePolyShapeHole parent child) children - 382 - 383 -- Merge - 384 mergePolyShapeHole :: PolyShape -> PolyShape -> PolyShape - 385 mergePolyShapeHole parent child = - 386 snd $ head $ - 387 sortOn fst - 388 [ cutSingleHole newParent child - 389 | newParent <- polyShapePermutations parent ] + 227 -- | Convert a polyshape to a list of line commands. + 228 plLineCommands :: PolyShape -> [LineCommand] + 229 plLineCommands pl = + 230 case curves of + 231 [] -> [] + 232 (CubicBezier start _ _ _:_) -> + 233 LineMove start : + 234 zipWith worker (drop 1 dstList ++ [start]) joinList ++ + 235 [LineEnd start] + 236 where + 237 ClosedPath closedPath = unPolyShape pl + 238 (dstList, joinList) = unzip closedPath + 239 curves = plCurves pl + 240 worker dst JoinLine = + 241 LineBezier [dst] + 242 worker dst (JoinCurve a b) = + 243 LineBezier [a,b,dst] + 244 + 245 -- | Extract all shapes from SVG nodes. Drawing attributes such + 246 -- as stroke and fill color are discarded. + 247 svgToPolyShapes :: Tree -> [PolyShape] + 248 svgToPolyShapes = cmdsToPolyShapes . toLineCommands . extractPath + 249 + 250 -- | Extract all polygons from SVG nodes. Curves are approximated to + 251 -- within the given tolerance. + 252 svgToPolygons :: Double -> SVG -> [Polygon] + 253 svgToPolygons tol = map (toPolygon . plPolygonify tol) . svgToPolyShapes + 254 where + 255 toPolygon :: [RPoint] -> Polygon + 256 toPolygon = mkPolygon . + 257 V.fromList . nub . map (fmap realToFrac) + 258 + 259 cmdsToPolyShapes :: [LineCommand] -> [PolyShape] + 260 cmdsToPolyShapes [] = [] + 261 cmdsToPolyShapes cmds = + 262 case cmds of + 263 (LineMove dst:cont) -> map PolyShape $ worker dst [] cont + 264 _ -> bad + 265 where + 266 bad = error $ "Reanimate.PolyShape: Invalid commands: " ++ show cmds + 267 finalize [] rest = rest + 268 finalize acc rest = ClosedPath (reverse acc) : rest + 269 worker _from acc [] = finalize acc [] + 270 worker _from acc (LineMove newStart : xs) = + 271 finalize acc $ + 272 worker newStart [] xs + 273 worker from acc (LineEnd orig:LineMove dst:xs) | from /= orig = + 274 finalize ((from, JoinLine):acc) $ + 275 worker dst [] xs + 276 worker _from acc (LineEnd{}:LineMove dst:xs) = + 277 finalize acc $ + 278 worker dst [] xs + 279 worker from acc [LineEnd orig] | from /= orig = + 280 finalize ((from, JoinLine):acc) [] + 281 worker _from acc [LineEnd{}] = + 282 finalize acc [] + 283 worker from acc (LineBezier [x]:xs) = + 284 worker x ((from, JoinLine) : acc) xs + 285 worker from acc (LineBezier [a,b]:xs) = + 286 let quad = QuadBezier from a b + 287 CubicBezier _ a' b' c' = quadToCubic quad + 288 in worker from acc (LineBezier [a',b',c']:xs) + 289 worker from acc (LineBezier [a,b,c]:xs) = + 290 worker c ((from, JoinCurve a b) : acc) xs + 291 worker _ _ _ = bad + 292 + 293 -- | Merge overlapping shapes. + 294 unionPolyShapes :: [PolyShape] -> [PolyShape] + 295 unionPolyShapes shapes = + 296 map PolyShape $ + 297 union (map unPolyShape shapes) FillNonZero (polyShapeTolerance/10000) + 298 + 299 -- | Merge overlapping shapes to within given tolerance. + 300 unionPolyShapes' :: Double -> [PolyShape] -> [PolyShape] + 301 unionPolyShapes' tol shapes = + 302 map PolyShape $ + 303 union (map unPolyShape shapes) FillNonZero tol + 304 + 305 -- | True iff lhs is inside of rhs. + 306 -- lhs and rhs may not overlap. + 307 -- Implementation: Trace a vertical line through the origin of A and check + 308 -- of this line intersects and odd number of times on both sides of A. + 309 isInsideOf :: PolyShape -> PolyShape -> Bool + 310 lhs `isInsideOf` rhs = + 311 odd (length upHits) && odd (length downHits) + 312 where + 313 (upHits, downHits) = polyIntersections origin rhs + 314 origin = polyShapeOrigin lhs + 315 + 316 polyIntersections :: RPoint -> PolyShape -> ([RPoint],[RPoint]) + 317 polyIntersections origin rhs = + 318 (nub $ concatMap (intersections rayUp) curves + 319 ,nub $ concatMap (intersections rayDown) curves) + 320 where + 321 curves = plCurves rhs + 322 + 323 intersections line bs = + 324 map (evalBezier bs . fst) (bezierIntersection bs line polyShapeTolerance) + 325 limit = 1000 + 326 rayUp = CubicBezier origin origin origin (V2 limit limit) + 327 rayDown = CubicBezier origin origin origin (V2 (-limit) (-limit)) + 328 + 329 polyShapeOrigin :: PolyShape -> V2 Double + 330 polyShapeOrigin (PolyShape closedPath) = + 331 case closedPath of + 332 ClosedPath [] -> V2 0 0 + 333 ClosedPath ((start,_):_) -> start + 334 + 335 -- | Find holes and group them with their parent. + 336 plGroupShapes :: [PolyShape] -> [PolyShapeWithHoles] + 337 plGroupShapes = worker + 338 where + 339 worker (s:rest) + 340 | null (parents s rest) = + 341 let isOnlyChild x = parents x (s:rest) == [s] + 342 (holes, nonHoles) = partition isOnlyChild rest + 343 prime = PolyShapeWithHoles + 344 { polyShapeParent = s + 345 , polyShapeHoles = holes } + 346 in prime : worker nonHoles + 347 | otherwise = worker (rest ++ [s]) + 348 worker [] = [] + 349 + 350 parents :: PolyShape -> [PolyShape] -> [PolyShape] + 351 parents self = filter (self `isInsideOf`) . filter (/=self) + 352 + 353 instance Eq PolyShape where + 354 a == b = plCurves a == plCurves b + 355 + 356 -- | Cut out holes. + 357 mergePolyShapeHoles :: PolyShapeWithHoles -> PolyShape + 358 mergePolyShapeHoles (PolyShapeWithHoles parent []) = parent + 359 mergePolyShapeHoles (PolyShapeWithHoles parent (child:children)) = + 360 mergePolyShapeHoles $ + 361 PolyShapeWithHoles (mergePolyShapeHole parent child) children + 362 + 363 -- Merge + 364 mergePolyShapeHole :: PolyShape -> PolyShape -> PolyShape + 365 mergePolyShapeHole parent child = + 366 snd $ head $ + 367 sortOn fst + 368 [ cutSingleHole newParent child + 369 | newParent <- polyShapePermutations parent ] + 370 + 371 {- + 372 parent: + 373 (a,b) + 374 (b,c) + 375 (c,a) + 376 + 377 child: + 378 (x,y) + 379 (y,z) + 380 (z,x) + 381 + 382 P = split (a,b) + 383 new: + 384 (P,b) p2b + 385 (b,c) pTail + 386 (c,a) pTail + 387 (a,P) a2p + 388 + 389 (P,x) p2x 390 - 391 {- - 392 parent: - 393 (a,b) - 394 (b,c) - 395 (c,a) + 391 (x,y) childCurves + 392 (y,z) childCurves + 393 (z,x) childCurves + 394 + 395 (x,P) x2p 396 - 397 child: - 398 (x,y) - 399 (y,z) - 400 (z,x) - 401 - 402 P = split (a,b) - 403 new: - 404 (P,b) p2b - 405 (b,c) pTail - 406 (c,a) pTail - 407 (a,P) a2p - 408 - 409 (P,x) p2x - 410 - 411 (x,y) childCurves - 412 (y,z) childCurves - 413 (z,x) childCurves - 414 - 415 (x,P) x2p - 416 - 417 -} - 418 cutSingleHole :: PolyShape -> PolyShape -> (Double, PolyShape) - 419 cutSingleHole parent child = - 420 (score, PolyShape $ curvesToClosed $ - 421 p2b:pTail ++ [a2p] ++ - 422 [p2x] ++ childCurves ++ - 423 [x2p] - 424 ) - 425 where - 426 -- vect = (childOrigin - p) * 0 -- 0.0001 - 427 vectL = 0 -- rotate90L $* vect - 428 vectR = 0 -- rotate90R $* vect - 429 score = vectorDistance childOrigin p - 430 childOrigin = polyShapeOrigin child - 431 childOrigin' = childOrigin - vectL - 432 (pHead:pTail) = plCurves parent - 433 childCurves = plCurves child - 434 - 435 pParam = closest pHead childOrigin polyShapeTolerance - 436 - 437 (a2p, p2b') = splitBezier pHead pParam - 438 p2b = case p2b' of - 439 CubicBezier a b c d -> CubicBezier (a - vectL) b c d - 440 - 441 p = evalBezier pHead pParam - 442 -- straight line to child origin - 443 p2x = lineBetween (p - vectR) childOrigin - 444 -- straight line from child origin - 445 x2p = lineBetween childOrigin' p - 446 - 447 lineBetween a = CubicBezier a a a - 448 - 449 -- | Destruct a polyshape into constituent curves. - 450 plCurves :: PolyShape -> [CubicBezier Double] - 451 plCurves = closedPathCurves . unPolyShape - 452 - 453 polyShapePermutations :: PolyShape -> [PolyShape] - 454 polyShapePermutations = - 455 map (PolyShape . curvesToClosed) . cycleList . plCurves - 456 where - 457 cycleList lst = - 458 let n = length lst in - 459 [ take n $ drop i $ cycle lst - 460 | i <- [0.. n-1] ] + 397 -} + 398 cutSingleHole :: PolyShape -> PolyShape -> (Double, PolyShape) + 399 cutSingleHole parent child = + 400 (score, PolyShape $ curvesToClosed $ + 401 p2b:pTail ++ [a2p] ++ + 402 [p2x] ++ childCurves ++ + 403 [x2p] + 404 ) + 405 where + 406 -- vect = (childOrigin - p) * 0 -- 0.0001 + 407 vectL = 0 -- rotate90L $* vect + 408 vectR = 0 -- rotate90R $* vect + 409 score = vectorDistance childOrigin p + 410 childOrigin = polyShapeOrigin child + 411 childOrigin' = childOrigin - vectL + 412 (pHead:pTail) = plCurves parent + 413 childCurves = plCurves child + 414 + 415 pParam = closest pHead childOrigin polyShapeTolerance + 416 + 417 (a2p, p2b') = splitBezier pHead pParam + 418 p2b = case p2b' of + 419 CubicBezier a b c d -> CubicBezier (a - vectL) b c d + 420 + 421 p = evalBezier pHead pParam + 422 -- straight line to child origin + 423 p2x = lineBetween (p - vectR) childOrigin + 424 -- straight line from child origin + 425 x2p = lineBetween childOrigin' p + 426 + 427 lineBetween a = CubicBezier a a a + 428 + 429 -- | Destruct a polyshape into constituent curves. + 430 plCurves :: PolyShape -> [CubicBezier Double] + 431 plCurves = closedPathCurves . unPolyShape + 432 + 433 polyShapePermutations :: PolyShape -> [PolyShape] + 434 polyShapePermutations = + 435 map (PolyShape . curvesToClosed) . cycleList . plCurves + 436 where + 437 cycleList lst = + 438 let n = length lst in + 439 [ take n $ drop i $ cycle lst + 440 | i <- [0.. n-1] ]