Fix issues in withTangent & strictlyParallel function

This commit is contained in:
sheaf 2024-04-29 19:35:53 +02:00
parent 63b9703faf
commit 2a21980ffc
4 changed files with 56 additions and 39 deletions

View file

@ -230,7 +230,8 @@ executable MetaBrush
, MetaBrush.Time
ghc-options:
-threaded -rtsopts
-threaded
-rtsopts
build-depends:
metabrushes

View file

@ -283,3 +283,7 @@ benchmark cusps
default-language:
Haskell2010
ghc-options:
-threaded
-rtsopts

View file

@ -56,8 +56,7 @@ import Data.Proxy
import Data.Semigroup
( sconcat )
import GHC.Exts
( newMutVar#, runRW#, inline
)
( newMutVar#, runRW# )
import GHC.STRef
( STRef(..), readSTRef, writeSTRef )
import GHC.Generics
@ -145,8 +144,6 @@ import Math.Orientation
import Math.Roots
import Math.Root.Isolation
--import Debug.Utils
--------------------------------------------------------------------------------
data Offset
@ -278,8 +275,8 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
endPt :: ptData
endPt = openCurveEnd lastCurve
startTgtFwd, startTgtBwd, endTgtFwd, endTgtBwd :: T ( 2 )
( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ 1 0
( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ 1 1
( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ 1 1e-4
( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ 1 (1 - 1e-4)
startBrush, endBrush :: SplinePts Closed
startBrush = brushShape spt0
endBrush = brushShape endPt
@ -342,8 +339,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 0
( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ 1 1
( ( _, startTgtFwd), ( _, startTgtBwd ) ) = outlineFn firstOutlineFn $ 1 1e-9
( ( _, endTgtFwd ), ( _, endTgtBwd ) ) = outlineFn lastOutlineFn $ 1 (1 - 1e-9)
fwdStartCap, bwdStartCap :: SplinePts Open
OutlineData ( fmap fst -> TwoSided fwdStartCap bwdStartCap ) _
= snd . runWriter
@ -368,7 +365,7 @@ computeStrokeOutline rootAlgo mbCuspOptions fitParams ptParams toBrushParams bru
where
outlineInfo :: ptData -> Curve Open crvData ptData -> OutlineInfo
outlineInfo = inline ( outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brush )
outlineInfo = outlineFunction rootAlgo mbCuspOptions ptParams toBrushParams brush
outlineFns :: Seq OutlineInfo
outlineFns = go spt0 ( openCurves $ splineCurves ( adjustSplineType @Open spline ) )
@ -407,8 +404,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 0
( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = outlineFn fwdBwd $ 1 1
( ( _, tgtFwd ), ( _, tgtBwd ) ) = outlineFn fwdBwd $ 1 1e-9
( ( _, next_tgtFwd ), ( _, next_tgtBwd ) ) = outlineFn fwdBwd $ 1 (1 - 1e-9)
lift $ tellBrushJoin ( prevTgt, prev_tgtFwd, tgtBwd ) ptData ( tgt, tgtFwd, prev_tgtBwd )
lift $ updateCurveData ( curveData curve ) fwdBwd
put ( next_tgt, next_tgtFwd, next_tgtBwd )
@ -857,7 +854,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
-- only allow well-defined query tangent vectors
, not (badTangent tgt_wanted)
= case runExcept . ( `runStateT` tgt_last ) $ ibifoldSpline go ( \ _ -> pure () ) $ adjustSplineType @Open spline of
Left off -> off
Left off ->
off
Right _ ->
error $
"withTangent: could not find any point with given tangent vector\n\
@ -910,12 +908,12 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
let
tgt1 :: T ( 2 )
tgt1 = p1 --> p2
in for_ ( convexCombination tgt0 tgt1 tgt_wanted ) \ t ->
in for_ ( convexCombination tgt0 tgt1 tgt_wanted ) \ s ->
throwE $
Offset
{ offsetIndex = i
, offsetParameter = Just t
, offset = T $ Quadratic.bezier @( T ( 2 ) ) ( Quadratic.Bezier {..} ) t
, offsetParameter = Just s
, offset = T $ Quadratic.bezier @( T ( 2 ) ) ( Quadratic.Bezier {..} ) s
}
handleSegment i p0 ( Bezier3To p1 p2 ( NextPoint p3 ) _ ) tgt0 =
let
@ -929,10 +927,10 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
c12 = tgt_wanted × tgt1
c23 = tgt_wanted × tgt2
correctTangentParam :: Double -> Maybe Double
correctTangentParam t
| t > -epsilon && t < 1 + epsilon
, tgt_wanted ^.^ Cubic.bezier' bez t > epsilon
= Just ( max 0 ( min 1 t ) )
correctTangentParam s
| s > -epsilon && s < 1 + epsilon
, tgt_wanted ^.^ Cubic.bezier' bez s > epsilon
= Just ( max 0 ( min 1 s ) )
| otherwise
= Nothing
in
@ -941,12 +939,12 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
mbParam = listToMaybe
. mapMaybe correctTangentParam
$ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 )
in for_ mbParam \ t ->
in for_ mbParam \ s ->
throwE $
Offset
{ offsetIndex = i
, offsetParameter = Just t
, offset = T $ Cubic.bezier @( T ( 2 ) ) bez t
, offsetParameter = Just s
, offset = T $ Cubic.bezier @( T ( 2 ) ) bez s
}
--------------------------------------------------------------------------------
@ -1101,7 +1099,7 @@ solveEnvelopeEquations :: RootSolvingAlgorithm
-> ( Offset, Offset )
-> Seq ( 1 -> StrokeDatum 2 () )
-> ( ( 2, T ( 2 ) ), ( 2, T ( 2 ) ) )
solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData
solveEnvelopeEquations rootAlgo ( 1 _t ) path_t path'_t ( fwdOffset, bwdOffset ) strokeData
= ( fwdSol, ( bwdPt, -1 *^ bwdTgt ) )
where
@ -1118,23 +1116,23 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
go is0 =
case sol desc strokeData is0 of
( goodSoln, pt, tgt )
| goodSoln && plausibleTangent tgt
| goodSoln
-> ( pt, tgt )
| otherwise
-> ( off path_t, path'_t )
plausibleTangent :: T ( 2 ) -> Bool
plausibleTangent tgt = path'_t ^.^ tgt > 0
sol :: String -> Seq ( 1 -> StrokeDatum 2 () ) -> Double -> ( Bool, 2, T ( 2 ) )
sol _desc f is0 =
let (solRes, _solSteps) = runSolveMethod ( eqn f ) is0
sol desc f is0 =
let ( solRes, _solSteps ) = runSolveMethod ( eqn f ) is0
( good, is ) =
case solRes of
Nothing -> ( False, is0 )
Just is1 -> ( True , is1 )
( ds, dcdt ) = finish f is
Just is1 -> ( if desc == "fwd"
then sgn >= 0
else sgn <= 0
, is1
)
( sgn, ds, dcdt ) = finish f is
in ( good, ds, dcdt )
runSolveMethod = case rootAlgo of
@ -1143,18 +1141,21 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
NewtonRaphson { maxIters, precision } ->
newtonRaphson maxIters precision domain
finish :: Seq ( 1 -> StrokeDatum 2 () ) -> Double -> ( 2, T ( 2 ) )
finish :: Seq ( 1 -> StrokeDatum 2 () ) -> Double -> ( Double, 2, T ( 2 ) )
finish fs is =
let (i, s) = fromDomain is in
case evalStrokeDatum fs is of -- TODO: a bit redundant to have to compute this again...
StrokeDatum
{ stroke
, mbRotation
, ee = D12 ( 1 _ee ) _ ( T ( 1 𝛿E𝛿s ) )
, du = D12 u _ _
, dv = D12 v _ _
, 𝛿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
let unrot_dcdt
| abs 𝛿E𝛿s < epsilon
, let s' = if s >= 0.5 then s - 1e-6 else s + 1e-6
= case ( fs `Seq.index` i ) ( 1 s' ) of
@ -1162,7 +1163,16 @@ solveEnvelopeEquations rootAlgo _t path_t path'_t ( fwdOffset, bwdOffset ) strok
-> recip 𝛿E𝛿s' *^ 𝛿E𝛿sdcdt'
| otherwise
= recip 𝛿E𝛿s *^ 𝛿E𝛿sdcdt
in ( stroke, dcdt )
dcdt = case mbRotation of
Nothing -> unrot_dcdt
Just ( D21 θ _ _ ) ->
let cosθ = cos θ
sinθ = sin θ
in rotate cosθ sinθ $ unrot_dcdt
in ( T u ^.^ T v, stroke, dcdt )
-- Compute the dot product of u and v (which are rotated versions of ∂c/∂t and ∂c/∂s).
-- The sign of this quantity determines which side of the envelope
-- we are on.
evalStrokeDatum :: Seq ( 1 -> StrokeDatum 2 () ) -> ( Double -> StrokeDatum 2 () )
evalStrokeDatum fs is =

View file

@ -129,11 +129,13 @@ instance Cross Double ( T ( 2 ) ) where
-- | Compute whether two vectors point in the same direction,
-- that is, whether each vector is a (strictly) positive multiple of the other.
--
-- Returns @False@ if either of the vectors is zero.
-- Returns @False@ if either of the vectors is zero (or very close to zero).
strictlyParallel :: T ( 2 ) -> T ( 2 ) -> Bool
strictlyParallel u v
= abs ( u × v ) < epsilon -- vectors are collinear
&& u ^.^ v > epsilon -- vectors point in the same direction (parallel and not anti-parallel)
= abs ( u × v ) < tol -- vectors are collinear
&& u ^.^ v > tol -- vectors point in the same direction (parallel and not anti-parallel)
where
tol = norm u * norm v * epsilon
-- | Finds whether the query vector @ u @ is a convex combination of the two provided vectors @ v0 @, @ v1 @.
--