From 15d50f8d76e3ff815569aaae972539fd4b45a89b Mon Sep 17 00:00:00 2001
From: sheaf <sam.derbyshire@gmail.com>
Date: Sun, 16 May 2021 00:33:58 +0200
Subject: [PATCH] add signed curvature computation

---
 src/lib/Math/Bezier/Cubic.hs     | 18 +++++++++++++-----
 src/lib/Math/Bezier/Quadratic.hs | 15 +++++++++++++--
 src/lib/Math/Module.hs           |  6 +++++-
 3 files changed, 31 insertions(+), 8 deletions(-)

diff --git a/src/lib/Math/Bezier/Cubic.hs b/src/lib/Math/Bezier/Cubic.hs
index 930fe21..0f392a8 100644
--- a/src/lib/Math/Bezier/Cubic.hs
+++ b/src/lib/Math/Bezier/Cubic.hs
@@ -18,7 +18,7 @@ module Math.Bezier.Cubic
   ( Bezier(..)
   , fromQuadratic
   , bezier, bezier', bezier''
-  , curvature, squaredCurvature
+  , curvature, squaredCurvature, signedCurvature
   , subdivide
   , ddist, closestPoint
   , drag, selfIntersectionParameters
@@ -70,12 +70,13 @@ import Math.Epsilon
 import Math.Module
   ( Module (..)
   , lerp
-  , Inner(..), squaredNorm
+  , Inner(..), norm, squaredNorm
+  , cross
   )
 import Math.Roots
   ( realRoots, solveQuadratic )
 import Math.Vector2D
-  ( Point2D(..) )
+  ( Point2D(..), Vector2D(..) )
 
 --------------------------------------------------------------------------------
 
@@ -125,11 +126,11 @@ bezier'' ( Bezier {..} ) t
       ( p1 --> p0 ^+^ p1 --> p2 )
       ( p2 --> p1 ^+^ p2 --> p3 )
 
--- | Curvature of a quadratic Bézier curve.
+-- | 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
 
--- | Square of curvature of a quadratic Bézier curve.
+-- | Square of curvature of a cubic Bézier curve.
 squaredCurvature :: forall v r p. ( Torsor v p, Inner r v, RealFloat r ) => Bezier p -> r -> r
 squaredCurvature bez t
   | sq_nm_g' < epsilon
@@ -144,6 +145,13 @@ squaredCurvature bez t
     sq_nm_g' :: r
     sq_nm_g' = squaredNorm @v g'
 
+-- | Signed curvature of a planar cubic Bézier curve.
+signedCurvature :: forall r. Floating r => Bezier ( Point2D r ) -> r -> r
+signedCurvature bez t = ( g' `cross` g'' ) / norm g'
+  where
+    g', g'' :: Vector2D r
+    g'  = bezier'  @( Vector2D r ) bez t
+    g'' = bezier'' @( Vector2D r ) bez t
 
 -- | Subdivide a cubic Bézier curve into two parts.
 subdivide :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> ( Bezier p, Bezier p )
diff --git a/src/lib/Math/Bezier/Quadratic.hs b/src/lib/Math/Bezier/Quadratic.hs
index 30e4b8d..f3f7f76 100644
--- a/src/lib/Math/Bezier/Quadratic.hs
+++ b/src/lib/Math/Bezier/Quadratic.hs
@@ -15,7 +15,7 @@
 module Math.Bezier.Quadratic
   ( Bezier(..)
   , bezier, bezier', bezier''
-  , curvature, squaredCurvature
+  , curvature, squaredCurvature, signedCurvature
   , subdivide
   , ddist, closestPoint
   , interpolate
@@ -65,10 +65,13 @@ import Math.Epsilon
 import Math.Module
   ( Module (..)
   , lerp
-  , Inner(..), squaredNorm
+  , Inner(..), norm, squaredNorm
+  , cross
   )
 import Math.Roots
   ( realRoots )
+import Math.Vector2D
+  ( Point2D(..), Vector2D(..) )
 
 --------------------------------------------------------------------------------
 
@@ -119,6 +122,14 @@ squaredCurvature bez t
     sq_nm_g' :: r
     sq_nm_g' = squaredNorm @v g'
 
+-- | Signed curvature of a planar quadratic Bézier curve.
+signedCurvature :: forall r. Floating r => Bezier ( Point2D r ) -> r -> r
+signedCurvature bez t = ( g' `cross` g'' ) / norm g'
+  where
+    g', g'' :: Vector2D r
+    g'  = bezier'  @( Vector2D r ) bez t
+    g'' = bezier'' @( Vector2D r ) bez
+
 -- | Subdivide a quadratic Bézier curve into two parts.
 subdivide :: forall v r p. ( Torsor v p, Module r v ) => Bezier p -> r -> ( Bezier p, Bezier p )
 subdivide ( Bezier {..} ) t = ( Bezier p0 q1 pt, Bezier pt r1 p2 )
diff --git a/src/lib/Math/Module.hs b/src/lib/Math/Module.hs
index 7e8ac89..dad2727 100644
--- a/src/lib/Math/Module.hs
+++ b/src/lib/Math/Module.hs
@@ -9,7 +9,7 @@
 module Math.Module
   ( Module(..), lerp
   , Inner(..)
-  , squaredNorm, quadrance, distance
+  , norm, squaredNorm, quadrance, distance
   , proj, projC, closestPointOnSegment
   , cross
   , strictlyParallel, convexCombination
@@ -75,6 +75,10 @@ infixl 8 ^.^
 class Module r m => Inner r m where
   (^.^) :: m -> m -> r
 
+-- | Norm of a vector, computed using the inner product.
+norm :: forall m r. ( Floating r, Inner r m ) => m -> r
+norm = sqrt . squaredNorm
+
 -- | Squared norm of a vector, computed using the inner product.
 squaredNorm :: forall m r. Inner r m => m -> r
 squaredNorm v = v ^.^ v