improve outline computation using divMod'

src/splines/Math/Bezier/Stroke.hs

@@ -30,6 +30,8 @@ import Data.Bifunctor
   ( Bifunctor(bimap) )
 import Data.Coerce
   ( coerce )
+import Data.Fixed
+  ( divMod' )
 import Data.Foldable
   ( for_ )
 import Data.Functor.Identity
@@ -67,11 +69,7 @@ import Data.Act
 import Data.Sequence
   ( Seq(..) )
 import qualified Data.Sequence as Seq
-  --( empty, reverse, singleton, index )
+  ( empty, index, length, reverse, singleton, zipWith )
-import Data.Set
-  ( Set )
-import qualified Data.Set as Set
-  ( insert, member, singleton )
 
 -- deepseq
 import Control.DeepSeq
@@ -427,6 +425,9 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline {
               OutlineData
                 ( TwoSided fwdData ( bimap reverseSpline Seq.reverse bwdData ) )
                 cusps
+          trace ( "bwd at t = 0.58: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.58 ) ) ( return () )
+          trace ( "bwd at t = 0.5966724346435021: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.5966724346435021 ) ) ( return () )
+          trace ( "bwd at t = 0.60: " ++ show ( ( snd . outlineFn fwdBwd ) $ ℝ1 0.60 ) ) ( return () )
           outlineData `deepseq` tell outlineData
           lift $ writeSTRef cachedStrokeRef ( Just outlineData )
 
@@ -978,94 +979,84 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
     fwdSol            = findSolFrom fwdOffset
     ( bwdPt, bwdTgt ) = findSolFrom bwdOffset
 
-    n :: Int
-    n = length strokeData
-
     findSolFrom :: Offset -> ( ℝ 2, T ( ℝ 2 ) )
     findSolFrom ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } )
-      = go ( Set.singleton i00 ) i00 ( fromMaybe 0.5 s00 )
+      = go ( fromIntegral i00 + fromMaybe 0.5 s00 )
       where
 
+        go :: Double -> ( ℝ 2, T ( ℝ 2 ) )
+        go is0 =
+          case sol strokeData is0 of
+            ( goodSoln, pt, tgt )
+              | goodSoln && plausibleTangent tgt
+              -> ( pt, tgt )
+              | otherwise
+              -> ( off • path_t, path'_t )
+
         plausibleTangent :: T ( ℝ 2 ) -> Bool
         plausibleTangent tgt = path'_t ^.^ tgt > 0
 
-        go :: Set Int -> Int -> Double -> ( ℝ 2, T ( ℝ 2 ) )
+    sol :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) )
-        go seen i0 s0 =
+    sol f is0 =
-          case sol s0 ( strokeData `Seq.index` i0 ) of
+      let ( good, is ) =
-            ( goodSoln, s, pt, tgt )
+            case newtonRaphson maxIters precision domain ( eqn f ) is0 of
-              | goodSoln && plausibleTangent tgt
+              Nothing  -> ( False, is0 )
-              -> ( pt, tgt )
+              Just is1 -> ( True , is1 )
-              | let ( i', s0' ) = mbNextPoint i0 ( unℝ1 s )
+          ( ds, dcdt ) = finish f is
-              , not ( i' `Set.member` seen )
+      in ( good, ds, dcdt )
-              -> go ( Set.insert i' seen ) i' s0'
-              | otherwise
-              -> ( off • path_t, path'_t )
 
-        mbNextPoint :: Int -> Double -> ( Int, Double )
+    finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) )
-        mbNextPoint i0 s0
+    finish f is =
-          | s0 <= 0.5
+      let (i, s) = fromDomain is in
-          = ( prev i0, 0.9 )
+      case ( f `Seq.index` i ) ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again...
-          | otherwise
+        StrokeDatum
-          = ( next i0, 0.1 )
+          { dstroke
+          , ee = D12 ( ℝ1 _ee ) ( T ( ℝ1 _𝛿E𝛿t ) ) ( T ( ℝ1 ee_s ) )
+          , 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt
+          } ->
+          -- The total derivative dc/dt is computed by dividing by ∂E/∂s,
+          -- so check it isn't zero first. This corresponds to cusps in the envelope.
+          let dcdt
+                | abs ee_s < epsilon
+                , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9
+                = case ( f `Seq.index` i ) ( ℝ1 s' ) of
+                   StrokeDatum { ee = D12 _ _ ( T ( ℝ1 ee_s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' }
+                     -> recip ee_s' *^ 𝛿E𝛿sdcdt'
+                | otherwise
+                = recip ee_s *^ 𝛿E𝛿sdcdt
+          in --trace
+             --  ( unlines
+             --    [ "solveEnvelopeEquations"
+             --    , "      t = " ++ show _t
+             --    , "      s = " ++ show s
+             --    , "      c = " ++ show dstroke
+             --    , "      E = " ++ show _ee
+             --    , "  ∂E/∂t = " ++ show _𝛿E𝛿t
+             --    , "  ∂E/∂s = " ++ show ee_s
+             --    , "  dc/dt = " ++ show dcdt
+             --    ] )
+              ( value @Double @2 @( ℝ 2 ) dstroke, dcdt )
 
-        prev, next :: Int -> Int
+    eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) )
-        prev i
+    eqn fs is =
-          | i == 0
+      let (i, s) = fromDomain is
-          = n - 1
+      in case ( fs `Seq.index` i ) ( ℝ1 s ) of
-          | otherwise
+           StrokeDatum { ee = D12 ee _ ee_s } ->
-          = i - 1
+             coerce ( ee, ee_s )
-        next i
-          | i == n - 1
-          = 0
-          | otherwise
-          = i + 1
-
-    sol :: Double -> ( ℝ 1 -> StrokeDatum 2 () ) -> ( Bool, ℝ 1, ℝ 2, T ( ℝ 2 ) )
-    sol initialGuess f =
-      let (good, s) = case newtonRaphson maxIters precision domain ( eqn f ) initialGuess of
-                        Nothing -> ( False, initialGuess )
-                        Just s0 -> ( True , s0 )
-      in case f ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again...
-          StrokeDatum
-            { dstroke
-            , ee = D12 ( ℝ1 _ee ) ( T ( ℝ1 _𝛿E𝛿t ) ) ( T ( ℝ1 ee_s ) )
-            , 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt
-            } ->
-            -- The total derivative dc/dt is computed by dividing by ∂E/∂s,
-            -- so check it isn't zero first. This corresponds to cusps in the envelope.
-            let dcdt
-                  | abs ee_s < epsilon
-                  , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9
-                  = case f ( ℝ1 s' ) of
-                     StrokeDatum { ee = D12 _ _ ( T ( ℝ1 ee_s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' }
-                       -> recip ee_s' *^ 𝛿E𝛿sdcdt'
-                  | otherwise
-                  = recip ee_s *^ 𝛿E𝛿sdcdt
-            in --trace
-               --  ( unlines
-               --    [ "solveEnvelopeEquations"
-               --    , "      t = " ++ show _t
-               --    , "      s = " ++ show s
-               --    , "      c = " ++ show dstroke
-               --    , "      E = " ++ show _ee
-               --    , "  ∂E/∂t = " ++ show _𝛿E𝛿t
-               --    , "  ∂E/∂s = " ++ show ee_s
-               --    , "  dc/dt = " ++ show dcdt
-               --    ] )
-                ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke, dcdt )
-
-    eqn :: ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) )
-    eqn f s =
-      case f ( ℝ1 s ) of
-        StrokeDatum { ee = D12 ee _ ee_s } ->
-          coerce ( ee, ee_s )
 
     maxIters :: Word
-    maxIters = 5 --30
+    maxIters = 20
     precision :: Int
-    precision = 10
+    precision = 8
+    n :: Int
+    n = Seq.length strokeData
     domain :: ( Double, Double )
-    domain = ( 0, 1 )
+    domain = ( 0, fromIntegral n )
+
+    fromDomain :: Double -> ( Int, Double )
+    fromDomain is =
+      let ( i0, s ) = is `divMod'` 1
+      in ( i0 `mod` n, s )
 
 newtype ZipSeq a = ZipSeq { getZipSeq :: Seq a }
   deriving stock Functor

src/splines/Math/Linear.hs

@@ -25,8 +25,6 @@ import Data.Kind
   ( Type )
 import GHC.Generics
   ( Generic, Generic1, Generically(..), Generically1(..) )
-import GHC.Show
-  ( showSpace )
 import GHC.TypeNats
   ( Nat, type (+) )
 import Unsafe.Coerce