From 0d4a5801debe26f1094a11a649e8d6a894646b41 Mon Sep 17 00:00:00 2001 From: David Himmelstrup Date: Fri, 7 Aug 2020 15:34:26 +0800 Subject: [PATCH] Morphology (#112) * Move hmatrix dependency out of reanimate's core. * Move chiphunk to a separate package. --- .azure/azure-linux-template.yml | 4 - .azure/azure-windows-template.yml | 15 +- examples/shatter.hs | 130 -------- examples/simulate_equation.hs | 189 ----------- examples/simulate_gravity.golden | 10 - examples/simulate_gravity.hs | 79 ----- examples/tut_glue_physics.hs | 59 ---- mkdocs.yml | 4 +- .../morphology_leastwork.hs | 0 .../morphology_linebend.hs | 0 .../morphology_linebend_intro.hs | 0 .../Math => reanimate-morph}/Compatible.hs | 238 +++++++++----- reanimate-morph/LeastDifference.hs | 195 +++++++++++ .../Morph => reanimate-morph}/LeastWork.hs | 0 .../Morph => reanimate-morph}/LineBend.hs | 0 .../Math => reanimate-morph}/Render.hs | 7 +- .../morphology_point_trajectory.hs | 0 reanimate.cabal | 77 +++-- src/Reanimate/Chiphunk.hs | 97 ------ src/Reanimate/Math/Common.hs | 39 ++- src/Reanimate/Math/DCEL.hs | 83 +++-- src/Reanimate/Math/Polygon.hs | 118 +++++-- src/Reanimate/Morph/Common.hs | 5 +- src/Reanimate/Morph/Linear.hs | 9 +- src/Reanimate/Morph/Rigid.hs | 105 ++++-- stack-lts-12.yaml | 13 +- videos/morph/dcel.hs | 42 ++- videos/morph/do_morph.hs | 286 ++++++++++++++++ videos/morph/least-difference.hs | 127 +++++++ videos/morph/morph.hs | 285 +--------------- videos/morph/rigid.hs | 309 ++++++++++++++++++ 31 files changed, 1408 insertions(+), 1117 deletions(-) delete mode 100755 examples/shatter.hs delete mode 100755 examples/simulate_equation.hs delete mode 100644 examples/simulate_gravity.golden delete mode 100755 examples/simulate_gravity.hs delete mode 100755 examples/tut_glue_physics.hs rename {examples => reanimate-examples}/morphology_leastwork.hs (100%) rename {examples => reanimate-examples}/morphology_linebend.hs (100%) rename {examples => reanimate-examples}/morphology_linebend_intro.hs (100%) rename {src/Reanimate/Math => reanimate-morph}/Compatible.hs (54%) create mode 100644 reanimate-morph/LeastDifference.hs rename {src/Reanimate/Morph => reanimate-morph}/LeastWork.hs (100%) rename {src/Reanimate/Morph => reanimate-morph}/LineBend.hs (100%) rename {src/Reanimate/Math => reanimate-morph}/Render.hs (97%) rename {examples => reanimate-morph}/morphology_point_trajectory.hs (100%) delete mode 100644 src/Reanimate/Chiphunk.hs create mode 100644 videos/morph/do_morph.hs create mode 100644 videos/morph/least-difference.hs create mode 100644 videos/morph/rigid.hs diff --git a/.azure/azure-linux-template.yml b/.azure/azure-linux-template.yml index 072f8cd..71e0c38 100644 --- a/.azure/azure-linux-template.yml +++ b/.azure/azure-linux-template.yml @@ -11,10 +11,6 @@ jobs: BUILD: stack STACK_YAML: stack.yaml ARGS: --pedantic - stack-no-hmatrix: - BUILD: stack - STACK_YAML: stack.yaml - ARGS: --flag reanimate:disable-hmatrix --flag reanimate:-test stack-lts-15: BUILD: stack STACK_YAML: stack-lts-15.yaml diff --git a/.azure/azure-windows-template.yml b/.azure/azure-windows-template.yml index 552ee15..7ab2cb0 100644 --- a/.azure/azure-windows-template.yml +++ b/.azure/azure-windows-template.yml @@ -34,23 +34,26 @@ jobs: displayName: Cache stack root - bash: | set -o xtrace - $CONDA/condabin/conda.bat install --yes --quiet -c msys2 m2w64-gcc-libgfortran m2w64-openblas + # $CONDA/condabin/conda.bat install --yes --quiet -c msys2 m2w64-gcc-libgfortran m2w64-openblas + + # cp $CONDA/Library/mingw-w64/bin/libgfortran-3.dll $CONDA/Library/mingw-w64/bin/libgfortran.dll - cp $CONDA/Library/mingw-w64/bin/libgfortran-3.dll $CONDA/Library/mingw-w64/bin/libgfortran.dll - curl -sSkL http://www.stackage.org/stack/windows-x86_64 -o /usr/bin/stack.zip unzip -o /usr/bin/stack.zip -d /usr/bin/ displayName: 'Setup' - bash: | set -o xtrace - stack --install-ghc build $ARGS --no-keep-going --fast --only-dependencies --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin - stack --install-ghc test $ARGS --no-keep-going --fast --flag reanimate:test --only-dependencies --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin + #stack --install-ghc build $ARGS --no-keep-going --fast --only-dependencies --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin + stack --install-ghc build $ARGS --no-keep-going --fast --only-dependencies + #stack --install-ghc test $ARGS --no-keep-going --fast --flag reanimate:test --only-dependencies --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin + stack --install-ghc test $ARGS --no-keep-going --fast --flag reanimate:test --only-dependencies displayName: 'Install dependencies' - bash: | set -o xtrace # stack ./examples/counter.hs check - stack test --flag reanimate:test --test-arguments="-j2" --fast --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin + # stack test --flag reanimate:test --test-arguments="-j2" --fast --flag reanimate:enable-hmatrix --flag hmatrix:openblas --extra-lib-dirs=$CONDA/Library/mingw-w64/bin + stack test --flag reanimate:test --test-arguments="-j2" --fast displayName: 'Build & Test' - bash: | # stack haddock --no-haddock-deps diff --git a/examples/shatter.hs b/examples/shatter.hs deleted file mode 100755 index d78ccce..0000000 --- a/examples/shatter.hs +++ /dev/null @@ -1,130 +0,0 @@ -#!/usr/bin/env stack --- stack runghc --package reanimate -{-# LANGUAGE OverloadedStrings #-} -module Main (main) where - -import Chiphunk.Low -import Control.Monad -import Graphics.SvgTree (Tree) -import Linear.V2 -import Reanimate -import Reanimate.Chiphunk -import Reanimate.PolyShape -import System.IO.Unsafe - -test :: Animation -test = unsafePerformIO $ do - bodyStore <- newBodyStore - let gravity = Vect 0 (-1) - - -- Create an empty space. - space <- spaceNew - spaceCollisionSlop space $= (screenWidth/2560) - spaceGravity space $= gravity - - -- Add a static line segment shape for the ground. - -- We'll make it slightly tilted so the ball will roll off. - -- We attach it to a static body to tell Chipmunk it shouldn't be movable. - static <- get $ spaceStaticBody space - ground <- segmentShapeNew static - (Vect (-screenWidth/2) 0) - (Vect (screenWidth/2) (-screenHeight/2)) 0 - -- ground <- polyShapeNewRaw static - -- [ Vect (-screenWidth/2) (screenHeight/2) - -- , Vect (screenWidth/2) (-screenHeight/2) - -- , Vect (-screenWidth/2) (-screenHeight/2) - -- , Vect (-screenWidth/2) (screenHeight/2) ] 0 - shapeFriction ground $= 1 - spaceAddShape space ground - - - - -- Now let's make a ball that falls onto the line and rolls off. - -- First we need to make a cpBody to hold the physical properties of the object. - -- These include the mass, position, velocity, angle, etc. of the object. - -- Then we attach collision shapes to the Body to give it a size and shape. - - let toVect (V2 x y) = Vect x y - - let svg = center $ scale 2 $ latex "\\LaTeX" - poly = svgToPolyShapes svg - vectGroup = plDecompose poly - --mkCircle (Num radius) - -- vects = svgToVects svg - -- vects' = fst $ convexHull vects 0 - -- svg' = - -- withFillOpacity 0 $ withStrokeColor "white" $ - -- withStrokeWidth 0.01 $ - -- renderPolyShapes (map plFromPolygon vectGroup) - -- - -- ballBody <- polyShapesToBody space poly - -- bodyPosition ballBody $= Vect (-screenWidth/4) (screenHeight/2) - -- - -- addToBodyStore bodyStore ballBody $ - -- withFillColor "white" $ - -- mkGroup - -- [ svg' ] - - -- let splitPolys = vectGroup - forM_ vectGroup $ \polygon -> do - bd <- polygonsToBody space [map toVect polygon] - bodyPosition bd $= Vect 0 (screenHeight/3) - addToBodyStore bodyStore bd $ - renderPolyShape $ plFromPolygon polygon - -- withFillColor "white" $ - -- mkGroup - -- [ --withStrokeWidth (Num 0.005) $ - -- -- withStrokeWidth (Num 0.00) $ - -- -- withStrokeColor "white" $ - -- -- withFillOpacity 1 $ - -- renderPolyShape $ plFromPolygon polygon - -- ] - - ani <- simulate space bodyStore 60 60 10 - spaceFreeRecursive space - return ani - --- data LineCommand --- = LineMove RPoint --- -- | LineDraw RPoint --- | LineBezier [RPoint] --- | LineEnd - --- vectsToSVG :: [Vect] -> Tree --- vectsToSVG (Vect x y:rest) = --- mkPath $ --- MoveTo OriginAbsolute [V2 x y] : --- [ LineTo OriginAbsolute [V2 a b] | Vect a b <- rest ] ++ --- [ EndPath ] --- where --- mkPath cmds = PathTree $ defaultSvg & pathDefinition .~ cmds - --- polygonsToSVG :: [[Vect]] -> Tree --- polygonsToSVG = merge . mkGroup . map vectsToSVG --- where --- merge svg = PathTree $ defaultSvg & pathDefinition .~ extractPath svg - --- svgToVects :: Tree -> [Vect] --- svgToVects svg = map worker (lineToPoints 200 cmds) --- where --- worker (V2 x y) = Vect x y --- cmds = toLineCommands $ wibble $ extractPath svg --- wibble xs = takeWhile (/=EndPath) xs ++ [EndPath] - -chunkPolyshapes :: Tree -> Tree -chunkPolyshapes t = - withStrokeColor "white" $ - withStrokeWidth 0.01 $ - withFillColor "white" $ t - -reorient :: Tree -> Tree -reorient = id -- scale 4 . translate 0 (-0.9) - -main :: IO () -main = reanimate $ bg `parA` mapA reorient (line `parA` mapA chunkPolyshapes test) - where - bg = animate $ const $ mkBackground "black" - line = animate $ const $ withStrokeColor "white" $ - withStrokeWidth 0.01 $ - mkLine (-screenWidth/2, 0) - (screenWidth/2, -screenHeight/2) diff --git a/examples/simulate_equation.hs b/examples/simulate_equation.hs deleted file mode 100755 index e0651cb..0000000 --- a/examples/simulate_equation.hs +++ /dev/null @@ -1,189 +0,0 @@ -#!/usr/bin/env stack --- stack runghc --package reanimate -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE RecordWildCards #-} -module Main (main) where - -import Chiphunk.Low -import Control.Monad -import Graphics.SvgTree (Tree) -import Linear.V2 -import Reanimate.Chiphunk -import Reanimate -import Reanimate.PolyShape -import System.IO.Unsafe - -test :: Animation -test = unsafePerformIO $ do - bodyStore <- newBodyStore - let gravity = Vect 0 (-1) - - -- Create an empty space. - space <- spaceNew - -- spaceCollisionSlop space $= (screenWidth/2560) - spaceGravity space $= gravity - - -- Add a static line segment shape for the ground. - -- We'll make it slightly tilted so the ball will roll off. - -- We attach it to a static body to tell Chipmunk it shouldn't be movable. - static <- get $ spaceStaticBody space - ground <- segmentShapeNew static - (Vect (-screenWidth/2) (screenHeight/2)) - (Vect (screenWidth/2) (-screenHeight/2)) 0 - shapeFriction ground $= 1 - spaceAddShape space ground - - -- Now let's make a ball that falls onto the line and rolls off. - -- First we need to make a cpBody to hold the physical properties of the object. - -- These include the mass, position, velocity, angle, etc. of the object. - -- Then we attach collision shapes to the Body to give it a size and shape. - - let toVect (V2 x y) = Vect x y - - let svg = center $ scale 2 $ latex "$\\Phi$" - poly = svgToPolyShapes svg - vectGroup = plDecompose' 1 poly - --mkCircle (Num radius) - -- vects = svgToVects svg - -- vects' = fst $ convexHull vects 0 - -- svg' = - -- withFillOpacity 0 $ withStrokeColor "white" $ - -- withStrokeWidth 0.01 $ - -- renderPolyShapes (map plFromPolygon vectGroup) - -- - -- ballBody <- polyShapesToBody space poly - -- bodyPosition ballBody $= Vect (-screenWidth/4) (screenHeight/2) - -- - -- addToBodyStore bodyStore ballBody $ - -- withFillColor "white" $ - -- mkGroup - -- [ svg' ] - - -- let splitPolys = vectGroup - forM_ vectGroup $ \polygon -> do - bd <- polygonsToBody space [map toVect polygon] - bodyPosition bd $= Vect 0 (screenHeight/3) - addToBodyStore bodyStore bd $ - renderPolyShape $ plFromPolygon polygon - -- withFillColor "white" $ - -- mkGroup - -- [ --withStrokeWidth (Num 0.005) $ - -- -- withStrokeWidth (Num 0.00) $ - -- -- withStrokeColor "white" $ - -- -- withFillOpacity 1 $ - -- renderPolyShape $ plFromPolygon polygon - -- ] - - ani <- simulate space bodyStore 60 10 10 - spaceFreeRecursive space - return ani - --- data LineCommand --- = LineMove RPoint --- -- | LineDraw RPoint --- | LineBezier [RPoint] --- | LineEnd - --- vectsToSVG :: [Vect] -> Tree --- vectsToSVG (Vect x y:rest) = --- mkPath $ --- MoveTo OriginAbsolute [V2 x y] : --- [ LineTo OriginAbsolute [V2 a b] | Vect a b <- rest ] ++ --- [ EndPath ] --- where --- mkPath cmds = PathTree $ defaultSvg & pathDefinition .~ cmds - --- polygonsToSVG :: [[Vect]] -> Tree --- polygonsToSVG = merge . mkGroup . map vectsToSVG --- where --- merge svg = PathTree $ defaultSvg & pathDefinition .~ extractPath svg - --- svgToVects :: Tree -> [Vect] --- svgToVects svg = map worker (lineToPoints 200 cmds) --- where --- worker (V2 x y) = Vect x y --- cmds = toLineCommands $ wibble $ extractPath svg --- wibble xs = takeWhile (/=EndPath) xs ++ [EndPath] - -chunkPolyshapes :: Tree -> Tree -chunkPolyshapes t = - mkGroup - [ - -- withFillColor "white" $ - -- withFillOpacity 1 $ - -- withStrokeWidth (Num 0.01) $ - -- withStrokeColor "blue" $ - -- mkGroup $ map renderPolyShape $ - -- map mergePolyShapeHoles $ - -- plGroupShapes $ - -- -- unionPolyShapes $ - -- traceUnion $ - -- svgToPolyShapes t - withStrokeWidth 0.001 $ - withStrokeColor "red" $ - withFillOpacity 0 $ - mkGroup $ map renderPolyShape $ - svgToPolyShapes t - -- withStrokeWidth 0.001 $ - -- withStrokeColor "red" $ - -- mkGroup $ map renderPolyShape [pl12,pl13] - ] - -reorient :: Tree -> Tree -reorient = scale 6 . translate 0 (-0.9) - --- pl1 = PolyShape {unPolyShape = ClosedPath [(Point 9.076152025049407e-2 0.8793158714390639,JoinLine),(Point 0.48865783997158957 0.5875446050720555,JoinLine),(Point 0.6291445682322502 0.560956966929786,JoinLine),(Point 0.684414274089106 0.6254399758525935,JoinLine),(Point 8.549729225029659e-2 0.9226054512958702,JoinLine)]} --- pl2 = PolyShape {unPolyShape = ClosedPath [(Point 0.3174796674934045 (-0.17830385091116296),JoinLine),(Point 0.3930521201476444 (-0.22384314349780887),JoinLine),(Point 0.42400439387727545 (-0.17026635397131717),JoinLine),(Point 0.367180526197663 (-0.13743818486413273),JoinLine)]} --- pl3 = PolyShape {unPolyShape = ClosedPath [(Point 3.818529643334073e-2 4.994390978944008e-2,JoinLine),(Point 3.5509598929490127e-3 3.0061273354143103e-2,JoinLine),(Point 2.6643683455205014e-2 (-1.3314171733262437e-2),JoinLine)]} --- pl4 = PolyShape {unPolyShape = ClosedPath [(Point 8.041123459068245e-2 (-4.592522065745819e-2),JoinLine),(Point 0.29675416515553454 (-0.17235110729572378),JoinLine),(Point 0.34590008699238584 (-0.13081971694800998),JoinLine),(Point 0.22737727571852562 4.513304505856075e-2,JoinLine),(Point 0.10498823211784264 (-3.373201779345303e-3),JoinLine)]} --- pl5 = PolyShape {unPolyShape = ClosedPath [(Point 4.1547772506933034e-2 0.10908491253751296,JoinLine),(Point 5.7556283592893165e-2 (-2.2936657573235042e-2),JoinLine),(Point 0.18645395693905423 2.887370173550552e-2,JoinLine),(Point 0.24063435970720265 0.12675372820097613,JoinLine),(Point 0.19807899284556882 0.17285980812856805,JoinLine),(Point 8.45668830191807e-2 0.18680141998908184,JoinLine)]} --- pl6 = PolyShape {unPolyShape = ClosedPath [(Point (-0.2574090456009813) 0.22352645168090446,JoinLine),(Point (-0.27139525022076605) 0.16073053060077958,JoinLine),(Point (-9.444562231853837e-2) 4.994785063168361e-2,JoinLine),(Point (-4.012849544398306e-2) 6.733655689772311e-2,JoinLine),(Point (-4.6550067413171026e-2) 0.20017002713448606,JoinLine)]} --- pl7 = PolyShape {unPolyShape = ClosedPath [(Point (-0.3714373753139388) 0.20874389306710772,JoinLine),(Point (-0.29452803703157854) 0.1655072120141602,JoinLine),(Point (-0.2835217605861194) 0.22889335799870397,JoinLine),(Point (-0.34012040364470725) 0.2621083279509192,JoinLine)]} --- pl8 = PolyShape {unPolyShape = ClosedPath [(Point 0.1712163207191154 0.2149766965620984,JoinLine),(Point 0.19616659305280632 0.1905738263277067,JoinLine),(Point 0.6315796735654418 0.5328570813993471,JoinLine),(Point 0.49895896036053317 0.5862908698518825,JoinLine)]} --- pl9 = PolyShape {unPolyShape = ClosedPath [(Point 0.10739670009207455 0.19390054402859663,JoinLine),(Point 0.16776933204268935 0.17681618959615353,JoinLine),(Point 0.19166035517967667 0.1970425575035583,JoinLine),(Point 0.16584071117957827 0.2205236594561943,JoinLine)]} --- pl10 = PolyShape {unPolyShape = ClosedPath [(Point (-0.3669271639958152) 1.1382616686769857,JoinLine),(Point (-0.30169633405804164) 1.1454435521784956,JoinLine),(Point (-0.2858326377769962) 1.2078015601462844,JoinLine),(Point (-0.3736986541543818) 1.199765022618315,JoinLine)]} --- pl11 = PolyShape {unPolyShape = ClosedPath [(Point (-3.0061173114158422e-2) 1.1779482589541495,JoinLine),(Point 5.4907529076574785e-3 1.2315286122957358,JoinLine),(Point (-3.4177530164954735e-2) 1.2269151735259667,JoinLine)]} --- pl12 = PolyShape {unPolyShape = ClosedPath --- [(Point (-0.301730083226429*10) 1.1457103198223786,JoinLine) --- ,(Point (-0.10299641326632633*10) 1.0714677466302471,JoinLine) --- ,(Point (-3.162369257155477e-2*10) 1.182092647608354,JoinLine) --- ,(Point (-3.636877035247056e-2*10) 1.2310026399913667,JoinLine) --- ,(Point (-0.28589249688644447*10) 1.208074964264014,JoinLine)]} --- pl13 = PolyShape {unPolyShape = ClosedPath --- [(Point (-9.39646610882567e-2*10) 0.96018603098516,JoinLine) --- ,(Point (-6.465077411411442e-2*10) 0.9514465910530718,JoinLine) --- ,(Point (0.14134955357291654*10) 1.2461928296570437,JoinLine) --- ,(Point (8.29425965486209e-3*10) 1.2330682613363688,JoinLine) --- ,(Point (-0.10287959610529117*10) 1.0717052650151784,JoinLine)]} --- pl14 = PolyShape {unPolyShape = ClosedPath [(Point (-0.6410657779553459) 0.5634659418454123,JoinLine),(Point (-0.641345175626591) 0.555247219163176,JoinLine),(Point (-0.4420444191709594) 0.5513358399243028,JoinLine),(Point (-0.10594072095664082) 0.9023816092645105,JoinLine),(Point (-7.577913627887163e-2) 0.9355404742795719,JoinLine),(Point (-0.10373629906690308) 0.9479533236181756,JoinLine)]} --- pl15 = PolyShape {unPolyShape = ClosedPath [(Point (-3.523485505752828e-2) 0.1998408875682416,JoinLine),(Point 1.2546982880291841e-2 0.15208851089779252,JoinLine),(Point 0.13018455413761065 0.208301106923033,JoinLine),(Point 8.433578523127425e-2 0.9418697065712005,JoinLine),(Point (-7.884702019521485e-2) 0.932428903501462,JoinLine),(Point (-8.084970167352624e-2) 0.9296667902794287,JoinLine)]} --- pl16 = PolyShape {unPolyShape = ClosedPath [(Point (-7.878197627394322e-2) 0.9319709148544943,JoinLine),(Point 8.438182768423004e-2 0.941734614325483,JoinLine),(Point 8.518191684165782e-2 0.9853357574063363,JoinLine),(Point 8.15663333974107e-2 0.9853939259545145,JoinLine),(Point (-4.936159754845773e-2) 0.9727171460725276,JoinLine)]} --- pl17 = PolyShape {unPolyShape = ClosedPath [(Point (-4.99911351869741e-2) 0.9726181699897635,JoinLine),(Point 8.094809532371453e-2 0.9851776987711782,JoinLine),(Point 7.384845354621125e-2 1.0968271981504014,JoinLine),(Point 4.400324523877765e-2 1.1025518925313644,JoinLine)]} --- pl18 = PolyShape {unPolyShape = ClosedPath [(Point 4.8482888984171635e-2 1.1105768787999895,JoinLine),(Point 7.512899062784217e-2 1.0959654499901599,JoinLine),(Point 0.2822087302816073 1.1420633097553252,JoinLine),(Point 0.27513287205836434 1.2060076125068004,JoinLine),(Point 0.20204313118994668 1.2231219927110881,JoinLine)]} --- pl19 = PolyShape {unPolyShape = ClosedPath [(Point 0.2748784478150325 1.2078902086770202,JoinLine),(Point 0.2820822061918643 1.1439601884266248,JoinLine),(Point 0.3457225261391361 1.127943000720636,JoinLine),(Point 0.36082444597621127 1.187946730956635,JoinLine)]} --- pl20 = PolyShape {unPolyShape = ClosedPath [(Point (-0.6864531908080556) 0.5032767356555317,JoinLine),(Point (-9.929294916146593e-2) 0.20829491627795255,JoinLine),(Point (-8.483230342755331e-2) 0.23523598173633797,JoinLine),(Point (-0.49217693932265094) 0.547917843618004,JoinLine)]} --- pl21 = PolyShape {unPolyShape = ClosedPath [(Point (-9.274301086931341e-2) 0.2153671046631685,JoinLine),(Point (-4.68781647762479e-2) 0.2005165600411279,JoinLine),(Point (-4.790067311086227e-2) 0.23865704773924823,JoinLine),(Point (-7.795810209993242e-2) 0.2330461435659677,JoinLine)]} - --- unionTest = parsed --- where --- svg = center $ scale 2 $ latex "I" --- poly = svgToPolyShapes svg --- vectGroup = plDecompose' 1 poly --- rendered = mkGroup $ take 100 $ map (renderPolyShape . plFromPolygon) vectGroup --- parsed = svgToPolyShapes rendered --- --- myDecompose tol = --- concatMap decomposePolygon . --- map (plPolygonify tol) . --- map mergePolyShapeHoles . --- plGroupShapes . --- unionPolyShapes - -main :: IO () -main = reanimate $ bg `parA` mapA reorient (line `parA` mapA chunkPolyshapes test) - where - bg = animate $ const $ mkBackground "black" - line = animate $ const $ withStrokeColor "white" $ - withStrokeWidth 0.01 $ - mkLine (-screenWidth/2, screenHeight/2) - (screenWidth/2, -screenHeight/2) diff --git a/examples/simulate_gravity.golden b/examples/simulate_gravity.golden deleted file mode 100644 index 7c8e128..0000000 --- a/examples/simulate_gravity.golden +++ /dev/null @@ -1,10 +0,0 @@ -0 -1 -2 -3 -4 -5 -6 -7 -8 -9 diff --git a/examples/simulate_gravity.hs b/examples/simulate_gravity.hs deleted file mode 100755 index 6af3fed..0000000 --- a/examples/simulate_gravity.hs +++ /dev/null @@ -1,79 +0,0 @@ -#!/usr/bin/env stack --- stack runghc --package reanimate -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE RecordWildCards #-} -module Main (main) where - -import Chiphunk.Low -import Reanimate.Chiphunk -import Reanimate -import System.IO.Unsafe - - -test :: Animation -test = unsafePerformIO $ do - bodyStore <- newBodyStore - let gravity = Vect 0 (-10) - - -- Create an empty space. - space <- spaceNew - spaceGravity space $= gravity - - -- Add a static line segment shape for the ground. - -- We'll make it slightly tilted so the ball will roll off. - -- We attach it to a static body to tell Chipmunk it shouldn't be movable. - static <- get $ spaceStaticBody space - ground <- segmentShapeNew static - (Vect (-screenWidth/2) 0) - (Vect (screenWidth/2) (-screenHeight/2)) 0 - shapeFriction ground $= 1 - spaceAddShape space ground - - -- Now let's make a ball that falls onto the line and rolls off. - -- First we need to make a cpBody to hold the physical properties of the object. - -- These include the mass, position, velocity, angle, etc. of the object. - -- Then we attach collision shapes to the Body to give it a size and shape. - - let radius = 1 - let mass = 1 - - -- The moment of inertia is like mass for rotation - -- Use the momentFor* functions to help you approximate it. - let moment = momentForCircle mass 0 radius (Vect 0 0) - - -- The spaceAdd* functions return the thing that you are adding. - ballBody <- bodyNew mass moment - spaceAddBody space ballBody - bodyPosition ballBody $= Vect 0 (screenHeight/2) - - -- Now we create the collision shape for the ball. - -- You can create multiple collision shapes that point to the same body. - -- They will all be attached to the body and move around to follow it. - ballShape <- circleShapeNew ballBody radius (Vect 0 0) - spaceAddShape space ballShape - shapeFriction ballShape $= 0.7 - - addToBodyStore bodyStore ballBody $ - withFillColor "white" $ - mkGroup - [ mkCircle radius - , withStrokeColor "black" $ - mkLine (0, 0) (0, radius) ] - - ani <- simulate space bodyStore 60 3 4 - - shapeFree ballShape - bodyFree ballBody - shapeFree ground - spaceFree space - return ani - - -main :: IO () -main = reanimate $ bg `parA` line `parA` test - where - bg = animate $ const $ mkBackground "black" - line = animate $ const $ withStrokeColor "white" $ - withStrokeWidth 0.01 $ - mkLine (-screenWidth/2, 0) - (screenWidth/2, -screenHeight/2) diff --git a/examples/tut_glue_physics.hs b/examples/tut_glue_physics.hs deleted file mode 100755 index 39e1dfc..0000000 --- a/examples/tut_glue_physics.hs +++ /dev/null @@ -1,59 +0,0 @@ -#!/usr/bin/env stack --- stack runghc --package reanimate -{-# LANGUAGE OverloadedStrings #-} -{-# LANGUAGE RecordWildCards #-} -module Main (main) where - -import Chiphunk.Low -import Codec.Picture (PixelRGBA8 (..)) -import Control.Lens -import Control.Monad (forM_) -import Linear.V2 (V2 (..)) -import Reanimate -import Reanimate.Chiphunk -import Reanimate.PolyShape -import System.IO.Unsafe (unsafePerformIO) - -shatter :: Animation -shatter = unsafePerformIO $ do - bodyStore <- newBodyStore - let gravity = Vect 0 (-1) -- Gravity points down (negative 1 y/s^2) - - space <- spaceNew - spaceGravity space $= gravity - - static <- get $ spaceStaticBody space - ground <- segmentShapeNew static - (Vect (-screenWidth/2) (-screenHeight/2)) - (Vect (screenWidth/2) (-screenHeight/2)) 0 - shapeFriction ground $= 1 - spaceAddShape space ground - - let toVect (V2 x y) = Vect x y - vectGroup = plDecompose $ svgToPolyShapes $ center $ scale 4 $ - latex "$F=ma$" - - forM_ vectGroup $ \polygon -> do - bd <- polygonsToBody space [map toVect polygon] - bodyPosition bd $= Vect 0 (screenHeight/3) - addToBodyStore bodyStore bd $ - renderPolyShape $ plFromPolygon polygon - - ani <- simulate space bodyStore fps stepsPerFrame shatterT - spaceFreeRecursive space - return $ mapA pp ani - where - shatterT = 10 - fps = 60 - stepsPerFrame = 10 - pp = withStrokeWidth 0.01 . withStrokeColor "black" . withFillColor "black" - -main :: IO () -main = reanimate $ parA bg $ sceneAnimation $ do - play $ shatter - play $ shatter - & reverseA - & setDuration 5 - & signalA (powerS 2) - where - bg = animate $ const $ mkBackgroundPixel (PixelRGBA8 252 252 252 0xFF) diff --git a/mkdocs.yml b/mkdocs.yml index fdec5b7..82a9742 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -19,8 +19,8 @@ nav: - Reanimate: index.md - Core concepts: introduction.md - Gluing together animations: glue_tut.md -# - Morphology (incomplete): morphology.md - - Voice control: voice.md + - Morphology (incomplete): morphology.md +# - Voice control: voice.md # - How to: # - Use fonts (incomplete): fonts.md - Showcase: diff --git a/examples/morphology_leastwork.hs b/reanimate-examples/morphology_leastwork.hs similarity index 100% rename from examples/morphology_leastwork.hs rename to reanimate-examples/morphology_leastwork.hs diff --git a/examples/morphology_linebend.hs b/reanimate-examples/morphology_linebend.hs similarity index 100% rename from examples/morphology_linebend.hs rename to reanimate-examples/morphology_linebend.hs diff --git a/examples/morphology_linebend_intro.hs b/reanimate-examples/morphology_linebend_intro.hs similarity index 100% rename from examples/morphology_linebend_intro.hs rename to reanimate-examples/morphology_linebend_intro.hs diff --git a/src/Reanimate/Math/Compatible.hs b/reanimate-morph/Compatible.hs similarity index 54% rename from src/Reanimate/Math/Compatible.hs rename to reanimate-morph/Compatible.hs index 90a1c84..72f0677 100644 --- a/src/Reanimate/Math/Compatible.hs +++ b/reanimate-morph/Compatible.hs @@ -1,31 +1,37 @@ -{-# LANGUAGE MultiWayIf #-} +{-# 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 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.Debug -import Reanimate.Math.Render import Reanimate.Svg -import Reanimate.Animation -- import Debug.Trace +-- import GHC.Stack -truncateP :: V2 Rational -> V2 Rational -truncateP = fmap (realToFrac . (realToFrac :: Rational -> Double)) +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 :: V2 Rational -> V2 Rational -> Int -> [V2 Rational] -mkSteinerPoints a b s_ = [ lerp (i / (s + 1)) b a | i <- [1 .. s] ] - where s = fromIntegral s_ +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 :: Polygon -> Int -> Int -> Int -> (Polygon, Polygon) +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) @@ -35,7 +41,7 @@ split1Link p i j s = left = map (pAccess p) [0 .. i] ++ sp ++ map (pAccess p) [j .. n - 1] right = map (pAccess p) [i .. j] ++ reverse sp -steiner2Link :: Polygon -> Int -> Int -> V2 Rational +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 @@ -46,57 +52,94 @@ steiner2Link 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 - = truncateP $ lerp 0.5 (fst vect) (intersects !! 0) + = 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 - distToV = approxDist (fst vect) + -- 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 vect - : [ u - | n <- [0 .. pSize p - 1] - , let edge = (pAccess p n, pAccess p $ pNext p n) - , u <- case rayIntersect vect edge of - Nothing -> [] - Just u -> [u] - , isBetween u edge - , u /= fst vect - , isForward vect u - ] - iP = pAdjustOffset p i - jP = pAdjustOffset p j - vect - | 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, p2) - | otherwise - = fromMaybe - (error $ "No window overlap: " ++ show - (isStraightLine, isGrandparent, oneBendBetween p i j) - ) - $ listToMaybe - $ [ (p1, p1 + (p2 - p1) + (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)) + -- 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 :: Polygon -> Int -> Int -> (Polygon, Polygon) +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 @@ -107,7 +150,7 @@ split2Link p i j = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right) data Link = OneLink | TwoLink -splitNLink :: Polygon -> Int -> [(Link, Int)] -> (Polygon, Polygon) +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 @@ -117,7 +160,7 @@ splitNLink p i js = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right) steiners = splitNLinks p i js -- 0, [(TwoLink,2),(OneLink,3)] -- 0, [(OneLink,5),(TwoLink,3)] -splitNLinks :: Polygon -> Int -> [(Link, Int)] -> [V2 Rational] +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)] = [] @@ -154,7 +197,7 @@ splitNLinks p i ((OneLink, j) : (TwoLink, j') : xs) = s' : s : splitNLinks p'' sIdx' ((OneLink, j') : xs) splitNLinks _p _i _ = error "splitNLinks: invalid input" -selectContains :: Polygon -> Polygon -> V2 Rational -> Polygon +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" @@ -218,7 +261,7 @@ data MeshPair = MeshPair Points Points Edges -- (a,b,1) -- (a,b,2) -- floyd -compatiblyTriangulateP :: Polygon -> Polygon -> [(Polygon, Polygon)] +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' @@ -226,39 +269,57 @@ compatiblyTriangulateP a b (pSetOffset a 0) (pSetOffset b 0) -showStep :: Polygon -> Polygon -> SVG -showStep a b = mkGroup - [ translate (-3) 0 - $ mkGroup [withFillColor "grey" $ polygonShape a, polygonNumDots a] - , translate 3 0 - $ mkGroup [withFillColor "grey" $ polygonShape b, polygonNumDots b] - ] +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] + -- ] -compatiblyTriangulateP' :: Polygon -> Polygon -> Polygon -> [(Polygon, Polygon)] +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 - (show - ("two link" :: String, toOriginIndex nodeL, toOriginIndex nodeR) - ) - $ -} - let (aL, aR) = if (nodeL, nodeR) `elem` aOneLink + 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 compatiblyTriangulateP' aOrigin aL bL - ++ compatiblyTriangulateP' aOrigin aR bR + 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) @@ -285,22 +346,21 @@ compatiblyTriangulateP' aOrigin a b bOneLink = polygonOneLinks b aTwoLink = polygonTwoLinks a bTwoLink = polygonTwoLinks b - nodeDist (i, j) = min (j - i) (n - j + i) + nodeDist (i, j) = (min (j - i) (n - j + i), hash (i,j)) -oneBendBetween :: Polygon -> Int -> Int -> Bool -oneBendBetween _p _a _b = False --- oneBendBetween p a b = --- direction --- (pAccess p b) --- (pAccess p (pParent p a b)) --- (pAccess p $ obstructedBy b) == 0 && --- (pParent p a b) /= obstructedBy b --- where --- obstructedBy n = --- case pParent p a n of --- i -> if i == a then n else obstructedBy i +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 :: Polygon -> [(Int, Int)] +polygonTwoLinks :: (Real a, Fractional a) => APolygon a -> [(Int, Int)] polygonTwoLinks p = [ (i, j) | i <- [0 .. n - 1] @@ -315,7 +375,7 @@ polygonTwoLinks p = ] where n = pSize p -polygonOneLinks :: Polygon -> [(Int, Int)] +polygonOneLinks :: APolygon a -> [(Int, Int)] polygonOneLinks p = [ (i, j) | i <- [0 .. pSize p - 1] diff --git a/reanimate-morph/LeastDifference.hs b/reanimate-morph/LeastDifference.hs new file mode 100644 index 0000000..1d71731 --- /dev/null +++ b/reanimate-morph/LeastDifference.hs @@ -0,0 +1,195 @@ +{-# LANGUAGE OverloadedStrings #-} +module Reanimate.Morph.LeastDifference where + +import Control.Applicative +import Control.Monad +import Data.List +import Data.Ord +import Data.Maybe +import Linear.V2 +import qualified Data.Vector as V +import Reanimate.Math.Common +import qualified Reanimate.Math.Compatible as Compat +import Reanimate.Math.Polygon +import Reanimate.Morph.Common + +import Debug.Trace + +leastDifference :: PointCorrespondence +leastDifference = undefined + +pCutsTop :: (Real a, Fractional a) => [(V2 a, V2 a)] -> APolygon a -> [(APolygon a, APolygon a, [Topology], [Topology], (V2 a, V2 a))] +pCutsTop edges p = + [ (l, r, topL, topR, edge) + | (l, r) <- pCuts p + , let edge = (pAccess l 0, pAccess r 1) + topL = pTopology (edge : edges) l + topR = pTopology (edge : edges) r + ] + +{-# INLINE pCompatibleCut #-} +pCompatibleCut :: (Real a, Fractional a) => [(Edge a,Edge a)] -> APolygon a -> APolygon a -> Maybe ((APolygon a, APolygon a),(APolygon a, APolygon a), [(Edge a,Edge a)]) +pCompatibleCut edges p1 p2 = listToMaybe $ sortOn pairSizeDifference $ do + (p1l, p1r, aTopL, aTopR, aEdge) <- pCutsTop (map fst edges) p1 + -- let aEdge = (pAccess p1l 0, pAccess p1r 1) + -- aTopL = pTopology (aEdge : map fst edges) p1l + -- aTopR = pTopology (aEdge : map fst edges) p1r + (p2l, p2r, bTopL, bTopR, bEdge) <- pCutsTop (map snd edges) p2 + -- let bEdge = (pAccess p2r 0, pAccess p2r 1) + -- bTopL = pTopology (bEdge : map snd edges) p2l + -- bTopR = pTopology (bEdge : map snd edges) p2r + let newEdges = (aEdge,bEdge) : edges + -- guard $ pTopology (map fst newEdges) p1l == pTopology (map snd newEdges) p2l + -- guard $ pTopology (map fst newEdges) p1r == pTopology (map snd newEdges) p2r + guard $ aTopL == bTopL + guard $ aTopR == bTopR + return ((p1l, p1r), (p2l, p2r), newEdges) + +sizeDifference :: (Fractional a) => (APolygon a, APolygon a) -> a +sizeDifference (a,b) = abs (pArea a - pArea b) + +pairSizeDifference :: (Fractional a, Ord a) => ((APolygon a, APolygon a),(APolygon a, APolygon a),b) -> a +pairSizeDifference (a,b,_) = max (sizeDifference a) (sizeDifference b) + +pEdges :: APolygon a -> [Edge a] +pEdges p = [ (pAccess p i, pAccess p $ i+1) | i <- [0 .. pSize p-1]] + +type Edge a = (V2 a, V2 a) + +data Topology = MutableSegment | ImmutableSegment Int + deriving (Eq, Show) + +{- INLINE pTopology -} +pTopology :: Eq a => [Edge a] -> APolygon a -> [Topology] +pTopology immutableEdges p = worker + [ maybe MutableSegment ImmutableSegment isImmutable + | i <- [0 .. pSize p-1] + , let isImmutable = + ((pAccess p i, pAccess p $ i+1) `elemIndex` immutableEdges) <|> + ((pAccess p $ i+1, pAccess p i) `elemIndex` immutableEdges) + ] + where + worker (MutableSegment:MutableSegment:xs) = worker (MutableSegment:xs) + worker (x:xs) = x : worker xs + worker [] = [] + +triangulate :: Polygon -> Polygon -> [(Polygon,Polygon)] +triangulate a b = + [ (castPolygon l, castPolygon r) + | (l,r) <- triangulate_ (castPolygon a) (castPolygon b) + ] + +triangulate_ :: FPolygon -> FPolygon -> [(FPolygon,FPolygon)] +triangulate_ a b + | distSquared (pCentroid al) (pCentroid bl) + + distSquared (pCentroid ar) (pCentroid br) < + distSquared (pCentroid al) (pCentroid br) + + distSquared (pCentroid ar) (pCentroid bl) + -- = [(al,bl), (ar,br)] + = triangulate' edges al bl ++ triangulate' edges ar br + | otherwise + -- = [(al,br), (ar,bl)] + = triangulate' edges al br ++ triangulate' edges ar bl + where + edges = [(aEdge,bEdge)] + aEdge = (pAccess ar 0, pAccess ar 1) + bEdge = (pAccess br 0, pAccess br 1) + (al, ar) = pCutEqual a + (bl, br) = pCutEqual b + +{-# SPECIALIZE triangulate' :: [(Edge Double,Edge Double)] -> FPolygon -> FPolygon -> [(FPolygon,FPolygon)] #-} +triangulate' :: (Real a, Fractional a) => [(Edge a,Edge a)] -> APolygon a -> APolygon a -> [(APolygon a,APolygon a)] +triangulate' edges a b + | pSize a == 3 || pSize b == 3 = giveUp + -- | trace (show (pSize a, pSize b)) False = undefined + | Just ((al,ar),(bl,br), newEdges) <- pCompatibleCut edges a b + = -- traceSVG (helper "topology") $ + triangulate' newEdges al bl ++ triangulate' newEdges ar br + -- = [(al,bl), (ar,br)] + | otherwise = giveUp + where + giveUp = [(a,b)] + -- let aNewEdges = max 0 (pSize b - pSize a) + -- bNewEdges = max 0 (pSize a - pSize b) + -- a' = pAddPointsRestricted (map fst edges) aNewEdges a + -- b' = pAddPointsRestricted (map snd edges) bNewEdges b + -- -- (a'', b'') = closestLinearCorrespondenceA a' b' + -- in [(a, b)] + --Compat.compatiblyTriangulateP a'' b'' + +polygonLengths :: Eq a => [APolygon a] -> APolygon a -> [(Int, Int, Int)] +polygonLengths polys origin' = worker 0 + where + Just firstIdx = V.elemIndex (pAccess (head polys) 0) (polygonPoints origin') + origin = pSetOffset origin' firstIdx + polyLength a b + | V.elem (pAccess a 1) (polygonPoints b) + = 1 + polyLength (pAdjustOffset a 1) b + | otherwise + = 0 + findPoly idx a b = fromMaybe (error $ "Missing: " ++ show idx) $ listToMaybe + [ (n, p, (i+polygonOffset p) `mod` pSize p) + | (n, p) <- zip [0..] polys + , let pts = polygonPoints p + , Just i <- [V.elemIndex a pts] + , V.elem b pts + ] + worker n | n == pSize origin = [] + worker n = + let (polyIdx, poly, polyKey) = findPoly n (pAccess origin n) (pAccess origin $ n+1) + originHere = pAdjustOffset origin n + len = polyLength originHere poly + in (polyIdx, polyKey, len) : worker (n+len) + + +alignPolygons :: (Fractional a, Real a, Eq a) => + [(APolygon a, APolygon a)] -> APolygon a -> APolygon a -> [(APolygon a, APolygon a)] +alignPolygons polys originL originR = + zip + (worker 0 (map fst polys) diffsL) + (worker 0 (map snd polys) diffsR) + where + lensL = polygonLengths (map fst polys) originL + lensR = polygonLengths (map snd polys) originR + diffsL = sortOn sortKey $ zipWith mkDiff lensL lensR + diffsR = sortOn sortKey $ zipWith mkDiff lensR lensL + worker n (y:ys) ((i,k,l,d):xs) | n < i + = y : worker (n+1) ys ((i,k,l,d):xs) + worker n (y:ys) ((_i,_k,_l,0):xs) + = worker n (y:ys) xs + worker n (y:ys) ((_i,k,l,d):xs) + = let y' = pAddPointsBetween (k,l) d y + in worker n (y':ys) xs + worker _ ys [] = ys + worker _ _ _ = error "bad worker input" + mkDiff (li, lk, ll) (_ri, _rk, rl) = + (li, lk, ll, max 0 (rl-ll)) + sortKey (i, k, _l, _d) = (i,Down k) + +compatTriagPairs :: (Real a, Fractional a) => [(APolygon a, APolygon a)] -> [(APolygon a, APolygon a)] +compatTriagPairs = concatMap (uncurry Compat.compatiblyTriangulateP) + +circumference :: (Real a, Fractional a) => [APolygon a] -> APolygon a -> APolygon a +circumference polys origin' = mkPolygon $ V.fromList $ worker 0 + where + Just firstIdx = V.elemIndex (pAccess (head polys) 0) (polygonPoints origin') + origin = pSetOffset origin' firstIdx + fromTo i j p + | i < j = [i .. j-1] + | otherwise = [i .. pSize p-1] ++ [0..j-1] + check [] = error "no results" + check [x] = x + check _ = error "multiple results" + findPoints n a b = check -- fromMaybe (error $ "Missing: " ++ show n) $ listToMaybe + [ map (pAccess p) (fromTo i j p) + | p <- map (flip pSetOffset 0) polys + , let pts = polygonPoints p + , Just i <- [V.elemIndex a pts] + , Just j <- [V.elemIndex b pts] + , if n == 2 then trace ("Found: " ++ show (i,j)) True else True + ] + worker n | n == pSize origin = [] + worker n = + let a = pAccess origin n + b = pAccess origin (n+1) + in findPoints n a b ++ worker (n+1) diff --git a/src/Reanimate/Morph/LeastWork.hs b/reanimate-morph/LeastWork.hs similarity index 100% rename from src/Reanimate/Morph/LeastWork.hs rename to reanimate-morph/LeastWork.hs diff --git a/src/Reanimate/Morph/LineBend.hs b/reanimate-morph/LineBend.hs similarity index 100% rename from src/Reanimate/Morph/LineBend.hs rename to reanimate-morph/LineBend.hs diff --git a/src/Reanimate/Math/Render.hs b/reanimate-morph/Render.hs similarity index 97% rename from src/Reanimate/Math/Render.hs rename to reanimate-morph/Render.hs index 1496daa..88ff3b7 100644 --- a/src/Reanimate/Math/Render.hs +++ b/reanimate-morph/Render.hs @@ -2,21 +2,16 @@ module Reanimate.Math.Render where import Codec.Picture.Types import Control.Monad --- import Data.List import qualified Data.Text as T import qualified Data.Vector as V import Linear.V2 import Reanimate.Animation import Reanimate.ColorMap --- import Reanimate.Constants import Reanimate.LaTeX --- import Reanimate.Math.Compatible --- import Reanimate.Math.EarClip -import Reanimate.Math.Polygon import Reanimate.Math.Common +import Reanimate.Math.Polygon import Reanimate.Math.SSSP import Reanimate.Math.Triangulate --- import Reanimate.Math.Visibility import Reanimate.Scene import Reanimate.Svg diff --git a/examples/morphology_point_trajectory.hs b/reanimate-morph/morphology_point_trajectory.hs similarity index 100% rename from examples/morphology_point_trajectory.hs rename to reanimate-morph/morphology_point_trajectory.hs diff --git a/reanimate.cabal b/reanimate.cabal index 38e35d8..6de9d5f 100644 --- a/reanimate.cabal +++ b/reanimate.cabal @@ -37,18 +37,6 @@ Source-Repository head Type: git Location: git://github.com/lemmih/reanimate.git --- If both 'enable-hmatrix' and 'disable-hmatrix' are false, then use --- hmatrix on Macos and Linux but not Windows. -Flag enable-hmatrix - Description: Enable hmatrix dependency on Windows (requires openblas) - Default: False - Manual: True - -Flag disable-hmatrix - Description: Force disable hmatrix dependency - Default: False - Manual: True - Flag test Description: Enable testing (requires stack) Default: False @@ -68,7 +56,6 @@ library Reanimate.Svg.Constructors Reanimate.Svg.BoundingBox Reanimate.Svg.LineCommand - -- Reanimate.Diagrams Reanimate.Transform Reanimate.Misc Reanimate.Math.Common @@ -78,13 +65,9 @@ library Reanimate.Math.Visibility Reanimate.Math.SSSP Reanimate.Math.Balloon - Reanimate.Math.Compatible - Reanimate.Math.Render Reanimate.Morph.Common Reanimate.Morph.Linear Reanimate.Morph.Rotational - Reanimate.Morph.LeastWork - Reanimate.Morph.LineBend Reanimate.Morph.Cache Reanimate.Raster Reanimate.ColorComponents @@ -102,7 +85,6 @@ library Reanimate.Builtin.Slide Reanimate.Constants Reanimate.Parameters - Reanimate.Chiphunk Reanimate.PolyShape Reanimate.GeoProjection Reanimate.Builtin.Documentation @@ -118,23 +100,48 @@ library Reanimate.Driver.Server Reanimate.Driver.Compile Paths_reanimate - build-depends: base >=4.10 && <5, - time, text, filepath, process, directory, - containers, reanimate-svg >= 0.10.0.0, xml, bytestring, lens, linear, mtl, matrix, - JuicyPixels >= 3.3.3, attoparsec, parallel, - cubicbezier, websockets >= 0.12.7.0, - hashable, fsnotify, open-browser, random, random-shuffle, base64-bytestring, - vector >= 0.12.0.0, colour, cassava, ansi-terminal, temporary, - optparse-applicative, chiphunk >= 0.1.2.1, - geojson, aeson >= 1.3.0.0, earcut >= 0.1.0.4, - split, array, cereal, vector-space, - neat-interpolation >= 0.3, ghcid, filelock - if (flag(enable-hmatrix) || !os(windows)) && !flag(disable-hmatrix) { - build-depends: hmatrix >= 0.20.0.0 - exposed-modules: Reanimate.Morph.Rigid - Reanimate.Math.Smooth - Reanimate.Math.DCEL - } + build-depends: + base >=4.10 && <5, + JuicyPixels >=3.3.3, + aeson >=1.3.0.0, + ansi-terminal >=0.8.0.4, + array >=0.5.2.0, + attoparsec >=0.13.2.0, + base64-bytestring >=1.0.0.1, + bytestring >=0.10.8.0, + cassava >=0.5.1.0, + cereal >=0.5.5.0, + colour >=2.3.4, + containers >=0.5.11.0, + cubicbezier >=0.6.0.5, + directory >=1.3.1.0, + earcut >=0.1.0.4, + filelock >=0.1.1.2, + filepath >=1.4.2, + fsnotify >=0.3.0.1, + geojson >=3.0.4, + ghcid >=0.7, + hashable >=1.2.7.0, + lens >=4.16.1, + linear >=1.20.8, + matrix >=0.3.6.1, + mtl >=2.2.2, + neat-interpolation >=0.3, + open-browser >=0.2.1.0, + optparse-applicative >=0.14.2.0, + parallel >=3.2.1.0, + process >=1.6.3.0, + random >=1.1, + random-shuffle >=0.0.4, + reanimate-svg >=0.10.0.0, + split >=0.2.3.3, + temporary >=1.3, + text >=1.2.3.0, + time >=1.8.0.0, + vector >=0.12.0.0, + vector-space >=0.13, + websockets >=0.12.7.0, + xml >=1.3.14 ghc-options: -Wall -fno-ignore-asserts test-suite spec diff --git a/src/Reanimate/Chiphunk.hs b/src/Reanimate/Chiphunk.hs deleted file mode 100644 index bccec97..0000000 --- a/src/Reanimate/Chiphunk.hs +++ /dev/null @@ -1,97 +0,0 @@ -module Reanimate.Chiphunk - ( simulate - , BodyStore - , newBodyStore - , addToBodyStore - , spaceFreeRecursive - , polyShapesToBody - , polygonsToBody - ) where - -import Chiphunk.Low -import Control.Monad -import Data.IORef -import Data.Map (Map) -import qualified Data.Map as Map -import qualified Data.Vector as V -import qualified Data.Vector.Mutable as V -import Foreign.Ptr -import Graphics.SvgTree (Tree) -import Linear.V2 (V2(..)) -import Reanimate.Animation -import Reanimate.PolyShape -import Reanimate.Svg.Constructors - -type BodyStore = IORef (Map WordPtr Tree) - -newBodyStore :: IO BodyStore -newBodyStore = newIORef Map.empty - -addToBodyStore :: BodyStore -> Body -> Tree -> IO () -addToBodyStore store body svg = do - key <- atomicModifyIORef' store $ \m -> - case Map.maxViewWithKey m of - Nothing -> (Map.singleton 1 svg, 1) - Just ((maxKey,_),_) -> - (Map.insert (maxKey+1) svg m, maxKey+1) - bodyUserData body $= wordPtrToPtr key - -renderBodyStore :: Space -> BodyStore -> IO Tree -renderBodyStore space store = do - m <- readIORef store - lst <- newIORef [] - spaceEachBody space (\body _dat -> do - key <- get (bodyUserData body) - case Map.lookup (ptrToWordPtr key) m of - Nothing -> putStrLn "Body doesn't have an associated SVG" - Just svg -> do - Vect posX posY <- get $ bodyPosition body - angle <- get $ bodyAngle body - let bodySvg = - translate posX posY $ - rotate (angle/pi*180) - svg - modifyIORef lst (bodySvg:) - ) nullPtr - result <- readIORef lst - return $ mkGroup result - - -simulate :: Space -> BodyStore -> Double -> Int -> Double -> IO Animation -simulate space store fps stepsPerFrame dur = do - let timeStep = 1/(fps*fromIntegral stepsPerFrame) - frames = round (dur * fps) - v <- V.new frames - forM_ [0..frames-1] $ \nth -> do - svg <- renderBodyStore space store - V.write v nth svg - replicateM_ stepsPerFrame $ spaceStep space timeStep - frozen <- V.unsafeFreeze v - return $ mkAnimation dur $ \t -> - let key = round (t * fromIntegral (frames-1)) - in frozen V.! key - -polyShapesToBody :: Space -> [PolyShape] -> IO Body -polyShapesToBody space poly = - polygonsToBody space (map (map toVect) $ plDecompose poly) - where - toVect (V2 x y) = Vect x y - -polygonsToBody :: Space -> [[Vect]] -> IO Body -polygonsToBody space polygons = do - plBody <- bodyNew 0 0 - spaceAddBody space plBody - - forM_ polygons $ \vects -> do - polyShape <- polyShapeNewRaw plBody vects 0.00 - shapeDensity polyShape $= 1 - spaceAddShape space polyShape - shapeFriction polyShape $= 0.7 - shapeElasticity polyShape $= 0.5 - return plBody - -spaceFreeRecursive :: Space -> IO () -spaceFreeRecursive space = do - spaceEachBody space (\body _ -> bodyFree body) nullPtr - spaceEachShape space (\shape _ -> shapeFree shape) nullPtr - spaceFree space diff --git a/src/Reanimate/Math/Common.hs b/src/Reanimate/Math/Common.hs index 0350bcd..b945c69 100644 --- a/src/Reanimate/Math/Common.hs +++ b/src/Reanimate/Math/Common.hs @@ -13,12 +13,13 @@ module Reanimate.Math.Common , area2X -- :: Fractional a => V2 a -> V2 a -> V2 a -> a , epsilon -- :: Fractional a => a , epsEq -- :: (Ord a, Fractional a) => a -> a -> Bool - , isLeftTurn -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool - , isLeftTurnOrLinear -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool - , isRightTurn -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool - , isRightTurnOrLinear -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool - , direction -- :: Fractional a => V2 a -> V2 a -> V2 a -> a + , isLeftTurn -- :: (Num a, Ord a) => V2 a -> V2 a -> V2 a -> Bool + , isLeftTurnOrLinear -- :: (Num a, Ord a) => V2 a -> V2 a -> V2 a -> Bool + , isRightTurn -- :: (Num a, Ord a) => V2 a -> V2 a -> V2 a -> Bool + , isRightTurnOrLinear -- :: (Num a, Ord a) => V2 a -> V2 a -> V2 a -> Bool + , direction -- :: Num a => V2 a -> V2 a -> V2 a -> a , isInside -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> V2 a -> Bool + , isInsideStrict -- :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> V2 a -> Bool , barycentricCoords -- :: Fractional a => V2 a -> V2 a -> V2 a -> V2 a -> (a, a, a) , rayIntersect -- :: (Fractional a, Ord a) => (V2 a,V2 a) -> (V2 a,V2 a) -> Maybe (V2 a) , isBetween -- :: (Ord a, Fractional a) => V2 a -> (V2 a, V2 a) -> Bool @@ -77,21 +78,26 @@ epsilon = 1e-13 epsEq :: (Ord a, Fractional a) => a -> a -> Bool epsEq a b = abs (a-b) < epsilon +compareEpsZero :: (Ord a, Fractional a) => a -> Ordering +compareEpsZero val + | abs val < epsilon = EQ + | otherwise = compare val 0 + {-# INLINE isLeftTurn #-} -- Left turn. isLeftTurn :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool isLeftTurn p1 p2 p3 = - case compare (direction p1 p2 p3) 0 of + case compareEpsZero (direction p1 p2 p3) of LT -> True - EQ -> False -- colnear + EQ -> False -- colinear GT -> False {-# INLINE isLeftTurnOrLinear #-} isLeftTurnOrLinear :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool isLeftTurnOrLinear p1 p2 p3 = - case compare (direction p1 p2 p3) 0 of + case compareEpsZero (direction p1 p2 p3) of LT -> True - EQ -> True -- colnear + EQ -> True -- colinear GT -> False {-# INLINE isRightTurn #-} @@ -103,15 +109,22 @@ isRightTurnOrLinear :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool isRightTurnOrLinear a b c = not (isLeftTurn a b c) {-# INLINE direction #-} -direction :: Fractional a => V2 a -> V2 a -> V2 a -> a +direction :: Num a => V2 a -> V2 a -> V2 a -> a direction p1 p2 p3 = crossZ (p3-p1) (p2-p1) {-# INLINE isInside #-} isInside :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> V2 a -> Bool isInside a b c d = - s >= 0 && s <= 1 && t >= 0 && t <= 1 + s >= 0 && s <= 1 && t >= 0 && t <= 1 && i >= 0 && i <= 1 where - (s, t, _) = barycentricCoords a b c d + (s, t, i) = barycentricCoords a b c d + +{-# INLINE isInsideStrict #-} +isInsideStrict :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> V2 a -> Bool +isInsideStrict a b c d = + s > 0 && s < 1 && t > 0 && t < 1 && i > 0 && i < 1 + where + (s, t, i) = barycentricCoords a b c d {-# INLINE barycentricCoords #-} barycentricCoords :: Fractional a => V2 a -> V2 a -> V2 a -> V2 a -> (a, a, a) @@ -153,7 +166,7 @@ lineIntersect a b = -- circleIntersect :: (Ord a, Fractional a) => (V2 a, V2 a) -> (V2 a, V2 a) -> [V2 a] -distSquared :: (Fractional a) => V2 a -> V2 a -> a +distSquared :: (Num a) => V2 a -> V2 a -> a distSquared a b = quadrance (a ^-^ b) approxDist :: (Real a, Fractional a) => V2 a -> V2 a -> a diff --git a/src/Reanimate/Math/DCEL.hs b/src/Reanimate/Math/DCEL.hs index bf0b703..cf7a187 100644 --- a/src/Reanimate/Math/DCEL.hs +++ b/src/Reanimate/Math/DCEL.hs @@ -8,22 +8,21 @@ import Control.Lens import Control.Monad.State import Control.Monad.Writer import Data.List -import qualified Data.List as L -import Data.Map (Map) -import qualified Data.Map as M +import qualified Data.List as L +import Data.Map (Map) +import qualified Data.Map as M import Data.Maybe -import qualified Data.Set as S -import qualified Data.Text as T -import qualified Data.Vector as V --- import Debug.Trace +import qualified Data.Set as S +import qualified Data.Text as T +import qualified Data.Vector as V import Linear.Metric import Linear.V2 -import Linear.Vector (lerp, (^/)) +import Linear.Vector (lerp, (^/)) import Reanimate -import Reanimate.Math.Common (isInside, lineIntersect, triangleAngles) -import Reanimate.Math.Polygon -import qualified Reanimate.Morph.Rigid as Rigid +import Reanimate.Math.Common (isInsideStrict, lineIntersect, + triangleAngles) import Text.Printf +-- import Debug.Trace type VertexId = Int type EdgeId = Int @@ -179,8 +178,8 @@ emptyMesh = Mesh 1 0 M.empty M.empty M.empty polygonsMesh :: Eq a => [a] -> [[a]] -> MeshM a () polygonsMesh outer trigs = do - polygonMeshOuter (reverse outer) modify $ meshFaces .~ M.empty + polygonMeshOuter (reverse outer) forM_ trigs $ \pts -> do innerFace <- createFace vIds <- mapM getVertex pts @@ -552,26 +551,43 @@ meshSmoothPosition :: Mesh (V2 Double) -> Mesh (V2 Double) meshSmoothPosition = execState worker where worker = gets steinerNodes >>= mapM_ smoothVertex + -- worker = mapM_ smoothVertex [221] smoothVertex :: VertexId -> MeshM (V2 Double) () smoothVertex steiner = do self <- _vertexPosition <$> requireVertex steiner es <- gets (outgoingEdges steiner) vs <- mapM (requireVertex . _edgeVertex) es - let ps = V.fromList (sortEdges self (map _vertexPosition vs)) + let ps = V.fromList (map _vertexPosition vs') + vs' = sortVertices self vs angleBased = angleSmooth self ps laplacian = sum ps ^/ (fromIntegral $ length ps) -- laplacian + -- trace ("self: " ++ show self) $ return () + -- trace ("angleBased: " ++ show angleBased) $ return () + -- trace ("Edges: " ++ show (map _vertexId vs')) $ return () + -- trace ("Valid angleBased: " ++ show (isValidLocation self ps angleBased)) $ return () + -- trace ("Valid laplacian: " ++ show (isValidLocation self ps laplacian)) $ return () + -- let checks = [ isInsideStrict self a b angleBased + -- | i <- [0 .. length ps - 1] + -- , let a = ps V.! i + -- b = ps V.! mod (i + 1) (length ps) + -- ] + -- trace ("Checks: " ++ show checks) $ return () if isValidLocation self ps angleBased then modifyVertex (vertexPosition .~ angleBased) steiner else if isValidLocation self ps laplacian then modifyVertex (vertexPosition .~ laplacian) steiner else return () where - sortEdges :: V2 Double -> [V2 Double] -> [V2 Double] - sortEdges = sortOn . dir + -- sortEdges :: V2 Double -> [V2 Double] -> [V2 Double] + -- sortEdges = sortOn . dir + sortVertices :: V2 Double -> [Vertex (V2 Double)] -> [Vertex (V2 Double)] + sortVertices self = sortOn (dirV self) -- Direction from south of 'a', to 'a', to 'b'. dir :: V2 Double -> V2 Double -> Double dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a))) + dirV :: V2 Double -> Vertex (V2 Double) -> Double + dirV a v = dir a (_vertexPosition v) angleSmooth :: V2 Double -> V.Vector (V2 Double) -> V2 Double @@ -604,7 +620,7 @@ angleSmooth origin js = V.sum (V.generate n nth) ^/ V.sum (V.generate n factor) isValidLocation :: V2 Double -> V.Vector (V2 Double) -> V2 Double -> Bool isValidLocation origin edges newLoc = or - [ isInside origin a b newLoc + [ isInsideStrict origin a b newLoc | i <- [0 .. length edges - 1] , let a = edges V.! i b = edges V.! mod (i + 1) (length edges) @@ -617,6 +633,11 @@ isValidLocation origin edges newLoc = dir :: V2 Double -> V2 Double -> Double dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a))) +isCCW :: (Ord a, Num a) => V2 a -> V2 a -> V2 a -> Bool +isCCW a b c = sum [fn a b, fn b c, fn c a] < 0 + where + fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1) + minAngle :: V2 Double -> V.Vector (V2 Double) -> Double minAngle origin edges = minimum $ concat [ [a1, a2, a3] @@ -703,7 +724,7 @@ delaunayFlip' eid m = ((f1',f2'), mAfter) = runState (flipEdge edge) m afterAng = min (faceMinAngle f1' mAfter) (faceMinAngle f2' mAfter) in - if (afterAng < beforeAng) || f1' == f1 + if (afterAng <= beforeAng) || f1' == f1 then Nothing else Just mAfter @@ -861,31 +882,3 @@ renderMeshStats mesh = mkGroup meanAngle = L.sort angs !! (length angs `div` 2) avgAngle = sum angs / (fromIntegral $ length angs) maxAngle = maximum angs - -toRigidMesh :: Mesh (V2 Double) -> Mesh (V2 Double) -> Rigid.Mesh -toRigidMesh meshA meshB = Rigid.Mesh - { meshPointsA = pointsA - , meshPointsB = pointsB - , meshOutline = outline - , meshTriangles = trigs } - where - pointsA = V.fromList $ map _vertexPosition $ M.elems (meshA^.meshVertices) - pointsB = V.fromList $ map _vertexPosition $ M.elems (meshB^.meshVertices) - outline = V.fromList - [ fromJust (V.elemIndex v pointsA) - | eid <- faceEdges (meshA^.meshOuterFace) meshA - , let edge = meshGetEdge eid meshA - v = _vertexPosition (meshGetVertex (edge^.edgeVertex) meshA) - ] - trigs = V.fromList - [ (aIdx, bIdx, cIdx) - | fid <- M.keys (meshA^.meshFaces) - , fid /= (meshA^.meshOuterFace) - , let edges = map (`meshGetEdge` meshA) (faceEdges fid meshA) - vs = map (`meshGetVertex` meshA) $ map _edgeVertex edges - ps = map _vertexPosition vs - [aIdx, bIdx, cIdx] = map (fromJust . (`V.elemIndex` pointsA)) ps - p = mkPolygon (V.fromList $ map (fmap realToFrac) ps) - , pIsSimple p || error "invalid polygon" - ] - diff --git a/src/Reanimate/Math/Polygon.hs b/src/Reanimate/Math/Polygon.hs index 8104c47..a0c35cd 100644 --- a/src/Reanimate/Math/Polygon.hs +++ b/src/Reanimate/Math/Polygon.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE BangPatterns #-} module Reanimate.Math.Polygon ( APolygon(..) , Polygon @@ -5,6 +6,7 @@ module Reanimate.Math.Polygon , P , mkPolygon -- :: (Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a , mkPolygonFromRing -- :: (Fractional a, Ord a) => Ring a -> APolygon a + , castPolygon -- :: (Real a, Fractional b, Ord a) => APolygon a -> APolygon b , pParent -- :: Polygon -> Int -> Int -> Int , pSetOffset -- :: APolygon a -> Int -> APolygon a , pAdjustOffset -- :: APolygon a -> Int -> APolygon a @@ -33,7 +35,10 @@ module Reanimate.Math.Polygon , pCircumference -- :: (Real a, Fractional a) => APolygon a -> a , pCircumference' -- :: (Real a, Fractional a) => APolygon a -> Double , pAddPoints -- :: Int -> Polygon -> Polygon + , pAddPointsRestricted -- :: [Int] -> Int -> Polygon -> Polygon + , pAddPointsBetween -- :: (Fractional a, Ord a, Real a) => (Int, Int) -> Int -> APolygon a -> APolygon a , pRayIntersect -- :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational) + , pOverlap -- :: Polygon -> Polygon -> Polygon , pCuts -- :: Polygon -> [(Polygon,Polygon)] , pCutEqual -- :: Polygon -> (Polygon, Polygon) -- * Triangulation @@ -119,14 +124,14 @@ instance Show a => Show (APolygon a) where instance Hashable a => Hashable (APolygon a) where hashWithSalt s p = V.foldl' hashWithSalt s (polygonPoints p) -instance (Real a, Fractional a, Ord a, Serialize a) => Serialize (APolygon a) where +instance (Real a, Fractional a, Serialize a) => Serialize (APolygon a) where put = put . V.toList . polygonPoints get = mkPolygon . V.fromList <$> get pRing :: APolygon a -> Ring a pRing = ringPack . polygonPoints -mkPolygon :: (Real a, Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a +mkPolygon :: (Real a, Fractional a) => V.Vector (V2 a) -> APolygon a mkPolygon points = Polygon { polygonPoints = points , polygonOffset = 0 @@ -139,6 +144,9 @@ mkPolygon points = Polygon trig = earCut ring -- earClip ring +castPolygon :: (Real a, Real b, Fractional b) => APolygon a -> APolygon b +castPolygon = mkPolygon . V.map (fmap realToFrac) . polygonPoints + mkPolygonFromRing :: (Real a, Fractional a, Ord a) => Ring a -> APolygon a mkPolygonFromRing = mkPolygon . ringUnpack @@ -146,7 +154,7 @@ pUnsafeMap :: (Ring a -> Ring a) -> APolygon a -> APolygon a pUnsafeMap fn p = p{ polygonPoints = ringUnpack (fn (pRing p)) } -- pParent p i j = shortest-path parent from j to i -pParent :: Polygon -> Int -> Int -> Int +pParent :: APolygon a -> Int -> Int -> Int pParent p i j = (sTree V.! mod (j + polygonOffset p) n - polygonOffset p) `mod` n where @@ -164,6 +172,7 @@ pAdjustOffset :: APolygon a -> Int -> APolygon a pAdjustOffset p offset = p { polygonOffset = (polygonOffset p + offset) `mod` pSize p } +{-# INLINE pSize #-} pSize :: APolygon a -> Int pSize = length . polygonPoints @@ -316,9 +325,10 @@ pIsInside p point = or -- findSmallest = minimumBy (comparing area2X) -- shareEdge p1 p2 = - +{-# INLINE pAccess #-} pAccess :: APolygon a -> Int -> V2 a -pAccess p i = polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p) +pAccess p i = -- polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p) + polygonPoints p `V.unsafeIndex` ((polygonOffset p + i) `mod` pSize p) triangle :: Polygon triangle = mkPolygon $ V.fromList [V2 1 1, V2 0 0, V2 2 0] @@ -503,7 +513,7 @@ pDeoverlap p = mkPolygon arr in lerp 0.99999 this prev else pAccess p n -pCycles :: Polygon -> [Polygon] +pCycles :: APolygon a -> [APolygon a] pCycles p = map (pAdjustOffset p) [0 .. pSize p-1] pCycle :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a @@ -526,7 +536,7 @@ pCycle p t = mkPolygon $ worker 0 0 len = pCircumference' p limit = t * len -pCentroid :: Polygon -> V2 Rational +pCentroid :: Fractional a => APolygon a -> V2 a pCentroid p = V2 cx cy where a = pArea p @@ -535,13 +545,29 @@ pCentroid p = V2 cx cy fnX (V2 x y) (V2 x' y') = (x+x')*(x*y' - x'*y) fnY (V2 x y) (V2 x' y') = (y+y')*(x*y' - x'*y) -pMapEdges :: (V2 Rational -> V2 Rational -> a) -> Polygon -> V.Vector a -pMapEdges fn p = V.generate (pSize p) $ \i -> - fn (pAccess p i) (pAccess p $ i+1) +{-# INLINE pMapEdges #-} +pMapEdges :: (V2 a -> V2 a -> b) -> APolygon a -> V.Vector b +pMapEdges fn p = V.generate n $ \i -> + if i == n-1 + then fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) + else fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1)) + where + n = pSize p + arr = polygonPoints p -pArea :: Polygon -> Rational +{-# SPECIALIZE pArea :: APolygon Double -> Double #-} +{-# SPECIALIZE pArea :: APolygon Rational -> Rational #-} +pArea :: (Fractional a) => APolygon a -> a pArea p = - 0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p) + -- 0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p) + 0.5 * worker 0 0 + where + fn (V2 x y) (V2 x' y') = x*y' - x'*y + arr = polygonPoints p + worker !acc i + | i == pSize p - 1 = acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` 0) + | otherwise = + worker (acc + fn (arr `V.unsafeIndex` i) (arr `V.unsafeIndex` (i+1))) (i+1) pCircumference :: (Real a, Fractional a) => APolygon a -> a pCircumference p = sum @@ -555,7 +581,7 @@ pCircumference' p = sum -- Add points by splitting the longest lines in half repeatedly. -pAddPoints :: Int -> Polygon -> Polygon +pAddPoints :: (Ord a, Real a, Fractional a) => Int -> APolygon a -> APolygon a pAddPoints n p | n <= 0 = p pAddPoints n p = pAddPoints (n-1) $ mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] @@ -572,6 +598,45 @@ pAddPoints n p = pAddPoints (n-1) $ distSquared (pAccess p a) (pAccess p $ a+1) `compare` distSquared (pAccess p b) (pAccess p $ b+1) +pAddPointsRestricted :: (Fractional a, Ord a, Real a) => [(V2 a, V2 a)] -> Int -> APolygon a -> APolygon a +pAddPointsRestricted _immutableEdges n p | n <= 0 = p +pAddPointsRestricted immutableEdges n p = pAddPointsRestricted immutableEdges (n-1) $ + mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] + where + isImmutable idx = + (pAccess p idx, pAccess p $ idx+1) `elem` immutableEdges || + (pAccess p $ idx+1, pAccess p idx) `elem` immutableEdges + worker idx + | idx == longestEdge && not (isImmutable idx) = + let start = pAccess p idx + end = pAccess p $ idx+1 + middle = lerp 0.5 end start + in [start, middle] + | otherwise = [pAccess p idx] + longestEdge = maximumBy cmpLength [0 .. pSize p-1] + cmpLength a _ | isImmutable a = LT + cmpLength _ b | isImmutable b = GT + cmpLength a b = + distSquared (pAccess p a) (pAccess p $ a+1) `compare` + distSquared (pAccess p b) (pAccess p $ b+1) + +pAddPointsBetween :: (Fractional a, Ord a, Real a) => (Int, Int) -> Int -> APolygon a -> APolygon a +pAddPointsBetween _ n p | n <= 0 = p +pAddPointsBetween (i,l) n p = pAddPointsBetween (i,l+1) (n-1) $ + mkPolygon $ V.fromList $ concatMap worker [0 .. pSize p-1] + where + worker idx + | idx == longestEdge = + let start = pAccess p idx + end = pAccess p $ idx+1 + middle = lerp 0.5 end start + in [start, middle] + | otherwise = [pAccess p idx] + longestEdge = maximumBy cmpLength [i .. i+l-1] + cmpLength a b = + distSquared (pAccess p a) (pAccess p $ a+1) `compare` + distSquared (pAccess p b) (pAccess p $ b+1) + -- addPoints :: Int -> Polygon -> Polygon -- addPoints n p = mkPolygon $ V.fromList $ worker n 0 (map (pAccess p) [0..s]) -- where @@ -603,11 +668,11 @@ pIsCCW p = V.sum (pMapEdges fn p) < 0 fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1) {-# INLINE pRayIntersect #-} -pRayIntersect :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational) +pRayIntersect :: (Fractional a, Ord a) => APolygon a -> (Int, Int) -> (Int,Int) -> Maybe (V2 a) pRayIntersect p (a,b) (c,d) = rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d) -pCuts :: Polygon -> [(Polygon,Polygon)] +pCuts :: (Real a, Fractional a) => APolygon a -> [(APolygon a,APolygon a)] pCuts p = [ pCutAt (pAdjustOffset p i) (j-i) | i <- [0 .. pSize p-1 ] @@ -615,26 +680,41 @@ pCuts p = , (j+1) `mod` pSize p /= i , pParent p i j == i ] -pCutEqual :: Polygon -> (Polygon, Polygon) +pCutEqual :: (Real a, Fractional a) => APolygon a -> (APolygon a, APolygon a) pCutEqual p = fromMaybe (p,p) $ listToMaybe $ sortOn f $ pCuts p where f (a,b) = abs (pArea a - pArea b) -- FIXME: This should be more efficient -pCutAt :: Polygon -> Int -> (Polygon, Polygon) +pCutAt :: (Real a, Fractional a) => APolygon a -> Int -> (APolygon a, APolygon a) pCutAt p i = (mkPolygon $ V.fromList left, mkPolygon $ V.fromList right) where n = pSize p left = map (pAccess p) [0 .. i] right = map (pAccess p) (0:[i..n-1]) - +pOverlap :: (Fractional a, Ord a, Real a) => APolygon a -> APolygon a -> APolygon a +pOverlap a b = mkPolygon $ V.fromList $ clearDups $ concatMap edgeIntersect [0 .. pSize a-1] + where + clearDups (x:y:xs) + | x == y = clearDups (y:xs) + | otherwise = x : clearDups (y:xs) + clearDups xs = xs + edgeIntersect edge = + sortOn (distSquared (pAccess a edge)) $ catMaybes + [ lineIntersect (aP, aP') (bP, bP') + | i <- [0 .. pSize b-1] + , let aP = pAccess a edge + aP' = pAccess a (edge+1) + bP = pAccess b i + bP' = pAccess b (i+1) + ] --------------------------------------------------------- -- SSSP visibility and SSSP windows -ssspVisibility :: Polygon -> Polygon +ssspVisibility :: (Real a, Fractional a) => APolygon a -> APolygon a ssspVisibility p = mkPolygon $ V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1]) where diff --git a/src/Reanimate/Morph/Common.hs b/src/Reanimate/Morph/Common.hs index 7732703..f698802 100644 --- a/src/Reanimate/Morph/Common.hs +++ b/src/Reanimate/Morph/Common.hs @@ -26,7 +26,7 @@ import Linear.V2 import Reanimate.Animation import Reanimate.ColorComponents import Reanimate.Ease -import Reanimate.Math.Polygon (Polygon, mkPolygon, pAddPoints, +import Reanimate.Math.Polygon (Polygon, APolygon, mkPolygon, pAddPoints, pCentroid, pCutEqual, pSize, polygonPoints) import Reanimate.PolyShape @@ -55,6 +55,7 @@ data Morph = Morph , morphObjectCorrespondence :: ObjectCorrespondence } +{-# INLINE morph #-} morph :: Morph -> SVG -> SVG -> Double -> SVG morph Morph{..} src dst = \t -> case t of @@ -77,7 +78,7 @@ morph Morph{..} src dst = \t -> , let arranged = morphPointCorrespondence srcPoly' dstPoly' ] -normalizePolygons :: Polygon -> Polygon -> (Polygon, Polygon) +normalizePolygons :: (Real a, Fractional a) => APolygon a -> APolygon a -> (APolygon a, APolygon a) normalizePolygons src dst = (pAddPoints (max 0 $ dstN-srcN) src ,pAddPoints (max 0 $ srcN-dstN) dst) diff --git a/src/Reanimate/Morph/Linear.hs b/src/Reanimate/Morph/Linear.hs index e6fbda0..266c350 100644 --- a/src/Reanimate/Morph/Linear.hs +++ b/src/Reanimate/Morph/Linear.hs @@ -2,6 +2,7 @@ module Reanimate.Morph.Linear ( linear, rawLinear , linearCorrespondence , closestLinearCorrespondence + , closestLinearCorrespondenceA , linearTrajectory ) where @@ -32,7 +33,10 @@ linearCorrespondence :: PointCorrespondence linearCorrespondence = normalizePolygons closestLinearCorrespondence :: PointCorrespondence -closestLinearCorrespondence src' dst' = +closestLinearCorrespondence = closestLinearCorrespondenceA + +closestLinearCorrespondenceA :: (Real a, Fractional a) => APolygon a -> APolygon a -> (APolygon a, APolygon a) +closestLinearCorrespondenceA src' dst' = (src, worker dst (score dst) options) where (src, dst) = normalizePolygons src' dst' @@ -44,7 +48,8 @@ closestLinearCorrespondence src' dst' = else worker bestP bestPScore xs options = pCycles dst score p = sum - [ approxDist (pAccess src n) (pAccess p n) + [ -- approxDist (pAccess src n) (pAccess p n) + distSquared (pAccess src n) (pAccess p n) | n <- [0 .. pSize src-1] ] linearTrajectory :: Trajectory diff --git a/src/Reanimate/Morph/Rigid.hs b/src/Reanimate/Morph/Rigid.hs index 1bc9290..ae4186c 100644 --- a/src/Reanimate/Morph/Rigid.hs +++ b/src/Reanimate/Morph/Rigid.hs @@ -66,28 +66,35 @@ on the shear. -} module Reanimate.Morph.Rigid where -import Data.Foldable (toList) +import Control.Lens +import Data.Foldable (toList) +import qualified Data.Map as M +import Data.Maybe import Data.Vector (Vector) import qualified Data.Vector as V import Linear.Quaternion -import Linear.Vector import Linear.V2 import Linear.V3 +import Linear.Vector import qualified Numeric.LinearAlgebra as Matrix -import Numeric.LinearAlgebra.HMatrix (GMatrix, Matrix, - toLists, (!), (><)) -import Reanimate.Animation -import Reanimate.Svg - +import Numeric.LinearAlgebra.HMatrix (GMatrix, Matrix, toLists, (!), + (><)) +import Reanimate.Animation +import Reanimate.Math.Compatible (compatiblyTriangulateP) +import qualified Reanimate.Math.DCEL as DCEL +import Reanimate.Math.Polygon (mkPolygon, pCopy, + polygonPoints) +import Reanimate.Morph.Common +import Reanimate.Svg type P = V2 Double type Trig = (P,P,P) type RelTrig = (Int,Int,Int) data Mesh = Mesh - { meshPointsA :: Vector P - , meshPointsB :: Vector P - , meshOutline :: Vector Int + { meshPointsA :: Vector P + , meshPointsB :: Vector P + , meshOutline :: Vector Int -- , meshSteiner :: Vector Int , meshTriangles :: Vector RelTrig } @@ -148,7 +155,7 @@ applyA :: Matrix Double -> Trig -> Trig applyA a p = case toLists (a <> matP) of [ [x1, x2], [y1, y2] ] -> (V2 0 0, V2 x1 y1, V2 x2 y2) - _ -> error "invalid matrix" + _ -> error "invalid matrix" where matP = trigToMatrix p @@ -178,12 +185,12 @@ trigToMatrix (p1,p2,p3) = matP data Prep = Prep - { prepPivot :: (P, P) - , prepPointsA :: Vector P - , prepPointsB :: Vector P - , prepRS :: Vector (Matrix Double, Matrix Double) - , prepRSRev :: Vector (Matrix Double, Matrix Double) - , prepUToB :: GMatrix + { prepPivot :: (P, P) + , prepPointsA :: Vector P + , prepPointsB :: Vector P + , prepRS :: Vector (Matrix Double, Matrix Double) + , prepRSRev :: Vector (Matrix Double, Matrix Double) + , prepUToB :: GMatrix } symmetric :: Bool @@ -204,7 +211,7 @@ prepare Mesh{..} = Prep where aOrigin = meshPointsA V.! pivotIdx bOrigin = meshPointsB V.! pivotIdx - pivotIdx = 0 + pivotIdx = 0 -- pivotIdx = case V.head meshTriangles of -- (a,_,_) -> a mkAbs p (a,b,c) = (p V.! a,p V.! b,p V.! c) @@ -252,7 +259,7 @@ interpolate Prep{..} t = V.fromList $ False 1e-9 1e-9 - 1000 + 10000 prepUToB b (Matrix.fromList $ concat [ [x,y] | V2 x y <- V.toList target ]) @@ -261,7 +268,7 @@ interpolate Prep{..} t = V.fromList $ then prepPointsA else prepPointsB worker (x:y:xs) = V2 x y ^+^ pivot : worker xs - worker _ = [] + worker _ = [] pivot = case prepPivot of (src, dst) -> lerp t dst src n = V.length prepRS @@ -279,3 +286,61 @@ interpolate Prep{..} t = V.fromList $ , let a = computeA_RSt r s (1-t) , j <- [0..3] ] + +toRigidMesh :: DCEL.Mesh (V2 Double) -> DCEL.Mesh (V2 Double) -> Mesh +toRigidMesh meshA meshB = Mesh + { meshPointsA = pointsA + , meshPointsB = pointsB + , meshOutline = outline + , meshTriangles = trigs } + where + pointsA = V.fromList $ map DCEL._vertexPosition $ M.elems (meshA^.DCEL.meshVertices) + pointsB = V.fromList $ map DCEL._vertexPosition $ M.elems (meshB^.DCEL.meshVertices) + outline = V.fromList + [ fromJust (V.elemIndex v pointsA) + | eid <- DCEL.faceEdges (meshA^.DCEL.meshOuterFace) meshA + , let edge = DCEL.meshGetEdge eid meshA + v = DCEL._vertexPosition (DCEL.meshGetVertex (edge^.DCEL.edgeVertex) meshA) + ] + trigs = V.fromList + [ (aIdx, bIdx, cIdx) + | fid <- M.keys (meshA^.DCEL.meshFaces) + , fid /= (meshA^.DCEL.meshOuterFace) + , let edges = map (`DCEL.meshGetEdge` meshA) (DCEL.faceEdges fid meshA) + vs = map (`DCEL.meshGetVertex` meshA) $ map DCEL._edgeVertex edges + ps = map DCEL._vertexPosition vs + [aIdx, bIdx, cIdx] = map (fromJust . (`V.elemIndex` pointsA)) ps + -- p = mkPolygon (V.fromList $ map (fmap realToFrac) ps) + -- , pIsSimple p || error "invalid polygon" + ] + + +{-# INLINE rigidMorph #-} +rigidMorph :: Trajectory +rigidMorph (p1', p2') = \t -> + let points = V.map (fmap realToFrac) $ interpolate optPrep t + in mkPolygon (V.map (\i -> points V.! i) (meshOutline optMesh)) + where + p1 = pCopy p1' + p2 = pCopy p2' + (p1s,p2s) = unzip (compatiblyTriangulateP p1 p2) + m2 = DCEL.buildMesh $ DCEL.polygonsMesh + (map (fmap realToFrac) $ V.toList $ polygonPoints p2) + (map (map (fmap realToFrac) . V.toList . polygonPoints) p2s) + + m1 = DCEL.buildMesh $ DCEL.polygonsMesh + (map (fmap realToFrac) $ V.toList $ polygonPoints p1) + (map (map (fmap realToFrac) . V.toList . polygonPoints) p1s) + + optPrep = prepare optMesh + optMesh = toRigidMesh m1final m2final + pipeline1 = last . take 20 . iterate + (uncurry DCEL.delaunayFlip . + uncurry DCEL.splitInternalEdges . + (\(a,b) -> (DCEL.meshSmoothPosition a, DCEL.meshSmoothPosition b))) + (m1final, m2final) = last $ take 20 $ iterate + (pipeline1 . + uncurry DCEL.splitLongestEdge . + pipeline1 + ) (m1,m2) + diff --git a/stack-lts-12.yaml b/stack-lts-12.yaml index 3b98243..3860cae 100644 --- a/stack-lts-12.yaml +++ b/stack-lts-12.yaml @@ -6,19 +6,12 @@ packages: - . extra-deps: -- reanimate-svg-0.10.0.0 -- chiphunk-0.1.2.1 -- cubicbezier-0.6.0.6@sha256:2191ff47144d9a13a2784651a33d340cd31be1926a6c188925143103eb3c8db3 -- fast-math-1.0.2@sha256:91181eb836e54413cc5a841e797c42b2264954e893ea530b6fc4da0dccf6a8b7 -- matrices-0.5.0@sha256:b2761813f6a61c84224559619cc60a16a858ac671c8436bbac8ec89e85473058 +- earcut-0.1.0.4@sha256:d5118b3eecf24d130263d81fb30f1ff56b1db43036582bfd1d8cc9ba3adae8be,1010 - geojson-4.0.1@sha256:276de5cb2aa3e07179a8d42c184b1a4e52d2c8d23cf2eb989ecfa6adbe726227 -- hmatrix-0.20.0.0@sha256:d79a9218e314f1a2344457c3851bd1d2536518ecb5f1a2fcd81daa45e46cd025,4870 +- JuicyPixels-3.3.5@sha256:5c67dc066d67f5045221a8202f701f4c77eab8220f44414c00fadcdc3e4f872b,3911 - websockets-0.12.7.1@sha256:4363a5ca794af094971266a1397807ac55d5b3124645060fb78c64ac4c88a1f1,7863 - clock-0.8@sha256:b4ae207e2d3761450060a0d0feb873269233898039c76fceef9cc1a544067767,4113 -- JuicyPixels-3.3.5@sha256:5c67dc066d67f5045221a8202f701f4c77eab8220f44414c00fadcdc3e4f872b,3911 -- th-printf-0.7@sha256:c6b2c08e62a32a327b29b1a3739939f0c85b4079ff2ae6d79fac639f25ee48b1,2297 -- earcut-0.1.0.4@sha256:d5118b3eecf24d130263d81fb30f1ff56b1db43036582bfd1d8cc9ba3adae8be,1010 +- reanimate-svg-0.10.0.0 - tasty-rerun-1.1.17@sha256:d4a3ccb0f63f499f36edc71b33c0f91c850eddb22dd92b928aa33b8459f3734a,1373 - tasty-1.3.1@sha256:01e35c97f7ee5ccbc28f21debea02a38cd010d53b4c3087f5677c5d06617a507,2520 -- optparse-applicative-0.15.1.0@sha256:29ff6146aabf54d46c4c8788e8d1eadaea27c94f6d360c690c5f6c93dac4b07e,4810 - ansi-terminal-0.10.3@sha256:e2fbcef5f980dc234c7ad8e2fa433b0e8109132c9e643bc40ea5608cd5697797,3226 diff --git a/videos/morph/dcel.hs b/videos/morph/dcel.hs index 10b85be..a4e3c58 100644 --- a/videos/morph/dcel.hs +++ b/videos/morph/dcel.hs @@ -52,8 +52,8 @@ p2 = pCopy p2' -- normalizePolygons closestLinearCorrespondence -- leastWork defaultStretchCosts defaultBendCosts - (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 4 $ latex "S") - (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 4 $ latex "C") + (pAtCenter $ unsafeSVGToPolygon 0.5 $ scale 6 $ latex "X") + (pAtCenter $ unsafeSVGToPolygon 0.5 $ scale 6 $ latex "S") (p1s,p2s) = unzip (compatiblyTriangulateP p1 p2) @@ -65,6 +65,10 @@ m1 = buildMesh $ polygonsMesh (map (fmap realToFrac) $ V.toList $ polygonPoints p1) (map (map (fmap realToFrac) . V.toList . polygonPoints) p1s) +doSmooth (a,b) = (meshSmoothPosition a, meshSmoothPosition b) +step0 = (m1,m2) +step1 = uncurry delaunayFlip step0 +step2 = doSmooth step1 main :: IO () main = reanimate $ sceneAnimation $ do @@ -76,7 +80,7 @@ main = reanimate $ sceneAnimation $ do -- DCEL.renderMesh m s <- newVar 1 mVar <- newVar (m1, m2) - let V2 centerX centerY = V2 0 0 -- realToFrac <$> pAccess p2 33 + let V2 centerX centerY = V2 0 0 -- V2 (-0.2) (-2)-- realToFrac <$> pAccess p2 33 adjustZ 2 $ newSprite_ $ do ~(m1, m2) <- unVar mVar pure $ mkGroup @@ -91,27 +95,33 @@ main = reanimate $ sceneAnimation $ do [ mkGroup [] , DCEL.renderMeshColored m1 , DCEL.renderMesh (0.05/sc) m1 - -- , renderMeshEdges m - -- , DCEL.renderMeshSimple (0.10/sc) m + -- , renderMeshEdges m1 + , withFillColor "grey" $ DCEL.renderMeshSimple (0.15/sc) m1 ] ,translate 2 0 $ withStrokeWidth (defaultStrokeWidth*1) $ lowerTransformations $ scale sc $ translate (negate centerX) (negate centerY) $ mkGroup [ mkGroup [] , DCEL.renderMeshColored m2 , DCEL.renderMesh (0.05/sc) m2 - -- , renderMeshEdges m - -- , DCEL.renderMeshSimple (0.10/sc) m + -- , renderMeshEdges m2 + -- , DCEL.renderMeshSimple (0.10/sc) m2 + , withFillColor "grey" $ DCEL.renderMeshSimple (0.15/sc) m2 ]] writeVar s 2 - wait (1/60) - let pipeline1 = last . take 20 . iterate - (uncurry delaunayFlip . - uncurry splitInternalEdges . - (\(a,b) -> (meshSmoothPosition a, meshSmoothPosition b))) - stages = take 30 $ iterate - (pipeline1 . - uncurry splitLongestEdge . - pipeline1 + -- wait (1/60) + let pipeline1 = last . take 20 . iterate ( + uncurry delaunayFlip . + -- uncurry splitInternalEdges . + doSmooth . + id + ) + stages = take 30 $ iterate ( + pipeline1 . + -- uncurry splitInternalEdges . + -- pipeline1 . + -- uncurry splitLongestEdge . + --pipeline1 + id ) (m1,m2) let pipeline = do modifyVar mVar (uncurry delaunayFlip) diff --git a/videos/morph/do_morph.hs b/videos/morph/do_morph.hs new file mode 100644 index 0000000..585dbab --- /dev/null +++ b/videos/morph/do_morph.hs @@ -0,0 +1,286 @@ +#!/usr/bin/env stack +-- stack runghc --package reanimate +module Main where + +import Codec.Picture.Types +import Control.Lens +import Control.Monad +import Data.List +import Data.Maybe +import qualified Data.Text as T +import Data.Tuple +import qualified Data.Vector as V +import Linear.V2 +import Linear.Vector +import qualified Numeric.LinearAlgebra as Matrix +import Numeric.LinearAlgebra.HMatrix (Matrix, linearSolve, toLists, + (><)) +import Reanimate +import Reanimate.Math.Common (barycentricCoords, isBetween, + rayIntersect) + +type Points = V.Vector (V2 Double) +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) +-- Linear interpolation on RelMesh gives smooth morph. +-- solveMesh :: RelMesh -> Mesh +-- mkRelative :: Mesh -> RelMesh +-- mkRelativePair :: MeshPair -> RelMeshPair +-- triangulate :: Polygon -> Polygon -> MeshPair ? +-- embed :: MeshPair -> MeshPair +-- compatible :: Mesh -> Mesh -> Maybe MeshPair +-- linearInterpolate :: MeshPair -> Double -> Mesh +-- convexInterpolate :: RelMeshPair -> Double -> RelMesh + +main :: IO () +main = reanimate morphAnimation + +morphAnimation :: Animation +morphAnimation = playThenReverseA $ pauseAround 1 1 $ addStatic (mkBackground "black") $ + signalA (curveS 2) $ mkAnimation 5 $ \t -> lowerTransformations $ scale 3 $ pathify $ center $ + mkGroup + [ translate (-1) 0 $ drawTrig (linearInterpolate meshPair t) + , translate 1 0 $ drawTrig $ solveMesh $ convexInterpolate relPair t + ] + where + relPair = mkRelativePair meshPair + meshPair = fromJust $ compatible example1 example2 + +mkLineP :: P -> P -> SVG +mkLineP (V2 x1 y1) (V2 x2 y2) = mkLine (x1,y1) (x2,y2) + +-- FIXME: Check that the triangles are all anticlockwise. +-- FIXME: Check that the edges connect all the points. +-- FIXME: Check that the edgse leave no gaps. +compatible :: Mesh -> Mesh -> Maybe MeshPair +compatible a b = + if meshEdges a == meshEdges b && V.length (meshPoints a) == V.length (meshPoints b) + then Just $ MeshPair (meshPoints a) (meshPoints b) (meshEdges a) + else Nothing + +linearInterpolate :: MeshPair -> Double -> Mesh +linearInterpolate (MeshPair aP bP edges) t = Mesh + { meshPoints = V.zipWith (lerp (1-t)) aP bP + , meshEdges = edges } + +example1 :: Mesh +example1 = Mesh points edges + where + points = V.fromList + [ V2 1 0 + , V2 (-1/2) (sqrt 3 / 2) + , V2 (-1/2) (-sqrt 3 / 2) + , 3 * points V.! 0 ^/ 4 + , 3 * points V.! 1 ^/ 4 + , 3 * points V.! 2 ^/ 4 + , points V.! 0 ^/ 2 + , points V.! 1 ^/ 2 + , points V.! 2 ^/ 2 + ] + edges = + [ (1,5,4), (1,2,5), (2,6,5), (2,3,6), (3,4,6), (3,1,4) + , (4,8,7), (4,5,8), (5,9,8), (5,6,9), (6,7,9), (6,4,7), (7,8,9)] + +example2 :: Mesh +example2 = Mesh points edges + where + points = V.fromList + [ V2 1 0 + , V2 (-1/2) (sqrt 3 / 2) + , V2 (-1/2) (-sqrt 3 / 2) + , 3 * points V.! 2 ^/ 4 + , 3 * points V.! 0 ^/ 4 + , 3 * points V.! 1 ^/ 4 + , points V.! 1 ^/ 2 + , points V.! 2 ^/ 2 + , points V.! 0 ^/ 2 + ] + edges = meshEdges example1 + +drawTrig (Mesh points gs) = withStrokeColor "grey" $ withFillColor "white" $ + withStrokeWidth (defaultStrokeWidth/2) $ mkGroup + [ mkGroup + [ mkGroup + [ mkLine (ax, ay) (bx, by) + , mkLine (bx, by) (cx, cy) + , mkLine (cx, cy) (ax, ay) + ] + | (a, b, c) <- gs + , let V2 ax ay = points V.! (a-1) + V2 bx by = points V.! (b-1) + V2 cx cy = points V.! (c-1) + ] + , mkGroup $ concat + [ [ colored v $ translate ax ay $ mkCircle circleRadius + , withStrokeWidth 0 $ + withStrokeColor "white" $ withFillColor "black" $ mkGroup + [ translate ax ay $ ppNum v ] + ] + | v <- nub $ concat [ [a,b,c] | (a,b,c) <- gs] + , let V2 ax ay = points V.! (v-1) + ]] + where + colored n = + let c = promotePixel $ turbo (fromIntegral n / fromIntegral (length gs-1)) + in withStrokeColorPixel c . withFillColorPixel c + ppNum n = cachedNumbers !! n + --scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}" + +cachedNumbers = + [ scaleToHeight (circleRadius*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}" + | n <- [0 .. ] ] + +circleRadius :: Double +circleRadius = 0.05 + +-- Anticlockwise. No duplicate vertices. length >= 3 +type Polygon = [V2 Double] +type P = V2 Double + +-- T = (U, G) +-- G = [Polygon] +-- U = nub $ concat G + +findStarNeighbours :: Eq a => [(a,a,a)] -> a -> [(a, a)] +findStarNeighbours allTrig self = + [ (b,c) + | (a,b,c) <- allTrig + , self == a + ] ++ + [ (c,a) + | (a,b,c) <- allTrig + , self == b + ] ++ + [ (a,b) + | (a,b,c) <- allTrig + , self == c + ] + +isInterior :: Eq a => [(a, a)] -> Bool +isInterior = isJust . getExteriorPoly + +getExteriorPoly :: Eq a => [(a, a)] -> Maybe [a] +getExteriorPoly [] = Nothing +getExteriorPoly ((a,b):rest) = worker [a] a b rest + where + worker acc start this [] = do + guard (start == this) + return (reverse acc) + worker acc start this xs = + case lookup this xs of + Just next -> worker (this:acc) start next (delete (this,next) xs) + Nothing -> + case lookup this (map swap xs) of + Just next -> worker (this:acc) start next (delete (next, this) xs) + Nothing -> Nothing + +convexInterpolate :: RelMeshPair -> Double -> RelMesh +convexInterpolate (RelMeshPair static edges leftM leftB rightM rightB) t = + RelMesh + { relMeshStatic = static + , relMeshEdges = edges + , relMeshA = Matrix.scale (1-t) leftM + + Matrix.scale t rightM + , relMeshB = Matrix.scale (1-t) leftB + + Matrix.scale t rightB + } + +solveMesh :: RelMesh -> Mesh +solveMesh (RelMesh static edges m b) = + case linearSolve m b of + Nothing -> error "Failed to solve mesh" + Just ret -> + Mesh (static <> V.fromList (worker (toLists ret))) edges + where + worker [] = [] + worker ([x]:[y]:rest) = V2 x y : worker rest + worker _ = error "invalid result" + +mkRelative :: Mesh -> RelMesh +mkRelative (Mesh points edges) = RelMesh (V.fromList exteriorPoints) edges mM bM + where + mM = (s><1) b + (s,exterior, (m, b)) = toParameters points edges + exteriorPoints = + [ points V.! (i-1) + | i <- exterior + ] + +mkRelativePair :: MeshPair -> RelMeshPair +mkRelativePair (MeshPair p1 p2 edges) = + let RelMesh static _ leftM leftB = mkRelative (Mesh p1 edges) + RelMesh _ _ rightM rightB = mkRelative (Mesh p2 edges) + in RelMeshPair static edges leftM leftB rightM rightB + +toParameters points groups = (length interior*2,exterior,unzip $ concat + [ let lst = [(if i == j then -1 else t) + | j <- interior + , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache + ] + pos = negate $ sum + [ pj ^* t + | j <- exterior + , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache + pj = points V.! (j-1) + ] + in [ (dupX lst, pos ^. _x) + , (dupY lst, pos ^. _y)] + | i <- interior ]) + where + lam_ij_cache = lam_ij points groups + dupX [] = [] + dupX (x:xs) = x:0:dupX xs + dupY [] = [] + dupY (x:xs) = 0:x:dupY xs + (interior, exterior) = + partition (isInterior . findStarNeighbours groups) [1 .. length points] + +lam_ij points groups = + [ ((i, j), t) + | i <- [1..V.length points] + , (j, t) <- lam_j points groups i + ] +lam_j points groups p = + [ (nP, sum [ t | (j,k,t) <- mu, j == nP ] / fromIntegral (length nPoints)) + | let n = findStarNeighbours groups p + nPoints = fromMaybe [] $ getExteriorPoly n + mu = calcMu points groups p + , nP <- nPoints ] +calcMu points groups p = concat + [ [ (nP, nP, t1) + , (a, nP, t2) + , (b, nP, t3) ] + -- (nP, a, b) + | {-p <- [1..length points]-} + let selfVert = points V.! (p-1) + n = findStarNeighbours groups p + nPoints :: [Int] + nPoints = fromMaybe [] $ getExteriorPoly n + , (i, nP) <- zip [1 .. ] nPoints + , let vert = points V.! (nP-1) + , let line = (vert, selfVert) + , let (a,b,aP,bP) = head $ + [ (a,b,aP,bP) + | (a,b) <- n + , let aP = points V.! (a-1) + bP = points V.! (b-1) + segment = (points V.! (a-1), points V.! (b-1)) + , case rayIntersect line segment of + Nothing -> False + Just u -> isBetween u segment + , a /= nP + , b /= nP ] + -- , b == (nPoints ++ nPoints) !! i + , let (t1,t2,t3) = barycentricCoords vert aP bP selfVert + ] diff --git a/videos/morph/least-difference.hs b/videos/morph/least-difference.hs new file mode 100644 index 0000000..0a24f40 --- /dev/null +++ b/videos/morph/least-difference.hs @@ -0,0 +1,127 @@ +#!/usr/bin/env stack +-- stack runghc --package reanimate +{-# LANGUAGE OverloadedStrings #-} +module Main where + +import Codec.Picture.Types +import Control.Lens () +import Control.Monad +import Data.Function +import Data.List +import Data.List.NonEmpty (NonEmpty) +import qualified Data.List.NonEmpty as NE +import Data.Maybe +import Data.Ratio +import qualified Data.Text as T +import Data.Tuple +import qualified Data.Vector as V +import Debug.Trace +import Linear.Matrix hiding (trace) +import Linear.Metric +import Linear.V2 +import Linear.V3 +import Linear.Vector +import Numeric.LinearAlgebra hiding (polar, scale, (<>)) +import qualified Numeric.LinearAlgebra as Matrix +import Numeric.LinearAlgebra.HMatrix hiding (polar, scale, (<>)) +import Reanimate +import Reanimate.Animation +import Reanimate.Math.Balloon +import Reanimate.Math.Common +import Reanimate.Math.Triangulate +import Reanimate.Math.Polygon +import Reanimate.Math.EarClip +import Reanimate.Math.SSSP +import Reanimate.Math.Render +import Reanimate.Math.Visibility +import Reanimate.Math.Compatible +import Reanimate.Morph.Common +import Reanimate.Morph.Linear +import Reanimate.Morph.LeastDifference +import Reanimate.PolyShape (svgToPolygons) +import Reanimate.Debug +import qualified Reanimate.Math.Compatible as Compat + +extraPoints = 0 + +-- p1 = pScale 2 $ pAtCenter $ pAddPoints extraPoints (pSetOffset shape13 0) +-- p1 = pSetOffset (addPoints 2 shape13) 0 +-- p2 = pScale 0.5 $ centerPolygon shape20 +-- p1 = centerPolygon shape2 +-- p2 = pScale 2 $ pAtCenter $ pAddPoints extraPoints (pSetOffset shape14 0 ) + +-- p1 = pAddPoints extraPoints $ (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "S") +-- p2 = pAddPoints extraPoints $ (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "C") +p1 = pAddPoints extraPoints $ (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "X") +p2 = pAddPoints extraPoints $ (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "I") + +p1_ = castPolygon p1 +p2_ = castPolygon p2 +polys = triangulate_ p1_ p2_ +p1_circ = circumference (map fst polys') p1_ +p2_circ = circumference (map snd polys') p2_ + +polys' = alignPolygons polys p1_ p2_ +polys'' = compatTriagPairs polys' + +-- (p1_circ, p2_circ) = closestLinearCorrespondence p1_circ' p2_circ' + +-- (p1s, p2s) = unzip $ triangulate p1 p2 + +main :: IO () +-- main = reanimate $ playTraces $ seq (last $ triangulate [] p1 p2) () +main = reanimate $ sceneAnimation $ do + bg <- newSpriteSVG $ mkBackground "black" + spriteZ bg (-1) + + -- let leastDiffTrig = triangulate p1 p2 + let lst = polys' + -- let lst = take 10 $ uncurry Compat.compatiblyTriangulateP (polys'!!3) -- polys'' + -- let lst = polys'' + forM_ (zip [0..] lst {-([polys''!!2]++p_snd)-}) $ \(n,(l, r)) -> do + let c = promotePixel $ turbo (n/fromIntegral (length lst-1)) + newSpriteSVG_ $ + translate (-2) 0 $ withFillColorPixel c $ mkGroup + [ polygonShape (castPolygon l) + -- , polygonNumDots (castPolygon l) + ] + newSpriteSVG_ $ + translate (2) 0 $ withFillColorPixel c $ mkGroup + [ polygonShape (castPolygon r) + -- , polygonNumDots (castPolygon r) + ] + nums <- newSpriteSVG $ mkGroup + [ translate (-2) 0 $ polygonNumDots (castPolygon l) + , translate (2) 0 $ polygonNumDots (castPolygon r) + ] + wait (1/60) + destroySprite nums + + -- fork $ play $ staticFrame (1/60) $ mkGroup + -- [ translate (2) 0 $ mkGroup + -- [ withFillColor "grey" $ polygonShape p2_circ + -- , polygonNumDots p2_circ + -- ] + -- , translate (6) 0 $ mkGroup + -- [ withFillColor "grey" $ polygonShape p2 + -- , polygonNumDots p2 + -- ] + -- ] + -- play $ staticFrame (1/60) $ mkGroup + -- [ translate (-2) 0 $ mkGroup + -- [ withFillColor "grey" $ polygonShape p1_circ + -- , polygonNumDots p1_circ + -- ] + -- , translate (-6) 0 $ mkGroup + -- [ withFillColor "grey" $ polygonShape p1 + -- , polygonNumDots p1 + -- ] + -- ] + + + return () + + +showP p = mkGroup + [ withFillColor "grey" $ polygonShape p + , polygonNumDots p ] diff --git a/videos/morph/morph.hs b/videos/morph/morph.hs index 585dbab..379b015 100644 --- a/videos/morph/morph.hs +++ b/videos/morph/morph.hs @@ -1,286 +1,3 @@ -#!/usr/bin/env stack --- stack runghc --package reanimate module Main where - -import Codec.Picture.Types -import Control.Lens -import Control.Monad -import Data.List -import Data.Maybe -import qualified Data.Text as T -import Data.Tuple -import qualified Data.Vector as V -import Linear.V2 -import Linear.Vector -import qualified Numeric.LinearAlgebra as Matrix -import Numeric.LinearAlgebra.HMatrix (Matrix, linearSolve, toLists, - (><)) -import Reanimate -import Reanimate.Math.Common (barycentricCoords, isBetween, - rayIntersect) - -type Points = V.Vector (V2 Double) -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) --- Linear interpolation on RelMesh gives smooth morph. --- solveMesh :: RelMesh -> Mesh --- mkRelative :: Mesh -> RelMesh --- mkRelativePair :: MeshPair -> RelMeshPair --- triangulate :: Polygon -> Polygon -> MeshPair ? --- embed :: MeshPair -> MeshPair --- compatible :: Mesh -> Mesh -> Maybe MeshPair --- linearInterpolate :: MeshPair -> Double -> Mesh --- convexInterpolate :: RelMeshPair -> Double -> RelMesh - main :: IO () -main = reanimate morphAnimation - -morphAnimation :: Animation -morphAnimation = playThenReverseA $ pauseAround 1 1 $ addStatic (mkBackground "black") $ - signalA (curveS 2) $ mkAnimation 5 $ \t -> lowerTransformations $ scale 3 $ pathify $ center $ - mkGroup - [ translate (-1) 0 $ drawTrig (linearInterpolate meshPair t) - , translate 1 0 $ drawTrig $ solveMesh $ convexInterpolate relPair t - ] - where - relPair = mkRelativePair meshPair - meshPair = fromJust $ compatible example1 example2 - -mkLineP :: P -> P -> SVG -mkLineP (V2 x1 y1) (V2 x2 y2) = mkLine (x1,y1) (x2,y2) - --- FIXME: Check that the triangles are all anticlockwise. --- FIXME: Check that the edges connect all the points. --- FIXME: Check that the edgse leave no gaps. -compatible :: Mesh -> Mesh -> Maybe MeshPair -compatible a b = - if meshEdges a == meshEdges b && V.length (meshPoints a) == V.length (meshPoints b) - then Just $ MeshPair (meshPoints a) (meshPoints b) (meshEdges a) - else Nothing - -linearInterpolate :: MeshPair -> Double -> Mesh -linearInterpolate (MeshPair aP bP edges) t = Mesh - { meshPoints = V.zipWith (lerp (1-t)) aP bP - , meshEdges = edges } - -example1 :: Mesh -example1 = Mesh points edges - where - points = V.fromList - [ V2 1 0 - , V2 (-1/2) (sqrt 3 / 2) - , V2 (-1/2) (-sqrt 3 / 2) - , 3 * points V.! 0 ^/ 4 - , 3 * points V.! 1 ^/ 4 - , 3 * points V.! 2 ^/ 4 - , points V.! 0 ^/ 2 - , points V.! 1 ^/ 2 - , points V.! 2 ^/ 2 - ] - edges = - [ (1,5,4), (1,2,5), (2,6,5), (2,3,6), (3,4,6), (3,1,4) - , (4,8,7), (4,5,8), (5,9,8), (5,6,9), (6,7,9), (6,4,7), (7,8,9)] - -example2 :: Mesh -example2 = Mesh points edges - where - points = V.fromList - [ V2 1 0 - , V2 (-1/2) (sqrt 3 / 2) - , V2 (-1/2) (-sqrt 3 / 2) - , 3 * points V.! 2 ^/ 4 - , 3 * points V.! 0 ^/ 4 - , 3 * points V.! 1 ^/ 4 - , points V.! 1 ^/ 2 - , points V.! 2 ^/ 2 - , points V.! 0 ^/ 2 - ] - edges = meshEdges example1 - -drawTrig (Mesh points gs) = withStrokeColor "grey" $ withFillColor "white" $ - withStrokeWidth (defaultStrokeWidth/2) $ mkGroup - [ mkGroup - [ mkGroup - [ mkLine (ax, ay) (bx, by) - , mkLine (bx, by) (cx, cy) - , mkLine (cx, cy) (ax, ay) - ] - | (a, b, c) <- gs - , let V2 ax ay = points V.! (a-1) - V2 bx by = points V.! (b-1) - V2 cx cy = points V.! (c-1) - ] - , mkGroup $ concat - [ [ colored v $ translate ax ay $ mkCircle circleRadius - , withStrokeWidth 0 $ - withStrokeColor "white" $ withFillColor "black" $ mkGroup - [ translate ax ay $ ppNum v ] - ] - | v <- nub $ concat [ [a,b,c] | (a,b,c) <- gs] - , let V2 ax ay = points V.! (v-1) - ]] - where - colored n = - let c = promotePixel $ turbo (fromIntegral n / fromIntegral (length gs-1)) - in withStrokeColorPixel c . withFillColorPixel c - ppNum n = cachedNumbers !! n - --scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}" - -cachedNumbers = - [ scaleToHeight (circleRadius*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}" - | n <- [0 .. ] ] - -circleRadius :: Double -circleRadius = 0.05 - --- Anticlockwise. No duplicate vertices. length >= 3 -type Polygon = [V2 Double] -type P = V2 Double - --- T = (U, G) --- G = [Polygon] --- U = nub $ concat G - -findStarNeighbours :: Eq a => [(a,a,a)] -> a -> [(a, a)] -findStarNeighbours allTrig self = - [ (b,c) - | (a,b,c) <- allTrig - , self == a - ] ++ - [ (c,a) - | (a,b,c) <- allTrig - , self == b - ] ++ - [ (a,b) - | (a,b,c) <- allTrig - , self == c - ] - -isInterior :: Eq a => [(a, a)] -> Bool -isInterior = isJust . getExteriorPoly - -getExteriorPoly :: Eq a => [(a, a)] -> Maybe [a] -getExteriorPoly [] = Nothing -getExteriorPoly ((a,b):rest) = worker [a] a b rest - where - worker acc start this [] = do - guard (start == this) - return (reverse acc) - worker acc start this xs = - case lookup this xs of - Just next -> worker (this:acc) start next (delete (this,next) xs) - Nothing -> - case lookup this (map swap xs) of - Just next -> worker (this:acc) start next (delete (next, this) xs) - Nothing -> Nothing - -convexInterpolate :: RelMeshPair -> Double -> RelMesh -convexInterpolate (RelMeshPair static edges leftM leftB rightM rightB) t = - RelMesh - { relMeshStatic = static - , relMeshEdges = edges - , relMeshA = Matrix.scale (1-t) leftM + - Matrix.scale t rightM - , relMeshB = Matrix.scale (1-t) leftB + - Matrix.scale t rightB - } - -solveMesh :: RelMesh -> Mesh -solveMesh (RelMesh static edges m b) = - case linearSolve m b of - Nothing -> error "Failed to solve mesh" - Just ret -> - Mesh (static <> V.fromList (worker (toLists ret))) edges - where - worker [] = [] - worker ([x]:[y]:rest) = V2 x y : worker rest - worker _ = error "invalid result" - -mkRelative :: Mesh -> RelMesh -mkRelative (Mesh points edges) = RelMesh (V.fromList exteriorPoints) edges mM bM - where - mM = (s><1) b - (s,exterior, (m, b)) = toParameters points edges - exteriorPoints = - [ points V.! (i-1) - | i <- exterior - ] - -mkRelativePair :: MeshPair -> RelMeshPair -mkRelativePair (MeshPair p1 p2 edges) = - let RelMesh static _ leftM leftB = mkRelative (Mesh p1 edges) - RelMesh _ _ rightM rightB = mkRelative (Mesh p2 edges) - in RelMeshPair static edges leftM leftB rightM rightB - -toParameters points groups = (length interior*2,exterior,unzip $ concat - [ let lst = [(if i == j then -1 else t) - | j <- interior - , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache - ] - pos = negate $ sum - [ pj ^* t - | j <- exterior - , let t = fromMaybe 0 $ lookup (i,j) lam_ij_cache - pj = points V.! (j-1) - ] - in [ (dupX lst, pos ^. _x) - , (dupY lst, pos ^. _y)] - | i <- interior ]) - where - lam_ij_cache = lam_ij points groups - dupX [] = [] - dupX (x:xs) = x:0:dupX xs - dupY [] = [] - dupY (x:xs) = 0:x:dupY xs - (interior, exterior) = - partition (isInterior . findStarNeighbours groups) [1 .. length points] - -lam_ij points groups = - [ ((i, j), t) - | i <- [1..V.length points] - , (j, t) <- lam_j points groups i - ] -lam_j points groups p = - [ (nP, sum [ t | (j,k,t) <- mu, j == nP ] / fromIntegral (length nPoints)) - | let n = findStarNeighbours groups p - nPoints = fromMaybe [] $ getExteriorPoly n - mu = calcMu points groups p - , nP <- nPoints ] -calcMu points groups p = concat - [ [ (nP, nP, t1) - , (a, nP, t2) - , (b, nP, t3) ] - -- (nP, a, b) - | {-p <- [1..length points]-} - let selfVert = points V.! (p-1) - n = findStarNeighbours groups p - nPoints :: [Int] - nPoints = fromMaybe [] $ getExteriorPoly n - , (i, nP) <- zip [1 .. ] nPoints - , let vert = points V.! (nP-1) - , let line = (vert, selfVert) - , let (a,b,aP,bP) = head $ - [ (a,b,aP,bP) - | (a,b) <- n - , let aP = points V.! (a-1) - bP = points V.! (b-1) - segment = (points V.! (a-1), points V.! (b-1)) - , case rayIntersect line segment of - Nothing -> False - Just u -> isBetween u segment - , a /= nP - , b /= nP ] - -- , b == (nPoints ++ nPoints) !! i - , let (t1,t2,t3) = barycentricCoords vert aP bP selfVert - ] +main = return () diff --git a/videos/morph/rigid.hs b/videos/morph/rigid.hs new file mode 100644 index 0000000..7676192 --- /dev/null +++ b/videos/morph/rigid.hs @@ -0,0 +1,309 @@ +#!/usr/bin/env stack +-- stack runghc --package reanimate +{-# LANGUAGE OverloadedStrings #-} +module Main where + +import Codec.Picture.Types +import Control.Exception +import Control.Lens () +import Control.Monad +import Data.Function +import Data.List +import Data.List.NonEmpty (NonEmpty) +import qualified Data.List.NonEmpty as NE +import Data.Maybe +import Data.Ratio +import qualified Data.Text as T +import Data.Tuple +import qualified Data.Vector as V +import Debug.Trace +import Linear.Matrix hiding (trace) +import Linear.Metric +import Linear.V2 +import Linear.V3 +import Linear.Vector +import Numeric.LinearAlgebra hiding (polar, scale, (<>)) +import qualified Numeric.LinearAlgebra as Matrix +import Numeric.LinearAlgebra.HMatrix hiding (polar, scale, (<>)) +import Reanimate +import Reanimate.Builtin.Documentation +import Reanimate.Debug +import Reanimate.Math.Balloon +import Reanimate.Math.Common +import Reanimate.Math.Compatible (compatiblyTriangulateP) +import qualified Reanimate.Math.DCEL as DCEL +import Reanimate.Math.EarClip +import Reanimate.Math.Polygon +import Reanimate.Math.Render +import Reanimate.Math.Smooth +import Reanimate.Math.SSSP +import Reanimate.Math.Visibility +import Reanimate.Morph.Common +import Reanimate.Morph.LeastDifference +import Reanimate.Morph.LeastWork +import Reanimate.Morph.Linear +import Reanimate.Morph.Rigid +import Reanimate.PolyShape (svgToPolygons) +import Text.Printf + +-- p1 = centerPolygon $ shape2 +--p1 = pScale 3.5 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0) +-- p1 = setOffset (addPoints 2 shape13) 0 +-- p2 = scalePolygon 0.5 $ centerPolygon shape20 +-- p1 = centerPolygon shape2 +--p2 = pScale 3.5 $ pAtCenter $ pAddPoints 0 (pSetOffset shape14 0) +-- p2 = setOffset shape14 0 +p1 = pCopy $ pSetOffset (pCopy p1') 0 +p2 = pCopy $ pSetOffset (pCopy p2') 0 +(p1', p2') = + -- normalizePolygons + -- closestLinearCorrespondence + (,) + -- leastWork zeroStretchCosts defaultBendCosts + -- leastWork defaultStretchCosts defaultBendCosts + -- (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 6 $ latex "S") + -- (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 6 $ latex "C") + (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "X") + (pAtCenter $ unsafeSVGToPolygon 0.1 $ scale 6 $ latex "I") + +p1_ = castPolygon p1 +p2_ = castPolygon p2 +polys = triangulate_ p1_ p2_ +p1_circ = castPolygon (circumference (map fst polys') p1_) +p2_circ = castPolygon (circumference (map snd polys') p2_) + +polys' = alignPolygons polys p1_ p2_ +(p1s, p2s) = unzip $ compatTriagPairs [ (castPolygon a, castPolygon b) | (a,b) <- polys' ] + +-- (p1s,p2s) = unzip (compatiblyTriangulateP p1 p2) + +m2 = DCEL.buildMesh $ DCEL.polygonsMesh + (map (fmap realToFrac) $ V.toList $ polygonPoints p2_circ) + (map (map (fmap realToFrac) . V.toList . polygonPoints) p2s) + +m1 = DCEL.buildMesh $ DCEL.polygonsMesh + (map (fmap realToFrac) $ V.toList $ polygonPoints p1_circ) + (map (map (fmap realToFrac) . V.toList . polygonPoints) p1s) + +pipeline1 = last . take 20 . iterate + ( uncurry DCEL.delaunayFlip . + uncurry DCEL.splitInternalEdges . + (\(a,b) -> (DCEL.meshSmoothPosition a, DCEL.meshSmoothPosition b))) +(m1final, m2final) = last $ take 30 $ iterate + (pipeline1 . + uncurry DCEL.splitLongestEdge . + pipeline1 + ) (m1,m2) + +-- (m1good, m2good) = last $ take 9 $ iterate +-- (pipeline1 . +-- -- uncurry DCEL.splitLongestEdge . +-- pipeline1 +-- ) (m1,m2) + +main :: IO () +-- main = do +-- evaluate (last $ take 900 $ compatiblyTriangulateP p1 p2) +-- return () +-- main = reanimate $ playTraces $ length $ take 29 p1s +-- main = reanimate $ playThenReverseA $ pauseAround 0.5 0.5 $ sceneAnimation $ do + +main = reanimate $ sceneAnimation $ do + newSpriteSVG_ $ mkBackgroundPixel rtfdBackgroundColor + -- play + -- $ playTraces + -- $ traceA (mapA (withFillOpacity 0 . withStrokeColor "white") drawCircle) + -- $ traceA (mapA (withFillOpacity 0 . withStrokeColor "white") drawBox) + -- $ traceSVG ((withFillOpacity 0 . withStrokeColor "white") $ mkCircle 2) + -- $ traceSVG ((withFillOpacity 0 . withStrokeColor "white") $ mkRect 2 2) + -- $ 20 + -- play $ playTraces $ last $ take 35 $ compatiblyTriangulateP p1 p2 + -- fork $ newSpriteA $ drawCompatible p1 p2 + -- newSpriteSVG_ $ translate (5) 0 $ lowerTransformations $ scale 2 $ + -- -- DCEL.renderMesh 0.02 m1final + -- DCEL.renderMeshColored m1final + -- renderMeshPair optMesh + -- newSpriteSVG_ $ translate (-5) 0 $ lowerTransformations $ scale 2 $ + -- DCEL.renderMesh 0.02 m1good + + -- newSpriteSVG_ $ + -- translate (-3) 0 $ scale 1 $ mkGroup + -- [ mkGroup [] + -- , withFillColor "grey" $ polygonShape p1 + -- , withFillColor "grey" $ polygonNumDots p1 + -- ] + -- newSpriteSVG_ $ + -- translate (3) 0 $ scale 1 $ mkGroup + -- [ mkGroup [] + -- , withFillColor "grey" $ polygonShape p2 + -- , withFillColor "grey" $ polygonNumDots p2 + -- ] + -- wait 1 + -- fork $ forM_ p1s $ \p1Piece -> do + -- newSpriteSVG_ $ + -- translate (-3) 0 $ scale 1 $ mkGroup + -- [ mkGroup [] + -- , withFillColor "white" $ polygonShape p1Piece + -- -- , withFillColor "grey" $ polygonNumDots p1Piece + -- ] + -- wait (1/60) + -- -- wait 1 + -- fork $ forM_ p2s $ \p1Piece -> do + -- newSpriteSVG_ $ + -- translate (3) 0 $ scale 1 $ mkGroup + -- [ mkGroup [] + -- , withFillColor "white" $ polygonShape p1Piece + -- -- , withFillColor "grey" $ polygonNumDots p1Piece + -- ] + -- wait (1/60) + -- -- wait 1 + + -- fork $ play $ mkAnimation 3 $ \t -> + -- let points = interpolate prep t + -- in translate (-3) 0 $ scale 1 $ mkGroup + -- [ mkGroup [] + -- , drawTrigs points (meshTriangles myMesh) + -- -- , DCEL.renderMeshColored m1 + -- ] + fork $ play $ mkAnimation 3 $ \t -> + let points = interpolate optPrep t + in translate (3) 0 $ lowerTransformations $ scale 1 $ mkGroup + [ mkGroup [] + , drawTrigs points (meshTriangles optMesh) + -- , DCEL.renderMeshColored m1final + -- , withGroupOpacity 0.5 $ DCEL.renderMesh 0.02 m1final + -- , withGroupOpacity 0.5 $ withFillOpacity 0 $ renderMeshPair optMesh + -- , let V2 x y = V.head points in + -- translate x y $ withFillColor "green" $ + -- mkCircle 0.03 + ] + -- , scale 2 $ polygonNumDots p1 + -- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t -> + -- let points = interpolate bestPrep t + -- in translate (3) 0 $ mkGroup [drawTrigs points (meshTriangles bestMesh) + -- -- , scale 2 $ polygonNumDots p2 + -- ] + -- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t -> + -- let points = interpolate prepRev t in + -- mkGroup + -- [ drawTrigs points (meshTriangles myMeshRev) + -- -- , polygonNumDots shape2 + -- ] + -- forM_ (myMesh : smoothMesh myMesh) $ \newMesh -> do + -- let (minAngA, minAngB) = meshMinAngle newMesh + -- txt <- newSpriteSVG $ withFillColor "white" $ + -- mkGroup + -- [translate (5) 2 $ center $ + -- latex $ T.pack $ printf "Min: %.1f" (minAngA/pi*180) + -- ,translate (5) 1 $ center $ + -- latex $ T.pack $ printf "Min: %.1f" (minAngB/pi*180) ] + -- s <- newSpriteSVG $ renderAMesh newMesh + -- wait (recip 60) + -- destroySprite s + -- destroySprite txt + -- newSpriteSVG $ translate (-3) 0 $ scale 2 $ renderMesh p1 (map fst trigs) + -- newSpriteSVG $ translate 3 0 $ scale 2 $ renderMesh p2 (map snd trigs) + -- wait 1 + where + prep = prepare myMesh + optPrep = prepare optMesh + -- bestPrep = prepare bestMesh + myMesh = mkMesh p1 p2 + optMesh = toRigidMesh m1final m2final + -- optMesh = DCEL.toRigidMesh m1good m2good + -- bestMesh = last $ smoothMesh myMesh + prepRev = prepare myMeshRev + myMeshRev = reverseMesh myMesh + + + +reverseMesh :: Mesh -> Mesh +reverseMesh mesh = + mesh { meshPointsA = meshPointsB mesh, meshPointsB = meshPointsA mesh } + +mkMesh :: Polygon -> Polygon -> Mesh +mkMesh a b = Mesh + { meshPointsA = V.map (fmap realToFrac) pointsA + , meshPointsB = V.map (fmap realToFrac) pointsB + , meshOutline = V.map (fromJust . flip V.elemIndex pointsA) (polygonPoints a) + , meshTriangles = relTrigs + } + where + pointsA = + V.fromList $ nub $ V.toList $ V.concat [ polygonPoints a | (a, b) <- trigs ] + pointsB = + V.fromList $ nub $ V.toList $ V.concat [ polygonPoints b | (a, b) <- trigs ] + trigs = compatiblyTriangulateP a b + mkRel arr p = + ( (fromJust $ V.elemIndex (pAccess p 0) arr) + , (fromJust $ V.elemIndex (pAccess p 1) arr) + , (fromJust $ V.elemIndex (pAccess p 2) arr) + ) + relTrigs = V.fromList [ (mkRel pointsA a) | (a, b) <- trigs ] + -- data Mesh = Mesh (Vector P) (Vector (RelTrig, RelTrig)) + +testMesh2 :: Mesh +testMesh2 = mkMesh (pTranslate (V2 0 0) shape2) shape20 + +testPrep = prepare testMesh2 + +drawTrigs :: V.Vector (V2 Double) -> V.Vector RelTrig -> SVG +drawTrigs points trigs = mkGroup + [ mkGroup + [ withFillOpacity 1 + $ withStrokeWidth (defaultStrokeWidth * 0) + $ withStrokeColor "grey" + $ withFillColor "black" + $ drawPolygon + $ map (points V.!) [a, b, c] + | (a, b, c) <- V.toList trigs + ] + -- , a == 0 + -- , mkGroup + -- [ withFillColor "red" $ mkGroup + -- [ drawPoint (points V.! a), drawPoint (points V.! b), drawPoint (points V.! c)] + -- | (a,b,c) <- V.toList trigs + -- -- , a == 0 + -- ] + ] + +drawTrigsLines :: V.Vector (V2 Double) -> V.Vector RelTrig -> SVG +drawTrigsLines points trigs = withFillOpacity 0 $ withStrokeColor "black" $ mkGroup + [ mkLinePathClosed + [ (aPx, aPy) + , (bPx, bPy) + , (cPx, cPy)] + | (a,b,c) <- V.toList trigs + -- , a==2 || b == 2 || c == 2 + , let V2 aPx aPy = points V.! a + V2 bPx bPy = points V.! b + V2 cPx cPy = points V.! c + ] + +drawPoint :: V2 Double -> SVG +drawPoint (V2 x y) = translate x y $ mkCircle 0.1 + +drawPolygon :: [V2 Double] -> SVG +drawPolygon lst = mkLinePathClosed [ (x, y) | V2 x y <- lst ] + +drawCompatible :: Polygon -> Polygon -> Animation +drawCompatible a b = sceneAnimation $ do + let left = -6 + right = 4 + newSpriteSVG $ translate left 0 $ mkGroup + [ withFillColor "grey" $ polygonShape a + , withFillColor "grey" $ polygonNumDots a + ] + newSpriteSVG $ translate right 0 $ mkGroup + [ withFillColor "grey" $ polygonShape b + , withFillColor "grey" $ polygonNumDots b + ] + return () + -- let compat = compatiblyTriangulateP a b + -- forM_ compat $ \(l, r) -> do + -- fork $ play $ staticFrame 1 $ + -- translate left 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $ + -- withFillOpacity 0 $ polygonShape l + -- fork $ play $ staticFrame 1 $ + -- translate right 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $ + -- withFillOpacity 0 $ polygonShape r