Offsetting Cubic segments.

offsetCubicSegment.hs

{-# LANGUAGE NoMonomorphismRestriction, ViewPatterns #-}
import Diagrams.Prelude
import Diagrams.Backend.Postscript.CmdLine
import Diagrams.Solve
import Diagrams.TwoD.Curvature

import Data.Monoid.PosInf


import qualified Debug.Trace as T

-- The basic plan here is to subdivide the segment until we have a reasonable approximation
-- of an arc, we can then scale offset handle lengths by the ratio of radii.  This means we
-- also need to approximate the radius of a segment.  One possibility for determining the
-- depth of recursion is to produce the offset curve, then compare points at some t-values.
-- If we have recursed far enough the distance between each pair with the same t-value
-- should be r.  There should be some good upper bound on the number of points we need to
-- check that follows from both curves being cubic.
--
-- There are some strange corner cases and I'm not sure how bad they are.  If there is some
-- sub-curve that we come to that approximates an arc with radius matching the offset radius
-- then we "should" end up with a degenerate sub-curve in the offset.  Really we should just
-- skip this segment, but that does mean that the result will not always be C^1 continuous.
-- On the other hand, we can't represent an arc exactly so that may not be precisely the
-- conditions that give us the degenerate curve.
--
offsetCubicSegment :: Double -> Double -> Segment R2 -> (Point R2, Trail R2)
offsetCubicSegment epsilon r s@(Cubic a b c) = (origin .+^ va, Trail (go (radiusOfCurvature s 0.5)) False)
  where
    -- Perpendiculars to handles.
    va = r *^ vperp (0 - a)
    vc = r *^ vperp (b - c)
    -- Split segments.
    ss = (\(a,b) -> [a,b]) $ splitAtParam s 0.5
    subdivided = concatMap (trailSegments . snd . offsetCubicSegment epsilon r) ss

    -- Offset with handles scaled based on curvature.
    offset factor = Cubic (a^*factor) ((b - c)^*factor + c + vc - va) (c + vc - va)


    -- We observe a corner.  Subdivide right away.
    go (Finite 0) = subdivided
    -- Some curvrature
    go roc
      | close     = [o]
      | otherwise = subdivided
      where
        -- We want the mulitplicative factor that takes us from the original
        -- segment's radius of curvature roc, to roc + r.
        --
        -- r + sr = x * sr
        --
        o = offset $ case roc of
              PosInfty  -> 1
              Finite sr -> 1 - r / sr -- TODO: I think my r's are backwards.

        close = and [epsilon > (magnitude (p o + va - p s - pp s))
--                    | t <- [0.01, 0.25, 0.5, 0.75, 0.99]
                    | t <- [0.25, 0.5, 0.75]
                    , let p = (`atParam` t)
                    , let pp = (`perpAtParam` t)
                    ]

---------------------------------------------------------------------

vperp :: R2 -> R2
vperp v = rotateBy (-1/4) (normalized v)

perpAtParam :: Segment R2 -> Double -> R2
perpAtParam s@(Cubic _ _ _) t = vperp (-a)
  where
    (Cubic a _ _) = snd $ splitAtParam s t

fromFixed' :: [FixedSegment R2] -> (Point R2, Trail R2)
fromFixed' [] = (p2 zeroV, Trail [] False) -- ???
fromFixed' (s:ss) = (a, Trail (b : map (snd . rel) ss) False)
  where
    (a, b) = rel s

rel (FLinear a b) = (a, Linear $ b .-. a)
rel (FCubic a b c d) = (a, Cubic (b .-. a) (c .-. a) (d .-. a))

---------------------------------------------------------------------

showExample :: Segment R2 -> Diagram Postscript R2
showExample s = pad 1.1 . centerXY $ d # lc blue # lw 0.1 <> d' # lw 0.1
  where
--    d  = stroke $ Path [(origin, Trail [s] False)]
    d  = mconcat . map (f blue) $ explodeTrail origin (Trail [s] False)
    d' = mconcat . zipWith f colors . uncurry explodeTrail $ offsetCubicSegment 0.1 1 s

    f c p@(Path [(a, Trail [Cubic vb vc vd] False)])
       = lw 0.01 (stroke (a          ~~ (a .+^ vb)))
      <> lw 0.01 (stroke ((a .+^ vc) ~~ (a .+^ vd)))
      <> (lc c . stroke $ p)
    f c p = lc c . stroke $ p

    colors = cycle [green, red]

showExample' :: Segment R2 -> Diagram Postscript R2
showExample' s = pad 1.1 . centerXY $ d # lc blue # lw 0.1 <> d' # lw 0.1
  where
      d  = stroke $ Path [(origin, Trail [s] False)]
      d' = mconcat . zipWith lc colors . map stroke . uncurry explodeTrail $ offsetCubicSegment 0.1 1 s

      colors = cycle [green, red]

----------------------------------------------------------------------

example :: Diagram Postscript R2
example = hcat . map showExample' $
        [ Cubic (10 &  0) (  5  & 18) (10 & 20)
        , Cubic ( 0 & 20) ( 10  & 10) ( 5 & 10)
        , Cubic (10 & 20) (  0  & 10) (10 &  0)
        , Cubic (10 & 20) ((-5) & 10) (10 &  0)
        ]

main = defaultMain example