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