{-# LANGUAGE ScopedTypeVariables #-}
module Reanimate.Math.Smooth
( steinerPoints
, renderMesh
, renderAMesh
, smoothMesh
, smoothStep
, angleSmooth
, meshMinAngle
, splitMeshEdges
)
where
import qualified Data.IntSet as ISet
import Control.Monad
import Control.Monad.ST
import Data.STRef
import Data.List
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import Linear.V2
import Linear.Metric
import Linear.Vector
import Reanimate.Animation
import Reanimate.Constants
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Math.Render
import Reanimate.Math.Triangulate
import Reanimate.Morph.Rigid
import Reanimate.Svg
steinerPoints :: Polygon -> [Polygon] -> [V2 Rational]
steinerPoints outerPolygon = concatMap
(V.toList . V.filter isSteiner . polygonPoints)
where isSteiner p = V.notElem p (polygonPoints outerPolygon)
renderMesh :: Polygon -> [Polygon] -> SVG
renderMesh outerPolygon innerPolygons = mkGroup
[ mkGroup
[ withFillOpacity 1
$ withStrokeWidth (defaultStrokeWidth * 0)
$ withStrokeColor "black"
$ withFillColor "lightgreen"
$ polygonShape p
| p <- innerPolygons
]
, mkGroup
[ withFillColor "red" $ aroundCenter (scale 0.2) $ drawPoint
(realToFrac <$> p)
| p <- steiners
]
]
where steiners = steinerPoints outerPolygon innerPolygons
renderAMesh :: Mesh -> SVG
renderAMesh m = mkGroup
[ translate (-1.5) 0 $ renderMesh
(mkPolygon outlineA)
[ mkPolygon $ V.fromList $ map (fmap realToFrac)
[pA V.! a, pA V.! b, pA V.! c]
| (a, b, c) <- V.toList (meshTriangles m)
]
, translate (1.5) 0 $ renderMesh
(mkPolygon outlineB)
[ mkPolygon $ V.fromList $ map (fmap realToFrac)
[pB V.! a, pB V.! b, pB V.! c]
| (a, b, c) <- V.toList (meshTriangles m)
]
]
where
pA = meshPointsA m
pB = meshPointsB m
outlineA = V.map (\i -> realToFrac <$> meshPointsA m V.! i) (meshOutline m)
outlineB = V.map (\i -> realToFrac <$> meshPointsB m V.! i) (meshOutline m)
drawPoint :: V2 Double -> SVG
drawPoint (V2 x y) = translate x y $ mkCircle 0.1
smoothMesh :: Mesh -> [Mesh]
smoothMesh m = runST $ do
trigs <- edgesToTriangulationM (length $ meshPointsA m) $ concat
[ [(a, b), (b, c), (c, a)] | (a, b, c) <- V.toList (meshTriangles m) ]
pA <- V.thaw $ V.map (fmap realToFrac) (meshPointsA m)
pB <- V.thaw $ V.map (fmap realToFrac) (meshPointsB m)
replicateM 120 $ do
smoothStep trigs pA steiner
smoothStep trigs pB steiner
newA <- V.freeze pA
newB <- V.freeze pB
return $ m { meshPointsA = newA, meshPointsB = newB }
where
steiner = V.fromList
[ i | i <- [0 .. length (meshPointsA m) - 1], V.notElem i (meshOutline m) ]
smoothStep
:: forall s
. V.MVector s [Int]
-> V.MVector s (V2 Double)
-> V.Vector Int
-> ST s ()
smoothStep triangulation pts steiner =
forM_ [0 .. length steiner - 1]
$ \s_i -> do
let i = steiner V.! s_i
pt <- MV.read pts i
edges <- V.fromList <$> (sortEdges pt =<< MV.read triangulation i)
let angleBased = angleSmooth pt edges
laplacian = sum edges ^/ (fromIntegral $ length edges)
if isValidLocation pt edges angleBased
then MV.write pts i angleBased
else if isValidLocation pt edges laplacian
then MV.write pts i laplacian
else return ()
where
sortEdges :: V2 Double -> [Int] -> ST s [V2 Double]
sortEdges pt edges = do
edgePoints <- mapM (MV.read pts) edges
return $ sortOn (dir pt) edgePoints
dir :: V2 Double -> V2 Double -> Double
dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a)))
angleSmooth :: V2 Double -> V.Vector (V2 Double) -> V2 Double
angleSmooth origin js = V.sum (V.generate n nth) ^/ V.sum (V.generate n factor)
where
n = length js
factor i =
let n_self = js V.! i
n_origin = origin - n_self
n_prev = js V.! mod (i - 1) n - n_self
n_next = js V.! mod (i + 1) n - n_self
a1 = acos (dot n_origin n_next / (norm n_origin * norm n_next))
a2 = acos (dot n_origin n_prev / (norm n_origin * norm n_prev))
alpha = a1 + a2
in recip (alpha * alpha)
nth i =
let V2 x y = origin
n_self@(V2 x_0 y_0) = js V.! i
n_origin = origin - n_self
n_prev = js V.! mod (i - 1) n - n_self
n_next = js V.! mod (i + 1) n - n_self
a1 = acos (dot n_origin n_next / (norm n_origin * norm n_next))
a2 = acos (dot n_origin n_prev / (norm n_origin * norm n_prev))
alpha = a1 + a2
b = (a2 - a1) / 2
x' = x_0 + (x - x_0) * cos b - (y - y_0) * sin b
y' = y_0 + (x - x_0) * sin b + (y - y_0) * cos b
in V2 x' y' ^/ (alpha * alpha)
isValidLocation :: V2 Double -> V.Vector (V2 Double) -> V2 Double -> Bool
isValidLocation origin edges newLoc =
or
[ isInside origin a b newLoc
| i <- [0 .. length edges - 1]
, let a = edges V.! i
b = edges V.! mod (i + 1) (length edges)
]
&& V.toList edges
== sortOn (dir newLoc) (V.toList edges)
&& minAngle origin edges
< minAngle newLoc edges
where
dir :: V2 Double -> V2 Double -> Double
dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a)))
minAngle :: V2 Double -> V.Vector (V2 Double) -> Double
minAngle origin edges = minimum $ concat
[ [a1, a2, a3]
| i <- [0 .. length edges - 1]
, let a = edges V.! i
b = edges V.! mod (i + 1) (length edges)
(a1, a2, a3) = triangleAngles origin a b
]
meshMinAngle :: Mesh -> (Double, Double)
meshMinAngle m =
( minimum $ concat
[ [a1, a2, a3]
| (a, b, c) <- V.toList (meshTriangles m)
, let (a1, a2, a3) = triangleAngles (meshPointsA m V.! a)
(meshPointsA m V.! b)
(meshPointsA m V.! c)
]
, minimum $ concat
[ [a1, a2, a3]
| (a, b, c) <- V.toList (meshTriangles m)
, let (a1, a2, a3) = triangleAngles (meshPointsB m V.! a)
(meshPointsB m V.! b)
(meshPointsB m V.! c)
]
)
splitMeshEdges :: Double -> Mesh -> Mesh
splitMeshEdges maxLen mesh = runST $ do
t <- MV.replicate (n * 2) []
forM_ (V.toList (meshTriangles mesh)) $ \(a, b, c) -> do
insertEdge t a b
insertEdge t b c
insertEdge t a c
forM_ [0 .. n - 1] $ \i -> MV.modify t (ISet.toList . ISet.fromList) i
triangulation <- V.freeze t
newPoints <- newSTRef []
offset <- newSTRef n
forM_ [0 .. n - 1] $ \i -> forM_ (triangulation V.! i) $ \j ->
when (i < j) $ when (maxEdgeDistanceSquared i j > maxLen * maxLen) $ do
p <- readSTRef offset
writeSTRef offset (p + 1)
let (p1, p2) = splitEdge i j
modifySTRef newPoints ((p1, p2) :)
let [c, d] = (triangulation V.! i) `intersect` (triangulation V.! j)
deleteEdge t i j
insertEdge t i p
insertEdge t j p
insertEdge t c p
insertEdge t d p
return undefined
where
splitEdge i j =
( lerp 0.5 (meshPointsA mesh V.! i) (meshPointsA mesh V.! j)
, lerp 0.5 (meshPointsB mesh V.! i) (meshPointsB mesh V.! j)
)
maxEdgeDistanceSquared i j =
distSquared (meshPointsA mesh V.! i) (meshPointsA mesh V.! j)
`max` distSquared (meshPointsB mesh V.! i) (meshPointsB mesh V.! j)
n = length (meshPointsA mesh)
deleteEdge v a b = do
MV.modify v (delete b) a
MV.modify v (delete a) b
insertEdge v a b = do
MV.modify v (b :) a
MV.modify v (a :) b