diff --git a/src/lib/Math/Bezier/Stroke.hs b/src/lib/Math/Bezier/Stroke.hs index 1552065..06db791 100644 --- a/src/lib/Math/Bezier/Stroke.hs +++ b/src/lib/Math/Bezier/Stroke.hs @@ -22,7 +22,6 @@ module Math.Bezier.Stroke , CachedStroke(..), discardCache, invalidateCache , computeStrokeOutline, joinWithBrush , withTangent - , between, parallel ) where @@ -32,7 +31,7 @@ import Prelude import Control.Arrow ( first, (***) ) import Control.Monad - ( guard, unless ) + ( unless ) import Control.Monad.ST ( RealWorld, ST ) import Data.Bifunctor @@ -42,7 +41,7 @@ import Data.Foldable import Data.List.NonEmpty ( unzip ) import Data.Maybe - ( fromMaybe, mapMaybe ) + ( fromMaybe, isJust, mapMaybe ) import GHC.Exts ( newMutVar#, runRW# ) import GHC.STRef @@ -103,17 +102,20 @@ import Math.Bezier.Spline , KnownSplineType ( bifoldSpline, ibifoldSpline ) , Spline(..), SplinePts, Curves(..), Curve(..) - , openCurveEnd, splitSplineAt, dropCurves + , openCurveStart, openCurveEnd, splitSplineAt, dropCurves ) import qualified Math.Bezier.Quadratic as Quadratic import Math.Epsilon ( epsilon ) import Math.Module - ( Module((^-^), (*^)), Inner((^.^)) + ( Module((*^)), Inner((^.^)) , lerp, squaredNorm, cross + , convexCombination, strictlyParallel ) import Math.Orientation - ( Orientation(..), convexOrientation, splineOrientation ) + ( Orientation(..), convexOrientation, splineOrientation + , between + ) import Math.Roots ( solveQuadratic ) import Math.Vector2D @@ -152,7 +154,7 @@ instance Monoid OutlineData where newtype CachedStroke s = CachedStroke { cachedStrokeRef :: STRef s ( Maybe OutlineData ) } instance Show ( CachedStroke s ) where - show _ = "CachedStroke..." + show _ = "<>" instance NFData ( CachedStroke s ) where rnf _ = () @@ -198,7 +200,8 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = SOpen | ( firstOutlineFwd, firstOutlineBwd ) :<| _ <- outlineFns , _ :|> ( lastOutlineFwd, lastOutlineBwd ) <- outlineFns - , _ :|> lastCurve <- openCurves $ splineCurves spline + , firstCurve :<| _ <- openCurves $ splineCurves spline + , prevCurves :|> lastCurve <- openCurves $ splineCurves spline , let endPt :: ptData endPt = openCurveEnd lastCurve @@ -210,9 +213,26 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = startBrush, endBrush :: SplinePts Closed startBrush = brushShape spt0 endBrush = brushShape endPt + + -- Computation of which brush segment to use for the end caps. + startTgt, endTgt :: Vector2D Double + startTgt = coords spt0 --> coords ( openCurveStart firstCurve ) + endTgt = case prevCurves of + Empty -> endTangent spt0 spt0 lastCurve + _ :|> lastPrev -> endTangent spt0 ( openCurveEnd lastPrev ) lastCurve + startTestTgt, endTestTgt :: Vector2D Double + ( startTestTgt, endTestTgt ) = + case brushOrientation of + CCW -> ( Vector2D sty (-stx), Vector2D ety (-etx) ) + CW -> ( Vector2D (-sty) stx , Vector2D (-ety) etx ) + where + stx, sty, etx, ety :: Double + Vector2D stx sty = startTgt + Vector2D etx ety = endTgt + startCap, endCap :: SplinePts Open startCap - | brushOrientation == convexOrientation [ startTgtBwd, startTgtFwd ] + | isJust $ between brushOrientation startTgtBwd startTgtFwd startTestTgt = fmap ( MkVector2D ( coords spt0 ) • ) $ joinWithBrush ( withTangent startTgtBwd startBrush ) ( withTangent startTgtFwd startBrush ) startBrush | otherwise @@ -220,7 +240,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = . reverseSpline $ joinWithBrush ( withTangent startTgtFwd startBrush ) ( withTangent startTgtBwd startBrush )startBrush endCap - | brushOrientation == convexOrientation [ endTgtFwd, endTgtBwd ] + | not . isJust $ between brushOrientation endTgtBwd endTgtFwd endTestTgt = fmap ( MkVector2D ( coords endPt ) • ) . reverseSpline $ joinWithBrush ( withTangent endTgtBwd endBrush ) ( withTangent endTgtFwd endBrush ) endBrush @@ -263,7 +283,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = where brushOrientation :: Orientation - brushOrientation = splineOrientation ( brushShape spt0 ) + brushOrientation = splineOrientation @Double ( brushShape spt0 ) outlineFns :: Seq @@ -349,7 +369,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = brush0 = brushShape sp0 fwdJoin, bwdJoin :: SplinePts Open fwdJoin - | tgtFwd `parallel` prevTgtFwd + | tgtFwd `strictlyParallel` prevTgtFwd = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty } | brushOrientation == convexOrientation [ prevTgtFwd, tgtFwd ] = fmap ( ptOffset • ) @@ -359,7 +379,7 @@ computeStrokeOutline fitParams ptParams brushFn spline@( Spline { splineStart = . reverseSpline $ joinWithBrush ( withTangent tgtFwd brush0 ) ( withTangent prevTgtFwd brush0 ) brush0 bwdJoin - | tgtBwd `parallel` prevTgtBwd + | tgtBwd `strictlyParallel` prevTgtBwd = Spline { splineStart = coords sp0, splineCurves = OpenCurves Empty } | brushOrientation == convexOrientation [ tgtBwd, prevTgtBwd ] = fmap ( ptOffset • ) @@ -433,33 +453,45 @@ outlineFunctions ptParams brushFn sp0 crv = -> ( brush3, Cubic.bezier @( Vector2D Double ) bez, Cubic.bezier' bez ) fwd, bwd :: Double -> ( Point2D Double, Vector2D Double ) fwd t - = ( off t - , if squaredNorm offTgt < epsilon then f' t else offTgt + = ( offset ( withTangent ( fwd' t ) ( brush t ) ) • f t + , fwd' t ) where off :: Double -> Point2D Double off u = offset ( withTangent ( f' u ) ( brush u ) ) • f u - offTgt :: Vector2D Double - offTgt - | t < 0.5 - = 1e9 *^ ( off t --> off (t + 1e-9) ) + offTgt :: Double -> Vector2D Double + offTgt u + | u < 0.5 + = 1e9 *^ ( off u --> off (u + 1e-9) ) | otherwise - = 1e9 *^ ( off (t - 1e-9) --> off t ) + = 1e9 *^ ( off (u - 1e-9) --> off u ) + fwd' :: Double -> Vector2D Double + fwd' u + | squaredNorm ( offTgt u ) < epsilon + = f' u + | otherwise + = offTgt u bwd t - = ( off s - , if squaredNorm offTgt < epsilon then (-1) *^ f' s else offTgt + = ( offset ( withTangent ( (-1) *^ bwd' s ) ( brush s ) ) • f s + , bwd' s ) where s :: Double s = 1 - t off :: Double -> Point2D Double off u = offset ( withTangent ( (-1) *^ f' u ) ( brush u ) ) • f u - offTgt :: Vector2D Double - offTgt - | s < 0.5 - = 1e9 *^ ( off s --> off (s + 1e-9) ) + offTgt :: Double -> Vector2D Double + offTgt u + | u < 0.5 + = 1e9 *^ ( off u --> off (u + 1e-9) ) | otherwise - = 1e9 *^ ( off (s - 1e-9) --> off s ) + = 1e9 *^ ( off (u - 1e-9) --> off u ) + bwd' :: Double -> Vector2D Double + bwd' u + | squaredNorm ( offTgt u ) < epsilon + = (-1) *^ f' u + | otherwise + = offTgt u in ( fwd, bwd ) ----------------------------------- @@ -647,6 +679,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli \tangent vector: " <> show tgt_wanted <> "\n\ \spline: " <> show spline <> "\n" where + ori :: Orientation + ori = splineOrientation @Double spline go :: Int -> ptData -> Curve Open crvData ptData -> StateT ( Vector2D Double ) ( Except Offset ) () go i cp cseg = do tgt_prev <- get @@ -659,8 +693,8 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli tgt_start = startTangent splineStart cp cseg tgt_end = endTangent splineStart cp cseg -- Handle corner. - unless ( parallel tgt_prev tgt_start ) do - for_ ( between tgt_wanted tgt_prev tgt_start ) \ _ -> + unless ( tgt_prev `strictlyParallel` tgt_start ) do + for_ ( between ori tgt_prev tgt_start tgt_wanted ) \ _ -> lift . throwE $ Offset { offsetIndex = i @@ -673,7 +707,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli handleSegment :: Int -> Point2D Double -> Curve Open crvData ( Point2D Double ) -> Vector2D Double -> Except Offset () handleSegment i p0 ( LineTo ( NextPoint p1 ) _ ) tgt0 - | parallel tgt_wanted tgt0 + | tgt_wanted `strictlyParallel` tgt0 , let offset :: Vector2D Double offset = MkVector2D $ lerp @( Vector2D Double ) 0.5 p0 p1 @@ -684,7 +718,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli let tgt1 :: Vector2D Double tgt1 = p1 --> p2 - in for_ ( between tgt_wanted tgt0 tgt1 ) \ t -> + in for_ ( convexCombination tgt0 tgt1 tgt_wanted ) \ t -> throwE $ Offset { offsetIndex = i @@ -715,7 +749,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli mbParam = case mapMaybe correctTangentParam $ solveQuadratic c01 ( 2 * ( c12 - c01 ) ) ( c01 + c23 - 2 * c12 ) of ( t : _ ) -> Just t - _ -> between tgt_wanted tgt0 tgt2 -- fallback in case we couldn't solve the quadratic for some reason + _ -> between ori tgt0 tgt2 tgt_wanted -- fallback in case we couldn't solve the quadratic for some reason in for_ mbParam \ t -> throwE $ Offset @@ -723,41 +757,3 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) = case lastTangent spli , offsetParameter = Just t , offset = MkVector2D $ Cubic.bezier @( Vector2D Double ) bez t } - --- | Finds whether the query vector @ u @ lies between the two provided vectors @ v0 @, @ v1 @. --- --- If so, returns @ t @ in @ [ 0, 1 ] @ such that @ ( 1 - t ) v0 + t v1 @ is a positive multiple of @ u @. -between - :: Vector2D Double -- ^ query vector - -> Vector2D Double -- ^ first vector - -> Vector2D Double -- ^ second vector - -> Maybe Double -between u v0 v1 - | abs c10 < epsilon - = if parallel u v0 - then Just 0 - else if parallel u v1 - then Just 1 - else Nothing - | otherwise - = do - let - t :: Double - t = c0 / c10 - guard ( t > - epsilon && t < 1 + epsilon ) - guard ( epsilon < u ^.^ ( lerp @( Vector2D Double ) t v0 v1 ) ) - Just $ min 1 ( max 0 t ) - - where - c0, c10 :: Double - c0 = v0 `cross` u - c10 = ( v0 ^-^ v1 ) `cross` u - --- | 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. -parallel :: Vector2D Double -> Vector2D Double -> Bool -parallel u v - = abs ( u `cross` v ) < epsilon -- vectors are collinear - && u ^.^ v > epsilon -- vectors point in the same direction (parallel and not anti-parallel) diff --git a/src/lib/Math/Module.hs b/src/lib/Math/Module.hs index f336631..7e8ac89 100644 --- a/src/lib/Math/Module.hs +++ b/src/lib/Math/Module.hs @@ -12,12 +12,15 @@ module Math.Module , squaredNorm, quadrance, distance , proj, projC, closestPointOnSegment , cross + , strictlyParallel, convexCombination ) where -- base import Control.Applicative ( liftA2 ) +import Control.Monad + ( guard ) import Data.Monoid ( Ap(..), Sum(..) ) @@ -34,6 +37,8 @@ import Data.Group ( invert ) -- MetaBrush +import Math.Epsilon + ( epsilon ) import Math.Vector2D ( Vector2D(..), Segment(..) ) @@ -107,6 +112,7 @@ closestPointOnSegment c ( Segment p0 p1 ) t :: r t = projC ( p0 --> c ) v01 +-------------------------------------------------------------------------------- instance Num a => Module a ( Sum a ) where @@ -136,6 +142,47 @@ instance Num a => Inner a ( Vector2D a ) where ( Vector2D x1 y1 ) ^.^ ( Vector2D x2 y2 ) = x1 * x2 + y1 * y2 +-- | Cross-product of two 2D vectors. cross :: Num a => Vector2D a -> Vector2D a -> a cross ( Vector2D x1 y1 ) ( Vector2D x2 y2 ) - = x1 * y2 - x2 * y1 \ No newline at end of file + = x1 * y2 - x2 * y1 + +-- | 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. +strictlyParallel :: RealFloat r => Vector2D r -> Vector2D r -> Bool +strictlyParallel u v + = abs ( u `cross` v ) < epsilon -- vectors are collinear + && u ^.^ v > epsilon -- vectors point in the same direction (parallel and not anti-parallel) + +-- | Finds whether the query vector @ u @ is a convex combination of the two provided vectors @ v0 @, @ v1 @. +-- +-- If so, returns @ t @ in @ [ 0, 1 ] @ such that @ ( 1 - t ) v0 + t v1 @ is a positive multiple of @ u @. +convexCombination + :: forall r + . RealFloat r + => Vector2D r -- ^ first vector + -> Vector2D r -- ^ second vector + -> Vector2D r -- ^ query vector + -> Maybe r +convexCombination v0 v1 u + | abs c10 < epsilon + = if strictlyParallel u v0 + then Just 0 + else if strictlyParallel u v1 + then Just 1 + else Nothing + | otherwise + = do + let + t :: r + t = c0 / c10 + guard ( t > - epsilon && t < 1 + epsilon ) + guard ( epsilon < u ^.^ ( lerp @( Vector2D r ) t v0 v1 ) ) + Just $ min 1 ( max 0 t ) + + where + c0, c10 :: r + c0 = v0 `cross` u + c10 = ( v0 ^-^ v1 ) `cross` u diff --git a/src/lib/Math/Orientation.hs b/src/lib/Math/Orientation.hs index 16fd8c3..bab531e 100644 --- a/src/lib/Math/Orientation.hs +++ b/src/lib/Math/Orientation.hs @@ -1,5 +1,7 @@ +{-# LANGUAGE AllowAmbiguousTypes #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE DerivingStrategies #-} +{-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE GADTs #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} @@ -7,9 +9,16 @@ module Math.Orientation ( Orientation(..), reverseOrientation , convexOrientation, splineOrientation, splineTangents + , between ) where +-- base +import Control.Monad + ( guard ) +import Data.Fixed + ( mod' ) + -- acts import Data.Act ( Torsor((-->)) ) @@ -18,6 +27,12 @@ import Data.Act import Data.Sequence ( Seq(..) ) +-- generic-lens +import Data.Generics.Product.Typed + ( HasType(typed) ) +import Data.GenericLens.Internal + ( view ) + -- MetaBrush import Math.Epsilon ( nearZero ) @@ -29,7 +44,7 @@ import Math.Bezier.Spline , ssplineType ) import Math.Vector2D - ( Point2D, Vector2D ) + ( Point2D, Vector2D(..) ) -------------------------------------------------------------------------------- @@ -58,34 +73,78 @@ convexOrientation _ = CCW -- default -- | Compute the orientation of a spline, assuming tangent vectors have a monotone angle. splineOrientation - :: ( KnownSplineType clo, RealFloat r ) - => Spline clo crvData ( Point2D r ) + :: forall r clo crvData ptData + . ( KnownSplineType clo, RealFloat r, HasType ( Point2D r ) ptData ) + => Spline clo crvData ptData -> Orientation -splineOrientation = convexOrientation . splineTangents +splineOrientation = convexOrientation . splineTangents @r -- | Compute the sequence of tangent vectors given by the control points of a Bézier spline. splineTangents - :: forall clo crvData r - . ( Num r, KnownSplineType clo ) - => Spline clo crvData ( Point2D r ) + :: forall r clo crvData ptData + . ( Num r, KnownSplineType clo, HasType ( Point2D r ) ptData ) + => Spline clo crvData ptData -> [ Vector2D r ] -splineTangents spline@( Spline { splineStart = p0, splineCurves = curves } ) = - case ssplineType @clo of - SOpen - | OpenCurves { openCurves = cs } <- curves - -> go p0 cs - SClosed - | OpenCurves cs@( c :<| _ ) <- splineCurves $ adjustSplineType @Open spline - -> go p0 ( cs :|> c ) - _ -> [] +splineTangents spline@( Spline { splineStart = sp0, splineCurves = curves } ) + | let + p0 :: Point2D r + p0 = view typed sp0 + = case ssplineType @clo of + SOpen + | OpenCurves { openCurves = cs } <- curves + -> go p0 cs + SClosed + | OpenCurves cs@( c :<| _ ) <- splineCurves $ adjustSplineType @Open spline + -> go p0 ( cs :|> c ) + _ -> [] where - go :: Point2D r -> Seq ( Curve Open crvData ( Point2D r ) ) -> [ Vector2D r ] + go :: Point2D r -> Seq ( Curve Open crvData ptData ) -> [ Vector2D r ] go _ Empty = [] go p ( crv :<| crvs ) = case crv of - LineTo { curveEnd = NextPoint q } -> - ( p --> q ) : go q crvs - Bezier2To { controlPoint = cp, curveEnd = NextPoint q } -> - ( p --> cp ) : ( cp --> q ) : go q crvs - Bezier3To { controlPoint1 = cp1, controlPoint2 = cp2, curveEnd = NextPoint q } -> - ( p --> cp1 ) : ( cp2 --> q ) : go q crvs + LineTo { curveEnd = NextPoint sq } + | let + q :: Point2D r + q = view typed sq + -> ( p --> q ) : go q crvs + Bezier2To { controlPoint = scp, curveEnd = NextPoint sq } + | let + cp, q :: Point2D r + cp = view typed scp + q = view typed sq + -> ( p --> cp ) : ( cp --> q ) : go q crvs + Bezier3To { controlPoint1 = scp1, controlPoint2 = scp2, curveEnd = NextPoint sq } + | let + cp1, cp2, q :: Point2D r + cp1 = view typed scp1 + cp2 = view typed scp2 + q = view typed sq + -> ( p --> cp1 ) : ( cp2 --> q ) : go q crvs + +-- | Checks whether a 2D vector lies "in between" two other vectors according to a given orientation, +-- i.e. whether the angle of the query vector lies in between the angles of the start and end vectors. +-- +-- Returns the proportion of the angle the vector is in between, or @Nothing@ if the query vector +-- is not in between. +-- +-- >>> between CCW ( Vector2D 1 0 ) ( Vector2D (-1) 1 ) ( Vector2D 1 1 ) +-- Just 0.3333333333333333 +between + :: forall r + . RealFloat r + => Orientation + -> Vector2D r -- ^ start vector + -> Vector2D r -- ^ end vector + -> Vector2D r -- ^ query vector: is in between the start and end vectors w.r.t. the provided orientation? + -> Maybe r +between CCW ( Vector2D x1 y1 ) ( Vector2D x2 y2 ) ( Vector2D a b ) = + let + τ, η, φ, θ :: r + τ = 2 * pi + η = atan2 y1 x1 + φ = ( atan2 y2 x2 - η ) `mod'` τ + θ = ( atan2 b a - η ) `mod'` τ + in do + guard ( θ < φ ) + pure ( θ / φ ) +between CW v1 v2 u = ( 1 - ) <$> between CCW v2 v1 u