diagrams-lib-1.3: src/Diagrams/TwoD/Curvature.hs
{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE GADTs #-}
-----------------------------------------------------------------------------
-- |
-- Module : Diagrams.TwoD.Curvature
-- Copyright : (c) 2013 diagrams-lib team (see LICENSE)
-- License : BSD-style (see LICENSE)
-- Maintainer : diagrams-discuss@googlegroups.com
--
-- Compute curvature for segments in two dimensions.
--
-----------------------------------------------------------------------------
module Diagrams.TwoD.Curvature
(
curvature
, radiusOfCurvature
, squaredCurvature
, squaredRadiusOfCurvature
) where
import Control.Lens (over)
import Control.Monad
import Data.Monoid.Inf
import Diagrams.Segment
import Diagrams.Tangent
import Diagrams.TwoD.Types
import Linear.Vector
-- | Curvature measures how curved the segment is at a point. One intuition
-- for the concept is how much you would turn the wheel when driving a car
-- along the curve. When the wheel is held straight there is zero curvature.
-- When turning a corner to the left we will have positive curvature. When
-- turning to the right we will have negative curvature.
--
-- Another way to measure this idea is to find the largest circle that we can
-- push up against the curve and have it touch (locally) at exactly the point
-- and not cross the curve. This is a tangent circle. The radius of that
-- circle is the \"Radius of Curvature\" and it is the reciprocal of curvature.
-- Note that if the circle is on the \"left\" of the curve, we have a positive
-- radius, and if it is to the right we have a negative radius. Straight
-- segments have an infinite radius which leads us to our representation. We
-- result in a pair of numerator and denominator so we can include infinity and
-- zero for both the radius and the curvature.
--
--
-- Lets consider the following curve:
--
-- <<diagrams/src_Diagrams_TwoD_Curvature_diagramA.svg#diagram=diagramA&height=200&width=400>>
--
-- The curve starts with positive curvature,
--
-- <<diagrams/src_Diagrams_TwoD_Curvature_diagramPos.svg#diagram=diagramPos&height=200&width=400>>
--
-- approaches zero curvature
--
-- <<diagrams/src_Diagrams_TwoD_Curvature_diagramZero.svg#diagram=diagramZero&height=200&width=400>>
--
-- then has negative curvature
--
-- <<diagrams/src_Diagrams_TwoD_Curvature_diagramNeg.svg#diagram=diagramNeg&height=200&width=400>>
--
-- > {-# LANGUAGE GADTs #-}
-- >
-- > import Diagrams.TwoD.Curvature
-- > import Data.Monoid.Inf
-- > import Diagrams.Coordinates
-- >
-- > segmentA :: Segment Closed V2 Double
-- > segmentA = Cubic (12 ^& 0) (8 ^& 10) (OffsetClosed (20 ^& 8))
-- >
-- > curveA = lw thick . strokeP . fromSegments $ [segmentA]
-- >
-- > diagramA = pad 1.1 . centerXY $ curveA
-- >
-- > diagramPos = diagramWithRadius 0.2
-- >
-- > diagramZero = diagramWithRadius 0.45
-- >
-- > diagramNeg = diagramWithRadius 0.8
-- >
-- > diagramWithRadius t = pad 1.1 . centerXY
-- > $ curveA
-- > <> showCurvature segmentA t
-- > # withEnvelope (curveA :: D V2 Double)
-- > # lc red
-- >
-- > showCurvature :: Segment Closed V2 Double -> Double -> Diagram SVG
-- > showCurvature bez@(Cubic b c (OffsetClosed d)) t
-- > | v == (0,0) = mempty
-- > | otherwise = go (radiusOfCurvature bez t)
-- > where
-- > v@(x,y) = unr2 $ firstDerivative b c d t
-- > vp = (-y) ^& x
-- >
-- > firstDerivative b c d t = let tt = t*t in (3*(3*tt-4*t+1))*^b + (3*(2-3*t)*t)*^c + (3*tt)*^d
-- >
-- > go Infinity = mempty
-- > go (Finite r) = (circle (abs r) # translate vpr
-- > <> strokeP (origin ~~ (origin .+^ vpr)))
-- > # moveTo (origin .+^ atParam bez t)
-- > where
-- > vpr = signorm vp ^* r
-- >
--
curvature :: RealFloat n
=> Segment Closed V2 n -- ^ Segment to measure on.
-> n -- ^ Parameter to measure at.
-> PosInf n -- ^ Result is a @PosInf@ value where @PosInfty@ represents
-- infinite curvature or zero radius of curvature.
curvature s = toPosInf . over _y sqrt . curvaturePair s
-- | With @squaredCurvature@ we can compute values in spaces that do not support
-- 'sqrt' and it is just as useful for relative ordering of curvatures or looking
-- for zeros.
squaredCurvature :: RealFloat n => Segment Closed V2 n -> n -> PosInf n
squaredCurvature s = toPosInf . over _x (join (*)) . curvaturePair s
-- | Reciprocal of @curvature@.
radiusOfCurvature :: RealFloat n
=> Segment Closed V2 n -- ^ Segment to measure on.
-> n -- ^ Parameter to measure at.
-> PosInf n -- ^ Result is a @PosInf@ value where @PosInfty@ represents
-- infinite radius of curvature or zero curvature.
radiusOfCurvature s = toPosInf . (\(V2 p q) -> V2 (sqrt q) p) . curvaturePair s
-- | Reciprocal of @squaredCurvature@
squaredRadiusOfCurvature :: RealFloat n => Segment Closed V2 n -> n -> PosInf n
squaredRadiusOfCurvature s = toPosInf . (\(V2 p q) -> (V2 q (p * p))) . curvaturePair s
-- Package up problematic values with the appropriate infinity.
toPosInf :: RealFloat a => V2 a -> PosInf a
toPosInf (V2 _ 0) = Infinity
toPosInf (V2 p q)
| isInfinite r || isNaN r = Infinity
| otherwise = Finite r
where r = p / q
-- Internal function that is not quite curvature or squaredCurvature but lets
-- us get there by either taking the square root of the numerator or squaring
-- the denominator respectively.
curvaturePair :: Num n
=> Segment Closed V2 n -> n -> V2 n
curvaturePair (Linear _) _ = V2 0 1 -- Linear segments always have zero curvature (infinite radius).
curvaturePair seg@(Cubic b c (OffsetClosed d)) t
= V2 (x'*y'' - y'*x'') ((x'*x' + y'*y')^(3 :: Int))
where
(V2 x' y' ) = seg `tangentAtParam` t
(V2 x'' y'') = secondDerivative
secondDerivative = (6*(3*t-2))*^b ^+^ (6-18*t)*^c ^+^ (6*t)*^d
-- TODO: We should be able to generalize this to higher dimensions. See
-- <http://en.wikipedia.org/wiki/Curvature>
--
-- TODO: I'm not sure what the best way to generalize squaredCurvature to other spaces is.
-- curvaturePair :: (Num t, Num (Scalar t), VectorSpace t)
-- => Segment Closed (t, t) -> Scalar t -> (t, t)
-- curvaturePair (Linear _) _ = (0,1) -- Linear segments always have zero curvature (infinite radius).
-- curvaturePair seg@(Cubic b c (OffsetClosed d)) t = ((x'*y'' - y'*x''), (x'*x' + y'*y')^(3 :: Integer))
-- where
-- (x' ,y' ) = seg `tangentAtParam` t
-- (x'',y'') = secondDerivative
-- secondDerivative = (6*(3*t-2))*^b ^+^ (6-18*t)*^c ^+^ (6*t)*^d