segmentSegment can fail to terminate on very specific and seemingly innocuous inputs

[bacchanalia]

Slight perturbations of d and r avoid the bug, as does increasing e (default 1e-8).

{-# LANGUAGE FlexibleContexts #-}
{-# LANGUAGE TypeFamilies #-}
import Diagrams.Prelude
import Diagrams.TwoD.Segment

d = 72       -- ~ 71.9 to 72.1
r = d/64     -- ~ 64 to 64.0000001
e = 6.969e-8 -- ~ <= 6.969e-8

path :: Path V2 Double
path = circle r # translateY d

trails :: [Located (Trail V2 Double)]
trails = head $ explodePath path

(s0:s1:_) = map (head . fixTrail) trails

bad = segmentSegment e s0 s1

-- Does not terminate or produce output.
main = print bad

[fryguybob]

I can reproduce the issue. Tracking it down further now.

[fryguybob]

We can make bezierClip in Diagrams.TwoD.Segment much more robust. It tries to narrow down potential intersections by intersecting a "fat line" of a portion of one curve with the convex hull of part of the other. In this case, there is no viable end condition as the end points of the two original segments match. There is no handling of the case where we have subdivided to a zero length line. It may be the case that the cubicbezier package has improved this code, but otherwise, we need to at least make lineEquation not produce NaN and bezierClip detect when the t and u intervals are equal (this might be as simple as changing the order of the cases, but that might have other implications).

[fryguybob]

It looks like there are significant changes in the cubicbezier package, but it
isn't immediately clear how they map to the code we have.

Here are some of the issues I see in Diagrams.TwoD.Segment as it is right now:

• The lineEquation function does not work if the two points are the same or very close.

  • diagrams-lib/src/Diagrams/TwoD/Segment.hs

    Lines 305 to 307 in 1f863bf

    	a = a' / d
    	b = b' / d
    	c = -(x1*a' + y1*b') / d

  • It is used by lineDistance. Logically that function should fallback
    to just being the distance to the point, but given it's type, we cannot
    check for NaN from lineEquation. Do we every rely on Floating not
    being RealFloat? How safe are zero checks?
  • We could avoid the division in lineEquation altogether in some cases
    or do it once, so it should probably return four values.
• Similar potential problems exist with intersectPt:

  • diagrams-lib/src/Diagrams/TwoD/Segment.hs

    Lines 319 to 321 in 1f863bf

    	intersectPt :: OrderedField n => n -> P2 n -> P2 n -> n
    	intersectPt d (P (V2 x1 y1)) (P (V2 x2 y2)) =
    	x1 + (d - y1) * (x2 - x1) / (y2 - y1)

  • The only code that uses this is chopHull and it does some (not
    sufficient I don't think) checking. It should be moved closer to those
    assumptions and the assumptions checked.
• The inRange function uses defEps and is used by functions that have an
  epsilon argument.
  • The semantics of the epsilons are not documented in general here. Are
    they in the vector space or the parameter space?

A broader issue with segmentSegment is the behavior when curves overlap. I
think we might be able do a reasonable job with that situation, but we need to
be really clear what the semantics are. For instance, there are Bezier curves
that have parts that overlap, don't overlap, and intersect at a single point.
For example:

import Diagrams.Prelude
import Diagrams.Backend.SVG.CmdLine

b = fromSegments [bezier3 (V2 4 0) (V2 1 3) (V2 1 (-1))] :: Located (Trail V2 Double)
b1 = stroke $ section b 0.0 0.6 :: Diagram B
b2 = stroke $ section b 0.4 1.0 :: Diagram B

main = mainWith . pad 1.1 . lw 7.0 $ mconcat
           [ b1 # lcA (blue `withOpacity` 0.5)
           , b2 # lcA (red `withOpacity` 0.5)
           ]

I think for this example it is clear what we would ideally want, segmentSegment to result in a list of values from a sum type that separates individual intersections from overlaps. What is less clear is the situation where there is a large part of the curve that "overlaps within epsilon". We probably need separate functions for the different meanings people might have.

[bacchanalia]

I think "overlaps within epsilon" aught to have to same behaviour as precise overlap because "precise" is going to be often be fuzzy do to floating point math anyway and otherwise we run into the problem of sectionSection trying to return a huge set of points along the overlap.

[fryguybob]

I imagine having two functions: segmentSegmentPoints and
segmentSegmentOverlaps for lack of better names. Both would have the same
result type:

data Intersection n = IPoint n n | IOverlap (n,n) (n,n)

segmentSegmentPoints :: OrderedField n
  => n -- epsilon in the parameter space
  -> FixedSegment V2 n
  -> FixedSegment V2 n
  -> [Intersection n]

segmentSegmentOverlaps :: OrderedField n
  => n -- epsilon in V2 Euclidean space
  -> FixedSegment V2 n
  -> FixedSegment V2 n
  -> [Intersection n]

The Points version will still include "exact" overlaps when they happen, but
will only include up to the maximum number of possible intersections for the
rest of the curve guided by the epsilon in parameter space (two intersections
within epsilon of each other may count as one). The Overlaps version would
try to find inexact overlaps guided by an epsilon in Euclidean space (as long
a stretch where the curves are within epsilon of each other).