Search Unity

Showcase Need a good 2D line segment intersector? Grab it while it's hot.

Discussion in 'Scripting' started by orionsyndrome, Sep 27, 2022.

  1. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    Code is tl;dr at the bottom

    It's a relatively simple thing (for most of us anyway) but here I've solved it somewhat unorthodoxly. And I'm thinking wow this should come bundled with Unity, with how complete it turned out, at least in 2D.

    Let's try something different: 2D line intersection via homogeneous coordinates.

    Hey I'm not an expert actually, but I know about homogeneous coordinates from working with barycentric coordinates, which are also homogeneous (you get them in HitInfo from Raycast for example). This is a pretty opaque topic in (common) geometry, one I won't bother you too much with (and psh I don't know what I'm talking about anyway), but what I know is that homogeneous coordinates live in a different space compared to Cartesian, where they nearly always satisfy some internally binding rule (i.e. you can usually scale them uniformly and this won't change their meaning, which is kind of a weird "zooming" property for a coordinate system). This gives them some unique properties which is why they were invented in the first place.

    Anyway, turns out you can intersect lines in 2D by using them and it's as easy as peanuts. It's not so easy to find an example though, so I had to scour through numpy freaks and AI research (sarcasm: yummy, my favorite activities). But after about a day, I got it working without too much hassle, and I wish you all to never have to go through that.

    So what it's all about? Well, you start by defining two lines with 2D points (a1, a2) and (b1, b2). So four points in total. Then you augment the coordinates by appending 1 as the third coordinate, as if they were in 3D, and whoa, we just got our homogeneous coordinates. Sounds so stupid and useless.

    Then you cross the point pairs to define two lines in space. That's how you do it, don't ask me why. And then you cross the lines again. Again, don't ask m.. because that's obviously a method for intersecting lines, alright, we ok now?
    And what you got in the end is the point of intersection in homogeneous space.
    Code (csharp):
    1. static Vector2 intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2) {
    2.   return (Vector2)(hp(a1).Cross(hp(a2))).Cross(hp(b1).Cross(hp(b2)));
    3. }
    4.  
    5. // hp = homogeneous point, erm or whatever it is
    6. static Vector3 hp(Vector2 pt) => new Vector3(pt.x, pt.y, 1f);
    Wait, WHAT?

    "You're kidding right? That can't possibly work"
    You're right, it can't. This is all a bad joke, I'm sorry. UNTIL

    you realize you need to scale this heavenly result down to proper legit down-to-earth glorious space we all know and love, by using that magical third component.
    Code (csharp):
    1. static Vector2 intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2) {
    2.   Vector3 cc = (hp(a1).Cross(hp(a2))).Cross(hp(b1).Cross(hp(b2)));
    3.   return ((Vector2)cc) / cc.z;
    4. }
    Man, this is so crazy! I love it.

    Btw I'm using my own libraries to type this down, and this Cross is the same as using Vector3.Cross, I just like a.Cross(b) better than Vector3.Cross(a, b). But a cross is a cross.

    If I did my math correctly (edit: I did not, but I fixed it), this is what actually happens (after simplifying multiplications with 1)
    Code (csharp):
    1. X: (ha2.x - ha1.x) * (hb1.x * hb2.y - hb1.y * hb2.x) - (hb2.x - hb1.x) * (ha1.x * ha2.y - ha1.y * ha2.x),
    2. Y: (hb1.y - hb2.y) * (ha1.x * ha2.y - ha1.y * ha2.x) - (ha1.y - ha2.y) * (hb1.x * hb2.y - hb1.y * hb2.x),
    3. Z: (ha1.y - ha2.y) * (hb2.x - hb1.x) - (ha2.x - ha1.x) * (hb1.y - hb2.y)
    which actually looks sensible. (You can also use this if you want the heart of it not to include any function calls, for performance.)

    It turns out that if Z ends up being 0, there is no intersection either due to lines being parallel or coincident (in 2D there are no skewed lines).

    You can find thorough explanations why this works (geometrically; numerically it looks the same as the algebraic solving-the-line-equations solution, only more robust). Here's one of the explanations with a picture if any of it makes any sense to you (it doesn't to me, I'm like "awkay" and it's all very awkward, don't be like me).

    Anyway, that was part I of this adventure. Let's now build a proper thing to round this off. Intersect methods typically return a short answer to a question "is there intersection at all?" and then shamelessly and unceremoniously plop the actual point sideways. Like so
    Code (csharp):
    1. static bool intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2, out Vector2 intersect) {
    2.   intersect = something;
    3.   return success;
    4. }
    Let's now try to make a segment intersection with the weird homogeneous stuff, because that's usually much more useful than having infinite lines intersect wherever.

    You can optimize segment intersection slightly by checking the bounding boxes first, and cancel the whole thing if the boxes do not overlap. But let's ignore the Bounds struct, that's for 3D, we can use Rect instead. The Rect is produced from the minimum corner (the one toward origin) and size. How can we easily determine what corner is the minimum for an arbitrary segment, especially if we consider all kinds of wild orientations for its two points?

    Well, we can find the component-wise minmax. This is one way to do it
    Code (csharp):
    1. static Rect rectFromSeg(Vector2 a, Vector2 b) {
    2.   var min = Vectx.LeastOf(a, b);
    3.   return new Rect(min, Vectx.MostOf(a, b) - min);
    4. }
    Pretty clean. LeastOf, just to illustrate the point, is functionally equivalent to
    Code (csharp):
    1. static Vector2 LeastOf(Vector2 a, Vector2 b) => new Vector2(Math.Min(a.x, b.x), Math.Min(a.y, b.y));
    Now we can attack the segment intersect method with more tools.
    Code (csharp):
    1. static bool intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2, out Vector2 intersect) {
    2.   intersect = new Vector2(float.NaN, float.NaN);
    3.  
    4.   var ar = rectFromSeg(a1, a2);
    5.   var br = rectFromSeg(b1, b2);
    6.   if(!ar.Overlaps(br)) return false; // quick and easy
    7.  
    8.   var cc = (hp(a1).Cross(hp(a2))).Cross(hp(b1).Cross(hp(b2)));
    9.   if(cc.z.Abs() < 1E-6f) return false; // skip the last piece (and avoid division by zero)
    10.  
    11.   intersect = ((Vector2)cc) * (1f / cc.z); // this replaces two divisions with one division and two multiplications
    12.   return true;
    13.  
    14.   // local functions
    15.   static Vector3 hp(Vector2 pt) => new Vector3(pt.x, pt.y, 1f);
    16. }
    But this doesn't work yet. Sure it finds the intersection like it did before, but doesn't treat our segments as of finite length. It can happen that the boxes overlap even though the segments don't.

    We have to check whether the resulting point actually lies within the interval of both segments. To do this, we again >shift< inside the "segment space" and there we claim that the starting point is as good as 0, and that the end point is as good as 1. The point has to lie inside this interval. If it's less than 0 or greater than 1, it's outside of the interval, and not included in the segment. This suddenly smells like a good ol' linear interpolation.

    Indeed, that's the lerp (or more precisely LerpUnclamped). We can always pretend that we have two points on a line and lerp between them with some interpolant t and this will yield us some point on this imaginary line, where 0 lands on the start, and 1 lands on the end. However, in this case, we got it backwards, we already have the point... and we'd like to get our t back! Yet there is no Vector2.InverseLerp, only Mathf.InverseLerp, so we'll have to think out of the box a little to figure it all out.

    The test boils down to this surprisingly non-pseudo code
    Code (csharp):
    1. static bool test(Vector2 p, Vector2 a, Vector2 b)
    2.   => p.LerpedBetween(a, b).InRange(0f, 1f);
    'Is this point within the interval of 0 and 1 as interpolated between some points a and b?' yes/no

    Code (csharp):
    1. static bool intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2, out Vector2 intersect) {
    2.   intersect = new Vector2(float.NaN, float.NaN);
    3.  
    4.   var ar = rectFromSeg(a1, a2);
    5.   var br = rectFromSeg(b1, b2);
    6.   if(!ar.Overlaps(br)) return false;
    7.  
    8.   var cc = (hp(a1).Cross(hp(a2))).Cross(hp(b1).Cross(hp(b2)));
    9.   if(cc.z.Abs() < 1E-6f) return false;
    10.  
    11.   var x = ((Vector2)cc) * (1f / cc.z);
    12.  
    13.   // both segments must contain the point
    14.   if(!test(x, a1, a2) || !test(x, b1, b2)) return false;
    15.  
    16.   intersect = x;
    17.   return true;
    18.  
    19.   // local functions
    20.   static Vector3 hp(Vector2 pt) => new Vector3(pt.x, pt.y, 1f);
    21.   static bool test(Vector2 p, Vector2 a, Vector2 b) => p.LerpedBetween(a, b).InRange(0f, 1f);
    22. }
    But we ought to figure out how to lower down the test function before we can try this out.

    Here's how to think about inverse lerps when there are more dimensions than 1. You can pick any one dimension and apply the inverse lerp there and it'll work. For example just X. But here's a caveat: imagine a line going straight horizontally, for example, its start is at (-1, 1) and its end at (5, 1), if you pick its X dimension, you're fine, because there is so much meat to extrapolate from, because it goes from -1 to 5 so there is a proper change in its attitude; but if you pick its Y dimension, nothing happens, there is nothing to lerp in between, no change at all. So that's the thing, you can pick any dimension, as long as there is some change recorded in it, so why not pick the one that records the greatest change (aka the delta)? That's one way to ensure the best possible precision.

    Lerp in general can be thought of as a mix of two interpolated values, where one grows at the expense of another. Thus
    x = a(1-t) + bt
    , whereas an inverse of that would look like this
    Code (csharp):
    1. x = a(1-t) + bt // we want t out
    2. x = a(1-t) + bt // apply a
    3. x = a - at + bt // rearrange
    4. x = a + bt - at // put t in front
    5. x = a + t(b - a) // move a to the left, then swap
    6. t(b - a) = x - a // divide everything by (b - a)
    7. t = (x - a) / (b - a)
    Here it becomes obvious that when there is no delta between a and b, the solution is undefined, because we end up dividing by zero. Another reason to stay clear from it. However in this context, if zero division happens, t is supposed to be 0 (because we know a == b). In other words, it's pretty well defined, it's just unsure whether it's 0 or 1, because it's kind of both at the same time.

    In the end, we're supposed to do three things: 1) find the maximum delta, 2) apply
    (x - a) / (b - a)
    to find
    t
    , 3) check the interval. Here are the first two (I choose not to do needless
    if
    blocks if I can help it, but you do you)
    Code (csharp):
    1. int i = (b.x - a.x).Abs() < (b.y - a.y).Abs()? 1 : 0; // pick the best dimension
    2. var num = p[i] - a[i]; // inverselerp numerator: (x - a)
    3. var denom = b[i] - a[i]; // inverselerp denominator: (b - a)
    To check the interval, we're supposed to first divide num and denom, but let's do something else instead. If we know 0 <= t <= 1 must return true, then a) any negative result must be false [because that would go below 0), and b) a num quantity larger then denom quantity must be false (quantity here denotes that we don't care about the sign) [because that would go over 1]. Finally, we avoid having to divide (and having to check against 0 denom).

    The (b) part is easy, we can do
    Code (csharp):
    1. if(num.Abs() > denom.Abs()) return false;
    but the (a) part with the negative result is "tricky" because we have 4 possible configurations (+/+, +/-, -/- and -/+), 2 of which will give us a negative result while two of the same will cancel out and give a positive. Now writing them all down is horrible practice imho, so let's just check if the signs are different. We can compare the two so-called SignSteps.
    Code (csharp):
    1. static bool signStep(float v) => (v >= 0f? 1 : -1);
    2. static bool sameSigns(float a, float b) => signStep(a) == signStep(b);
    SignStep is similar to the ordinary Sign function, but lacks zero. This is very useful when you need a sign stepping increment (who wants an increment of 0?), but also in this case, where we don't actually care about zero, we just want to see if num and denom are bearing the same signs.

    But then, who needs SignSteps, they're only useful as a concept. Why not simplify to
    Code (csharp):
    1. a >= 0f == b >= 0f
    Now that's an interesting-looking expression.

    We can finish this off
    Code (csharp):
    1. static bool test(Vector2 p, Vector2 a, Vector2 b) {
    2.   int i = (b.x - a.x).Abs() < (b.y - a.y).Abs()? 1 : 0;
    3.   float n = p[i] - a[i], d = b[i] - a[i];
    4.   return n >= 0f == d >= 0f && n.Abs() <= d.Abs();
    5. }
    And that's it.
    .
    .
    .
    But what about rays? *cue VSauce music*

    Wouldn't it be nice if we had all three categories elegantly supported by this one method: lines, rays, and segments? In fact, I didn't even change the name of the method to properly reflect that it's now computing segments. Let's add this, for funzies and, I guess, completeness.
    Code (csharp):
    1. enum Mode {
    2.   Lines,
    3.   Rays,
    4.   RayLine,
    5.   RaySegment,
    6.   Segments
    7. }
    With rays, we want to test only if that inverse lerp leads to a negative result (< 0). This means we should split and selectively exclude the other (> 1) test.
    Code (csharp):
    1. // when bidi(rection) is true, we test in two directions
    2. static bool test(Vector2 p, Vector2 a, Vector2 b, bool bidi = false) {
    3.   int i = (b.x - a.x).Abs() < (b.y - a.y).Abs()? 1 : 0;
    4.   float n = p[i] - a[i], d = b[i] - a[i];
    5.   if(bidi && n.Abs() > d.Abs()) return false;
    6.   return n >= 0f == d >= 0f;
    7. }

    tl;dr
    I'm sorry about the code not being directly copy/pasteable, I will switch this to Unity API asap, I just don't have time to do it right now.
    Code (csharp):
    1. // compute intersects like a boss you most certainly are™
    2. static bool intersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2, out Vector2 intersect, Mode mode = Mode.Segments) {
    3.   intersect = new Vector2(float.NaN, float.NaN);
    4.  
    5.   if(mode == Mode.Segments) {
    6.     var ar = rectFromSeg(a1, a2);
    7.     var br = rectFromSeg(b1, b2);
    8.     if(!ar.Overlaps(br)) return false;
    9.   }
    10.  
    11.   var cc = (hp(a1).Cross(hp(a2))).Cross(hp(b1).Cross(hp(b2)));
    12.   if(cc.z.Abs() < 1E-6f) return false;
    13.  
    14.   var x = ((Vector2)cc) * (1f / cc.z);
    15.  
    16.   if(mode switch {
    17.     Mode.Rays => !test(x, a1, a2) || !test(x, b1, b2),
    18.     Mode.RayLine => !test(x, a1, a2),
    19.     Mode.RaySegment => !test(x, a1, a2) || !test(x, b1, b2, bidi: true),
    20.     Mode.Segments => !test(x, a1, a2, bidi: true) || !test(x, b1, b2, bidi: true),
    21.     _ => false // lines can't fail at this
    22.   }) return false;
    23.  
    24.   intersect = x;
    25.   return true;
    26.  
    27.   // local functions
    28.   static Vector3 hp(Vector2 p) => new Vector3(p.x, p.y, 1f);
    29.  
    30.   static bool test(Vector2 p, Vector2 a, Vector2 b, bool bidi = false) {
    31.     int i = (b.x - a.x).Abs() < (b.y - a.y).Abs()? 1 : 0;
    32.     float n = p[i] - a[i], d = b[i] - a[i];
    33.     if(bidi && n.Abs() > d.Abs()) return false;
    34.     return n >= 0f == d >= 0f;
    35.   }
    36.  
    37. }
    38.  
    39. static Rect rectFromSeg(Vector2 a, Vector2 b) {
    40.   var min = Vectx.LeastOf(a, b); // build your own, example is in the post
    41.   return new Rect(min, Vectx.MostOf(a, b) - min);
    42. }
    43.  
    44. enum Mode {
    45.   Lines,
    46.   Rays,
    47.   RayLine,
    48.   RaySegment,
    49.   Segments
    50. }
    Usage example
    Code (csharp):
    1. if(IntersectLines(a1, a2, b1, b2, out var point, IntersectMode.RaySegment)) {
    2.   drawGizmo(point);
    3. }

    "Can I like the post?"
    Only if you insist lol
     
    Last edited: Sep 28, 2022
  2. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    Performance considerations and extended usage.

    If you want to lower down the 'crossing heart' of this method, you can introduce the following change (in place of lines 11 and 12)
    Code (csharp):
    1. var z = (a1.y - a2.y) * (b2.x - b1.x) - (a2.x - a1.x) * (b1.y - b2.y);
    2. if(z.Abs() < 1E-6f) return false;
    3.  
    4. var ip = (1f / z) * new Vector2(
    5.   (a2.x - a1.x) * (b1.x * b2.y - b1.y * b2.x) - (b2.x - b1.x) * (a1.x * a2.y - a1.y * a2.x),
    6.   (b1.y - b2.y) * (a1.x * a2.y - a1.y * a2.x) - (a1.y - a2.y) * (b1.x * b2.y - b1.y * b2.x)
    7. );
    This way we do a minimum of work to test the waters before we go on with the full thing.

    The whole solution has ~20 subtractions, 16 multiplications, 1 division, has no trigonometry, and takes zero roots. Frankly with a method in this performance class, the biggest bottleneck is the actual branching. But here you should note that we've avoided all kinds of edge-case detection thanks to homogeneous coordinates, and that nearly all branching is highly predictable because it mostly relies on the input arguments themselves. The only place where this isn't true is the
    test
    itself, and that's pretty much unavoidable (unless you set the mode to Lines).

    All local functions are very likely to be inlined by the compiler, so they serve only to reduce the source clutter. On a Mono compiler set to Debug, accessing vectors via indexer used to be extremely slow, but this isn't the case any more in the last 3-4 years or so, especially with IL2CPP set to Release, I don't know what was the problem before or if anyone had even noticed it.

    You can also expand this to return a meaningful failure evaluation. For example, once you know z is (almost) 0, you can immediately tell that the lines must be parallel/collinear. Each point of failure can turn into a mini-query.

    This is one example where you can mask certain results
    Code (csharp):
    1. [System.Flags]
    2. public enum IntersectResult {
    3.   Intersect = 0,
    4.   NoOverlap = 0x1, // test bounds failure
    5.   JustParallel = 0x2, // lines are parallel if z turns out 0
    6.   Collinear = 0x6, // lines are collinear if they are parallel AND the space between them is too small
    7.   NotInIntervalA = 0x10, // if a (a1, a2) test fails
    8.   NotInIntervalB = 0x20, // if a (b1, b2) test fails; you must split the switch expression in two parts
    9.   Uncontained = 0x30
    10. }
    I will cover the collinearity test and a visual test-driver component later.

    You can now test for broader states, for example, if any interval was violated
    Code (csharp):
    1. var res = IntersectLines(a1, a2, b1, b2, out var point, Mode.Segments);
    2. if(((int)res & (int)IntersectResult.Uncontained) > 0) { // we can use this as a mask
    3.   var msg = "Lines intersect but the point is outside the following segments: ".
    4.   if(res.HasFlag(IntersectResult.NotInIntervalA)) msg += "A ";
    5.   if(res.HasFlag(IntersectResult.NotInIntervalB)) msg += "B";
    6.   Debug.Log(msg);
    7. }

    Finally, as an extra, here's a suite of functions that find the nearest point on line/ray/segment (and optionally spit out distance and direction)
    (a and b define the line/ray/segment, point is some arbitrary point in space)
    Code (csharp):
    1. using System.Runtime.CompilerServices;
    2.  
    3. [MethodImpl(MethodImplOptions.AggressiveInlining)]
    4. static private Vector2 nearestPoint_core(this Vector2 point, Vector2 a, Vector2 b, Func<float, float, float> limiter = null) {
    5.   var dir = a.Toward(b, out var mag);
    6.   var dot = a.To(point).Dot(dir);
    7.   dot = limiter is null? dot : limiter(dot, mag);
    8.   return a.Directed(dir, dot);
    9. }
    10.  
    11. static public Vector2 NearestPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance, out Vector2 dir, Func<float, float, float> limiter = null) {
    12.   var np = point.nearestPoint_core(a, b, limiter);
    13.   dir = point.Toward(np, out distance);
    14.   return np;
    15. }
    16.  
    17. static public Vector2 NearestPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance, Func<float, float, float> limiter = null) {
    18.   var np = point.nearestPoint_core(a, b, limiter);
    19.   point.Toward(np, out distance);
    20.   return np;
    21. }
    22.  
    23. static public Vector2 NearestPoint(this Vector2 point, Vector2 a, Vector2 b, Func<float, float, float> limiter = null)
    24.   => point.nearestPoint_core(a, b, limiter);
    25.  
    26. static public Vector2 NearestSegPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance, out Vector2 dir)
    27.   => point.NearestPoint(a, b, out distance, out dir, (l,d) => l.Clamp(0f, d) );
    28.  
    29. static public Vector2 NearestSegPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance)
    30.   => point.NearestPoint(a, b, out distance, out _, (l,d) => l.Clamp(0f, d) );
    31.  
    32. static public Vector2 NearestSegPoint(this Vector2 point, Vector2 a, Vector2 b)
    33.   => point.NearestPoint(a, b, out _, out _, (l,d) => l.Clamp(0f, d) );
    34.  
    35. static public Vector2 NearestRayPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance, out Vector2 dir)
    36.   => point.NearestPoint(a, b, out distance, out dir, (l,_) => l.MinimumOf(0f) );
    37.  
    38. static public Vector2 NearestRayPoint(this Vector2 point, Vector2 a, Vector2 b, out float distance)
    39.   => point.NearestPoint(a, b, out distance, out _, (l,_) => l.MinimumOf(0f) );
    40.  
    41. static public Vector2 NearestRayPoint(this Vector2 point, Vector2 a, Vector2 b)
    42.   => point.NearestPoint(a, b, out _, out _, (l,_) => l.MinimumOf(0f) );
    Again, my stupid libraries that stand in the way of proper copypasta, but here's a translation
    a.Dot(b) <=> Vector2.Dot(a, b) <=> a.x * b.x + a.y * b.y
    a.To(b) <=> b - a <=> (b.x - a.x, b.y - a.y)
    p.Directed(d, m) <=> p + m * d <=> (p.x + m * d.x, p.y + m * d.y)
    a.Toward(b) <=> (b - a).normalized
    a.magnitude <=> MathF.Sqrt(a.x * a.x + a.y * a.y)
    a.normalized <=> { var invd = 1f / a.magnitude; return (a.x * invd, a.y * invd) }
    a.Toward(b, out distance) <=> { var dlt = b - a; distance = dlt.magnitude; return dlt / distance; }
    v.MinimumOf(min) <=> Math.Max(min, v) <=> v < min? min : v
    v.Clamp(min, max) <=> Mathf.Clamp(v, min, max) <=> Math.Min(Math.Max(v, min), max)

    Edit: I messed up distance and direction computation, it's fixed now. Though it's less optimal now and I can see what pushed me to make this error in the first place. I'll optimize this in the future.
     
    Last edited: Sep 28, 2022
    Terraya and Bunny83 like this.
  3. Bunny83

    Bunny83

    Joined:
    Oct 18, 2010
    Posts:
    4,006
    homogeneous coordinates are really powerful indeed and once you get the gist, they aren't even that complicated. Sure in 3d (where we use 4d homogeneous coordinates) it's a bit difficult to "visualize" them. However the 2d analog works the same way. If you understand the 2d case, 3d should be no problem.
     
    orionsyndrome likes this.
  4. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    Let's tackle collinearity. So this is when two lines practically both lie on the exact same line.

    Not only the lines are parallel to each other, but they also seem to be the same thing. Why's that a problem? Well, it's not easy to tell if they're on the same line just by looking at the points they were defined with.

    You might have two completely disjointed segments, miles apart, that in fact, belong to the same line, and you want to be able to check whether that's true. Or you have two segments which overlap in such a manner that you should only keep one of them (or merge them together?). Or two rays looking away from each other, etc.

    So it's useful to be able to tell parallels from collinears apart. We already know if a line is parallel when z is 0. It makes sense now to add extra something to resolve whether it's collinear too, just in that case.

    For example
    Code (csharp):
    1. var z = (a1.y - a2.y) * (b2.x - b1.x) - (a2.x - a1.x) * (b1.y - b2.y);
    2.  
    3. if(z.Abs() < 1E-6f) {
    4.   if(collinear()) return IntersectResult.Collinear;
    5.   return IntersectResult.Parallel;
    6. }
    But how exactly to implement this?
    A regular analytical way would be to define a line as y = mx + b where m is the slope, and b is the y-intercept. Let's try this. A slope is "rise over run" so how much y changes per x step. This means that m = (a2.y - a1.y) / (a2.x - a1.x)
    Ohnoes, we got an explosive division here, one that blows up with perfectly vertical lines because, well, tan(90) is undefined. However, once again, it's not really undefined, it's just confused about the kind of infinity it should represent. So let's give it a bit of support, and it'll turn out fine.

    Code (csharp):
    1. static float slope(Vector2 a, Vector2 b)
    2.   => ((b.y - a.y) / (b.x - a.x)).SubstNaN( _ => (b.y - a.y) > 0f? float.PositiveInfinity : float.NegativeInfinity );
    This is one ad-hoc method one can use, to conditionally inject a fixed value in a non-pre-evaluated manner. What do I mean by this?
    Well this particular SubstNaN method is fitted with a lambda that is called only after we evaluate NaN. This means we evaluate that subtraction twice only when (b.x - a.x) == 0.
    Code (csharp):
    1. static public float SubstNaN(this float v, Func<float, float> subst) => float.IsNaN(v)? subst(v) : v;
    This way, we don't care about the result any more. It just works. Next we want to get that y-intercept, and this is why we needed the slope in the first place.
    Code (csharp):
    1. static float yIntercept(Vector2 a, Vector2 b) => a.y - slope(a, b) * a.x;
    Now if that slope blows-up to infinity, the y-intercept will blow up as well. And it does make sense for vertical lines to have undefined y-intercept.

    What was that thing we were after? Ah, yes, collinearity can be defined by saying that the two parallel lines must be at no distance apart. And this is why this is tricky, how do we tell the distance between two parallel lines?

    Well, if we know that they share the same slope m, here's one handy formula.
    Code (csharp):
    1. d = |(b2 - b1) / sqrt(1 + m^2)|
    And sure enough, we're equipped to do this
    Code (csharp):
    1. static float parallelDistance(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2)
    2.   => ( (yIntercept(b1, b2) - yIntercept(a1, a2)) / Mathx.Sqrt(1f + slope(a1, a2).Sqr()) ).Abs();
    3.  
    4. static bool areCollinear(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2) => parallelDistance(...) < 1E-5f;
    But it's not very good, in fact it's horrible. Edge cases might turn up like this
    Code (csharp):
    1. |(INF - INF) / sqrt(1 + INF^2)| => |0 / sqrt(INF)| => NaN
    and this
    Code (csharp):
    1. |(-INF - INF) / sqrt(1 + INF^2)| => |INF / sqrt(INF)| => NaN
    But you get the general idea behind it, you can see the costs involved, and how we can approach it. Now let's do something that works in line with homogeneous coordinates, not against them. Remember NearestPoint functions I listed in the previous post? Well you sample a point from one segment, and use that to get the distance to another, but treat that other segment as a line instead.

    And if you've seen that code for NearestPoint method, it basically has the same projective geometry heart, that dot lets us do the actual projection, and also there are some similarities with the analytical method, namely in both cases we can't get rid of the square root (y'know ... distance is a distance, it's always like this).

    The second overload is what we care about, and we're supposed to ignore the limiter.
    Edit: Had to rewrite this because of a stupid mistake I made, and I've decided to introduce a more optimal variant in this particular case. Edit2: and a similar mistake yet again grr, nvm practice makes perfect
    Code (csharp):
    1. static public Vector2 NearestLinePoint(Vector2 point, Vector2 a, Vector2 b, out float distance) {
    2.   var dir = a.Toward(b); // here lies the square root and division
    3.   var dot = a.To(point).Dot(dir); // just some subtractions, additions and multiplications
    4.   var np = a.Directed(dir, dot); // not much going on in this line either
    5.   distance = point.DistanceTo(np); // another square root and division, trying to optimize this away is why I'm making mistakes
    6.   return np;
    7. }
    // Edit: updated this as well
    Code (csharp):
    1. static public Vector2 NearestLinePoint(this Vector2 point, Vector2 a, Vector2 b, out float distance) {
    2.   var dir = (b - a).normalized;
    3.   var dot = Vector2.Dot(point - a, dir);
    4.   var np = a + dot * dir;
    5.   distance = Vector3.Distance(point, np);
    6.   return np;
    7. }

    So the collinear test isn't exactly cheap, but honestly, we only do it if we already decided to skip doing everything else. Regardless it's nice to include this test only as an option.

    So for example
    Code (csharp):
    1. [Flags]
    2. public enum IntersectOptions {
    3.   None = 0,
    4.   SegmentBoundsTest = 0x1,
    5.   DetectCollinearity = 0x2
    6. }
    7.  
    8. static public IntersectResult IntersectLines(Vector2 a1, Vector2 a2, Vector2 b1, Vector2 b2, out Vector2 intersect,
    9.   IntersectMode mode = IntersectMode.Segments, IntersectOptions options = IntersectOptions.SegmentBoundsTest) {
    10.   intersect = new Vector2(float.NaN, float.NaN);
    11.  
    12.   if(mode == IntersectMode.Segments && options.HasFlag(IntersectOptions.SegmentBoundsTest)) {
    13.     var ar = rectFromSeg(a1, a2);
    14.     var br = rectFromSeg(b1, b2);
    15.     if(!ar.Overlaps(br)) return IntersectResult.NoOverlap;
    16.   }
    17.  
    18.   var z = (a1.y - a2.y) * (b2.x - b1.x) - (a2.x - a1.x) * (b1.y - b2.y);
    19.  
    20.   if(z.Abs() < 1E-6f) {
    21.     if(options.HasFlag(IntersectOptions.DetectCollinearity)) {
    22.       a1.NearestLinePoint(b1, b2, out var distance);
    23.       if(distance < 1E-5f) return IntersectResult.Collinear;
    24.     }
    25.     return IntersectResult.Parallel;
    26.   }
    27.  
    28.   // ....
    29.  
    30. }
     
    Last edited: Sep 28, 2022
  5. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    Whoops, I messed up NearestPoints :)
    The core is ok, I just rushed the wrong answers to distance and direction.
    Fixing it. Could be better, but should be fixed. I'll try to revisit this at some point in the future.
     
    Last edited: Sep 28, 2022
  6. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    Thanks for the encouragement. Right now I don't need an intersector for 3D, but it does seem like it should be a walk in the park. I've tested this a lot and I quite like the overall stability of this, compared to my earlier solutions. I started to resent excessive branching for the edge cases in geometry, nowadays it just looks like a code smell to me, and there is always some NaN or a floating-point artifact lurking behind the corner. So this is part of an effort to upgrade the libraries for the 2D project I'm currently on. Plus I learned something new.

    I also managed to round up the Ellipse code I was working on before. It's quite a beast now. Also figured out how to cheaply tell if a point lies inside an arbitrary ellipse or not. And I'm not sure, the tangent solution isn't an exact one, but if I'd flatten this out, idk, it's very similar to Inigo's second sdf method. In any case, it's more than perfect for that celestial bodies experiment I did before.

    Edit:
    A 2D ellipse tutorial is located here.
     
    Last edited: Apr 30, 2023
  7. orionsyndrome

    orionsyndrome

    Joined:
    May 4, 2014
    Posts:
    3,116
    For the record: the solution from the post #1, modified so that it works with vanilla API, can be found here.
    At some point I might also upload the Segment class struct that I've developed with the aforementioned concepts, fully tested and optimized.