diff --git a/brush-strokes/src/cusps/bench/Bench/Cases.hs b/brush-strokes/src/cusps/bench/Bench/Cases.hs index 8aae46c..5d13e3e 100644 --- a/brush-strokes/src/cusps/bench/Bench/Cases.hs +++ b/brush-strokes/src/cusps/bench/Bench/Cases.hs @@ -69,10 +69,4 @@ trickyCusp2BrushStroke = defaultStartBoxes :: [ Int ] -> [ ( Int, [ Box 2 ] ) ] defaultStartBoxes is = - [ ( i, [ 𝕀ℝ2 ( 𝕀 zero one ) ( 𝕀 zero one ) ] ) | i <- is ] - -zero, one :: Double -zero = 1e-6 -one = 1 - zero -{-# INLINE zero #-} -{-# INLINE one #-} + [ ( i, [ 𝕀ℝ2 ( 𝕀 0 1 ) ( 𝕀 0 1 ) ] ) | i <- is ] diff --git a/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs b/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs index 89ec016..8ae7c95 100644 --- a/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs +++ b/brush-strokes/src/lib/Math/Bezier/Cubic/Fit.hs @@ -1,7 +1,7 @@ {-# LANGUAGE ScopedTypeVariables #-} module Math.Bezier.Cubic.Fit - ( FitParameters(..), FitPoint(..) + ( FitParameters(..), FitPoint(..), FwdBwd(..) , fitSpline, fitPiece ) where @@ -92,11 +92,19 @@ data FitParameters } deriving stock Show +data FwdBwd + = Fwd + | Bwd + deriving stock ( Eq, Show, Generic ) + deriving anyclass NFData data FitPoint = FitPoint - { fitPoint :: !( ℝ 2 ) } + { fitDir :: !FwdBwd + , fitPoint :: !( ℝ 2 ) + } | FitTangent - { fitPoint :: !( ℝ 2 ) + { fitDir :: !FwdBwd + , fitPoint :: !( ℝ 2 ) , fitTangent :: !( T ( ℝ 2 ) ) } | JoinPoint @@ -118,11 +126,12 @@ data FitPoint -- See 'fitPiece' for more information on the fitting process, -- including the meaning of \( \texttt{t_tol} \) and \( \texttt{maxIters} \). fitSpline - :: FitParameters + :: FwdBwd + -> FitParameters -> ( ℝ 1 -> ( ℝ 2, T ( ℝ 2 ) ) ) -- ^ curve \( t \mapsto ( C(t), C'(t) ) \) to fit -> ( ℝ 1, ℝ 1 ) -- ^ interval \( [t_min, t_max] \) -> ( SplinePts Open, Seq FitPoint ) -fitSpline ( FitParameters {..} ) curveFn = go 0 +fitSpline fwdOrBwd ( FitParameters {..} ) curveFn = go 0 where dt :: Double dt = recip ( fromIntegral nbSegments ) @@ -130,13 +139,13 @@ fitSpline ( FitParameters {..} ) curveFn = go 0 :: Int -> ( ℝ 1, ℝ 1 ) -> ( SplinePts Open, Seq FitPoint ) - go subdiv (t_min, t_max) = + go subdiv ( t_min, t_max ) = let p, r :: ℝ 2 tp, tr :: T ( ℝ 2 ) qs :: [ ℝ 2 ] - (p, tp) = curveFn $ T ( ℝ1 0.0001 ) β€’ t_min - (r, tr) = curveFn $ T ( ℝ1 -0.0001 ) β€’ t_max + (p, tp) = curveFn t_min + (r, tr) = curveFn t_max qs = [ fst $ curveFn ( lerp @( T ( ℝ 1 ) ) ( dt * fromIntegral j ) t_min t_max ) | j <- [ 1 .. nbSegments - 1 ] ] in @@ -145,7 +154,7 @@ fitSpline ( FitParameters {..} ) curveFn = go 0 | subdiv >= maxSubdiv || max_sq_error <= dist_tol ^ ( 2 :: Int ) -> -- trace ( unlines [ "fitSpline: piece is OK", "t_min = " ++ show t_min, "start = " ++ show p, "start tgt = " ++ show tp, "t_max = " ++ show t_max, "end = " ++ show r, "end tgt = " ++ show tr ] ) $ - ( openCubicBezierCurveSpline () bez, ( FitTangent p tp :<| Seq.fromList ( map FitPoint qs ) ) :|> FitTangent r tr ) + ( openCubicBezierCurveSpline () bez, ( FitTangent fwdOrBwd p tp :<| Seq.fromList ( map ( FitPoint fwdOrBwd ) qs ) ) :|> FitTangent fwdOrBwd r tr ) | let t_split :: ℝ 1 t_split = ℝ1 $ min ( 1 - dt ) $ max dt t_split_0 diff --git a/brush-strokes/src/lib/Math/Bezier/Spline.hs b/brush-strokes/src/lib/Math/Bezier/Spline.hs index 3715e43..1fcecc2 100644 --- a/brush-strokes/src/lib/Math/Bezier/Spline.hs +++ b/brush-strokes/src/lib/Math/Bezier/Spline.hs @@ -36,7 +36,7 @@ import qualified Data.Bifunctor.Tannen as Biff import Data.Sequence ( Seq(..) ) import qualified Data.Sequence as Seq - ( singleton, drop, splitAt ) + ( drop, length, singleton, splitAt ) -- deepseq import Control.DeepSeq @@ -196,6 +196,10 @@ data instance Curves Closed crvData ptData deriving stock ( Show, Generic, Generic1, Functor, Foldable, Traversable ) deriving anyclass NFData +nbCurves :: Curves Closed crvData ptData -> Int +nbCurves NoCurves = 0 +nbCurves ( ClosedCurves { prevOpenCurves } ) = 1 + Seq.length prevOpenCurves + instance Bifunctor ( Curves Closed ) where bimap _ _ NoCurves = NoCurves bimap f g ( ClosedCurves p l ) = ClosedCurves ( fmap ( bimap f g ) p ) ( bimap f g l ) diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index 1f31184..4a87d0c 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -113,7 +113,9 @@ import Calligraphy.Brushes import Math.Algebra.Dual import qualified Math.Bezier.Cubic as Cubic import Math.Bezier.Cubic.Fit - ( FitPoint (..), FitParameters, fitSpline ) + ( FwdBwd(..), FitPoint (..) + , FitParameters, fitSpline + ) import Math.Bezier.Spline import qualified Math.Bezier.Quadratic as Quadratic import Math.Bezier.Stroke.EnvelopeEquation @@ -136,6 +138,8 @@ import Math.Orientation import Math.Roots import Math.Root.Isolation +import Debug.Trace + -------------------------------------------------------------------------------- data Offset @@ -273,8 +277,8 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru endPt :: ptData endPt = openCurveEnd lastCurve startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( ℝ 2 ) - ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ ℝ1 1e-4 - ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ ℝ1 (1 - 1e-4) + ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ ℝ1 0 + ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ ℝ1 1 startBrush, endBrush :: SplinePts Closed startBrush = brushShape spt0 endBrush = brushShape endPt @@ -351,8 +355,8 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru endTgt = case prevCurves of Empty -> endTangent spt0 spt0 lastCurve _ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve - ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ ℝ1 1e-9 - ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ ℝ1 (1 - 1e-9) + ( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ ℝ1 0 + ( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ ℝ1 1 fwdStartCap, bwdStartCap :: SplinePts Open OutlineData ( fmap fst -> TwoSided fwdStartCap bwdStartCap ) _ = snd . runWriter @@ -416,8 +420,8 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru tgt, next_tgt, tgtFwd, next_tgtFwd, tgtBwd, next_tgtBwd :: T ( ℝ 2 ) tgt = startTangent spt0 ptData curve next_tgt = endTangent spt0 ptData curve - ( ( _, tgtFwd ), ( _, tgtBwd ) ) = outlineFn fwdBwd $ ℝ1 1e-9 - ( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = outlineFn fwdBwd $ ℝ1 (1 - 1e-9) + ( ( _, tgtFwd ), ( _, tgtBwd ) ) = outlineFn fwdBwd $ ℝ1 0 + ( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = outlineFn fwdBwd $ ℝ1 1 lift $ tellBrushJoin ( prevTgt, prev_tgtFwd, tgtBwd ) ptData ( tgt, tgtFwd, prev_tgtBwd ) lift $ updateCurveData ( curveData curve ) fwdBwd put ( next_tgt, next_tgtFwd, next_tgtBwd ) @@ -447,8 +451,9 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru fwdData, bwdData :: ( SplinePts Open, Seq FitPoint ) ( fwdData, bwdData ) = - ( sconcat $ fmap ( fitSpline fitParams ( fst . outlineFn fwdBwd ) ) intervals - , sconcat $ fmap ( fitSpline fitParams ( snd . outlineFn fwdBwd ) ) intervals ) + ( sconcat $ fmap ( fitSpline Fwd fitParams ( fst . outlineFn fwdBwd ) ) intervals + , sconcat $ fmap ( fitSpline Bwd fitParams ( snd . outlineFn fwdBwd ) ) intervals + ) -- TODO: use foldMap1 once that's in base. `Strats.using` ( Strats.parTuple2 Strats.rdeepseq Strats.rdeepseq ) @@ -899,7 +904,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) -- only allow non-empty splines | Just tgt_last <- lastTangent spline -- only allow well-defined query tangent vectors - , not (badTangent tgt_wanted) + , not $ badTangent tgt_wanted = case runExcept . ( `runStateT` tgt_last ) $ ibifoldSpline go ( \ _ -> pure () ) $ adjustSplineType @Open spline of Left off -> off @@ -931,13 +936,31 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) tgt_end = endTangent splineStart cp cseg -- Handle corner. unless ( tgt_prev `strictlyParallel` tgt_start ) do - for_ ( between ori tgt_prev tgt_start tgt_wanted ) \ _ -> - lift . throwE $ - Offset + for_ ( between ori tgt_prev tgt_start tgt_wanted ) \ a -> do + let + -- Decide whether the offset should be at the start of this + -- segment or at the end of the previous segment, depending + -- on which tangent is closer. + -- + -- This is important to get the best initial guess possible + -- for Newton's method when solving the envelope equation. + off + | a < 0.5 + = Offset + { offsetIndex = + if i == 0 + then nbCurves ( splineCurves spline ) - 1 + else i - 1 + , offsetParameter = Just 1 + , offset = T p + } + | otherwise + = Offset { offsetIndex = i , offsetParameter = Just 0 , offset = T p } + lift $ throwE $ off -- Handle segment. lift $ handleSegment i p seg tgt_start put tgt_end @@ -1215,11 +1238,6 @@ data RootSolvingAlgorithm -- | Use the modified Halley M2 method. | HalleyM2 { maxIters :: Word, precision :: Int } -data FwdBwd - = Fwd - | Bwd - deriving stock ( Eq, Show ) - -- | Solve the envelope equations at a given point \( t = t_0 \), to find -- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke. solveEnvelopeEquations :: RootSolvingAlgorithm @@ -1243,8 +1261,10 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse mbBwdCornerSol = getFirst $ foldMap ( First . cornerSol Bwd ) corners findSolFrom :: FwdBwd -> Offset -> ( ℝ 2, T ( ℝ 2 ) ) - findSolFrom fwdOrBwd ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } ) - = go ( fromIntegral i00 + fromMaybe 0.5 s00 ) + findSolFrom fwdOrBwd ( Offset { offsetIndex = i00, offsetParameter = s00, offset = off } ) = + go ( fromIntegral i00 + maybe 0.5 ( min 0.99999999999999989 ) s00 ) + -- NB: the 'fromDomain' function requires values in the + -- half-open interval [0,1[, hence the @min 0.99999999999999989@ above. where go :: Double -> ( ℝ 2, T ( ℝ 2 ) ) @@ -1254,7 +1274,20 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse | goodSoln -> ( pt, tgt ) | otherwise - -> ( off β€’ path_t, path'_t ) + -> trace + ( unlines + [ "solveEnvelopeEquations: bad solution" + , "falling back to naive guess; this will lead to inaccurate fitting" + , " t: " ++ show _t + , " dir: " ++ show fwdOrBwd + , " p: " ++ show path_t + , " p': " ++ show path'_t + , " is0: " ++ show is0 + , " BAD pt: " ++ show pt + , " BAD tgt: " ++ show tgt + ] + ) + ( off β€’ path_t, path'_t ) sol :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Bool, ℝ 2, T ( ℝ 2 ) ) sol fwdOrBwd f is0 = @@ -1267,8 +1300,8 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse Bwd -> sgn <= 0 , is1 ) - ( sgn, ds, dcdt ) = finish f is - in ( good, ds, dcdt ) + ( sgn, c, dcdt ) = finish fwdOrBwd f is + in ( good, c, dcdt ) runSolveMethod = case rootAlgo of HalleyM2 { maxIters, precision } -> @@ -1276,8 +1309,8 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse NewtonRaphson { maxIters, precision } -> newtonRaphson maxIters precision domain - finish :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) ) - finish fs is = + finish :: FwdBwd -> Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> Double -> ( Double, ℝ 2, T ( ℝ 2 ) ) + finish _fwdOrBwd fs is = let (i, s) = fromDomain is in case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again... StrokeDatum @@ -1292,7 +1325,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse -- so check it isn't zero first. This corresponds to cusps in the envelope. let unrot_dcdt | abs 𝛿E𝛿s < epsilon - , let s' = if s >= 0.5 then s - 1e-6 else s + 1e-6 + , let s' = if s >= 0.5 then s - 1e-7 else s + 1e-7 = case ( fs `Seq.index` i ) ( ℝ1 s' ) of StrokeDatum { ee = D12 _ _ ( T ( ℝ1 𝛿E𝛿s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' } -> recip 𝛿E𝛿s' *^ 𝛿E𝛿sdcdt' @@ -1313,7 +1346,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse evalStrokeDatum :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> ( Double -> StrokeDatum 2 NonIV ) evalStrokeDatum fs is = let (i, s) = fromDomain is - in ( fs `Seq.index` i ) ( ℝ1 $ max 1e-6 $ min (1 - 1e-6) $ s ) + in ( fs `Seq.index` i ) ( ℝ1 $ max 0 $ min 1 $ s ) eqn :: Seq ( ℝ 1 -> StrokeDatum 2 NonIV ) -> ( Double -> (# Double, Double #) ) eqn fs is = @@ -1336,7 +1369,6 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse Fwd -> id Bwd -> ( -1 *^ ) b_s = flipWhenBwd ( T u ) - ori = CCW --case fwdOrBwd of { Fwd -> CCW ; Bwd -> CW } res@( _pt, _tgt ) = case mbRotation of Nothing -> ( unT $ T p ^+^ T b , T u ) @@ -1346,22 +1378,7 @@ solveEnvelopeEquations rootAlgo ( ℝ1 _t ) path_t path'_t ( fwdOffset, bwdOffse in ( unT $ T p ^+^ rotate cosΞΈ sinΞΈ ( T b ) , rotate cosΞΈ sinΞΈ ( T u ) ) in -{- - trace (unlines - [ "cornerSol" - , "fwdOrBwd: " ++ show fwdOrBwd - , "ori: " ++ show ori - , "dp: " ++ show _dp - , "startTgt: " ++ show startTgt - , "endTgt: " ++ show endTgt - , "u: " ++ show u - , "b_s: " ++ show b_s - , "ok: " ++ show ( between ori startTgt endTgt b_s ) - , "pt: " ++ show _pt - , "tgt: " ++ show _tgt - ]) $ --} - if isNothing $ between ori startTgt endTgt b_s + if isNothing $ between CCW startTgt endTgt b_s then Nothing else Just res @@ -1441,13 +1458,8 @@ findCusps ( opts1, opts2 ) boxStrokeData = ) where unit :: 𝕀 - unit = 𝕀 zero one + unit = 𝕀 0 1 {-# INLINE unit #-} - zero, one :: Double - zero = 1e-6 - one = 1 - zero - {-# INLINE zero #-} - {-# INLINE one #-} -- | Like 'findCusps', but parametrised over the initial boxes for the -- root isolation method. diff --git a/src/app/MetaBrush/Application.hs b/src/app/MetaBrush/Application.hs index bbea452..4f75912 100644 --- a/src/app/MetaBrush/Application.hs +++ b/src/app/MetaBrush/Application.hs @@ -221,7 +221,7 @@ runApplication application = do maxHistorySizeTVar <- STM.newTVarIO 1000 fitParametersTVar <- STM.newTVarIO $ FitParameters - { maxSubdiv = 2 --5 --2 --3 -- 6 + { maxSubdiv = 1 --2 --5 --2 --3 -- 6 , nbSegments = 3 , dist_tol = 5e-3 , t_tol = 1e-4 diff --git a/src/app/MetaBrush/Render/Document.hs b/src/app/MetaBrush/Render/Document.hs index a6c2792..49ca653 100644 --- a/src/app/MetaBrush/Render/Document.hs +++ b/src/app/MetaBrush/Render/Document.hs @@ -66,7 +66,7 @@ import Math.Algebra.Dual import qualified Math.Bezier.Cubic as Cubic ( Bezier(..), fromQuadratic ) import Math.Bezier.Cubic.Fit - ( FitPoint(..), FitParameters ) + ( FitPoint(..), FitParameters, FwdBwd(..) ) import qualified Math.Bezier.Quadratic as Quadratic ( Bezier(..) ) import Math.Bezier.Spline @@ -691,7 +691,7 @@ drawDebugInfo cols zoom@( Zoom { zoomFactor } ) ( fitPts, cusps ) = do for_ cusps ( drawCusp cols zoom ) drawFitPoint :: Colours -> Zoom -> FitPoint -> StateT Double Cairo.Render () -drawFitPoint _ ( Zoom { zoomFactor } ) ( FitPoint { fitPoint = ℝ2 x y } ) = do +drawFitPoint _ zoom ( FitPoint { fitDir = fwdOrBwd, fitPoint = ℝ2 x y } ) = do hue <- get put ( hue + 0.01 ) @@ -701,9 +701,7 @@ drawFitPoint _ ( Zoom { zoomFactor } ) ( FitPoint { fitPoint = ℝ2 x y } ) = do lift do Cairo.save Cairo.translate x y - Cairo.arc 0 0 ( 4 / zoomFactor ) 0 ( 2 * pi ) - Cairo.setSourceRGBA r g b 0.8 - Cairo.fill + drawFitPointHelper zoom fwdOrBwd r g b Cairo.restore drawFitPoint _ ( Zoom { zoomFactor } ) ( JoinPoint { joinPoint = ℝ2 x y } ) = lift do @@ -717,7 +715,7 @@ drawFitPoint _ ( Zoom { zoomFactor } ) ( JoinPoint { joinPoint = ℝ2 x y } ) = Cairo.stroke Cairo.restore -drawFitPoint _ ( Zoom { zoomFactor } ) ( FitTangent { fitPoint = ℝ2 x y, fitTangent = V2 tx ty } ) = do +drawFitPoint _ zoom@( Zoom { zoomFactor } ) ( FitTangent { fitDir = fwdOrBwd, fitPoint = ℝ2 x y, fitTangent = V2 tx ty } ) = do hue <- get put ( hue + 0.01 ) @@ -728,14 +726,30 @@ drawFitPoint _ ( Zoom { zoomFactor } ) ( FitTangent { fitPoint = ℝ2 x y, fitT Cairo.save Cairo.translate x y Cairo.moveTo 0 0 - Cairo.lineTo ( 0.05 * tx ) ( 0.05 * ty ) + Cairo.lineTo ( min 0.05 ( 1 / zoomFactor ) * tx ) ( min 0.05 ( 1 / zoomFactor ) * ty ) Cairo.setLineWidth ( 4 / zoomFactor ) - Cairo.setSourceRGBA r g b 1.0 + Cairo.setSourceRGBA r g b 0.8 Cairo.stroke - Cairo.arc 0 0 ( 2 / zoomFactor ) 0 ( 2 * pi ) - Cairo.fill + drawFitPointHelper zoom fwdOrBwd r g b Cairo.restore +drawFitPointHelper :: Zoom -> FwdBwd -> Double -> Double -> Double -> Cairo.Render () +drawFitPointHelper ( Zoom { zoomFactor } ) fwdOrBwd r g b = do + Cairo.moveTo 0 0 + Cairo.arc 0 0 ( 4 / zoomFactor ) 0 ( 2 * pi ) + Cairo.setSourceRGBA r g b 0.8 + Cairo.fill + Cairo.moveTo 0 0 + case fwdOrBwd of + Fwd -> do + Cairo.arc 0 ( -1 / zoomFactor ) ( 1 / zoomFactor ) 0 ( 2 * pi ) + Cairo.setSourceRGBA 255 255 255 0.8 + Cairo.fill + Bwd -> do + Cairo.arc 0 ( 1 / zoomFactor ) ( 1 / zoomFactor ) 0 ( 2 * pi ) + Cairo.setSourceRGBA 0 0 0 0.8 + Cairo.fill + drawCusp :: Colours -> Zoom -> Cusp -> Cairo.Render () drawCusp _ ( Zoom { zoomFactor } ) ( Cusp { cuspPathCoords = D21 { _D21_v = ℝ2 px py