```{-# LANGUAGE GADTs #-}
-----------------------------------------------------------------------------
-- |
-- Module      :  Diagrams.TwoD.Offset
--
-- Compute offsets to segments in two dimensions.
--
-----------------------------------------------------------------------------
module Diagrams.TwoD.Offset
( offsetSegment
) where

import Data.AffineSpace
import Data.Monoid.Inf
import Data.VectorSpace

import Diagrams.Core

import Diagrams.Parametric
import Diagrams.Path
import Diagrams.Segment
import Diagrams.Trail
import Diagrams.TwoD.Curvature
import Diagrams.TwoD.Transform
import Diagrams.TwoD.Types

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

unitPerp :: R2 -> R2
unitPerp = normalized . perp

perpAtParam :: Segment Closed R2 -> Double -> R2
perpAtParam   (Linear (OffsetClosed a))  t = unitPerp a
perpAtParam s@(Cubic _ _ _)              t = unitPerp a
where
(Cubic a _ _) = snd \$ splitAtParam s t

-- | Compute the offset of a segment.  Given a segment compute the offset
--   curve that is a fixed distance from the original curve.  For linear
--   segments nothing special happens, the same linear segment is returned
--   with a point that is offset by a perpendicular vector of the given offset
--   length.
--
--   Cubic segments require a search for a subdivision of cubic segments that
--   gives an approximation of the offset within the given epsilon tolerance.
--   We must do this because the offset of a cubic is not a cubic itself (the
--   degree of the curve increases).  Cubics do, however, approach constant
--   curvature as we subdivide.  In light of this we scale the handles of
--   the offset cubic segment in proportion to the radius of curvature difference
--   between the original subsegment and the offset which will have a radius
--   increased by the offset parameter.
--
--   In the following example the blue lines are the original segments and
--   the alternating green and red lines are the resulting offset trail segments.
--
--   <<diagrams/cubicOffsetExample.svg#diagram=cubicOffsetExample&width=600>>
--
--   Note that when the original curve has a cusp, the offset curve forms a
--   radius around the cusp, and when there is a loop in the original curve,
--   there can be two cusps in the offset curve.
--
offsetSegment :: Double     -- ^ Epsilon value that represents the maximum
--   allowed deviation from the true offset.  In
--   the current implementation each result segment
--   should be bounded by arcs that are plus or
--   minus epsilon from the radius of curvature of
--   the offset.
-> Double     -- ^ Offset from the original segment, positive is
--   on the right of the curve, negative is on the
--   left.
-> Segment Closed R2 -- ^ Original segment
-> (Point R2, Trail R2) -- ^ Resulting offset point and trail.
offsetSegment _       r s@(Linear (OffsetClosed a))    = (origin .+^ va, trailFromSegments [s])
where va = r *^ unitPerp a

offsetSegment epsilon r s@(Cubic a b (OffsetClosed c)) = (origin .+^ va, t)
where
t = trailFromSegments (go (radiusOfCurvature s 0.5))
-- Perpendiculars to handles.
va = r *^ unitPerp a
vc = r *^ unitPerp (c - b)
-- Split segments.
ss = (\(a,b) -> [a,b]) \$ splitAtParam s 0.5
subdivided = concatMap (trailSegments . snd . offsetSegment epsilon r) ss

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

-- We observe a corner.  Subdivide right away.
go (Finite 0) = subdivided
-- We have some curvature
go roc
| close     = [o]
| otherwise = subdivided
where
-- We want the multiplicative 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
Infinity  -> 1          -- Do the right thing.
Finite sr -> 1 + r / sr

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

-- > import Diagrams.TwoD.Offset
-- >
-- > showExample :: Segment R2 -> Diagram SVG 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
-- >          \$ offsetSegment 0.1 (-1) s
-- >
-- >       colors = cycle [green, red]
-- >
-- > cubicOffsetExample :: Diagram SVG R2
-- > cubicOffsetExample = 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)
-- >         ]
```