Fix default monoidal annotation for closed 2D paths

[fryguybob]

(Imported from http://code.google.com/p/diagrams/issues/detail?id=17. Original issue from byor...@gmail.com on April 2, 2011, 06:13:46 PM UTC)

The stroke and strokeT methods create diagrams from paths and trails respectively. Currently the monoidal annotation they assign is (const (Any False)). However, for closed 2D paths the default annotation ought to test whether a given point lies inside the path. This is non-trivial but well-studied so there ought to be some existing algorithms that can be easily ported.

This is quite critical since pretty much every shape is ultimately based on paths, and without this the whole monoidal annotation scheme is quite useless.

[fryguybob]

(Imported. Original comment by byor...@gmail.com on April 9, 2011, 12:22:23 AM UTC)

http://paulbourke.net/geometry/insidepoly/ may be useful.

[fryguybob]

(Imported. Original comment by fryguy...@gmail.com on April 25, 2011, 02:13:08 PM UTC)

Do we also need to know what "fill rule" is going to be used to decide?

http://hackage.haskell.org/packages/archive/cairo/latest/doc/html/Graphics-Rendering-Cairo.html#t:FillRule

Another reference:

http://softsurfer.com/Archive/algorithm_0103/algorithm_0103.htm

[fryguybob]

(Imported. Original comment by fryguy...@gmail.com on April 25, 2011, 03:01:15 PM UTC)

It didn't verify this much, but the simple tests work:

{{{

!hs

import Data.VectorSpace import Control.Applicative

cross (x,y) (x',y') = x * y' - y * x'

segments = zip <*> tail

isInsideWinding p = (/= 0) . crossings p isInsideEvenOdd p = odd . crossings p

crossings p@(,y) = sum . map test . segments where
test l@((,ay),(_,by)) | ay <= y = if by <= y then isLeft l else 0 | by <= y = negate (isLeft l) | otherwise = 0

isLeft (a,b) = floor $ signum (cross (b ^-^ a) (p ^-^ a))

example :: [(Double,Double)] example = [(0,0),(4,0),(4,3),(1,3),(1,2),(3,2),(3,1),(2,1),(2,4),(0,4),(0,0)]

cases :: [(Double,Double)] cases = [(1,1),(5,5),(1.5,2.5),(2.5,1.5)]

testW = map (isInsideWinding example) cases == [True,False,True,False] testE = map (isInsideEvenOdd example) cases == [True,False,False,False] }}}

[fryguybob]

(Imported. Original comment by byor...@gmail.com on April 25, 2011, 04:30:15 PM UTC)

This is a great start, thanks! We'll also need to take Bezier segments into account which makes things a bit more interesting. But in principle it shouldn't be too hard to compute intersection points with Bezier curves, figure out whether they are "up" or "down", and so on, although I suppose caution is needed in dealing with points where Bezier curves self-intersect.

As for the fill rule, I suppose we should make this an attribute of paths (like open/closed is now)?

[fryguybob]

(Imported. Original comment by fryguy...@gmail.com on April 25, 2011, 09:11:02 PM UTC)

Yeah fill rule would make sense as an attribute. Humm, I know you can get nine intersection points out of just two non-self intersecting beziers (I don't know if that is the upper-bound or not). Maybe something like subdividing at horizontal inflection points would simplify things. For now we could just use bezier interpolation. I wrote this somewhat related code a while ago:

{{{

!hs

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

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 }}}

[fryguybob]

(Imported. Original comment by byor...@gmail.com on May 4, 2011, 07:55:04 PM UTC)

Whew, it works! Mostly. I think. Someone still needs to generalize it so the user can choose the fill rule to use (I'll add another ticket for this), and there are probably still lingering numerical stability problems with the cubic solver, but I'm happy to finally have this working. =)