improve corner initial guess for Newton iteration

This commit improves the initial guess we use for the Newton iteration
to solve the envelope equation when the best tangent we found arises
from a corner: instead of starting at (i = i0, s = 0), i.e. at the start
of the current curve in the spline, we might instead start at
(i = i0 - 1, s = 1).

This commit also improves the debug output to make it clearer which part
of the path goes forwards (little white dot above, inside fit points)
or backwards (little black dot below, inside fit points).
This commit is contained in:
sheaf 2024-11-17 15:34:55 +01:00
parent aeb7a425f2
commit 13cc649319
6 changed files with 111 additions and 78 deletions

View file

@ -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 ]

View file

@ -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

View file

@ -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 )

View file

@ -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.

View file

@ -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

View file

@ -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