[diagrams] Segment splitting and arc length

[Ryan Yates]

Hi,

I was looking through the code and saw a TODO for computing arc length of a
segment and for splitting segments and thought I would share some of the
successes I have had in this area.  For splitting at a parameter Bézier
curves couldn't be nicer (I'm pretty new to Haskell so I'm sure this code
could be nicer):

segmentSplitAtParam :: (InnerSpace v) => Segment v -> Scalar v -> (Segment
v, Segment v)
segmentSplitAtParam (Linear x1) t = (left, right)
      where left  = Linear p
            right = Linear (x1 ^-^ p)
            p = lerp zeroV x1 t
segmentSplitAtParam (Cubic c1 c2 x2) t = (left, right)
      where left  = Cubic a b e
            right = Cubic (c ^-^ e) (d ^-^ e) (x2 ^-^ e)
            p = lerp c1    c2 t
            a = lerp zeroV c1 t
            b = lerp a     p  t
            d = lerp c2    x2 t
            c = lerp p     d  t
            e = lerp b     c  t

Wikipedia has some nice animations of this:
http://en.wikipedia.org/wiki/Bézier_curve#Constructing_B.C3.A9zier_curves<http://en.wikipedia.org/wiki/B%C3%A9zier_curve#Constructing_B.C3.A9zier_curves>.
 You can write AtParam in terms of segmentSplitAtParam but I'm not sure if
that gives better numeric results or not.  Arc length (in my experience) is
the bane of Bézier curves.  What I have done in the past is to write ArcLength
:: Bezier -> Double -> Double where the second argument is the desired error
bound.  I would then recurse over splitting and use the distance to the end
point as a lower bound and the sum of the length of the distances of zeroV
to c1, c1 to c2, and c2 to x2.  I think it follows from the subdivision
definition that this is an upper bound, but I'm not completely sure
(particularly in spaces above R^2).

arcLength :: (InnerSpace v, Floating (Scalar v), Ord (Scalar v))
          => Segment v -> Scalar v -> Scalar v
arcLength (Linear x1) _ = magnitude x1
arcLength s@(Cubic c1 c2 x2) m
  | ub - lb < m = (ub + lb) / 2
  | otherwise   = arcLength l m + arcLength r m
 where (l,r) = segmentSplitAtParam s 0.5
       ub    = sum (map magnitude [c1, c2 ^-^ c1, x2 ^-^ c2])
       lb    = magnitude x2

I'm sure this implementation can be improved upon, but one nice thing is I
think it works fine with degenerate Béziers.  It will probably fail for
small enough values of m (but isn't floating arithmetic always a fail at
some point?).

Anyway, hope this is useful,

Ryan Yates
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.haskell.org/pipermail/diagrams/attachments/20100528/317e95d0/attachment.htm