From 570a593f4f59147c38111ce8d8ac934006224db1 Mon Sep 17 00:00:00 2001 From: sheaf Date: Wed, 12 May 2021 15:16:26 +0200 Subject: [PATCH] fix stupid error in quadratic solver --- src/lib/Math/Bezier/Stroke.hs | 50 ++++++++++++++++++++++------------- src/lib/Math/Roots.hs | 4 +-- 2 files changed, 34 insertions(+), 20 deletions(-) diff --git a/src/lib/Math/Bezier/Stroke.hs b/src/lib/Math/Bezier/Stroke.hs index 06db791..068649f 100644 --- a/src/lib/Math/Bezier/Stroke.hs +++ b/src/lib/Math/Bezier/Stroke.hs @@ -41,7 +41,7 @@ import Data.Foldable import Data.List.NonEmpty ( unzip ) import Data.Maybe - ( fromMaybe, isJust, mapMaybe ) + ( fromMaybe, isJust, listToMaybe, mapMaybe ) import GHC.Exts ( newMutVar#, runRW# ) import GHC.STRef @@ -212,7 +212,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = endTgtBwd = (-1) *^ snd ( lastOutlineBwd 0 ) startBrush, endBrush :: SplinePts Closed startBrush = brushShape spt0 - endBrush = brushShape endPt + endBrush = brushShape endPt -- Computation of which brush segment to use for the end caps. startTgt, endTgt :: Vector2D Double @@ -229,24 +229,24 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = stx, sty, etx, ety :: Double Vector2D stx sty = startTgt Vector2D etx ety = endTgt - + startCap, endCap :: SplinePts Open startCap | isJust $ between brushOrientation startTgtBwd startTgtFwd startTestTgt = fmap ( MkVector2D ( coords spt0 ) • ) - $ joinWithBrush ( withTangent startTgtBwd startBrush ) ( withTangent startTgtFwd startBrush ) startBrush + $ joinWithBrush startBrush startTgtBwd startTgtFwd | otherwise = fmap ( MkVector2D ( coords spt0 ) • ) . reverseSpline - $ joinWithBrush ( withTangent startTgtFwd startBrush ) ( withTangent startTgtBwd startBrush )startBrush + $ joinWithBrush startBrush startTgtFwd startTgtBwd endCap | not . isJust $ between brushOrientation endTgtBwd endTgtFwd endTestTgt = fmap ( MkVector2D ( coords endPt ) • ) . reverseSpline - $ joinWithBrush ( withTangent endTgtBwd endBrush ) ( withTangent endTgtFwd endBrush ) endBrush + $ joinWithBrush endBrush endTgtBwd endTgtFwd | otherwise = fmap ( MkVector2D ( coords endPt ) • ) - $ joinWithBrush ( withTangent endTgtFwd endBrush ) ( withTangent endTgtBwd endBrush ) endBrush + $ joinWithBrush endBrush endTgtFwd endTgtBwd -> do TwoSided ( fwdPts, fwdFits ) ( bwdPts, bwdFits ) <- updateSpline ( startTgtFwd, startTgtBwd ) pure @@ -373,21 +373,21 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty } | brushOrientation == convexOrientation [ prevTgtFwd, tgtFwd ] = fmap ( ptOffset • ) - $ joinWithBrush ( withTangent prevTgtFwd brush0 ) ( withTangent tgtFwd brush0 ) brush0 + $ joinWithBrush brush0 prevTgtFwd tgtFwd | otherwise = fmap ( ptOffset • ) . reverseSpline - $ joinWithBrush ( withTangent tgtFwd brush0 ) ( withTangent prevTgtFwd brush0 ) brush0 + $ joinWithBrush brush0 tgtFwd prevTgtFwd bwdJoin | tgtBwd `strictlyParallel` prevTgtBwd = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty } | brushOrientation == convexOrientation [ tgtBwd, prevTgtBwd ] = fmap ( ptOffset • ) - $ joinWithBrush ( withTangent tgtBwd brush0 ) ( withTangent prevTgtBwd brush0 ) brush0 + $ joinWithBrush brush0 tgtBwd prevTgtBwd | otherwise = fmap ( ptOffset • ) . reverseSpline - $ joinWithBrush ( withTangent prevTgtBwd brush0 ) ( withTangent tgtBwd brush0 ) brush0 + $ joinWithBrush brush0 prevTgtBwd tgtBwd brushJoin :: OutlineData brushJoin = TwoSided ( fwdJoin, Empty ) ( bwdJoin, Empty ) @@ -499,7 +499,6 @@ outlineFunctions ptParams brushFn sp0 crv = -- used in the "stroke" function. ----- - startTangent, endTangent :: ( SplineTypeI clo, HasType ( Point2D Double ) ptData ) => ptData -> ptData -> Curve clo crvData ptData -> Vector2D Double startTangent sp p0 ( LineTo mp1 _ ) = coords p0 --> coords ( fromNextPoint sp mp1 ) startTangent _ p0 ( Bezier2To p1 _ _ ) = coords p0 --> coords p1 @@ -516,11 +515,27 @@ lastTangent ( Spline { splineStart, splineCurves = ClosedCurves ( _ :|> prev ) l -------------------------------------------------------------------------------- -- | Compute the join at a point of discontinuity of the tangent vector direction (G1 discontinuity). -joinWithBrush :: forall crvData ptData. HasType ( Point2D Double ) ptData => Offset -> Offset -> Spline Closed crvData ptData -> SplinePts Open joinWithBrush + :: ( HasType ( Point2D Double ) ptData + -- debugging + , Show ptData, Show crvData + ) + => Spline Closed crvData ptData + -> Vector2D Double + -> Vector2D Double + -> SplinePts Open +joinWithBrush brush startTgt endTgt = joinBetweenOffsets brush startOffset endOffset + where + startOffset, endOffset :: Offset + startOffset = withTangent startTgt brush + endOffset = withTangent endTgt brush + +-- | Select the section of a spline in between two offsets. +joinBetweenOffsets :: forall crvData ptData. HasType ( Point2D Double ) ptData => Spline Closed crvData ptData -> Offset -> Offset -> SplinePts Open +joinBetweenOffsets + spline ( Offset { offsetIndex = i1, offsetParameter = mb_t1 } ) ( Offset { offsetIndex = i2, offsetParameter = mb_t2 } ) - spline | i2 > i1 = let pcs, lastAndRest :: Maybe ( SplinePts Open ) @@ -746,10 +761,9 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli in let mbParam :: Maybe Double - mbParam = - case mapMaybe correctTangentParam $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 ) of - ( t : _ ) -> Just t - _ -> between ori tgt0 tgt2 tgt_wanted -- fallback in case we couldn't solve the quadratic for some reason + mbParam = listToMaybe + . mapMaybe correctTangentParam + $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 ) in for_ mbParam \ t -> throwE $ Offset diff --git a/src/lib/Math/Roots.hs b/src/lib/Math/Roots.hs index e26894c..232f726 100644 --- a/src/lib/Math/Roots.hs +++ b/src/lib/Math/Roots.hs @@ -69,8 +69,8 @@ solveQuadratic a0 a1 a2 r = if a1 >= 0 then 2 * a0 / ( - a1 - sqrt disc ) - else 0.5 * ( - a1 + sqrt disc) / a2 - in [ r, -r - a1 ] + else 0.5 * ( - a1 + sqrt disc ) / a2 + in [ r, -r - a1 / a2 ] where disc :: a disc = a1 * a1 - 4 * a0 * a2