lineQuadratic() not finding one specific intersection

[hamoid]

Hi,

I'm working with intersections in OPENRNDR and I have the following design made out of 4 segments (two lines and two curves):

There should be 5 intersections detected but only 4 are found. The original program simplified to only use the two involved segments looks like this:

fun doTest() {
    val l = Line2.line(Vec2(512.0, 96.0), Vec2(128.0, 384.0))
    val c = Bezier2.curve(Vec2(128.0, 384.0), Vec2(320.0, 168.0), Vec2(512.0, 384.0)) as Bezier2.QuadraticBezier2?
    println(l)
    println(c)
    val i = Intersections.lineQuadratic(l, c)
    println(i.size)
    i.forEach { println(it) }
}

This is the output of the program:

a=[x=512.0, y=96.0], b=[x=128.0, y=384.0]
p0=[x=128.0, y=384.0], p1=[x=320.0, y=168.0], p2=[x=512.0, y=384.0]
0

Zero intersections are found.

This is a visualization of those vertices using Inkscape. Even if the blue line should be a curve, I think it shows that there should be an intersection:

Any ideas why this happens?

[hamoid]

The order of the vertices seems to affect the result, when it probably shouldn't.

fun doTest() {
    val l0 = Vec2(512.0, 96.0)
    val l1 = Vec2(128.0, 384.0)
    val c0 = Vec2(128.0, 384.0)
    val c1 = Vec2(320.0, 168.0)
    val c2 = Vec2(512.0, 384.0)

    // fails to find intersection
    println(Intersections.lineQuadratic(
        Line2.line(l0, l1),
        Bezier2.curve(c0, c1, c2) as Bezier2.QuadraticBezier2?
    ).size)

    // reversed line works
    println(Intersections.lineQuadratic(
        Line2.line(l1, l0),
        Bezier2.curve(c0, c1, c2) as Bezier2.QuadraticBezier2?
    ).size)

    // reversed curve works
    println(Intersections.lineQuadratic(
        Line2.line(l0, l1),
        Bezier2.curve(c2, c1, c0) as Bezier2.QuadraticBezier2?
    ).size)

    // both reversed works
    println(Intersections.lineQuadratic(
        Line2.line(l1, l0),
        Bezier2.curve(c2, c1, c0) as Bezier2.QuadraticBezier2?
    ).size)
}

Output:

0
2
2
2

Equations.solveQuadratic() returns zero results in the first case.

[hamoid]

Adding a small offset to one of the vertices also makes it returns the expected intersection.
For example replacing (128.0, 384.0) with (128.01, 384.03)

[hamoid]

Here a java program demonstrating the issue which can be dropped into src/io/lacuna/artifex

package io.lacuna.artifex;

import io.lacuna.artifex.utils.Intersections;

public class issue3 {
    public static void main(String[] args) {
        Vec2 l0 = new Vec2(512.0, 96.0);
        Vec2 l1 = new Vec2(128.0, 384.0);
        Vec2 c0 = new Vec2(128.0, 384.0);
        Vec2 c1 = new Vec2(320.0, 168.0);
        Vec2 c2 = new Vec2(512.0, 384.0);
        calculate(l0, l1, c0, c1, c2);
        calculate(l1, l0, c0, c1, c2); // reverse line
        calculate(l0, l1, c2, c1, c0); // reverse curve
        calculate(l1, l0, c2, c1, c0); // reverse both
    }

    private static void calculate(Vec2 l0, Vec2 l1, Vec2 c0, Vec2 c1, Vec2 c2) {
        Line2 l = Line2.line(l0, l1);
        Bezier2.QuadraticBezier2 c =
                (Bezier2.QuadraticBezier2) Bezier2.curve(c0, c1, c2);

        System.out.println("\n# Intersections between");
        System.out.println("line: " + l);
        System.out.println("curve: " + c);

        Vec2[] intersections = Intersections.lineQuadratic(l, c);
        System.out.println(intersections.length + " intersections");
        for (Vec2 intersection : intersections) {
            System.out.println(intersection);
        }
    }
}

and the output of the program

# Intersections between
line: a=[x=512.0, y=96.0], b=[x=128.0, y=384.0]
curve: p0=[x=128.0, y=384.0], p1=[x=320.0, y=168.0], p2=[x=512.0, y=384.0]
0 intersections

# Intersections between
line: a=[x=128.0, y=384.0], b=[x=512.0, y=96.0]
curve: p0=[x=128.0, y=384.0], p1=[x=320.0, y=168.0], p2=[x=512.0, y=384.0]
2 intersections
[x=0.33333333333333326, y=0.3333333333333333]
[x=0.0, y=0.0]

# Intersections between
line: a=[x=512.0, y=96.0], b=[x=128.0, y=384.0]
curve: p0=[x=512.0, y=384.0], p1=[x=320.0, y=168.0], p2=[x=128.0, y=384.0]
2 intersections
[x=1.0, y=1.0]
[x=0.6666666666666666, y=0.6666666666666666]

# Intersections between
line: a=[x=128.0, y=384.0], b=[x=512.0, y=96.0]
curve: p0=[x=512.0, y=384.0], p1=[x=320.0, y=168.0], p2=[x=128.0, y=384.0]
2 intersections
[x=0.0, y=1.0]
[x=0.3333333333333333, y=0.6666666666666666]

[hamoid]

I keep digging :-)

Something goes wrong in solveQuadratic(). See how k is a very huge number (around the max double) which results in some variables becoming infinite.

And the reason seems to be that among a, b and c, there's two negative values and one is zero. The max value of the three is therefore zero.

getExponent() has two special cases, one of them is if the argument is zero or subnormal, then the result is Double.MIN_EXPONENT - 1.

Next I'll try comparing the failing case with the one that works in which one vertex is slightly displaced.

Update: the reason it works replacing (128.0, 384.0) with (128.01, 384.03) is because a, b and c are all non-zero, which doesn't help me yet.

[hamoid]

I don't know if this is a good fix, but it solves the issue for me:

diff --git a/src/io/lacuna/artifex/utils/Scalars.java b/src/io/lacuna/artifex/utils/Scalars.java
index 6341851..5aa7221 100644
--- a/src/io/lacuna/artifex/utils/Scalars.java
+++ b/src/io/lacuna/artifex/utils/Scalars.java
@@ -59,7 +59,11 @@ public class Scalars {
   }
 
   public static double normalizationFactor(double a, double b, double c) {
-    double exponent = getExponent(max(max(a, b), c));
+    double maxValue = max(a, b, c);
+    if(maxValue == 0.0) {
+      maxValue = min(a, b, c);
+    }
+    double exponent = getExponent(maxValue);
     return (exponent < -8 || exponent > 8) ? Math.pow(2, -exponent) : 1;
   }
 
@@ -75,4 +79,12 @@ public class Scalars {
   public static double max(double a, double b, double c) {
     return max(a, max(b, c));
   }
+
+  public static double min(double a, double b) {
+    return a > b ? b : a;
+  }
+
+  public static double min(double a, double b, double c) {
+    return min(a, min(b, c));
+  }
 }

[hamoid]

Is getExponent(max(abs(a), abs(b), abs(c))) better?

... Maybe not, as that might change the behavior by ruling out negative exponents...