improve intervallic Bézier evaluation

Now we evaluate Bézier curves using an AABB computation. This results
in tighter intervals, which means that the cusp-finding algorithm
is better behaved.
This commit is contained in:
sheaf 2023-04-25 23:07:18 +02:00
parent 7db3cbef33
commit 50d99e1e4b
9 changed files with 193 additions and 59 deletions

View file

@ -13,9 +13,9 @@ import Control.Exception
--------------------------------------------------------------------------------
{-# INLINE assert #-}
assert :: String -> a -> a
assert :: Bool -> String -> a -> a
#ifdef ASSERTS
assert message _ = throw ( AssertionFailed message )
assert False message _ = throw ( AssertionFailed message )
#else
assert _ a = a
assert _ _ a = a
#endif

View file

@ -220,7 +220,7 @@ addToAnchor anchor newSpline = over ( field' @"documentContent" . field' @"strok
setBrushData = set ( field @"brushParams" ) ( brushParams ( splineEnd prevSpline ) )
in prevSpline <> fmap setBrushData newSpline
| otherwise
= assert ( "addToAnchor: trying to add to closed spline " <> show strokeUnique )
= assert False ( "addToAnchor: trying to add to closed spline " <> show strokeUnique )
prevSpline -- should never add to a closed spline
= overStrokeSpline updateSpline stroke
| otherwise

View file

@ -6,10 +6,12 @@ module Math.Bezier.Cubic
( Bezier(..)
, fromQuadratic
, bezier, bezier', bezier'', bezier'''
, derivative
, curvature, squaredCurvature, signedCurvature
, subdivide
, ddist, closestPoint
, drag, selfIntersectionParameters
, extrema
)
where
@ -105,6 +107,11 @@ bezier ( Bezier {..} ) t =
( Quadratic.bezier @v ( Quadratic.Bezier p0 p1 p2 ) t )
( Quadratic.bezier @v ( Quadratic.Bezier p1 p2 p3 ) t )
-- | The derivative of a Cubic Bézier curve, as a quadratic Bézier curve.
derivative :: ( Group v, Module r v ) => Bezier v -> Quadratic.Bezier v
derivative ( Bezier {..} ) = ( Ring.fromInteger 3 *^ )
<$> Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 )
-- | Derivative of a cubic Bézier curve.
bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
bezier' ( Bezier {..} )
@ -255,3 +262,11 @@ selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2
c0 = f2
c1 = f3 - 2 * f2
c2 = f1 + f2 - f3
-- | Extremal values of the Bézier parameter for a cubic Bézier curve.
extrema :: RealFloat r => Bezier r -> [ r ]
extrema ( Bezier {..} ) = solveQuadratic c b a
where
a = p3 - 3 * p2 + 3 * p1 - p0
b = 2 * ( p0 - 2 * p1 + p2 )
c = p1 - p0

View file

@ -9,6 +9,7 @@ module Math.Bezier.Quadratic
, subdivide
, ddist, closestPoint
, interpolate
, extrema
)
where
@ -188,3 +189,9 @@ interpolate p0 p2 t q = Bezier {..}
p1 = ( ( 0.5 * ( t - 1 ) / t ) *^ ( q --> p0 :: v )
^+^ ( 0.5 * t / ( t - 1 ) ) *^ ( q --> p2 :: v )
) q
-- | Extremal values of the Bézier parameter for a quadratic Bézier curve.
extrema :: Fractional r => Bezier r -> [ r ]
extrema ( Bezier {..} ) = [ t ]
where
t = ( p0 - p1 ) / ( p0 - 2 * p1 + p2 )

View file

@ -235,6 +235,7 @@ computeStrokeOutline ::
, HasChainRule Double 2 brushParams
, HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams )
, Traversable ( D 2 brushParams )
, Representable Double usedParams
-- Debugging.
, Show ptData, Show brushParams
@ -513,6 +514,9 @@ outlineFunction
, HasChainRule ( 𝕀 Double ) 3 ( 𝕀 brushParams )
, Traversable ( D 2 brushParams )
-- Computing AABBs
, Representable Double usedParams
-- Debugging.
, Show ptData, Show brushParams
)
@ -600,13 +604,16 @@ outlineFunction ptParams toBrushParams brushFromParams = \ sp0 crv ->
pathAndUsedParams :: forall k i arr crvData ptData usedParams
. ( HasType ( 2 ) ptData
, CurveOrder k
, HasBézier k i
, arr ~ C k
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
, Module ( I i Double ) ( T ( I i usedParams ) )
, Torsor ( T ( I i usedParams ) ) ( I i usedParams )
, Module Double ( T usedParams )
, Representable Double usedParams
, Torsor ( T usedParams ) usedParams
)
=> ( I i ( 1 ) -> I i Double )
-> ( forall a. a -> I i a )
@ -618,16 +625,16 @@ pathAndUsedParams co toI ptParams sp0 crv =
case crv of
LineTo { curveEnd = NextPoint sp1 }
| let seg = Segment sp0 sp1
-> ( line @k @i co ( fmap ( toI . coords ) seg )
, line @k @i co ( fmap ( toI . ptParams ) seg ) )
-> ( line @k @i @( 2 ) co ( fmap ( toI . coords ) seg )
, line @k @i @usedParams co ( fmap ( toI . ptParams ) seg ) )
Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 }
| let bez2 = Quadratic.Bezier sp0 sp1 sp2
-> ( bezier2 @k @i co ( fmap ( toI . coords ) bez2 )
, bezier2 @k @i co ( fmap ( toI . ptParams ) bez2 ) )
-> ( bezier2 @k @i @( 2 ) co ( fmap ( toI . coords ) bez2 )
, bezier2 @k @i @usedParams co ( fmap ( toI . ptParams ) bez2 ) )
Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 }
| let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3
-> ( bezier3 @k @i co ( fmap ( toI . coords ) bez3 )
, bezier3 @k @i co ( fmap ( toI . ptParams ) bez3 ) )
-> ( bezier3 @k @i @( 2 ) co ( fmap ( toI . coords ) bez3 )
, bezier3 @k @i @usedParams co ( fmap ( toI . ptParams ) bez3 ) )
-----------------------------------
-- Various utility functions
@ -923,7 +930,7 @@ withTangent tgt_wanted spline@( Spline { splineStart } )
}
splineCurveFns :: forall k i
. ( CurveOrder k
. ( HasBézier k i
, D k ( I i ( 1 ) ) ~ D k ( 1 )
, Module ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) ) )
@ -940,14 +947,14 @@ splineCurveFns co spls
where
curveFn :: I i ( 2 )
-> Curve Open () ( I i ( 2 ) )
-> ( C k ( I i ( 1 ) ) ( I i ( 2 ) ) )
-> C k ( I i ( 1 ) ) ( I i ( 2 ) )
curveFn p0 = \case
LineTo { curveEnd = NextPoint p1 }
-> line @k @i co $ Segment p0 p1
-> line @k @i @( 2 ) co $ Segment p0 p1
Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 }
-> bezier2 @k @i co $ Quadratic.Bezier p0 p1 p2
-> bezier2 @k @i @( 2 ) co $ Quadratic.Bezier p0 p1 p2
Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 }
-> bezier3 @k @i co $ Cubic.Bezier p0 p1 p2 p3
-> bezier3 @k @i @( 2 ) co $ Cubic.Bezier p0 p1 p2 p3
-- | Solve the envelope equations at a given point \( t = t_0 \), to find
-- \( s_0 \) such that \( c(t_0, s_0) \) is on the envelope of the brush stroke.
@ -1072,7 +1079,8 @@ instance Applicative ZipSeq where
{-# INLINE liftA2 #-}
brushStrokeData :: forall k brushParams i arr
. ( CurveOrder k, arr ~ C k
. ( arr ~ C k
, HasBézier k i, HasEnvelopeEquation k
, Differentiable k i brushParams
, HasChainRule ( I i Double ) k ( I i ( 1 ) )
, Applicative ( D k ( 1 ) )
@ -1084,7 +1092,6 @@ brushStrokeData :: forall k brushParams i arr
, D k ( I i ( 2 ) ) ~ D k ( 2 )
, Cross ( I i Double ) ( T ( I i ( 2 ) ) )
, Torsor ( T ( I i ( 2 ) ) ) ( I i ( 2 ) )
-- , Coercible ( I i ( 1 ) ) ( I i Double )
, Show brushParams
)
=> ( I i Double -> I i ( 1 ) )

View file

@ -5,17 +5,16 @@
module Math.Bezier.Stroke.EnvelopeEquation
( StrokeDatum(..)
, CurveOrder(..)
, HasBézier(..)
, HasEnvelopeEquation(..)
) where
-- base
import Prelude hiding ( Num(..), (^) )
import Data.Coerce
( Coercible, coerce )
import Data.Functor.Identity
( Identity(..) )
import Data.Kind
( Type, Constraint )
import Data.List.NonEmpty
( NonEmpty(..) )
import GHC.TypeNats
( Nat, type (-) )
@ -72,32 +71,40 @@ deriving stock instance Show ( StrokeDatum 3 𝕀 )
--------------------------------------------------------------------------------
type CurveOrder :: Nat -> Constraint
class CurveOrder k where
type HasBézier :: forall {t}. Nat -> t -> Constraint
class HasBézier k i where
-- | Linear interpolation, as a differentiable function.
line :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
line :: forall b
. ( Module Double ( T b ), Torsor ( T b ) b
, Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b )
, D k ( I i ( 1 ) ) ~ D k ( 1 )
)
=> ( I i ( 1 ) -> I i Double )
-> Segment b -> C k ( I i ( 1 ) ) b
-> Segment ( I i b ) -> C k ( I i ( 1 ) ) ( I i b )
-- | A quadratic Bézier curve, as a differentiable function.
bezier2 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
bezier2 :: forall b
. ( Module Double ( T b ), Torsor ( T b ) b
, Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b )
, Representable Double b
, D k ( I i ( 1 ) ) ~ D k ( 1 )
)
=> ( I i ( 1 ) -> I i Double )
-> Quadratic.Bezier b -> C k ( I i ( 1 ) ) b
-> Quadratic.Bezier ( I i b ) -> C k ( I i ( 1 ) ) ( I i b )
-- | A cubic Bézier curve, as a differentiable function.
bezier3 :: forall i b
. ( Module ( I i Double ) ( T b ), Torsor ( T b ) b
bezier3 :: forall b
. ( Module Double ( T b ), Torsor ( T b ) b
, Module ( I i Double ) ( T (I i b ) ), Torsor ( T ( I i b ) ) ( I i b )
, Representable Double b
, D k ( I i ( 1 ) ) ~ D k ( 1 )
)
=> ( I i ( 1 ) -> I i Double )
-> Cubic.Bezier b -> C k ( I i ( 1 ) ) b
-> Cubic.Bezier ( I i b ) -> C k ( I i ( 1 ) ) ( I i b )
type HasEnvelopeEquation :: Nat -> Constraint
class HasEnvelopeEquation k where
uncurryD :: D k a ~ D k ( 1 )
=> D k ( 1 ) ( C k a b ) -> a -> D k ( 2 ) b
@ -137,7 +144,7 @@ class CurveOrder k where
-> ( D ( k - 1 ) ( I i ( 2 ) ) ( I i ( 1 ) )
, D ( k - 2 ) ( I i ( 2 ) ) ( T ( I i ( 2 ) ) ) )
instance CurveOrder 2 where
instance HasBézier 2 () where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D21 ( lerp @( T b ) t a b )
@ -156,6 +163,8 @@ instance CurveOrder 2 where
( Cubic.bezier' bez t )
( Cubic.bezier'' bez t )
instance HasEnvelopeEquation 2 where
uncurryD = uncurryD2
brushStroke ( D21 p p_t p_tt ) ( D22 b b_t b_s b_tt b_ts b_ss ) =
@ -185,7 +194,7 @@ instance CurveOrder 2 where
--
-- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s.
instance CurveOrder 3 where
instance HasBézier 3 () where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
@ -208,6 +217,8 @@ instance CurveOrder 3 where
( Cubic.bezier'' bez t )
( Cubic.bezier''' bez )
instance HasEnvelopeEquation 3 where
uncurryD = uncurryD3
brushStroke
@ -246,3 +257,79 @@ instance CurveOrder 3 where
( T $ co ee_t ) ( T $ co ee_s )
( T $ co ee_tt) ( T $ co ee_ts ) ( T $ co ee_ss )
, D12 𝛿E𝛿sdcdt ( T 𝛿E𝛿sdcdt_t ) ( T 𝛿E𝛿sdcdt_s ) )
instance HasBézier 3 𝕀 where
line co ( Segment a b :: Segment b ) =
D \ ( co -> t ) ->
D31 ( lerp @( T b ) t a b )
( a --> b )
origin
origin
bezier2 co ( bez :: Quadratic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( aabb bez ( `evaluateQuadratic` t ) )
( Quadratic.bezier' bez t )
( Quadratic.bezier'' bez )
origin
bezier3 co ( bez :: Cubic.Bezier b ) =
D \ ( co -> t ) ->
D31 ( aabb bez ( `evaluateCubic` t ) )
( T $ aabb ( fmap unT $ Cubic.derivative ( fmap T bez ) ) ( `evaluateQuadratic` t ) )
( Cubic.bezier'' bez t )
( Cubic.bezier''' bez )
{- Note [Computing Béziers over intervals]
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Here's how we evaluate a Bézier curve when the coefficients are intervals.
As we are using axis-aligned interval arithmetic, we reduce to the 1D situation.
Now, the formulas for Bézier curves (with value of the t parameter in [0,1])
are convex combinations
b_0(t) [p_0,q_0] + b_1(t) [p_1,q_1] + ... + b_n(t) [p_n,q_n]
-- Here b_1, ..., b_n are Bernstein polynomials.
This means that the minimum value attained by the Bézier curve as we vary
both the time parameter and the values of the points within their respective
intervals will occur when we take
b_0(t) p_0 + b_1(t) p_1 + ... + b_n(t) p_n
and the maximum when we take
b_0(t) q_0 + b_1(t) q_1 + ... + b_n(t) q_n
For each of these, we can compute a minimum/maximum by computing an axis-aligned
bounding box for the Bézier curve. This is done by computing the derivative
with respect to t and root-finding.
-}
-- | Evaluate a cubic Bézier curve, when both the coefficients and the
-- parameter are intervals.
evaluateCubic :: Cubic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
evaluateCubic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = fmap inf bez
sup_bez = fmap sup bez
mins = fmap (Cubic.bezier @( T Double ) inf_bez)
$ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema inf_bez ) )
maxs = fmap (Cubic.bezier @( T Double ) sup_bez)
$ inf t :| ( sup t : filter ( inside t ) ( Cubic.extrema sup_bez ) )
in 𝕀 ( minimum mins ) ( maximum maxs )
-- | Evaluate a quadratic Bézier curve, when both the coefficients and the
-- parameter are intervals.
evaluateQuadratic :: Quadratic.Bezier ( 𝕀 Double ) -> 𝕀 Double -> 𝕀 Double
evaluateQuadratic bez t =
-- assert (inf t >= 0 && sup t <= 1) "evaluateCubic: t ⊊ [0,1]" $ -- Requires t ⊆ [0,1]
let inf_bez = fmap inf bez
sup_bez = fmap sup bez
mins = fmap (Quadratic.bezier @( T Double ) inf_bez)
$ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema inf_bez ) )
maxs = fmap (Quadratic.bezier @( T Double ) sup_bez)
$ inf t :| ( sup t : filter ( inside t ) ( Quadratic.extrema sup_bez ) )
in 𝕀 ( minimum mins ) ( maximum maxs )

View file

@ -12,11 +12,15 @@ module Math.Interval
, scaleInterval
, 𝕀
, singleton, nonDecreasing
, inside
, aabb
)
where
-- base
import Prelude hiding ( Num(..) )
import Data.List.NonEmpty
( NonEmpty(..) )
-- acts
import Data.Act
@ -32,20 +36,18 @@ import Data.Group.Generics
-- splines
import Math.Algebra.Dual
import Math.Algebra.Dual.Internal
( chainRule1NQ )
import Math.Interval.Internal
( 𝕀(𝕀), inf, sup, scaleInterval )
import Math.Linear
( (..), T(..)
, RepresentableQ(..)
, RepresentableQ(..), Representable(..)
)
import Math.Module
import Math.Monomial
import Math.Ring
--------------------------------------------------------------------------------
-- Interval arithmetic using rounded-hw library.
-- Interval arithmetic.
type 𝕀 n = 𝕀 ( n )
type instance D k ( 𝕀 v ) = D k v
@ -57,6 +59,8 @@ singleton a = 𝕀 a a
nonDecreasing :: ( a -> b ) -> 𝕀 a -> 𝕀 b
nonDecreasing f ( 𝕀 lo hi ) = 𝕀 ( f lo ) ( f hi )
inside :: Ord a => 𝕀 a -> a -> Bool
inside ( 𝕀 lo hi ) x = x >= lo && x <= hi
deriving via ViaAbelianGroup ( T ( 𝕀 Double ) )
instance Semigroup ( T ( 𝕀 Double ) )
@ -71,12 +75,19 @@ instance Torsor ( T ( 𝕀 Double ) ) ( 𝕀 Double ) where
a --> b = T $ b - a
-------------------------------------------------------------------------------
-- Lattices.
aabb :: ( Representable r v, Ord r, Functor f )
=> f ( 𝕀 v ) -> ( f ( 𝕀 r ) -> 𝕀 r ) -> 𝕀 v
aabb fv extrema = tabulate \ i -> extrema ( fmap ( `index` i ) fv )
--------------------------------------------------------------------------------
instance Module ( 𝕀 Double ) ( T ( 𝕀 0 ) ) where
origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )
T a ^+^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^+^ T $$( indexQ [|| b ||] i ) ||] )
T a ^-^ T b = T $$( tabulateQ \ i -> [|| unT $ T $$( indexQ [|| a ||] i ) ^-^ T $$( indexQ [|| b ||] i ) ||] )
k *^ T a = T $$( tabulateQ \ i -> [|| unT $ k *^ T $$( indexQ [|| a ||] i ) ||] )
origin = T $ singleton 0
_ ^+^ _ = T $ singleton 0
_ ^-^ _ = T $ singleton 0
_ *^ _ = T $ singleton 0
instance Module ( 𝕀 Double ) ( T ( 𝕀 1 ) ) where
origin = T $$( tabulateQ \ _ -> [|| unT $ origin ||] )

View file

@ -1,7 +1,8 @@
{-# LANGUAGE ScopedTypeVariables #-}
module Math.Ring
( AbelianGroup(..), Signed(..), Ring(..), Field(..), Transcendental(..)
( AbelianGroup(..), Signed(..), Ring(..), Field(..)
, Floating(..), Transcendental(..)
, ViaPrelude(..), ViaAbelianGroup(..)
@ -11,7 +12,7 @@ module Math.Ring
-- base
import Prelude ( Num, Fractional )
import Prelude hiding ( Num(..), Fractional(..) )
import Prelude hiding ( Num(..), Fractional(..), Floating(..) )
import qualified Prelude
import Control.Applicative
( liftA2 )
@ -60,11 +61,13 @@ class Ring a => Field a where
recip x = fromInteger 1 / x
x / y = x * recip y
class Field a => Floating a where
sqrt :: a -> a
class Field a => Transcendental a where
pi :: a
cos :: a -> a
sin :: a -> a
-- sqrt :: a -> a
--------------------------------------------------------------------------------
@ -126,11 +129,13 @@ instance Fractional a => Field ( ViaPrelude a ) where
fromRational = coerce $ Prelude.fromRational @a
(/) = coerce $ (Prelude./) @a
instance Floating a => Transcendental ( ViaPrelude a ) where
instance Prelude.Floating a => Floating ( ViaPrelude a ) where
sqrt = coerce $ Prelude.sqrt @a
instance Prelude.Floating a => Transcendental ( ViaPrelude a ) where
pi = coerce $ Prelude.pi @a
sin = coerce $ Prelude.sin @a
cos = coerce $ Prelude.cos @a
-- sqrt = coerce $ Prelude.sqrt @a
--------------------------------------------------------------------------------
@ -160,12 +165,14 @@ deriving via ViaPrelude Float instance AbelianGroup Float
deriving via ViaPrelude Float instance Ring Float
deriving via ViaPrelude Float instance Signed Float
deriving via ViaPrelude Float instance Field Float
deriving via ViaPrelude Float instance Floating Float
deriving via ViaPrelude Float instance Transcendental Float
deriving via ViaPrelude Double instance AbelianGroup Double
deriving via ViaPrelude Double instance Ring Double
deriving via ViaPrelude Double instance Signed Double
deriving via ViaPrelude Double instance Field Double
deriving via ViaPrelude Double instance Floating Double
deriving via ViaPrelude Double instance Transcendental Double
--------------------------------------------------------------------------------

View file

@ -62,26 +62,26 @@ solveQuadratic
-> a -- ^ linear coefficient
-> a -- ^ quadratic coefficient
-> [ a ]
solveQuadratic a0 a1 a2
| nearZero a1 && nearZero a2
= if nearZero a0
solveQuadratic c b a
| nearZero b && nearZero a
= if nearZero c
then [ 0, 0.5, 1 ] -- convention
else []
| nearZero ( a0 * a0 * a2 / ( a1 * a1 ) )
= [ -a0 / a1 ]
| nearZero ( c * c * a / ( b * b ) )
= [ -c / b ]
| disc < 0
= [] -- non-real solutions
| otherwise
= let
r :: a
r =
if a1 >= 0
then 2 * a0 / ( -a1 - sqrt disc )
else 0.5 * ( -a1 + sqrt disc ) / a2
in [ r, -r - a1 / a2 ]
if b >= 0
then 2 * c / ( -b - sqrt disc )
else 0.5 * ( -b + sqrt disc ) / a
in [ r, c / ( a * r ) ]
where
disc :: a
disc = a1 * a1 - 4 * a0 * a2
disc = b * b - 4 * a * c
--------------------------------------------------------------------------------
-- Root finding using Laguerre's method