From 2a21980ffc9b028414efac787b83c9081b328abd Mon Sep 17 00:00:00 2001 From: sheaf Date: Mon, 29 Apr 2024 19:35:53 +0200 Subject: [PATCH] Fix issues in withTangent & strictlyParallel function --- MetaBrush.cabal | 3 +- brush-strokes/brush-strokes.cabal | 4 ++ brush-strokes/src/lib/Math/Bezier/Stroke.hs | 80 ++++++++++++--------- brush-strokes/src/lib/Math/Module.hs | 8 ++- 4 files changed, 56 insertions(+), 39 deletions(-) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index c73d9ef..41752a3 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -230,7 +230,8 @@ executable MetaBrush , MetaBrush.Time ghc-options: - -threaded -rtsopts + -threaded + -rtsopts build-depends: metabrushes diff --git a/brush-strokes/brush-strokes.cabal b/brush-strokes/brush-strokes.cabal index e63d12b..f38f9a6 100644 --- a/brush-strokes/brush-strokes.cabal +++ b/brush-strokes/brush-strokes.cabal @@ -283,3 +283,7 @@ benchmark cusps default-language: Haskell2010 + + ghc-options: + -threaded + -rtsopts diff --git a/brush-strokes/src/lib/Math/Bezier/Stroke.hs b/brush-strokes/src/lib/Math/Bezier/Stroke.hs index dbccb1a..acfc766 100644 --- a/brush-strokes/src/lib/Math/Bezier/Stroke.hs +++ b/brush-strokes/src/lib/Math/Bezier/Stroke.hs @@ -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 = diff --git a/brush-strokes/src/lib/Math/Module.hs b/brush-strokes/src/lib/Math/Module.hs index 6e97568..fc78230 100644 --- a/brush-strokes/src/lib/Math/Module.hs +++ b/brush-strokes/src/lib/Math/Module.hs @@ -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 @. --