fix stupid error in quadratic solver

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 =
+          mbParam = listToMaybe
-            case mapMaybe correctTangentParam $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 ) of
+                  . mapMaybe correctTangentParam
-              ( t : _ ) -> Just t
+                  $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 )
-              _         -> between ori tgt0 tgt2 tgt_wanted -- fallback in case we couldn't solve the quadratic for some reason
         in for_ mbParam \ t ->
             throwE $
               Offset

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
+        else 0.5 * ( - a1 + sqrt disc ) / a2
-    in [ r, -r - a1 ]
+    in [ r, -r - a1 / a2 ]
   where
     disc :: a  
     disc = a1 * a1 - 4 * a0 * a2