From 136882510340e16d7eb416546ab6b57ec17220de Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 8 Jan 2024 11:31:18 +0100 Subject: [PATCH] improve outline computation using divMod' --- src/splines/Math/Bezier/Stroke.hs | 151 ++++++++++++++---------------- src/splines/Math/Linear.hs | 2 - 2 files changed, 71 insertions(+), 82 deletions(-) diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 77bf8c5..b48b3ee 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/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 ) -import Data.Set - ( Set ) -import qualified Data.Set as Set - ( insert, member, singleton ) + ( empty, index, length, reverse, singleton, zipWith ) -- 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 ) ) - go seen i0 s0 = - case sol s0 ( strokeData `Seq.index` i0 ) of - ( goodSoln, s, pt, tgt ) - | goodSoln && plausibleTangent tgt - -> ( pt, tgt ) - | let ( i', s0' ) = mbNextPoint i0 ( unℝ1 s ) - , not ( i' `Set.member` seen ) - -> go ( Set.insert i' seen ) i' s0' - | otherwise - -> ( off • path_t, path'_t ) + sol :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) ) + sol f is0 = + let ( good, is ) = + case newtonRaphson maxIters precision domain ( eqn f ) is0 of + Nothing -> ( False, is0 ) + Just is1 -> ( True , is1 ) + ( ds, dcdt ) = finish f is + in ( good, ds, dcdt ) - mbNextPoint :: Int -> Double -> ( Int, Double ) - mbNextPoint i0 s0 - | s0 <= 0.5 - = ( prev i0, 0.9 ) - | otherwise - = ( next i0, 0.1 ) + finish :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> Double -> ( ℝ 2, T ( ℝ 2 ) ) + finish f is = + let (i, s) = fromDomain is in + case ( f `Seq.index` i ) ( ℝ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 `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 - prev i - | i == 0 - = n - 1 - | otherwise - = i - 1 - 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 ) + eqn :: Seq ( ℝ 1 -> StrokeDatum 2 () ) -> ( Double -> ( Double, Double ) ) + eqn fs is = + let (i, s) = fromDomain is + in case ( fs `Seq.index` i ) ( ℝ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 diff --git a/src/splines/Math/Linear.hs b/src/splines/Math/Linear.hs index e6a6863..08cb5d9 100644 --- a/src/splines/Math/Linear.hs +++ b/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