{-# 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: -- -- <> -- -- The curve starts with positive curvature, -- -- <> -- -- approaches zero curvature -- -- <> -- -- then has negative curvature -- -- <> -- -- > {-# 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 -- -- -- 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