diff --git a/.azure/azure-windows-template.yml b/.azure/azure-windows-template.yml
index 3ff941f..05d84ef 100644
--- a/.azure/azure-windows-template.yml
+++ b/.azure/azure-windows-template.yml
@@ -62,8 +62,8 @@ jobs:
     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
+      stack --install-ghc build $ARGS --no-keep-going --fast --flag reanimate:-hmatrix --extra-lib-dirs="$CONDA/Library/bin/" --only-dependencies
+      stack --install-ghc test $ARGS --no-keep-going --fast --flag reanimate:-hmatrix --flag reanimate:test --extra-lib-dirs="$CONDA/Library/bin/" --only-dependencies
     env:
       OS_NAME: ${{ parameters.os }}
     displayName: 'Install dependencies'
@@ -71,7 +71,7 @@ jobs:
       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/"
+      stack test --flag reanimate:test --flag reanimate:-hmatrix --test-arguments="-j2" --fast --extra-lib-dirs="$CONDA/Library/bin/"
     env:
       OS_NAME: ${{ parameters.os }}
     displayName: 'Build & Test'
diff --git a/ChangeLog.md b/ChangeLog.md
index 9f361a8..3cbc532 100644
--- a/ChangeLog.md
+++ b/ChangeLog.md
@@ -18,6 +18,8 @@
 * Fix bug that forced GIF width to 320 pixels.
 * Add helper for creating custom viewboxes, withViewBox.
 * Add newSprite_, newSpriteSVG_.
+* Drop support for lts-11 and lts-12. Lts-15, lts-14 and lts-13 are supported.
+* Add common interface for polygon morphing.
 
 ## 0.3.1.0 -- 2020-05-12
 
diff --git a/reanimate.cabal b/reanimate.cabal
index ac344fb..ee1796b 100644
--- a/reanimate.cabal
+++ b/reanimate.cabal
@@ -66,7 +66,6 @@ library
                       Reanimate.Math.Common
                       Reanimate.Math.Polygon
                       Reanimate.Math.Triangulate
-                      Reanimate.Math.DCEL
                       Reanimate.Math.EarClip
                       Reanimate.Math.Visibility
                       Reanimate.Math.SSSP
@@ -100,6 +99,7 @@ library
                       Reanimate.GeoProjection
                       Reanimate.Builtin.Documentation
                       Reanimate.Builtin.Images
+                      Reanimate.Debug
                       Reanimate.Voice
   other-modules:      Reanimate.Cache
                       Reanimate.Driver
@@ -122,6 +122,7 @@ library
     build-depends: hmatrix >= 0.20.0.0
     exposed-modules:  Reanimate.Morph.Rigid
                       Reanimate.Math.Smooth
+                      Reanimate.Math.DCEL
   }
   ghc-options: -Wall -fno-ignore-asserts
 
diff --git a/src/Reanimate/Debug.hs b/src/Reanimate/Debug.hs
new file mode 100644
index 0000000..a2c13bf
--- /dev/null
+++ b/src/Reanimate/Debug.hs
@@ -0,0 +1,58 @@
+module Reanimate.Debug
+  ( traceSVG
+  , traceA
+  , playTraces
+  )
+where
+
+import           Data.IORef
+import           System.IO.Unsafe
+import           Control.Exception
+import           Reanimate.Animation
+import           Reanimate.Svg
+import           Reanimate.LaTeX
+import           Reanimate.Constants
+import           Text.Printf
+import qualified Data.Text                     as T
+
+{-# NOINLINE traceBuffer #-}
+traceBuffer :: IORef [Animation]
+traceBuffer = unsafePerformIO (newIORef [])
+
+{-# NOINLINE traceSVG #-}
+traceSVG :: SVG -> a -> a
+traceSVG = traceA . staticFrame (recip 60)
+
+{-# NOINLINE traceA #-}
+traceA :: Animation -> a -> a
+traceA a v = unsafePerformIO $ do
+  modifyIORef' traceBuffer (a :)
+  evaluate v
+
+{-# NOINLINE playTraces #-}
+playTraces :: a -> Animation
+playTraces v = unsafePerformIO $ do
+  _   <- evaluate v
+  lst <- atomicModifyIORef' traceBuffer (\x -> ([], reverse x))
+  let n = length lst :: Int
+  return $ foldr
+    seqA
+    (pause 0)
+    [ f `parA` staticFrame (duration f) (counter i n) | (i, f) <- zip [1 :: Int ..] lst ]
+ where
+  counter a b = mkGroup
+    [ withStrokeWidth defaultStrokeWidth
+    $ withStrokeColor "black"
+    $ translate 6.5 4
+    $ center
+    $ latex
+    $ T.pack
+    $ printf "%d/%d" a b
+    , withStrokeWidth 0
+    $ withFillColor "white"
+    $ translate 6.5 4
+    $ center
+    $ latex
+    $ T.pack
+    $ printf "%d/%d" a b
+    ]
diff --git a/src/Reanimate/Driver/Server.hs b/src/Reanimate/Driver/Server.hs
index 3fdffe7..cdb1271 100644
--- a/src/Reanimate/Driver/Server.hs
+++ b/src/Reanimate/Driver/Server.hs
@@ -108,15 +108,15 @@ slaveHandler verbose mbGHCPath extraGHCOpts conn self svgDir =
     sendTextData conn (T.pack "status\nCompiling")
     ret <- case mbGHCPath of
       Nothing -> do
-        let opts = ["ghc", "--"] ++ ghcOptions tmpDir ++ extraGHCOpts ++ [takeFileName self, "-o", tmpExecutable]
+        let args = ["ghc", "--"] ++ ghcOptions tmpDir ++ extraGHCOpts ++ [takeFileName self, "-o", tmpExecutable]
         when verbose $
-          putStrLn $ "Running: " ++ showCommandForUser "stack" opts
-        runCmd_ "stack" opts
+          putStrLn $ "Running: " ++ showCommandForUser "stack" args
+        runCmd_ "stack" args
       Just ghc -> do
-        let opts = ghcOptions tmpDir ++ extraGHCOpts ++ [takeFileName self, "-o", tmpExecutable]
+        let args = ghcOptions tmpDir ++ extraGHCOpts ++ [takeFileName self, "-o", tmpExecutable]
         when verbose $
-          putStrLn $ "Running: " ++ showCommandForUser ghc opts
-        runCmd_ ghc opts
+          putStrLn $ "Running: " ++ showCommandForUser ghc args
+        runCmd_ ghc args
     case ret of
       Left err ->
         sendTextData conn $ T.pack $ "error\n" ++ unlines (drop 3 (lines err))
diff --git a/src/Reanimate/Math/Common.hs b/src/Reanimate/Math/Common.hs
index 942411b..0350bcd 100644
--- a/src/Reanimate/Math/Common.hs
+++ b/src/Reanimate/Math/Common.hs
@@ -72,43 +72,48 @@ area2X (V2 a1 a2) (V2 b1 b2) (V2 c1 c2) =
             (V3 c1 c2 1))
 
 epsilon :: Fractional a => a
-epsilon = 1e-9
+epsilon = 1e-13
 
 epsEq :: (Ord a, Fractional a) => a -> a -> Bool
 epsEq a b = abs (a-b) < epsilon
 
+{-# INLINE isLeftTurn #-}
 -- 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
+  case compare (direction p1 p2 p3) 0 of
     LT -> True
     EQ -> False -- colnear
     GT -> False
 
+{-# INLINE isLeftTurnOrLinear #-}
 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
+  case compare (direction p1 p2 p3) 0 of
     LT -> True
     EQ -> True -- colnear
     GT -> False
 
+{-# INLINE isRightTurn #-}
 isRightTurn :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
 isRightTurn a b c = not (isLeftTurnOrLinear a b c)
 
+{-# INLINE isRightTurnOrLinear #-}
 isRightTurnOrLinear :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> Bool
 isRightTurnOrLinear a b c = not (isLeftTurn a b c)
 
+{-# INLINE direction #-}
 direction :: Fractional a => V2 a -> V2 a -> V2 a -> a
 direction p1 p2 p3 = crossZ (p3-p1) (p2-p1)
 
+{-# INLINE isInside #-}
 isInside :: (Fractional a, Ord a) => V2 a -> V2 a -> V2 a -> V2 a -> Bool
 isInside a b c d =
     s >= 0 && s <= 1 && t >= 0 && t <= 1
   where
     (s, t, _) = barycentricCoords a b c d
 
+{-# INLINE barycentricCoords #-}
 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)
diff --git a/src/Reanimate/Math/Compatible.hs b/src/Reanimate/Math/Compatible.hs
index 4d4441e..90a1c84 100644
--- a/src/Reanimate/Math/Compatible.hs
+++ b/src/Reanimate/Math/Compatible.hs
@@ -9,140 +9,155 @@ import           Linear.V2
 import           Linear.Vector
 import           Reanimate.Math.Common
 import           Reanimate.Math.Polygon
+import           Reanimate.Debug
+import           Reanimate.Math.Render
+import           Reanimate.Svg
+import           Reanimate.Animation
 
-import           Debug.Trace
+-- import           Debug.Trace
+
+truncateP :: V2 Rational -> V2 Rational
+truncateP = fmap (realToFrac . (realToFrac :: Rational -> Double))
 
 mkSteinerPoints :: V2 Rational -> V2 Rational -> Int -> [V2 Rational]
-mkSteinerPoints a b s_ = [ lerp (i/(s+1)) b a | i <- [1..s]]
-  where
-    s = fromIntegral s_
+mkSteinerPoints 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
+  (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 -> []
+  | 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
+  = truncateP $ 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
         ]
-    isForward (a,b) v =
-      not (isBetween a (b,v))
+  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]
+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
+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 :: 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"
+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.
@@ -206,46 +221,71 @@ data MeshPair = MeshPair Points Points Edges
 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)
+  | otherwise = traceSVG (showStep a b) $ compatiblyTriangulateP'
+    (pSetOffset a 0)
+    (pSetOffset a 0)
+    (pSetOffset b 0)
+
+showStep :: Polygon -> Polygon -> SVG
+showStep a b = mkGroup
+  [ translate (-3) 0
+    $ mkGroup [withFillColor "grey" $ polygonShape a, polygonNumDots a]
+  , translate 3    0
+    $ mkGroup [withFillColor "grey" $ polygonShape b, polygonNumDots b]
+  ]
 
 compatiblyTriangulateP' :: Polygon -> Polygon -> Polygon -> [(Polygon, Polygon)]
 compatiblyTriangulateP' aOrigin a b
-  | trace ("Polygon size: " ++ show (pSize a)) False = undefined
-  | n == 3 = trace ("Done") $
-    [(a,b)]
+  | n == 3 = traceSVG (showStep a b) $ {- 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)
+    -- trace ("aOneLink: " ++ show aOneLink) $
+    -- trace ("bOneLink: " ++ show bOneLink) $
+    -- trace ("aTwoLink: " ++ show aTwoLink) $
+    -- trace ("bTwoLink: " ++ show bTwoLink) $
+    traceSVG (showStep a b) $ case bestOneLink of
+    Nothing -> case bestTwoLink of
+      Nothing -> error $ "no 2-links"
+      Just (nodeL, nodeR) ->
+        {-trace
+            (show
+              ("two link" :: String, toOriginIndex nodeL, toOriginIndex nodeR)
+            )
+          $ -}
+            let (aL, aR) = if (nodeL, nodeR) `elem` aOneLink
+                  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
@@ -260,25 +300,26 @@ oneBendBetween _p _a _b = False
 --       case pParent p a n of
 --         i -> if i == a then n else obstructedBy i
 
-polygonTwoLinks :: Polygon -> [(Int,Int)]
+polygonTwoLinks :: Polygon -> [(Int, Int)]
 polygonTwoLinks p =
   [ (i, j)
-  | i <- [0 .. n-1]
-  , j <- [i+2 .. n-1]
-  , not (i==0 && j == n-1)
+  | 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
+  , 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.
+-- Points on a straight line should be 2-link even though they are not
+-- direct grandparents.
   , isTwoLink || isStraightLine
-  ] where n = pSize p
+  ]
+  where n = pSize p
 
-polygonOneLinks :: Polygon -> [(Int,Int)]
+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)
+  | i <- [0 .. pSize p - 1]
+  , j <- [i + 2 .. pSize p - 1]
+  , not (i == 0 && j == pSize p - 1)
   , pParent p i j == i
   ]
diff --git a/src/Reanimate/Math/DCEL.hs b/src/Reanimate/Math/DCEL.hs
index 780ca72..72bab3f 100644
--- a/src/Reanimate/Math/DCEL.hs
+++ b/src/Reanimate/Math/DCEL.hs
@@ -1,136 +1,897 @@
 {-# LANGUAGE RecordWildCards #-}
+{-# LANGUAGE TemplateHaskell #-}
+-- Code originally from: https://github.com/wereHamster/haskell-half-edge-mesh
 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           Codec.Picture.Types
+import           Control.Lens
+import           Control.Monad.State
+import           Control.Monad.Writer
+import           Data.List
+import qualified Data.List              as L
+import           Data.Map               (Map)
+import qualified Data.Map               as M
+import           Data.Maybe
+import qualified Data.Set               as S
+import qualified Data.Text              as T
+import qualified Data.Vector            as V
+import           Debug.Trace
+import           Linear.Metric
 import           Linear.V2
-import           Linear.V3
+import           Linear.Vector
+import           Reanimate
+import           Reanimate.Math.Common  (isInside, lineIntersect, triangleAngles)
+import           Reanimate.Math.Polygon
+import qualified Reanimate.Morph.Rigid  as Rigid
 import           Text.Printf
 
-import           Reanimate.Math.Common
+type VertexId = Int
+type EdgeId = Int
+type FaceId = Int
 
-type EdgeIdx = Int
-type VertexIdx = Int
-type FaceIdx = Int
+-- Invariants:
+--   No duplicate positions
+--   All edges have an opposite twin
+--   Edge prev/next keep the same direction (CW or CCW)
+--   twin (twin edge) = edge
+--   All next/prev links have the same face
+data Mesh a = Mesh
+    { _meshIdCounter :: Int
+    , _meshOuterFace :: FaceId
+    , _meshVertices  :: Map VertexId (Vertex a)
+    , _meshEdges     :: Map EdgeId Edge
+    , _meshFaces     :: Map FaceId Face
+    } deriving (Show)
 
-data DCEL = DCEL
-  { dcelFaces    :: Vector EdgeIdx
-  , dcelVertices :: Vector EdgeIdx
-  , dcelEdges    :: Vector Edge
-  , dcelPoints   :: Vector P
-  } deriving (Show)
+data Vertex a = Vertex
+    { _vertexId       :: VertexId
+    , _vertexPosition :: a
+    -- , _vertexEdge     :: EdgeId -- Index to CCW edge pointing to this vertex.
+    } deriving (Show)
 
 data Edge = Edge
-  { edgeVertex :: VertexIdx
-  , edgeFace   :: FaceIdx
-  , edgeTwin   :: EdgeIdx
-  , edgeNext   :: EdgeIdx
-  , edgePrev   :: EdgeIdx
-  } deriving (Show)
+    { _edgeId     :: EdgeId
+    , _edgeVertex :: VertexId
+    , _edgeTwin   :: EdgeId
+    , _edgeNext   :: EdgeId
+    , _edgePrev   :: EdgeId
+    , _edgeFace   :: FaceId
+    } deriving (Show)
 
-getPoint :: Int -> DCEL -> P
-getPoint i DCEL{..} = dcelPoints V.! (i `mod` V.length dcelPoints)
+data Face = Face
+    { _faceId   :: FaceId
+    , _faceEdge :: EdgeId
+    } deriving (Show)
 
-dcelSize :: DCEL -> Int
-dcelSize DCEL{..} = V.length dcelPoints
+type MeshM position a = State (Mesh position) a
 
+validMesh :: Mesh a -> [String]
+validMesh Mesh{..} = execWriter $ do
+  return ()
 
+makeLenses ''Mesh
+makeLenses ''Vertex
+makeLenses ''Edge
+makeLenses ''Face
 
-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"
-  ]
+meshGetEdge :: EdgeId -> Mesh a -> Edge
+meshGetEdge eid Mesh{..} = M.findWithDefault err eid _meshEdges
   where
-    infix 0 ~>
-    True ~> _    = Nothing
-    False ~> msg = Just msg
+    err = error $ "Edge not found: " ++ show eid
 
-ppDCEL :: DCEL -> String
-ppDCEL DCEL{..} =
-  unlines
-  [ "Faces: Edge"
+meshGetVertex :: VertexId -> Mesh a -> Vertex a
+meshGetVertex vid Mesh{..} = M.findWithDefault err vid _meshVertices
+  where
+    err = error $ "Vertex not found: " ++ show vid
+
+meshGetFace :: FaceId -> Mesh a -> Face
+meshGetFace fid Mesh{..} = M.findWithDefault err fid _meshFaces
+  where
+    err = error $ "Face not found: " ++ show fid
+
+meshAngles :: Mesh (V2 Double) -> [Double]
+meshAngles mesh@Mesh{..} =
+  [ faceMinAngle fid mesh
+  | fid <- M.keys _meshFaces
+  ]
+
+meshEdgeAngle :: Mesh (V2 Double) -> EdgeId -> Double
+meshEdgeAngle m eId =
+    ang
+  where
+    edge = meshGetEdge eId m
+    next = meshGetEdge (edge ^. edgeNext) m
+    prev = meshGetEdge (edge ^. edgePrev) m
+    v1 = meshGetVertex (prev ^. edgeVertex) m ^. vertexPosition
+    v2 = meshGetVertex (edge ^. edgeVertex) m ^. vertexPosition
+    v3 = meshGetVertex (next ^. edgeVertex) m ^. vertexPosition
+    (_, ang, _) = triangleAngles v1 v2 v3
+
+facePositions :: FaceId -> Mesh a -> [a]
+facePositions fid m =
+  [ meshGetVertex (edge^.edgeVertex) m ^. vertexPosition
+  | eid <- faceEdges fid m
+  , let edge = meshGetEdge eid m
+  ]
+
+faceEdges :: FaceId -> Mesh a -> [EdgeId]
+faceEdges fid m = worker (meshGetEdge lastEdge m ^.edgeNext)
+  where
+    face = meshGetFace fid m
+    lastEdge = face^.faceEdge
+    worker eid
+      | eid == lastEdge = [eid]
+      | otherwise = eid : worker (meshGetEdge eid m ^. edgeNext)
+
+faceMinAngle :: FaceId -> Mesh (V2 Double) -> Double
+faceMinAngle fid m = minimum (map (meshEdgeAngle m) (faceEdges fid m))
+
+ppMesh :: Show a => Mesh a -> String
+ppMesh Mesh{..} = unlines
+  [ printf "Outer face: %d" _meshOuterFace
+  , ""
+  , "Vertices:"
+  , printf "%4s %5s %10s" ("ID"::String) ("Edge"::String) ("Position"::String)
   , unlines
-    [ printf "%3d:  %3d" face edge
-    | (face, edge) <- zip [0::Int .. ] (V.toList dcelFaces)
+    [ printf "%4d %10s" _vertexId  (show _vertexPosition)
+    | Vertex{..} <- M.elems _meshVertices
     ]
-  , "Vertices: Edge"
+  , "Edges:"
+  , printf "%4s %3s %4s %4s %4s %4s"
+    ("ID"::String) ("Vtx"::String) ("Twin"::String) ("Next"::String)
+    ("Prev"::String) ("Face"::String)
   , unlines
-    [ printf "%3d:     %3d" vertex edge
-    | (vertex, edge) <- zip [0::Int ..] (V.toList dcelVertices)
+    [ printf "%4d %3d %4d %4d %4d %4d" _edgeId _edgeVertex _edgeTwin _edgeNext _edgePrev _edgeFace
+    | Edge{..} <- M.elems _meshEdges
     ]
-  , "Edges: Vertex Face Twin Next Prev"
+  , "Faces:"
   , 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))
+    [ printf "%4d %4d" _faceId _faceEdge
+    | Face{..} <- M.elems _meshFaces
     ]
   ]
-  where chunkSize = 5
 
-emptyDCEL :: DCEL
-emptyDCEL = undefined
+emptyMesh :: Mesh a
+emptyMesh = Mesh 1 0 M.empty M.empty M.empty
+
+-- addFace :: FaceId -> [a] -> MeshM a ()
+-- addFace parentFace subPoly = do
+--   innerFace <- createFace
+--   vIds <- mapM getVertex pts
+--   edges     <- forM (zip vIds (tail vIds ++ take 1 vIds)) $ \(v0, v1) -> do
+--     e <- do
+--       mbEdge <- findEdge' v0 v1
+--       case mbEdge of
+--         Just e  -> pure (e ^. edgeId)
+--         Nothing -> do
+--           e <- createEdge v1
+--           twin <- createEdge v0
+--           setTwinEdge e twin
+--           setFace twin parentFace
+--           return e
+--     setFace e innerFace
+--     modifyVertex (vertexEdge .~ e) v0
+--     pure e
+
+--   setFaceEdge innerFace (head edges)
+--   forM_ (zip edges (tail edges ++ take 1 edges)) $ \(e0, e1) ->
+--     linkEdges e0 e1
+
+polygonsMesh :: Eq a => [a] -> [[a]] -> MeshM a ()
+polygonsMesh outer trigs = do
+  polygonMeshOuter (reverse outer)
+  modify $ meshFaces .~ M.empty
+  forM_ trigs $ \pts -> do
+    innerFace <- createFace
+    vIds <- mapM getVertex pts
+    edges     <- forM (zip vIds (tail vIds ++ take 1 vIds)) $ \(v0, v1) -> do
+      e <- createEdge v1
+      setFace e innerFace
+      -- modifyVertex (vertexEdge .~ e) v0
+      pure e
+
+    setFaceEdge innerFace (head edges)
+    forM_ (zip edges (tail edges ++ take 1 edges)) $ \(e0, e1) ->
+      linkEdges e0 e1
+    forM edges $ \e -> do
+      edge <- getEdge e
+      prev <- getEdge (edge ^. edgePrev)
+      mbTwin <- findEdge' (edge ^. edgeVertex) (prev ^. edgeVertex)
+      case mbTwin of
+        Nothing -> return ()
+        Just twin ->
+          setTwinEdge e (twin^.edgeId)
+
+-- Positions must be in CW order
+polygonMesh :: [a] -> MeshM a ()
+polygonMesh vs = do
+  outerFace <- gets _meshOuterFace
+  innerFace <- createFace
+  vIds      <- mapM createVertex vs
+  edges     <- forM (zip vIds (tail vIds ++ take 1 vIds)) $ \(v0, v1) -> do
+    e  <- createEdge v1 -- Edge from v0 to v1
+    e' <- createEdge v0 -- Twin edge from v1 to v0
+    setFace e  innerFace
+    setFace e' outerFace
+    -- modifyVertex (vertexEdge .~ e) v0
+    setTwinEdge e e'
+    pure e
+
+  outerEdge <- _edgeTwin <$> getEdge (head edges)
+  setFaceEdge innerFace (head edges)
+  setFaceEdge outerFace outerEdge
+
+  forM_ (zip edges (tail edges ++ take 1 edges)) $ \(e0, e1) -> do
+    linkEdges e0 e1
+    e0' <- _edgeTwin <$> getEdge e0
+    e1' <- _edgeTwin <$> getEdge e1
+    linkEdges e1' e0'
+
+polygonMeshOuter :: [a] -> MeshM a ()
+polygonMeshOuter vs = do
+  outerFace <- gets _meshOuterFace
+  vIds      <- mapM createVertex vs
+  edges     <- forM vIds $ \v0 -> do
+    e' <- createEdge v0
+    setFace e' outerFace
+    pure e'
+
+  setFaceEdge outerFace (head edges)
+
+  let f = Face outerFace (head edges)
+  modify $ meshFaces %~ M.insert outerFace f
+
+  forM_ (zip edges (tail edges ++ take 1 edges)) $ \(e0, e1) -> do
+    linkEdges e0 e1
+
+
+newId :: MeshM a Int
+newId = do
+  counter <- gets _meshIdCounter
+  modify $ meshIdCounter %~ succ
+  return counter
+
+createVertex :: a -> State (Mesh a) VertexId
+createVertex position = do
+  k <- newId
+  let v = Vertex k position -- (error "_vertexEdge not set")
+  modify $ meshVertices %~ M.insert k v
+  return k
+
+getVertex :: Eq a => a -> State (Mesh a) VertexId
+getVertex position = do
+  vs <- gets _meshVertices
+  case L.find comparingPosition (M.assocs vs) of
+    Nothing     -> createVertex position
+    Just (k, _) -> return k
+  where comparingPosition (_k, v) = (_vertexPosition v) == position
+
+requireVertex :: VertexId -> MeshM a (Vertex a)
+requireVertex vid = do
+  ret <- gets (M.lookup vid . _meshVertices)
+  case ret of
+    Nothing -> error "Invalid vertex id"
+    Just v  -> pure v
+
+modifyVertex :: (Vertex a -> Vertex a) -> VertexId -> State (Mesh a) ()
+modifyVertex f k = modify $ meshVertices %~ M.adjust f k
+
+getEdge :: EdgeId -> State (Mesh a) Edge
+getEdge e = do
+  edges <- gets _meshEdges
+  maybe (bug edges) return $ M.lookup e edges
+    --return $ fromJust $ M.lookup e edges
+  where bug _edges = error $ "Can't find edge with id " ++ (show e)
+
+withEdge :: EdgeId -> (Edge -> State (Mesh a) ()) -> State (Mesh a) ()
+withEdge e f = f =<< getEdge e
+
+-- findEdge :: VertexId -> VertexId -> State (Mesh a) EdgeId
+-- findEdge v0 v1 = do
+--   edges <- gets _meshEdges
+--   case L.find comparingVertices (M.assocs edges) of
+--     Nothing     -> error $ "No edge between: " ++ show (v0, v1)
+--     Just (k, _) -> return k
+--   where comparingVertices (k, e) = (k, _edgeNext e) == (v0, v1)
+
+createEdge :: VertexId -> State (Mesh a) EdgeId
+createEdge v = do
+  k <- newId
+  let e = Edge k
+               v
+               (error $ "_edgeTwin not set " ++ show (k,v))
+               (error $ "_edgeNext not set " ++ show (k,v))
+               (error $ "_edgePrev not set " ++ show (k,v))
+               (error $ "_edgeFace not set " ++ show (k,v))
+  modify $ meshEdges %~ M.insert k e
+  return k
+
+modifyEdge :: (Edge -> Edge) -> EdgeId -> State (Mesh a) ()
+modifyEdge f k = modify $ meshEdges %~ M.adjust f k
+
+linkEdges :: EdgeId -> EdgeId -> State (Mesh a) ()
+linkEdges e0 e1 = do
+  setNextEdge e0 e1
+  setPreviousEdge e1 e0
+
+setNextEdge :: EdgeId -> EdgeId -> State (Mesh a) ()
+setNextEdge e0 e1 = modifyEdge (edgeNext .~ e1) e0
+
+setPreviousEdge :: EdgeId -> EdgeId -> State (Mesh a) ()
+setPreviousEdge e0 e1 = modifyEdge (edgePrev .~ e1) e0
+
+setTwinEdge :: EdgeId -> EdgeId -> MeshM a ()
+setTwinEdge e0 e1 = do
+  modifyEdge (edgeTwin .~ e1) e0
+  modifyEdge (edgeTwin .~ e0) e1
+
+setFace :: EdgeId -> FaceId -> State (Mesh a) ()
+setFace e0 f0 = modifyEdge (edgeFace .~ f0) e0
+
+updateFaces :: EdgeId -> FaceId -> State (Mesh a) ()
+updateFaces e0 f0 = do
+  e <- getEdge e0
+  setFace e0 f0
+  worker (e ^. edgeNext)
+ where
+  worker e1
+    | e0 == e1 = return ()
+    | otherwise = do
+      e <- getEdge e1
+      setFace e1 f0
+      worker (e ^. edgeNext)
+
+createFace :: State (Mesh a) FaceId
+createFace = do
+  k <- newId
+  let f = Face k (error "_faceEdge not set")
+  modify $ meshFaces %~ M.insert k f
+  return k
+
+modifyFace :: (Face -> Face) -> FaceId -> State (Mesh a) ()
+modifyFace f k = modify $ meshFaces %~ M.adjust f k
+
+setFaceEdge :: FaceId -> EdgeId -> MeshM a ()
+setFaceEdge f0 e0 = modifyFace (\f -> f { _faceEdge = e0 }) f0
+
+deleteFace :: FaceId -> MeshM a ()
+deleteFace f0 = modify $ meshFaces %~ M.delete f0
+
+getFace :: FaceId -> State (Mesh a) Face
+getFace f = do
+  faces <- gets _meshFaces
+  maybe bug return (M.lookup f faces)
+  where bug = error $ "Can't find face with id " ++ (show f)
+
+
+buildMesh :: State (Mesh a) b -> Mesh a
+buildMesh f = execState f emptyMesh
+
+numVertices :: Mesh a -> Int
+numVertices = M.size . _meshVertices
+
+numEdges :: Mesh a -> Int
+numEdges = M.size . _meshEdges
+
+numFaces :: Mesh a -> Int
+numFaces = M.size . _meshFaces
+
+vertices :: Mesh a -> [Vertex a]
+vertices = M.elems . _meshVertices
+
+outgoingEdges :: VertexId -> Mesh a -> [Edge]
+outgoingEdges v mesh =
+  [ meshGetEdge (edge ^. edgeTwin) mesh
+  | edge <- M.elems (_meshEdges mesh)
+  , _edgeVertex edge == v
+  ]
+
+
+-- for each inner node:
+--    get all neighbour nodes. sort by CCW
+--    compute laplacian
+--    compute angle-based
+--    check correctness
+--    move vertex
+
+--   p    e
+-- _ -> u -> v
+--   p'   e'
+-- _ <- u <- v
+--
+--   p    e1   e
+-- _ -> u -> t -> v
+--   p'   e1'  e'
+-- _ <- u <- t <- v
+--
+-- faces stay the same
+splitEdge :: VertexId -> EdgeId -> MeshM a EdgeId
+splitEdge t e = do
+  eEdge     <- getEdge e
+  eEdgeTwin <- getEdge (eEdge ^. edgeTwin)
+
+  e1        <- createEdge t
+  setFace e1 (eEdge ^. edgeFace)
+  -- modifyVertex (vertexEdge .~ e1) t
+
+  e1' <- createEdge (eEdgeTwin ^. edgeVertex)
+  setFace e1' (eEdgeTwin ^. edgeFace)
+  setTwinEdge e1 e1'
+
+  let oldPrev = eEdge ^. edgePrev
+      oldNext = eEdgeTwin ^. edgeNext
+
+  linkEdges oldPrev             e1
+  linkEdges (eEdge ^. edgeTwin) e1'
+  linkEdges e1                  e
+  linkEdges e1'                 oldNext
+
+  modifyEdge (edgeVertex .~ t) (eEdge ^. edgeTwin)
+
+  return e1
+
+splitTriangle :: V2 Double -> EdgeId -> MeshM (V2 Double) (FaceId, FaceId, FaceId, FaceId)
+splitTriangle vertex e = do
+  outer <- gets _meshOuterFace
+  edge <- getEdge e
+  next <- getEdge (edge^.edgeNext)
+  twin <- getEdge (edge^.edgeTwin)
+  next' <- getEdge (twin^.edgeNext)
+
+  t         <- createVertex vertex
+  splitEdge t e
+
+  smoothVertex t
+
+  (f1,f2) <- insertEdge (edge^.edgeFace) t (next^.edgeVertex)
+  (f3,f4) <- insertEdge (twin^.edgeFace) t (next'^.edgeVertex)
+  return (f1,f2,f3,f4)
+
+splitOuterTriangle :: V2 Double -> EdgeId -> MeshM (V2 Double) (FaceId, FaceId)
+splitOuterTriangle vertex e = do
+  outer <- gets _meshOuterFace
+  edge <- getEdge e
+  twin <- getEdge (edge^.edgeTwin)
+  next' <- getEdge (twin^.edgeNext)
+
+  t         <- createVertex vertex
+  splitEdge t e
+
+  (f1,f2) <- insertEdge (twin^.edgeFace) t (next'^.edgeVertex)
+  return (f1,f2)
+  -- return (outer, outer)
+
+findEdge :: FaceId -> VertexId -> MeshM a Edge
+findEdge f v = do
+  edges <- gets (M.elems . _meshEdges)
+  case
+      listToMaybe [ e | e <- edges, e ^. edgeFace == f, e ^. edgeVertex == v ]
+    of
+      Nothing -> error $ "Edge not found: " ++ show (f, v)
+      Just e  -> pure e
+
+findEdge' :: VertexId -> VertexId -> MeshM a (Maybe Edge)
+findEdge' v0 v1 = do
+  m <- get
+  pure $ listToMaybe
+    [ e
+    | e <- M.elems (m^.meshEdges)
+    , e ^. edgeVertex == v1
+    , let prev = meshGetEdge (e ^. edgePrev) m
+    , prev ^. edgeVertex == v0
+    ]
+
+findEdge'' :: [FaceId] -> VertexId -> MeshM a (Maybe Edge)
+findEdge'' fs v = do
+  edges <- gets (M.elems . _meshEdges)
+  pure $
+      listToMaybe [ e | e <- edges, e ^. edgeFace `notElem` fs, e ^. edgeVertex == v ]
 
 {-
-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
 -}
+deleteEdge :: EdgeId -> MeshM a ()
+deleteEdge e0 = do
+  e <- getEdge e0
+  e' <- getEdge (e^.edgeTwin)
+  linkEdges (e^.edgePrev) (e'^.edgeNext)
+  linkEdges (e'^.edgePrev) (e^.edgeNext)
+  updateFaces (e^.edgeNext) (e^.edgeFace)
+  setFaceEdge (e^.edgeFace) (e^.edgeNext)
+  deleteFace (e'^.edgeFace)
+  modify $ meshEdges %~ M.delete (e^.edgeId)
+  modify $ meshEdges %~ M.delete (e'^.edgeId)
+
+insertEdge :: FaceId -> VertexId -> VertexId -> MeshM a (FaceId, FaceId)
+insertEdge f0 v0 v1 = do
+  hv0 <- findEdge f0 v0
+  hv1 <- findEdge f0 v1
+  f1  <- createFace
+  f2  <- createFace
+  h1  <- createEdge v1
+  h2  <- createEdge v0
+  setFaceEdge f1 h1
+  setFaceEdge f2 h2
+  setTwinEdge h1 h2
+  linkEdges (hv0 ^. edgeId) h1
+  linkEdges (hv1 ^. edgeId) h2
+  linkEdges h1              (hv1 ^. edgeNext)
+  linkEdges h2              (hv0 ^. edgeNext)
+  updateFaces (hv0 ^. edgeId) f1
+  updateFaces (hv1 ^. edgeId) f2
+  modify $ meshFaces %~ M.delete f0
+  return (f1,f2)
+
+steinerNodes :: Mesh a -> [VertexId]
+steinerNodes Mesh{..} =
+  [ v
+  | v <- M.keys _meshVertices
+  , S.notMember v notSteiner
+  ]
+  where
+    notSteiner = S.fromList
+      [ _edgeVertex edge
+      | edge <- M.elems _meshEdges
+      , _edgeFace edge == _meshOuterFace ]
+
+-- findNeighbours :: VertexId -> Mesh a -> [Edge]
+-- findNeighbours vertex Mesh{..} =
+--   [ edge
+--   | edge <- M.elems _meshEdges
+--   , _edgeVertex edge == vertex
+--   ]
+
+{-
+For each internal node:
+  * Find neighbours
+  * Compute angle-based position
+  * Compute laplacian position
+  * Update if possible
+-}
+meshSmoothPosition :: Mesh (V2 Double) -> Mesh (V2 Double)
+meshSmoothPosition = execState worker
+  where
+    worker = gets steinerNodes >>= mapM_ smoothVertex
+
+smoothVertex :: VertexId -> MeshM (V2 Double) ()
+smoothVertex steiner = do
+    self <- _vertexPosition <$> requireVertex steiner
+    es <- gets (outgoingEdges steiner)
+    vs <- mapM (requireVertex . _edgeVertex) es
+    let ps = V.fromList (sortEdges self (map _vertexPosition vs))
+        angleBased = angleSmooth self ps
+        laplacian  = sum ps ^/ (fromIntegral $ length ps) -- laplacian
+    if isValidLocation self ps angleBased
+      then modifyVertex (vertexPosition .~ angleBased) steiner
+      else if isValidLocation self ps laplacian
+        then modifyVertex (vertexPosition .~ laplacian) steiner
+        else return ()
+  where
+    sortEdges :: V2 Double -> [V2 Double] -> [V2 Double]
+    sortEdges = sortOn . dir
+    -- 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
+  ]
+
+flipEdge :: Edge -> MeshM (V2 Double) (FaceId, FaceId)
+flipEdge e = do
+  e' <- getEdge (e^.edgeTwin)
+
+  v0 <- _vertexPosition <$> requireVertex (e^.edgeVertex)
+  v1 <- _vertexPosition <$> requireVertex (e'^.edgeVertex)
+
+  a' <- getEdge (e^.edgeNext)
+  b' <- getEdge (e'^.edgeNext)
+
+  v0' <- _vertexPosition <$> requireVertex (a'^.edgeVertex)
+  v1' <- _vertexPosition <$> requireVertex (b'^.edgeVertex)
+
+  case lineIntersect (v0, v1) (v0', v1') of
+    Nothing -> pure (e^.edgeFace, e'^.edgeFace)
+    Just{} -> do
+      deleteEdge (e^.edgeId)
+      insertEdge (e^.edgeFace) (a'^.edgeVertex) (b'^.edgeVertex)
+
+internalEdges :: Mesh a -> [EdgeId]
+internalEdges m =
+  [ edge^.edgeId
+  | edge <- M.elems (m^.meshEdges)
+  , edge^.edgeFace /= m^.meshOuterFace
+  , let twin = meshGetEdge (edge^.edgeTwin) m
+  , edge^.edgeVertex < twin^.edgeVertex
+  , twin^.edgeFace /= m^.meshOuterFace
+  ]
+
+outerEdges :: Mesh a -> [EdgeId]
+outerEdges m =
+  [ edge^.edgeId
+  | edge <- M.elems (m^.meshEdges)
+  , edge^.edgeFace == m^.meshOuterFace
+  ]
+
+longestEdge :: Mesh (V2 Double) -> Mesh (V2 Double) -> EdgeId
+longestEdge m1 m2 =
+    maximumBy cmp (internalEdges m1)
+  where
+    cmp e1 e2 =
+      compare
+        (max (edgeLength e1 m1) (edgeLength e1 m2))
+        (max (edgeLength e2 m1) (edgeLength e2 m2))
+
+edgeLength :: EdgeId -> Mesh (V2 Double) -> Double
+edgeLength eid m =
+  let edge = meshGetEdge eid m
+      twin = meshGetEdge (edge^.edgeTwin) m
+      p1 = _vertexPosition $ meshGetVertex (edge^.edgeVertex) m
+      p2 = _vertexPosition $ meshGetVertex (twin^.edgeVertex) m
+  in distance p1 p2
+
+applyCompatible :: (t -> Mesh a -> Maybe (Mesh a)) -> (Mesh a -> [t]) -> Mesh a -> Mesh a -> (Mesh a, Mesh a)
+applyCompatible fn vs m1 m2
+  | _meshIdCounter m1 /= _meshIdCounter m2 = error "invalid ID counter"
+applyCompatible fn vs m1 m2 = worker m1 m2 (vs m1)
+  where
+    worker l r [] = (l,r)
+    worker l r (eid:rest) =
+      case (,) <$> fn eid l <*> fn eid r of
+        Nothing       -> worker l r rest
+        Just (l', r') -> worker l' r' rest
+
+delaunayFlip :: Mesh (V2 Double) -> Mesh (V2 Double) -> (Mesh (V2 Double), Mesh (V2 Double))
+delaunayFlip = applyCompatible delaunayFlip' internalEdges
+
+delaunayFlip' :: EdgeId -> Mesh (V2 Double) -> Maybe (Mesh (V2 Double))
+delaunayFlip' eid m =
+  let edge = meshGetEdge eid m
+      twin = meshGetEdge (edge^.edgeTwin) m
+      f1 = edge^.edgeFace
+      f2 = twin^.edgeFace
+      beforeAng = min (faceMinAngle f1 m) (faceMinAngle f2 m)
+      ((f1',f2'), mAfter) = runState (flipEdge edge) m
+      afterAng = min (faceMinAngle f1' mAfter) (faceMinAngle f2' mAfter)
+  in
+  if (afterAng < beforeAng) || f1' == f1
+    then Nothing
+    else Just mAfter
+
+splitInternalEdges :: Mesh (V2 Double) -> Mesh (V2 Double) -> (Mesh (V2 Double), Mesh (V2 Double))
+splitInternalEdges = applyCompatible splitInternalEdge internalEdges
+
+splitInternalEdge :: EdgeId -> Mesh (V2 Double) -> Maybe (Mesh (V2 Double))
+splitInternalEdge eid m = evalState worker m
+  where
+    edgeFaces edge = do
+      twin <- getEdge (edge^.edgeTwin)
+      return (edge^.edgeFace, twin^.edgeFace)
+    worker = do
+      edge <- getEdge eid
+      twin <- getEdge (edge^.edgeTwin)
+      v0 <- _vertexPosition <$> requireVertex (edge^.edgeVertex)
+      v1 <- _vertexPosition <$> requireVertex (twin^.edgeVertex)
+      let middle = lerp 0.5 v0 v1
+      (f1,f2) <- edgeFaces edge
+      mBefore <- get
+      let beforeAng = min (faceMinAngle f1 mBefore) (faceMinAngle f2 mBefore)
+      (f3,f4,f5,f6) <- splitTriangle middle eid
+      mAfter <- get
+      let afterAng = minimum
+            [ faceMinAngle f3 mAfter
+            , faceMinAngle f4 mAfter
+            , faceMinAngle f5 mAfter
+            , faceMinAngle f6 mAfter ]
+      return ()
+      if (afterAng < beforeAng)
+        then return Nothing
+        else return (Just mAfter)
+
+splitLongestEdge :: Mesh (V2 Double) -> Mesh (V2 Double) -> (Mesh (V2 Double), Mesh (V2 Double))
+splitLongestEdge m1 m2 =
+    (splitInternalEdgeForced longest m1, splitInternalEdgeForced longest m2)
+  where
+    longest = longestEdge m1 m2
+
+splitInternalEdgeForced :: EdgeId -> Mesh (V2 Double) -> Mesh (V2 Double)
+splitInternalEdgeForced eid m = execState worker m
+  where
+    edgeFaces edge = do
+      twin <- getEdge (edge^.edgeTwin)
+      return (edge^.edgeFace, twin^.edgeFace)
+    worker = do
+      edge <- getEdge eid
+      twin <- getEdge (edge^.edgeTwin)
+      v0 <- _vertexPosition <$> requireVertex (edge^.edgeVertex)
+      v1 <- _vertexPosition <$> requireVertex (twin^.edgeVertex)
+      let middle = lerp 0.5 v0 v1
+      splitTriangle middle eid
+
+splitOuterEdges :: Mesh (V2 Double) -> Mesh (V2 Double) -> (Mesh (V2 Double), Mesh (V2 Double))
+splitOuterEdges = applyCompatible splitOuterEdge outerEdges
+
+splitOuterEdge :: EdgeId -> Mesh (V2 Double) -> Maybe (Mesh (V2 Double))
+splitOuterEdge eid m = evalState worker m
+  where
+    edgeFaces edge = do
+      twin <- getEdge (edge^.edgeTwin)
+      return (edge^.edgeFace, twin^.edgeFace)
+    worker = do
+      edge <- getEdge eid
+      twin <- getEdge (edge^.edgeTwin)
+      v0 <- _vertexPosition <$> requireVertex (edge^.edgeVertex)
+      v1 <- _vertexPosition <$> requireVertex (twin^.edgeVertex)
+      let middle = lerp 0.5 v0 v1
+      let f1 = twin^.edgeFace
+      mBefore <- get
+      let beforeAng = faceMinAngle f1 mBefore
+      (f2,f3) <- splitOuterTriangle middle eid
+      mAfter <- get
+      let afterAng = minimum
+            [ faceMinAngle f2 mAfter
+            , faceMinAngle f3 mAfter ]
+      if (afterAng < beforeAng)
+        then return Nothing
+        else return (Just mAfter)
+
+renderMesh :: Double -> Mesh (V2 Double) -> SVG
+renderMesh radius m = mkGroup
+  [ mkGroup
+    [ withStrokeColor "black" $
+      mkLine (x1, y1) (x2, y2)
+    | edge <- M.elems (_meshEdges m)
+    , _edgeFace edge /= _meshOuterFace m
+    , let twin = meshGetEdge (edge ^. edgeTwin) m
+    , let V2 x1 y1 = _vertexPosition (meshGetVertex (_edgeVertex edge) m)
+          V2 x2 y2 = _vertexPosition (meshGetVertex (_edgeVertex twin) m)
+    ]
+  -- , mkGroup
+  --   [ mkLinePathClosed
+  --     [
+  --     |
+  --     ] [(Double, Double)] -> Tree
+  --   | face <- M.elems _meshFaces
+  --   ]
+  , mkGroup
+    [ translate x y $ withFillColor "red" $ mkCircle radius
+    | vertex <- M.elems (_meshVertices m)
+    , let V2 x y = _vertexPosition vertex
+    ]
+  ]
+
+renderMeshEdges :: Mesh (V2 Double) -> SVG
+renderMeshEdges mesh@Mesh {..} = mkGroup
+  [ mkGroup
+    [ mkGroup
+      [ withStrokeColorPixel (promotePixel $ viridis $ ang/pi) $
+        mkLine (x1, y1) (x2, y2)
+      ]
+    | edge <- M.elems _meshEdges
+    , _edgeFace edge /= _meshOuterFace
+    , let twin = meshGetEdge (edge ^. edgeTwin) mesh
+          V2 x1 y1 = _vertexPosition (meshGetVertex (_edgeVertex edge) mesh)
+          V2 x2 y2 = _vertexPosition (meshGetVertex (_edgeVertex twin) mesh)
+          ang = meshEdgeAngle mesh (edge^.edgeId)
+    , ang/pi*180 < 2
+    ]
+  ]
+
+renderMeshSimple :: Double -> Mesh (V2 Double) -> SVG
+renderMeshSimple radius Mesh {..} = mkGroup
+  [ mkGroup
+    [ translate x y $ mkGroup
+      [ withFillColor "red" $ mkCircle 0.01
+      , scaleToHeight (radius*2) $ center $ latex $ T.pack $ show (_vertexId vertex) ]
+    | vertex <- M.elems _meshVertices
+    , let V2 x y = _vertexPosition vertex
+    ]
+  ]
+
+renderMeshColored :: Mesh (V2 Double) -> SVG
+renderMeshColored m@Mesh{..} = mkGroup
+  [ withFillColorPixel c $ mkLinePathClosed $
+    [ (x,y) | V2 x y <- facePositions fid m ]
+  | fid <- M.keys _meshFaces
+  , let ang = faceMinAngle fid m
+        bestAng = pi/3
+        score = ang/bestAng
+        c = promotePixel $ viridis score
+  ]
+-- facePositions :: FaceId -> Mesh a -> [a]
+-- faceMinAngle :: FaceId -> Mesh (V2 Double) -> Double
+
+renderMeshStats :: Mesh (V2 Double) -> SVG
+renderMeshStats mesh = mkGroup
+    [ latex $ T.pack $ printf "Min:  %.1f" (minAngle/pi*180)
+    , translate 0 (sep*1) $
+      latex $ T.pack $ printf "Mean: %.1f" (meanAngle/pi*180)
+    , translate 0 (sep*2) $
+      latex $ T.pack $ printf "Avg:  %.1f" (avgAngle/pi*180)
+    , translate 0 (sep*3) $
+      latex $ T.pack $ printf "Max:  %.1f" (maxAngle/pi*180)
+    ]
+  where
+    sep = -1
+    angs = meshAngles mesh
+    minAngle = minimum angs
+    meanAngle = L.sort angs !! (length angs `div` 2)
+    avgAngle = sum angs / (fromIntegral $ length angs)
+    maxAngle = maximum angs
+
+toRigidMesh :: Mesh (V2 Double) -> Mesh (V2 Double) -> Rigid.Mesh
+toRigidMesh meshA meshB = Rigid.Mesh
+    { meshPointsA = pointsA
+    , meshPointsB = pointsB
+    , meshOutline = outline
+    , meshTriangles = trigs }
+  where
+    pointsA = V.fromList $ map _vertexPosition $ M.elems (meshA^.meshVertices)
+    pointsB = V.fromList $ map _vertexPosition $ M.elems (meshB^.meshVertices)
+    outline = V.fromList
+      [ fromJust (V.elemIndex v pointsA)
+      | eid <- faceEdges (meshA^.meshOuterFace) meshA
+      , let edge = meshGetEdge eid meshA
+            v = _vertexPosition (meshGetVertex (edge^.edgeVertex) meshA)
+      ]
+    trigs = V.fromList
+      [ (aIdx, bIdx, cIdx)
+      | fid <- M.keys (meshA^.meshFaces)
+      , fid /= (meshA^.meshOuterFace)
+      , let edges = map (`meshGetEdge` meshA) (faceEdges fid meshA)
+            vs = map (`meshGetVertex` meshA) $ map _edgeVertex edges
+            ps = map _vertexPosition vs
+            [aIdx, bIdx, cIdx] = map (fromJust . (`V.elemIndex` pointsA)) ps
+            p = mkPolygon (V.fromList $ map (fmap realToFrac) ps)
+      , pIsSimple p || error "invalid polygon"
+      ]
+
diff --git a/src/Reanimate/Math/Polygon.hs b/src/Reanimate/Math/Polygon.hs
index a467c7f..21e3053 100644
--- a/src/Reanimate/Math/Polygon.hs
+++ b/src/Reanimate/Math/Polygon.hs
@@ -597,6 +597,7 @@ pIsCCW p = V.sum (pMapEdges fn p) < 0
   where
     fn (V2 x1 y1) (V2 x2 y2) = (x2-x1)*(y2+y1)
 
+{-# INLINE pRayIntersect #-}
 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)
diff --git a/src/Reanimate/Math/Smooth.hs b/src/Reanimate/Math/Smooth.hs
index e77f33a..cdee482 100644
--- a/src/Reanimate/Math/Smooth.hs
+++ b/src/Reanimate/Math/Smooth.hs
@@ -266,6 +266,8 @@ splitMeshEdges maxLen mesh = runST $ do
       insertEdge t j p
       insertEdge t c p
       insertEdge t d p
+  -- Add 'newPoints' to mesh points
+  -- Convert edges to triangles
   return undefined
  where
   splitEdge i j =
@@ -282,8 +284,8 @@ splitMeshEdges maxLen mesh = runST $ do
   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))
+  -- sortEdges :: V.Vector (V2 Double) -> V2 Double -> [Int] -> [Int]
+  -- sortEdges v pt = sortOn (\i -> dir 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))
+  -- dir :: V2 Double -> V2 Double -> Double
+  -- dir a b = (atan2 (crossZ (V2 0 1) (b - a)) (dot (V2 0 1) (b - a)))
diff --git a/src/Reanimate/Morph/Linear.hs b/src/Reanimate/Morph/Linear.hs
index 6c9a6cb..e6fbda0 100644
--- a/src/Reanimate/Morph/Linear.hs
+++ b/src/Reanimate/Morph/Linear.hs
@@ -43,7 +43,9 @@ closestLinearCorrespondence src' dst' =
         then worker x newScore xs
         else worker bestP bestPScore xs
     options = pCycles dst
-    score = V.sum . V.zipWith approxDist (polygonPoints src) . polygonPoints
+    score p = sum
+      [ approxDist (pAccess src n) (pAccess p n)
+      | n <- [0 .. pSize src-1] ]
 
 linearTrajectory :: Trajectory
 linearTrajectory (src,dst)
diff --git a/src/Reanimate/Morph/Rigid.hs b/src/Reanimate/Morph/Rigid.hs
index 381322f..1bc9290 100644
--- a/src/Reanimate/Morph/Rigid.hs
+++ b/src/Reanimate/Morph/Rigid.hs
@@ -76,6 +76,9 @@ import           Linear.V3
 import qualified Numeric.LinearAlgebra         as Matrix
 import           Numeric.LinearAlgebra.HMatrix (GMatrix, Matrix,
                                                 toLists, (!), (><))
+import Reanimate.Animation
+import Reanimate.Svg
+
 
 type P = V2 Double
 type Trig = (P,P,P)
@@ -88,6 +91,18 @@ data Mesh = Mesh
   -- , meshSteiner :: Vector Int
   , meshTriangles :: Vector RelTrig }
 
+renderMeshPair :: Mesh -> SVG
+renderMeshPair Mesh{..} = withStrokeColor "black" $ mkGroup
+  [ mkLinePathClosed
+    [ (aPx, aPy)
+    , (bPx, bPy)
+    , (cPx, cPy)]
+  | (a,b,c) <- V.toList meshTriangles
+  , let V2 aPx aPy = meshPointsA V.! a
+        V2 bPx bPy = meshPointsA V.! b
+        V2 cPx cPy = meshPointsA V.! c
+  ]
+
 -- applyA (computeA a b) a = b + some_constant_translation
 -- A = Q P_inv
 computeA :: Trig -> Trig -> Matrix Double
@@ -189,8 +204,9 @@ prepare Mesh{..} = Prep
   where
     aOrigin = meshPointsA V.! pivotIdx
     bOrigin = meshPointsB V.! pivotIdx
-    pivotIdx = case V.head meshTriangles of
-        (a,_,_) -> a
+    pivotIdx = 0 
+    -- 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
diff --git a/src/Reanimate/Svg.hs b/src/Reanimate/Svg.hs
index ea094c1..f8c912f 100644
--- a/src/Reanimate/Svg.hs
+++ b/src/Reanimate/Svg.hs
@@ -21,7 +21,7 @@ import           Reanimate.Svg.Unuse
 import qualified Reanimate.Transform          as Transform
 
 lowerTransformations :: Tree -> Tree
-lowerTransformations = worker Transform.identity
+lowerTransformations = worker False Transform.identity
   where
     updLineCmd m cmd =
       case cmd of
@@ -30,17 +30,30 @@ lowerTransformations = worker Transform.identity
         LineBezier ps -> LineBezier $ map (Transform.transformPoint m) ps
         LineEnd p     -> LineEnd $ Transform.transformPoint m p
     updPath m = lineToPath . map (updLineCmd m) . toLineCommands
-    worker m t =
+    updPoint m (Num a,Num b) =
+      case Transform.transformPoint m (V2 a b) of
+        V2 x y -> (Num x, Num y)
+    updPoint _ other = other -- XXX: Can we do better here?
+    worker hasPathified m t =
       let m' = m * Transform.mkMatrix (t^.transform) in
       case t of
         PathTree path -> PathTree $
           path & pathDefinition %~ updPath m'
                & transform .~ Nothing
         GroupTree g -> GroupTree $
-          g & groupChildren %~ map (worker m')
+          g & groupChildren %~ map (worker hasPathified m')
             & transform .~ Nothing
+        LineTree line ->
+          LineTree $
+            line & linePoint1 %~ updPoint m
+                 & linePoint2 %~ updPoint m
         ClipPathTree{} -> t
-        _ -> mkGroup [t] & transform ?~ [ Transform.toTransformation m ]
+        -- If we encounter an unknown node and we've already tried to convert
+        -- to paths, give up and insert an explicit transformation.
+        _ | hasPathified ->
+          mkGroup [t] & transform ?~ [ Transform.toTransformation m ]
+        -- If we haven't tried to pathify, run pathify only once.
+        _ -> worker True m (pathify t)
 
 lowerIds :: Tree -> Tree
 lowerIds = mapTree worker
diff --git a/videos/morph/dcel.hs b/videos/morph/dcel.hs
new file mode 100644
index 0000000..10b85be
--- /dev/null
+++ b/videos/morph/dcel.hs
@@ -0,0 +1,168 @@
+#!/usr/bin/env stack
+-- stack runghc --package reanimate
+{-# LANGUAGE OverloadedStrings #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ApplicativeDo #-}
+module Main where
+
+import           Control.Lens
+import           Control.Monad
+import           Control.Monad.State
+import           Data.Function
+import           Data.List
+import           Text.Printf
+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           Reanimate
+import           Reanimate.Builtin.Documentation
+import           Reanimate.Math.Balloon
+import           Reanimate.Math.Common
+import           Reanimate.Math.Polygon
+import           Reanimate.Math.EarClip
+import           Reanimate.Math.Smooth
+import           Reanimate.Math.Render
+import           Reanimate.Math.Compatible
+import           Reanimate.Morph.Common
+import           Reanimate.Morph.Linear
+import           Reanimate.Morph.Rigid
+import           Reanimate.Debug
+import Reanimate.Math.DCEL as DCEL
+
+-- p1 = centerPolygon $ shape2
+--p1 = pScale 3.5 $ pAtCenter $ pAddPoints (0+2) (pSetOffset shape13 0)
+-- p1 = setOffset (addPoints 2 shape13) 0
+-- p2 = scalePolygon 0.5 $ centerPolygon shape20
+-- p1 = centerPolygon shape2
+--p2 = pScale 3.5 $ pAtCenter $ pAddPoints 0 (pSetOffset shape14 0)
+-- p2 = setOffset shape14 0
+p1 = pCopy p1'
+p2 = pCopy p2'
+(p1', p2') =
+  -- normalizePolygons
+             closestLinearCorrespondence
+  -- leastWork defaultStretchCosts defaultBendCosts
+  (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 4 $ latex "S")
+  (pAtCenter $ unsafeSVGToPolygon 0.01 $ scale 4 $ latex "C")
+
+(p1s,p2s) = unzip (compatiblyTriangulateP p1 p2)
+
+m2 = buildMesh $ polygonsMesh 
+        (map (fmap realToFrac) $ V.toList $ polygonPoints p2)
+        (map (map (fmap realToFrac) . V.toList . polygonPoints) p2s)
+
+m1 = buildMesh $ polygonsMesh 
+        (map (fmap realToFrac) $ V.toList $ polygonPoints p1)
+        (map (map (fmap realToFrac) . V.toList . polygonPoints) p1s)
+
+
+main :: IO ()
+main = reanimate $ sceneAnimation $ do
+  -- newSpriteSVG_ $ mkBackground "black"
+  newSpriteSVG_ $ mkBackgroundPixel rtfdBackgroundColor
+  -- let m = buildMesh $ polygonMesh $ map (fmap realToFrac) $ V.toList $ polygonPoints p2
+  
+  -- newSpriteSVG_ $ withStrokeWidth (defaultStrokeWidth*0.5) $ scale 2 $
+  --   DCEL.renderMesh m
+  s <- newVar 1
+  mVar <- newVar (m1, m2)
+  let V2 centerX centerY = V2 0 0 -- realToFrac <$> pAccess p2 33
+  adjustZ 2 $ newSprite_ $ do
+    ~(m1, m2) <- unVar mVar
+    pure $ mkGroup
+      [ translate 5 3.5 $ scale 0.5 $ renderMeshStats m1
+      , translate 5 0.5 $ scale 0.5 $ renderMeshStats m2 ]
+  newSprite_ $ do
+    ~(m1,m2) <- unVar mVar
+    sc <- unVar s
+    pure $ mkGroup
+      [translate (-4) 0 $ withStrokeWidth (defaultStrokeWidth*1) $ lowerTransformations $ scale sc $
+        translate (negate centerX) (negate centerY) $ mkGroup
+        [ mkGroup []
+        , DCEL.renderMeshColored m1
+        , DCEL.renderMesh (0.05/sc) m1
+        -- , renderMeshEdges m
+        -- , DCEL.renderMeshSimple (0.10/sc) m
+        ]
+      ,translate 2 0 $ withStrokeWidth (defaultStrokeWidth*1) $ lowerTransformations $ scale sc $
+        translate (negate centerX) (negate centerY) $ mkGroup
+        [ mkGroup []
+        , DCEL.renderMeshColored m2
+        , DCEL.renderMesh (0.05/sc) m2
+        -- , renderMeshEdges m
+        -- , DCEL.renderMeshSimple (0.10/sc) m
+        ]]
+  writeVar s 2
+  wait (1/60)
+  let pipeline1 = last . take 20 . iterate 
+        (uncurry delaunayFlip .
+         uncurry splitInternalEdges . 
+         (\(a,b) -> (meshSmoothPosition a, meshSmoothPosition b)))
+      stages = take 30 $ iterate 
+        (pipeline1 .
+         uncurry splitLongestEdge .
+         pipeline1
+         ) (m1,m2)
+  let pipeline = do
+        modifyVar mVar (uncurry delaunayFlip)
+        modifyVar mVar (uncurry splitInternalEdges)
+        modifyVar mVar $ \(a,b) -> (meshSmoothPosition a, meshSmoothPosition b)
+  forM_ stages $ \newM -> do
+    writeVar mVar newM
+    wait (1/60)
+  -- replicateM_ 30 $ do
+  --   pipeline
+  --   wait (1/60)
+
+  -- replicateM_ 200 $ do
+  --   modifyVar mVar (uncurry splitOuterEdges)
+  --   modifyVar mVar (uncurry splitLongestEdge)
+  --   wait (1/60)
+
+  --   replicateM_ 5 $ do
+  --     pipeline
+  --     wait (1/60)
+
+  -- let m' = execState (flipEdge 39 17) m
+  -- writeVar mVar m'
+  -- wait 1
+  -- tweenVar s 5 $ \v -> fromToS v 10
+    
+  -- newSpriteSVG_ $ scale 3 $ 
+  --   mkGroup $ map polygonShape p2s
+  -- newSpriteSVG_ $ scale 0.5 $ withFillColor "red" $
+  --   mkGroup $ map (polygonDots . pScale 6) p2s
+  -- newSpriteSVG_ $ withFillColor "red" $ scale 3 $ polygonNumDots p2
+  -- wait 1
+  return ()
+ where
+   
+
+drawTrigs :: V.Vector (V2 Double) -> V.Vector RelTrig -> SVG
+drawTrigs points trigs = mkGroup
+  [ mkGroup
+      [ withFillOpacity 1
+        $ withStrokeWidth (defaultStrokeWidth * 0)
+        $ withStrokeColor "grey"
+        $ withFillColor "black"
+        $ drawPolygon
+        $ map (points V.!) [a, b, c]
+      | (a, b, c) <- V.toList trigs
+      ]
+  ]
+
+drawPoint :: V2 Double -> SVG
+drawPoint (V2 x y) = translate x y $ mkCircle 0.1
+
+drawPolygon :: [V2 Double] -> SVG
+drawPolygon lst = mkLinePathClosed [ (x, y) | V2 x y <- lst ]
diff --git a/videos/morph/morph.hs b/videos/morph/morph.hs
index ef4cf2d..585dbab 100644
--- a/videos/morph/morph.hs
+++ b/videos/morph/morph.hs
@@ -1,5 +1,5 @@
 #!/usr/bin/env stack
--- stack --resolver lts-15.04 runghc --package reanimate
+-- stack runghc --package reanimate
 module Main where
 
 import           Codec.Picture.Types
diff --git a/videos/morph/vis.hs b/videos/morph/vis.hs
index 30b3d5e..97f60dc 100644
--- a/videos/morph/vis.hs
+++ b/videos/morph/vis.hs
@@ -1,5 +1,5 @@
 #!/usr/bin/env stack
--- stack --resolver lts-15.04 runghc --package reanimate
+-- stack runghc --package reanimate
 {-# LANGUAGE OverloadedStrings #-}
 module Main where