add cubic Bézier crunode computation

This commit is contained in:
sheaf 2021-05-15 21:24:58 +02:00
parent 9f16cff978
commit df2469487d

View file

@ -21,7 +21,7 @@ module Math.Bezier.Cubic
, curvature, squaredCurvature
, subdivide
, ddist, closestPoint
, drag
, drag, selfIntersectionParameters
)
where
@ -73,7 +73,9 @@ import Math.Module
, Inner(..), squaredNorm
)
import Math.Roots
( realRoots )
( realRoots, solveQuadratic )
import Math.Vector2D
( Point2D(..) )
--------------------------------------------------------------------------------
@ -111,11 +113,9 @@ bezier ( Bezier {..} ) t =
-- | Derivative of cubic Bézier curve.
bezier' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
bezier' ( Bezier {..} ) t
bezier' ( Bezier {..} )
= ( 3 *^ )
$ lerp @v t
( lerp @v t ( p0 --> p1 ) ( p1 --> p2 ) )
( lerp @v t ( p1 --> p2 ) ( p2 --> p3 ) )
. Quadratic.bezier @v ( Quadratic.Bezier ( p0 --> p1 ) ( p1 --> p2 ) ( p2 --> p3 ) )
-- | Second derivative of a cubic Bézier curve.
bezier'' :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> v
@ -179,7 +179,7 @@ ddist ( Bezier {..} ) c = [ a5, a4, a3, a2, a1, a0 ]
closestPoint
:: forall v r p. ( Torsor v p, Inner r v, RealFloat r, Prim r, NFData r )
=> Bezier p -> p -> ArgMin r ( r, p )
closestPoint pts c = pickClosest ( 0 :| 1 : roots ) -- todo: also include the self-intersection point if one exists
closestPoint pts c = pickClosest ( 0 :| 1 : roots )
where
roots :: [ r ]
roots = filter ( \ r -> r > 0 && r < 1 ) ( realRoots 50 $ ddist @v pts c )
@ -222,3 +222,29 @@ drag ( Bezier {..} ) t q = Bezier { p0, p1 = p1', p2 = p2', p3 }
p1', p2' :: p
p1' = ( ( 1 - t ) *^ delta ) p1
p2' = ( t *^ delta ) p2
-- | Compute parameter values for the self-intersection of a planar cubic Bézier curve, if such exist.
--
-- The parameter values might lie outside the interval [0,1],
-- indicating a self-intersection of the extended curve.
--
-- Formula taken from:
-- "A Basis for the Implicit Representation of Planar Rational Cubic Bézier Curves"
-- Oliver J. D. Barrowclough, 2016
selfIntersectionParameters :: forall r. RealFloat r => Bezier ( Point2D r ) -> [ r ]
selfIntersectionParameters ( Bezier {..} ) = solveQuadratic c0 c1 c2
where
areaConstant :: Point2D r -> Point2D r -> Point2D r -> r
areaConstant ( Point2D x1 y1 ) ( Point2D x2 y2 ) ( Point2D x3 y3 ) =
x1 * ( y2 - y3 ) + x2 * ( y3 - y1 ) + x3 * ( y1 - y2 )
l0, l1, l2, l3, f1, f2, f3, c0, c1, c2 :: r
l0 = areaConstant p3 p2 p1
l1 = areaConstant p2 p3 p0
l2 = areaConstant p1 p0 p3
l3 = areaConstant p0 p1 p2
f1 = 3 * ( l1 * l1 - 3 * l0 * l2 )
f2 = 3 * ( l2 * l2 - 3 * l1 * l3 )
f3 = 3 * ( 9 * l0 * l3 - l1 * l2 )
c0 = f2
c1 = f3 - 2 * f2
c2 = f1 + f2 - f3