From 236055b4ca58ced124447d81ed4475b408b645b6 Mon Sep 17 00:00:00 2001 From: sheaf Date: Sat, 21 Jan 2023 15:24:08 +0100 Subject: [PATCH] do the interval brush stroking at degree 3 --- MetaBrush.cabal | 1 + src/metabrushes/MetaBrush/Asset/Brushes.hs | 49 ++- src/metabrushes/MetaBrush/Brush.hs | 13 +- src/metabrushes/MetaBrush/Records.hs | 3 + src/splines/Math/Algebra/Dual.hs | 14 +- src/splines/Math/Bezier/Cubic.hs | 14 +- src/splines/Math/Bezier/Quadratic.hs | 4 +- src/splines/Math/Bezier/Stroke.hs | 360 ++++++++------------- src/splines/Math/Differentiable.hs | 8 +- src/splines/Math/Interval.hs | 2 +- src/splines/Math/Module.hs | 12 +- src/splines/Math/Orientation.hs | 4 +- 12 files changed, 205 insertions(+), 279 deletions(-) diff --git a/MetaBrush.cabal b/MetaBrush.cabal index 62b376f..c776dd5 100644 --- a/MetaBrush.cabal +++ b/MetaBrush.cabal @@ -179,6 +179,7 @@ library splines , Math.Bezier.Quadratic , Math.Bezier.Spline , Math.Bezier.Stroke + , Math.Bezier.Stroke.EnvelopeEquation , Math.Differentiable , Math.Epsilon , Math.Interval diff --git a/src/metabrushes/MetaBrush/Asset/Brushes.hs b/src/metabrushes/MetaBrush/Asset/Brushes.hs index b2e6637..27c7af1 100644 --- a/src/metabrushes/MetaBrush/Asset/Brushes.hs +++ b/src/metabrushes/MetaBrush/Asset/Brushes.hs @@ -29,8 +29,6 @@ import qualified Data.HashMap.Strict as HashMap -- MetaBrush import Math.Algebra.Dual import Math.Bezier.Spline -import Math.Differentiable - ( DiffInterp ) import Math.Interval ( type I ) import Math.Linear @@ -94,46 +92,65 @@ circleSpline p = sequenceA $ lastCrv = Bezier3To ( p κ -1 ) ( p 1 -κ ) BackToStart () -circleBrush :: forall i k - . ( k ~ 2, DiffInterp i ( Record CircleBrushFields ) ) +circleBrush :: forall i k irec + . ( irec ~ I i ( Record CircleBrushFields ) + , Module + ( D k irec ( I i Double ) ) + ( D k irec ( I i ( ℝ 2 ) ) ) + , Module ( I i Double ) ( T ( I i Double ) ) + , HasChainRule ( I i Double ) k irec + , Representable ( I i Double ) irec + , Applicative ( D k irec ) + ) => Proxy# i -> ( forall a. a -> I i a ) - -> C k ( I i ( Record CircleBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) + -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) circleBrush _ mkI = D \ params -> - let r :: D k ( I i ( Record CircleBrushFields ) ) ( I i Double ) + let r :: D k irec( I i Double ) r = runD ( var @_ @k ( Fin 1 ) ) params - mkPt :: Double -> Double -> D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) + mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) ) mkPt ( kon -> x ) ( kon -> y ) = ( x * r ) *^ e_x ^+^ ( y * r ) *^ e_y in circleSpline @i @k @( Record CircleBrushFields ) @( ℝ 2 ) mkPt where - e_x, e_y :: D k ( I i ( Record CircleBrushFields ) ) ( I i ( ℝ 2 ) ) + e_x, e_y :: D k irec ( I i ( ℝ 2 ) ) e_x = pure $ mkI $ ℝ2 1 0 e_y = pure $ mkI $ ℝ2 0 1 - kon = konst @( I i Double ) @k @( I i ( Record CircleBrushFields ) ) . mkI + kon = konst @( I i Double ) @k @irec . mkI -ellipseBrush :: forall i k - . ( k ~ 2, DiffInterp i ( Record EllipseBrushFields ) ) +ellipseBrush :: forall i k irec + . ( irec ~ I i ( Record EllipseBrushFields ) + , Module + ( D k irec ( I i Double ) ) + ( D k irec ( I i ( ℝ 2 ) ) ) + , Module ( I i Double ) ( T ( I i Double ) ) + , HasChainRule ( I i Double ) k irec + , Representable ( I i Double ) irec + , Applicative ( D k irec ) + , Transcendental ( D k irec ( I i Double ) ) + -- TODO: make a synonym for the above... + -- it seems DiffInterp isn't exactly right + ) => Proxy# i -> ( forall a. a -> I i a ) - -> C k ( I i ( Record EllipseBrushFields ) ) ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) + -> C k irec ( Spline 'Closed () ( I i ( ℝ 2 ) ) ) ellipseBrush _ mkI = D \ params -> - let a, b, phi :: D k ( I i ( Record EllipseBrushFields ) ) ( I i Double ) + let a, b, phi :: D k irec ( I i Double ) a = runD ( var @_ @k ( Fin 1 ) ) params b = runD ( var @_ @k ( Fin 2 ) ) params phi = runD ( var @_ @k ( Fin 3 ) ) params - mkPt :: Double -> Double -> D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) + mkPt :: Double -> Double -> D k irec ( I i ( ℝ 2 ) ) mkPt ( kon -> x ) ( kon -> y ) = ( x * a * cos phi - y * b * sin phi ) *^ e_x ^+^ ( y * b * cos phi + x * a * sin phi ) *^ e_y in circleSpline @i @k @( Record EllipseBrushFields ) @( ℝ 2 ) mkPt where - e_x, e_y :: D k ( I i ( Record EllipseBrushFields ) ) ( I i ( ℝ 2 ) ) + e_x, e_y :: D k irec ( I i ( ℝ 2 ) ) e_x = pure $ mkI $ ℝ2 1 0 e_y = pure $ mkI $ ℝ2 0 1 - kon = konst @( I i Double ) @k @( I i ( Record EllipseBrushFields ) ) . mkI + kon = konst @( I i Double ) @k @irec . mkI diff --git a/src/metabrushes/MetaBrush/Brush.hs b/src/metabrushes/MetaBrush/Brush.hs index e4b3996..8ea2d0c 100644 --- a/src/metabrushes/MetaBrush/Brush.hs +++ b/src/metabrushes/MetaBrush/Brush.hs @@ -40,15 +40,14 @@ import qualified Data.Text as Text -- MetaBrush import Math.Algebra.Dual - ( type (~>) ) -import Math.Linear -import Math.Interval - ( type I, Extent(Point, Interval) ) - + ( C ) import Math.Bezier.Spline ( SplineType(Closed), Spline ) import Math.Differentiable - ( DiffInterp ) + ( DiffInterp, ExtentOrder ) +import Math.Interval + ( type I, Extent(Point, Interval) ) +import Math.Linear import MetaBrush.Records ( KnownSymbols, Length, Record ) import MetaBrush.Serialisable @@ -67,7 +66,7 @@ data WithParams params f = . ( DiffInterp i params ) => Proxy# i -> ( forall a. a -> I i a ) - -> I i params ~> f ( I i ( ℝ 2 ) ) + -> C ( ExtentOrder i ) ( I i params ) ( f ( I i ( ℝ 2 ) ) ) } -------------------------------------------------------------------------------- diff --git a/src/metabrushes/MetaBrush/Records.hs b/src/metabrushes/MetaBrush/Records.hs index 58ea2b1..49711fc 100644 --- a/src/metabrushes/MetaBrush/Records.hs +++ b/src/metabrushes/MetaBrush/Records.hs @@ -150,6 +150,9 @@ deriving newtype instance HasChainRule Double 2 ( ℝ ( Length ks ) ) deriving via 𝕀ℝ ( Length ks ) instance HasChainRule ( 𝕀 Double ) 2 ( 𝕀ℝ ( Length ks ) ) => HasChainRule ( 𝕀 Double ) 2 ( 𝕀 ( Record ks ) ) +deriving via 𝕀ℝ ( Length ks ) + instance HasChainRule ( 𝕀 Double ) 3 ( 𝕀ℝ ( Length ks ) ) + => HasChainRule ( 𝕀 Double ) 3 ( 𝕀 ( Record ks ) ) -------------------------------------------------------------------------------- diff --git a/src/splines/Math/Algebra/Dual.hs b/src/splines/Math/Algebra/Dual.hs index cc5e841..f6f341d 100644 --- a/src/splines/Math/Algebra/Dual.hs +++ b/src/splines/Math/Algebra/Dual.hs @@ -10,7 +10,7 @@ -dsuppress-unfoldings -dsuppress-coercions #-} module Math.Algebra.Dual - ( C(..), D, type (~>), type (~~>) + ( C(..), D , HasChainRule(..), chainRule , uncurryD2, uncurryD3 , linear, fun, var @@ -49,11 +49,6 @@ type C :: Nat -> Type -> Type -> Type newtype C k u v = D { runD :: u -> D k u v } deriving stock instance Functor ( D k u ) => Functor ( C k u ) --- | \( C^2 \)-differentiable mappings. -type (~>) = C 2 --- | \( C^3 \)-differentiable mappings. -type (~~>) = C 3 - -- | @D k u v@ is the space of @k@-th order germs of functions from @u@ to @v@, -- represented by the algebra: -- @@ -64,6 +59,10 @@ type D :: Nat -> Type -> Type -> Type type family D k u type instance D k ( ℝ 0 ) = D𝔸0 +type instance D 0 ( ℝ 1 ) = D𝔸0 +type instance D 0 ( ℝ 2 ) = D𝔸0 +type instance D 0 ( ℝ 3 ) = D𝔸0 +type instance D 0 ( ℝ 4 ) = D𝔸0 type instance D 1 ( ℝ 1 ) = D1𝔸1 type instance D 1 ( ℝ 2 ) = D1𝔸2 @@ -93,6 +92,9 @@ instance ( Applicative ( D k u ), Module r ( T v ) ) -------------------------------------------------------------------------------- +-- TODO: split up this class into the chain rule operation +-- and all the other operations. + -- | @HasChainRule r k v@ means we have a chain rule -- with @D k v w@ in the middle, for any @r@-module @w@. class HasChainRule r k v where diff --git a/src/splines/Math/Bezier/Cubic.hs b/src/splines/Math/Bezier/Cubic.hs index b40a76f..6bf3c50 100644 --- a/src/splines/Math/Bezier/Cubic.hs +++ b/src/splines/Math/Bezier/Cubic.hs @@ -5,7 +5,7 @@ module Math.Bezier.Cubic ( Bezier(..) , fromQuadratic - , bezier, bezier', bezier'' + , bezier, bezier', bezier'', bezier''' , curvature, squaredCurvature, signedCurvature , subdivide , ddist, closestPoint @@ -56,8 +56,8 @@ import Math.Epsilon import Math.Module ( Module (..) , lerp - , Inner(..), norm, squaredNorm - , cross + , Inner((^.^)), norm, squaredNorm + , Cross((×)) ) import Math.Roots ( realRoots, solveQuadratic ) @@ -119,6 +119,12 @@ bezier'' ( Bezier {..} ) t ( p1 --> p0 ^+^ p1 --> p2 ) ( p2 --> p1 ^+^ p2 --> p3 ) +-- | Third derivative of a cubic Bézier curve. +bezier''' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> v +bezier''' ( Bezier {..} ) + = ( Ring.fromInteger 6 *^ ) + $ ( ( p0 --> p3 ) ^+^ Ring.fromInteger 3 *^ ( p2 --> p1 ) ) + -- | Curvature of a cubic Bézier curve. curvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r curvature bez t = sqrt $ squaredCurvature @v bez t @@ -140,7 +146,7 @@ squaredCurvature bez t -- | Signed curvature of a planar cubic Bézier curve. signedCurvature :: Bezier ( ℝ 2 ) -> Double -> Double -signedCurvature bez t = ( g' `cross` g'' ) / norm g' ^ ( 3 :: Int ) +signedCurvature bez t = ( g' × g'' ) / norm g' ^ ( 3 :: Int ) where g', g'' :: T ( ℝ 2 ) g' = bezier' @( T ( ℝ 2 ) ) bez t diff --git a/src/splines/Math/Bezier/Quadratic.hs b/src/splines/Math/Bezier/Quadratic.hs index c8c5c72..3da728d 100644 --- a/src/splines/Math/Bezier/Quadratic.hs +++ b/src/splines/Math/Bezier/Quadratic.hs @@ -54,7 +54,7 @@ import Math.Module ( Module (..) , lerp , Inner(..), norm, squaredNorm - , cross + , Cross((×)) ) import Math.Roots ( realRoots ) @@ -119,7 +119,7 @@ squaredCurvature bez t -- | Signed curvature of a planar quadratic Bézier curve. signedCurvature :: Bezier ( ℝ 2 ) -> Double -> Double -signedCurvature bez t = ( g' `cross` g'' ) / norm g' ^ ( 3 :: Int ) +signedCurvature bez t = ( g' × g'' ) / norm g' ^ ( 3 :: Int ) where g', g'' :: T ( ℝ 2 ) g' = bezier' @( T ( ℝ 2 ) ) bez t diff --git a/src/splines/Math/Bezier/Stroke.hs b/src/splines/Math/Bezier/Stroke.hs index 85cae93..c8c82bc 100644 --- a/src/splines/Math/Bezier/Stroke.hs +++ b/src/splines/Math/Bezier/Stroke.hs @@ -1,5 +1,4 @@ {-# LANGUAGE AllowAmbiguousTypes #-} -{-# LANGUAGE PartialTypeSignatures #-} {-# LANGUAGE QuantifiedConstraints #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE UndecidableInstances #-} @@ -32,7 +31,7 @@ import Control.Monad.ST import Data.Bifunctor ( Bifunctor(bimap) ) import Data.Coerce - ( Coercible, coerce ) + ( Coercible ) import Data.Foldable ( for_, toList ) import Data.Functor.Identity @@ -49,6 +48,8 @@ import GHC.STRef ( STRef(..), readSTRef, writeSTRef ) import GHC.Generics ( Generic, Generic1 ) +import GHC.TypeNats + ( type (-) ) -- acts import Data.Act @@ -113,14 +114,17 @@ import Math.Bezier.Spline , showSplinePoints ) import qualified Math.Bezier.Quadratic as Quadratic +import Math.Bezier.Stroke.EnvelopeEquation import Math.Differentiable - ( Differentiable, DiffInterp ) + ( Differentiable, DiffInterp + , type ExtentOrder + ) import Math.Epsilon ( epsilon ) import Math.Interval import Math.Linear import Math.Module - ( Module(..), Inner((^.^)), Cross(cross), Interpolatable + ( Module(..), Inner((^.^)), Cross((×)), Interpolatable , lerp, convexCombination, strictlyParallel ) import Math.Orientation @@ -197,15 +201,15 @@ computeStrokeOutline :: , NFData ptData, NFData crvData -- Differentiability. - , DiffInterp 'Point brushParams - , DiffInterp 'Interval brushParams , Interpolatable Double usedParams , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) - , HasChainRule Double 2 brushParams - , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams ) - , HasChainRule Double 2 usedParams - , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams ) - , Traversable ( D 2 brushParams ) + , DiffInterp 'Point brushParams + , DiffInterp 'Interval brushParams + , HasChainRule Double ( ExtentOrder 'Point ) usedParams + , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams ) + , HasChainRule Double ( ExtentOrder 'Point ) brushParams + , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams ) + , Traversable ( D ( ExtentOrder 'Point ) brushParams ) -- Debugging. , Show ptData, Show brushParams @@ -218,7 +222,9 @@ computeStrokeOutline :: . DiffInterp i brushParams => Proxy# i -> ( forall a. a -> I i a ) - -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) ) + -> C ( ExtentOrder i ) + ( I i brushParams ) + ( Spline Closed () ( I i ( ℝ 2 ) ) ) ) -> Spline clo crvData ptData -> ST s @@ -398,7 +404,7 @@ computeStrokeOutline fitParams ptParams toBrushParams brushFn spline@( Spline { ori = splineOrientation brush0 fwdCond, bwdCond :: Bool ( fwdCond, bwdCond ) - | prevTgt `cross` tgt < 0 && prevTgt ^.^ tgt < 0 + | prevTgt × tgt < 0 && prevTgt ^.^ tgt < 0 = ( isJust $ between ori prevTgtFwd tgtFwd testTgt1 , isJust $ between ori prevTgtBwd tgtBwd ( -1 *^ testTgt1 ) ) @@ -448,14 +454,11 @@ outlineFunction , Interpolatable ( 𝕀 Double ) ( 𝕀 usedParams ) , DiffInterp 'Point brushParams , DiffInterp 'Interval brushParams - , HasChainRule Double 2 usedParams - , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 usedParams ) - , HasChainRule Double 2 brushParams - , HasChainRule ( 𝕀 Double ) 2 ( 𝕀 brushParams ) - , Traversable ( D 2 brushParams ) - - -- , Diffy Double usedParams - -- , Diffy ( 𝕀 Double ) ( 𝕀 usedParams ) + , HasChainRule Double ( ExtentOrder 'Point ) usedParams + , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 usedParams ) + , HasChainRule Double ( ExtentOrder 'Point ) brushParams + , HasChainRule ( 𝕀 Double ) ( ExtentOrder 'Interval ) ( 𝕀 brushParams ) + , Traversable ( D ( ExtentOrder 'Point ) brushParams ) -- Debugging. , Show ptData, Show brushParams @@ -466,15 +469,19 @@ outlineFunction . DiffInterp i brushParams => Proxy# i -> ( forall a. a -> I i a ) - -> I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) ) + -> C ( ExtentOrder i ) + ( I i brushParams ) + ( Spline Closed () ( I i ( ℝ 2 ) ) ) ) -> ptData -> Curve Open crvData ptData -> OutlineFn outlineFunction ptParams toBrushParams brushFromParams sp0 crv = let - pathAndUsedParams :: forall i - . ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) + pathAndUsedParams :: forall i k arr + . ( k ~ ExtentOrder i, CurveOrder k + , arr ~ C k + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) , Module ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) , Torsor ( T ( I i ( ℝ 2 ) ) ) ( I i ( ℝ 2 ) ) @@ -482,34 +489,37 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = , Torsor ( T ( I i usedParams ) ) ( I i usedParams ) ) => ( forall a. a -> I i a ) - -> ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ), I i ( ℝ 1 ) ~> I i usedParams ) + -> ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ), I i ( ℝ 1 ) `arr` I i usedParams ) pathAndUsedParams toI = case crv of LineTo { curveEnd = NextPoint sp1 } | let seg = Segment sp0 sp1 - -> ( line @i ( fmap ( toI . coords ) seg ) - , line @i ( fmap ( toI . ptParams ) seg ) ) + -> ( line @k @i ( fmap ( toI . coords ) seg ) + , line @k @i ( fmap ( toI . ptParams ) seg ) ) Bezier2To { controlPoint = sp1, curveEnd = NextPoint sp2 } | let bez2 = Quadratic.Bezier sp0 sp1 sp2 - -> ( bezier2 @i ( fmap ( toI . coords ) bez2 ) - , bezier2 @i ( fmap ( toI . ptParams ) bez2 ) ) + -> ( bezier2 @k @i ( fmap ( toI . coords ) bez2 ) + , bezier2 @k @i ( fmap ( toI . ptParams ) bez2 ) ) Bezier3To { controlPoint1 = sp1, controlPoint2 = sp2, curveEnd = NextPoint sp3 } | let bez3 = Cubic.Bezier sp0 sp1 sp2 sp3 - -> ( bezier3 @i ( fmap ( toI . coords ) bez3 ) - , bezier3 @i ( fmap ( toI . ptParams ) bez3 ) ) + -> ( bezier3 @k @i ( fmap ( toI . coords ) bez3 ) + , bezier3 @k @i ( fmap ( toI . ptParams ) bez3 ) ) - usedParams :: ℝ 1 ~> usedParams - path :: ℝ 1 ~> ℝ 2 + usedParams :: C ( ExtentOrder 'Point ) ( ℝ 1 ) usedParams + path :: C ( ExtentOrder 'Point ) ( ℝ 1 ) ( ℝ 2 ) ( path, usedParams ) = pathAndUsedParams @Point id curvesI :: 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) - curvesI = brushStrokeData @'Interval @brushParams + curvesI = brushStrokeData @'Interval @( ExtentOrder 'Interval ) @brushParams pathI - ( chainRule @( 𝕀 Double ) @2 usedParamsI $ linear ( nonDecreasing toBrushParams ) ) + ( chainRule @( 𝕀 Double ) @( ExtentOrder 'Interval ) + usedParamsI + ( linear ( nonDecreasing toBrushParams ) ) + ) ( brushFromParams @'Interval proxy# singleton ) - usedParamsI :: 𝕀ℝ 1 ~> 𝕀 usedParams - pathI :: 𝕀ℝ 1 ~> 𝕀ℝ 2 + usedParamsI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀 usedParams ) + pathI :: C ( ExtentOrder 'Interval ) ( 𝕀ℝ 1 ) ( 𝕀ℝ 2 ) ( pathI, usedParamsI ) = pathAndUsedParams @'Interval singleton fwdBwd :: OutlineFn @@ -520,9 +530,12 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = where curves :: Seq ( ℝ 1 -> StrokeDatum Point ) - curves = brushStrokeData @Point @brushParams + curves = brushStrokeData @Point @( ExtentOrder 'Point ) @brushParams path - ( chainRule @Double usedParams $ linear toBrushParams ) + ( chainRule @Double @( ExtentOrder 'Point ) + usedParams + ( linear toBrushParams ) + ) ( brushFromParams @Point proxy# id ) t @@ -537,11 +550,11 @@ outlineFunction ptParams toBrushParams brushFromParams sp0 crv = bisSols = bisection 0.0001 curvesI - in trace - ( unlines $ - ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) : - "" : - map show bisSols ) + in --trace + -- ( unlines $ + -- ( "bisectionMethod: #(possible zeroes) = " ++ show ( length bisSols ) ) : + -- "" : + -- map show bisSols ) fwdBwd ----------------------------------- @@ -799,9 +812,9 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) bez :: Cubic.Bezier ( ℝ 2 ) bez = Cubic.Bezier {..} c01, c12, c23 :: Double - c01 = tgt_wanted `cross` tgt0 - c12 = tgt_wanted `cross` tgt1 - c23 = tgt_wanted `cross` tgt2 + c01 = tgt_wanted × tgt0 + c12 = tgt_wanted × tgt1 + c23 = tgt_wanted × tgt2 correctTangentParam :: Double -> Maybe Double correctTangentParam t | t > -epsilon && t < 1 + epsilon @@ -823,109 +836,15 @@ withTangent tgt_wanted spline@( Spline { splineStart } ) , offset = T $ Cubic.bezier @( T ( ℝ 2 ) ) bez t } --------------------------------------------------------------------------------- --- | A brush stroke, as described by the equation --- --- \[ c(t,s) = p(t) + b(t,s) \] --- --- where: --- --- - \( p(t) \) is the path that the brush follows, and --- - \( b(t,s) \) is the brush shape, as it varies along the path. -brushStroke :: Module r ( T v ) - => D 2 ( ℝ 1 ) v -- ^ stroke path \( p(t) \) - -> D 2 ( ℝ 2 ) v -- ^ brush \( b(t,s) \) - -> D 2 ( ℝ 2 ) v -brushStroke ( D21 p dpdt d2pdt2 ) ( D22 b dbdt dbds d2bdt2 d2bdtds d2bds2 ) = - D22 ( unT $ T p ^+^ T b ) - -- c = p + b - ( dpdt ^+^ dbdt ) dbds - -- ∂c/∂t = dp/dt + ∂b/∂t - -- ∂c/∂s = ∂b/∂s - - ( d2pdt2 ^+^ d2bdt2 ) d2bdtds d2bds2 - -- ∂²c/∂t² = d²p/dt² + ∂²b/∂t² - -- ∂²c/∂t∂s = ∂²b/∂t∂s - -- ∂²c/∂s² = ∂²b/∂s² - --- | The envelope equation --- --- \[ E = \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} = 0, \] --- --- together with the total derivative --- --- \[ \frac{\mathrm{d} c}{\mathrm{d} t}, \] --- --- and the partial derivatives --- --- \[ \frac{\partial E}{\partial s}, \qquad \frac{\partial E}{\partial s}. \] --- --- NB: if \( \frac{\partial E}{\partial s} \) is zero, the total derivative is ill-defined. -envelopeEquation :: forall i - . ( D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 ) - , Cross ( I i Double ) ( T ( I i ( ℝ 2 ) ) ) - , Fractional ( I i Double ) - ) - => D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) - -> ( I i Double, T ( I i ( ℝ 2 ) ), T ( I i ( ℝ 2 ) ), I i Double, I i Double ) -envelopeEquation ( D22 _ dcdt dcds d2cdt2 d2cdtds d2cds2 ) = - let ee = dcdt `cross` dcds - dEdt = d2cdt2 `cross` dcds + dcdt `cross` d2cdtds - dEds = d2cdtds `cross` dcds + dcdt `cross` d2cds2 - tot = dcdt -- ^-^ ( dEdt / dEds ) *^ dcds - dEdsTot = dEds *^ dcdt ^-^ dEdt *^ dcds - in ( ee, tot, dEdsTot, dEdt, dEds ) - -- Computation of total derivative dc/dt: - -- - -- dc/dt = ∂c/∂t + ∂c/∂s ∂s/∂t - -- ∂s/∂t = - ∂E/∂t / ∂E/∂s - -- - -- ∂E/∂s dc/dt = ∂E/∂s ∂c/∂t - ∂E/∂t ∂c/∂s. - --- | Linear interpolation, as a differentiable function. -line :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) - , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) - ) - => Segment b -> I i ( ℝ 1 ) ~> b -line ( Segment a b ) = D \ ( coerce -> t ) -> - D21 ( lerp @( T b ) t a b ) - ( a --> b ) - origin - --- | A quadratic Bézier curve, as a differentiable function. -bezier2 :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) - , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) - ) - => Quadratic.Bezier b -> I i ( ℝ 1 ) ~> b -bezier2 bez = D \ ( coerce -> t ) -> - D21 ( Quadratic.bezier @( T b ) bez t ) - ( Quadratic.bezier' bez t ) - ( Quadratic.bezier'' bez ) - --- | A cubic Bézier curve, as a differentiable function. -bezier3 :: forall i b - . ( Module ( I i Double ) ( T b ), Torsor ( T b ) b - , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) - , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) - ) - => Cubic.Bezier b -> I i ( ℝ 1 ) ~> b -bezier3 bez = D \ ( coerce -> t ) -> - D21 ( Cubic.bezier @( T b ) bez t ) - ( Cubic.bezier' bez t ) - ( Cubic.bezier'' bez t ) - -splineCurveFns :: forall i - . ( D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) +splineCurveFns :: forall i k + . ( k ~ ExtentOrder i, CurveOrder 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 ) ) , Coercible ( I i ( ℝ 1 ) ) ( I i Double ) ) - => Spline Closed () ( I i ( ℝ 2 ) ) -> Seq ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) + => Spline Closed () ( I i ( ℝ 2 ) ) -> Seq ( C k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) ) splineCurveFns spls = runIdentity . bifoldSpline @@ -936,14 +855,14 @@ splineCurveFns spls where curveFn :: I i ( ℝ 2 ) -> Curve Open () ( I i ( ℝ 2 ) ) - -> ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) + -> ( C k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) ) curveFn p0 = \case LineTo { curveEnd = NextPoint p1 } - -> line @i $ Segment p0 p1 + -> line @k @i $ Segment p0 p1 Bezier2To { controlPoint = p1, curveEnd = NextPoint p2 } - -> bezier2 @i $ Quadratic.Bezier p0 p1 p2 + -> bezier2 @k @i $ Quadratic.Bezier p0 p1 p2 Bezier3To { controlPoint1 = p1, controlPoint2 = p2, curveEnd = NextPoint p3 } - -> bezier3 @i $ Cubic.Bezier p0 p1 p2 p3 + -> bezier3 @k @i $ 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. @@ -1019,15 +938,21 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData Nothing -> ( False, initialGuess ) Just s0 -> ( True , s0 ) in case f ( ℝ1 s ) of -- TODO: a bit redundant to have to compute this again... - StrokeDatum { ee = _ee, dstroke, 𝛿E𝛿t = _𝛿E𝛿t, 𝛿E𝛿s, dcdt } -> + StrokeDatum + { dstroke + , ee = D12 ( ℝ1 _ee ) ( T ( ℝ1 _𝛿E𝛿t ) ) ( T ( ℝ1 ee_s ) ) + , 𝛿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 totDeriv - | abs 𝛿E𝛿s < epsilon + let dcdt + | abs ee_s < epsilon , let s' = if s >= 0.5 then s - 1e-9 else s + 1e-9 - = case f ( ℝ1 s' ) of { StrokeDatum { dcdt = dcdt_s' } -> dcdt_s' } + = case f ( ℝ1 s' ) of + StrokeDatum { ee = D12 _ _ ( T ( ℝ1 ee_s' ) ), 𝛿E𝛿sdcdt = D0 𝛿E𝛿sdcdt' } + -> recip ee_s' *^ 𝛿E𝛿sdcdt' | otherwise - = dcdt + = recip ee_s *^ 𝛿E𝛿sdcdt in --trace -- ( unlines -- [ "solveEnvelopeEquations" @@ -1036,17 +961,16 @@ solveEnvelopeEquations _t path_t path'_t ( fwdOffset, bwdOffset ) strokeData -- , " c = " ++ show dstroke -- , " E = " ++ show _ee -- , " ∂E/∂t = " ++ show _𝛿E𝛿t - -- , " ∂E/∂s = " ++ show 𝛿E𝛿s - -- , " dc/dt = " ++ show totDeriv + -- , " ∂E/∂s = " ++ show ee_s + -- , " dc/dt = " ++ show dcdt -- ] ) - ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke - , totDeriv ) + ( good, ℝ1 s, value @Double @2 @( ℝ 2 ) dstroke, dcdt ) eqn :: ( ℝ 1 -> StrokeDatum Point ) -> ( Double -> ( Double, Double ) ) eqn f s = case f ( ℝ1 s ) of - StrokeDatum { ee, 𝛿E𝛿s } -> - ( ee, 𝛿E𝛿s ) + StrokeDatum { ee = D12 ( ℝ1 ee ) _ ( T ( ℝ1 ee_s ) ) } -> + ( ee, ee_s ) maxIters :: Word maxIters = 5 --30 @@ -1061,48 +985,55 @@ instance Applicative ZipSeq where pure _ = error "only use Applicative ZipSeq with non-empty Traversable functors" liftA2 f ( ZipSeq xs ) ( ZipSeq ys ) = ZipSeq ( Seq.zipWith f xs ys ) -brushStrokeData :: forall i brushParams - . ( Differentiable i brushParams +brushStrokeData :: forall i k brushParams arr + . ( k ~ ExtentOrder i, CurveOrder k, arr ~ C k + , Differentiable i brushParams , Fractional ( I i Double ) - , D 2 ( I i ( ℝ 1 ) ) ~ D 2 ( ℝ 1 ) - , D 2 ( I i ( ℝ 2 ) ) ~ D 2 ( ℝ 2 ) + , HasChainRule ( I i Double ) k ( I i ( ℝ 1 ) ) + , Applicative ( D k ( ℝ 1 ) ) + + , D ( k - 2 ) ( I i ( ℝ 2 ) ) ~ D ( k - 2 ) ( ℝ 2 ) + , D ( k - 1 ) ( I i ( ℝ 2 ) ) ~ D ( k - 1 ) ( ℝ 2 ) + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , D k ( I i ( ℝ 1 ) ) ~ D k ( ℝ 1 ) + , 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 ( ℝ 1 ) ~> I i ( ℝ 2 ) ) + => ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) -- ^ path - -> ( I i ( ℝ 1 ) ~> I i brushParams ) + -> ( I i ( ℝ 1 ) `arr` I i brushParams ) -- ^ brush parameters - -> ( I i brushParams ~> Spline Closed () ( I i ( ℝ 2 ) ) ) + -> ( I i brushParams `arr` Spline Closed () ( I i ( ℝ 2 ) ) ) -- ^ brush from parameters -> ( I i ( ℝ 1 ) -> Seq ( I i ( ℝ 1 ) -> StrokeDatum i ) ) brushStrokeData path params brush = \ t -> let - dpath_t :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + dpath_t :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) !dpath_t = runD path t - dparams_t :: D 2 ( I i ( ℝ 1 ) ) ( I i brushParams ) - !dparams_t@( D21 { _D21_v = params_t } ) = runD params t - dbrush_params :: D 2 ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) - !dbrush_params = runD brush params_t - splines :: Seq ( D 2 ( I i brushParams ) ( I i ( ℝ 1 ) ~> I i ( ℝ 2 ) ) ) - !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i ) dbrush_params - dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) - !dbrushes_t = force $ fmap ( uncurryD2 . ( chain @(I i Double) @2 dparams_t ) ) splines + dparams_t :: D k ( I i ( ℝ 1 ) ) ( I i brushParams ) + !dparams_t = runD params t + dbrush_params :: D k ( I i brushParams ) ( Spline Closed () ( I i ( ℝ 2 ) ) ) + !dbrush_params = runD brush $ value @( I i Double ) @k @( I i ( ℝ 1 ) ) dparams_t + splines :: Seq ( D k ( I i brushParams ) ( I i ( ℝ 1 ) `arr` I i ( ℝ 2 ) ) ) + !splines = getZipSeq $ traverse ( ZipSeq . splineCurveFns @i @k ) dbrush_params + dbrushes_t :: Seq ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + !dbrushes_t = force $ fmap ( uncurryD @k . ( chain @(I i Double) @k dparams_t ) ) splines -- This is the crucial use of the chain rule. in fmap ( mkStrokeDatum dpath_t ) dbrushes_t where - mkStrokeDatum :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -> ( I i ( ℝ 1 ) -> D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) + mkStrokeDatum :: D k ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) + -> ( I i ( ℝ 1 ) -> D k ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) ) -> ( I i ( ℝ 1 ) -> StrokeDatum i ) mkStrokeDatum dpath_t dbrush_t s = let dbrush_t_s = dbrush_t s - dstroke@( D22 _c _𝛿c𝛿t _𝛿c𝛿s _ _ _ ) = brushStroke dpath_t dbrush_t_s - ( ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s ) = envelopeEquation @i dstroke + dstroke = brushStroke @k dpath_t dbrush_t_s + ( ee, 𝛿E𝛿sdcdt ) = envelopeEquation @k @i dstroke in -- trace -- ( unlines -- [ "envelopeEquation:" @@ -1113,62 +1044,29 @@ brushStrokeData path params brush = -- , " ∂c/∂s = " ++ show _𝛿c𝛿s -- , " E = " ++ show ee -- , " ∂E/∂t = " ++ show _𝛿E𝛿t - -- , " ∂E/∂s = " ++ show 𝛿E𝛿s + -- , " ∂E/∂s = " ++ show ee_s -- , " dc/dt = " ++ show dcdt ] ) $ StrokeDatum { dpath = dpath_t , dbrush = dbrush_t_s , dstroke - , ee, dcdt, 𝛿E𝛿sdcdt, 𝛿E𝛿t, 𝛿E𝛿s } - - --- | The value and derivative of a brush stroke at a given coordinate --- \( (t_0, s_0) \), together with the value of the envelope equation at that --- point. -data StrokeDatum i - = StrokeDatum - { -- | Path \( p(t_0) \) (with its 1st and 2nd derivatives). - dpath :: D 2 ( I i ( ℝ 1 ) ) ( I i ( ℝ 2 ) ) - -- | Brush shape \( b(t_0, s_0) \) (with its 1st and 2nd derivatives). - , dbrush :: D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) - - -- Everything below can be computed in terms of the first two fields. - - -- | Stroke \( c(t_0,s_0) = p(t_0) + b(t_0,s_0) \) (with its 1st and 2nd derivatives). - , dstroke :: D 2 ( I i ( ℝ 2 ) ) ( I i ( ℝ 2 ) ) - -- | Envelope - -- - -- \[ E(t_0,s_0) = \left ( \frac{\partial c}{\partial t} \times \frac{\partial c}{\partial s} \right )_{(t_0,s_0)}. \] - , ee :: I i Double - -- | \( \left ( \frac{\partial E}{\partial s} \right )_{(t_0,s_0)}. \) - , 𝛿E𝛿s :: I i Double - -- | \( \left ( \frac{\partial E}{\partial t} \right )_{(t_0,s_0)}. \) - , 𝛿E𝛿t :: I i Double - -- | Total derivative - -- - -- \[ \left ( \frac{\mathrm{d} c}{\mathrm{d} t} \right )_{(t_0,s_0)}. \] - -- - -- This is ill-defined when \( \frac{\partial E}{\partial s} = 0 \). - , dcdt, 𝛿E𝛿sdcdt :: T ( I i ( ℝ 2 ) ) - - - } - -deriving stock instance Show ( StrokeDatum 'Point ) -deriving stock instance Show ( StrokeDatum 'Interval ) + , ee + , 𝛿E𝛿sdcdt + } -------------------------------------------------------------------------------- bisection :: Double -> ( 𝕀ℝ 1 -> Seq ( 𝕀ℝ 1 -> StrokeDatum 'Interval ) ) - -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -bisection minWidth eqs = bisect initialCands [] [] + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀ℝ 1, 𝕀ℝ 2 ) ] +bisection minWidth eqs = + bisect initialCands [] [] where - bisect :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, need bisection to refine + bisect :: [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀ℝ 1, 𝕀ℝ 2 ) ] -- have solutions, need bisection to refine -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1 ) ] -- have been bisected, don't know if they contain solutions yet - -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] -- have solutions, don't bisect further - -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀 Double, 𝕀ℝ 2 ) ] + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀ℝ 1, 𝕀ℝ 2 ) ] -- have solutions, don't bisect further + -> [ ( 𝕀ℝ 1, Int, 𝕀ℝ 1, 𝕀ℝ 1, 𝕀ℝ 2 ) ] bisect [] [] sols = sols bisect cands ( ( t, i, s ) : toTry ) sols | Just ( ee, 𝛿E𝛿sdcdt ) <- isCand t i s @@ -1211,17 +1109,19 @@ bisection minWidth eqs = bisect initialCands [] [] [ (t, i, s, ee, 𝛿E𝛿sdcdt ) | let !eqs_t = eqs t , ( eq_t, i ) <- zip ( toList eqs_t ) ( [0,1..] :: [Int] ) - , let !( StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } ) = eq_t s - , Interval.inf ( ival ee ) < 0 && Interval.sup ( ival ee ) > 0 - , cmpℝ2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) - && cmpℝ2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + , let !( StrokeDatum { ee = D22 ee _ _ _ _ _, 𝛿E𝛿sdcdt = D12 ( T 𝛿E𝛿sdcdt ) _ _ } ) = eq_t s + , Interval.inf ( ival ee ) < Rounded ( ℝ1 0 ) + , Interval.sup ( ival ee ) > Rounded ( ℝ1 0 ) + , cmpℝ2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) + , cmpℝ2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) ] - isCand :: 𝕀ℝ 1 -> Int -> 𝕀ℝ 1 -> Maybe ( 𝕀 Double, 𝕀ℝ 2 ) + isCand :: 𝕀ℝ 1 -> Int -> 𝕀ℝ 1 -> Maybe ( 𝕀ℝ 1, 𝕀ℝ 2 ) isCand t i s = case ( ( eqs t ) `Seq.index` i ) s of - StrokeDatum { ee, 𝛿E𝛿sdcdt = T 𝛿E𝛿sdcdt } -> + StrokeDatum { ee = D22 ee _ _ _ _ _, 𝛿E𝛿sdcdt = D12 ( T 𝛿E𝛿sdcdt ) _ _ } -> do guard $ - Interval.inf ( ival ee ) < 0 && Interval.sup ( ival ee ) > 0 + Interval.inf ( ival ee ) < Rounded ( ℝ1 0 ) + && Interval.sup ( ival ee ) > Rounded ( ℝ1 0 ) && cmpℝ2 (<) ( getRounded ( Interval.inf $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) && cmpℝ2 (>) ( getRounded ( Interval.sup $ ival 𝛿E𝛿sdcdt ) ) ( ℝ2 0 0 ) return ( ee, 𝛿E𝛿sdcdt ) diff --git a/src/splines/Math/Differentiable.hs b/src/splines/Math/Differentiable.hs index 9c77026..433f5b8 100644 --- a/src/splines/Math/Differentiable.hs +++ b/src/splines/Math/Differentiable.hs @@ -1,7 +1,7 @@ {-# LANGUAGE UndecidableInstances #-} module Math.Differentiable - ( Differentiable, DiffInterp ) + ( ExtentOrder, Differentiable, DiffInterp ) where -- base @@ -24,12 +24,10 @@ import Math.Ring type ExtentOrder :: Extent -> Nat type family ExtentOrder e where - ExtentOrder i = 2 - --ExtentOrder 'Point = 2 - --ExtentOrder 'Interval = 2 + ExtentOrder 'Point = 2 + ExtentOrder 'Interval = 3 -- Currently we're doing order 2 derivatives for the brush stroke fitting, -- but order 3 derivatives for the interval Newton method to find cusps. --- TODO: using 2 for both until migration finishes. type Differentiable :: Extent -> Type -> Constraint class diff --git a/src/splines/Math/Interval.hs b/src/splines/Math/Interval.hs index 8d6e706..84c0fc2 100644 --- a/src/splines/Math/Interval.hs +++ b/src/splines/Math/Interval.hs @@ -120,7 +120,7 @@ instance Inner ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where in x1x2 + y1y2 instance Cross ( 𝕀 Double ) ( T ( 𝕀ℝ 2 ) ) where - T ( 𝕀 ( ℝ2 x1_lo y1_lo ) ( ℝ2 x1_hi y1_hi ) ) `cross` + T ( 𝕀 ( ℝ2 x1_lo y1_lo ) ( ℝ2 x1_hi y1_hi ) ) × T ( 𝕀 ( ℝ2 x2_lo y2_lo ) ( ℝ2 x2_hi y2_hi ) ) = let !x1y2 = 𝕀 x1_lo x1_hi * 𝕀 y2_lo y2_hi !y2x1 = 𝕀 x2_lo x2_hi * 𝕀 y1_lo y1_hi diff --git a/src/splines/Math/Module.hs b/src/splines/Math/Module.hs index dea92d3..f5512dc 100644 --- a/src/splines/Math/Module.hs +++ b/src/splines/Math/Module.hs @@ -63,7 +63,7 @@ lerp :: forall v r p. ( Module r v, Torsor v p ) => r -> p -> p -> p lerp t p0 p1 = ( t *^ ( p0 --> p1 :: v ) ) • p0 class Module r m => Cross r m where - cross :: m -> m -> r + (×) :: m -> m -> r -- | Norm of a vector, computed using the inner product. norm :: forall m r. ( Floating r, Inner r m ) => m -> r @@ -122,7 +122,7 @@ instance Inner Double ( T ( ℝ 2 ) ) where V2 x1 y1 ^.^ V2 x2 y2 = x1 Ring.* x2 + y1 Ring.* y2 instance Cross Double ( T ( ℝ 2 ) ) where - cross ( V2 x1 y1 ) ( V2 x2 y2 ) = x1 Ring.* y2 Ring.- x2 Ring.* y1 + V2 x1 y1 × V2 x2 y2 = x1 Ring.* y2 Ring.- x2 Ring.* y1 -- | Compute whether two vectors point in the same direction, -- that is, whether each vector is a (strictly) positive multiple of the other. @@ -130,8 +130,8 @@ instance Cross Double ( T ( ℝ 2 ) ) where -- Returns @False@ if either of the vectors is zero. strictlyParallel :: T ( ℝ 2 ) -> T ( ℝ 2 ) -> 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) + = abs ( u × 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 @. -- @@ -159,8 +159,8 @@ convexCombination v0 v1 u where c0, c10 :: Double - c0 = v0 `cross` u - c10 = ( v0 ^-^ v1 ) `cross` u + c0 = v0 × u + c10 = ( v0 ^-^ v1 ) × u -------------------------------------------------------------------------------- -- Not sure how to set things up to automate the following... diff --git a/src/splines/Math/Orientation.hs b/src/splines/Math/Orientation.hs index 91d1adf..81719fe 100644 --- a/src/splines/Math/Orientation.hs +++ b/src/splines/Math/Orientation.hs @@ -32,7 +32,7 @@ import Data.Generics.Internal.VL import Math.Epsilon ( nearZero ) import Math.Module - ( cross ) + ( (×) ) import Math.Bezier.Spline ( Spline(..), Curves(..), Curve(..), NextPoint(..) , SplineType(..), KnownSplineType(..), SSplineType(..) @@ -63,7 +63,7 @@ convexOrientation ( v1 : v2 : vs ) = CW where crossProduct :: Double - crossProduct = v1 `cross` v2 + crossProduct = v1 × v2 convexOrientation _ = CCW -- default -- | Compute the orientation of a spline, assuming tangent vectors have a monotone angle.