never executed always true always false
1 {-# LANGUAGE BangPatterns #-}
2 module Reanimate.Math.Polygon
3 ( APolygon(..)
4 , Polygon
5 , FPolygon
6 , P
7 , mkPolygon -- :: (Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a
8 , mkPolygonFromRing -- :: (Fractional a, Ord a) => Ring a -> APolygon a
9 , castPolygon -- :: (Real a, Fractional b, Ord a) => APolygon a -> APolygon b
10 , pParent -- :: Polygon -> Int -> Int -> Int
11 , pSetOffset -- :: APolygon a -> Int -> APolygon a
12 , pAdjustOffset -- :: APolygon a -> Int -> APolygon a
13 , pSize -- :: APolygon a -> Int
14 , pNull -- :: APolygon a -> Bool
15 , pNext -- :: APolygon a -> Int -> Int
16 , pPrev -- :: APolygon a -> Int -> Int
17 , pIsSimple -- :: Polygon -> Bool
18 , pIsConvex -- :: Polygon -> Bool
19 , pIsCCW -- :: Polygon -> Bool
20 , pScale -- :: Rational -> Polygon -> Polygon
21 , pAtCentroid -- :: Polygon -> Polygon
22 , pAtCenter -- :: Polygon -> Polygon
23 , pTranslate -- :: V2 Rational -> Polygon -> Polygon
24 , pCenter -- :: Polygon -> V2 Rational
25 , pBoundingBox -- :: Polygon -> (Rational, Rational, Rational, Rational)
26 , pIsInside -- :: Polygon -> V2 Rational -> Bool
27 , pAccess -- :: APolygon a -> Int -> V2 a
28 , pMkWinding -- :: Int -> Polygon
29 , pDeoverlap -- :: Polygon -> Polygon
30 , pCycles -- :: Polygon -> [Polygon]
31 , pCycle -- :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a
32 , pCentroid -- :: Polygon -> V2 Rational
33 , pMapEdges -- :: (V2 Rational -> V2 Rational -> a) -> Polygon -> V.Vector a
34 , pArea -- :: Polygon -> Rational
35 , pCircumference -- :: (Real a, Fractional a) => APolygon a -> a
36 , pCircumference' -- :: (Real a, Fractional a) => APolygon a -> Double
37 , pAddPoints -- :: Int -> Polygon -> Polygon
38 , pAddPointsRestricted -- :: [Int] -> Int -> Polygon -> Polygon
39 , pAddPointsBetween -- :: (Fractional a, Ord a, Real a) => (Int, Int) -> Int -> APolygon a -> APolygon a
40 , pRayIntersect -- :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
41 , pOverlap -- :: Polygon -> Polygon -> Polygon
42 , pCuts -- :: Polygon -> [(Polygon,Polygon)]
43 , pCutEqual -- :: Polygon -> (Polygon, Polygon)
44 -- * Triangulation
45 , isValidTriangulation -- :: Polygon -> Triangulation -> Bool
46 , triangulationsToPolygons -- :: Polygon -> Triangulation -> [Polygon]
47 -- * Single-Source-Shortest-Path
48 , ssspVisibility -- :: Polygon -> Polygon
49 , ssspWindows -- :: Polygon -> [(V2 Rational, V2 Rational)]
50 -- * Duals
51 , pdualPolygons -- :: Polygon -> PDual -> [Polygon]
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 ) where
86
87 -- import Control.Exception
88 import Data.Hashable
89 import Data.List (intersect, maximumBy, sort, sortOn,
90 tails)
91 import Data.Maybe
92 import Data.Ratio
93 import Data.Serialize
94 import Data.Vector (Vector)
95 import qualified Data.Vector as V
96 import Linear.V2
97 import Linear.Vector
98 import Reanimate.Math.Common
99 import Reanimate.Math.EarClip
100 import Reanimate.Math.SSSP
101 import Reanimate.Math.Triangulate
102
103 -- import Debug.Trace
104
105 -- Generate random polygons, options:
106 -- 1. put corners around a circle. Vary the radius.
107 -- 2. close a hilbert curve
108 type FPolygon = APolygon Double
109 -- Optimize representation?
110 -- Polygon = (Vector XNumerator, Vector XDenominator
111 -- ,Vector YNumerator, Vector YDenominator)
112 data APolygon a = Polygon
113 { polygonPoints :: Vector (V2 a)
114 , polygonOffset :: Int
115 , polygonTriangulation :: Triangulation
116 , polygonSSSP :: Vector SSSP
117 }
118 type Polygon = APolygon Rational
119 type P = V2 Double
120
121 instance Show a => Show (APolygon a) where
122 show = show . V.toList . polygonPoints
123
124 instance Hashable a => Hashable (APolygon a) where
125 hashWithSalt s p = V.foldl' hashWithSalt s (polygonPoints p)
126
127 instance (Real a, Fractional a, Serialize a) => Serialize (APolygon a) where
128 put = put . V.toList . polygonPoints
129 get = mkPolygon . V.fromList <$> get
130
131 pRing :: APolygon a -> Ring a
132 pRing = ringPack . polygonPoints
133
134 mkPolygon :: (Real a, Fractional a) => V.Vector (V2 a) -> APolygon a
135 mkPolygon points = Polygon
136 { polygonPoints = points
137 , polygonOffset = 0
138 , polygonTriangulation = trig
139 , polygonSSSP = V.generate n $ \i -> sssp ring (dual i trig)
140 }
141 where
142 n = length points
143 ring = ringPack points
144 trig = earCut ring
145 -- earClip ring
146
147 castPolygon :: (Real a, Real b, Fractional b) => APolygon a -> APolygon b
148 castPolygon = mkPolygon . V.map (fmap realToFrac) . polygonPoints
149
150 mkPolygonFromRing :: (Real a, Fractional a, Ord a) => Ring a -> APolygon a
151 mkPolygonFromRing = mkPolygon . ringUnpack
152
153 pUnsafeMap :: (Ring a -> Ring a) -> APolygon a -> APolygon a
154 pUnsafeMap fn p = p{ polygonPoints = ringUnpack (fn (pRing p)) }
155
156 -- pParent p i j = shortest-path parent from j to i
157 pParent :: APolygon a -> Int -> Int -> Int
158 pParent p i j =
159 (sTree V.! mod (j + polygonOffset p) n - polygonOffset p) `mod` n
160 where
161 sTree = polygonSSSP p V.! mod (i + polygonOffset p) n
162 n = pSize p
163
164 pCopy :: Polygon -> Polygon
165 pCopy p = mkPolygon $ V.generate (pSize p) $ pAccess p
166
167 pSetOffset :: APolygon a -> Int -> APolygon a
168 pSetOffset p offset =
169 p { polygonOffset = offset `mod` pSize p }
170
171 pAdjustOffset :: APolygon a -> Int -> APolygon a
172 pAdjustOffset p offset =
173 p { polygonOffset = (polygonOffset p + offset) `mod` pSize p }
174
175 {-# INLINE pSize #-}
176 pSize :: APolygon a -> Int
177 pSize = length . polygonPoints
178
179 pNull :: APolygon a -> Bool
180 pNull = V.null . polygonPoints
181
182 pNext :: APolygon a -> Int -> Int
183 pNext p i = (i+1) `mod` pSize p
184
185 pPrev :: APolygon a -> Int -> Int
186 pPrev p i = (i-1) `mod` pSize p
187
188 -- When is a polygon valid/simple?
189 -- It is counter-clockwise.
190 -- No edges intersect.
191 -- O(n^2)
192 -- 'checkEdge' takes 90% of the time.
193 pIsSimple :: Polygon -> Bool
194 pIsSimple p | pSize p < 3 = False
195 pIsSimple p = pIsCCW p && noDups && checkEdge 0 2
196 where
197 noDups = checkForDups (sort (V.toList (polygonPoints p)))
198 checkForDups (x:y:xs)
199 = x /= y && checkForDups (y:xs)
200 checkForDups _ = True
201 len = pSize p
202 -- check i,i+1 against j,j+1
203 -- j > i+1
204 checkEdge i j
205 | j >= len = (i > len-3) || checkEdge (i+1) (i+3)
206 | otherwise =
207 case lineIntersect (pAccess p i, pAccess p $ i+1)
208 (pAccess p j, pAccess p $ j+1) of
209 Just u | u /= pAccess p i -> False
210 _nothing -> checkEdge i (j+1)
211
212 pScale :: Rational -> Polygon -> Polygon
213 pScale s = pUnsafeMap (ringMap (^* s))
214
215 pAtCentroid :: Polygon -> Polygon
216 pAtCentroid p = pTranslate (negate c) p
217 where c = pCentroid p ^/ 2
218
219 pAtCenter :: Polygon -> Polygon
220 pAtCenter p = pTranslate (negate $ pCenter p) p
221
222 pTranslate :: V2 Rational -> Polygon -> Polygon
223 pTranslate v = pUnsafeMap (ringMap (+v))
224
225 pCenter :: Polygon -> V2 Rational
226 pCenter p = V2 (x+w/2) (y+h/2)
227 where
228 (x,y,w,h) = pBoundingBox p
229
230 -- Returns (min-x, min-y, width, height)
231 pBoundingBox :: Polygon -> (Rational, Rational, Rational, Rational)
232 pBoundingBox = \p ->
233 let V2 x y = pAccess p 0 in
234 case V.foldl' worker (x, y, 0, 0) (polygonPoints p) of
235 (xMin, yMin, xMax, yMax) ->
236 (xMin, yMin, xMax-xMin, yMax-yMin)
237 where
238 worker (xMin,yMin,xMax,yMax) (V2 thisX thisY) =
239 (min xMin thisX, min yMin thisY
240 ,max xMax thisX, max yMax thisY)
241
242 -- Place n points on a circle, use one parameter to slide the points back and forth.
243 -- Use second parameter to move points closer to center circle.
244 pGenerate :: [(Double, Double)] -> Polygon
245 pGenerate points
246 | len < 4 = error "pGenerate: require at least four points"
247 | otherwise = mkPolygon $ V.fromList
248 [ V2 (realToFrac $ cos ang * rMod)
249 (realToFrac $ sin ang * rMod)
250 | (i,(angMod,rMod)) <- zip [0..] points
251 , let minAngle = tau / len * i - pi
252 maxAngle = tau / len * (i+1) - pi
253 ang = minAngle + (maxAngle-minAngle)*angMod
254 ]
255 where
256 tau = 2*pi
257 len = fromIntegral (length points)
258
259 pUnGenerate :: Polygon -> [(Double, Double)]
260 pUnGenerate p =
261 [ worker i (fmap realToFrac e)
262 | (i,e) <- zip [0..] (V.toList $ polygonPoints p) ]
263 where
264 len = fromIntegral (pSize p)
265 worker i (V2 x y) =
266 let ang = atan2 y x
267 minAngle = tau / len * i - pi
268 maxAngle = tau / len * (i+1) - pi
269 in ((ang-minAngle)/(maxAngle-minAngle), sqrt (x*x+y*y))
270 tau = 2*pi
271
272 -- When is a triangulation valid?
273 -- Intersection: No internal edges intersect.
274 -- Completeness: All edge neighbours share a single internal edge.
275 isValidTriangulation :: Polygon -> Triangulation -> Bool
276 isValidTriangulation p t = isComplete && intersectionFree
277 where
278 o = polygonOffset p
279 isComplete = all isProper [0 .. pSize p-1]
280 isProper i =
281 let j = pNext p i in
282 length ((pPrev p i : (t V.! i)) `intersect` (pNext p j : t V.! j)) == 1
283 intersectionFree = and
284 [ case lineIntersect (pAccess p (a-o), pAccess p (b-o)) (pAccess p (c-o), pAccess p (d-o)) of
285 Nothing -> True
286 Just u -> u == pAccess p (a-o) || u == pAccess p (b-o) ||
287 u == pAccess p (c-o) || u == pAccess p (d-o)
288 | ((a,b),(c,d)) <- edgePairs ]
289 edgePairs = [ (e1, e2) | (e1, rest) <- zip edges (drop 1 $ tails edges), e2 <- rest]
290 edges =
291 [ (n, i)
292 | (n, lst) <- zip [0..] (V.toList t)
293 , i <- lst
294 , n < i
295 ]
296
297 triangulationsToPolygons :: Polygon -> Triangulation -> [Polygon]
298 triangulationsToPolygons p t =
299 [ mkPolygon $ V.fromList
300 [ pAccess p g, pAccess p i, pAccess p j ]
301 | i <- [0 .. pSize p-1]
302 , let js = filter (i<) $ t V.! i
303 , (g, j) <- zip (i-1:js) js
304 ]
305
306 pIsInside :: Polygon -> V2 Rational -> Bool
307 pIsInside p point = or
308 [ isInside (rawAccess g) (rawAccess i) (rawAccess j) point
309 | i <- [0 .. pSize p-1]
310 , let js = filter (i<) $ polygonTriangulation p V.! i
311 , (g, j) <- zip (i-1:js) js
312 ]
313 where
314 rawAccess x = polygonPoints p V.! x
315
316 -- reducePolygons :: Int -> [Polygon] -> [Polygon]
317 -- reducePolygons n ps
318 -- | length ps <= n = ps
319 -- | otherwise =
320 -- let p = findSmallest ps
321 -- es = edges p
322 -- e = findSmallest es
323 -- in reducePolygons n (merge p e : delete p (delete e ps))
324 -- where
325 -- findSmallest = minimumBy (comparing area2X)
326 -- shareEdge p1 p2 =
327
328 {-# INLINE pAccess #-}
329 pAccess :: APolygon a -> Int -> V2 a
330 pAccess p i = -- polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p)
331 polygonPoints p `V.unsafeIndex` ((polygonOffset p + i) `mod` pSize p)
332
333 triangle :: Polygon
334 triangle = mkPolygon $ V.fromList [V2 1 1, V2 0 0, V2 2 0]
335
336 triangle' :: [P]
337 triangle' = reverse [V2 1 1, V2 0 0, V2 2 0]
338
339 shape1 :: Polygon
340 shape1 = mkPolygon $ V.fromList
341 [ V2 0 0, V2 2 0
342 , V2 2 1, V2 2 2, V2 2 3, V2 2 4, V2 2 5, V2 2 6
343 , V2 1 1, V2 0 1 ]
344
345 shape2 :: Polygon
346 shape2 = mkPolygon $ V.fromList
347 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 0 (-1), V2 0 (-2)
348 , V2 3 (-2), V2 3 2, V2 0 2]
349
350 shape3 :: Polygon
351 shape3 = mkPolygon $ V.fromList
352 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 2, V2 0 2]
353
354 shape4 :: Polygon
355 shape4 = mkPolygon $ V.fromList
356 [ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 3 (-1),V2 3 2, V2 0 2]
357
358 shape5 :: Polygon
359 shape5 = pCycles shape4 !! 2
360
361 -- square
362 shape6 :: Polygon
363 shape6 = mkPolygon $ V.fromList [ V2 0 0, V2 1 0, V2 1 1, V2 0 1 ]
364
365 shape7 :: Polygon
366 shape7 = pScale 6 $ mkPolygon $ V.fromList
367 [V2 ((-1567171105775771) % 144115188075855872) ((-7758063241391039) % 1152921504606846976)
368 ,V2 ((-2711114907999263) % 18014398509481984) ((-3561889280168807) % 18014398509481984)
369 ,V2 ((-6897139157863177) % 72057594037927936) ((-1632144794297397) % 4503599627370496)
370 ,V2 (5592137945106423 % 36028797018963968) ((-71351641856107) % 281474976710656)
371 ,V2 (2568147525079071 % 4503599627370496) ((-4312925637247687) % 18014398509481984)
372 ,V2 (1291079014395023 % 2251799813685248) (321513444515769 % 2251799813685248)
373 ,V2 (2071709221627247 % 4503599627370496) (4019115966736491 % 9007199254740992)
374 ,V2 ((-1589087869859839) % 144115188075855872) (4904023654354179 % 9007199254740992)
375 ,V2 ((-2328090886101149) % 36028797018963968) (2587887893460759 % 36028797018963968)
376 ,V2 ((-7990199074159871) % 18014398509481984) (1301850651537745 % 4503599627370496)]
377
378 shape8 :: Polygon
379 shape8 = pScale 10 $ pGenerate
380 [(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)]
381
382 shape9 :: Polygon
383 shape9 = pScale 5 $ pGenerate
384 [(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)]
385
386 shape10 :: Polygon
387 shape10 = pGenerate
388 [(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)]
389
390 shape11 :: Polygon
391 shape11 = pGenerate
392 [(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)]
393
394 shape12 :: Polygon
395 shape12 = mkPolygon $ V.fromList
396 [ V2 0 0, V2 0.5 1.5, V2 2 2, V2 (-2) 2, V2 (-0.5) 1.5 ]
397
398 -- F shape
399 shape13 :: Polygon
400 shape13 = pCycles (mkPolygon $ V.reverse (V.fromList
401 [ V2 0 0, V2 0 2
402 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1
403 , V2 1 1, V2 1 0.7
404 , V2 0.3 0.7, V2 0.3 0 ])) !! 7
405
406 -- E shape
407 shape14 :: Polygon
408 shape14 = pCycles (mkPolygon $ V.reverse $ V.fromList
409 [ V2 0 0, V2 0 2 -- up
410 , V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 -- first prong
411 , V2 1 1, V2 1 0.7, V2 0.3 0.7, V2 0.3 0.3 -- second prong
412 , V2 1 0.3, V2 1 0 -- last prong
413 ]) !! 9
414
415 --
416 shape15 :: Polygon
417 shape15 = mkPolygon $ V.fromList
418 [ V2 0 0, V2 2 0
419 , V2 2 2, V2 1 2
420 , V2 1 1, V2 0 1]
421
422 shape16 :: Polygon
423 shape16 = mkPolygon $ V.fromList
424 [ V2 0 0, V2 2 0
425 , V2 2 1, V2 1 1
426 , V2 1 2, V2 0 2]
427
428 shape17 :: Polygon
429 shape17 = mkPolygon $ V.fromList
430 [ V2 2 0, V2 2 1
431 , V2 1 1, V2 1 2
432 , V2 0 2, V2 0 1, V2 0 0 ]
433
434 shape18 :: Polygon
435 shape18 = mkPolygon $ V.fromList
436 [ V2 2 0, V2 2 1, V2 2 2
437 , V2 1 2, V2 1 1
438 , V2 0 1, V2 0 0 ]
439
440 shape19 :: Polygon
441 shape19 = mkPolygon $ V.fromList
442 [ V2 (-3) (-3), V2 0 (-1)
443 , V2 3 (-3), V2 1 0
444 , V2 3 3, V2 0 1
445 , V2 (-3) 3, V2 (-1) 0 ]
446
447 shape20 :: Polygon
448 shape20 = mkPolygon $ V.fromList
449 [ V2 (-3) (-3)
450 , V2 0 (-1)
451 , V2 3 (-3)
452 , V2 5 0
453 , V2 2.5 (-2)
454 , V2 1 0
455 , V2 3 3
456 , V2 0 1
457 , V2 (-3) 3
458 , V2 (-1) 0 ]
459
460 shape21 :: Polygon
461 shape21 = mkPolygon $ V.fromList
462 [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)
463 ,V2 3.0 2.0,V2 0.0 2.0]
464
465 shape22 :: Polygon
466 shape22 = pScale 2 $ mkPolygon $ V.fromList
467 [V2 (-0.17) (-0.08)
468 ,V2 (-0.34) (-0.21)
469 ,V2 0.0 0.0
470 ,V2 (-0.10) 0.60
471 ,V2 (-0.14) 0.19
472 ,V2 (-0.05) 0.03
473 ]
474
475 shape23 :: Polygon
476 shape23 = mkPolygon $ V.fromList
477 [ V2 0 0, V2 4 0
478 , V2 4 3, V2 2 3
479 , V2 2 2, V2 3 2
480 , V2 3 1, V2 1 1
481 , V2 1 2, V2 2 2
482 , V2 2 3, V2 0 3 ]
483
484 concave :: Polygon
485 concave = mkPolygon $
486 V.fromList [V2 0 0, V2 2 0, V2 2 2, V2 1 1, V2 0 2]
487
488 pMkWinding :: Int -> Polygon
489 pMkWinding n | n < 1 = error "Polygon must have at least one winding."
490 pMkWinding n = mkPolygon $
491 V.fromList $ p0 : p1 : walkTo p1 1 n (V2 1 0) ++ reverse (walkTo p0 1 (n+2) (V2 (-1) 0))
492 where
493 p0 = V2 0 0
494 p1 = V2 0 1
495 walkTo at a b dir
496 | a == b = []
497 | otherwise =
498 let newAt = at + (dir ^* toRational a)
499 in newAt : walkTo newAt (a+1) b (rot dir)
500 rot (V2 x y) =
501 V2 y (-x)
502
503 pDeoverlap :: Polygon -> Polygon
504 pDeoverlap p = mkPolygon arr
505 where
506 arr = V.generate (pSize p) worker
507 worker 0 = pAccess p 0
508 worker n =
509 if length (V.elemIndices (pAccess p n) (polygonPoints p)) /= 1
510 then
511 let prev = arr V.! (n-1)
512 this = pAccess p n
513 in lerp 0.99999 this prev
514 else pAccess p n
515
516 pCycles :: APolygon a -> [APolygon a]
517 pCycles p = map (pAdjustOffset p) [0 .. pSize p-1]
518
519 pCycle :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a
520 pCycle p 0 = p
521 pCycle p t = mkPolygon $ worker 0 0
522 where
523 worker acc i
524 | segment + acc > limit =
525 V.singleton (lerp (realToFrac $ (segment + acc - limit)/segment) x y) <>
526 -- V.drop (i+1) (polygonPoints p) <>
527 V.fromList (map (pAccess p) [i+1..pSize p-1]) <>
528 V.fromList (map (pAccess p) [0 .. i])
529 -- V.take (i+1) (polygonPoints p)
530 | i == pSize p-1 = V.fromList (map (pAccess p) [0 .. pSize p-1])
531 | otherwise = worker (acc+segment) (i+1)
532 where
533 x = pAccess p i
534 y = pAccess p $ i+1
535 segment = distance' x y
536 len = pCircumference' p
537 limit = t * len
538
539 pCentroid :: Fractional a => APolygon a -> V2 a
540 pCentroid p = V2 cx cy
541 where
542 a = pArea p
543 cx = recip (6*a) * V.sum (pMapEdges fnX p)
544 cy = recip (6*a) * V.sum (pMapEdges fnY p)
545 fnX (V2 x y) (V2 x' y') = (x+x')*(x*y' - x'*y)
546 fnY (V2 x y) (V2 x' y') = (y+y')*(x*y' - x'*y)
547
548 {-# INLINE pMapEdges #-}
549 pMapEdges :: (V2 a -> V2 a -> b) -> APolygon a -> V.Vector b
550 pMapEdges fn p = V.generate n $ \i ->
551 if i == n-1
552 then fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0)
553 else fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1))
554 where
555 n = pSize p
556 arr = polygonPoints p
557
558 {-# SPECIALIZE pArea :: APolygon Double -> Double #-}
559 {-# SPECIALIZE pArea :: APolygon Rational -> Rational #-}
560 pArea :: (Fractional a) => APolygon a -> a
561 pArea p =
562 -- 0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p)
563 0.5 * worker 0 0
564 where
565 fn (V2 x y) (V2 x' y') = x*y' - x'*y
566 arr = polygonPoints p
567 worker !acc i
568 | i == pSize p - 1 = acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0)
569 | otherwise =
570 worker (acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1))) (i+1)
571
572 pCircumference :: (Real a, Fractional a) => APolygon a -> a
573 pCircumference p = sum
574 [ approxDist (pAccess p i) (pAccess p $ i+1)
575 | i <- [0 .. pSize p-1]]
576
577 pCircumference' :: (Real a, Fractional a) => APolygon a -> Double
578 pCircumference' p = sum
579 [ distance' (pAccess p i) (pAccess p $ i+1)
580 | i <- [0 .. pSize p-1]]
581
582
583 -- Add points by splitting the longest lines in half repeatedly.
584 pAddPoints :: (Ord a, Real a, Fractional a) => Int -> APolygon a -> APolygon a
585 pAddPoints n p | n <= 0 = p
586 pAddPoints n p = pAddPoints (n-1) $
587 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1]
588 where
589 worker idx
590 | idx == longestEdge =
591 let start = pAccess p idx
592 end = pAccess p $ idx+1
593 middle = lerp 0.5 end start
594 in [start, middle]
595 | otherwise = [pAccess p idx]
596 longestEdge = maximumBy cmpLength [0 .. pSize p-1]
597 cmpLength a b =
598 distSquared (pAccess p a) (pAccess p $ a+1) `compare`
599 distSquared (pAccess p b) (pAccess p $ b+1)
600
601 pAddPointsRestricted :: (Fractional a, Ord a, Real a) => [(V2 a, V2 a)] -> Int -> APolygon a -> APolygon a
602 pAddPointsRestricted _immutableEdges n p | n <= 0 = p
603 pAddPointsRestricted immutableEdges n p = pAddPointsRestricted immutableEdges (n-1) $
604 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1]
605 where
606 isImmutable idx =
607 (pAccess p idx, pAccess p $ idx+1) `elem` immutableEdges ||
608 (pAccess p $ idx+1, pAccess p idx) `elem` immutableEdges
609 worker idx
610 | idx == longestEdge && not (isImmutable idx) =
611 let start = pAccess p idx
612 end = pAccess p $ idx+1
613 middle = lerp 0.5 end start
614 in [start, middle]
615 | otherwise = [pAccess p idx]
616 longestEdge = maximumBy cmpLength [0 .. pSize p-1]
617 cmpLength a _ | isImmutable a = LT
618 cmpLength _ b | isImmutable b = GT
619 cmpLength a b =
620 distSquared (pAccess p a) (pAccess p $ a+1) `compare`
621 distSquared (pAccess p b) (pAccess p $ b+1)
622
623 pAddPointsBetween :: (Fractional a, Ord a, Real a) => (Int, Int) -> Int -> APolygon a -> APolygon a
624 pAddPointsBetween _ n p | n <= 0 = p
625 pAddPointsBetween (i,l) n p = pAddPointsBetween (i,l+1) (n-1) $
626 mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1]
627 where
628 worker idx
629 | idx == longestEdge =
630 let start = pAccess p idx
631 end = pAccess p $ idx+1
632 middle = lerp 0.5 end start
633 in [start, middle]
634 | otherwise = [pAccess p idx]
635 longestEdge = maximumBy cmpLength [i .. i+l-1]
636 cmpLength a b =
637 distSquared (pAccess p a) (pAccess p $ a+1) `compare`
638 distSquared (pAccess p b) (pAccess p $ b+1)
639
640 -- addPoints :: Int -> Polygon -> Polygon
641 -- addPoints n p = mkPolygon $ V.fromList $ worker n 0 (map (pAccess p) [0..s])
642 -- where
643 -- worker 0 _ rest = init rest
644 -- worker i acc (x:y:xs) =
645 -- let xy = approxDist x y in
646 -- if acc + xy > limit
647 -- then x : worker (i-1) 0 (lerp ((limit-acc)/xy) y x : y:xs)
648 -- else x : worker i (acc+xy) (y:xs)
649 -- worker _ _ [_] = []
650 -- worker _ _ _ = error "addPoints: invalid polygon"
651 -- s = pSize p
652 -- len = polygonLength p
653 -- limit = len / fromIntegral (n+1)
654
655 pIsConvex :: Polygon -> Bool
656 pIsConvex p = and
657 [ area2X (pAccess p i) (pAccess p j) (pAccess p k) > 0
658 | i <- [0..n-1]
659 , j <- [i+1..n-1]
660 , k <- [j+1..n-1]
661 ]
662 where n = pSize p
663
664 pIsCCW :: Polygon -> Bool
665 pIsCCW p | pNull p = False
666 pIsCCW p = V.sum (pMapEdges fn p) < 0
667 where
668 fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1)
669
670 {-# INLINE pRayIntersect #-}
671 pRayIntersect :: (Fractional a, Ord a) => APolygon a -> (Int, Int) -> (Int,Int) -> Maybe (V2 a)
672 pRayIntersect p (a,b) (c,d) =
673 rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d)
674
675 pCuts :: (Real a, Fractional a) => APolygon a -> [(APolygon a,APolygon a)]
676 pCuts p =
677 [ pCutAt (pAdjustOffset p i) (j-i)
678 | i <- [0 .. pSize p-1 ]
679 , j <- [i+2 .. pSize p-1 ]
680 , (j+1) `mod` pSize p /= i
681 , pParent p i j == i ]
682
683 pCutEqual :: (Real a, Fractional a) => APolygon a -> (APolygon a, APolygon a)
684 pCutEqual p =
685 fromMaybe (p,p) $ listToMaybe $ sortOn f $ pCuts p
686 where
687 f (a,b) = abs (pArea a - pArea b)
688
689 -- FIXME: This should be more efficient
690 pCutAt :: (Real a, Fractional a) => APolygon a -> Int -> (APolygon a, APolygon a)
691 pCutAt p i = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right)
692 where
693 n = pSize p
694 left = map (pAccess p) [0 .. i]
695 right = map (pAccess p) (0:[i..n-1])
696
697 pOverlap :: (Fractional a, Ord a, Real a) => APolygon a -> APolygon a -> APolygon a
698 pOverlap a b = mkPolygon $ V.fromList $ clearDups $ concatMap edgeIntersect [0 .. pSize a-1]
699 where
700 clearDups (x:y:xs)
701 | x == y = clearDups (y:xs)
702 | otherwise = x : clearDups (y:xs)
703 clearDups xs = xs
704 edgeIntersect edge =
705 sortOn (distSquared (pAccess a edge)) $ catMaybes
706 [ lineIntersect (aP, aP') (bP, bP')
707 | i <- [0 .. pSize b-1]
708 , let aP = pAccess a edge
709 aP' = pAccess a (edge+1)
710 bP = pAccess b i
711 bP' = pAccess b (i+1)
712 ]
713
714 ---------------------------------------------------------
715 -- SSSP visibility and SSSP windows
716
717 ssspVisibility :: (Real a, Fractional a) => APolygon a -> APolygon a
718 ssspVisibility p = mkPolygon $
719 V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1])
720 where
721 clearDups (x:y:xs)
722 | x == y = clearDups (y:xs)
723 | otherwise = x : clearDups (y:xs)
724 clearDups xs = xs
725 obstructedBy n =
726 case pParent p 0 n of
727 0 -> n
728 i -> obstructedBy i
729 go [] = []
730 go [x] = [pAccess p x]
731 go (x:y:xs) =
732 let xO = obstructedBy x
733 yO = obstructedBy y
734 in case () of
735 ()
736 -- Both ends are visible.
737 | xO == x && yO == y -> pAccess p x : go (y:xs)
738 -- X is visible, x to intersect (0,yO) (x,y)
739 | xO == x ->
740 pAccess p x : fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y)) : go (y:xs)
741 -- Y is visible
742 | yO == y -> fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) : pAccess p y : go (y:xs)
743 -- Neither is visible and they've obstructed by the same point
744 -- so the entire edge is hidden.
745 | xO == yO -> go (y:xs)
746 -- Neither is visible. Cast shadow from obstruction points to
747 -- find if a subsection of the edge is visible.
748 | otherwise ->
749 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
750 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
751 in if a /= b
752 then a : b : go (y:xs)
753 else go (y:xs)
754
755 ssspWindows :: Polygon -> [(V2 Rational, V2 Rational)]
756 ssspWindows p = clearDups $ go (pAccess p 0) [0..pSize p-1]
757 where
758 clearDups (x:y:xs)
759 | x == y = clearDups (y:xs)
760 | otherwise = x : clearDups (y:xs)
761 clearDups xs = xs
762 obstructedBy n =
763 case pParent p 0 n of
764 0 -> n
765 i -> obstructedBy i
766 go _ [] = []
767 go _ [_] = []
768 go l (x:y:xs) =
769 let xO = obstructedBy x
770 yO = obstructedBy y
771 in case () of
772 ()
773 -- Both ends are visible.
774 | xO == x && yO == y -> go (pAccess p x) (y:xs)
775 -- X is visible, x to intersect (0,yO) (x,y)
776 | xO == x ->
777 go (fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y))) (y:xs)
778 -- Y is visible
779 | yO == y ->
780 let newL = fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) in
781 (l, newL) :
782 go newL (y:xs)
783 -- Neither is visible and they've obstructed by the same point
784 -- so the entire edge is hidden.
785 | xO == yO -> go l (y:xs)
786 -- Neither is visible. Cast shadow from obstruction points to
787 -- find if a subsection of the edge is visible.
788 | otherwise ->
789 let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
790 b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
791 in if a /= b
792 then (l, a) : (b, pAccess p yO) : go (pAccess p yO) (y:xs)
793 else go l (y:xs)
794
795 pdualPolygons :: Polygon -> PDual -> [Polygon]
796 pdualPolygons p pdual = map mkPolygonFromRing (pdualRings (pRing p) pdual)