mirror of
https://github.com/reanimate/reanimate.git
synced 2026-09-13 09:02:54 +00:00
0d4a5801de* Move hmatrix dependency out of reanimate's core. * Move chiphunk to a separate package.
385 lines
14 KiB
Haskell
385 lines
14 KiB
Haskell
{-# LANGUAGE MultiWayIf #-}
|
|
{-# LANGUAGE ScopedTypeVariables #-}
|
|
module Reanimate.Math.Compatible where
|
|
|
|
import Data.Hashable
|
|
import Data.List
|
|
import Data.Maybe
|
|
import Data.Ord
|
|
import qualified Data.Vector as V
|
|
import Linear.V2
|
|
import Linear.Vector
|
|
import Reanimate.Animation
|
|
import Reanimate.Debug
|
|
import Reanimate.Math.Common
|
|
import Reanimate.Math.Polygon
|
|
import Reanimate.Svg
|
|
|
|
-- import Debug.Trace
|
|
-- import GHC.Stack
|
|
|
|
truncateP :: forall a. (Fractional a, Real a) => V2 a -> V2 a
|
|
truncateP = fmap (fromDouble . toDouble)
|
|
where
|
|
toDouble :: a -> Double
|
|
toDouble = fromRational . toRational
|
|
fromDouble :: Double -> a
|
|
fromDouble = fromRational . toRational
|
|
|
|
mkSteinerPoints :: (Fractional a) => V2 a -> V2 a -> Int -> [V2 a]
|
|
mkSteinerPoints a b s = [ lerp (fromIntegral i / (fromIntegral s + 1)) b a | i <- [1 .. s] ]
|
|
|
|
-- 0..i,j..n-1
|
|
-- i..j
|
|
split1Link :: (Real a, Fractional a) => APolygon a -> Int -> Int -> Int -> (APolygon a, APolygon a)
|
|
split1Link p i j s | j < i = split1Link p j i s
|
|
split1Link p i j s =
|
|
(mkPolygon $ V.fromList left, mkPolygon $ V.fromList right)
|
|
where
|
|
n = pSize p
|
|
sp = mkSteinerPoints (pAccess p i) (pAccess p j) s
|
|
left = map (pAccess p) [0 .. i] ++ sp ++ map (pAccess p) [j .. n - 1]
|
|
right = map (pAccess p) [i .. j] ++ reverse sp
|
|
|
|
steiner2Link :: (Real a, Fractional a) => APolygon a -> Int -> Int -> V2 a
|
|
steiner2Link p i j | j < i = steiner2Link p j i
|
|
steiner2Link p i j
|
|
| isNeighbour
|
|
= error "steiner2Link: Points are neighbours"
|
|
| isParent
|
|
= error "steiner2Link: Points can directly see each other."
|
|
| not (isStraightLine || isGrandparent || oneBendBetween p i j)
|
|
= error
|
|
$ "steiner2Link: Cannot construct 2-link chain between points: "
|
|
++ show (i, j, pParent p i j, pParent p i (pParent p i j))
|
|
| oneBendBetween p i j
|
|
= let p1 = lerp 0.5 (pAccess p i) (pAccess p j)
|
|
p2 = case p1 - pAccess p i of
|
|
V2 x y -> p1 + V2 (-y) x -- rotate 90 degrees.
|
|
intersects =
|
|
sortOn (approxDist p1)
|
|
$ [ u
|
|
| n <- [0 .. pSize p - 1]
|
|
, n < i || n >= j
|
|
, let edge = (pAccess p n, pAccess p $ pNext p n)
|
|
, u <- case rayIntersect (p1,p2) edge of
|
|
Nothing -> []
|
|
Just u -> [u]
|
|
, isBetween u edge
|
|
]
|
|
in lerp 0.5 p1 (intersects !! 0)
|
|
| otherwise
|
|
= let pI = ssspVisibility $ pSetOffset p i
|
|
pJ = ssspVisibility $ pSetOffset p j
|
|
in truncateP $ pCentroid $ pOverlap pI pJ
|
|
-- | otherwise
|
|
-- = truncateP $ lerp 0.5 windowCommon (intersects !! 0)
|
|
where
|
|
-- trigger = True -- i == 4 && j == 8
|
|
-- info = mkGroup
|
|
-- [ scale 2 $ mkGroup
|
|
-- [ mkGroup [withFillColor "grey" $ polygonShape p, polygonNumDots p]
|
|
-- , let V2 x1 y1 = realToFrac <$> windowCommon
|
|
-- in translate x1 y1 $ withFillColor "red" $ mkCircle 0.08
|
|
-- , let V2 x1 y1 = realToFrac <$> iWindow
|
|
-- in translate x1 y1 $ withFillColor "blue" $ mkCircle 0.04
|
|
-- , let V2 x1 y1 = realToFrac <$> jWindow
|
|
-- in translate x1 y1 $ withFillColor "green" $ mkCircle 0.04
|
|
-- ]
|
|
-- , translate 5 0 $ withFillColor "white" $ latex $ T.pack $ show (i, j) ]
|
|
-- distToV = approxDist windowCommon
|
|
isNeighbour = i == pNext p j || i == pPrev p j
|
|
isParent = pParent p i j == i
|
|
isGrandparent = pParent p i (pParent p i j) == i
|
|
isStraightLine =
|
|
direction (pAccess p j) (pAccess p $ pParent p i j) (pAccess p i) == 0
|
|
-- intersects =
|
|
-- sortOn distToV
|
|
-- $ snd windowDirection
|
|
-- : [ u
|
|
-- | n <- [0 .. pSize p - 1]
|
|
-- , let edge = (pAccess p n, pAccess p $ pNext p n)
|
|
-- , u <- case rayIntersect windowDirection edge of
|
|
-- Nothing -> []
|
|
-- Just u -> [u]
|
|
-- , isBetween u edge
|
|
-- , u /= windowCommon
|
|
-- , isForward windowDirection u
|
|
-- ] ++
|
|
-- [ u
|
|
-- | not isStraightLine
|
|
-- , ray <- [(pAccess p i, jWindow), (pAccess p j, iWindow)]
|
|
-- , u <- maybeToList $ rayIntersect windowDirection ray ]
|
|
-- iP = pAdjustOffset p i
|
|
-- jP = pAdjustOffset p j
|
|
-- windowDirection = (windowCommon, windowOpposite)
|
|
-- (windowCommon, iWindow, jWindow, windowOpposite)
|
|
-- | isStraightLine
|
|
-- = let p1 = lerp 0.5 (pAccess p i) (pAccess p j)
|
|
-- p2 = case p1 - pAccess p i of
|
|
-- V2 x y -> p1 + V2 (-y) x -- rotate 90 degrees.
|
|
-- in (p1, undefined, undefined, p2)
|
|
-- | otherwise
|
|
-- = fromMaybe
|
|
-- (error $ "No window overlap: " ++ show
|
|
-- (isStraightLine, isGrandparent, oneBendBetween p i j)
|
|
-- )
|
|
-- $ listToMaybe
|
|
-- $ [ (p1, p2, p3, p2+p3-p1)
|
|
-- | (a, b) <- ssspWindows iP
|
|
-- , (c, d) <- ssspWindows jP
|
|
-- , (p1, p2, p3) <- if
|
|
-- | a == c -> pure (a, b, d)
|
|
-- | a == d -> pure (a, b, c)
|
|
-- | b == c -> pure (b, a, d)
|
|
-- | b == d -> pure (b, a, c)
|
|
-- | otherwise -> []
|
|
-- ]
|
|
-- isForward (a, b) v = not (isBetween a (b, v))
|
|
|
|
-- 0..i,s,j..n-1
|
|
-- i..j,s
|
|
split2Link :: (Real a, Fractional a) => APolygon a -> Int -> Int -> (APolygon a, APolygon a)
|
|
split2Link p i j | j < i = split2Link p j i
|
|
split2Link p i j = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right)
|
|
where
|
|
s = steiner2Link p i j
|
|
n = pSize p
|
|
left = map (pAccess p) [0 .. i] ++ [s] ++ map (pAccess p) [j .. n - 1]
|
|
right = map (pAccess p) [i .. j] ++ [s]
|
|
|
|
data Link = OneLink | TwoLink
|
|
|
|
splitNLink :: (Real a, Fractional a) => APolygon a -> Int -> [(Link, Int)] -> (APolygon a, APolygon a)
|
|
splitNLink p i js = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right)
|
|
where
|
|
n = pSize p
|
|
left = map (pAccess p) [0 .. i] ++ steiners ++ map (pAccess p) [j .. n - 1]
|
|
right = map (pAccess p) [i .. j] ++ reverse steiners
|
|
j = snd (last js)
|
|
steiners = splitNLinks p i js
|
|
-- 0, [(TwoLink,2),(OneLink,3)]
|
|
-- 0, [(OneLink,5),(TwoLink,3)]
|
|
splitNLinks :: (Real a, Fractional a) => APolygon a -> Int -> [(Link, Int)] -> [V2 a]
|
|
splitNLinks _p _i [] = []
|
|
splitNLinks p i [(TwoLink, j )] = [steiner2Link p i j]
|
|
splitNLinks _p _i [(OneLink, _j)] = []
|
|
splitNLinks p i ((TwoLink, j) : (OneLink, j') : xs) =
|
|
let
|
|
(l, r) = split2Link p i j
|
|
p' = selectContains l r (pAccess p j')
|
|
s = steiner2Link p i j
|
|
sIdx =
|
|
fromMaybe (error "missing steiner") $ V.elemIndex s (polygonPoints p')
|
|
in
|
|
s : splitNLinks p' sIdx ((TwoLink, j') : xs)
|
|
splitNLinks p i ((TwoLink, j) : (TwoLink, j') : xs) =
|
|
let
|
|
(l, r) = split2Link p i j
|
|
p' = selectContains l r (pAccess p j')
|
|
s = steiner2Link p i j
|
|
sIdx =
|
|
fromMaybe (error "missing steiner") $ V.elemIndex s (polygonPoints p')
|
|
in
|
|
s : splitNLinks p' sIdx ((OneLink, j) : (TwoLink, j') : xs)
|
|
splitNLinks p i ((OneLink, j) : (TwoLink, j') : xs) =
|
|
let
|
|
(l, r) = split2Link p j j'
|
|
p' = selectContains l r (pAccess p i)
|
|
p'' = selectContains l r (pAccess p j')
|
|
s = steiner2Link p j j'
|
|
s' = steiner2Link p' i sIdx
|
|
sIdx = fromMaybe (error "missing steiner sIdx")
|
|
$ V.elemIndex s (polygonPoints p')
|
|
sIdx' = fromMaybe (error "missing steiner sIdx'")
|
|
$ V.elemIndex s' (polygonPoints p'')
|
|
in
|
|
s' : s : splitNLinks p'' sIdx' ((OneLink, j') : xs)
|
|
splitNLinks _p _i _ = error "splitNLinks: invalid input"
|
|
|
|
selectContains :: Eq a => APolygon a -> APolygon a -> V2 a -> APolygon a
|
|
selectContains p1 p2 elt | V.elem elt (polygonPoints p1) = p1
|
|
| V.elem elt (polygonPoints p2) = p2
|
|
| otherwise = error "elt not member of either polygons"
|
|
-- [(P, Vis)] -> Polygon -> [P]
|
|
{-
|
|
If there is a 1-link path, connect it.
|
|
Recurse into the two new polygon pairs.
|
|
|
|
If there is a 2-link path, connect it, possibly add bogus steiner point.
|
|
Recurse into the two new polygon pairs.
|
|
|
|
If there is a n-link path, connect it, add bogus steiner points.
|
|
Recurse into the two new polygon pairs.
|
|
-}
|
|
|
|
-- a 2 b 1 c
|
|
-- ab = 2-link of (a,b)
|
|
-- abc = 2-link of (cut-a-ab-b, c)
|
|
-- steps: a ab abc c
|
|
-- a, [(2,b), (1,c)]
|
|
-- [(ab, p1,p2)]
|
|
|
|
{-
|
|
, mkGroup
|
|
[ withStrokeColor "green" $mkGroup
|
|
[ mkLine (x1,y1) (x2,y2)
|
|
, mkLine (x1',y1') (x2',y2') ]
|
|
| (a,b) <- mWins
|
|
, (i,j) <- oWins
|
|
, a == i || a == j || b == i || b == j
|
|
, not (sort [a,b] == sort [i,j])
|
|
, let V2 x1 y1 = realToFrac <$> a
|
|
V2 x2 y2 = realToFrac <$> b
|
|
V2 x1' y1' = realToFrac <$> i
|
|
V2 x2' y2' = realToFrac <$> j
|
|
]
|
|
-}
|
|
|
|
|
|
|
|
type Points = V.Vector (V2 Rational)
|
|
type Edges = [(Int, Int, Int)]
|
|
data Mesh = Mesh { meshPoints :: Points, meshEdges :: Edges }
|
|
data MeshPair = MeshPair Points Points Edges
|
|
-- The points in a RelMesh are:
|
|
-- relMeshStatic ++ x where Ax = B
|
|
-- data RelMesh = RelMesh
|
|
-- { relMeshStatic :: Points
|
|
-- , relMeshEdges :: Edges
|
|
-- , relMeshA :: Matrix Double
|
|
-- , relMeshB :: Matrix Double
|
|
-- }
|
|
-- data RelMeshPair = RelMeshPair Points Edges (Matrix Double) (Matrix Double) (Matrix Double) (Matrix Double)
|
|
--
|
|
--
|
|
|
|
-- Use SSSP to find 1-link connections
|
|
-- Compute visibility windows for each node
|
|
-- Check windows for ovelap, giving 2-link connections.
|
|
-- Use floyd to generate all other connections
|
|
-- (a,b,1)
|
|
-- (a,b,2)
|
|
-- floyd
|
|
compatiblyTriangulateP :: (Real a, Fractional a) => APolygon a -> APolygon a -> [(APolygon a, APolygon a)]
|
|
compatiblyTriangulateP a b
|
|
| pSize a /= pSize b = error "polygon size mismatch"
|
|
| otherwise = traceSVG (showStep a b) $ compatiblyTriangulateP'
|
|
(pSetOffset a 0)
|
|
(pSetOffset a 0)
|
|
(pSetOffset b 0)
|
|
|
|
showStep :: APolygon a -> APolygon a -> SVG
|
|
showStep _a _b = mkGroup []
|
|
-- [ translate (-3) 0 $ scale 2
|
|
-- $ mkGroup [withFillColor "grey" $ polygonShape a, polygonNumDots a]
|
|
-- , translate 3 0 $ scale 2
|
|
-- $ mkGroup [withFillColor "grey" $ polygonShape b, polygonNumDots b]
|
|
-- ]
|
|
|
|
showSolution :: APolygon a -> APolygon a -> (Int, Int) -> SVG
|
|
showSolution _a _b (_nodeL, _nodeR) = mkGroup []
|
|
-- [ showStep a b
|
|
-- -- , translate (-3) 0 $ scale 2 $
|
|
-- -- translate xA yA $ withFillColor "red" $
|
|
-- -- mkCircle 0.05
|
|
-- , withFillColor "white" $ latex $ T.pack $ show (nodeL, nodeR)
|
|
-- ]
|
|
where
|
|
-- V2 xA yA = realToFrac <$> steiner2Link a nodeL nodeR
|
|
-- V2 xB yB = steiner2Link b nodeL nodeR
|
|
|
|
compatiblyTriangulateP' :: (Real a, Fractional a) => APolygon a -> APolygon a -> APolygon a -> [(APolygon a, APolygon a)]
|
|
compatiblyTriangulateP' aOrigin a b
|
|
-- | trace (show (pSize a, pSize b)) False = undefined
|
|
| n == 3 = traceSVG (showStep a b) $ {- trace ("Done") $ -} [(a, b)]
|
|
| otherwise =
|
|
-- trace ("PolygonA " ++ show a) $
|
|
-- trace ("PolygonB " ++ show b) $
|
|
-- trace ("pSize: " ++ show (pSize a)) $
|
|
-- trace ("aOneLink: " ++ show aOneLink) $
|
|
-- trace ("bOneLink: " ++ show bOneLink) $
|
|
-- trace ("aTwoLink: " ++ show aTwoLink) $
|
|
-- trace ("bTwoLink: " ++ show bTwoLink) $
|
|
-- trace ("bestOneLink: " ++ show bestOneLink) $
|
|
-- trace ("bestTwoLink: " ++ show bestTwoLink) $
|
|
traceSVG (showStep a b) $ case bestOneLink of
|
|
Nothing -> case bestTwoLink of
|
|
Nothing -> error $ "no 2-links"
|
|
Just (nodeL, nodeR) -> --trace ("Creating steiner node: " ++ show (nodeL, nodeR)) $
|
|
let isTwolink =
|
|
(nodeL, nodeR) `notElem` aOneLink ||
|
|
(nodeL, nodeR) `notElem` bOneLink
|
|
(aL, aR) = if (nodeL, nodeR) `elem` aOneLink
|
|
then split1Link a nodeL nodeR 1
|
|
else split2Link a nodeL nodeR
|
|
(bL, bR) = if (nodeL, nodeR) `elem` bOneLink
|
|
then split1Link b nodeL nodeR 1
|
|
else split2Link b nodeL nodeR
|
|
in (if isTwolink then traceSVG (showSolution a b (nodeL, nodeR)) else id) $
|
|
(compatiblyTriangulateP' aOrigin aL bL
|
|
++ compatiblyTriangulateP' aOrigin aR bR
|
|
)
|
|
Just (nodeL, nodeR) ->
|
|
{-trace
|
|
(show ("one link" :: String, toOriginIndex nodeL, toOriginIndex nodeR)
|
|
)
|
|
$ -}
|
|
let (aL, aR) = split1Link a nodeL nodeR 0
|
|
(bL, bR) = split1Link b nodeL nodeR 0
|
|
in compatiblyTriangulateP' aOrigin aL bL
|
|
++ compatiblyTriangulateP' aOrigin aR bR
|
|
where
|
|
-- toOriginIndex idx =
|
|
-- ( idx
|
|
-- , fromMaybe (-1) (V.elemIndex (pAccess a idx) (polygonPoints aOrigin))
|
|
-- + polygonOffset aOrigin
|
|
-- )
|
|
n = pSize a
|
|
bestOneLink =
|
|
listToMaybe (sortOn (Down . nodeDist) (aOneLink `intersect` bOneLink))
|
|
bestTwoLink = listToMaybe
|
|
(sortOn (Down . nodeDist)
|
|
((aOneLink ++ aTwoLink) `intersect` (bOneLink ++ bTwoLink))
|
|
)
|
|
aOneLink = polygonOneLinks a
|
|
bOneLink = polygonOneLinks b
|
|
aTwoLink = polygonTwoLinks a
|
|
bTwoLink = polygonTwoLinks b
|
|
nodeDist (i, j) = (min (j - i) (n - j + i), hash (i,j))
|
|
|
|
oneBendBetween :: (Real a, Fractional a) => APolygon a -> Int -> Int -> Bool
|
|
-- oneBendBetween _p _a _b = False
|
|
oneBendBetween p a b =
|
|
abs (direction
|
|
(pAccess p a)
|
|
(pAccess p (pParent p a b))
|
|
(pAccess p b)) < epsilon
|
|
where
|
|
-- obstructedBy n =
|
|
-- case pParent p a n of
|
|
-- i -> if i == a then n else obstructedBy i
|
|
|
|
polygonTwoLinks :: (Real a, Fractional a) => APolygon a -> [(Int, Int)]
|
|
polygonTwoLinks p =
|
|
[ (i, j)
|
|
| i <- [0 .. n - 1]
|
|
, j <- [i + 2 .. n - 1]
|
|
, not (i == 0 && j == n - 1)
|
|
, pParent p i j /= i -- check for 1-link
|
|
, let isTwoLink = pParent p i (pParent p i j) == i
|
|
isStraightLine = oneBendBetween p i j
|
|
-- Points on a straight line should be 2-link even though they are not
|
|
-- direct grandparents.
|
|
, isTwoLink || isStraightLine
|
|
]
|
|
where n = pSize p
|
|
|
|
polygonOneLinks :: APolygon a -> [(Int, Int)]
|
|
polygonOneLinks p =
|
|
[ (i, j)
|
|
| i <- [0 .. pSize p - 1]
|
|
, j <- [i + 2 .. pSize p - 1]
|
|
, not (i == 0 && j == pSize p - 1)
|
|
, pParent p i j == i
|
|
]
|