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