never executed always true always false
1 {-# LANGUAGE BangPatterns #-}
2 {-# LANGUAGE ConstraintKinds #-}
3 {-# OPTIONS_HADDOCK hide #-}
4 module Reanimate.Math.Polygon
5 ( APolygon(..)
6 , Polygon
7 , FPolygon
8 , P
9 , mkPolygon -- :: (Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a
10 , mkPolygonFromRing -- :: (Fractional a, Ord a) => Ring a -> APolygon a
11 , castPolygon -- :: (Real a, Fractional b, Ord a) => APolygon a -> APolygon b
12 , pParent -- :: Polygon -> Int -> Int -> Int
13 , pSetOffset -- :: APolygon a -> Int -> APolygon a
14 , pAdjustOffset -- :: APolygon a -> Int -> APolygon a
15 , pSize -- :: APolygon a -> Int
16 , pNull -- :: APolygon a -> Bool
17 , pNext -- :: APolygon a -> Int -> Int
18 , pPrev -- :: APolygon a -> Int -> Int
19 , pIsSimple -- :: Polygon -> Bool
20 , pIsConvex -- :: Polygon -> Bool
21 , pIsCCW -- :: Polygon -> Bool
22 , pScale -- :: Rational -> Polygon -> Polygon
23 , pAtCentroid -- :: Polygon -> Polygon
24 , pAtCenter -- :: Polygon -> Polygon
25 , pTranslate -- :: V2 Rational -> Polygon -> Polygon
26 , pCenter -- :: Polygon -> V2 Rational
27 , pBoundingBox -- :: Polygon -> (Rational, Rational, Rational, Rational)
28 , pIsInside -- :: Polygon -> V2 Rational -> Bool
29 , pAccess -- :: APolygon a -> Int -> V2 a
30 , pMkWinding -- :: Int -> Polygon
31 , pDeoverlap -- :: Polygon -> Polygon
32 , pCycles -- :: Polygon -> [Polygon]
33 , pCycle -- :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a
34 , pCentroid -- :: Polygon -> V2 Rational
35 , pMapEdges -- :: (V2 Rational -> V2 Rational -> a) -> Polygon -> V.Vector a
36 , pArea -- :: Polygon -> Rational
37 , pCircumference -- :: (Real a, Fractional a) => APolygon a -> a
38 , pCircumference' -- :: (Real a, Fractional a) => APolygon a -> Double
39 , pAddPoints -- :: Int -> Polygon -> Polygon
40 , pAddPointsRestricted -- :: [Int] -> Int -> Polygon -> Polygon
41 , pAddPointsBetween -- :: (Fractional a, Ord a, Real a) => (Int, Int) -> Int -> APolygon a -> APolygon a
42 , pRayIntersect -- :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
43 , pOverlap -- :: Polygon -> Polygon -> Polygon
44 , pCuts -- :: Polygon -> [(Polygon,Polygon)]
45 , pCutEqual -- :: Polygon -> (Polygon, Polygon)
46 -- * Triangulation
47 , isValidTriangulation -- :: Polygon -> Triangulation -> Bool
48 , triangulationsToPolygons -- :: Polygon -> Triangulation -> [Polygon]
49 -- * Single-Source-Shortest-Path
50 , ssspVisibility -- :: Polygon -> Polygon
51 , ssspWindows -- :: Polygon -> [(V2 Rational, V2 Rational)]
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 -- 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 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 = pAddPointsRestricted []
589
590 pAddPointsRestricted :: PolyCtx a => [(V2 a, V2 a)] -> Int -> APolygon a -> APolygon a
591 pAddPointsRestricted _immutableEdges n p | n <= 0 = p
592 pAddPointsRestricted immutableEdges n p = pAddPointsRestricted immutableEdges (n-1) $
593 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1]
594 where
595 isImmutable idx =
596 (pAccess p idx, pAccess p $ idx+1) `elem` immutableEdges ||
597 (pAccess p $ idx+1, pAccess p idx) `elem` immutableEdges
598 worker idx
599 | idx == longestEdge && not (isImmutable idx) =
600 [pAccess p idx, pMiddlePoint p idx]
601 | otherwise = [pAccess p idx]
602 longestEdge = maximumBy cmpLength [0 .. pSize p-1]
603 cmpLength a _ | isImmutable a = LT
604 cmpLength _ b | isImmutable b = GT
605 cmpLength a b =
606 distSquared (pAccess p a) (pAccess p $ a+1) `compare`
607 distSquared (pAccess p b) (pAccess p $ b+1)
608
609 pMiddlePoint :: PolyCtx a => APolygon a -> Int -> V2 a
610 pMiddlePoint p idx
611 = lerp 0.5 (pAccess p $ idx+1) (pAccess p idx)
612
613 pAddPointsBetween :: PolyCtx a => (Int, Int) -> Int -> APolygon a -> APolygon a
614 pAddPointsBetween _ n p | n <= 0 = p
615 pAddPointsBetween (i,l) n p = pAddPointsBetween (i,l+1) (n-1) $
616 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1]
617 where
618 worker idx
619 | idx == longestEdge =
620 [pAccess p idx, pMiddlePoint p idx]
621 | otherwise = [pAccess p idx]
622 longestEdge = maximumBy cmpLength [i .. i+l-1]
623 cmpLength a b =
624 distSquared (pAccess p a) (pAccess p $ a+1) `compare`
625 distSquared (pAccess p b) (pAccess p $ b+1)
626
627 -- addPoints :: Int -> Polygon -> Polygon
628 -- addPoints n p = mkPolygon $ V.fromList $ worker n 0 (map (pAccess p) [0..s])
629 -- where
630 -- worker 0 _ rest = init rest
631 -- worker i acc (x:y:xs) =
632 -- let xy = approxDist x y in
633 -- if acc + xy > limit
634 -- then x : worker (i-1) 0 (lerp ((limit-acc)/xy) y x : y:xs)
635 -- else x : worker i (acc+xy) (y:xs)
636 -- worker _ _ [_] = []
637 -- worker _ _ _ = error "addPoints: invalid polygon"
638 -- s = pSize p
639 -- len = polygonLength p
640 -- limit = len / fromIntegral (n+1)
641
642 pIsConvex :: Polygon -> Bool
643 pIsConvex p = and
644 [ area2X (pAccess p i) (pAccess p j) (pAccess p k) > 0
645 | i <- [0..n-1]
646 , j <- [i+1..n-1]
647 , k <- [j+1..n-1]
648 ]
649 where n = pSize p
650
651 pIsCCW :: Polygon -> Bool
652 pIsCCW p | pNull p = False
653 pIsCCW p = V.sum (pMapEdges fn p) < 0
654 where
655 fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1)
656
657 {-# INLINE pRayIntersect #-}
658 pRayIntersect :: PolyCtx a => APolygon a -> (Int, Int) -> (Int,Int) -> Maybe (V2 a)
659 pRayIntersect p (a,b) (c,d) =
660 rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d)
661
662 pCuts :: (Real a, Fractional a, Epsilon a) => APolygon a -> [(APolygon a,APolygon a)]
663 pCuts p =
664 [ pCutAt (pAdjustOffset p i) (j-i)
665 | i <- [0 .. pSize p-1 ]
666 , j <- [i+2 .. pSize p-1 ]
667 , (j+1) `mod` pSize p /= i
668 , pParent p i j == i ]
669
670 pCutEqual :: PolyCtx a => APolygon a -> (APolygon a, APolygon a)
671 pCutEqual p =
672 fromMaybe (p,p) $ listToMaybe $ sortOn f $ pCuts p
673 where
674 f (a,b) = abs (pArea a - pArea b)
675
676 -- FIXME: This should be more efficient
677 pCutAt :: PolyCtx a => APolygon a -> Int -> (APolygon a, APolygon a)
678 pCutAt p i = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right)
679 where
680 n = pSize p
681 left = map (pAccess p) [0 .. i]
682 right = map (pAccess p) (0:[i..n-1])
683
684 pOverlap :: PolyCtx a => APolygon a -> APolygon a -> APolygon a
685 pOverlap a b = mkPolygon $ V.fromList $ clearDups $ concatMap edgeIntersect [0 .. pSize a-1]
686 where
687 clearDups (x:y:xs)
688 | x == y = clearDups (y:xs)
689 | otherwise = x : clearDups (y:xs)
690 clearDups xs = xs
691 edgeIntersect edge =
692 sortOn (distSquared (pAccess a edge)) $ catMaybes
693 [ lineIntersect (aP, aP') (bP, bP')
694 | i <- [0 .. pSize b-1]
695 , let aP = pAccess a edge
696 aP' = pAccess a (edge+1)
697 bP = pAccess b i
698 bP' = pAccess b (i+1)
699 ]
700
701 ---------------------------------------------------------
702 -- SSSP visibility and SSSP windows
703
704 ssspVisibility :: PolyCtx a => APolygon a -> APolygon a
705 ssspVisibility p = mkPolygon $
706 V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1])
707 where
708 clearDups (x:y:xs)
709 | x == y = clearDups (y:xs)
710 | otherwise = x : clearDups (y:xs)
711 clearDups xs = xs
712 obstructedBy n =
713 case pParent p 0 n of
714 0 -> n
715 i -> obstructedBy i
716 go [] = []
717 go [x] = [pAccess p x]
718 go (x:y:xs) =
719 let xO = obstructedBy x
720 yO = obstructedBy y
721 in case () of
722 ()
723 -- Both ends are visible.
724 | xO == x && yO == y -> pAccess p x : go (y:xs)
725 -- X is visible, x to intersect (0,yO) (x,y)
726 | xO == x ->
727 pAccess p x : fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y)) : go (y:xs)
728 -- Y is visible
729 | yO == y -> fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) : pAccess p y : go (y:xs)
730 -- Neither is visible and they've obstructed by the same point
731 -- so the entire edge is hidden.
732 | xO == yO -> go (y:xs)
733 -- Neither is visible. Cast shadow from obstruction points to
734 -- find if a subsection of the edge is visible.
735 | otherwise ->
736 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
737 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
738 in if a /= b
739 then a : b : go (y:xs)
740 else go (y:xs)
741
742 ssspWindows :: Polygon -> [(V2 Rational, V2 Rational)]
743 ssspWindows p = clearDups $ go (pAccess p 0) [0..pSize p-1]
744 where
745 clearDups (x:y:xs)
746 | x == y = clearDups (y:xs)
747 | otherwise = x : clearDups (y:xs)
748 clearDups xs = xs
749 obstructedBy n =
750 case pParent p 0 n of
751 0 -> n
752 i -> obstructedBy i
753 go _ [] = []
754 go _ [_] = []
755 go l (x:y:xs) =
756 let xO = obstructedBy x
757 yO = obstructedBy y
758 in case () of
759 ()
760 -- Both ends are visible.
761 | xO == x && yO == y -> go (pAccess p x) (y:xs)
762 -- X is visible, x to intersect (0,yO) (x,y)
763 | xO == x ->
764 go (fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y))) (y:xs)
765 -- Y is visible
766 | yO == y ->
767 let newL = fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) in
768 (l, newL) :
769 go newL (y:xs)
770 -- Neither is visible and they've obstructed by the same point
771 -- so the entire edge is hidden.
772 | xO == yO -> go l (y:xs)
773 -- Neither is visible. Cast shadow from obstruction points to
774 -- find if a subsection of the edge is visible.
775 | otherwise ->
776 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
777 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
778 in if a /= b
779 then (l, a) : (b, pAccess p yO) : go (pAccess p yO) (y:xs)
780 else go l (y:xs)