hashirama/reanimate
Watch
1
0
Fork
You've already forked reanimate
 0
mirror of https://github.com/reanimate/reanimate.git synced 2026-09-09 15:12:21 +00:00
Code Issues Projects Releases Wiki Activity

Feature morph (#62)

Browse source 
* Disable windows CI (because of hmatrix)
* Implement ear-clip triangulation.
* Implement common framework for polygon morphing.
* Draft implementations of compatible triangulation and mesh smoothing.
* Implement least-work morphing.
* Implement as-rigid-as-possible morphing.
* Implement SSSP.
* Draft article about polygon morphing.

Former-commit-id: f6f19f58b9dd58c655bcca0ed2d65f72726f0a06
This commit is contained in:
David Himmelstrup 2020-05-14 15:14:29 +08:00 • committed by GitHub
commit 4356efbc7e
70 changed files with 6504 additions and 195 deletions
Download patch file Download diff file Expand all files Collapse all files

19
.azure/azure-linux-template.yml
View file

@ -10,12 +10,12 @@ jobs:
BUILD: stack
STACK_YAML: stack.yaml
ARGS: --pedantic
stack-lts-12:
stack-lts-14:
BUILD: stack
STACK_YAML: stack-lts-12.yaml
stack-lts-11:
STACK_YAML: stack-lts-14.yaml
stack-lts-13:
BUILD: stack
STACK_YAML: stack-lts-11.yaml
STACK_YAML: stack-lts-13.yaml
#stack-nightly:
# BUILD: stack
# ARGS: --resolver nightly
@ -28,11 +28,12 @@ jobs:
cacheHitVar: CACHE_RESTORED
displayName: Cache stack root
- script: |
sudo add-apt-repository ppa:jonathonf/ffmpeg-4
#sudo add-apt-repository ppa:jonathonf/ffmpeg-4
sudo apt-get update
sudo apt-get install ffmpeg
sudo apt-get install libblas-dev liblapack-dev
#sudo apt-get install ffmpeg
#sudo apt-get install texlive texlive-latex-base texlive-latex-extra texlive-fonts-extra texlive-science texlive-xetex texlive-latex-recommended texlive-lang-english texlive-lang-chinese
sudo apt-get install povray
#sudo apt-get install povray
mkdir -p ~/.local/bin
curl -L https://get.haskellstack.org/stable/linux-x86_64.tar.gz | tar xz --wildcards --strip-components=1 -C ~/.local/bin '*/stack'
case "$BUILD" in
@ -68,7 +69,7 @@ jobs:
;;
*)
export PATH=$HOME/.local/bin:$PATH
stack --install-ghc test $ARGS --only-dependencies
stack --install-ghc test $ARGS --flag reanimate:test --only-dependencies
;;
esac
set -ex
@ -102,7 +103,7 @@ jobs:
;;
*)
stack ./examples/counter.hs check
stack test $ARGS --flag reanimate:test
stack test --flag reanimate:test --test-arguments="-j2"
;;
esac
set +ex

14
.azure/azure-osx-template.yml
View file

@ -10,12 +10,12 @@ jobs:
BUILD: stack
STACK_YAML: stack.yaml
ARGS: --pedantic
stack-lts-12:
stack-lts-14:
BUILD: stack
STACK_YAML: stack-lts-12.yaml
stack-lts-11:
STACK_YAML: stack-lts-14.yaml
stack-lts-13:
BUILD: stack
STACK_YAML: stack-lts-11.yaml
STACK_YAML: stack-lts-13.yaml
maxParallel: 6
steps:
- task: Cache@2
@ -25,12 +25,14 @@ jobs:
cacheHitVar: CACHE_RESTORED
displayName: Cache stack root
- script: |
# brew install gcc
# brew install lapack
mkdir -p ~/.local/bin
curl -skL https://get.haskellstack.org/stable/osx-x86_64.tar.gz | tar xz --strip-components=1 --include '*/stack' -C ~/.local/bin;
export PATH=$HOME/.local/bin:$PATH
stack --install-ghc test $ARGS --only-dependencies
stack --install-ghc test $ARGS --flag reanimate:test --only-dependencies
stack ./examples/counter.hs check
stack test $ARGS --flag reanimate:test
stack test --flag reanimate:test --test-arguments="-j2"
env:
OS_NAME: ${{ parameters.os }}
displayName: 'Installation ${{parameters.os}} & Test'

61
.azure/azure-windows-template.yml
View file

@ -9,28 +9,69 @@ jobs:
matrix:
stack:
BUILD: stack
STACK_YAML: stack.yaml
STACK_YAML: stack-lts-14.yaml # lts-15 is broken. "Access violation in generated code" :-/
ARGS: --pedantic
stack-lts-12:
stack-lts-14:
BUILD: stack
STACK_YAML: stack-lts-12.yaml
stack-lts-11:
STACK_YAML: stack-lts-14.yaml
stack-lts-13:
BUILD: stack
STACK_YAML: stack-lts-11.yaml
STACK_YAML: stack-lts-13.yaml
maxParallel: 6
steps:
- task: Cache@2
inputs:
key: ${{ parameters.name }} | ${{ parameters.vmImage }} | $(STACK_YAML) | stack-root | $(Agent.OS)
key: ${{ parameters.name }} | ${{ parameters.vmImage }} | $(STACK_YAML) | stack-root | $(Agent.OS) | version 2
path: $(STACK_ROOT)
cacheHitVar: CACHE_RESTORED
displayName: Cache stack root
- bash: |
set -o xtrace
#curl -sSkL https://icl.cs.utk.edu/lapack-for-windows/libraries/VisualStudio/3.7.0/Dynamic-MINGW/Win64/libblas.dll -o libblas.dll
# curl -sSkL https://icl.cs.utk.edu/lapack-for-windows/libraries/VisualStudio/3.7.0/Dynamic-MINGW/Win64/libblas.lib -o libblas.lib
#curl -sSkL https://icl.cs.utk.edu/lapack-for-windows/libraries/VisualStudio/3.7.0/Dynamic-MINGW/Win64/liblapack.dll -o liblapack.dll
# curl -sSkL https://icl.cs.utk.edu/lapack-for-windows/libraries/VisualStudio/3.7.0/Dynamic-MINGW/Win64/liblapack.lib -o liblapack.lib
# cp libblas.dll libblas.lib liblapack.dll liblapack.lib /usr/lib/
# cp libblas.dll libblas.lib liblapack.dll liblapack.lib /lib/
# cp libblas.dll libblas.lib liblapack.dll liblapack.lib /usr/local/lib/
# cp libblas.dll libblas.lib liblapack.dll liblapack.lib /c/mingw810/prerequisites/x86_64-binutils-nomulti/x86_64-w64-mingw32/lib
# cp libblas.dll libblas.lib liblapack.dll liblapack.lib /c/mingw810/prerequisites/x86_64-binutils-nomulti/lib
echo $CONDA
echo $PATH
CWD=`pwd`
cd "$CONDA/condabin/"
export PATH=`pwd`:$PATH
cd $CWD
echo $PATH
echo `which conda.bat`
echo `which conda`
conda info --envs
ls $CONDA/condabin
# $CONDA/condabin/conda.bat create --yes --name reanimate
# $CONDA/condabin/conda.bat activate reanimate
$CONDA/condabin/conda.bat
$CONDA/condabin/conda.bat install --yes --quiet -c conda-forge libblas
$CONDA/condabin/conda.bat install --yes --quiet -c conda-forge liblapack
find "$CONDA" -name '*libblas*'
find "$CONDA" -name '*lapack*'
# export LIBRARY_PATH="`pwd`;$LIBRARY_PATH"
curl -sSkL http://www.stackage.org/stack/windows-x86_64 -o /usr/bin/stack.zip
unzip -o /usr/bin/stack.zip -d /usr/bin/
stack --install-ghc test $ARGS --only-dependencies
stack ./examples/counter.hs check
stack test $ARGS --flag reanimate:test
env:
OS_NAME: ${{ parameters.os }}
displayName: 'Installation ${{parameters.os}} & Test'
displayName: 'Setup'
- bash: |
set -o xtrace
stack --install-ghc build $ARGS --no-keep-going --fast --extra-lib-dirs="$CONDA/Library/bin/" --only-dependencies
stack --install-ghc test $ARGS --no-keep-going --fast --flag reanimate:test --extra-lib-dirs="$CONDA/Library/bin/" --only-dependencies
env:
OS_NAME: ${{ parameters.os }}
displayName: 'Install dependencies'
- bash: |
set -o xtrace
stack ./examples/counter.hs check
rm examples/*.golden
#stack test --flag reanimate:test --test-arguments="-j2" --fast --extra-lib-dirs="$CONDA/Library/bin/"
env:
OS_NAME: ${{ parameters.os }}
displayName: 'Build & Test'

1
README.md
View file

@ -73,6 +73,7 @@ animation source code, the browser window will automatically reload and show the
# Documentation
* API reference: http://hackage.haskell.org/package/reanimate/docs/Reanimate.html
* Core concepts: https://reanimate.readthedocs.io/en/latest/introduction/
* Design overview: https://reanimate.readthedocs.io/en/latest/glue_tut/
* Gallery with source code: https://reanimate.readthedocs.io/en/latest/gallery/

85
docs/index.md
View file

@ -1,3 +1,10 @@
# Reanimate is:
* **An animation library** - a system for turning code into movies, similar to [3b1b's manim](https://github.com/3b1b/manim) of [YouTube fame](https://www.youtube.com/channel/UCYO_jab_esuFRV4b17AJtAw).
* **User friendly** - extensively documented with example animations covering the [entire API](http://hackage.haskell.org/package/reanimate/docs/Reanimate.html).
* **Cross-platform** - works on Linux, Mac and Windows.
# Getting started
## Viewing a basic animation
@ -6,19 +13,87 @@
$ git clone https://github.com/Lemmih/reanimate.git
$ cd reanimate/
$ stack build
$ stack ./examples/raster.hs
$ stack ./examples/doc_andThen.hs
```
Running the above code should open a browser window containing this animation:
<details>
<summary>Toggle doc_andThen.hs source code.</summary>
<pre><code class="haskell">
{!examples/doc_andThen.hs!}
</code></pre>
</details>
<br/>
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/bvLqalX.mp4">
</video>
## Checking dependencies
Reanimate includes a run-time for optional dependencies.
If something doesn't work, you can run a check for run-time dependencies like
this:
```console
$ stack ./examples/raster.hs check
$ stack ./examples/doc_andThen.hs check
reanimate checks:
Has ffmpeg: 4.1.3-0ubuntu1
Has dvisvgm: /usr/bin/dvisvgm
Has povray: /usr/bin/povray
Has blender: 2.82
Has rsvg-convert: 2.44.10
Has inkscape: 0.92.4
Has convert: 6.9.10-14
Has LaTeX: /usr/bin/latex
Has LaTeX package 'babel': OK
Has LaTeX package 'preview': OK
Has LaTeX package 'amsmath': OK
Has XeLaTeX: /usr/bin/xelatex
Has XeLaTeX package 'ctex': OK
```
## Rendering movies
None of these dependencies are vital but the functionality of **reanimate**
will be reduced if they are missing. For example, without 'ffmpeg', you won't
be able to generate stand-alone video files (but can still view animations
in a browser window). Without 'LaTeX' and 'dvisvgm', you won't be able use
LaTeX for typesetting.
## Rendering animations
**reanimate** has builtin support for rendering animations to video files and
gifs. To render an animation with default options, use the 'render' command:
```console
$ stack ./examples/raster.hs render
$ stack ./examples/doc_andThen.hs render
```
For computationally intense animations, use the `--compile` flag to compile
your animation script and run it using all available cores:
```console
$ stack ./examples/doc_andThen.hs render --compile
```
If the defaults are too aggressive for you tastes, change them directly or use
a different preset:
```console
$ stack ./examples/doc_andThen.hs render --preset high --fps 10
```
Rendering to gifs automatically does two passes and uses a color palette to
improve the quality and size of the gif. GIFs are, by default, rendered at
24 fps with a resolution of 320x180:
```console
$ stack ./examples/doc_andThen.hs render --format gif
```
# Learning resources
* [API documentation](http://hackage.haskell.org/package/reanimate/docs/Reanimate.html)
* [Core concepts](introduction)
* [Gallery](gallery)
* [Examples](https://github.com/Lemmih/reanimate/tree/master/examples)
* [Design overview](glue_tut)

2
docs/introduction.md
View file

@ -1,4 +1,4 @@
# Core concepts
# Introduction
This document will introduce the foundamental concepts used in **reanimate**.
It is assumed that you are already familiar with Haskell. After reading, you will be able to understand all the concepts used in the [API reference documentation](http://hackage.haskell.org/package/reanimate/docs/Reanimate.html). This is not a tutorial, though, and you may have to look at the examples before you can turn these concepts into beautiful animations.

200
docs/morphology.md Normal file
View file

@ -0,0 +1,200 @@
# 2D morphology
Morphing 2D shapes is conceptually simple. We can imagine it as a function that
takes two SVG images and produces the intermediate steps that smoothly
transforms the source image to the target image. If the SVG images are simple
shapes then it looks like this:
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/4Uwayst.mp4">
</video>
Often it is quite easy for a human to tell if a morph looks right or not. This,
unfortunately, does not mean the problem is easy to solve for a computer. As
it turns out, there is not a single, perfect algorithm for morphing between 2D
shapes. Instead, there are many different algorithms that each have their own
benefits and trade-offs. This article will specify the details of the
morphology-problem and show a handful of different solutions.
# Breaking it down
Conceptually, a morphing function should take two SVG images and give a new
SVG image over time. In Haskell, it would look like this:
```haskell
morph :: SVG → SVG → (Time → SVG)
```
However, morphing has several sub-problems, each of which that can be solved
indepently in a variety of ways. These solutions have different benefits and
drawbacks, and wildly different performance. Since no single morphing algorithm
is obviously superior in all cases, let's identify the othogonal sub-problems
and use their solutions as parameters to our morphing function.
## 1. Point correspondence
```haskell
type PointCorrespondence = Polygon → Polygon → (Polygon, Polygon)
```
SVG images can be simplified to polygons which are made out of points. To morph
between two shapes, we somehow have to transform the points in a source
polygon into the points in a target polygon. Finding a good point correspondence
between two polygons is a tricky thing, though, and it can have a huge impact
on the quality of the morph as seen in this illustration:
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/9WJ6mC6.mp4">
</video>
Point correspondence algorithms take two polygons (which may have a different
number of points) and align their points so there's an exact one-to-one
correspondence. This may require shuffling, adding, splitting, or merging points.
## 2. Point trajectory
```haskell
type Trajectory = (Polygon, Polygon) → (Double → Polygon)
```
Once the points in the source and target polygons have been aligned, the next
step is to figure out which path they should take. The simplest approach would
be to move the points in a straight line but this doesn't always look right.
The points in the leftmost figure of the following illustration moves in a
straight line while the points in the rightmost figure take a curved path
designed to avoid self-intersections:
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/jXOR7Ij.mp4">
</video>
The design space for trajectory algorithms is quite large and this article
will cover five different approaches in addition to the simplest straight-line
path.
## 3. Object correspondence
```haskell
type ObjectCorrespondence = [Polygon] → [Polygon] → [(Polygon, Polygon)]
```
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/GoiZxgo.mp4">
</video>
## 4. Color interpolation
Morphing between colors is something many libraries don't pay a lot of attention
but there is a fair amount of complexity here and it does deserve consideration.
In short, humans (usually) have three types light sensitive cells that respond
to different ranges of light wavelengths. So, even though the light we see may
contain hundreds of different wavelengths, the eye only notices how much each
of the three types of cells are stimulated and will therefore always reduce
colors to just three numbers.
So far, so good. With three coordinates necessary to specify a color, each color
has a position in a three dimensional space and we can morph between two colors
by going in a straight line. This is where things get tricky, though, because
going in a straight line in the XYZ space (the 3D colorspace defined by the
sensitivity of a typical human's eye) is physically correct but often a bit
counter intuitive. For example, morphing between **blue** and **yellow** will go
through **grey**. Also, the XYZ doesn't take into account how humans perceive
colors (ie. what happens in the brain rather than the eye) and
morphing between two colors with the same perceived brightness will not keep the
brightness a constant. These concerns, and many more, have lead to the
development of different 3D colorspaces, such as CIE LAB, which take human color perception into
account. **Reanimate** has built-in support for several spaces and this
illustration shows how they stack up against each other:
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/AUqFi7L.mp4">
</video>
By default, **reanimate** uses the LAB colorspace for morphing colors.
# Linear interpolation
TBD.
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/ZBh1ena.mp4">
</video>
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/2jWIXLL.mp4">
</video>
# Rotational interpolation
```haskell
rotationalTrajectory :: Origin -> Trajectory
```
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/dFY0IZz.mp4">
</video>
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/nAk8xJ1.mp4">
</video>
# Edge+angle interpolation
```haskell
lineBend :: Trajectory
```
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/cmG1Wwr.mp4">
</video>
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/hoSncAC.mp4">
</video>
TBD.
# Least-work correspondence
<video width="640" height="360" muted autoplay loop>
<source src="https://i.imgur.com/C8bCnea.mp4">
</video>
$$
work_{stretch} = k_s\frac{\delta^2}{2L_0}
$$
$$
work_{bend} =
\left\{
\begin{array}{ll}
k_b(\Delta \theta + m_b \Delta \theta^\star)^{e_b}
& \text{if $\theta(t)$ never goes to zero} \\
k_b(\Delta \theta + m_b \Delta \theta^\star)^{e_b} + p_b
& \text{if $\theta(t)$ does go to zero}
\end{array}
\right.
$$
TBD.
# Skeleton
TBD.
# As-rigid-as-possible
TBD.
# Guaranteed intersection-free
TBD.
# Related work
* Flubber
* Polymorph
* raphael
* Kute
* svg_morph
* D3

1
docs/requirements.txt
View file

@ -1 +1,2 @@
markdown_include==0.5.1
pymdown-extensions==7.0

1
examples/bug_play1.hs
View file

@ -3,7 +3,6 @@
module Main(main) where
import Reanimate
import Reanimate.Scene
-- This should give circle.
main :: IO ()

1
examples/bug_play2.hs
View file

@ -3,7 +3,6 @@
module Main(main) where
import Reanimate
import Reanimate.Scene
-- This should give rect.
main :: IO ()

1
examples/bug_play3.hs
View file

@ -3,7 +3,6 @@
module Main(main) where
import Reanimate
import Reanimate.Scene
-- This should give rect.
main :: IO ()

1
examples/doc_queryNow.hs
View file

@ -5,7 +5,6 @@ module Main(main) where
import Reanimate
import Reanimate.Builtin.Documentation
import Data.Monoid ((<>))
import qualified Data.Text as T
main :: IO ()

1
examples/fourier.hs
View file

@ -8,7 +8,6 @@ import qualified Data.Text as T
import Graphics.SvgTree
import Reanimate
import Reanimate.LaTeX (latex)
waveMultiplier :: Int
-- waveMultiplier = 1 -- Sawtooth wave

61
examples/morphology_closest.hs Executable file
View file

@ -0,0 +1,61 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Control.Monad
import Graphics.SvgTree (LineJoin (..))
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-4) 4 $
center $ latex "no-op"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (4) 4 $
center $ latex "closest"
forM_ pairs $ uncurry showPair
where
showPair from to =
waitOn $ do
fork $ play $ mkAnimation 4 (morph rawLinear from to)
# mapA (translate (-4) (-0.5))
# signalA (curveS 4)
fork $ play $ mkAnimation 4 (morph linear from to)
# mapA (translate (4) (-0.5))
# signalA (curveS 4)
pairs = zip stages (tail stages ++ [head stages])
stages = map (lowerTransformations . scale 8 . pathify . center) $ colorize
[ latex "X"
, latex "$\\aleph$"
, latex "Y"
, latex "$\\infty$"
, latex "I"
, latex "$\\pi$"
, latex "1"
, latex "S"
, mkRect 0.5 0.5
]
colorize :: [SVG] -> [SVG]
colorize lst =
[ withFillColorPixel (promotePixel $ parula (n/fromIntegral (length lst-1))) elt
| elt <- lst
| n <- [0..]
]

45
examples/morphology_color.hs Executable file
View file

@ -0,0 +1,45 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Reanimate
import Reanimate.Interpolate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
setDuration 5 $
mapA (withStrokeWidth 0) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
doMorph yellow blue
doMorph blue yellow
where
doMorph from to =
play $ pauseAtEnd 0.5 $ animate $ \t ->
mkGroup
[ translate (-screenWidth/2 + screenWidth/6*(n*2+1)) 0 $ mkGroup
[ translate 0 3 $ scale 1.5 $ center $ latex label
, translate 0 (-0.5) $ morph conf from to t
]
| conf <- [ linear
, linear{morphColorComponents=xyzComponents}
, linear{morphColorComponents=hsvComponents}]
| label <- ["LAB", "XYZ", "HSV"]
| n <- [0..]
]
radius = 2.5
yellow = withFillColor "yellow" $ mkCircle radius
blue = withFillColor "blue" $ mkCircle radius
-- green = withFillColor "green" $ mkCircle radius
-- pink = withFillColor "pink" $ mkCircle radius
-- red = withFillColor "red" $ mkCircle radius
-- purple = withFillColor "purple" $ mkCircle radius

36
examples/morphology_intro.hs Executable file
View file

@ -0,0 +1,36 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth 0) $
mapA (withStrokeColor "black") $
mapA (withFillOpacity 1) $
sceneAnimation $ do
play $ step stage1 stage2
play $ step stage2 stage3
play $ step stage3 stage4
play $ step stage4 stage1
# setDuration 5
where
radius = 2.5
step from to =
signalA (curveS 4) $ animate $ morph linear from to
stage1 = translate (-3) 0 $ withFillColor "red" $ mkCircle radius
stage2 = translate 3 0 $ withFillColor "blue" $ mkRect (radius*2) (radius*2)
stage3 = mkGroup
[translate (-1) (-1) $ withFillColor "green" $ mkRect (radius*0.5) (radius*0.5)
,translate 1 (1) $ withFillColor "black" $ mkRect (radius*0.5) (radius*0.5) ]
stage4 = translate (-3) 0 $ withFillColor "purple" $ mkCircle (radius*0.24)

69
examples/morphology_leastwork.hs Executable file
View file

@ -0,0 +1,69 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Control.Monad
import Graphics.SvgTree (LineJoin (..))
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.LeastWork
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-3) 4 $
center $ latex "linear"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (3) 4 $
center $ latex "least-work"
forM_ pairs $ uncurry showPair
where
showPair from to =
waitOn $ do
fork $ play $ mkAnimation 4 (morph linear from to)
# mapA (translate (-3) (-0.5))
# signalA (curveS 4)
fork $ play $ mkAnimation 4 (morph myMorph from to)
# mapA (translate (3) (-0.5))
# signalA (curveS 4)
stretchCosts = defaultStretchCosts
{ stretchStiffness = 2 }
bendCosts = defaultBendCosts
myMorph = linear{morphPointCorrespondence = leastWork stretchCosts bendCosts }
pairs = zip stages (tail stages ++ [head stages])
stages = map (lowerTransformations . scale 6 . pathify . center) $ colorize
[ latex "X"
, latex "$\\aleph$"
, latex "Y"
, latex "$\\infty$"
, latex "I"
, latex "$\\pi$"
, latex "1"
, latex "T"
, latex "I"
, latex "L"
, latex "S"
, mkRect 0.5 0.5
]
colorize :: [SVG] -> [SVG]
colorize lst =
[ withFillColorPixel (promotePixel $ parula (n/fromIntegral (length lst-1))) elt
| elt <- lst
| n <- [0..]
]

60
examples/morphology_linear.hs Executable file
View file

@ -0,0 +1,60 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Graphics.SvgTree (LineJoin (..))
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $
showPair (stages ++ take 1 stages)
where
showPair (from:to:rest) =
waitOn $ do
toS <- newSpriteSVG $
translate (4) 0 to
spriteE toS $ overBeginning 0.2 fadeInE
spriteE toS $ overEnding 0.2 fadeOutE
m <- fork $ newSpriteA $ animate (morph linear from to)
# signalA (curveS 4)
# mapA (translate (-4) 0)
wait 4
destroySprite m
destroySprite toS
showPair (to:rest)
return ()
showPair _ = return ()
stages = map (lowerTransformations . scale 8 . pathify . center) $ colorize
[ latex "X"
, latex "$\\aleph$"
, latex "Y"
, latex "$\\infty$"
, latex "I"
, latex "$\\pi$"
, latex "1"
, latex "S"
, mkRect 0.5 0.5
]
colorize :: [SVG] -> [SVG]
colorize lst =
[ withFillColorPixel (promotePixel $ parula (n/fromIntegral (length lst-1))) elt
| elt <- lst
| n <- [0..]
]

63
examples/morphology_linebend.hs Executable file
View file

@ -0,0 +1,63 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Control.Monad
import Graphics.SvgTree (LineJoin (..))
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
import Reanimate.Morph.LineBend
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-3) 4 $
center $ latex "linear"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (3) 4 $
center $ latex "line bend"
forM_ pairs $ uncurry showPair
where
showPair from to =
waitOn $ do
fork $ play $ mkAnimation 4 (morph linear from to)
# mapA (translate (-3) (-0.5))
# signalA (curveS 4)
fork $ play $ mkAnimation 4 (morph myMorph from to)
# mapA (translate (3) (-0.5))
# signalA (curveS 4)
myMorph = linear{morphTrajectory = lineBend }
pairs = zip stages (tail stages ++ [head stages])
stages = map (lowerTransformations . scale 6 . pathify . center) $ colorize
[ latex "X"
, latex "$\\aleph$"
, latex "Y"
, latex "$\\infty$"
, latex "I"
, latex "$\\pi$"
, latex "1"
, latex "S"
, mkRect 0.5 0.5
]
colorize :: [SVG] -> [SVG]
colorize lst =
[ withFillColorPixel (promotePixel $ parula (n/fromIntegral (length lst-1))) elt
| elt <- lst
| n <- [0..]
]

100
examples/morphology_linebend_intro.hs Executable file
View file

@ -0,0 +1,100 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Control.Monad
import Data.List (transpose)
import qualified Data.Vector as V
import Graphics.SvgTree (LineJoin (..))
import Linear.V2
import Reanimate
import Reanimate.Math.Polygon
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
import Reanimate.Morph.LineBend
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
spike1 :: Polygon
spike1 = pScale 2.5 $ mkPolygon $ V.map (\x -> x - V2 1 1) $ V.fromList
[ V2 0 0, V2 2 0
, V2 2 1 ]
spike2 :: Polygon
spike2 = pScale 2.5 $ mkPolygon $ V.map (\x -> x - V2 1 1) $ V.fromList
[ V2 0 0, V2 2 0
, V2 0 1 ]
spike3 :: Polygon
spike3 = pScale 2.5 $ mkPolygon $ V.map (\x -> x - V2 1 1) $ V.fromList
[ V2 0.5 0, V2 1.5 0
, V2 1 2 ]
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-4) 4 $
center $ latex "linear"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (4) 4 $
center $ latex "line bend"
showTrails
forM_ pairs $ uncurry showPair
where
showTrails = do
_ <- newSpriteSVG $
translate (-4) (-0.5) $
genTrails (map linearTrajectory pairs)
_ <- newSpriteSVG $
translate (4) (-0.5) $
genTrails (map lineBend pairs)
return ()
showPair from to =
waitOn $ do
fork $ play $ mkAnimation 4 (morph linear (polygonShape from) (polygonShape to))
# mapA (translate (-4) (-0.5))
# mapA (withFillColor "lightgreen")
# signalA (curveS 4)
fork $ play $ mkAnimation 4 (morph myMorph (polygonShape from) (polygonShape to))
# mapA (translate (4) (-0.5))
# mapA (withFillColor "cyan")
# signalA (curveS 4)
myMorph = linear{morphTrajectory = lineBend }
pairs = zip stages (tail stages ++ [head stages])
stages =
[ spike1
, spike2
, spike3
]
genTrails :: [(Double -> Polygon)] -> SVG
genTrails plotters =
withFillOpacity 0 $
withStrokeWidth (defaultStrokeWidth*0.5) $
withStrokeColor "black" $
withStrokeDashArray [0.1,0.05] $
mkGroup $ map mkTrail $ transpose $ concat
[
[ V.toList $ V.map (fmap realToFrac) (polygonPoints poly)
| n <- [0..steps]
, let poly = plotter (fromIntegral n / fromIntegral steps)
]
| plotter <- plotters ]
where
steps = 100 :: Int
mkTrail lst = mkLinePath [ (x,y) | V2 x y <- lst ]

59
examples/morphology_object_correspondence.hs Executable file
View file

@ -0,0 +1,59 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth 0) $
mapA (withStrokeColor "black") $
mapA (withFillOpacity 1) $
sceneAnimation $ do
let showLabel label = do
fork $ play $ staticFrame 4 (center $ latex label)
# mapA (translate 0 4)
# applyE (overBeginning 0.1 fadeInE)
# applyE (overEnding 0.1 fadeOutE)
let sides = mkGroup
[ translate (-4) 0 $ withFillColor "blue" $
mkCircle 2
, translate 4 0 $ withFillColor "green" $
mkCircle 2]
middle =
withFillColor "red" $
mkRect 4 4
showLabel "Cut"
play $ step sides middle
# pauseAtEnd 1
play $ step middle sides
# pauseAtEnd 1
showLabel "Overlap"
play $ stepDup sides middle
# pauseAtEnd 1
play $ stepDup middle sides
# pauseAtEnd 1
showLabel "Obliterate"
play $ stepGenesis sides middle
# pauseAtEnd 1
play $ stepGenesis middle sides
# pauseAtEnd 1
where
step from to =
signalA (curveS 2) $ animate $ morph linear from to
stepDup from to =
signalA (curveS 2) $ animate $ morph linear{morphObjectCorrespondence=dupObjectCorrespondence} from to
stepGenesis from to =
signalA (curveS 2) $ animate $ morph linear{morphObjectCorrespondence=genesisObjectCorrespondence} from to

100
examples/morphology_point_correspondence.hs Executable file
View file

@ -0,0 +1,100 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Control.Monad (replicateM_)
import qualified Data.Text as T
import qualified Data.Vector as V
import Linear.V2
import Linear.Vector
import Reanimate
import Reanimate.Math.Polygon
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth 0) $
mapA (withStrokeColor "black") $
mapA (withFillOpacity 1) $
sceneAnimation $ do
let pl1 = translate (-4) 0 $ mkGroup
[ lowerTransformations $ scale 3 $ withFillOpacity 1 $
withStrokeColor "black" $
withStrokeWidth defaultStrokeWidth $
withFillColor "lightgreen" $
polygonShape octogon
, lowerTransformations $ scale 3 $ polygonNumDots octogon 0
]
pl2 n t = translate 4 0 $ lowerTransformations $ scale 3 $ mkGroup
[ withFillOpacity 1 $
withStrokeColor "black" $
withStrokeWidth defaultStrokeWidth $
withFillColor "cyan" $
polygonShape (pSetOffset star n)
, polygonNumDots (pSetOffset star n) t ]
offset <- newVar 0
slide <- newVar 0
_ <- newSprite $ pl2 <$> unVar offset <*> unVar slide
_ <- newSpriteSVG $ pl1
let slideLeft = do
tweenVar slide 1 $ \v -> fromToS v 1 . curveS 4
writeVar slide 0
modifyVar offset succ
replicateM_ 8 $ do
offsetVal <- readVar offset
play $ step pl1 (pl2 offsetVal 0)
# setDuration 3
slideLeft
where
step from to =
signalA (curveS 2) $ animate $ morph rawLinear from to
octogon :: Polygon
octogon = mkPolygon $ V.fromList
[ realToFrac <$> V2 (cos phi) (sin phi)
| n <- [0..7]
, let phi = n*2*pi/8 + pi/8 :: Double
]
star :: Polygon
star = pScale 0.5 $ mkPolygon $ V.fromList
[ V2 0 1, V2 (-2) 2, V2 (-1) 0
, V2 (-2) (-2), V2 0 (-1), V2 2 (-2)
, V2 1 0, V2 2 2 ]
genColor :: Int -> Int -> PixelRGBA8
genColor n m =
promotePixel $ parula (fromIntegral (n+offset) / fromIntegral (m+offset))
where
offset = 5
polygonNumDots :: Polygon -> Double -> SVG
polygonNumDots p t = mkGroup $ reverse
[ mkGroup
[ colored n $ withStrokeWidth (defaultStrokeWidth*0.5) $ withStrokeColor "black" $
translate x y $ pathify $ mkCircle circR
, withFillColor "black" $
translate x y $ ppNum n ]
| n <- [0..pSize p-1]
, let a = realToFrac <$> pAccess p n
b = realToFrac <$> pAccess p (pNext p n)
V2 x y = lerp t b a ]
where
circR = 0.1
colored n =
let c = genColor n (pSize p-1)
in withFillColorPixel c
ppNum n = scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]

320
examples/morphology_point_trajectory.hs Normal file
View file

@ -0,0 +1,320 @@
#!/usr/bin/env stack
-- stack --resolver lts-15.04 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
import Reanimate.Math.Polygon
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
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
meshToPolygon :: Mesh -> Polygon
meshToPolygon mesh = pScale 2 $ mkPolygon $ V.fromList
[ realToFrac <$> (meshPoints mesh V.! (n-1))
| n <- polygonNodes ]
where
polygonNodes :: [Int]
polygonNodes = [5,8,9,7,4,6]
main :: IO ()
main = reanimate morphAnimation
genTrails :: (Double -> Mesh) -> SVG
genTrails mkMesh =
withFillOpacity 0 $
withStrokeWidth (defaultStrokeWidth*0.5) $
withStrokeColor "black" $
withStrokeDashArray [0.1,0.1] $
mkGroup $ map mkTrail $ transpose
[ V.toList $ V.map (fmap realToFrac) $ polygonPoints $ meshToPolygon mesh
| n <- [0..steps]
, let mesh = mkMesh (fromIntegral n / fromIntegral steps)
]
where
steps = 100 :: Int
mkTrail lst =mkLinePath [ (x,y) | V2 x y <- lst ]
morphAnimation :: Animation
morphAnimation = playThenReverseA $ pauseAround 1 1 $ addStatic (mkBackgroundPixel bgColor) $
signalA (curveS 2) $ mkAnimation 5 $ \t -> lowerTransformations $ scale 3 $ pathify $
mkGroup
[ translate (-1.5) 0 $ withFillColor "lightgreen" $ mkGroup
[ withStrokeWidth defaultStrokeWidth $
withStrokeColor "black" $
polygonShape $ meshToPolygon $ linearInterpolate meshPair t
, linearTrails
, polygonNumDots (meshToPolygon $ linearInterpolate meshPair t) 0
]
, translate 1 0 $ withFillColor "cyan" $ mkGroup
[ withStrokeWidth defaultStrokeWidth $
withStrokeColor "black" $
polygonShape $ meshToPolygon $ solveMesh $ convexInterpolate relPair t
, curveTrails
, polygonNumDots (meshToPolygon $ solveMesh $ convexInterpolate relPair t) 0
]
]
where
linearTrails = genTrails (linearInterpolate meshPair)
curveTrails = genTrails (solveMesh . convexInterpolate relPair)
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
-- 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><s) (concat m)
bM = (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 :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> (Int, [Int], ([[b]], [b]))
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 :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> [((Int, Int), b)]
lam_ij points groups =
[ ((i, j), t)
| i <- [1..V.length points]
, (j, t) <- lam_j points groups i
]
lam_j :: (Ord b, Fractional b) => V.Vector (V2 b) -> [(Int, Int, Int)] -> Int -> [(Int, b)]
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 :: (Ord c, Fractional c) => V.Vector (V2 c) -> [(Int, Int, Int)] -> Int -> [(Int, Int, c)]
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
, nP <- 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
]
genColor :: Int -> Int -> PixelRGBA8
genColor n m =
promotePixel $ parula (fromIntegral (n+offset) / fromIntegral (m+offset))
where
offset = 5
polygonNumDots :: Polygon -> Double -> SVG
polygonNumDots p t = mkGroup $ reverse
[ mkGroup
[ colored n $ withStrokeWidth (defaultStrokeWidth*0.5) $ withStrokeColor "black" $
translate x y $ mkCircle circR
, withFillColor "black" $
translate x y $ ppNum n ]
| n <- [0..pSize p-1]
, let a = realToFrac <$> pAccess p n
b = realToFrac <$> pAccess p (pNext p n)
V2 x y = lerp t b a ]
where
circR = 0.1
colored n =
let c = genColor n (pSize p-1)
in withFillColorPixel c
ppNum n = scaleToHeight (circR*1.5) $ center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]

63
examples/morphology_rotational.hs Executable file
View file

@ -0,0 +1,63 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Codec.Picture.Types
import Control.Monad
import Graphics.SvgTree (LineJoin(..))
import Reanimate
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
import Reanimate.Morph.Rotational
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-3) 4 $
center $ latex "linear"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (3) 4 $
center $ latex "rotational"
forM_ pairs $ uncurry showPair
where
showPair from to =
waitOn $ do
fork $ play $ mkAnimation 4 (morph linear from to)
# mapA (translate (-3) (-0.5))
# signalA (curveS 4)
fork $ play $ mkAnimation 4 (morph myMorph from to)
# mapA (translate (3) (-0.5))
# signalA (curveS 4)
myMorph = linear{morphTrajectory = rotationalTrajectory origin }
origin = (0.5, 0.5)
pairs = zip stages (tail stages ++ [head stages])
stages = map (lowerTransformations . scale 6 . pathify . center) $ colorize
[ latex "X"
, latex "$\\aleph$"
, latex "Y"
, latex "$\\infty$"
, latex "I"
, latex "$\\pi$"
, latex "1"
, latex "S"
, mkRect 0.5 0.5
]
colorize :: [SVG] -> [SVG]
colorize lst =
[ withFillColorPixel (promotePixel $ parula (n/fromIntegral (length lst-1))) elt
| elt <- lst
| n <- [0..]
]

125
examples/morphology_rotational_intro.hs Executable file
View file

@ -0,0 +1,125 @@
#!/usr/bin/env stack
-- stack runghc --package reanimate
{-# LANGUAGE ApplicativeDo #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE ParallelListComp #-}
module Main(main) where
import Codec.Picture
import Data.List (transpose)
import qualified Data.Vector as V
import Graphics.SvgTree (LineJoin (..))
import Linear.V2
import Linear.Vector
import Reanimate
import Reanimate.Math.Polygon
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
import Reanimate.Morph.Rotational
bgColor :: PixelRGBA8
bgColor = PixelRGBA8 252 252 252 0xFF
spike1 :: Polygon
spike1 = pScale 2.5 $ mkPolygon $ V.map (\x -> x - V2 1 1) $ V.fromList
[ V2 0 0, V2 2 0
, V2 2 1 ]
spike2 :: Polygon
spike2 = pScale 2.5 $ mkPolygon $ V.map (\x -> x - V2 1 1) $ V.fromList
[ V2 0 0, V2 2 0
, V2 0 1 ]
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]
main :: IO ()
main = reanimate $
addStatic (mkBackgroundPixel bgColor) $
mapA (withStrokeWidth defaultStrokeWidth) $
mapA (withStrokeColor "black") $
mapA (withStrokeLineJoin JoinRound) $
mapA (withFillOpacity 1) $
sceneAnimation $ do
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (-4) 4 $
center $ latex "linear"
_ <- newSpriteSVG $
withStrokeWidth 0 $ translate (4) 4 $
center $ latex "rotational"
originVar <- newVar (0.5, 0.5)
showTrails originVar
let pushLeft newOrigin = do
fork $ tweenVar originVar 0.5 $ \(x,y) t ->
let t' = curveS 2 t in
(fromToS x (fst newOrigin) t', fromToS y (snd newOrigin) t')
showPair originVar spike1 spike2
pushRight newOrigin = do
fork $ tweenVar originVar 0.5 $ \(x,y) t ->
let t' = curveS 2 t in
(fromToS x (fst newOrigin) t', fromToS y (snd newOrigin) t')
showPair originVar spike2 spike1
pushLeft (1, 0.5)
pushRight (0.1, 0.95)
pushLeft (0.5, 1)
pushRight (0.5, 2)
pushLeft (-0.3, -0.2)
pushRight (0.5, 0.5)
where
showTrails originVar = do
_ <- newSpriteSVG $
translate (-4) (-0.5) $
genTrails (map linearTrajectory [(spike1, spike2)])
_ <- newSprite $ do
origin <- unVar originVar
pure $
translate (4) (-0.5) $
genTrails (map (rotationalTrajectory origin) [(spike1, spike2)])
return ()
showPair originVar from to =
waitOn $ do
fork $ play $ mkAnimation 2 (morph linear (polygonShape from) (polygonShape to))
# mapA (translate (-4) (-0.5))
# mapA (withFillColor "lightgreen")
# signalA (curveS 4)
s <- fork $ newSprite $ do
origin <- unVar originVar
t <- spriteT
d <- spriteDuration
pure $
let localTime = curveS 4 (t/d)
myMorph = linear{morphTrajectory = rotationalTrajectory origin }
originSrc = polygonOrigin from origin
originDst = polygonOrigin to origin
V2 originX originY = lerp localTime originDst originSrc
in
translate (4) (-0.5) $ mkGroup
[ withFillColor "cyan" $
morph myMorph (polygonShape from) (polygonShape to) localTime
, withFillColor "red" $ translate originX originY $
mkCircle 0.1
]
wait 2
destroySprite s
genTrails :: [(Double -> Polygon)] -> SVG
genTrails plotters =
withFillOpacity 0 $
withStrokeWidth (defaultStrokeWidth*0.5) $
withStrokeColor "black" $
withStrokeDashArray [0.1,0.05] $
mkGroup $ map mkTrail $ transpose $ concat
[
[ V.toList $ V.map (fmap realToFrac) (polygonPoints poly)
| n <- [0..steps]
, let poly = plotter (fromIntegral n / fromIntegral steps)
]
| plotter <- plotters ]
where
steps = 100 :: Int
mkTrail lst = mkLinePath [ (x,y) | V2 x y <- lst ]

1
examples/signals.hs
View file

@ -3,7 +3,6 @@
{-# LANGUAGE OverloadedStrings #-}
module Main (main) where
import Data.Monoid ((<>))
import Data.Text (Text, pack)
import Graphics.SvgTree hiding (Text)
import Numeric

10
hie.yaml Normal file
View file

@ -0,0 +1,10 @@
cradle:
multi:
- path: "./src"
config: {cradle: {stack: {component: "reanimate:lib"}} }
- path: "./videos"
config: {cradle: {stack: {component: "reanimate:lib"}} }
- path: "./examples"
config: {cradle: {stack: {component: "reanimate:lib"}} }
- path: "./test"
config: {cradle: {none: }}

9
mkdocs.yml
View file

@ -4,5 +4,14 @@ theme:
name: readthedocs
markdown_extensions:
- markdown_include.include
- pymdownx.arithmatex
extra_javascript:
- https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.0/MathJax.js?config=TeX-MML-AM_CHTML
extra_css:
- fix.css
nav:
- Reanimate: index.md
- Core concepts: introduction.md
- Gluing together animations: glue_tut.md
- Morphology (incomplete): morphology.md
- Gallery (outdated): gallery.md

41
reanimate.cabal
View file

@ -36,6 +36,11 @@ Source-Repository head
Type: git
Location: git://github.com/lemmih/reanimate.git
Flag hmatrix
Description: Enable hmatrix dependency (requires blas and lapack)
Default: True
Manual: True
Flag test
Description: Enable testing (requires stack)
Default: False
@ -58,6 +63,24 @@ library
-- Reanimate.Diagrams
Reanimate.Transform
Reanimate.Misc
Reanimate.Math.Common
Reanimate.Math.Polygon
Reanimate.Math.Triangulate
Reanimate.Math.DCEL
Reanimate.Math.EarClip
Reanimate.Math.Visibility
Reanimate.Math.SSSP
Reanimate.Math.Balloon
Reanimate.Math.Compatible
Reanimate.Math.Render
Reanimate.Math.Smooth
Reanimate.Morph.Common
Reanimate.Morph.Linear
Reanimate.Morph.Rotational
Reanimate.Morph.LeastWork
Reanimate.Morph.LineBend
Reanimate.Morph.Cache
Reanimate.Morph.Rigid
Reanimate.Raster
Reanimate.ColorMap
Reanimate.ColorSpace
@ -89,12 +112,16 @@ library
time, text, filepath, process, directory,
containers, reanimate-svg >= 0.9.7.0, xml, bytestring, lens, linear, mtl, matrix,
JuicyPixels, attoparsec, parallel,
cubicbezier, websockets,
cubicbezier, websockets >= 0.12.7.0,
hashable, fsnotify, open-browser, random-shuffle, base64-bytestring,
vector >= 0.12.0.0, colour, cassava, ansi-wl-pprint, here, temporary,
optparse-applicative, chiphunk >= 0.1.2.1,
geojson, aeson >= 1.3.0.0
ghc-options: -Wall
geojson, aeson >= 1.3.0.0,
split, fingertree, array, cereal, vector-space
if flag(hmatrix) {
build-depends: hmatrix >= 0.20.0.0
}
ghc-options: -Wall -fno-ignore-asserts
test-suite spec
type: exitcode-stdio-1.0
@ -102,16 +129,20 @@ test-suite spec
default-language: Haskell2010
other-modules:
UnitTests
Properties
Helpers
hs-source-dirs: test
if !flag(test) {
Buildable: False
}
ghc-options: -w
build-depends:
base,
directory,
filepath,
bytestring, process,
reanimate,
QuickCheck >= 2.1.0,
reanimate, vector, linear,
QuickCheck >= 2.1.0, text,
tasty, tasty-golden, tasty-hunit,
tasty-quickcheck,
ansi-wl-pprint

25
src/Reanimate/Cache.hs
View file

@ -5,10 +5,12 @@ module Reanimate.Cache
, cacheDiskSvg
, cacheDiskKey
, cacheDiskLines
, encodeInt
) where
import Control.Exception
import Control.Monad (unless)
import Data.Bits
import Data.Hashable
import Data.IORef
import Data.Map (Map)
@ -40,11 +42,11 @@ cacheFile template gen = do
renameOrCopyFile tmp path
evaluate path
cacheDisk :: (T.Text -> Maybe a) -> (a -> T.Text) -> (Text -> IO a) -> (Text -> IO a)
cacheDisk parse render gen key = do
cacheDisk :: String -> (T.Text -> Maybe a) -> (a -> T.Text) -> (Text -> IO a) -> (Text -> IO a)
cacheDisk cacheType parse render gen key = do
root <- getXdgDirectory XdgCache "reanimate"
createDirectoryIfMissing True root
let path = root </> show (hash key)
let path = root </> encodeInt (hash key) <.> cacheType
hit <- doesFileExist path
if hit
then do
@ -55,7 +57,7 @@ cacheDisk parse render gen key = do
else genCache root path
where
genCache root path = do
(tmpPath, tmpHandle) <- openTempFile root (show (hash key))
(tmpPath, tmpHandle) <- openTempFile root (encodeInt (hash key))
new <- gen key
T.hPutStr tmpHandle (render new)
hClose tmpHandle
@ -66,7 +68,7 @@ cacheDiskKey :: Text -> IO Tree -> IO Tree
cacheDiskKey key gen = cacheDiskSvg (const gen) key
cacheDiskSvg :: (Text -> IO Tree) -> (Text -> IO Tree)
cacheDiskSvg = cacheDisk parse render
cacheDiskSvg = cacheDisk "svg" parse render
where
parse txt = case parseXML txt of
[Elem t] -> Just (unparse t)
@ -74,7 +76,7 @@ cacheDiskSvg = cacheDisk parse render
render = T.pack . renderTree
cacheDiskLines :: (Text -> IO [Text]) -> (Text -> IO [Text])
cacheDiskLines = cacheDisk parse render
cacheDiskLines = cacheDisk "txt" parse render
where
parse = Just . T.lines
render = T.unlines
@ -96,3 +98,14 @@ cacheMem gen key = do
-- don't store the result.
None -> pure None
_ -> atomicModifyIORef cache (\m -> (Map.insert key svg m, svg))
encodeInt :: Int -> String
encodeInt i = worker (fromIntegral i) 60
where
worker :: Word -> Int -> String
worker key sh
| sh < 0 = []
| otherwise =
case (key `shiftR` sh) `mod` 64 of
idx -> alphabet !! fromIntegral idx : worker key (sh-6)
alphabet = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+$"

102
src/Reanimate/Driver.hs
View file

@ -1,17 +1,23 @@
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE MultiWayIf #-}
module Reanimate.Driver ( reanimate ) where
module Reanimate.Driver
( reanimate
)
where
import Control.Applicative ((<|>))
import Control.Applicative ( (<|>) )
import Data.Maybe
import Reanimate.Animation (Animation)
import Reanimate.Animation ( Animation )
import Reanimate.Driver.Check
import Reanimate.Driver.CLI
import Reanimate.Driver.Compile
import Reanimate.Driver.Server
import Reanimate.Parameters
import Reanimate.Render (FPS, Format (..), Height, Width,
render, renderSnippets, renderSvgs, selectRaster)
import Reanimate.Render ( render
, renderSnippets
, renderSvgs
, selectRaster
)
import System.Directory
import System.FilePath
import Text.Printf
@ -95,18 +101,19 @@ Rendering animation can be controlled with these arguments:
-}
reanimate :: Animation -> IO ()
reanimate animation = do
Options{..} <- getDriverOptions
Options {..} <- getDriverOptions
case optsCommand of
Raw -> setFPS 60 >> renderSvgs animation
Test -> do
Raw -> setFPS 60 >> renderSvgs animation
Test -> do
setNoExternals True
-- hSetBinaryMode stdout True
renderSnippets animation
Check -> checkEnvironment
View -> serve
Render{..} -> do
let fmt = guessParameter renderFormat (fmap presetFormat renderPreset) $
case renderTarget of
Check -> checkEnvironment
View -> serve
Render {..} -> do
let fmt =
guessParameter renderFormat (fmap presetFormat renderPreset)
$ case renderTarget of
-- Format guessed from output
Just target -> case takeExtension target of
".mp4" -> RenderMp4
@ -125,38 +132,54 @@ reanimate animation = do
RenderWebm -> replaceExtension self "webm"
Just target -> makeAbsolute target
let fps = guessParameter renderFPS (fmap presetFPS renderPreset) $ formatFPS fmt
(width, height) =
fromMaybe
( maybe (formatWidth fmt) presetWidth renderPreset
, maybe (formatHeight fmt) presetHeight renderPreset )
(userPreferredDimensions renderWidth renderHeight)
let
fps =
guessParameter renderFPS (fmap presetFPS renderPreset) $ formatFPS fmt
(width, height) = fromMaybe
( maybe (formatWidth fmt) presetWidth renderPreset
, maybe (formatHeight fmt) presetHeight renderPreset
)
(userPreferredDimensions renderWidth renderHeight)
if renderCompile
then
compile
["render"
,"--fps", show fps
,"--width", show width
,"--height", show height
,"--format", showFormat fmt
,"--raster", showRaster renderRaster
,"--target", target
,"+RTS", "-N", "-RTS"]
then compile
[ "render"
, "--fps"
, show fps
, "--width"
, show width
, "--height"
, show height
, "--format"
, showFormat fmt
, "--raster"
, showRaster renderRaster
, "--target"
, target
, "+RTS"
, "-N"
, "-RTS"
]
else do
raster <- selectRaster renderRaster
setRaster raster
setFPS fps
setWidth width
setHeight height
printf "Animation options:\n\
printf
"Animation options:\n\
\ fps: %d\n\
\ width: %d\n\
\ height: %d\n\
\ fmt: %s\n\
\ target: %s\n\
\ raster: %s\n"
fps width height (showFormat fmt) target (show raster)
fps
width
height
(showFormat fmt)
target
(show raster)
render animation target raster fmt width height fps
@ -166,13 +189,14 @@ guessParameter a b def = fromMaybe def (a <|> b)
-- If user specifies exactly one dimension explicitly, calculate the other
userPreferredDimensions :: Maybe Width -> Maybe Height -> Maybe (Width, Height)
userPreferredDimensions (Just width) (Just height) = Just (width, height)
userPreferredDimensions (Just width) Nothing = Just (width, makeEven $ width * 9 `div` 16)
userPreferredDimensions Nothing (Just height) = Just (makeEven $ height * 16 `div` 9, height)
userPreferredDimensions Nothing Nothing = Nothing
userPreferredDimensions (Just width) (Just height) = Just (width, height)
userPreferredDimensions (Just width) Nothing =
Just (width, makeEven $ width * 9 `div` 16)
userPreferredDimensions Nothing (Just height) =
Just (makeEven $ height * 16 `div` 9, height)
userPreferredDimensions Nothing Nothing = Nothing
-- Avoid ffmpeg failures "height not divisible by 2"
makeEven :: Int -> Int
makeEven x
| even x = x
| otherwise = x - 1
makeEven x | even x = x
| otherwise = x - 1

7
src/Reanimate/Driver/Server.hs
View file

@ -47,7 +47,12 @@ serve = withManager $ \watch -> do
unless hasConn openViewer
putStrLn "Listening..."
runServerWith "127.0.0.1" 9161 opts $ \pending -> do
let options = ServerOptions
{ serverHost = "127.0.0.1"
, serverPort = 9161
, serverConnectionOptions = opts
, serverRequirePong = Nothing }
runServerWithOptions options $ \pending -> do
putStrLn "New connection received."
hasConn <- swapMVar hasConnectionVar True
if hasConn

29
src/Reanimate/Interpolate.hs
View file

@ -1,18 +1,20 @@
{-# LANGUAGE RecordWildCards #-}
module Reanimate.Interpolate where
import Codec.Picture
import Data.Colour
import Data.Colour.CIE
import Data.Colour.CIE.Illuminant (d65)
import Data.Colour.SRGB
import Data.Colour.RGBSpace.HSV
import Data.Colour.RGBSpace
import Data.Fixed
import Codec.Picture
import Codec.Picture.Types
import Data.Colour
import Data.Colour.CIE
import Data.Colour.CIE.Illuminant (d65)
import Data.Colour.RGBSpace
import Data.Colour.RGBSpace.HSV
import Data.Colour.SRGB
import Data.Fixed
import Reanimate.Signal
data ColorComponents = ColorComponents
{ colorUnpack :: Colour Double -> (Double, Double, Double)
, colorPack :: Double -> Double -> Double -> Colour Double }
, colorPack :: Double -> Double -> Double -> Colour Double }
-- | > interpolate rgbComponents yellow blue
--
@ -79,6 +81,15 @@ interpolate ColorComponents{..} from to = \d ->
interpolateRGB8 :: ColorComponents -> PixelRGB8 -> PixelRGB8 -> (Double -> PixelRGB8)
interpolateRGB8 comps from to = toRGB8 . interpolate comps (fromRGB8 from) (fromRGB8 to)
interpolateRGBA8 :: ColorComponents -> PixelRGBA8 -> PixelRGBA8 -> (Double -> PixelRGBA8)
interpolateRGBA8 comps from to = \t ->
case interp t of
PixelRGB8 r g b ->
let alpha = fromToS (fromIntegral $ pixelOpacity from) (fromIntegral $ pixelOpacity to) t
in PixelRGBA8 r g b (round alpha)
where
interp = interpolateRGB8 comps (dropTransparency from) (dropTransparency to)
toRGB8 :: Colour Double -> PixelRGB8
toRGB8 c = PixelRGB8 r g b
where

1
src/Reanimate/LaTeX.hs
View file

@ -3,7 +3,6 @@
module Reanimate.LaTeX (latex,xelatex,latexAlign) where
import qualified Data.ByteString as B
import Data.Monoid ((<>))
import Data.Text (Text)
import qualified Data.Text as T
import qualified Data.Text.IO as T

167
src/Reanimate/Math/Balloon.hs Normal file
View file

@ -0,0 +1,167 @@
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE OverloadedStrings #-}
module Reanimate.Math.Balloon where
import Control.Lens
import qualified Data.Vector as V
import Graphics.SvgTree (drawAttributes)
import Linear.V2
import Linear.Vector
import Reanimate.Animation
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Morph.Common (toShapes)
import Reanimate.Svg.Constructors
-- import Debug.Trace
balloon :: SVG -> (Double -> SVG)
balloon = balloon' 0.01
balloon' :: Double -> SVG -> (Double -> SVG)
balloon' tol svg = \t ->
mkGroup
[ polygonShape (gen t) & drawAttributes .~ attr
| (attr, gen) <- lst ]
where
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]
lst =
[ (attr, balloonP $ shiftLongestDiameter poly)
| (attr, poly) <- toShapes tol svg
]
-- x <= 1
-- diameter (balloonP x p) = diameter p * x
balloonP :: Polygon -> Double -> Polygon
balloonP p = \t ->
let targetLength = d * t
nodeVisible x = ds V.! x <= targetLength
-- Move 'a' closer to 'target' such that the length from point 0 to 'a'
-- is targetLength.
moveCloser :: Int -> V2 Rational -> V2 Rational
moveCloser target a =
let targetDist = ds V.! target
aDist = distance' (pAccess p target) a
frac = min 1 $ realToFrac $ (targetLength - targetDist) / aDist
in lerp frac a (pAccess p target)
worker 0 = [pAccess p 0]
worker a =
let b = pNext p a in
if nodeVisible a && nodeVisible b
then (pAccess p a : pAccess p b : [])
else
chunkRight a b (pAccess p a) (pAccess p b) (fst $ getFunnel a b) ++
chunkCenter a b ++
chunkLeft a b (pAccess p a) (pAccess p b) (snd $ getFunnel a b)
chunkRight ai bi a b (x:y:xs) =
case rayIntersect (a,b) (pAccess p x,pAccess p y) of
Just u ->
if nodeVisible x
then
map (moveCloser x) (split a u) ++
chunkRight ai bi u b (y:xs)
else chunkRight ai bi u b (y:xs)
_ -> -- error $ "chunkRight: urk: " ++ show (ai,bi,x,y)
if nodeVisible x
then map (moveCloser x) [a]
else []
chunkRight _ai _bi _a _b _ = []
chunkLeft ai bi a b (x:y:xs) =
case rayIntersect (a,b) (pAccess p x,pAccess p y) of
Just u ->
if nodeVisible x
then
chunkLeft ai bi a u (y:xs) ++
map (moveCloser x) (split u b)
else chunkLeft ai bi a u (y:xs)
_ -> -- error $ "chunkLeft: urk: " ++ show (ai,bi, x,y)
if nodeVisible x
then map (moveCloser x) [b]
else []
chunkLeft _ai _bi _a _b _ = []
chunkCenter a b =
let (aF, bF) = getFunnel a b
aP = pAccess p a
bP = pAccess p b in
case (reverse aF, reverse bF) of
([x], [_]) | nodeVisible x ->
map (moveCloser x) (split aP bP)
([x], _:left:_) | nodeVisible x ->
case rayIntersect (aP,bP) (pAccess p x,pAccess p left) of
Just v ->
map (moveCloser x) (split aP v)
Nothing -> map (moveCloser x) [aP,bP]
(x:right:_, [_]) | nodeVisible x ->
case rayIntersect (aP,bP) (pAccess p x,pAccess p right) of
Just u -> map (moveCloser x) (split u bP)
Nothing -> map (moveCloser x) [aP,bP] -- error $ "urk: " ++ show (a,b, right)
(x:right:_, _:left:_) | nodeVisible x ->
case rayIntersect (aP,bP) (pAccess p x,pAccess p right) of
Just u ->
case rayIntersect (aP,bP) (pAccess p x,pAccess p left) of
Just v -> map (moveCloser x) (split u v)
Nothing -> map (moveCloser x) [aP,bP]
Nothing -> map (moveCloser x) [aP,bP]
_ -> []
in mkPolygon $ V.fromList $ clearDups $
concatMap worker [0..pSize p-1]
where
clearDups (x:y:xs)
| x == y = clearDups (x:xs)
clearDups (x:xs) = x : clearDups xs
clearDups [] = []
getParents 0 = []
getParents x =
let parent = pParent p 0 x
in parent : getParents parent
getFunnel a b =
let aP = getParents a
bP = getParents b in
(takeUntil (`elem` bP) aP
,takeUntil (`elem` aP) bP)
split aP bP =
let steps = 50 in
[ lerp (t/steps) bP aP
| t <- [0 .. steps]
]
d = V.maximum ds
ds = ssspDistances p
takeUntil :: (a -> Bool) -> [a] -> [a]
takeUntil _fn [] = []
takeUntil fn (x:xs)
| fn x = [x]
| otherwise = x : takeUntil fn xs
diameter :: Polygon -> Double
diameter p = V.maximum (ssspDistances p)
shiftLongestDiameter :: Polygon -> Polygon
shiftLongestDiameter p = findBest 0 p (pCycles p)
where
margin = 0.01
findBest _score elt [] = elt
findBest score elt (x:xs) =
let newScore = diameter x in
if
| newScore-score > score * margin -> findBest newScore x xs
| score-newScore > newScore * margin -> findBest score elt xs
| isTopLeft x elt -> findBest newScore x xs
| otherwise -> findBest score elt xs
isTopLeft a b =
case pAccess a 0-pAccess b 0 of
V2 x y -> y > x
-- Shortest distances from point 0 to all other points.
ssspDistances :: Polygon -> V.Vector Double
ssspDistances p = arr
where
arr = V.generate (pSize p) $ \i ->
case i of
0 -> 0
_ ->
let parent = pParent p 0 i in
arr V.! parent + distance' (pAccess p i) (pAccess p parent)

168
src/Reanimate/Math/Common.hs Normal file
View file

@ -0,0 +1,168 @@
module Reanimate.Math.Common
( -- * Ring
Ring(..)
, ringSize -- :: Ring a -> Int
, ringAccess -- :: Ring a -> Int -> V2 a
, ringClamp -- :: Ring a -> Int -> Int
, ringUnpack -- :: Ring a -> Vector (V2 a)
, ringPack -- :: Vector (V2 a) -> Ring a
, ringMap -- :: (V2 a -> V2 b) -> Ring a -> Ring b
, ringRayIntersect -- :: Ring Rational -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
-- * Math
, area -- :: Fractional a => V2 a -> V2 a -> V2 a -> a
, 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
, isInside -- :: (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
, lineIntersect -- :: (Ord a, Fractional a) => (V2 a, V2 a) -> (V2 a, V2 a) -> Maybe (V2 a)
, distSquared -- :: (Fractional a) => V2 a -> V2 a -> a
, approxDist -- :: (Real a, Fractional a) => V2 a -> V2 a -> a
, distance' -- :: (Real a, Fractional a) => V2 a -> V2 a -> Double
, triangleAngles -- :: V2 Double -> V2 Double -> V2 Double -> (Double, Double, Double)
) where
import Data.Vector (Vector)
import qualified Data.Vector as V
import Linear.Matrix (det33)
import Linear.Metric
import Linear.V2
import Linear.V3
import Linear.Vector
newtype Ring a = Ring (Vector (V2 a))
ringSize :: Ring a -> Int
ringSize (Ring v) = length v
ringAccess :: Ring a -> Int -> V2 a
ringAccess (Ring v) i = v V.! mod i (length v)
ringClamp :: Ring a -> Int -> Int
ringClamp (Ring v) i = mod i (length v)
ringUnpack :: Ring a -> Vector (V2 a)
ringUnpack (Ring v) = v
ringPack :: Vector (V2 a) -> Ring a
ringPack = Ring
ringMap :: (V2 a -> V2 b) -> Ring a -> Ring b
ringMap fn (Ring v) = Ring (V.map fn v)
ringRayIntersect :: Ring Rational -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
ringRayIntersect p (a,b) (c,d) =
rayIntersect (ringAccess p a, ringAccess p b) (ringAccess p c, ringAccess p d)
area :: Fractional a => V2 a -> V2 a -> V2 a -> a
area a b c = 1/2 * area2X a b c
area2X :: Fractional a => V2 a -> V2 a -> V2 a -> a
area2X (V2 a1 a2) (V2 b1 b2) (V2 c1 c2) =
det33 (V3 (V3 a1 a2 1)
(V3 b1 b2 1)
(V3 c1 c2 1))
epsilon :: Fractional a => a
epsilon = 1e-9
epsEq :: (Ord a, Fractional a) => a -> a -> Bool
epsEq a b = abs (a-b) < epsilon
-- Left turn.
isLeftTurn :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
isLeftTurn p1 p2 p3 =
let d = direction p1 p2 p3 in
case compare d 0 of
LT -> True
EQ -> False -- colnear
GT -> False
isLeftTurnOrLinear :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
isLeftTurnOrLinear p1 p2 p3 =
let d = direction p1 p2 p3 in
case compare d 0 of
LT -> True
EQ -> True -- colnear
GT -> False
isRightTurn :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
isRightTurn a b c = not (isLeftTurnOrLinear a b c)
isRightTurnOrLinear :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
isRightTurnOrLinear a b c = not (isLeftTurn a b c)
direction :: Fractional a => V2 a -> V2 a -> V2 a -> a
direction p1 p2 p3 = crossZ (p3-p1) (p2-p1)
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
where
(s, t, _) = barycentricCoords a b c d
barycentricCoords :: Fractional a => V2 a -> V2 a -> V2 a -> V2 a -> (a, a, a)
barycentricCoords (V2 x1 y1) (V2 x2 y2) (V2 x3 y3) (V2 x y) =
(lam1, lam2, lam3)
where
lam1 = ((y2-y3)*(x-x3) + (x3 - x2)*(y-y3)) /
((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
lam2 = ((y3-y1)*(x-x3) + (x1-x3)*(y-y3)) /
((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
lam3 = 1 - lam1 - lam2
{-# INLINE rayIntersect #-}
rayIntersect :: (Fractional a, Ord a) => (V2 a,V2 a) -> (V2 a,V2 a) -> Maybe (V2 a)
rayIntersect (V2 x1 y1,V2 x2 y2) (V2 x3 y3, V2 x4 y4)
| yBot == 0 = Nothing
| otherwise = Just $
V2 (xTop/xBot) (yTop/yBot)
where
xTop = (x1*y2 - y1*x2)*(x3-x4) - (x1 - x2)*(x3*y4-y3*x4)
xBot = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4)
yTop = (x1*y2 - y1*x2)*(y3-y4) - (y1-y2)*(x3*y4-y3*x4)
yBot = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
{-# INLINE isBetween #-}
isBetween :: (Ord a, Fractional a) => V2 a -> (V2 a, V2 a) -> Bool
isBetween (V2 x y) (V2 x1 y1, V2 x2 y2) =
((y1 > y) /= (y2 > y) || y == y1 || y == y2) && -- y is between y1 and y2
((x1 > x) /= (x2 > x) || x == x1 || x == x2)
{-# INLINE lineIntersect #-}
lineIntersect :: (Ord a, Fractional a) => (V2 a, V2 a) -> (V2 a, V2 a) -> Maybe (V2 a)
lineIntersect a b =
case rayIntersect a b of
Just u
| isBetween u a && isBetween u b -> Just u
_ -> Nothing
-- circleIntersect :: (Ord a, Fractional a) => (V2 a, V2 a) -> (V2 a, V2 a) -> [V2 a]
distSquared :: (Fractional a) => V2 a -> V2 a -> a
distSquared a b = quadrance (a ^-^ b)
approxDist :: (Real a, Fractional a) => V2 a -> V2 a -> a
approxDist a b = realToFrac (sqrt (realToFrac (distSquared a b) :: Double))
distance' :: (Real a, Fractional a) => V2 a -> V2 a -> Double
distance' a b = sqrt (realToFrac (distSquared a b))
-- sum of angles is always pi.
triangleAngles :: V2 Double -> V2 Double -> V2 Double -> (Double, Double, Double)
triangleAngles a b c =
(findAngle (b-a) (c-a)
,findAngle (c-b) (a-b)
,findAngle (a-c) (b-c))
where
findAngle v1 v2 = abs (atan2 (crossZ v1 v2) (dot v1 v2))
-- findAngle v1 v2 = acos (dot v1 v2 / (norm v1 * norm v2))

284
src/Reanimate/Math/Compatible.hs Normal file
View file

@ -0,0 +1,284 @@
{-# LANGUAGE MultiWayIf #-}
module Reanimate.Math.Compatible where
import Data.List
import Data.Maybe
import Data.Ord
import qualified Data.Vector as V
import Linear.V2
import Linear.Vector
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Debug.Trace
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_
-- 0..i,j..n-1
-- i..j
split1Link :: Polygon -> Int -> Int -> Int -> (Polygon, Polygon)
split1Link p i j s | j < i = split1Link p j i s
split1Link p i j s =
(mkPolygon $ V.fromList left
,mkPolygon $ V.fromList right)
where
n = pSize p
sp = mkSteinerPoints (pAccess p i) (pAccess p j) s
left = map (pAccess p) [0..i] ++ sp ++ map (pAccess p) [j..n-1]
right = map (pAccess p) [i..j] ++ reverse sp
steiner2Link :: Polygon -> Int -> Int -> V2 Rational
steiner2Link p i j | j < i = steiner2Link p j i
steiner2Link p i j
| isNeighbour = error "steiner2Link: Points are neighbours"
| isParent = error "steiner2Link: Points can directly see each other."
| not (isStraightLine || isGrandparent || oneBendBetween p i j) =
error $ "steiner2Link: Cannot construct 2-link chain between points: " ++ show (i,j,pParent p i j,pParent p i (pParent p i j))
| otherwise =
lerp 0.5 (fst vect) (intersects!!0)
where
distToV = approxDist (fst vect)
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))
-- 0..i,s,j..n-1
-- i..j,s
split2Link :: Polygon -> Int -> Int -> (Polygon, Polygon)
split2Link p i j | j < i = split2Link p j i
split2Link p i j =
(mkPolygon $ V.fromList left
,mkPolygon $ V.fromList right)
where
s = steiner2Link p i j
n = pSize p
left = map (pAccess p) [0..i] ++ [s] ++ map (pAccess p) [j..n-1]
right = map (pAccess p) [i..j] ++ [s]
data Link = OneLink | TwoLink
splitNLink :: Polygon -> Int -> [(Link,Int)] -> (Polygon, Polygon)
splitNLink p i js =
(mkPolygon $ V.fromList left
,mkPolygon $ V.fromList right)
where
n = pSize p
left = map (pAccess p) [0..i] ++ steiners ++ map (pAccess p) [j..n-1]
right = map (pAccess p) [i..j] ++ reverse steiners
j = snd (last js)
steiners = splitNLinks p i js
-- 0, [(TwoLink,2),(OneLink,3)]
-- 0, [(OneLink,5),(TwoLink,3)]
splitNLinks :: Polygon -> Int -> [(Link,Int)] -> [V2 Rational]
splitNLinks _p _i [] = []
splitNLinks p i [(TwoLink,j)] = [steiner2Link p i j]
splitNLinks _p _i [(OneLink,_j)] = []
splitNLinks p i ((TwoLink,j):(OneLink,j'):xs) =
let (l,r) = split2Link p i j
p' = selectContains l r (pAccess p j')
s = steiner2Link p i j
sIdx = fromMaybe (error "missing steiner") $ V.elemIndex s (polygonPoints p')
in s : splitNLinks p' sIdx ((TwoLink,j') : xs)
splitNLinks p i ((TwoLink,j):(TwoLink,j'):xs) =
let (l,r) = split2Link p i j
p' = selectContains l r (pAccess p j')
s = steiner2Link p i j
sIdx = fromMaybe (error "missing steiner") $ V.elemIndex s (polygonPoints p')
in s : splitNLinks p' sIdx ((OneLink,j):(TwoLink,j') : xs)
splitNLinks p i ((OneLink,j):(TwoLink,j'):xs) =
let (l,r) = split2Link p j j'
p' = selectContains l r (pAccess p i)
p'' = selectContains l r (pAccess p j')
s = steiner2Link p j j'
s' = steiner2Link p' i sIdx
sIdx = fromMaybe (error "missing steiner sIdx") $ V.elemIndex s (polygonPoints p')
sIdx' = fromMaybe (error "missing steiner sIdx'") $ V.elemIndex s' (polygonPoints p'')
in s' : s : splitNLinks p'' sIdx' ((OneLink,j') : xs)
splitNLinks _p _i _ = error "splitNLinks: invalid input"
selectContains :: Polygon -> Polygon -> V2 Rational -> Polygon
selectContains p1 p2 elt
| V.elem elt (polygonPoints p1) = p1
| V.elem elt (polygonPoints p2) = p2
| otherwise = error "elt not member of either polygons"
-- [(P, Vis)] -> Polygon -> [P]
{-
If there is a 1-link path, connect it.
Recurse into the two new polygon pairs.
If there is a 2-link path, connect it, possibly add bogus steiner point.
Recurse into the two new polygon pairs.
If there is a n-link path, connect it, add bogus steiner points.
Recurse into the two new polygon pairs.
-}
-- a 2 b 1 c
-- ab = 2-link of (a,b)
-- abc = 2-link of (cut-a-ab-b, c)
-- steps: a ab abc c
-- a, [(2,b), (1,c)]
-- [(ab, p1,p2)]
{-
, mkGroup
[ withStrokeColor "green" $mkGroup
[ mkLine (x1,y1) (x2,y2)
, mkLine (x1',y1') (x2',y2') ]
| (a,b) <- mWins
, (i,j) <- oWins
, a == i || a == j || b == i || b == j
, not (sort [a,b] == sort [i,j])
, let V2 x1 y1 = realToFrac <$> a
V2 x2 y2 = realToFrac <$> b
V2 x1' y1' = realToFrac <$> i
V2 x2' y2' = realToFrac <$> j
]
-}
type Points = V.Vector (V2 Rational)
type Edges = [(Int, Int, Int)]
data Mesh = Mesh { meshPoints :: Points, meshEdges :: Edges }
data MeshPair = MeshPair Points Points Edges
-- The points in a RelMesh are:
-- relMeshStatic ++ x where Ax = B
-- data RelMesh = RelMesh
-- { relMeshStatic :: Points
-- , relMeshEdges :: Edges
-- , relMeshA :: Matrix Double
-- , relMeshB :: Matrix Double
-- }
-- data RelMeshPair = RelMeshPair Points Edges (Matrix Double) (Matrix Double) (Matrix Double) (Matrix Double)
--
--
-- Use SSSP to find 1-link connections
-- Compute visibility windows for each node
-- Check windows for ovelap, giving 2-link connections.
-- Use floyd to generate all other connections
-- (a,b,1)
-- (a,b,2)
-- floyd
compatiblyTriangulateP :: Polygon -> Polygon -> [(Polygon, Polygon)]
compatiblyTriangulateP a b
| pSize a /= pSize b = error "polygon size mismatch"
| otherwise = compatiblyTriangulateP' (pSetOffset a 0) (pSetOffset a 0) (pSetOffset b 0)
compatiblyTriangulateP' :: Polygon -> Polygon -> Polygon -> [(Polygon, Polygon)]
compatiblyTriangulateP' aOrigin a b
| trace ("Polygon size: " ++ show (pSize a)) False = undefined
| n == 3 = trace ("Done") $
[(a,b)]
| otherwise =
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
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
Just (nodeL,nodeR) -> trace (show ("one link"::String, toOriginIndex nodeL, toOriginIndex nodeR)) $
let (aL, aR) = split1Link a nodeL nodeR 0
(bL, bR) = split1Link b nodeL nodeR 0
in compatiblyTriangulateP' aOrigin aL bL ++ compatiblyTriangulateP' aOrigin aR bR
where
toOriginIndex idx =
(idx,fromMaybe (-1) (V.elemIndex (pAccess a idx) (polygonPoints aOrigin)) +
polygonOffset aOrigin)
n = pSize a
bestOneLink = listToMaybe (sortOn (Down . nodeDist)
(aOneLink `intersect` bOneLink))
bestTwoLink = listToMaybe (sortOn (Down . nodeDist)
((aOneLink++aTwoLink) `intersect` (bOneLink++bTwoLink)))
aOneLink = polygonOneLinks a
bOneLink = polygonOneLinks b
aTwoLink = polygonTwoLinks a
bTwoLink = polygonTwoLinks b
nodeDist (i,j) = min (j-i) (n-j+i)
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
polygonTwoLinks :: Polygon -> [(Int,Int)]
polygonTwoLinks p =
[ (i, j)
| i <- [0 .. n-1]
, j <- [i+2 .. n-1]
, not (i==0 && j == n-1)
, pParent p i j /= i -- check for 1-link
, let isTwoLink = pParent p i (pParent p i j) == i
isStraightLine = oneBendBetween p i j
-- Points on a straight line should be 2-link even though they are not
-- direct grandparents.
, isTwoLink || isStraightLine
] where n = pSize p
polygonOneLinks :: Polygon -> [(Int,Int)]
polygonOneLinks p =
[ (i, j)
| i <- [0 .. pSize p-1]
, j <- [i+2 .. pSize p-1]
, not (i==0 && j == pSize p-1)
, pParent p i j == i
]

136
src/Reanimate/Math/DCEL.hs Normal file
View file

@ -0,0 +1,136 @@
{-# LANGUAGE RecordWildCards #-}
module Reanimate.Math.DCEL where
{-
import Control.Applicative
import Data.List.Split
import Data.Vector (Vector)
import qualified Data.Vector as V
import Linear.Matrix
import Linear.V2
import Linear.V3
import Text.Printf
import Reanimate.Math.Common
type EdgeIdx = Int
type VertexIdx = Int
type FaceIdx = Int
data DCEL = DCEL
{ dcelFaces :: Vector EdgeIdx
, dcelVertices :: Vector EdgeIdx
, dcelEdges :: Vector Edge
, dcelPoints :: Vector P
} deriving (Show)
data Edge = Edge
{ edgeVertex :: VertexIdx
, edgeFace :: FaceIdx
, edgeTwin :: EdgeIdx
, edgeNext :: EdgeIdx
, edgePrev :: EdgeIdx
} deriving (Show)
getPoint :: Int -> DCEL -> P
getPoint i DCEL{..} = dcelPoints V.! (i `mod` V.length dcelPoints)
dcelSize :: DCEL -> Int
dcelSize DCEL{..} = V.length dcelPoints
isValidDCEL :: DCEL -> Maybe String
isValidDCEL DCEL{..} = foldr (<|>) empty $
[ V.all (< V.length dcelEdges) dcelFaces ~> "Bad face"
, V.all (< V.length dcelEdges) dcelVertices ~> "Bad vertice"
, V.all (< V.length dcelEdges) (V.map edgeNext dcelEdges) ~> "Bad next"
, V.all (< V.length dcelEdges) (V.map edgePrev dcelEdges) ~> "Bad prev"
, V.all (< V.length dcelEdges) (V.map edgeTwin dcelEdges) ~> "Bad twin"
, V.all (< V.length dcelPoints) (V.map edgeVertex dcelEdges) ~> "Bad edge vertex"
, V.all (< V.length dcelFaces) (V.map edgeFace dcelEdges) ~> "Bad edge face"
, V.length dcelVertices == V.length dcelPoints ~> "Vertex/Point mismatch"
]
where
infix 0 ~>
True ~> _ = Nothing
False ~> msg = Just msg
ppDCEL :: DCEL -> String
ppDCEL DCEL{..} =
unlines
[ "Faces: Edge"
, unlines
[ printf "%3d: %3d" face edge
| (face, edge) <- zip [0::Int .. ] (V.toList dcelFaces)
]
, "Vertices: Edge"
, unlines
[ printf "%3d: %3d" vertex edge
| (vertex, edge) <- zip [0::Int ..] (V.toList dcelVertices)
]
, "Edges: Vertex Face Twin Next Prev"
, unlines
[ printf "%3d: %3d %3d %3d %3d %3d"
edge edgeVertex edgeFace edgeTwin edgeNext edgePrev
| (edge, Edge{..}) <- zip [0::Int ..] (V.toList dcelEdges)
]
, "Points:"
, unlines
[ printf "%3d: " (idx*chunkSize) ++
unwords [ printf "<%.1f,%.1f>" x y | V2 x y <- chunk ]
| (idx, chunk) <- zip [0::Int ..] (chunksOf chunkSize (V.toList dcelPoints))
]
]
where chunkSize = 5
emptyDCEL :: DCEL
emptyDCEL = undefined
{-
faces: [0,1] -- 0 is outside, 1 is inside.
vertices: [0..n-1]
edges:
i: vertex=i
face=1
twin=i+n
next=(i+1)%n
prev=(i-1)%n
i: j=i-n
vertex=edge[j].prev.vertex
face=0
twin=j
next=(j-1)%n + i
prev=(j+1)%n +i
points: pts
-}
fromSimplePolygon :: [V2 Double] -> DCEL
fromSimplePolygon [] = error "empty polygon"
fromSimplePolygon pts | not (isCCW $ V.fromList pts) = error "Polygon not counter-clockwise"
fromSimplePolygon pts = DCEL
{ dcelFaces = V.fromList [0,1]
, dcelVertices = V.fromList [0..n-1]
, dcelEdges = V.fromList $
[ Edge
{ edgeVertex = i
, edgeFace = 1
, edgeTwin = i+n
, edgeNext = (i+1) `mod` n
, edgePrev = (i-1) `mod` n }
| i <- [0..n-1] ] ++
[ Edge
{ edgeVertex = (i-1) `mod` n
, edgeFace = 0
, edgeTwin = i
, edgeNext = (i-1) `mod` n + i
, edgePrev = (i-1) `mod` n + i }
| i <- [0..n-1] ]
, dcelPoints = V.fromList pts
}
where
n = length pts
-- addEdge :: VertexIdx -> VertexIdx -> DCEL -> DCEL ->
-- addEdge a b dcel =
-- add edge, vertex=a, face=a
-- add twin
-}

98
src/Reanimate/Math/EarClip.hs Normal file
View file

@ -0,0 +1,98 @@
module Reanimate.Math.EarClip
( earClip
, earClip'
, isEarCorner
) where
import Data.List
import qualified Data.Set as Set
import Reanimate.Math.Common
import Reanimate.Math.Triangulate
-- import Debug.Trace
-- Triangulation by ear clipping. O(n^2)
earClip :: (Fractional a, Ord a) => Ring a -> Triangulation
earClip = last . earClip'
earClip' :: (Fractional a, Ord a) => Ring a -> [Triangulation]
earClip' p = map (edgesToTriangulation $ ringSize p) $ inits $
let ears = Set.fromList [ i
| i <- elts
, isEarCorner p elts (mod (i-1) n) i (mod (i+1) n) ]
in worker ears (mkQueue elts)
where
n = ringSize p
elts = [0 .. n-1]
-- worker :: Set.Set Int -> PolyQueue Int -> [(P,P)]
worker _ears queue | isSimpleQ queue = []
worker ears queue
-- | trace (show (x, Set.member x ears)) False = undefined
| x `Set.member` ears =
let dq = dropQ queue
v0 = prevQ 1 queue
v1 = prevQ 0 queue
v3 = peekQ dq
v4 = peekQ (nextQ dq)
e1 = if isEarCorner p (toList dq) v0 v1 v3
then Set.insert v1 ears
else Set.delete v1 ears
e2 = if isEarCorner p (toList dq) v1 v3 v4
then Set.insert v3 e1
else Set.delete v3 e1
in (v1,v3) : worker e2 dq
| otherwise = worker ears (nextQ queue)
where
x = peekQ queue
data PolyQueue a = PolyQueue a [a] [a] [a]
-- sizeQ :: PolyQueue a -> Int
-- sizeQ (PolyQueue _ a b _) = 1 + length a + length b
mkQueue :: [a] -> PolyQueue a
mkQueue (x:xs) = PolyQueue x xs [] (reverse (x:xs))
mkQueue [] = error "mkQueue: empty"
toList :: PolyQueue a -> [a]
toList (PolyQueue e a b _) = e : a ++ b
isSimpleQ :: PolyQueue a -> Bool
isSimpleQ (PolyQueue _ xs ys _) =
case xs ++ ys of
[_,_] -> True
_ -> False
peekQ :: PolyQueue a -> a
peekQ (PolyQueue e _ _ _) = e
nextQ :: PolyQueue a -> PolyQueue a
nextQ (PolyQueue x [] ys p) =
let (y:xs) = reverse (x:ys)
in PolyQueue y xs [] (x:p)
nextQ (PolyQueue x (y:xs) ys p) = PolyQueue y xs (x:ys) (x:p)
dropQ :: PolyQueue a -> PolyQueue a
dropQ (PolyQueue _ [] ys p) =
let (x:xs) = reverse ys
in PolyQueue x xs [] p
dropQ (PolyQueue _ (x:xs) ys p) = PolyQueue x xs ys p
prevQ :: Int -> PolyQueue a -> a
prevQ nth (PolyQueue _ _ _ p) = p!!nth
-- O(n)
-- Returns true if ac can be cut from polygon. That is, true if 'b' is an ear.
-- isEarCorner polygon a b c = True iff ac can be cut
isEarCorner :: (Fractional a, Ord a) => Ring a -> [Int] -> Int -> Int -> Int -> Bool
isEarCorner p polygon a b c =
isLeftTurn aP bP cP &&
-- If it is a right turn then the line ac will be outside the polygon
and [ not (isInside aP bP cP (ringAccess p k))
| k <- polygon, k /= a && k /= b && k /= c
]
where
aP = ringAccess p a
bP = ringAccess p b
cP = ringAccess p c

47
src/Reanimate/Math/FloydWarshall.hs Normal file
View file

@ -0,0 +1,47 @@
module Reanimate.Math.FloydWarshall where
import Control.Monad (join)
import Data.List
import Data.Maybe
import qualified Data.Map as Map
import Data.Map (Map)
data Shortest b a = Shortest { distance :: a, path :: [b] }
deriving Show
instance (Ord a, Eq b) => Semigroup (Shortest b a) where
a <> b = case distance a `compare` distance b of
GT -> b
LT -> a
EQ -> a { path = path a `union` path b }
floydWarshall :: (Semigroup b, Ord b, Ord p) =>
[p]
-> Map (p, p) (Shortest p b)
-> Map (p, p) (Maybe (Shortest p b))
floydWarshall v dist = foldr innerCycle (Just <$> dist) v
where
innerCycle k dist = (newDist <$> v <*> v) `setTo` dist
where
newDist i j =
((i,j), do a <- join $ Map.lookup (i, k) dist
b <- join $ Map.lookup (k, j) dist
return $ Shortest (distance a <> distance b) (path a))
setTo = Map.unionWith (<>) . Map.fromList
buildPaths :: Ord p => Map (p,p) (Shortest p b) -> Map (p,p) (Shortest [p] b)
buildPaths d = Map.mapWithKey (\pair s -> s { path = buildPath pair}) d
where
buildPath (i,j)
| i == j = [[j]]
| otherwise = do
k <- path $ fromJust $ Map.lookup (i,j) d
p <- buildPath (k, j)
[i : p]
findMinDistances :: (Monoid b, Ord p, Ord b) => [p] -> Map (p,p) b -> Map (p,p) (Shortest [p] b)
findMinDistances v g =
let weights = Map.mapWithKey (\(_,j) w -> Shortest w [j]) g
trivial = Map.fromList [ ((i,i), Shortest mempty []) | i <- v ]
clean d = fromJust <$> Map.filter isJust (d Map.\\ trivial)
in buildPaths $ clean $ floydWarshall v (weights <> trivial)

693
src/Reanimate/Math/Polygon.hs Normal file
View file

@ -0,0 +1,693 @@
module Reanimate.Math.Polygon
( APolygon(..)
, Polygon
, FPolygon
, P
, mkPolygon -- :: (Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a
, mkPolygonFromRing -- :: (Fractional a, Ord a) => Ring a -> APolygon a
, pParent -- :: Polygon -> Int -> Int -> Int
, pSetOffset -- :: APolygon a -> Int -> APolygon a
, pAdjustOffset -- :: APolygon a -> Int -> APolygon a
, pSize -- :: APolygon a -> Int
, pNull -- :: APolygon a -> Bool
, pNext -- :: APolygon a -> Int -> Int
, pPrev -- :: APolygon a -> Int -> Int
, pIsSimple -- :: Polygon -> Bool
, pIsConvex -- :: Polygon -> Bool
, pIsCCW -- :: Polygon -> Bool
, pScale -- :: Rational -> Polygon -> Polygon
, pAtCentroid -- :: Polygon -> Polygon
, pAtCenter -- :: Polygon -> Polygon
, pTranslate -- :: V2 Rational -> Polygon -> Polygon
, pCenter -- :: Polygon -> V2 Rational
, pBoundingBox -- :: Polygon -> (Rational, Rational, Rational, Rational)
, pIsInside -- :: Polygon -> V2 Rational -> Bool
, pAccess -- :: APolygon a -> Int -> V2 a
, pMkWinding -- :: Int -> Polygon
, pDeoverlap -- :: Polygon -> Polygon
, pCycles -- :: Polygon -> [Polygon]
, pCycle -- :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a
, pCentroid -- :: Polygon -> V2 Rational
, pMapEdges -- :: (V2 Rational -> V2 Rational -> a) -> Polygon -> V.Vector a
, pArea -- :: Polygon -> Rational
, pCircumference -- :: (Real a, Fractional a) => APolygon a -> a
, pCircumference' -- :: (Real a, Fractional a) => APolygon a -> Double
, pAddPoints -- :: Int -> Polygon -> Polygon
, pRayIntersect -- :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
-- * Triangulation
, isValidTriangulation -- :: Polygon -> Triangulation -> Bool
, triangulationsToPolygons -- :: Polygon -> Triangulation -> [Polygon]
-- * Single-Source-Shortest-Path
, ssspVisibility -- :: Polygon -> Polygon
, ssspWindows -- :: Polygon -> [(V2 Rational, V2 Rational)]
-- * Duals
, pdualPolygons -- :: Polygon -> PDual -> [Polygon]
-- * Built-in shapes for testing
, triangle -- :: Polygon
, triangle' -- :: [P]
, shape1 -- :: Polygon
, shape2 -- :: Polygon
, shape3 -- :: Polygon
, shape4 -- :: Polygon
, shape5 -- :: Polygon
, shape6 -- :: Polygon
, shape7 -- :: Polygon
, shape8 -- :: Polygon
, shape9 -- :: Polygon
, shape10 -- :: Polygon
, shape11 -- :: Polygon
, shape12 -- :: Polygon
, shape13 -- :: Polygon
, shape14 -- :: Polygon
, shape15 -- :: Polygon
, shape16 -- :: Polygon
, shape17 -- :: Polygon
, shape18 -- :: Polygon
, shape19 -- :: Polygon
, shape20 -- :: Polygon
, shape21 -- :: Polygon
, shape22 -- :: Polygon
, shape23 -- :: Polygon
, concave -- :: Polygon
-- * Internals
, pRing -- :: APolygon a -> Ring a
, pUnsafeMap -- :: (Ring a -> Ring a) -> APolygon a -> APolygon a
, pCopy -- :: Polygon -> Polygon
, pGenerate -- :: [(Double, Double)] -> Polygon
, pUnGenerate -- :: Polygon -> [(Double, Double)]
) where
import Data.Hashable
import Data.List (intersect, tails, sort, maximumBy)
import Data.Maybe
import Data.Ratio
import Data.Serialize
import Data.Vector (Vector)
import qualified Data.Vector as V
import Linear.V2
import Linear.Vector
import Reanimate.Math.Common
import Reanimate.Math.EarClip
import Reanimate.Math.SSSP
import Reanimate.Math.Triangulate
-- import Debug.Trace
-- Generate random polygons, options:
-- 1. put corners around a circle. Vary the radius.
-- 2. close a hilbert curve
type FPolygon = APolygon Double
-- Optimize representation?
-- Polygon = (Vector XNumerator, Vector XDenominator
-- ,Vector YNumerator, Vector YDenominator)
data APolygon a = Polygon
{ polygonPoints :: Vector (V2 a)
, polygonOffset :: Int
, polygonTriangulation :: Triangulation
, polygonSSSP :: Vector SSSP
}
type Polygon = APolygon Rational
type P = V2 Double
instance Show a => Show (APolygon a) where
show = show . V.toList . polygonPoints
instance Hashable a => Hashable (APolygon a) where
hashWithSalt s p = V.foldl' hashWithSalt s (polygonPoints p)
instance (Fractional a, Ord 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 :: (Fractional a, Ord a) => V.Vector (V2 a) -> APolygon a
mkPolygon points = Polygon
{ polygonPoints = points
, polygonOffset = 0
, polygonTriangulation = trig
, polygonSSSP = V.generate n $ \i -> sssp ring (dual i trig)
}
where
n = length points
ring = ringPack points
trig = earClip ring
mkPolygonFromRing :: (Fractional a, Ord a) => Ring a -> APolygon a
mkPolygonFromRing = mkPolygon . ringUnpack
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 p i j =
(sTree V.! mod (j + polygonOffset p) n - polygonOffset p) `mod` n
where
sTree = polygonSSSP p V.! mod (i + polygonOffset p) n
n = pSize p
pCopy :: Polygon -> Polygon
pCopy p = mkPolygon $ V.generate (pSize p) $ pAccess p
pSetOffset :: APolygon a -> Int -> APolygon a
pSetOffset p offset =
p { polygonOffset = offset `mod` pSize p }
pAdjustOffset :: APolygon a -> Int -> APolygon a
pAdjustOffset p offset =
p { polygonOffset = (polygonOffset p + offset) `mod` pSize p }
pSize :: APolygon a -> Int
pSize = length . polygonPoints
pNull :: APolygon a -> Bool
pNull = V.null . polygonPoints
pNext :: APolygon a -> Int -> Int
pNext p i = (i+1) `mod` pSize p
pPrev :: APolygon a -> Int -> Int
pPrev p i = (i-1) `mod` pSize p
-- When is a polygon valid/simple?
-- It is counter-clockwise.
-- No edges intersect.
-- O(n^2)
-- 'checkEdge' takes 90% of the time.
pIsSimple :: Polygon -> Bool
pIsSimple p | pSize p < 3 = False
pIsSimple p = pIsCCW p && noDups && checkEdge 0 2
where
noDups = checkForDups (sort (V.toList (polygonPoints p)))
checkForDups (x:y:xs)
= x /= y && checkForDups (y:xs)
checkForDups _ = True
len = pSize p
-- check i,i+1 against j,j+1
-- j > i+1
checkEdge i j
| j >= len = (i > len-3) || checkEdge (i+1) (i+3)
| otherwise =
case lineIntersect (pAccess p i, pAccess p $ i+1)
(pAccess p j, pAccess p $ j+1) of
Just u | u /= pAccess p i -> False
_nothing -> checkEdge i (j+1)
pScale :: Rational -> Polygon -> Polygon
pScale s = pUnsafeMap (ringMap (^* s))
pAtCentroid :: Polygon -> Polygon
pAtCentroid p = pTranslate (negate c) p
where c = pCentroid p ^/ 2
pAtCenter :: Polygon -> Polygon
pAtCenter p = pTranslate (negate $ pCenter p) p
pTranslate :: V2 Rational -> Polygon -> Polygon
pTranslate v = pUnsafeMap (ringMap (+v))
pCenter :: Polygon -> V2 Rational
pCenter p = V2 (x+w/2) (y+h/2)
where
(x,y,w,h) = pBoundingBox p
-- Returns (min-x, min-y, width, height)
pBoundingBox :: Polygon -> (Rational, Rational, Rational, Rational)
pBoundingBox = \p ->
let V2 x y = pAccess p 0 in
case V.foldl' worker (x, y, 0, 0) (polygonPoints p) of
(xMin, yMin, xMax, yMax) ->
(xMin, yMin, xMax-xMin, yMax-yMin)
where
worker (xMin,yMin,xMax,yMax) (V2 thisX thisY) =
(min xMin thisX, min yMin thisY
,max xMax thisX, max yMax thisY)
-- Place n points on a circle, use one parameter to slide the points back and forth.
-- Use second parameter to move points closer to center circle.
pGenerate :: [(Double, Double)] -> Polygon
pGenerate points
| len < 4 = error "pGenerate: require at least four points"
| otherwise = mkPolygon $ V.fromList
[ V2 (realToFrac $ cos ang * rMod)
(realToFrac $ sin ang * rMod)
| (i,(angMod,rMod)) <- zip [0..] points
, let minAngle = tau / len * i - pi
maxAngle = tau / len * (i+1) - pi
ang = minAngle + (maxAngle-minAngle)*angMod
]
where
tau = 2*pi
len = fromIntegral (length points)
pUnGenerate :: Polygon -> [(Double, Double)]
pUnGenerate p =
[ worker i (fmap realToFrac e)
| (i,e) <- zip [0..] (V.toList $ polygonPoints p) ]
where
len = fromIntegral (pSize p)
worker i (V2 x y) =
let ang = atan2 y x
minAngle = tau / len * i - pi
maxAngle = tau / len * (i+1) - pi
in ((ang-minAngle)/(maxAngle-minAngle), sqrt (x*x+y*y))
tau = 2*pi
-- When is a triangulation valid?
-- Intersection: No internal edges intersect.
-- Completeness: All edge neighbours share a single internal edge.
isValidTriangulation :: Polygon -> Triangulation -> Bool
isValidTriangulation p t = isComplete && intersectionFree
where
o = polygonOffset p
isComplete = all isProper [0 .. pSize p-1]
isProper i =
let j = pNext p i in
length ((pPrev p i : (t V.! i)) `intersect` (pNext p j : t V.! j)) == 1
intersectionFree = and
[ case lineIntersect (pAccess p (a-o), pAccess p (b-o)) (pAccess p (c-o), pAccess p (d-o)) of
Nothing -> True
Just u -> u == pAccess p (a-o) || u == pAccess p (b-o) ||
u == pAccess p (c-o) || u == pAccess p (d-o)
| ((a,b),(c,d)) <- edgePairs ]
edgePairs = [ (e1, e2) | (e1, rest) <- zip edges (drop 1 $ tails edges), e2 <- rest]
edges =
[ (n, i)
| (n, lst) <- zip [0..] (V.toList t)
, i <- lst
, n < i
]
triangulationsToPolygons :: Polygon -> Triangulation -> [Polygon]
triangulationsToPolygons p t =
[ mkPolygon $ V.fromList
[ pAccess p g, pAccess p i, pAccess p j ]
| i <- [0 .. pSize p-1]
, let js = filter (i<) $ t V.! i
, (g, j) <- zip (i-1:js) js
]
pIsInside :: Polygon -> V2 Rational -> Bool
pIsInside p point = or
[ isInside (rawAccess g) (rawAccess i) (rawAccess j) point
| i <- [0 .. pSize p-1]
, let js = filter (i<) $ polygonTriangulation p V.! i
, (g, j) <- zip (i-1:js) js
]
where
rawAccess x = polygonPoints p V.! x
-- reducePolygons :: Int -> [Polygon] -> [Polygon]
-- reducePolygons n ps
-- | length ps <= n = ps
-- | otherwise =
-- let p = findSmallest ps
-- es = edges p
-- e = findSmallest es
-- in reducePolygons n (merge p e : delete p (delete e ps))
-- where
-- findSmallest = minimumBy (comparing area2X)
-- shareEdge p1 p2 =
pAccess :: APolygon a -> Int -> V2 a
pAccess p i = polygonPoints p V.! ((polygonOffset p + i) `mod` pSize p)
triangle :: Polygon
triangle = mkPolygon $ V.fromList [V2 1 1, V2 0 0, V2 2 0]
triangle' :: [P]
triangle' = reverse [V2 1 1, V2 0 0, V2 2 0]
shape1 :: Polygon
shape1 = mkPolygon $ V.fromList
[ V2 0 0, V2 2 0
, V2 2 1, V2 2 2, V2 2 3, V2 2 4, V2 2 5, V2 2 6
, V2 1 1, V2 0 1 ]
shape2 :: Polygon
shape2 = mkPolygon $ V.fromList
[ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 0 (-1), V2 0 (-2)
, V2 3 (-2), V2 3 2, V2 0 2]
shape3 :: Polygon
shape3 = mkPolygon $ V.fromList
[ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 2, V2 0 2]
shape4 :: Polygon
shape4 = mkPolygon $ V.fromList
[ V2 0 0, V2 1 0, V2 1 1, V2 2 1, V2 2 (-1), V2 3 (-1),V2 3 2, V2 0 2]
shape5 :: Polygon
shape5 = pCycles shape4 !! 2
-- square
shape6 :: Polygon
shape6 = mkPolygon $ V.fromList [ V2 0 0, V2 1 0, V2 1 1, V2 0 1 ]
shape7 :: Polygon
shape7 = pScale 6 $ mkPolygon $ V.fromList
[V2 ((-1567171105775771) % 144115188075855872) ((-7758063241391039) % 1152921504606846976)
,V2 ((-2711114907999263) % 18014398509481984) ((-3561889280168807) % 18014398509481984)
,V2 ((-6897139157863177) % 72057594037927936) ((-1632144794297397) % 4503599627370496)
,V2 (5592137945106423 % 36028797018963968) ((-71351641856107) % 281474976710656)
,V2 (2568147525079071 % 4503599627370496) ((-4312925637247687) % 18014398509481984)
,V2 (1291079014395023 % 2251799813685248) (321513444515769 % 2251799813685248)
,V2 (2071709221627247 % 4503599627370496) (4019115966736491 % 9007199254740992)
,V2 ((-1589087869859839) % 144115188075855872) (4904023654354179 % 9007199254740992)
,V2 ((-2328090886101149) % 36028797018963968) (2587887893460759 % 36028797018963968)
,V2 ((-7990199074159871) % 18014398509481984) (1301850651537745 % 4503599627370496)]
shape8 :: Polygon
shape8 = pScale 10 $ pGenerate
[(0.36,0.4),(0.7,1.8e-2),(0.7,0.2),(0.1,0.4),(0.2,0.2),(0.7,0.1),(0.4,8.0e-2)]
shape9 :: Polygon
shape9 = pScale 5 $ pGenerate
[(0.5,0.2),(0.7,0.6),(0.4,0.3),(0.1,0.7),(0.3,1.0e-2),(0.5,0.3),(0.2,0.8),(0.1,0.8),(0.7,6.0e-2),(0.1,0.6)]
shape10 :: Polygon
shape10 = pGenerate
[(0.4,0.7),(0.2,0.2),(0.3,0.9),(5.0e-2,0.1),(0.7,1.0e-2),(0.7,0.9),(0.2,0.1),(0.5,6.0e-2),(0.6,9.0e-2)]
shape11 :: Polygon
shape11 = pGenerate
[(0.1,0.8),(0.7,0.6),(0.7,0.4),(0.3,0.5),(0.8,0.9),(0.8,6.0e-2),(1.0e-2,4.0e-2),(0.8,0.1)]
shape12 :: Polygon
shape12 = mkPolygon $ V.fromList
[ V2 0 0, V2 0.5 1.5, V2 2 2, V2 (-2) 2, V2 (-0.5) 1.5 ]
-- F shape
shape13 :: Polygon
shape13 = pCycles (mkPolygon $ V.reverse (V.fromList
[ V2 0 0, V2 0 2
, V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1
, V2 1 1, V2 1 0.7
, V2 0.3 0.7, V2 0.3 0 ])) !! 7
-- E shape
shape14 :: Polygon
shape14 = pCycles (mkPolygon $ V.reverse $ V.fromList
[ V2 0 0, V2 0 2 -- up
, V2 1 2, V2 1 1.7, V2 0.3 1.7, V2 0.3 1 -- first prong
, V2 1 1, V2 1 0.7, V2 0.3 0.7, V2 0.3 0.3 -- second prong
, V2 1 0.3, V2 1 0 -- last prong
]) !! 9
--
shape15 :: Polygon
shape15 = mkPolygon $ V.fromList
[ V2 0 0, V2 2 0
, V2 2 2, V2 1 2
, V2 1 1, V2 0 1]
shape16 :: Polygon
shape16 = mkPolygon $ V.fromList
[ V2 0 0, V2 2 0
, V2 2 1, V2 1 1
, V2 1 2, V2 0 2]
shape17 :: Polygon
shape17 = mkPolygon $ V.fromList
[ V2 2 0, V2 2 1
, V2 1 1, V2 1 2
, V2 0 2, V2 0 1, V2 0 0 ]
shape18 :: Polygon
shape18 = mkPolygon $ V.fromList
[ V2 2 0, V2 2 1, V2 2 2
, V2 1 2, V2 1 1
, V2 0 1, V2 0 0 ]
shape19 :: Polygon
shape19 = mkPolygon $ V.fromList
[ V2 (-3) (-3), V2 0 (-1)
, V2 3 (-3), V2 1 0
, V2 3 3, V2 0 1
, V2 (-3) 3, V2 (-1) 0 ]
shape20 :: Polygon
shape20 = mkPolygon $ V.fromList
[ V2 (-3) (-3)
, V2 0 (-1)
, V2 3 (-3)
, V2 5 0
, V2 2.5 (-2)
, V2 1 0
, V2 3 3
, V2 0 1
, V2 (-3) 3
, V2 (-1) 0 ]
shape21 :: Polygon
shape21 = mkPolygon $ V.fromList
[V2 0.0 0.0,V2 1.0 0.0,V2 1.0 1.0,V2 2.0 1.0,V2 2.0 (-1.0),V2 3.0 (-1.0)
,V2 3.0 2.0,V2 0.0 2.0]
shape22 :: Polygon
shape22 = pScale 2 $ mkPolygon $ V.fromList
[V2 (-0.17) (-0.08)
,V2 (-0.34) (-0.21)
,V2 0.0 0.0
,V2 (-0.10) 0.60
,V2 (-0.14) 0.19
,V2 (-0.05) 0.03
]
shape23 :: Polygon
shape23 = mkPolygon $ V.fromList
[ V2 0 0, V2 4 0
, V2 4 3, V2 2 3
, V2 2 2, V2 3 2
, V2 3 1, V2 1 1
, V2 1 2, V2 2 2
, V2 2 3, V2 0 3 ]
concave :: Polygon
concave = mkPolygon $
V.fromList [V2 0 0, V2 2 0, V2 2 2, V2 1 1, V2 0 2]
pMkWinding :: Int -> Polygon
pMkWinding n | n < 1 = error "Polygon must have at least one winding."
pMkWinding n = mkPolygon $
V.fromList $ p0 : p1 : walkTo p1 1 n (V2 1 0) ++ reverse (walkTo p0 1 (n+2) (V2 (-1) 0))
where
p0 = V2 0 0
p1 = V2 0 1
walkTo at a b dir
| a == b = []
| otherwise =
let newAt = at + (dir ^* toRational a)
in newAt : walkTo newAt (a+1) b (rot dir)
rot (V2 x y) =
V2 y (-x)
pDeoverlap :: Polygon -> Polygon
pDeoverlap p = mkPolygon arr
where
arr = V.generate (pSize p) worker
worker 0 = pAccess p 0
worker n =
if length (V.elemIndices (pAccess p n) (polygonPoints p)) /= 1
then
let prev = arr V.! (n-1)
this = pAccess p n
in lerp 0.99999 this prev
else pAccess p n
pCycles :: Polygon -> [Polygon]
pCycles p = map (pAdjustOffset p) [0 .. pSize p-1]
pCycle :: (Real a, Fractional a, Ord a) => APolygon a -> Double -> APolygon a
pCycle p 0 = p
pCycle p t = mkPolygon $ worker 0 0
where
worker acc i
-- | segment + acc == limit =
-- V.drop i p <>
-- V.take i p
| segment + acc > limit =
V.singleton (lerp (realToFrac $ (segment + acc - limit)/segment) x y) <>
-- V.drop (i+1) (polygonPoints p) <>
V.fromList (map (pAccess p) [i+1..pSize p-1]) <>
V.fromList (map (pAccess p) [0 .. i])
-- V.take (i+1) (polygonPoints p)
| i == pSize p-1 = V.fromList (map (pAccess p) [0 .. pSize p-1])
| otherwise = worker (acc+segment) (i+1)
where
x = pAccess p i
y = pAccess p $ i+1
segment = distance' x y
len = pCircumference' p
limit = t * len
pCentroid :: Polygon -> V2 Rational
pCentroid p = V2 cx cy
where
a = pArea p
cx = recip (6*a) * V.sum (pMapEdges fnX p)
cy = recip (6*a) * V.sum (pMapEdges fnY p)
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)
pArea :: Polygon -> Rational
pArea p =
0.5 * V.sum (pMapEdges (\(V2 x y) (V2 x' y') -> x*y' - x'*y) p)
pCircumference :: (Real a, Fractional a) => APolygon a -> a
pCircumference p = sum
[ approxDist (pAccess p i) (pAccess p $ i+1)
| i <- [0 .. pSize p-1]]
pCircumference' :: (Real a, Fractional a) => APolygon a -> Double
pCircumference' p = sum
[ distance' (pAccess p i) (pAccess p $ i+1)
| i <- [0 .. pSize p-1]]
-- Add points by splitting the longest lines in half repeatedly.
pAddPoints :: Int -> Polygon -> Polygon
pAddPoints n p | n <= 0 = p
pAddPoints n p = pAddPoints (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 [0 .. pSize p-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
-- worker 0 _ rest = init rest
-- worker i acc (x:y:xs) =
-- let xy = approxDist x y in
-- if acc + xy > limit
-- then x : worker (i-1) 0 (lerp ((limit-acc)/xy) y x : y:xs)
-- else x : worker i (acc+xy) (y:xs)
-- worker _ _ [_] = []
-- worker _ _ _ = error "addPoints: invalid polygon"
-- s = pSize p
-- len = polygonLength p
-- limit = len / fromIntegral (n+1)
pIsConvex :: Polygon -> Bool
pIsConvex p = and
[ area2X (pAccess p i) (pAccess p j) (pAccess p k) > 0
| i <- [0..n-1]
, j <- [i+1..n-1]
, k <- [j+1..n-1]
]
where n = pSize p
pIsCCW :: Polygon -> Bool
pIsCCW p | pNull p = False
pIsCCW p = V.sum (pMapEdges fn p) < 0
where
fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1)
pRayIntersect :: Polygon -> (Int, Int) -> (Int,Int) -> Maybe (V2 Rational)
pRayIntersect p (a,b) (c,d) =
rayIntersect (pAccess p a, pAccess p b) (pAccess p c, pAccess p d)
---------------------------------------------------------
-- SSSP visibility and SSSP windows
ssspVisibility :: Polygon -> Polygon
ssspVisibility p = mkPolygon $
V.fromList $ clearDups $ go [0 .. pSize p-1] -- ([root..pSize p-1] ++ [0 .. root-1])
where
clearDups (x:y:xs)
| x == y = clearDups (y:xs)
| otherwise = x : clearDups (y:xs)
clearDups xs = xs
obstructedBy n =
case pParent p 0 n of
0 -> n
i -> obstructedBy i
go [] = []
go [x] = [pAccess p x]
go (x:y:xs) =
let xO = obstructedBy x
yO = obstructedBy y
in case () of
()
-- Both ends are visible.
| xO == x && yO == y -> pAccess p x : go (y:xs)
-- X is visible, x to intersect (0,yO) (x,y)
| xO == x ->
pAccess p x : fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y)) : go (y:xs)
-- Y is visible
| yO == y -> fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) : pAccess p y : go (y:xs)
-- Neither is visible and they've obstructed by the same point
-- so the entire edge is hidden.
| xO == yO -> go (y:xs)
-- Neither is visible. Cast shadow from obstruction points to
-- find if a subsection of the edge is visible.
| otherwise ->
let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
in if a /= b
then a : b : go (y:xs)
else go (y:xs)
ssspWindows :: Polygon -> [(V2 Rational, V2 Rational)]
ssspWindows p = clearDups $ go (pAccess p 0) [0..pSize p-1]
where
clearDups (x:y:xs)
| x == y = clearDups (y:xs)
| otherwise = x : clearDups (y:xs)
clearDups xs = xs
obstructedBy n =
case pParent p 0 n of
0 -> n
i -> obstructedBy i
go _ [] = []
go _ [_] = []
go l (x:y:xs) =
let xO = obstructedBy x
yO = obstructedBy y
in case () of
()
-- Both ends are visible.
| xO == x && yO == y -> go (pAccess p x) (y:xs)
-- X is visible, x to intersect (0,yO) (x,y)
| xO == x ->
go (fromMaybe (pAccess p y) (pRayIntersect p (0,yO) (x,y))) (y:xs)
-- Y is visible
| yO == y ->
let newL = fromMaybe (pAccess p x) (pRayIntersect p (0,xO) (x,y)) in
(l, newL) :
go newL (y:xs)
-- Neither is visible and they've obstructed by the same point
-- so the entire edge is hidden.
| xO == yO -> go l (y:xs)
-- Neither is visible. Cast shadow from obstruction points to
-- find if a subsection of the edge is visible.
| otherwise ->
let a = fromMaybe (error "a") (pRayIntersect p (0,xO) (x,y))
b = fromMaybe (error "b") (pRayIntersect p (0,yO) (x,y))
in if a /= b
then (l, a) : (b, pAccess p yO) : go (pAccess p yO) (y:xs)
else go l (y:xs)
pdualPolygons :: Polygon -> PDual -> [Polygon]
pdualPolygons p pdual = map mkPolygonFromRing (pdualRings (pRing p) pdual)

231
src/Reanimate/Math/Render.hs Normal file
View file

@ -0,0 +1,231 @@
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.SSSP
import Reanimate.Math.Triangulate
-- import Reanimate.Math.Visibility
import Reanimate.Scene
import Reanimate.Svg
-- drawVisibility :: Polygon -> Animation
-- drawVisibility p' = mkAnimation 5 $ \t ->
-- let p = cyclePolygon p' (t::Double) in
-- centerUsing (polygonShape p) $
-- mkGroup
-- [ withFillColor "grey" $ polygonShape p
-- , withFillColor "grey" $ polygonDots p
-- , withFillColor "white" $ mkLinePathClosed
-- [ (x,y) | V2 x y <- visibility (map (fmap realToFrac) $ V.toList $ polygonPoints p) ]
-- , let V2 x y = fmap realToFrac $ pAccess p 0 in
-- translate x y $ withFillColor "red" $ mkCircle 0.1
-- -- , withFillColor "blue" $ latex $ T.pack $ show (t)
-- ]
-- drawSSSPVisibility :: Polygon -> Animation
-- drawSSSPVisibility p' = mkAnimation 5 $ \t ->
-- let p = pSetOffset p' (round $ t*(fromIntegral $ pSize p'-1))
-- vis = ssspVisibility p in
-- centerUsing (polygonShape p) $
-- mkGroup
-- [ withFillColor "grey" $ polygonShape p
-- -- , withFillColor "grey" $ polygonDots p
-- -- , withFillColor "white" $ polygonShape vis
-- , let V2 x y = fmap realToFrac $ pAccess p 0 in
-- translate x y $ withFillColor "red" $ mkCircle 0.1
-- -- , withFillColor "blue" $ latex $ T.pack $ show (t)
-- ]
-- drawSSSPVisibilityFast :: Polygon -> Animation
-- drawSSSPVisibilityFast p' = mkAnimation 5 $ \t ->
-- let root = min (pSize p-1) $ (floor $ t*(fromIntegral $ pSize p))
-- p = pSetOffset p' root
-- vis = ssspVisibility p in
-- -- centerUsing (polygonShape p) $
-- mkGroup
-- [ withFillColor "grey" $ polygonShape p
-- , withFillColor "white" $ polygonShape vis
-- , withFillColor "grey" $ polygonNumDots p
-- , let V2 x y = fmap realToFrac $ pAccess p 0 in
-- translate x y $ withFillColor "red" $ mkCircle 0.09
-- -- , withFillColor "blue" $ latex $ T.pack $ show (t)
-- ]
-- drawCompatible :: Polygon -> Polygon -> Animation
-- drawCompatible a b = sceneAnimation $ do
-- newSpriteSVG $ translate (-3) 0 $ mkGroup
-- [ withFillColor "grey" $ polygonShape a
-- , withFillColor "grey" $ polygonNumDots a
-- ]
-- newSpriteSVG $ translate (3) 0 $ mkGroup
-- [ withFillColor "grey" $ polygonShape b
-- , withFillColor "grey" $ polygonNumDots b
-- ]
-- let compat = compatiblyTriangulateP a b
-- forM_ compat $ \(l, r) -> do
-- fork $ play $ staticFrame 1 $
-- translate (-3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
-- withFillOpacity 0 $ polygonShape l
-- fork $ play $ staticFrame 1 $
-- translate (3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
-- withFillOpacity 0 $ polygonShape r
-- drawWindow :: Polygon -> SVG
-- drawWindow p =
-- let mWins = ssspWindows p
-- in
-- mkGroup
-- [ mkGroup
-- [ mkLine (x1,y1) (x2,y2)
-- | (a,b) <- mWins
-- , let V2 x1 y1 = realToFrac <$> a
-- V2 x2 y2 = realToFrac <$> b
-- ]
-- ]
-- drawWindowOverlap :: Polygon -> Int -> Int -> SVG
-- drawWindowOverlap p a b =
-- let aWins = ssspWindows (modOffset p a)
-- bWins = ssspWindows (modOffset p b)
-- in
-- mkGroup
-- [ mkGroup $
-- [ withStrokeColor "green" $mkGroup
-- [ mkLine (x1,y1) (x2,y2)
-- , mkLine (x1',y1') (x2',y2') ]
-- | (a,b) <- aWins
-- , (i,j) <- bWins
-- , a == i || a == j || b == i || b == j
-- , not (sort [a,b] == sort [i,j])
-- , let V2 x1 y1 = realToFrac <$> a
-- V2 x2 y2 = realToFrac <$> b
-- V2 x1' y1' = realToFrac <$> i
-- V2 x2' y2' = realToFrac <$> j
-- ]
-- ]
polygonShape :: Polygon -> SVG
polygonShape p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList (polygonPoints p) ]
polygonDots :: Polygon -> SVG
polygonDots p = mkGroup
[ translate x y $ mkCircle 0.1 | V2 x y <- V.toList $ V.map (fmap realToFrac) $ polygonPoints p ]
polygonNumDots :: Polygon -> SVG
polygonNumDots p = mkGroup $ reverse
[ mkGroup
[ colored n $
translate x y $ mkCircle circR
, withFillColor "black" $
translate x y $ ppNum n ]
| n <- [0..pSize p-1]
, let V2 x y = realToFrac <$> pAccess p n ]
where
circR = 0.05
colored n =
let c = promotePixel $ turbo (fromIntegral (n+2) / fromIntegral (pSize p-1+2))
in withStrokeColorPixel c . withFillColorPixel c
ppNum n =
scaleToSize (circR*1.7) (circR*1.5) $
center $ latex $ T.pack $ "\\texttt{" ++ show n ++ "}"
-- drawSSSP :: Polygon -> (Polygon -> SSSP) -> Animation
-- drawSSSP p gen = mkAnimation 5 $ \t -> centerUsing outline $
-- let p' = cyclePolygons p !! (round $ t*(fromIntegral $ pSize p-1)) in
-- mkGroup
-- [ outline
-- , renderSSSP p' (gen p')
-- , withFillColor "grey" $ polygonNumDots $ p'
-- ]
-- where
-- outline =
-- withFillColor "grey" $ mkLinePathClosed
-- [ (x,y) | V2 x y <- map (fmap realToFrac) (V.toList (polygonPoints p) ++ [pAccess p 0]) ]
drawSSSP :: Polygon -> Animation
drawSSSP p = mkAnimation 5 $ \t -> centerUsing (polygonShape p) $
let root = round $ t*fromIntegral (pSize p-1)
sTree = polygonSSSP p V.! root in
mkGroup
[ polygonShape p
, renderSSSP p sTree
, polygonNumDots p
]
-- Inline such that 'trees' is only computed once.
{-# INLINE drawSSSPNaive #-}
drawSSSPNaive :: Polygon -> Animation
drawSSSPNaive p = mkAnimation 5 $ \t -> centerUsing (polygonShape p) $
let root = round $ t*fromIntegral (pSize p-1) in
mkGroup
[ polygonShape p
, renderSSSP p (trees !! root)
, polygonNumDots p
]
where
trees =
[ naive $ pRing $ pCopy $ pSetOffset p i
| i <- [0 .. pSize p-1]]
renderSSSP :: Polygon -> SSSP -> SVG
renderSSSP p s = withFillOpacity 0 $ withStrokeColor "white" $ mkGroup
[ mkLinePath (lineFrom i)
| i <- [0 .. length s-1] ]
where
lineFrom i =
let V2 ax ay = realToFrac <$> pAccess p i
next = (s V.! i)
in (ax,ay) : if next == i then [] else lineFrom next
drawTriangulation :: Polygon -> (Ring Rational -> [Triangulation]) -> Animation
drawTriangulation p gen = sceneAnimation $
forM_ (gen $ pRing p) $ \t -> play $ staticFrame 1 $ renderTriangulation p t
renderTriangulation :: Polygon -> Triangulation -> SVG
renderTriangulation p t = center $ mkGroup
[ withFillColor "grey" $ polygonShape p
, withStrokeColor "white" $ mkGroup $ concat
[ [ mkLine (ax,ay) (bx,by) ]
| i <- [0..pSize p-1]
, y <- t V.! i
, let V2 ax ay = realToFrac <$> rawAccess i
V2 bx by = realToFrac <$> rawAccess y
]
, withFillColor "grey" $ polygonNumDots p
]
where
rawAccess x = polygonPoints p V.! x
renderDual :: Ring Rational -> Dual -> SVG
renderDual ring d = case d of
Dual (a,b,c) l r -> mkGroup
[ withFillColor "blue" $ mkTrig a b c
, worker c a l
, worker b c r
]
where
mkTrig a b c =
let V2 x1 y1 = realToFrac <$> ringAccess ring a
V2 x2 y2 = realToFrac <$> ringAccess ring b
V2 x3 y3 = realToFrac <$> ringAccess ring c
in mkLinePathClosed [ (x1, y1), (x2,y2), (x3,y3) ]
worker _p1 _p2 EmptyDual = mkGroup []
worker p1 p2 (NodeDual x l r) = mkGroup
[ mkTrig p1 p2 x
, worker x p2 l
, worker p1 x r
]

341
src/Reanimate/Math/SSSP.hs Normal file
View file

@ -0,0 +1,341 @@
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE MultiParamTypeClasses #-}
{-# OPTIONS_GHC -fno-warn-orphans #-}
module Reanimate.Math.SSSP
( -- * Single-Source-Shortest-Path
SSSP
, sssp -- :: (Fractional a, Ord a) => Ring a -> Dual -> SSSP
, dual -- :: Int -> Triangulation -> Dual
, Dual(..)
, DualTree(..)
, PDual
, toPDual -- :: Ring Rational -> Dual -> PDual
, pdualRings -- :: Ring Rational -> PDual -> [Ring Rational]
-- * Misc
, dualToTriangulation -- :: Ring Rational -> Dual -> Triangulation
, pdualReduce -- :: Ring Rational -> PDual -> Int -> PDual
, visibilityArray -- :: Ring Rational -> V.Vector [Int]
, naive -- :: Ring Rational -> SSSP
, naive2 -- :: Ring Rational -> SSSP
, drawDual -- :: Dual -> String
) where
import Control.Monad
import Control.Monad.ST
-- import Data.FingerTree (SearchResult (..), (|>))
-- import qualified Data.FingerTree as F
import Data.Foldable
import Data.List
import qualified Data.Map as Map
import Data.Maybe
import Data.Ord
import Data.STRef
import Data.Tree
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import Reanimate.Math.Common
import Reanimate.Math.Triangulate
-- import Debug.Trace
type SSSP = V.Vector Int
-- ssspParent :: Polygon -> SSSP -> Int -> Int
-- ssspParent p sTree x =
-- (sTree V.! ((x - polygonOffset p) `mod` n) + polygonOffset p) `mod` n
-- where
-- n = polygonSize p
visibilityArray :: Ring Rational -> V.Vector [Int]
visibilityArray p = arr
where
n = ringSize p
arr = V.fromList
[ visibility y
| y <- [0..n-1]
]
visibility y =
[ i
| i <- [0..y-1]
, y `elem` arr V.! i ] ++
[ i
| i <- [y+1 .. n-1]
, let pI = ringAccess p i
isOpen = isRightTurn pYp pY pYn
, ringClamp p (y+1) == i || ringClamp p (y-1) == i || if isOpen
then isLeftTurnOrLinear pY pYn pI ||
isLeftTurnOrLinear pYp pY pI
else not $ isRightTurn pY pYn pI ||
isRightTurn pYp pY pI
, let myEdges = [(e1,e2) | (e1,e2) <- edges, e1/=y, e1/=i, e2/=y,e2/=i]
, all (isNothing . lineIntersect (pY,pI))
[ (ringAccess p e1, ringAccess p e2) | (e1,e2) <- myEdges ]]
where
pY = ringAccess p y
pYn = ringAccess p $ y+1
pYp = ringAccess p $ y-1
edges = zip [0..n-1] (tail [0..n-1] ++ [0])
-- Iterative Single Source Shortest Path solver. Quite slow.
naive :: Ring Rational -> SSSP
naive p =
V.fromList $ Map.elems $
Map.map snd $
worker initial
where
initial = Map.singleton 0 (0,0)
visibility = visibilityArray p
worker :: Map.Map Int (Rational, Int) -> Map.Map Int (Rational, Int)
worker m
| m==newM = newM
| otherwise = worker newM
where
ms' = [ Map.fromList
[ case Map.lookup v m of
Nothing -> (v, (distThroughI, i))
Just (otherDist,parent)
| otherDist > distThroughI -> (v, (distThroughI, i))
| otherwise -> (v, (otherDist, parent))
| v <- visibility V.! i
, let distThroughI = dist + approxDist (ringAccess p i) (ringAccess p v) ]
| (i,(dist,_)) <- Map.toList m
]
newM = Map.unionsWith g (m:ms') :: Map.Map Int (Rational,Int)
g a b = if fst a < fst b then a else b
naive2 :: Ring Rational -> SSSP
naive2 p = runST $ do
parents <- MV.replicate (ringSize p) (-1)
costs <- MV.replicate (ringSize p) (-1)
MV.write parents 0 0
MV.write costs 0 0
changedRef <- newSTRef False
let loop i
| i == ringSize p = do
changed <- readSTRef changedRef
when changed $ do
writeSTRef changedRef False
loop 0
| otherwise = do
myCost <- MV.read costs i
unless (myCost < 0) $
forM_ (visibility V.! i) $ \n -> do
-- n is visible from i.
theirCost <- MV.read costs n
let throughCost = myCost + approxDist (ringAccess p i) (ringAccess p n)
when (throughCost < theirCost || theirCost < 0) $ do
MV.write parents n i
MV.write costs n throughCost
writeSTRef changedRef True
loop (i+1)
loop 0
V.unsafeFreeze parents
where
visibility = visibilityArray p
data PDual = PDual (V.Vector Int) Rational [PDual]
deriving (Show)
toPDual :: Ring Rational -> Dual -> PDual
toPDual p d =
case d of
Dual (a,b,c) l r ->
PDual (V.fromList [a,b,c])
(area2X (ringAccess p a) (ringAccess p b) (ringAccess p c))
(catMaybes [ worker c a l, worker b c r])
where
worker _ _ EmptyDual = Nothing
worker a b (NodeDual x l r) = Just $
PDual (V.fromList [a,x,b])
(area2X (ringAccess p a) (ringAccess p x) (ringAccess p b))
(catMaybes [ worker x b l, worker a x r])
pdualSize :: PDual -> Int
pdualSize (PDual _ _ children) = 1 + sum (map pdualSize children)
pdualArea :: PDual -> Rational
pdualArea (PDual _ faceArea _) = faceArea
-- FIXME: 'origin' isn't used. Remove.
pdualReduce :: Ring Rational -> PDual -> Int -> PDual
pdualReduce origin pdual n
| pdualSize pdual <= n = pdual
| otherwise =
let smallest = minimum $ pAreas pdual
in pdualReduce origin (merge smallest pdual) n
where
merge _s (PDual p faceArea []) = PDual p faceArea []
merge s (PDual p faceArea children)
| faceArea == s =
let (PDual p2 area2 children2:xs) = sortBy (comparing pdualArea) children
in PDual (joinP p p2) (faceArea+area2) (children2++xs)
| otherwise =
let (PDual p2 area2 children2:xs) = sortBy (comparing pdualArea) children
in if area2 == s
then PDual (joinP p p2) (faceArea+area2) (children2++xs)
else PDual p faceArea (map (merge s) children)
pAreas (PDual _ faceArea children) = faceArea : concatMap pAreas children
joinP a b = V.fromList (sort (V.toList a ++ V.toList b))
pdualRings :: Ring Rational -> PDual -> [Ring Rational]
pdualRings p (PDual pts _area children) =
ringPack (V.map (ringAccess p) pts) : concatMap (pdualRings p) children
-- Dual of triangulated polygon
data Dual = Dual (Int,Int,Int) -- (a,b,c)
DualTree -- borders ca
DualTree -- borders bc
deriving (Show)
data DualTree
= EmptyDual
| NodeDual Int -- axb triangle, a and b are from parent.
DualTree -- borders xb
DualTree -- borders ax
deriving (Show)
drawDual :: Dual -> String
drawDual d = drawTree $
case d of
Dual (a,b,c) l r -> Node (show (a,b,c)) [worker c a l, worker b c r]
where
worker _a _b EmptyDual = Node "Leaf" []
worker a b (NodeDual x l r) =
Node (show (b,a,x)) [worker x b l, worker a x r]
dualToTriangulation :: Ring Rational -> Dual -> Triangulation
dualToTriangulation p d = edgesToTriangulation (ringSize p) $ filter goodEdge $
case d of
Dual (a,b,c) l r ->
(a,b):(a,c):(b,c):worker c a l ++ worker b c r
where
goodEdge (a,b)
= a /= ringClamp p (b+1) && a /= ringClamp p (b-1)
worker _a _b EmptyDual = []
worker a b (NodeDual x l r) =
(a,x) : (x, b) : worker x b l ++ worker a x r
-- Dual path:
-- (Int,Int,Int) + V.Vector Int + V.Vector LeftOrRight
-- simplifyDual :: DualTree -> DualTree
-- -- simplifyDual (NodeDual x EmptyDual EmptyDual) = NodeLeaf x
-- -- simplifyDual (NodeDual x l EmptyDual) = NodeDualL x l
-- -- simplifyDual (NodeDual x EmptyDual r) = NodeDualR x r
-- simplifyDual d = d
dual :: Int -> Triangulation -> Dual
dual root t =
case hasTriangle of
[] -> error "weird triangulation"
-- [] -> Dual (0,1,V.length t-1) EmptyDual (dualTree t (1, (V.length t-1)) 0)
(x:_) -> Dual (root,rootNext,x) (dualTree t (x,root) rootNext) (dualTree t (rootNext,x) root)
where
rootNext = idx (root+1)
rootPrev = idx (root-1)
rootNNext = idx (root+2)
idx i = i `mod` n
hasTriangle = (rootPrev : t V.! root) `intersect` (rootNNext : t V.! rootNext)
n = V.length t
-- a=6, b=0, e=1
dualTree :: Triangulation -> (Int,Int) -> Int -> DualTree
dualTree t (a,b) e = -- simplifyDual $
case hasTriangle of
[] -> EmptyDual
[(ab)] ->
NodeDual ab
(dualTree t (ab,b) a)
(dualTree t (a,ab) b)
_ -> error $ "Invalid triangulation: " ++ show (a,b,e,hasTriangle)
where
hasTriangle = (prev a : next a : t V.! a) `intersect` (prev b : next b : t V.! b)
\\ [e]
n = V.length t
next x = (x+1) `mod` n
prev x = (x-1) `mod` n
-- data MinMax = MinMax Int Int | MinMaxEmpty deriving (Show)
-- instance Semigroup MinMax where
-- MinMaxEmpty <> b = b
-- a <> MinMaxEmpty = a
-- MinMax a b <> MinMax c d
-- = MinMax (min a c) (max b d)
-- -- = MinMax c b
-- instance Monoid MinMax where
-- mempty = MinMaxEmpty
--
-- instance F.Measured MinMax Int where
-- measure i = MinMax i i
-- dualRoot :: Dual -> Int
-- dualRoot (Dual (a,_,_) _ _) = a
-- O(n*ln n), could be O(n) if I could figure out how to use fingertrees...
sssp :: (Fractional a, Ord a) => Ring a -> Dual -> SSSP
sssp p d = toSSSP $
case d of
Dual (a,b,c) l r ->
(a, a) :
(b, a) :
(c, a) :
worker [c] [b] a r ++
loopLeft a c l
where
toSSSP edges =
(V.fromList . map snd . sortOn fst) edges
loopLeft a outer l =
case l of
EmptyDual -> []
NodeDual x l' r' ->
(x,a) :
worker [x] [outer] a r' ++
loopLeft a x l'
searchFn _checkStep _cusp _x [] = Nothing
searchFn checkStep cusp x (y:ys)
| not (checkStep (ringAccess p cusp) (ringAccess p y) (ringAccess p x))
= Just $ helper [] y ys
| otherwise = Nothing
where
helper acc v [] = (v, [], reverse acc)
helper acc v1 (v2:vs)
| checkStep (ringAccess p v1) (ringAccess p v2) (ringAccess p x) =
(v1, v2:vs, reverse acc)
| otherwise = helper (v1:acc) v2 vs
searchRight = searchFn isLeftTurn
searchLeft = searchFn isRightTurn
-- adj x = x -- ringClamp p (x-dualRoot d)
-- optTrace msg =
-- if False -- dualRoot d == 1 || dualRoot d == 0
-- then trace msg
-- else id
worker _ _ _ EmptyDual = []
worker f1 f2 cusp (NodeDual x l r) =
-- (optTrace ("Funnel: " ++ show
-- (map adj $ toList f1
-- ,adj cusp
-- ,map adj $ toList f2
-- ,adj x
-- , dualRoot d))
-- ) $
case searchLeft cusp x (toList f1) of
Just (v, f1Hi, f1Lo) ->
-- optTrace (" Visble from left: " ++ show (adj x,adj v)) $
(x, v::Int) :
worker f1Hi [x] v l ++
worker (f1Lo ++ [v, x]) f2 cusp r
Nothing ->
case searchRight cusp x (toList f2) of
Just (v, f2Hi, f2Lo) ->
-- optTrace (" Visble from right: " ++ show (adj x,adj v)) $
(x, v::Int) :
worker f1 (f2Lo ++ [v, x]) cusp l ++
worker [x] f2Hi v r
Nothing ->
-- optTrace (" Visble from cusp: " ++ show (adj x,adj cusp)) $
(x, cusp::Int) :
worker f1 [x] cusp l ++
worker [x] f2 cusp r

289
src/Reanimate/Math/Smooth.hs Normal file
View file

@ -0,0 +1,289 @@
{-# LANGUAGE ScopedTypeVariables #-}
module Reanimate.Math.Smooth
( steinerPoints
, renderMesh
, renderAMesh
, smoothMesh
, smoothStep
, angleSmooth
, meshMinAngle
, splitMeshEdges
)
where
import qualified Data.IntSet as ISet
import Control.Monad
import Control.Monad.ST
import Data.STRef
import Data.List
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import Linear.V2
import Linear.Metric
import Linear.Vector
import Reanimate.Animation
import Reanimate.Constants
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Math.Render
import Reanimate.Math.Triangulate
import Reanimate.Morph.Rigid
import Reanimate.Svg
steinerPoints :: Polygon -> [Polygon] -> [V2 Rational]
steinerPoints outerPolygon = concatMap
(V.toList . V.filter isSteiner . polygonPoints)
where isSteiner p = V.notElem p (polygonPoints outerPolygon)
-- compatiblyTriangulateP :: Polygon -> Polygon -> [(Polygon, Polygon)]
renderMesh :: Polygon -> [Polygon] -> SVG
renderMesh outerPolygon innerPolygons = mkGroup
[ mkGroup
[ withFillOpacity 1
$ withStrokeWidth (defaultStrokeWidth * 0)
$ withStrokeColor "black"
$ withFillColor "lightgreen"
$ polygonShape p
| p <- innerPolygons
]
, mkGroup
[ withFillColor "red" $ aroundCenter (scale 0.2) $ drawPoint
(realToFrac <$> p)
| p <- steiners
]
]
where steiners = steinerPoints outerPolygon innerPolygons
renderAMesh :: Mesh -> SVG
renderAMesh m = mkGroup
[ translate (-1.5) 0 $ renderMesh
(mkPolygon outlineA)
[ mkPolygon $ V.fromList $ map (fmap realToFrac)
[pA V.! a, pA V.! b, pA V.! c]
| (a, b, c) <- V.toList (meshTriangles m)
]
, translate (1.5) 0 $ renderMesh
(mkPolygon outlineB)
[ mkPolygon $ V.fromList $ map (fmap realToFrac)
[pB V.! a, pB V.! b, pB V.! c]
| (a, b, c) <- V.toList (meshTriangles m)
]
]
where
pA = meshPointsA m
pB = meshPointsB m
outlineA = V.map (\i -> realToFrac <$> meshPointsA m V.! i) (meshOutline m)
outlineB = V.map (\i -> realToFrac <$> meshPointsB m V.! i) (meshOutline m)
drawPoint :: V2 Double -> SVG
drawPoint (V2 x y) = translate x y $ mkCircle 0.1
-- Plot angles for a mesh.
-- Plot min-angles for meshes as they are iteratively smoothed.
-- Smoothing algorithm:
-- For each steiner point:
-- Find desired position from each vertex, take average, update position.
-- For each edge:
-- Does flipping it increase min-angle in both meshes?
-- If yes, do it.
-- For long edges:
-- Cut in half:
-- Cut the two connected faces in half.
-- all points
-- steiner points
--
{-
data Mesh = Mesh
{ meshPointsA :: Vector P
, meshPointsB :: Vector P
, meshOutline :: Vector Int
-- , meshSteiner :: Vector Int
, meshTriangles :: Vector RelTrig }
-}
smoothMesh :: Mesh -> [Mesh]
smoothMesh m = runST $ do
trigs <- edgesToTriangulationM (length $ meshPointsA m) $ concat
[ [(a, b), (b, c), (c, a)] | (a, b, c) <- V.toList (meshTriangles m) ]
pA <- V.thaw $ V.map (fmap realToFrac) (meshPointsA m)
pB <- V.thaw $ V.map (fmap realToFrac) (meshPointsB m)
replicateM 120 $ do
smoothStep trigs pA steiner
smoothStep trigs pB steiner
newA <- V.freeze pA
newB <- V.freeze pB
return $ m { meshPointsA = newA, meshPointsB = newB }
where
steiner = V.fromList
[ i | i <- [0 .. length (meshPointsA m) - 1], V.notElem i (meshOutline m) ]
smoothStep
:: forall s
. V.MVector s [Int]
-> V.MVector s (V2 Double)
-> V.Vector Int
-> ST s ()
smoothStep triangulation pts steiner =
forM_ [0 .. length steiner - 1] {-length steiner-1-}
$ \s_i -> do
let i = steiner V.! s_i
pt <- MV.read pts i
edges <- V.fromList <$> (sortEdges pt =<< MV.read triangulation i)
let angleBased = angleSmooth pt edges
laplacian = sum edges ^/ (fromIntegral $ length edges) -- laplacian
-- unsafeIOToST $ hPrint stderr (pt, edges, newX)
if isValidLocation pt edges angleBased
then MV.write pts i angleBased
else if isValidLocation pt edges laplacian
then MV.write pts i laplacian
else return ()
where
sortEdges :: V2 Double -> [Int] -> ST s [V2 Double]
sortEdges pt edges = do
edgePoints <- mapM (MV.read pts) edges
return $ sortOn (dir pt) edgePoints
-- 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)))
angleSmooth :: V2 Double -> V.Vector (V2 Double) -> V2 Double
angleSmooth origin js = V.sum (V.generate n nth) ^/ V.sum (V.generate n factor)
where
n = length js
factor i =
let n_self = js V.! i
n_origin = origin - n_self
n_prev = js V.! mod (i - 1) n - n_self
n_next = js V.! mod (i + 1) n - n_self
a1 = acos (dot n_origin n_next / (norm n_origin * norm n_next))
a2 = acos (dot n_origin n_prev / (norm n_origin * norm n_prev))
alpha = a1 + a2
in recip (alpha * alpha)
nth i =
let V2 x y = origin
n_self@(V2 x_0 y_0) = js V.! i
n_origin = origin - n_self
n_prev = js V.! mod (i - 1) n - n_self
n_next = js V.! mod (i + 1) n - n_self
a1 = acos (dot n_origin n_next / (norm n_origin * norm n_next))
a2 = acos (dot n_origin n_prev / (norm n_origin * norm n_prev))
alpha = a1 + a2
b = (a2 - a1) / 2
x' = x_0 + (x - x_0) * cos b - (y - y_0) * sin b
y' = y_0 + (x - x_0) * sin b + (y - y_0) * cos b
in V2 x' y' ^/ (alpha * alpha)
isValidLocation :: V2 Double -> V.Vector (V2 Double) -> V2 Double -> Bool
isValidLocation origin edges newLoc =
or
[ isInside origin a b newLoc
| i <- [0 .. length edges - 1]
, let a = edges V.! i
b = edges V.! mod (i + 1) (length edges)
]
&& V.toList edges
== sortOn (dir newLoc) (V.toList edges)
&& minAngle origin edges
< minAngle newLoc edges
where
dir :: V2 Double -> V2 Double -> Double
dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a)))
minAngle :: V2 Double -> V.Vector (V2 Double) -> Double
minAngle origin edges = minimum $ concat
[ [a1, a2, a3]
| i <- [0 .. length edges - 1]
, let a = edges V.! i
b = edges V.! mod (i + 1) (length edges)
(a1, a2, a3) = triangleAngles origin a b
]
meshMinAngle :: Mesh -> (Double, Double)
meshMinAngle m =
( minimum $ concat
[ [a1, a2, a3]
| (a, b, c) <- V.toList (meshTriangles m)
, let (a1, a2, a3) = triangleAngles (meshPointsA m V.! a)
(meshPointsA m V.! b)
(meshPointsA m V.! c)
]
, minimum $ concat
[ [a1, a2, a3]
| (a, b, c) <- V.toList (meshTriangles m)
, let (a1, a2, a3) = triangleAngles (meshPointsB m V.! a)
(meshPointsB m V.! b)
(meshPointsB m V.! c)
]
)
-- Smooth:
-- Find all edges
-- Sort counter-clockwise from 12 o'clock. (ie direction from south)
-- Compute new position.
-- Update position.
{-
data Mesh = Mesh
{ meshPointsA :: Vector P
, meshPointsB :: Vector P
, meshOutline :: Vector Int
-- , meshSteiner :: Vector Int
, meshTriangles :: Vector RelTrig }
-}
-- edge: a to b
-- opposite points: c and d
-- delete edge (a,b)
-- insert point p
-- insert edge (a,p)
-- insert edge (p,b)
-- insert edge (c,p)
-- insert edge (p,d)
splitMeshEdges :: Double -> Mesh -> Mesh
splitMeshEdges maxLen mesh = runST $ do
t <- MV.replicate (n * 2) []
forM_ (V.toList (meshTriangles mesh)) $ \(a, b, c) -> do
insertEdge t a b
insertEdge t b c
insertEdge t a c
forM_ [0 .. n - 1] $ \i -> MV.modify t (ISet.toList . ISet.fromList) i
triangulation <- V.freeze t
newPoints <- newSTRef []
offset <- newSTRef n
forM_ [0 .. n - 1] $ \i -> forM_ (triangulation V.! i) $ \j ->
when (i < j) $ when (maxEdgeDistanceSquared i j > maxLen * maxLen) $ do
p <- readSTRef offset
writeSTRef offset (p + 1)
let (p1, p2) = splitEdge i j
modifySTRef newPoints ((p1, p2) :)
let [c, d] = (triangulation V.! i) `intersect` (triangulation V.! j)
deleteEdge t i j
insertEdge t i p
insertEdge t j p
insertEdge t c p
insertEdge t d p
return undefined
where
splitEdge i j =
( lerp 0.5 (meshPointsA mesh V.! i) (meshPointsA mesh V.! j)
, lerp 0.5 (meshPointsB mesh V.! i) (meshPointsB mesh V.! j)
)
maxEdgeDistanceSquared i j =
distSquared (meshPointsA mesh V.! i) (meshPointsA mesh V.! j)
`max` distSquared (meshPointsB mesh V.! i) (meshPointsB mesh V.! j)
n = length (meshPointsA mesh)
deleteEdge v a b = do
MV.modify v (delete b) a
MV.modify v (delete a) b
insertEdge v a b = do
MV.modify v (b :) a
MV.modify v (a :) b
-- sortEdges :: V.Vector (V2 Double) -> V2 Double -> V2 Double -> [Int] -> [Int]
-- sortEdges v prev pt = sortOn (\i -> dir prev pt (v V.! i))
-- dir :: V2 Double -> V2 Double -> V2 Double -> Double
-- dir l a b = atan2 (crossZ (l - a) (b - a)) (dot (l - a) (b - a))

54
src/Reanimate/Math/Triangulate.hs Normal file
View file

@ -0,0 +1,54 @@
module Reanimate.Math.Triangulate
( Triangulation
, edgesToTriangulation
, edgesToTriangulationM
, trianglesToTriangulation
, trianglesToTriangulationM
)
where
import Control.Monad
import qualified Data.IntSet as ISet
import qualified Data.Vector as V
import qualified Data.Vector.Mutable as MV
import Control.Monad.ST
-- Max edges: n-2
-- Each edge is represented twice: 2n-4
-- Flat structure:
-- edges :: V.Vector Int -- max length (2n-4)
-- offsets :: V.Vector Int -- length n
-- Combine the two vectors? < n => offsets, >= n => edges?
type Triangulation = V.Vector [Int]
-- FIXME: Move to Common or a Triangulation module
-- O(n)
edgesToTriangulation :: Int -> [(Int, Int)] -> Triangulation
edgesToTriangulation size edges = runST $ do
v <- edgesToTriangulationM size edges
V.unsafeFreeze v
edgesToTriangulationM :: Int -> [(Int, Int)] -> ST s (V.MVector s [Int])
edgesToTriangulationM size edges = do
v <- MV.replicate size []
forM_ edges $ \(e1, e2) -> do
MV.modify v (e1 :) e2
MV.modify v (e2 :) e1
forM_ [0 .. size - 1] $ \i -> MV.modify v (ISet.toList . ISet.fromList) i
return v
trianglesToTriangulation :: Int -> V.Vector (Int, Int, Int) -> Triangulation
trianglesToTriangulation size edges = runST $ do
v <- trianglesToTriangulationM size edges
V.unsafeFreeze v
trianglesToTriangulationM
:: Int -> V.Vector (Int, Int, Int) -> ST s (V.MVector s [Int])
trianglesToTriangulationM size trigs = do
v <- MV.replicate size []
forM_ (V.toList trigs) $ \(a, b, c) -> do
MV.modify v (\x -> b : c : x) a
MV.modify v (\x -> a : c : x) b
MV.modify v (\x -> a : b : x) c
forM_ [0 .. size - 1] $ \i -> MV.modify v (ISet.toList . ISet.fromList) i
return v

118
src/Reanimate/Math/Visibility.hs Normal file
View file

@ -0,0 +1,118 @@
module Reanimate.Math.Visibility where
-- import qualified Data.Set as Set
import Data.Maybe
import Linear.V2
import Reanimate.Math.Common
import Reanimate.Math.Polygon
-- import Debug.Trace
visibility :: [P] -> [P]
visibility (z:v:vs) = reverse $ go z [v,z] vs
visibility _ = undefined
-- visibility (z:v0:v1:rest)
-- | isLeftTurn z v0 v1 = left z [v1,v0] rest
-- | otherwise = scanA
-- Three cases:
-- v is visible: zsv is a left-turn
-- v moves in front of the stack: s'sv is right turn
-- v moves behind stack:
go :: (Ord a, Fractional a) => V2 a -> [V2 a] -> [V2 a] -> [V2 a]
go _z stack [] = stack
go z stack@(s:s':_ss) (v:vs)
-- | isLeftTurn z s v && isRightTurn s' s v = trace ("FF: " ++ show (z,s,s',v)) $ fastForward z stack s (v:vs)
| isLeftTurn z s v = {-trace ("Left: " ++ show (z,s,v)) $ -}go z (v:stack) vs
| isLeftTurn s' s v = {-trace ("Right: " ++ show (s',s,v,vs)) $ -}rightTurn z stack v vs
| otherwise = {-trace ("FF: " ++ show (z,s,s',v)) $ -}fastForward z stack s (v:vs)
go _ _ _ = undefined
{-
z: 2,2
stack: [1,1 2,6]
v: 0,1
-}
rightTurn :: (Ord a, Fractional a) => V2 a -> [V2 a] -> V2 a -> [V2 a] -> [V2 a]
rightTurn z stack' v (v1:vs)
| isRightTurn z v v1 = {-trace ("Double right: " ++ show (v,v1)) $ -}rightTurn z stack v1 vs
| isLeftTurn z v v1 && isRightTurn (head stack') v v1
= {-trace ("Right->Left: " ++ show (v,v1)) $ -}go z (v:stack) (v1:vs)
| otherwise
= {-trace ("Scan: " ++ show (v, stack')) $ -}scanc z stack v (v1:vs)
where
stack@(_s1:_ss) = unwindStack z stack' v (v1:vs)
rightTurn z stack v [] = unwindStack z stack v []
-- scan forwards until edge intersects zv ray
scanc :: (Ord a, Fractional a) => V2 a -> [V2 a] -> V2 a -> [V2 a] -> [V2 a]
scanc z stack v (v1:v2:vs)
| isBetween u (v1,v2) = -- trace ("Found: " ++ show (u, stack)) $
go z (u:stack) (v2:vs)
--rightTurn z stack u (v2:vs)
| otherwise = scanc z stack v (v2:vs)
where
Just u = rayIntersect (z,v) (v1,v2)
scanc _z stack _v _vs = stack
unwindStack :: (Ord a, Fractional a) => V2 a -> [V2 a] -> V2 a -> t -> [V2 a]
unwindStack z (s1:s2:ss) v vs
| isRightTurn z s1 v && isLeftTurn z s2 v = (u:s2:ss)
| otherwise = unwindStack z (s2:ss) v vs
where
Just u = rayIntersect (z,v) (s1,s2)
unwindStack _z stack _v _vs = stack
-- We've moved into shadow. There are three ways out:
-- 1. Find edge that crosses zv from right to left.
-- In this case, push intersection point to the stack and continue as normal.
-- 2. Find edge that crosses zv from left to right and is above v
-- unwind stack
-- 3. Find edge that crosses zv from left to right and is below v
fastForward :: (Ord a, Fractional a) => V2 a -> [V2 a] -> V2 a -> [V2 a] -> [V2 a]
fastForward z stack v (v1:v2:vs)
| isNothing i || not (isBetween u (v1, v2)) = {-trace ("FF past: " ++ show (z,v,v1,v2)) $ -}fastForward z stack v (v2:vs)
| distSquared v u > distSquared z u = {-trace ("FF skip: " ++ show (z,v,v1,v2)) $ -}fastForward z stack v (v2:vs)
| isLeftTurn z u v2 = {-trace ("FF to: " ++ show u) $ -}go z (v2:u:stack) vs
| distSquared z v < distSquared z u = {-trace ("FF unwind: " ++ show (z,v,v1,v2,u)) $ -}unwindStack z stack v2 (vs)
| otherwise = fastForward z stack v (v2:vs)
where
i = rayIntersect (z, v) (v1, v2)
Just u = rayIntersect (z, v) (v1, v2)
fastForward _z stack _v _vs = stack
{-
v2: 2,1
v1: 1,1
v: 1,0
-}
{-
left z [] [] = finish
left z stack@(s1:s2:ss) (v1:vs)
| isLeftTurn z s1 v1
= left z (v1:stack) rest
| isRightTurn z s1 v1 && isRightTurn s2 s1 v1
= scanA
| otherwise
= right
-- We've made a right turn and need to find the edge in the stack that
-- intersects the zv ray.
-- zv intersects the segment ab iff zav is right-turn and zbv is a left-turn.
right z [] v vs = []
right z (s1:s2:ss) v1 (v2:vs)
| isRightTurn z s1 v1 && isLeftTurn z s2 v1
= let u = lineIntersect (z,v) (s1,s2)
in if isRightTurn z v1 v2
then right
else if isLeftTurn z v1 v2 && isRightTurn
-}
-- Joe and Simpson.
-- vispol (z:v0:v1:vs)
-- | isLeftTurn z v0 v1 = left z [v1,v0] vs
-- | otherwise = scana z [v0] (v1:vs)

1
src/Reanimate/Misc.hs
View file

@ -16,7 +16,6 @@ import GHC.IO.Exception
import System.Directory (copyFile, createDirectory, findExecutable,
getTemporaryDirectory, removeFile,
renameFile)
import System.Exit (ExitCode (..))
import System.FilePath ((<.>), (</>))
import System.IO (hClose, hGetContents, hIsEOF, openTempFile)
import System.Process (readProcessWithExitCode,

42
src/Reanimate/Morph/Cache.hs Normal file
View file

@ -0,0 +1,42 @@
module Reanimate.Morph.Cache
( cachePointCorrespondence -- :: Int -> PointCorrespondence -> PointCorrespondence
) where
import Control.Exception
import qualified Data.ByteString as B
import Data.Hashable
import Data.Serialize
import Reanimate.Cache (encodeInt)
import Reanimate.Morph.Common
import System.Directory
import System.FilePath
import System.IO
import System.IO.Temp
import System.IO.Unsafe
-- type PointCorrespondence = Polygon → Polygon → (Polygon, Polygon)
cachePointCorrespondence :: Int -> PointCorrespondence -> PointCorrespondence
cachePointCorrespondence ident fn src dst = unsafePerformIO $ do
root <- getXdgDirectory XdgCache "reanimate"
createDirectoryIfMissing True root
let path = root </> template
hit <- doesFileExist path
if hit
then do
inp <- B.readFile path
case decode inp of
Left{} -> do
removeFile path
gen path
Right out -> return out
else gen path
where
gen path = do
correspondence <- evaluate (fn src dst)
withSystemTempFile template $ \tmp h -> do
hClose h
B.writeFile tmp (encode correspondence)
renameFile tmp path
return correspondence
template = encodeInt key <.> "morph"
key = hashWithSalt ident (src,dst)

188
src/Reanimate/Morph/Common.hs Normal file
View file

@ -0,0 +1,188 @@
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE TupleSections #-}
{-# LANGUAGE UnicodeSyntax #-}
module Reanimate.Morph.Common
( PointCorrespondence
, Trajectory
, ObjectCorrespondence
, Morph(..)
, morph
, splitObjectCorrespondence
, dupObjectCorrespondence
, genesisObjectCorrespondence
, toShapes
, normalizePolygons
, annotatePolygons
, unsafeSVGToPolygon
) where
import Control.Lens
import qualified Data.Vector as V
import Graphics.SvgTree (DrawAttributes, Texture (..),
drawAttributes, fillColor, fillOpacity,
groupOpacity, strokeColor,
strokeOpacity)
import Linear.V2
import Reanimate.Animation
import Reanimate.Interpolate
import Reanimate.Math.EarClip
import Reanimate.Math.Polygon (Polygon, mkPolygon, pAddPoints,
pCentroid, pRing, pSize, pdualPolygons,
polygonPoints)
import Reanimate.Math.SSSP
import Reanimate.PolyShape
import Reanimate.Signal
import Reanimate.Svg
-- Correspondence
-- Trajectory
-- Color interpolation
-- Polygon holes
-- Polygon splitting
-- Graphical polygon? FIXME: Come up with a better name.
type GPolygon = (DrawAttributes, Polygon)
type PointCorrespondence = Polygon → Polygon → (Polygon, Polygon)
type Trajectory = (Polygon, Polygon) → (Double → Polygon)
type ObjectCorrespondence = [GPolygon] → [GPolygon] → [(GPolygon, GPolygon)]
data Morph = Morph
{ morphTolerance :: Double
, morphColorComponents :: ColorComponents
, morphPointCorrespondence :: PointCorrespondence
, morphTrajectory :: Trajectory
, morphObjectCorrespondence :: ObjectCorrespondence
}
morph :: Morph -> SVG -> SVG -> Double -> SVG
morph Morph{..} src dst = \t ->
case t of
-- 0 -> lowerTransformations src
-- 1 -> lowerTransformations dst
_ -> mkGroup
[ render (genPoints t)
& drawAttributes .~ genAttrs t
| (genAttrs, genPoints) <- gens
]
where
render p = mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) $ V.toList $ polygonPoints p ]
srcShapes = toShapes morphTolerance src
dstShapes = toShapes morphTolerance dst
pairs = morphObjectCorrespondence srcShapes dstShapes
gens =
[ (interpolateAttrs morphColorComponents srcAttr dstAttr, morphTrajectory arranged)
| ((srcAttr, srcPoly'), (dstAttr, dstPoly')) <- pairs
, let arranged = morphPointCorrespondence srcPoly' dstPoly'
]
normalizePolygons :: Polygon -> Polygon -> (Polygon, Polygon)
normalizePolygons src dst =
(pAddPoints (max 0 $ dstN-srcN) src
,pAddPoints (max 0 $ srcN-dstN) dst)
where
srcN = pSize src
dstN = pSize dst
interpolateAttrs :: ColorComponents -> DrawAttributes -> DrawAttributes -> Double -> DrawAttributes
interpolateAttrs colorComps src dst t =
src & fillColor .~ (interpColor <$> src^.fillColor <*> dst^.fillColor)
& strokeColor .~ (interpColor <$> src^.strokeColor <*> dst^.strokeColor)
& fillOpacity .~ (interpOpacity <$> src^.fillOpacity <*> dst^.fillOpacity)
& groupOpacity .~ (interpOpacity <$> src^.groupOpacity <*> dst^.groupOpacity)
& strokeOpacity .~ (interpOpacity <$> src^.strokeOpacity <*> dst^.strokeOpacity)
where
interpColor (ColorRef a) (ColorRef b) =
ColorRef $ interpolateRGBA8 colorComps a b t
-- interpolateColor (ColorRef a) FillNone = ColorRef a
interpColor a _ = a
interpOpacity a b = realToFrac (fromToS (realToFrac a) (realToFrac b) t)
genesisObjectCorrespondence :: ObjectCorrespondence
genesisObjectCorrespondence left right =
case (left, right) of
([] , []) -> []
([], (y1,y2):ys) ->
((y1,y2), (y1, emptyFrom y2 y2)) : genesisObjectCorrespondence [] ys
((x1,x2):xs, []) ->
((x1,x2), (x1, emptyFrom x2 x2)) : genesisObjectCorrespondence xs []
(x:xs, y:ys) ->
(x,y) : genesisObjectCorrespondence xs ys
where
emptyFrom a b = mkPolygon $ V.map (const $ pCentroid a) (polygonPoints b)
dupObjectCorrespondence :: ObjectCorrespondence
dupObjectCorrespondence left right =
case (left, right) of
(_, []) -> []
([], _) -> []
([x], [y]) ->
[(x,y)]
([(x1,x2)], yShapes) ->
let x2s = replicate (length yShapes) x2
in dupObjectCorrespondence (map (x1,) x2s) yShapes
(xShapes, [(y1,y2)]) ->
let y2s = replicate (length xShapes) y2
in dupObjectCorrespondence xShapes (map (y1,) y2s)
(x:xs, y:ys) ->
(x, y) : dupObjectCorrespondence xs ys
splitObjectCorrespondence :: ObjectCorrespondence
splitObjectCorrespondence left right =
case (left, right) of
(_, []) -> []
([], _) -> []
([x], [y]) ->
[(x,y)]
([(x1,x2)], yShapes) ->
let x2s = splitPolygon (length yShapes) x2
in splitObjectCorrespondence (map (x1,) x2s) yShapes
(xShapes, [(y1,y2)]) ->
let y2s = splitPolygon (length xShapes) y2
in splitObjectCorrespondence xShapes (map (y1,) y2s)
(x:xs, y:ys) ->
(x,y) : splitObjectCorrespondence xs ys
splitPolygon :: Int -> Polygon -> [Polygon]
splitPolygon n polygon =
let trig = earClip $ pRing polygon
d = dual 0 trig
pd = toPDual (pRing polygon) d
reduced = pdualReduce (pRing polygon) pd n
polygons = pdualPolygons polygon reduced
in polygons
-- joinPairs :: Correspondence -> [(DrawAttributes, PolyShape)] -> [(DrawAttributes, PolyShape)]
-- -> [(DrawAttributes, DrawAttributes, [(RPoint, RPoint)])]
-- joinPairs _ _ [] = []
-- joinPairs _ [] _ = []
-- joinPairs corr [(x1,x2)] [(y1,y2)] =
-- [(x1,y1, corr x2 y2)]
-- joinPairs corr [(x1,x2)] yShapes =
-- let x2s = splitPolyShape 0.001 (length yShapes) x2
-- in joinPairs corr (map (x1,) x2s) yShapes
-- joinPairs corr xShapes [(y1,y2)] =
-- let y2s = reverse $ splitPolyShape 0.001 (length xShapes) y2
-- in joinPairs corr xShapes (map (y1,) y2s)
-- joinPairs corr ((x1,x2):xs) ((y1,y2):ys) =
-- (x1,y1, corr x2 y2) : joinPairs corr xs ys
-- joinPairs _ _ _ = []
-- FIXME: sort by size, smallest to largest
toShapes :: Double -> SVG -> [(DrawAttributes, Polygon)]
toShapes tol src =
[ (attrs, plToPolygon tol shape)
| (_, attrs, glyph) <- svgGlyphs $ lowerTransformations $ pathify src
, shape <- map mergePolyShapeHoles $ plGroupShapes $ svgToPolyShapes glyph
]
unsafeSVGToPolygon :: Double -> SVG -> Polygon
unsafeSVGToPolygon tol src = snd $ head $ toShapes tol src
annotatePolygons :: (Polygon -> SVG) -> SVG -> SVG
annotatePolygons fn svg = mkGroup
[ fn poly & drawAttributes .~ attr
| (attr, poly) <- toShapes 0.001 svg
]

449
src/Reanimate/Morph/LeastWork.hs Normal file
View file

@ -0,0 +1,449 @@
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE OverloadedStrings #-}
{-# LANGUAGE RecordWildCards #-}
module Reanimate.Morph.LeastWork
( StretchCosts(..)
, defaultStretchCosts
, zeroStretchCosts
, stretchWork
, BendCosts(..)
, defaultBendCosts
, bendWork
, bendInfo
, leastWork
, anyLeastWork
, rawLeastWork
, leastWork'
) where
import Control.Monad
import Control.Monad.ST
import Control.Monad.ST.Unsafe
import Data.Array.Base
import Data.Hashable
import Data.Maybe
import qualified Data.Vector as V
import qualified Data.Vector.Unboxed as VU
import Linear.Metric
import Linear.V2
import Linear.Vector
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Morph.Cache
import Reanimate.Morph.Common
import Reanimate.Morph.Linear
import System.IO
-- import Text.Printf
-- import Debug.Trace
data StretchCosts = StretchCosts
{ stretchStiffness :: Double
, stretchCollapsePenalty :: Double
, stretchElasticity :: Double
}
instance Hashable StretchCosts where
hashWithSalt salt StretchCosts{..} =
hashWithSalt salt (stretchStiffness, stretchCollapsePenalty, stretchElasticity)
defaultStretchCosts :: StretchCosts
defaultStretchCosts = StretchCosts
{ stretchStiffness = 1
, stretchCollapsePenalty = 2
, stretchElasticity = 2 }
zeroStretchCosts :: StretchCosts
zeroStretchCosts = defaultStretchCosts{ stretchStiffness = 0 }
-- stretch work paramters:
-- stretch stiffness
-- collapse penalty
-- exponent
stretchWork :: StretchCosts -> Double -> Double -> Double
stretchWork _ 0 0 = 0
stretchWork StretchCosts{..} lenStart lenEnd = checkNaN $
stretchStiffness * (delta**stretchElasticity)/
((1-c)*lenMin + c*lenMax)
where
checkNaN v
| isNaN v = 0
| otherwise = v
c = recip stretchCollapsePenalty
lenMin = min lenStart lenEnd
lenMax = max lenStart lenEnd
delta = lenEnd-lenStart
data BendCosts = BendCosts
{ bendStiffness :: Double
, bendElasticity :: Double
, bendDeviationPenalty :: Double
, bendZeroPenalty :: Double
}
instance Hashable BendCosts where
hashWithSalt salt BendCosts{..} =
hashWithSalt salt
(bendStiffness
,bendElasticity
,bendDeviationPenalty
,bendZeroPenalty)
defaultBendCosts :: BendCosts
defaultBendCosts = BendCosts
{ bendStiffness = 1
, bendElasticity = 2
, bendDeviationPenalty = 10
, bendZeroPenalty = 1000 }
type PType = Double
bendWork :: BendCosts -> V2 PType -> V2 PType -> V2 PType -> V2 PType -> Double
bendWork BendCosts{..} startA endA startB endB =
case bendInfo startA endA startB endB of
(deltaAngle, deviation, throughZero) ->
(bendStiffness * (deltaAngle + bendDeviationPenalty * deviation))**bendElasticity +
if throughZero then bendZeroPenalty else 0
-- deltaX :: Fractional a => V2 a
-- deltaX = V2 (recip 1000000000) 0
-- testBendInfo :: Polygon -> Polygon -> [Int] -> [Int] -> (Double, Double, Bool)
-- testBendInfo a b ~[a1,a2,a3] ~[b1,b2,b3] =
-- let startA = realToFrac <$> (pAccess a a1 - pAccess a a2)
-- endA = realToFrac <$> (pAccess b b1 - pAccess b b2)
-- endB = realToFrac <$> (pAccess b b3 - pAccess b b2)
-- startB = realToFrac <$> (pAccess a a3 - pAccess a a2)
-- in bendInfo startA endA startB endB
-- zeroBendFactor :: Fractional a => a
-- zeroBendFactor = 1
vZero :: (Ord a, Fractional a) => V2 a -> Bool
vZero (V2 x y) = abs x < epsilon && abs y < epsilon
bendInfo :: V2 PType -> V2 PType -> V2 PType -> V2 PType -> (Double, Double, Bool)
-- bendInfo startA endA startB endB
-- | vZero startA, not (vZero endA) = bendInfo (endA ^* zeroBendFactor) endA startB endB
-- | vZero endA, not (vZero startA) = bendInfo startA (startA ^* zeroBendFactor) startB endB
-- | vZero startB, not (vZero endB) = bendInfo startA endA (endB ^* zeroBendFactor) endB
-- | vZero endB, not (vZero startB) = bendInfo startA endA startB (startB ^* zeroBendFactor)
bendInfo startA endA startB endB =
-- trace (show $ map (fmap realToFrac) [startA, endA, startB, endB]) $
-- trace (show (realToFrac <$> q0, realToFrac <$> q1, realToFrac <$> q2)) $
-- (ang4, alpha+beta, throughZero)
(ang4, alpha+beta, throughZero)
where
!x0 = dot startA startB
!y0 = crossZ startA startB
!x1 = (dot endA startB + dot startA endB) / 2
!y1 = (crossZ endA startB + crossZ startA endB) / 2
!x2 = dot endA endB
!y2 = crossZ endA endB
fudge x | vZero x = V2 (-0.1) 0
fudge x = x
q0,q1, q2 :: V2 PType
!q0 = fudge $ V2 x0 y0
!q1 = V2 x1 y1
!q2 = fudge $ V2 x2 y2
curve n = q0^*squared (1-n) + q1^*(2*n*(1-n)) + q2^*squared n
squared a = a*a
ang3 = abs (angleR q0 q2)
ang4 = if isWrap then 2*pi - ang3 else ang3
d0 = crossZ q0 q1
d1 = crossZ q0 q2 / 2
d2 = crossZ q1 q2
roots = filter (>0) $ filter (<1) $ quadraticRoot (d0+d2-2*d1) (2*d1-2*d0) d0
alpha = fromMaybe 0 $ listToMaybe
[ min q0Angle q2Angle - ang
| root <- roots
, let ang = angleR (V2 1 0) (curve root)
, ang < min q0Angle q2Angle ]
beta = fromMaybe 0 $ listToMaybe
[ ang - max q0Angle q2Angle
| root <- roots
, let ang = angleR (V2 1 0) (curve root)
, ang > max q0Angle q2Angle ]
q0Angle = angleR (V2 1 0) q0
q2Angle = angleR (V2 1 0) q2
isWrap = d1*d1 - d0*d2 < 0 && isInside q0 q1 q2 (V2 0 0)
throughZero = quadThroughZero q0 q1 q2
angleR :: V2 Double -> V2 Double -> Double
angleR a b = atan2 (realToFrac $ crossZ a b) (realToFrac $ dot a b)
leastWork :: StretchCosts -> BendCosts -> PointCorrespondence
leastWork stretchCosts bendCosts =
cachePointCorrespondence key (leastWork_ stretchCosts bendCosts)
where
key = hash ("leastWork v0"::String, stretchCosts, bendCosts)
leastWork_ :: StretchCosts -> BendCosts -> PointCorrespondence
leastWork_ stretchCosts bendCosts src dst
| pSize src > pSize dst =
case leastWork stretchCosts bendCosts dst src of
(dst', src') -> (src', dst')
leastWork_ stretchCosts bendCosts src dst =
worker undefined (-1) options
where
worker bestP _bestPScore [] =
--trace ("\nBest score: " ++ show _bestPScore)
bestP
worker bestP bestPScore (x:xs) =
let (src',dst', newScore) = leastWork' stretchCosts bendCosts src x in
if newScore < bestPScore || bestPScore < 0
then worker (src',dst') newScore xs
else worker bestP bestPScore xs
options = pCycles dst
-- options = [pCycles dst !! 79 ]
anyLeastWork :: StretchCosts -> BendCosts -> PointCorrespondence
anyLeastWork stretchCosts bendCosts src dst =
case closestLinearCorrespondence src dst of
(src', dst') -> case leastWork' stretchCosts bendCosts src' dst' of
(src'', dst'', _score) -> (src'', dst'')
rawLeastWork :: StretchCosts -> BendCosts -> PointCorrespondence
rawLeastWork stretchCosts bendCosts src dst =
case leastWork' stretchCosts bendCosts src dst of
(src', dst', _score) -> (src', dst')
-- type PointCorrespondence = Polygon → Polygon → (Polygon, Polygon)
leastWork' :: StretchCosts -> BendCosts -> Polygon -> Polygon -> (Polygon, Polygon, Double)
leastWork' stretchCosts bendCosts src dst =
leastWork'' stretchCosts bendCosts src dst
srcV dstV (mkDistances srcV) (mkDistances dstV)
where
srcV = V.map (fmap realToFrac) $ polygonPoints src
dstV = V.map (fmap realToFrac) $ polygonPoints dst
mkDistances poly = VU.generate (length poly) $ \i ->
distance (V.unsafeIndex poly i) (V.unsafeIndex poly ((i+1) `mod` length poly))
leastWork'' :: StretchCosts -> BendCosts -> Polygon -> Polygon
-> V.Vector (V2 Double) -> V.Vector (V2 Double)
-> VU.Vector Double -> VU.Vector Double
-> (Polygon, Polygon, Double)
leastWork'' stretchCosts bendCosts src dst srcV dstV srcDist dstDist = runST $ do
-- unsafeIOToST $ hPutStrLn stderr $ "Size: " ++ show (nSrc, nDst)
-- unsafeIOToST $ hFlush stderr
work <- asArray $ newArray ((0,0),(nSrc,nDst)) (0::Double)
north <- asArray $ newArray ((0,0),(nSrc,nDst)) (0::Int)
west <- asArray $ newArray ((0,0),(nSrc,nDst)) (0::Int)
let -- calc :: Int -> Int -> ST s ()
calc = \i j -> do
let cond = False -- i==14 && j==26
w0 <- if i==0 then return 0 else do
prevWork <- readArray work (i-1,j)
prevWest <- readArray west (i-1,j)
prevNorth <- readArray north (i-1,j)
let !len = VU.unsafeIndex srcDist (i-1)
-- !len = distance (realToFrac <$> pAccess src (i-1)) (realToFrac <$> pAccess src i)
!sWork = stretchWork stretchCosts len 0
-- !srcPrev = realToFrac <$> pAccess src (i-1-prevWest)
!srcPrev = V.unsafeIndex srcV (i-1-prevWest)
-- !srcMiddle = realToFrac <$> pAccess src (i-1)
!srcMiddle = V.unsafeIndex srcV (i-1)
-- !srcNext = realToFrac <$> pAccess src i
!srcNext = V.unsafeIndex srcV i
-- !dstPrev = realToFrac <$> pAccess dst (j-prevNorth)
!dstPrev = V.unsafeIndex dstV (j-prevNorth)
-- !dstMiddle = realToFrac <$> pAccess dst j
!dstMiddle = V.unsafeIndex dstV j
!bWork =
bendWork bendCosts
(srcPrev-srcMiddle) (dstPrev-dstMiddle)
(srcNext-srcMiddle) (V2 0 0)
-- when (len' /= len) $ do
-- error $ "w0: Bad length: " ++ show (len', len, i, j)
when cond $ do
unsafeIOToST $ hPutStrLn stderr $ "West: " ++ show (i-1-prevWest, i-1, i)
unsafeIOToST $ hPutStrLn stderr $ "West: " ++ show (j-prevNorth, j, j)
unsafeIOToST $ hPutStrLn stderr $ "bWork: " ++ show (bWork)
unsafeIOToST $ hPutStrLn stderr $ "len: " ++ show (len)
unsafeIOToST $ hPutStrLn stderr $ "sWork: " ++ show (sWork)
return (prevWork + sWork + bWork + if prevWest==0&&prevNorth==1 then 10000 else 0)
w1 <- if j==0 then return 0 else do
prevWork <- readArray work (i,j-1)
prevWest <- readArray west (i,j-1)
prevNorth <- readArray north (i,j-1)
let len = VU.unsafeIndex dstDist (j-1)
--len = distance (realToFrac <$> pAccess dst (j-1)) (realToFrac <$> pAccess dst j)
sWork = stretchWork stretchCosts 0 len
-- srcPrev = realToFrac <$> pAccess src (i-prevWest)
-- srcMiddle = realToFrac <$> pAccess src i
-- dstPrev = realToFrac <$> pAccess dst (j-1-prevNorth)
-- dstMiddle = realToFrac <$> pAccess dst (j-1)
-- dstNext = realToFrac <$> pAccess dst j
srcPrev = V.unsafeIndex srcV (i-prevWest)
srcMiddle = V.unsafeIndex srcV i
dstPrev = V.unsafeIndex dstV (j-1-prevNorth)
dstMiddle = V.unsafeIndex dstV (j-1)
dstNext = V.unsafeIndex dstV j
bWork =
bendWork bendCosts
(srcPrev-srcMiddle) (dstPrev-dstMiddle)
(V2 0 0) (dstNext-dstMiddle)
-- when (len' /= len) $ do
-- error $ "w1: Bad length: " ++ show (len', len, j, realToFrac <$> pAccess dst (j-1), realToFrac <$> pAccess dst j)
when (isNaN bWork) $
error $ "bWork NaN: " ++ show (i,j)
when (isNaN sWork) $
error $ "sWork NaN: " ++ show (i,j, len)
when cond $ do
unsafeIOToST $ hPutStrLn stderr $ "North: " ++ show (i-prevWest, i, i)
unsafeIOToST $ hPutStrLn stderr $ "North: " ++ show (j-1-prevNorth, j-1, j)
unsafeIOToST $ hPutStrLn stderr $ "bWork: " ++ show (bWork)
unsafeIOToST $ hPutStrLn stderr $ "len: " ++ show (len)
unsafeIOToST $ hPutStrLn stderr $ "sWork: " ++ show (sWork)
return (prevWork + sWork + bWork + if prevWest==1&&prevNorth==0 then 10000 else 0)
w2 <- if j==0 || i==0 then return 0 else do
prevWork <- readArray work (i-1,j-1)
prevWest <- readArray west (i-1,j-1)
prevNorth <- readArray north (i-1,j-1)
-- let findPrev key = do
-- keyNorth <- readArray north
let -- startLen = distance (realToFrac <$> pAccess src (i-1)) (realToFrac <$> pAccess src i)
-- endLen = distance (realToFrac <$> pAccess dst (j-1)) (realToFrac <$> pAccess dst j)
startLen = VU.unsafeIndex srcDist (i-1)
endLen = VU.unsafeIndex dstDist (j-1)
sWork = stretchWork stretchCosts startLen endLen
-- srcPrev = realToFrac <$> pAccess src (i-1-prevWest)
-- srcMiddle = realToFrac <$> pAccess src (i-1)
-- srcNext = realToFrac <$> pAccess src i
-- dstPrev = realToFrac <$> pAccess dst (j-1-prevNorth)
-- dstMiddle = realToFrac <$> pAccess dst (j-1)
-- dstNext = realToFrac <$> pAccess dst j
srcPrev = V.unsafeIndex srcV (i-1-prevWest)
srcMiddle = V.unsafeIndex srcV (i-1)
srcNext = V.unsafeIndex srcV i
dstPrev = V.unsafeIndex dstV (j-1-prevNorth)
-- dstPrev' = V.unsafeIndex dstV (j-2)
dstMiddle = V.unsafeIndex dstV (j-1)
dstNext = V.unsafeIndex dstV j
bWork =
bendWork bendCosts
(srcPrev-srcMiddle) (dstPrev-dstMiddle)
(srcNext-srcMiddle) (dstNext-dstMiddle)
-- bWork' =
-- bendWork bendCosts
-- (srcPrev-srcMiddle) (dstPrev'-dstMiddle)
-- (srcNext-srcMiddle) (dstNext-dstMiddle)
when cond $ do
unsafeIOToST $ hPutStrLn stderr $ "Diag: " ++ show (i-1-prevWest, i-1, i)
unsafeIOToST $ hPutStrLn stderr $ "Diag: " ++ show (j-1-prevNorth, j-1, j)
unsafeIOToST $ hPutStrLn stderr $ "bWork: " ++ show (bWork)
-- unsafeIOToST $ hPutStrLn stderr $ "bWork': " ++ show (bWork')
return (prevWork + sWork + bWork)
let goNorth = do
when cond $ do
unsafeIOToST $ hPutStrLn stderr "Go north"
unsafeIOToST $ hPrint stderr (w0, w1, w2)
writeArray work (i,j) w1
writeArray north (i,j) 1
goWest = do
when cond $ do
unsafeIOToST $ hPutStrLn stderr "Go west"
unsafeIOToST $ hPrint stderr (w0, w1, w2)
writeArray work (i,j) w0
writeArray west (i,j) 1
goNorthWest = do
when cond $
unsafeIOToST $ hPutStrLn stderr "Go northwest"
writeArray work (i,j) w2
writeArray north (i,j) 1
writeArray west (i,j) 1
-- when (j==2 && i==2) $
-- unsafeIOToST $ hPrint stderr (i,j,w0, w1, w2)
when (isNaN w0 || isNaN w1 || isNaN w2) $
error $ "NaN: " ++ show (w0,w1,w2, i, j)
if | i==0 -> goNorth
| j==0 -> goWest
| w2 <= w0 && w2 <= w1 -> goNorthWest
| w0 <= w1 && w0 <= w2 -> goWest
| w1 <= w0 && w1 <= w2 -> goNorth
| otherwise -> error $ "urk: " ++ show (w0, w1, w2, i, j)
forM_ [1.. max nSrc nDst-1] $ \idx -> do
forM_ [0..idx-1] $ \i -> do
when (idx < nSrc && i < nDst) $
calc idx i
when (i < nSrc && idx < nDst) $
calc i idx
when (idx < nSrc && idx < nDst) $
calc idx idx
-- workMatrix <- unsafeFreeze work
-- unsafeIOToST $ do
-- hPutStr stderr " "
-- forM_ [0..nSrc-1] $ \i -> do
-- hPutStr stderr (printf "%6d " i)
-- hPutStr stderr "\n"
-- forM_ [0..nDst-1] $ \j -> do
-- hPutStr stderr (printf "%2d " j)
-- forM_ [0..nSrc-1] $ \i -> do
-- hPutStr stderr (printf "%6.1f " (workMatrix ! (i,j)))
-- hPutStr stderr "\n"
finalWork <- readArray work (nSrc-1,nDst-1)
let walk = \i j acc -> do
-- w <- readArray work (i,j)
-- when (i==14) $
-- unsafeIOToST $ hPrint stderr (i,j)
let acc' = (i,j):acc
if (i==0 && j==0)
then return acc'
else do
isNorth <- (==1) <$> readArray north (i,j)
isWest <- (==1) <$> readArray west (i,j)
if | isNorth && isWest -> walk (i-1) (j-1) acc'
| isNorth -> walk i (j-1) acc'
| otherwise -> walk (i-1) j acc'
pairs <- walk (nSrc-1) (nDst-1) []
return
( mkPolygon $ V.fromList [ pAccess src idx | idx <- map fst pairs ]
, mkPolygon $ V.fromList [ pAccess dst idx | idx <- map snd pairs ]
, finalWork
)
where
nSrc = pSize src
nDst = pSize dst
asArray :: ST s (STUArray s i e) -> ST s (STUArray s i e)
asArray = id
{-
quadToCubic (QuadBezier a b c) =
CubicBezier a ((1/3)*^(a ^+^ 2*^b)) ((1/3)*^(2*^b ^+^ c)) c
cubicRoot (p3 - 3*p2 + 3*p1 - p0) (3*p2 - 6*p1 + 3*p0) (3*p1 - 3*p0) p0
where (CubicBezier (Point p0 _) (Point p1 _) (Point p2 _) (Point p3 _)) = b
p3 = a
p2 = ((1/3)*^(a ^+^ 2*^b))
p1 = ((1/3)*^(2*^b ^+^ c))
p0 = c
cubicRoot 0 (a - 2b + c) (2b - 2c) c
-}
quadThroughZero :: V2 Double -> V2 Double -> V2 Double -> Bool
quadThroughZero a@(V2 _ q0) b@(V2 _ q1) c@(V2 _ q2) =
any (xPositive . eval) $ filter (\x -> x > 0 && x < 1) extrema
where
xPositive (V2 x _) = x>0
eval n = a^*squared (1-n) + b^*(2*n*(1-n)) + c^*squared n
squared x = x*x
extrema = quadraticRoot (q2-2*q1+q0) (2*q1-2*q0) q0
quadraticRoot :: Double -> Double -> Double -> [Double]
quadraticRoot a b c
| a == 0 && b == 0 = []
| a == 0 = [-c/b]
| otherwise = result
where
d = b*b - 4*a*c
bSign | b < 0 = -1
| otherwise = 1
q = - (b + bSign * sqrt d) / 2
x1 = q/a
x2 = c/q
result | d < 0 = []
| d == 0 = [x1]
| otherwise = [x1, x2]

108
src/Reanimate/Morph/LineBend.hs Normal file
View file

@ -0,0 +1,108 @@
module Reanimate.Morph.LineBend
( lineBend
, lineBendRaw
) where
import qualified Data.Vector as V
import Linear.Matrix (det22)
import Linear.Metric
import Linear.V2
import Linear.Vector
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Morph.Common
lineBend :: Trajectory
lineBend = lineBend' True
lineBendRaw :: Trajectory
lineBendRaw = lineBend' False
-- Note: Returns polygon with n+1 elements.
lineBend' :: Bool -> Trajectory
lineBend' corrected (a,b) = \t ->
let u = 1 - t
-- phi = V.zipWith (\l r -> u*l + t*r) phi_a phi_b
phi = V.zipWith (\l r -> lerpAngle l r t) phi_a phi_b
-- alpha_zero = u * alpha_a_zero + t * alpha_b_zero
alpha_zero = lerpAngle alpha_a_zero alpha_b_zero t
alpha = V.scanl (+) alpha_zero phi
bigE =
V.sum $ V.zipWith (\diff alp -> squared diff * squared (cos alp)) lengths_diff alpha
bigF =
V.sum $ V.zipWith (\diff alp -> squared diff * sin alp * cos alp) lengths_diff alpha
bigG =
V.sum $ V.zipWith (\diff alp -> squared diff * squared (sin alp)) lengths_diff alpha
bigU =
2 * V.sum
(V.zipWith3 (\len_a len_b alp -> (u*len_a + t*len_b) * cos alp) lengths_a lengths_b alpha)
bigV =
2 * V.sum
(V.zipWith3 (\len_a len_b alp -> (u*len_a + t*len_b) * sin alp) lengths_a lengths_b alpha)
lam1 = det22 (V2 (V2 bigU bigF) (V2 bigV bigG)) /
det22 (V2 (V2 bigE bigF) (V2 bigF bigG))
lam2 = det22 (V2 (V2 bigE bigU) (V2 bigF bigV)) /
det22 (V2 (V2 bigE bigF) (V2 bigF bigG))
s_vect = V.fromList
[ -0.5 * squared (lengths_diff V.! n) *
(lam1 * cos (alpha V.! n) +
lam2 * sin (alpha V.! n))
| n <- [0 .. pSize a-1]]
lengths = V.zipWith3 (\l r s -> u*l + t*r + if corrected then s else 0) lengths_a lengths_b s_vect
movingCenter :: V2 Double
movingCenter = lerp t (realToFrac <$> pCentroid b) (realToFrac <$> pCentroid a)
centerAng :: Double
centerAng = lerpAngle centoid_angle_a centoid_angle_b t
centerLen = u*centoid_len_a + t*centoid_len_b
-- V2 x_zero y_zero = lerp t (realToFrac <$> b V.! 0) (realToFrac <$> a V.! 0)
V2 x_zero y_zero = movingCenter + V2 (cos centerAng * centerLen) (sin centerAng * centerLen)
x_modifiers = V.zipWith (*) lengths (V.map cos alpha)
y_modifiers = V.zipWith (*) lengths (V.map sin alpha)
xs = V.scanl (+) x_zero x_modifiers
ys = V.scanl (+) y_zero y_modifiers
in mkPolygon $ V.map (fmap realToFrac) $ V.zipWith V2 xs ys
where
squared x = x*x
centoid_angle_a :: Double
centoid_angle_a = lineAngle (pCentroid a + V2 1 0) (pCentroid a) (pAccess a 0)
centoid_angle_b = lineAngle (pCentroid b + V2 1 0) (pCentroid b) (pAccess b 0)
centoid_len_a :: Double
centoid_len_a = realToFrac $ approxDist (pCentroid a) (pAccess a 0)
centoid_len_b = realToFrac $ approxDist (pCentroid b) (pAccess b 0)
alpha_a_zero = lineAngle (pAccess a 0 + V2 1 0) (pAccess a 0) (pAccess a 1)
alpha_b_zero = lineAngle (pAccess b 0 + V2 1 0) (pAccess b 0) (pAccess b 1)
lengths_a = computeLength a
lengths_b = computeLength b
lengths_diff' = V.zipWith (\l r -> abs (l-r)) lengths_a lengths_b
lengths_tol = 0.0001 * V.maximum lengths_diff'
lengths_diff = V.map (max lengths_tol) lengths_diff'
phi_a = computePhi a
phi_b = computePhi b
computeLength :: Polygon -> V.Vector Double
computeLength poly = V.fromList
[ realToFrac $ approxDist this next
| n <- [0 .. pSize poly-1]
, let this = pAccess poly n
next = pAccess poly (pNext a n)
]
computePhi poly = V.fromList $
[ negate $ angle' (next-this) ((prev-this) ^* (-1))
| n <- [1 .. pSize poly-1]
, let prev = pAccess poly (n-1)
this = pAccess poly n
next = pAccess poly (pNext a n)
]
lerpAngle :: Double -> Double -> (Double -> Double)
lerpAngle fromAng toAng t
| abs (fromAng - (toAng+2*pi)) < abs (fromAng - toAng) = (1-t)*fromAng + t*(toAng+2*pi)
| abs (fromAng - (toAng-2*pi)) < abs (fromAng - toAng) = (1-t)*fromAng + t*(toAng-2*pi)
| otherwise = (1-t)*fromAng + t*toAng
-- Angle from a through b to c. Measured on the right-hand side, from 0 to tau
lineAngle :: V2 Rational -> V2 Rational -> V2 Rational -> Double
lineAngle a b c = angle' (a-b) (c-b)
angle' :: V2 Rational -> V2 Rational -> Double
angle' a b = atan2 (realToFrac $ crossZ a b) (realToFrac $ dot a b)

52
src/Reanimate/Morph/Linear.hs Normal file
View file

@ -0,0 +1,52 @@
module Reanimate.Morph.Linear
( linear, rawLinear
, linearCorrespondence
, closestLinearCorrespondence
, linearTrajectory
) where
import Data.Hashable
import qualified Data.Vector as V
import Linear.Vector
import Reanimate.Interpolate
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Morph.Cache
import Reanimate.Morph.Common
linear :: Morph
linear = rawLinear
{ morphPointCorrespondence =
cachePointCorrespondence (hash ("closest"::String))
closestLinearCorrespondence }
rawLinear :: Morph
rawLinear = Morph
{ morphTolerance = 0.001
, morphColorComponents = labComponents
, morphPointCorrespondence = linearCorrespondence
, morphTrajectory = linearTrajectory
, morphObjectCorrespondence = splitObjectCorrespondence }
linearCorrespondence :: PointCorrespondence
linearCorrespondence = normalizePolygons
closestLinearCorrespondence :: PointCorrespondence
closestLinearCorrespondence src' dst' =
(src, worker dst (score dst) options)
where
(src, dst) = normalizePolygons src' dst'
worker bestP _bestPScore [] = bestP
worker bestP bestPScore (x:xs) =
let newScore = score x in
if newScore < bestPScore
then worker x newScore xs
else worker bestP bestPScore xs
options = pCycles dst
score = V.sum . V.zipWith approxDist (polygonPoints src) . polygonPoints
linearTrajectory :: Trajectory
linearTrajectory (src,dst)
| pSize src == pSize dst = \t -> mkPolygon $
V.zipWith (lerp $ realToFrac t) (polygonPoints dst) (polygonPoints src)
| otherwise = error $ "Invalid lengths: " ++ show (pSize src, pSize dst)

265
src/Reanimate/Morph/Rigid.hs Normal file
View file

@ -0,0 +1,265 @@
{-# LANGUAGE RecordWildCards #-}
{-
The shape of a triangle can be represented by two explicit points (ie vectors)
and a third implicit point at (0,0). A triangle defined by three points can be
converted by subtracting one of the points from the others:
shape (p1,p2,p3) = (p2-p1, p3-p1)
This keeps the shape intact but drops the position of the triangle.
If P and Q are two such triangles, there exists a matrix A that transforms
P into Q: AP = Q
The 'A' matrix can be calculated by multiplying Q and the inverse of P:
A = QP_inv
So, if we have a set of triangles, we can drop their positional information
and compute their shape transformation matrices. But can we do it in reverse?
Can we start with a set of transformation matrices (which don't give us the
positions of the triangles) and compute a set of points that represent triangles
with the right transformation?
In other words, if we have a set of triangles (with absolute position) and
a transformation matrix for each triangle, can we compute new absolute positions
such that AP=Q?
Solving this is not straightforward as moving a single point will change the
transformation of multiple triangles. Fortunately, though, this problem can be
expression as linear equations and solved quickly with LAPACK.
To express the problem as a linear system, we need 4 equations for each
triangle and 2 unknowns for each point.
The system has this shape: Mx=B
We're trying to find 'x' which are the positions of the new points.
'M' is a matrix that transforms points into transformation matrices.
'B' is the target transformation matrices.
With a bit of rewriting, we can reduce 'A = QP_inv' to 4 linear equations
with 6 unknowns. The 6 unknowns are the x and y positions for the three corners
in a triangle.
This is just another way of calculating the transformation matrix 'A' and we
can verify it by looking at the output of:
* computeA src dst
* applyCoeff (coeffOfB src) dst
The output should be the same because they're both computing 'A'.
'applyCoeff' simply multiplies the coefficients by a 6x1 matrix containing
the x and y positions for 'dst'.
So, this means we have roughly the right shape: Mx=B. M is the coefficients
we found for 'src', 'x' is dst, and 'B' is the transformation matrix that
turns 'src' into 'dst'.
We can get the coefficients for each triangle and place them together in a
single large matrix. Then we can put the desired transformation matrices in 'B'
and ask a linear solver to find 'x' which satisfies the equations.
Phew, step one is done. We now have the machinery for turning arbitrary
transformation matrices into proper, connected triangles. Next step is figuring
out which rotational matrices to use. At the two extremes, we have the identity
matrices (representing our starting shape) and the matrices computed from the
target triangles (represetning out target shape). How can these two matrices be
interpolated for a smooth morph?
There are several ways of interpolating matrices. Linear interpolation is
possible. But, it looks prettier if we separate rotation from shear. Then we
can use spherical linear interpolation on the rotation and linear interpolation
on the shear.
-}
module Reanimate.Morph.Rigid where
import Data.Foldable (toList)
import Data.Vector (Vector)
import qualified Data.Vector as V
import Linear.Quaternion
import Linear.Vector
import Linear.V2
import Linear.V3
import qualified Numeric.LinearAlgebra as Matrix
import Numeric.LinearAlgebra.HMatrix (GMatrix, Matrix,
toLists, (!), (><))
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
-- , meshSteiner :: Vector Int
, meshTriangles :: Vector RelTrig }
-- applyA (computeA a b) a = b + some_constant_translation
-- A = Q P_inv
computeA :: Trig -> Trig -> Matrix Double
computeA p q = matQ <> Matrix.inv matP
where
matP = trigToMatrix p
matQ = trigToMatrix q
-- A = UDV = U(VV_t)DV = (UV)(V_tDV) = RS
-- R = UV
-- S = V_tDV
computeRS :: Matrix Double -> (Matrix Double, Matrix Double)
computeRS a = (r, s)
where
(u,d,vt) = Matrix.svd a
v = Matrix.tr vt
r = u <> v
s = vt <> Matrix.diag d <> v
matrixToQuaternion :: Matrix Double -> Quaternion Double
matrixToQuaternion r = q
where
w = 0.5 * sqrt (1 + r!0!0 + r!1!1 + 1)
z = 1/(4*w) * (r!0!1 - r!1!0)
q = Quaternion w (V3 0 0 z)
quaternionToMatrix :: Quaternion Double -> Matrix Double
quaternionToMatrix q = (2><2)
[ 1 - 2*(qj*qj + qk*qk), 2*(qi*qj+qk*qr)
, 2*(qi*qj - qk*qr), 1 - 2*(qi*qi + qk*qk) ]
where
[qr,qi,qj,qk] = toList q
-- A = R((1-t)I + tS)
computeA_RSt :: Matrix Double -> Matrix Double -> Double -> Matrix Double
computeA_RSt r s t = r_t <> (realToFrac (1-t) * Matrix.ident 2 + realToFrac t * s)
where
i = Quaternion 1 0
q = slerp i (matrixToQuaternion r) t
r_t = quaternionToMatrix q
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"
where
matP = trigToMatrix p
coeffOfB :: Trig -> Matrix Double
coeffOfB p = (4><6)
[ -a0-a2, 0, a0, 0, a2, 0
, -a1-a3, 0, a1, 0, a3, 0
, 0, -a0-a2, 0, a0, 0, a2
, 0, -a1-a3, 0, a1, 0, a3 ]
where
[[a0,a1],[a2,a3]] = toLists (Matrix.inv (trigToMatrix p))
applyCoeff :: Matrix Double -> Trig -> Matrix Double
applyCoeff b (V2 x1 y1, V2 x2 y2, V2 x3 y3) = b <> matQ
where
matQ = (6><1) [ x1, y1, x2, y2, x3, y3]
trigToMatrix :: Trig -> Matrix Double
trigToMatrix (p1,p2,p3) = matP
where
V2 p12x p12y = p2-p1
V2 p13x p13y = p3-p1
matP = (2><2)
[p12x, p13x
,p12y, p13y]
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
}
symmetric :: Bool
symmetric = True
prepare :: Mesh -> Prep
prepare Mesh{..} = Prep
{ prepPivot = (aOrigin, bOrigin)
, prepPointsA =
V.map (subtract aOrigin) $
V.take pivotIdx meshPointsA <> V.drop (pivotIdx+1) meshPointsA
, prepPointsB =
V.map (subtract bOrigin) $
V.take pivotIdx meshPointsB <> V.drop (pivotIdx+1) meshPointsB
, prepRS = rsList
, prepRSRev = rsRevList
, prepUToB = Matrix.mkSparse uToB }
where
aOrigin = meshPointsA V.! pivotIdx
bOrigin = meshPointsB V.! pivotIdx
pivotIdx = case V.head meshTriangles of
(a,_,_) -> a
mkAbs p (a,b,c) = (p V.! a,p V.! b,p V.! c)
absATrigs = V.map (mkAbs meshPointsA) meshTriangles
absBTrigs = V.map (mkAbs meshPointsB) meshTriangles
aList = V.zipWith computeA absATrigs absBTrigs
rsList = V.map computeRS aList
revList = V.zipWith computeA absBTrigs absATrigs
rsRevList = V.map computeRS revList
n = length meshTriangles
-- nT = if symmetric then n*2 else n
-- pivotToB = ((nT*4)><2) $ concat
-- [ [ fromMaybe 0 (lookup (x,pivotIdx*2) bigM)
-- , fromMaybe 0 (lookup (x,pivotIdx*2+1) bigM)]
-- | x <- [0..(nT*4)-1] ]
uToB =
[ ((x,y-2),key) | ((x,y),key) <- bigM, y /= pivotIdx*2 && y /= pivotIdx*2+1 ]
bigM =
concat (zipWith worker [0..] (V.toList meshTriangles)) ++
if symmetric
then concat $ zipWith workerRev [n..] (V.toList meshTriangles)
else []
worker i src@(a,b,c) = concat $
let effs = coeffOfB (mkAbs meshPointsA src) in
[ [((i*4+h, e*2), effs!h!(j*2))
,((i*4+h, e*2+1), effs!h!(j*2+1))]
| h <- [0..3]
, (e,j) <- zip [a,b,c] [0..]
]
workerRev i dst@(a,b,c) = concat $
let effs = coeffOfB (mkAbs meshPointsB dst) in
[ [((i*4+h, e*2), effs!h!(j*2))
,((i*4+h, e*2+1), effs!h!(j*2+1))]
| h <- [0..3]
, (e,j) <- zip [a,b,c] [0..]
]
interpolate :: Prep -> Double -> Vector P
interpolate Prep{..} t = V.fromList $
pivot : worker (Matrix.toList solution)
where
-- solution = Matrix.cgSolve False prepUToB b
solution = Matrix.cgx $ last solutions
solutions = Matrix.cgSolve'
False
1e-9
1e-9
1000
prepUToB
b
(Matrix.fromList $ concat [ [x,y] | V2 x y <- V.toList target ])
-- 0
target = if t < 0.5
then prepPointsA
else prepPointsB
worker (x:y:xs) = V2 x y ^+^ pivot : worker xs
worker _ = []
pivot = case prepPivot of
(src, dst) -> lerp t dst src
n = V.length prepRS
b = Matrix.vector $
[ concat (toLists a)!!j
| i <- [0..n-1]
, let (r,s) = prepRS V.! i
, let a = computeA_RSt r s t
, j <- [0..3]
] ++
[ concat (toLists a)!!j
| symmetric
, i <- [0..n-1]
, let (r,s) = prepRSRev V.! i
, let a = computeA_RSt r s (1-t)
, j <- [0..3]
]

56
src/Reanimate/Morph/Rotational.hs Normal file
View file

@ -0,0 +1,56 @@
module Reanimate.Morph.Rotational
( rotationalTrajectory
, polygonOrigin
) where
import qualified Data.Vector as V
import Linear.Vector
import Linear.V2
import Linear.Metric
import Reanimate.Signal
import Reanimate.Morph.Common
import Reanimate.Math.Polygon
type Origin = (Double, Double)
rotationalTrajectory :: Origin -> Trajectory
rotationalTrajectory origin (src,dst) =
\t ->
let thisOrigin = lerp t dstOrigin srcOrigin in
mkPolygon $
V.generate (pSize src) $ \i ->
let len = fromToS (srcLengths V.! i) (dstLengths V.! i) t
ang = lerpAngle (srcAngles V.! i) (dstAngles V.! i) t
in realToFrac <$> (thisOrigin + V2 (cos ang * len) (sin ang * len))
where
srcOrigin = polygonOrigin src origin
dstOrigin = polygonOrigin dst origin
srcLengths :: V.Vector Double
srcLengths = V.map (distance srcOrigin . fmap realToFrac) $ polygonPoints src
dstLengths = V.map (distance dstOrigin . fmap realToFrac) $ polygonPoints dst
srcAngles = V.map (originAngle srcOrigin . fmap realToFrac) $ polygonPoints src
dstAngles = V.map (originAngle dstOrigin . fmap realToFrac) $ polygonPoints dst
originAngle o = lineAngle (o + V2 1 0) o
polygonOrigin :: Polygon -> Origin -> V2 Double
polygonOrigin poly (originX, originY) =
case pBoundingBox poly of
(polyX, polyY, polyWidth, polyHeight) ->
V2 (realToFrac polyX + realToFrac polyWidth * originX)
(realToFrac polyY + realToFrac polyHeight * originY)
lerpAngle :: Double -> Double -> Double -> Double
lerpAngle fromAng toAng t
| abs (fromAng - (toAng+2*pi)) < abs (fromAng - toAng) = (1-t)*fromAng + t*(toAng+2*pi)
| abs (fromAng - (toAng-2*pi)) < abs (fromAng - toAng) = (1-t)*fromAng + t*(toAng-2*pi)
| otherwise = (1-t)*fromAng + t*toAng
-- Angle from a through b to c.
lineAngle :: V2 Double -> V2 Double -> V2 Double -> Double
lineAngle a b c = angle' (a-b) (c-b)
angle' :: V2 Double -> V2 Double -> Double
angle' a b = atan2 (crossZ a b) (dot a b)

53
src/Reanimate/Morph/Triangulate.hs Normal file
View file

@ -0,0 +1,53 @@
module Reanimate.Morph.Triagulate where
-- import Graphics.SvgTree
import Linear.V2
import Linear.Vector
-- import Reanimate
import Numeric.LinearAlgebra
-- import Numeric.LinearAlgebra.Sparse
-- import Data.Sparse.SpMatrix
-- import Data.Sparse.SpVector
type Poly = [V2 Double]
-- renderPolygon :: Poly -> SVG
-- renderPolygon = undefined
triangulate :: Poly -> Poly -> ([Poly],[Poly])
triangulate = undefined
type Point = V2 Double
p1,p2,p3,p4 :: Point
p1 = V2 0 0
p2 = V2 2 0
p3 = V2 1 2
p4 = V2 1 1
-- let (t1,t2,t3) = barycentricCoordinates p1 p2 p3 p4
-- p4 = t1*p1 + t2*p2 + t3*p3
--
--
-- [ t1
-- , t2
-- , t3 ]
-- [ [p1_x, p2_x, p3_x ] [ p4_x
-- , [p1_y, p2_y, p3_y ] , p4_y
-- , [1, 1, 1] ] , 1 ]
-- barycentricCoordinates :: Point -> Point -> Point -> Point -> (Double, Double, Double)
barycentricCoordinates (V2 p1x p1y) (V2 p2x p2y) (V2 p3x p3y) (V2 p4x p4y) =
linearSolve m v
-- (m, v)
where
m = (3><3) [p1x, p2x, p3x, p1y, p2y, p3y, 1, 1, 1 :: Double]
v = (3><1) [ p4x, p4y, 1 :: Double ]
bCoords (V2 x1 y1) (V2 x2 y2) (V2 x3 y3) (V2 x y) =
(lam1, lam2, lam3)
where
lam1 = ((y2-y3)*(x-x3) + (x3 - x2)*(y-y3)) /
((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
lam2 = ((y3-y1)*(x-x3) + (x1-x3)*(y-y3)) /
((y2-y3)*(x1-x3) + (x3-x2)*(y1-y3))
lam3 = 1 - lam1 - lam2

83
src/Reanimate/PolyShape.hs
View file

@ -2,6 +2,7 @@ module Reanimate.PolyShape
( PolyShape(..)
, PolyShapeWithHoles(..)
, svgToPolyShapes -- :: Tree -> [PolyShape]
, svgToPolygons -- :: Double -> Svg -> [Polygon]
, renderPolyShape -- :: PolyShape -> Tree
, renderPolyShapes -- :: [PolyShape] -> Tree
@ -16,6 +17,7 @@ module Reanimate.PolyShape
, isInsideOf -- :: PolyShape -> PolyShape -> Bool
, plFromPolygon -- :: [RPoint] -> PolyShape
, plToPolygon -- :: Double -> PolyShape -> Polygon
, plPolygonify -- :: Double -> PolyShape -> [Point Double]
, plDecompose -- :: [PolyShape] -> [[RPoint]]
, unionPolyShapes -- :: [PolyShape] -> [PolyShape]
@ -26,26 +28,35 @@ module Reanimate.PolyShape
, mergePolyShapeHoles -- :: PolyShapeWithHoles -> PolyShape
, polyShapeTolerance
, plPartial, plPartialGroup, plPartial'
, splitPolyShape -- :: Double -> Int -> PolyShape -> [PolyShape]
, plGroupTouching
) where
import Chiphunk.Low
import Control.Lens ((&), (.~))
import Data.List (minimumBy, nub, partition, sortOn)
import Control.Lens ((&), (.~))
import Data.AdditiveGroup
import Data.List (minimumBy, nub, partition, sortOn)
import Data.Ord
import qualified Data.Vector as V
import Debug.Trace
import Geom2D.CubicBezier (ClosedPath (..), CubicBezier (..), DPoint,
FillRule (..), PathJoin (..), Point (..),
arcLength, arcLengthParam,
bezierIntersection, bezierSubsegment,
closedPathCurves, closest, colinear,
curvesToClosed, evalBezier,
interpolateVector, reorient, splitBezier,
union, vectorDistance)
import Graphics.SvgTree (PathCommand (..), RPoint, Tree (..),
defaultSvg, pathDefinition)
import Geom2D (rotate90L, rotate90R, ($*), (^*))
import Geom2D.CubicBezier (ClosedPath (..), CubicBezier (..),
DPoint, FillRule (..), PathJoin (..),
Point (..), arcLength, arcLengthParam,
bezierIntersection, bezierSubsegment,
closedPathCurves, closest, colinear,
curvesToClosed, evalBezier,
interpolateVector, reorient,
splitBezier, union, vectorDistance)
import Graphics.SvgTree (PathCommand (..), RPoint, Tree (..),
defaultSvg, pathDefinition)
import Linear.V2
import Reanimate.Animation
import Reanimate.Constants
import Reanimate.Math.EarClip
import Reanimate.Math.Polygon (Polygon, mkPolygon, pIsCCW, pRing,
pdualPolygons, polygonPoints)
import Reanimate.Math.SSSP
import Reanimate.Svg
-- | Shape drawn by continuous line. May have overlap, may be convex.
@ -91,6 +102,13 @@ plFromPolygon = PolyShape . ClosedPath . map worker
where
worker (V2 x y) = (Point x y, JoinLine)
plToPolygon :: Double -> PolyShape -> Polygon
plToPolygon tol pl =
let p = V.init . V.fromList . map (\(Point x y) -> realToFrac <$> V2 x y) .
plPolygonify tol $ pl
in if pIsCCW (mkPolygon p) then mkPolygon p else mkPolygon (V.reverse p)
plPartial :: Double -> PolyShape -> PolyShape
plPartial delta pl | delta >= 1 = pl
plPartial delta pl = PolyShape $ curvesToClosed (lineOut ++ [joinB] ++ lineIn)
@ -110,6 +128,26 @@ plPartial delta pl = PolyShape $ curvesToClosed (lineOut ++ [joinB] ++ lineIn)
then [bezierSubsegment c 0 (arcLengthParam c l polyShapeTolerance)]
else c : takeLen (l-cLen) cs
-- earClip :: Polygon -> Triangulation
-- dual :: Triangulation -> Dual
-- toPDual :: Polygon -> Dual -> PDual
-- pdualReduce :: Polygon -> PDual -> Int -> PDual
-- pdualPolygons :: Polygon -> PDual -> [Polygon]
splitPolyShape :: Double -> Int -> PolyShape -> [PolyShape]
splitPolyShape tol n poly =
let polygon = toPolygon (plPolygonify tol poly)
trig = earClip $ pRing polygon
d = dual 0 trig
pd = toPDual (pRing polygon) d
reduced = pdualReduce (pRing polygon) pd n
polygons = pdualPolygons polygon reduced
in map toPolyShape polygons
where
toPolygon :: [Point Double] -> Polygon
toPolygon = mkPolygon . V.fromList . nub . map (\(Point x y) -> V2 (realToFrac x) (realToFrac y))
toPolyShape :: Polygon -> PolyShape
toPolyShape = plFromPolygon . map (fmap realToFrac) . V.toList . polygonPoints
plPartialGroup :: Double -> [PolyShape] -> [PolyShape]
plPartialGroup _delta [] = []
plPartialGroup delta pls =
@ -186,7 +224,7 @@ plPolygonify tol shape =
worker c | endPoint c == startPoint c =
[] -- error $ "Bad bezier: " ++ show c
worker c =
if colinear c tol
if colinear c tol -- && arcLength c 1 tol < 1
then [endPoint c]
else
let (lhs,rhs) = splitBezier c 0.5
@ -220,6 +258,13 @@ plLineCommands pl =
svgToPolyShapes :: Tree -> [PolyShape]
svgToPolyShapes = cmdsToPolyShapes . toLineCommands . extractPath
svgToPolygons :: Double -> SVG -> [Polygon]
svgToPolygons tol = map (toPolygon . plPolygonify tol) . svgToPolyShapes
where
toPolygon :: [Point Double] -> Polygon
toPolygon = mkPolygon .
V.fromList . nub . map (\(Point x y) -> V2 (realToFrac x) (realToFrac y))
cmdsToPolyShapes :: [LineCommand] -> [PolyShape]
cmdsToPolyShapes [] = []
cmdsToPolyShapes cmds =
@ -362,20 +407,26 @@ cutSingleHole parent child =
[x2p]
)
where
vect = (childOrigin ^-^ p) ^* 0.0001
vectL = rotate90L $* vect
vectR = rotate90R $* vect
score = vectorDistance childOrigin p
childOrigin = polyShapeOrigin child
childOrigin' = childOrigin ^-^ vectL
(pHead:pTail) = plCurves parent
childCurves = plCurves child
pParam = closest pHead childOrigin polyShapeTolerance
(a2p, p2b) = splitBezier pHead pParam
(a2p, p2b') = splitBezier pHead pParam
p2b = case p2b' of
CubicBezier a b c d -> CubicBezier (a ^-^ vectL) b c d
p = evalBezier pHead pParam
-- straight line to child origin
p2x = lineBetween p childOrigin
p2x = lineBetween (p ^-^ vectR) childOrigin
-- straight line from child origin
x2p = lineBetween childOrigin p
x2p = lineBetween childOrigin' p
lineBetween a = CubicBezier a a a

1
src/Reanimate/Raster.hs
View file

@ -13,7 +13,6 @@ module Reanimate.Raster
) where
import Codec.Picture
import Codec.Picture.Types (dynamicMap)
import Control.Lens ((&), (.~))
import Control.Monad
import qualified Data.ByteString as B

3
src/Reanimate/Svg.hs
View file

@ -12,7 +12,6 @@ import Control.Monad.State
import Graphics.SvgTree hiding (height, line, path, use,
width)
import Linear.V2 hiding (angle)
import Data.Monoid ((<>))
import Reanimate.Constants
import Reanimate.Animation (SVG)
import Reanimate.Svg.Constructors
@ -179,7 +178,7 @@ svgGlyphs = worker id defaultSvg
let acc' sub = acc (GroupTree $ g & groupChildren .~ [sub])
attr' = (g^.drawAttributes) `mappend` attr
in concatMap (worker acc' attr') (g ^. groupChildren)
t -> [(acc, attr, t)]
t -> [(acc, (t^.drawAttributes) `mappend` attr, t)]
{-| Convert primitive SVG shapes (like those created by 'mkCircle', 'mkRect', 'mkLine' or
'mkEllipse') into SVG path. This can be useful for creating animations of these shapes being

37
src/Reanimate/Svg/Constructors.hs
View file

@ -9,6 +9,7 @@ module Reanimate.Svg.Constructors
, mkPathString
, mkPathText
, mkLinePath
, mkLinePathClosed
, mkClipPath
, mkText
-- * Grouping shapes and definitions
@ -19,6 +20,7 @@ module Reanimate.Svg.Constructors
, withId
, withStrokeColor
, withStrokeColorPixel
, withStrokeDashArray
, withStrokeLineJoin
, withFillColor
, withFillColorPixel
@ -30,6 +32,7 @@ module Reanimate.Svg.Constructors
, center
, centerX
, centerY
, centerUsing
, translate
, rotate
, rotateAroundCenter
@ -168,9 +171,7 @@ flipYAxis = scaleXY 1 (-1)
-- | Translate given image so that the center of its bouding box coincides with coordinates
-- @(0, 0)@.
center :: Tree -> Tree
center t = translate (-x-w/2) (-y-h/2) t
where
(x, y, w, h) = boundingBox t
center t = centerUsing t t
-- | Translate given image so that the X-coordinate of the center of its bouding box is 0.
centerX :: Tree -> Tree
@ -184,6 +185,11 @@ centerY t = translate 0 (-y-h/2) t
where
(_x, y, _w, h) = boundingBox t
centerUsing :: Tree -> Tree -> Tree
centerUsing a = translate (-x-w/2) (-y-h/2)
where
(x, y, w, h) = boundingBox a
-- | Create 'Texture' based on SVG color name.
-- See <https://en.wikipedia.org/wiki/Web_colors#X11_color_names> for the list of available names.
-- If the provided name doesn't correspond to valid SVG color name, white-ish color is used.
@ -200,6 +206,9 @@ withStrokeColor color = strokeColor .~ pure (mkColor color)
withStrokeColorPixel :: PixelRGBA8 -> Tree -> Tree
withStrokeColorPixel color = strokeColor .~ pure (ColorRef color)
withStrokeDashArray :: [Double] -> Tree -> Tree
withStrokeDashArray arr = strokeDashArray .~ pure (map Num arr)
-- | See <https://developer.mozilla.org/en-US/docs/Web/SVG/Attribute/stroke-linejoin>
withStrokeLineJoin :: LineJoin -> Tree -> Tree
withStrokeLineJoin ljoin = strokeLineJoin .~ pure ljoin
@ -317,6 +326,16 @@ mkLinePath ((startX, startY):rest) =
cmds = [ MoveTo OriginAbsolute [V2 startX startY]
, LineTo OriginAbsolute [ V2 x y | (x, y) <- rest ] ]
-- | Create a path from a list of @(x, y)@ coordinates of points along the path.
mkLinePathClosed :: [(Double, Double)] -> Tree
mkLinePathClosed [] = mkGroup []
mkLinePathClosed ((startX, startY):rest) =
PathTree $ defaultSvg & pathDefinition .~ cmds
where
cmds = [ MoveTo OriginAbsolute [V2 startX startY]
, LineTo OriginAbsolute [ V2 x y | (x, y) <- rest ]
, EndPath ]
-- | Rectangle with a uniform color and the same size as the screen.
--
-- Example:
@ -340,15 +359,15 @@ mkBackgroundPixel pixel =
-- rows. Each row can contain different number of cells.
gridLayout :: [[Tree]] -> Tree
gridLayout rows = mkGroup
[ translate (-screenWidth/2+colSep*nCol)
(screenHeight/2-rowSep*nRow)
[ translate (-screenWidth/2+colSep*nCol + colSep*0.5)
(screenHeight/2-rowSep*nRow - rowSep*0.5)
elt
| (nRow, row) <- zip [1..] rows
| (nRow, row) <- zip [0..] rows
, let nCols = length row
colSep = screenWidth / fromIntegral (nCols+1)
, (nCol, elt) <- zip [1..] row ]
colSep = screenWidth / fromIntegral nCols
, (nCol, elt) <- zip [0..] row ]
where
rowSep = screenHeight / fromIntegral (nRows+1)
rowSep = screenHeight / fromIntegral nRows
nRows = length rows
-- | Insert a native text object anchored at the middle.

16
stack-lts-13.yaml Normal file
View file

@ -0,0 +1,16 @@
resolver: lts-13.19
allow-newer: false
packages:
- .
extra-deps:
- reanimate-svg-0.9.8.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
- hmatrix-0.20.0.0@sha256:d79a9218e314f1a2344457c3851bd1d2536518ecb5f1a2fcd81daa45e46cd025,4870
- websockets-0.12.7.0@sha256:fc96169cc8268d3efb93bdc1e0b060dd88ce932237b837904118777a90d6db80,7863
- clock-0.8@sha256:b4ae207e2d3761450060a0d0feb873269233898039c76fceef9cc1a544067767,4113

14
stack-lts-14.yaml Normal file
View file

@ -0,0 +1,14 @@
resolver: lts-14.27
allow-newer: false
packages:
- .
extra-deps:
- reanimate-svg-0.9.8.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
- hmatrix-0.20.0.0@sha256:d79a9218e314f1a2344457c3851bd1d2536518ecb5f1a2fcd81daa45e46cd025,4870

2
stack.yaml
View file

@ -1,4 +1,4 @@
resolver: lts-13.19
resolver: lts-15.4
allow-newer: false

47
stack.yaml.lock
View file

@ -1,47 +0,0 @@
# This file was autogenerated by Stack.
# You should not edit this file by hand.
# For more information, please see the documentation at:
# https://docs.haskellstack.org/en/stable/lock_files
packages:
- completed:
hackage: reanimate-svg-0.9.8.0@sha256:541ea5d9677d9ad055ea0b47ce3ebdd585d7197e5bdec248a53a5d53a281519f,2507
pantry-tree:
size: 1117
sha256: c7a66d694d4486c908018875515206e598946cbf2fba38fa72738231aaf3cf3e
original:
hackage: reanimate-svg-0.9.8.0
- completed:
hackage: chiphunk-0.1.2.1@sha256:1e17e0e3d2fc91317e2cfcb8fd9f7add4a5038d76a4ae6b444641c0a2adec881,4619
pantry-tree:
size: 6333
sha256: 505d40929bd792b8003c0d03ec93d99ecba5e7927fbc6193824c0dab92c54dec
original:
hackage: chiphunk-0.1.2.1
- completed:
hackage: cubicbezier-0.6.0.6@sha256:2191ff47144d9a13a2784651a33d340cd31be1926a6c188925143103eb3c8db3,2243
pantry-tree:
size: 1029
sha256: 1a6e3a16c3d33c8baac53b4345b29695ff90cf2b559cfed349ab1bc454949632
original:
hackage: cubicbezier-0.6.0.6@sha256:2191ff47144d9a13a2784651a33d340cd31be1926a6c188925143103eb3c8db3
- completed:
hackage: fast-math-1.0.2@sha256:91181eb836e54413cc5a841e797c42b2264954e893ea530b6fc4da0dccf6a8b7,1387
pantry-tree:
size: 421
sha256: f22194924e26d3a27faef9cc9943ebb703c508d8ee9c9c7507618e39d62853af
original:
hackage: fast-math-1.0.2@sha256:91181eb836e54413cc5a841e797c42b2264954e893ea530b6fc4da0dccf6a8b7
- completed:
hackage: matrices-0.5.0@sha256:b2761813f6a61c84224559619cc60a16a858ac671c8436bbac8ec89e85473058,1656
pantry-tree:
size: 1202
sha256: f5f030a2c4641dfbee6e403590493a30d7de7bca15d1c304616ff186f64860c8
original:
hackage: matrices-0.5.0@sha256:b2761813f6a61c84224559619cc60a16a858ac671c8436bbac8ec89e85473058
snapshots:
- completed:
size: 498155
url: https://raw.githubusercontent.com/commercialhaskell/stackage-snapshots/master/lts/13/19.yaml
sha256: b9367a80d4393d02e58a46b8a9fdfbd7bc19f59c0c2bbf90034ba15cf52cf213
original: lts-13.19

97
test/Helpers.hs Normal file
View file

@ -0,0 +1,97 @@
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -w #-}
-- {-# LANGUAGE TemplateHaskell #-}
module Helpers where
import qualified Data.Vector as V
import Linear.V2
import Linear.Vector
import Numeric
import Data.List
import Test.QuickCheck
import Test.Tasty
import Test.Tasty.QuickCheck
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Math.EarClip
import Reanimate.Math.SSSP
-- [(0.8,1),(0.04,0.1),(0.05,1)]
newtype PolyParam = PolyParam { unParam :: (Double, Double) }
deriving (Show)
instance Arbitrary PolyParam where
arbitrary = fmap PolyParam $
(,) <$> choose (0.01,0.99)
<*> choose (0.01,1)
shrink (PolyParam (a,b)) =
[ PolyParam (a', b')
| a' <- a : shrinkDouble a
, b' <- if a==a' then b : shrinkDouble b else b : shrinkDouble b
, a' >= 0.01, b' >= 0.01
]
data Parameters = Parameters [(Double, Double)]
deriving (Show)
instance Arbitrary Parameters where
arbitrary = do
PolyParam e1 <- arbitrary
PolyParam e2 <- arbitrary
PolyParam e3 <- arbitrary
PolyParam e4 <- arbitrary
rest <- arbitrary
return $ Parameters $ e1:e2:e3:e4:map unParam rest
shrink (Parameters xs) =
map Parameters (delete xs $ shrinkListByOne xs) ++
map Parameters (delete xs $
sequence $ map (map unParam . shrink . PolyParam) xs)
shrinkListByOne :: [a] -> [[a]]
shrinkListByOne x | length x <= 4 = [x]
shrinkListByOne x = x : zipWith (++) (inits x) (tails $ drop 1 x)
shrinkDouble :: Double -> [Double]
shrinkDouble x =
let s = showFFloat Nothing x "" in
if length s < 4
then []
else [read (take (length s-1) s)]
instance Arbitrary Polygon where
arbitrary =
frequency
[ (1, elements premade)
, (2, pMkWinding <$> choose (1,100))
, (99, randomPoly)
]
where
premade = [shape1, shape2, shape3, shape4, shape5, shape6
,shape7, shape8, shape9, shape11 ]
randomPoly = do
let genParameters = do
angMod <- choose (0.01, 0.99)
rMod <- choose (0.01, 0.99)
return (angMod, rMod)
PolyParam e1 <- arbitrary
PolyParam e2 <- arbitrary
PolyParam e3 <- arbitrary
PolyParam e4 <- arbitrary
es <- listOf genParameters
pure $ pGenerate (e1:e2:e3:e4:es)
-- shrink p = filter isSimple $
-- map V.fromList (sequence (map shrinkV2 (V.toList p)))
shrinkV2 (V2 a b) = V2 <$> take 1 (shrinkRealFrac a) <*> take 1 (shrinkRealFrac b)
instance (RealFrac a, Arbitrary a) => Arbitrary (V2 a) where
arbitrary = do
(a,b) <- arbitrary
pure $ V2 a b
shrink = shrinkV2
-- shrink (V2 a b) =
-- [ V2 a b | (a',b') <- shrink (a,b) ]
percent :: (Testable prop) => (Rational -> prop) -> Property
percent gen = forAll (choose (0,1::Double)) $ gen . realToFrac

72
test/Properties.hs Normal file
View file

@ -0,0 +1,72 @@
{-# LANGUAGE FlexibleInstances #-}
{-# LANGUAGE TypeSynonymInstances #-}
{-# OPTIONS_GHC -w #-}
{-# LANGUAGE TemplateHaskell #-}
module Properties where
import qualified Data.Vector as V
import Linear.V2
import Linear.Vector
import Test.QuickCheck
import Test.Tasty
import Test.Tasty.QuickCheck
import Reanimate.Math.Common
import Reanimate.Math.Polygon
import Reanimate.Math.EarClip
import Reanimate.Math.SSSP
import Helpers
import Debug.Trace
prop_pGenerate (PolyParam a) (PolyParam b) (PolyParam c) (PolyParam d) =
pIsSimple $ pGenerate [a,b,c,d]
prop_pIsSimple p = pIsSimple p
prop_isBetween a b = percent $ \t ->
a /= b ==>
isBetween (lerp t a b) (a,b) &&
isBetween a (a,b) &&
isBetween b (a,b)
prop_isBetweenAxis a = percent $ \t ->
let b = a + V2 0 1
c = a + V2 1 0 in
isBetween (lerp t a b) (a,b) &&
isBetween (lerp t a c) (a,c)
prop_isNotBetween a b = forAll (choose (1.1,100)) $ \t ->
a /= 0 && b /= 0 ==>
not (isBetween (lerp (realToFrac (t::Double)) a b :: V2 Rational) (a,b))
prop_winding = forAll (choose (1,100)) $ \n -> pIsSimple (pMkWinding n)
prop_ccw p = label (if pIsConvex p then "convex" else "concave") $ pIsCCW p
prop_rev_ccw p = not $ pIsCCW (mkPolygon $ V.reverse $ polygonPoints p)
prop_pCycle p = forAll (choose (0,1)) $ \t ->
pIsSimple (pCycle p t)
prop_validEarClip p = isValidTriangulation p (earClip $ pRing p)
-- dualToTrangulation . dual = id
prop_dualInv p =
let t = polygonTriangulation p in
dualToTriangulation (pRing p) (dual (polygonOffset p) t) == t
prop_ssspEq p =
pSize p < 20 ==>
naive (pRing p) == polygonSSSP p V.! 0
prop_ssspVisibilityLength (Parameters xs) =
let p = pGenerate xs in
pSize (ssspVisibility p) <= pSize p
prop_dualCycle p = forAll (choose (0,pSize p-1)) $ \n ->
let p1 = pSetOffset p n
p2 = pCopy p1
in and [ pParent p1 0 i == pParent p2 0 i | i <- [0.. pSize p-1] ]
return []
all_props :: TestTree
all_props = testProperties "properties" $allProperties

3
test/Spec.hs
View file

@ -3,10 +3,11 @@ module Main (main) where
import Test.Tasty
import UnitTests
import Properties
main :: IO ()
main = do
tests1 <- unitTestFolder "examples/"
tests2 <- compileTestFolder "examples/"
tests3 <- compileVideoFolder "videos/"
defaultMain $ testGroup "tests" [tests1, tests2, tests3]
defaultMain $ testGroup "tests" [tests1, tests2, tests3, all_props]

25
test/UnitTests.hs
View file

@ -5,9 +5,12 @@ module UnitTests
) where
import Control.Exception
import qualified Data.ByteString as BS
import qualified Data.ByteString.Lazy as LBS
import Data.List (sort)
import Reanimate.Misc (withTempDir, withTempFile)
import qualified Data.Text as T
import qualified Data.Text.IO as T
import Reanimate.Misc (runCmd, withTempDir, withTempFile)
import System.Directory
import System.Exit
import System.FilePath
@ -22,7 +25,7 @@ unitTestFolder path = do
files <- sort <$> getDirectoryContents path
mbWDiff <- findExecutable "wdiff"
let diff = case mbWDiff of
Nothing -> ["diff", "--strip-trailing-cr"]
Nothing -> ["diff", "--strip-trailing-cr"]
Just wdiff -> [wdiff, "--no-common"]
return $ testGroup "animate"
[ goldenVsStringDiff file (\ref new -> diff ++ [ref, new]) fullPath (genGolden hsPath)
@ -34,21 +37,23 @@ unitTestFolder path = do
genGolden :: FilePath -> IO LBS.ByteString
genGolden path = withTempDir $ \tmpDir -> withTempFile ".exe" $ \tmpExecutable -> do
let ghcOpts = ["-rtsopts", "--make", "-O2"] ++
let ghcOpts = ["-rtsopts", "--make", "-O0", "-Werror", "-Wall"] ++
["-odir", tmpDir, "-hidir", tmpDir, "-o", tmpExecutable]
runOpts = ["+RTS", "-M1G"]
-- XXX: Check for errors.
_ <- readProcessWithExitCode "stack" (["ghc","--", path] ++ ghcOpts) ""
-- ret <- runCmd_ "stack" $ ["ghc", "--"] ++ ghcOptions tmpDir ++ [self, "-o", tmpExecutable]
-- ["-rtsopts", "--make", "-threaded", "-O2"] ++
-- ["-odir", tmpDir, "-hidir", tmpDir]
(inh, outh, errh, _pid) <- runInteractiveProcess tmpExecutable (["test"] ++ runOpts)
runCmd "stack" $ ["ghc","--", path] ++ ghcOpts
(inh, outh, errh, pid) <- runInteractiveProcess tmpExecutable (["test"] ++ runOpts)
Nothing Nothing
-- hSetBinaryMode outh True
-- hSetNewlineMode outh universalNewlineMode
hClose inh
hClose errh
LBS.hGetContents outh
out <- BS.hGetContents outh
err <- T.hGetContents errh
code <- waitForProcess pid
case code of
ExitSuccess -> return $ LBS.fromChunks [out]
ExitFailure{} -> error $ "Failed to run: " ++ T.unpack err
compileTestFolder :: FilePath -> IO TestTree
compileTestFolder path = do

27
videos/map-projection/distortion.hs
View file

@ -23,19 +23,21 @@ import Graphics.SvgTree (PathCommand (..), Tree (None))
import Reanimate
import Reanimate.Animation
import Reanimate.GeoProjection
import Reanimate.Raster
import Reanimate.Scene
import System.IO.Unsafe
main :: IO ()
main = seq equirectangular $ reanimate $ sceneAnimation $ do
main = seq equirectangular $ reanimate $ setDuration 59 $ sceneAnimation $ do
bg <- newSpriteSVG $ mkBackground "white"
spriteZ bg (-1)
prevProj <- newVar equirectangularP
txtVar <- newVar "Equirectangular"
txtS <- newSprite $ renderLabel <$> unVar txtVar
txtFade <- spriteVar txtS 1 withGroupOpacity
spriteZ txtS 2
let pushInterp = pushMerge' (\a b t -> interpP src a b t)
pushMerge' fn label proj = do
let pushInterp label proj = do
fork $ do
tweenVar txtFade 0.2 $ \v -> fromToS v 0 . curveS 2
writeVar txtVar label
@ -43,10 +45,15 @@ main = seq equirectangular $ reanimate $ sceneAnimation $ do
prev <- readVar prevProj
play $ pauseAtEnd waitT $ signalA (curveS 2) $
mkAnimation morphT $ \t ->
mkGroup $
[ scaleToSize screenWidth screenHeight $
embedImage $ fn prev proj t
, grid $ mergeP prev proj t ]
let imgKey = (label, projectionLabel prev, projectionLabel proj, t
,"distort"::String)
imgFile = cacheImage imgKey $
interpP src prev proj t
in mkGroup $
[ --scaleToSize screenWidth screenHeight $
-- embedImage $ interP src prev proj t
mkImage screenWidth screenHeight imgFile
, grid $ mergeP prev proj t ]
writeVar prevProj proj
pushT mkLabel mkProj = do
fork $ tweenVar txtVar morphT $ \v t -> if t > 0 then mkLabel t else v
@ -114,12 +121,12 @@ renderLabel label =
translate (-screenWidth/2) (-screenHeight/2) $
translate 0 (svgHeight ref) svgTxt
equirectangular :: Image PixelRGB8
equirectangular :: Image PixelRGBA8
equirectangular = unsafePerformIO $ do
dat <- BS.readFile "earth.jpg"
dat <- BS.readFile "earth-extreme.jpg"
case decodeJpeg dat of
Left err -> error err
Right img -> return $ convertRGB8 img
Right img -> return $ convertRGBA8 img
toRads :: Double -> Double
toRads dec = dec/180 * pi

286
videos/morph/morph.hs Normal file
View file

@ -0,0 +1,286 @@
#!/usr/bin/env stack
-- stack --resolver lts-15.04 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><s) (concat m)
bM = (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
]

434
videos/morph/vis.hs Normal file
View file

@ -0,0 +1,434 @@
#!/usr/bin/env stack
-- stack --resolver lts-15.04 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.PolyShape (svgToPolygons)
p1 = pScale 2 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
-- p1 = pSetOffset (addPoints 2 shape13) 0
-- p2 = pScale 0.5 $ centerPolygon shape20
-- p1 = centerPolygon shape2
p2 = pScale 2 $ pAtCenter $ pAddPoints 0 (pSetOffset shape14 0)
main :: IO ()
main = reanimate $ sceneAnimation $ do
bg <- newSpriteSVG $ mkBackground "black"
spriteZ bg (-1)
-- play $ drawVisibleFrom triangle
-- play $ drawVisibleFrom shape1
-- play $ drawVisibleFrom shape2
-- play $ drawVisibleFrom shape3
-- play $ drawVisibleFrom shape4
-- play $ drawVisibleFrom shape5
-- play $ drawVisibleFrom shape6
-- play $ drawOverlap shape20
-- play $ drawSSSP triangle naive
-- play $ drawSSSPFast shape5
-- fork $ play $ mapA (translate (-3) 0) $ drawSSSPVisibilityFast shape2
-- fork $ play $ mapA (translate (3) 0) $ drawSSSPVisibilityFast shape7
-- play $ drawOverlap $ fst $ split1Link (fst (split1Link shape7 3 7)) 0 2
-- fork $ play $ mapA (translate (-4) 0) $ drawSSSPVisibilityFast $ pSetOffset (addPoints 0 shape14) 0
-- fork $ play $ mapA (translate (4) 0) $ drawSSSPVisibilityFast $ pSetOffset (addPoints 2 shape13) 0
fork $ play $ staticFrame 1 $
translate (4) 0 $ mkGroup
[ withFillColor "grey" $ polygonShape p1
, polygonNumDots p1 ]
fork $ play $ staticFrame 1 $
let -- pOrigin = (addPoints 10 $ pSetOffset shape14 0)
pOrigin = pScale 2 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
--pOrigin = (addPoints 2 shape13)
p1 = pScale 1 $ pSetOffset pOrigin 0
p2 = snd $ split1Link p1 1 11 0
p3 = fst $ split1Link p2 0 2 1
-- p4 = fst $ split2Link p3 2 5
p4 = snd $ split1Link p3 0 9 1
p5 = fst $ split1Link p4 1 3 0
p6 = fst $ split2Link p5 0 2
p7 = fst $ split1Link p6 1 3 0
p8 = fst $ split2Link p7 0 2
p9 = fst $ split1Link p7 1 3 0
p = p4
in mkGroup
[ withFillColor "grey" $ polygonShape p
, polygonNumDots p ]
-- newSpriteSVG $
-- let p1 = pSetOffset shape14 0
-- V2 x y = realToFrac <$> steiner2Link p1 2 6
-- in translate (-4) 0 $
-- translate x y $ withFillColor "red" $
-- mkCircle 0.1
-- newSpriteSVG $
-- let p1 = pSetOffset shape14 0
-- in translate (-4) 0 $
-- mkGroup
-- [ mkGroup
-- [ withStrokeColor "purple" $
-- mkLinePath [(x1,y1),(x2,y2)]
-- , translate x2 y2 $ withFillColor "red" $ mkCircle 0.1 ]
-- | (eA, eB) <- take 1 $ steiner2Edges p1 4 11
-- , let V2 x1 y1 = realToFrac <$> eA
-- V2 x2 y2 = realToFrac <$> eB
-- ]
-- play $ staticFrame 1 $
-- let p1 = pSetOffset shape14 0
-- in mkGroup
-- [ withFillColor "grey" $ polygonShape p1
-- , polygonNumDots p1
-- , drawWindowOverlap p1 1 6
-- -- , withStrokeColor "red" $ drawWindow (pSetOffset p1 2)
-- -- , withStrokeColor "blue" $ drawWindow (pSetOffset p1 6)
-- ]
-- play $ drawCompatible (pSetOffset (addPoints 2 shape13) 0) (pSetOffset shape14 0)
-- play $ drawSSSP shape2 naive
-- play $ drawSSSP shape3 naive
-- play $ drawSSSP shape4 naive
-- play $ drawSSSP shape5 naive
-- play $ drawSSSP (pScale 0.5 $ winding 10) (\p -> sssp p (dual (earClip p)))
-- play $ drawSSSP shape1 (\p -> sssp p (dual (earClip p)))
-- play $ drawTriangulation (pCycle shape5 0.2910962834555265) earClip'
-- play $ drawTriangulation shape5 earClip'
-- play $ mkAnimation 1 $ \t ->
-- let p = shape4
-- in polygonNumDots (pCycle p t)
-- play $ setDuration 20 $ drawSSSPVisibility $ pScale 1 $ shape7
-- let shapeI = head $ svgToPolygons 0.1 $ scale 8 $ center $ latex "I"
-- play $ animate $ \_ -> withFillColor "grey" $ polygonNumDots shapeI
-- play $ drawTriangulation (pScale 0.5 $ winding 10) earClip'
-- play $ staticFrame 1 $
-- mkGroup
-- [ withFillColor "grey" $ polygonShape shape13
-- , withFillColor "grey" $ polygonDots shape13 ]
-- let p = balloonP origin
-- origin = centerPolygon $ pScale 1 $ shiftLongestDiameter shapeI
-- mkB = balloon (scale 8 $ center $ latex "C")
-- inf = unsafeSVGToPolygon 0.1 $ (scale 8 $ center $ latex "$\\infty$")
-- cShape = shiftLongestDiameter $ unsafeSVGToPolygon 0.01 $ (scale 8 $ center $ latex "C")
-- _ <- newSpriteSVG $
-- translate (0) (3) $ withFillColor "white" $
-- center $ latex $ T.pack $ show (isSimple inf)
-- play $ pauseAtEnd 1 $ mkAnimation 3 $ \t ->
-- let inflated = p t in
-- translate (0) (0) $ withFillColor "white" $ withStrokeColor "red" $
-- withStrokeWidth (defaultStrokeWidth*0) $
-- polygonShape inf
-- renderTriangulation inf (earClip inf)
-- mkB t
-- scale 3 $ -- translate (-2.2) (-1) $
-- mkGroup
-- [ mkGroup []
-- , withStrokeWidth (defaultStrokeWidth*0.05) $
-- renderDual (pRing shape22) (dual (polygonOffset shape22) $ polygonTriangulation shape22)
-- , polygonNumDots shape22
-- -- renderTriangulation cShape (polygonTriangulation cShape)
-- ]
-- mkGroup
-- [ -- translate (-2) 0 $ withFillColor "white" $ polygonShape $ balloonP (min 1 $ t+0.1) origin
-- if False then mkGroup [] else translate (0) 0 $ mkGroup
-- [ withFillColor "white" $ polygonShape inflated
-- -- , polygonNumDots inflated
-- ]
-- -- , translate (2) 0 $ mkGroup
-- -- [ withFillColor "grey" $ polygonShape origin
-- -- , polygonNumDots origin
-- -- ]
-- ]
-- play $ mapA (withStrokeWidth (defaultStrokeWidth*0.2)) $ drawTriangulation cShape earClip'
-- play $ staticFrame 1 $ renderTriangulation shape3 earClip
-- play $ staticFrame 1 $ renderTriangulation shape4 earClip
-- play $ staticFrame 1 $ renderTriangulation shape5 earClip
-- play $ staticFrame 1 $ renderTriangulation shape6 earClip
-- let p = deoverlapPolygon shape23
-- newSpriteSVG $ translate (-3) 0 $ scale 2 $ center $ mkGroup
-- [ withFillColor "grey" $ polygonShape p
-- , polygonNumDots p
-- ]
-- wait 1
-- play $ mapA (scale 2 . withStrokeWidth (defaultStrokeWidth*0.2)) $
-- drawTriangulation p (earClip')
return ()
drawVisibility :: Polygon -> Animation
drawVisibility p' = mkAnimation 5 $ \t ->
let p = pCycle p' (t::Double) in
centerUsing (polygonShape p) $
mkGroup
[ withFillColor "grey" $ polygonShape p
, withFillColor "grey" $ polygonDots p
, withFillColor "white" $ mkLinePathClosed
[ (x,y) | V2 x y <- visibility (map (fmap realToFrac) $ V.toList $ polygonPoints p) ]
, let V2 x y = fmap realToFrac $ pAccess p 0 in
translate x y $ withFillColor "red" $ mkCircle 0.1
-- , withFillColor "blue" $ latex $ T.pack $ show (t)
]
drawSSSPVisibility :: Polygon -> Animation
drawSSSPVisibility p' = mkAnimation 5 $ \t ->
let p = pSetOffset p' (round $ t*(fromIntegral $ pSize p'-1))
vis = ssspVisibility p in
centerUsing (polygonShape p) $
mkGroup
[ withFillColor "grey" $ polygonShape p
-- , withFillColor "grey" $ polygonDots p
-- , withFillColor "white" $ polygonShape vis
, let V2 x y = fmap realToFrac $ pAccess p 0 in
translate x y $ withFillColor "red" $ mkCircle 0.1
-- , withFillColor "blue" $ latex $ T.pack $ show (t)
]
drawSSSPVisibilityFast :: Polygon -> Animation
drawSSSPVisibilityFast p' = mkAnimation 5 $ \t ->
let root = min (pSize p-1) $ (floor $ t*(fromIntegral $ pSize p))
p = pSetOffset p' root
vis = ssspVisibility p in
-- centerUsing (polygonShape p) $
mkGroup
[ withFillColor "grey" $ polygonShape p
, withFillColor "white" $ polygonShape vis
, withFillColor "grey" $ polygonNumDots p
, let V2 x y = fmap realToFrac $ pAccess p 0 in
translate x y $ withFillColor "red" $ mkCircle 0.09
-- , withFillColor "blue" $ latex $ T.pack $ show (t)
]
drawCompatible :: Polygon -> Polygon -> Animation
drawCompatible a b = sceneAnimation $ do
newSpriteSVG $ translate (-3) 0 $ mkGroup
[ withFillColor "grey" $ polygonShape a
, withFillColor "grey" $ polygonNumDots a
]
newSpriteSVG $ translate (3) 0 $ mkGroup
[ withFillColor "grey" $ polygonShape b
, withFillColor "grey" $ polygonNumDots b
]
let compat = compatiblyTriangulateP a b
forM_ compat $ \(l, r) -> do
fork $ play $ staticFrame 1 $
translate (-3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
withFillOpacity 0 $ polygonShape l
fork $ play $ staticFrame 1 $
translate (3) 0 $ withStrokeColor "white" $ withStrokeWidth (defaultStrokeWidth*0.2) $
withFillOpacity 0 $ polygonShape r
-- forM_ compat $ \(l, r) -> waitOn $ do
-- fork $ play $ staticFrame 1 $
-- translate (-3) 0 $
-- withFillColor "white" $ polygonShape l
-- fork $ play $ staticFrame 1 $
-- translate (3) 0 $
-- withFillColor "white" $ polygonShape r
drawOverlap :: Polygon -> Animation
drawOverlap p' = mkAnimation 5 $ \t ->
let p = pCycle p' (t::Double)
vis = ssspVisibility p
vis' = ssspVisibility p'
mWins = ssspWindows p
oWins = ssspWindows p'
-- (left, right) = split1Link p' 0 3
-- (left, right) = split2Link p' 0 2
-- (left, right) = splitNLink p' 0 [(TwoLink,2),(OneLink,3)]
-- (left, right) = splitNLink p' 0 [(OneLink,5),(TwoLink,3)]
-- sPoly = if t < 0.5 then left else right
in
centerUsing (polygonShape p) $
mkGroup
[ withFillColor "grey" $ polygonShape p
-- , withFillColor "grey" $ polygonDots p
, withFillOpacity 0.5 $ withFillColor "lightgreen" $ polygonShape vis
-- , withFillOpacity 0.5 $ withFillColor "cyan" $ polygonShape vis'
, let V2 x y = fmap realToFrac $ pAccess p 0 in
translate x y $ withFillColor "red" $ mkCircle 0.1
, let V2 x y = fmap realToFrac $ pAccess p' 0 in
translate x y $ withFillColor "red" $ mkCircle 0.1
-- , withFillColor "blue" $ latex $ T.pack $ show (t)
, mkGroup
[ withStrokeColor "red" $
mkLine (x1,y1) (x2,y2)
| (a,b) <- mWins
, let V2 x1 y1 = realToFrac <$> a
V2 x2 y2 = realToFrac <$> b
]
-- , mkGroup
-- [ withStrokeColor "blue" $
-- mkLine (x1,y1) (x2,y2)
-- | (a,b) <- oWins
-- , let V2 x1 y1 = realToFrac <$> a
-- V2 x2 y2 = realToFrac <$> b
-- ]
, mkGroup
[ withStrokeColor "green" $mkGroup
[ mkLine (x1,y1) (x2,y2)
, mkLine (x1',y1') (x2',y2') ]
| (a,b) <- mWins
, (i,j) <- oWins
, a == i || a == j || b == i || b == j
, not (sort [a,b] == sort [i,j])
, let V2 x1 y1 = realToFrac <$> a
V2 x2 y2 = realToFrac <$> b
V2 x1' y1' = realToFrac <$> i
V2 x2' y2' = realToFrac <$> j
]
-- , mkGroup
-- [ let V2 x y = realToFrac <$> link in
-- translate x y $
-- withFillColor "red" $
-- mkCircle 0.1
-- | link <- [steiner2Link p' 0 2, steiner2Link p' 5 3] ]
-- , withFillColor "blue" $ polygonShape sPoly
]
drawWindow :: Polygon -> SVG
drawWindow p =
let mWins = ssspWindows p
in
mkGroup
[ mkGroup
[ mkLine (x1,y1) (x2,y2)
| (a,b) <- mWins
, let V2 x1 y1 = realToFrac <$> a
V2 x2 y2 = realToFrac <$> b
]
]
drawWindowOverlap :: Polygon -> Int -> Int -> SVG
drawWindowOverlap p a b =
let aWins = ssspWindows (pAdjustOffset p a)
bWins = ssspWindows (pAdjustOffset p b)
in
mkGroup
[ mkGroup $
[ withStrokeColor "green" $mkGroup
[ mkLine (x1,y1) (x2,y2)
, mkLine (x1',y1') (x2',y2') ]
| (a,b) <- aWins
, (i,j) <- bWins
, a == i || a == j || b == i || b == j
, not (sort [a,b] == sort [i,j])
, let V2 x1 y1 = realToFrac <$> a
V2 x2 y2 = realToFrac <$> b
V2 x1' y1' = realToFrac <$> i
V2 x2' y2' = realToFrac <$> j
]
]
drawSSSP :: Polygon -> (Polygon -> SSSP) -> Animation
drawSSSP p gen = mkAnimation 5 $ \t -> centerUsing outline $
let p' = pCycles p !! (round $ t*(fromIntegral $ pSize p-1)) in
mkGroup
[ outline
, renderSSSP p' (gen p')
-- , let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
-- translate x y $ withFillColor "red" $ mkCircle 0.1
, withFillColor "grey" $ polygonNumDots $ p'
]
where
outline =
withFillColor "grey" $ mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) (V.toList (polygonPoints p) ++ [pAccess p 0]) ]
{-# INLINE drawSSSPFast #-}
drawSSSPFast :: Polygon -> Animation
drawSSSPFast p = mkAnimation 5 $ \t -> centerUsing outline $
let root = (round $ t*(fromIntegral $ pSize p-1))
d = dual root triangulation
sTree = sssp (pRing p) d in
mkGroup
[ outline
, renderSSSP p sTree
-- , let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
-- translate x y $ withFillColor "red" $ mkCircle 0.1
, withFillColor "grey" $ polygonNumDots $ p
]
where
triangulation = earClip $ pRing p
outline =
withFillColor "grey" $ mkLinePathClosed
[ (x,y) | V2 x y <- map (fmap realToFrac) (V.toList (polygonPoints p) ++ [pAccess p 0]) ]
drawVisibleFrom :: Polygon -> Animation
drawVisibleFrom p = mkAnimation 5 $ \t -> centerUsing (polygonShape p) $ mkGroup
[ withFillColor "grey" $ polygonShape p
, withFillColor "grey" $ polygonDots p
, renderVisibleFrom (pCycle p t)
, let V2 x y = fmap realToFrac $ pAccess (pCycle p t) 0 in
translate x y $ withFillColor "red" $ mkCircle 0.1
]
renderVisibleFrom :: Polygon -> SVG
renderVisibleFrom p = withStrokeColor "white" $ withFillColor "white" $ mkGroup
[ mkGroup
[ mkLine (ax,ay) (bx,by)
, translate bx by $ mkCircle 0.1 ]
| i <- visibilityArray (pRing p) V.! 0
, let V2 ax ay = fmap realToFrac $ pAccess p 0
V2 bx by = fmap realToFrac $ pAccess p i ]
drawTriangulation :: Polygon -> (Ring Rational -> [Triangulation]) -> Animation
drawTriangulation p gen = sceneAnimation $ do
forM_ (gen $ pRing p) $ \t -> play $ staticFrame 1 $ renderTriangulation p t
renderDual :: Ring Rational -> Dual -> SVG
renderDual ring d = case d of
Dual (a,b,c) l r -> mkGroup
[ withFillColor "blue" $ mkTrig a b c
, worker c a l
, worker b c r
]
where
mkTrig a b c =
let V2 x1 y1 = realToFrac <$> ringAccess ring a
V2 x2 y2 = realToFrac <$> ringAccess ring b
V2 x3 y3 = realToFrac <$> ringAccess ring c
in mkLinePathClosed [ (x1, y1), (x2,y2), (x3,y3) ]
worker p1 p2 EmptyDual = mkGroup []
worker p1 p2 (NodeDual x l r) = mkGroup
[ mkTrig p1 p2 x
, worker x p2 l
, worker p1 x r
]
-- data Dual = Dual (Int,Int,Int) -- (a,b,c)
-- DualTree -- borders ca
-- DualTree -- borders bc
-- deriving (Show)
--
-- data DualTree
-- = EmptyDual
-- | NodeDual Int -- axb triangle, a and b are from parent.
-- DualTree -- borders xb
-- DualTree -- borders ax
-- deriving (Show)