The function cannot detect intersections of two curves

Question

The below source code is supposed to find an intersection between two curves.

However, the function is unable to detect the intersection point.

How can I fix it?

using MathNet.Numerics.Interpolation;
using System.Collections.Generic;
using System.Linq;

public static Vec2 FindIntersectionOfTwoCurves(List<double> xList1, List<double> yList1,
                                               List<double> xList2, List<double> yList2)
{
    if (xList1 == null || yList1 == null || xList2 == null || yList2 == null ||
        xList1.Count != yList1.Count || xList2.Count != yList2.Count)
        return null;

    IInterpolation interpolation1 = Interpolate.Linear(xList1, yList1);
    IInterpolation interpolation2 = Interpolate.Linear(xList2, yList2);

    double lowerBound = Math.Max(xList1.Min(), xList2.Min());
    double upperBound = Math.Min(xList1.Max(), xList2.Max());

    double step = (upperBound - lowerBound) / 1000; // adjust the step size as needed
    for (double x = lowerBound; x <= upperBound; x += step)
    {
        double y1 = interpolation1.Interpolate(x);
        double y2 = interpolation2.Interpolate(x);

        if (Math.Abs(y1 - y2) < 1e-7)
        {
            return new Vec2(x, y1);
        }
    }

    return null;
}

public class Vec2
{
    public double X { get; set; }
    public double Y { get; set; }

    public Vec2(double x, double y)
    {
        X = x;
        Y = y;
    }
}

Answer 1

since you are trying to find an intersection of two curves with a non-continous first derivative you should use the Bisection method.

a quickly written C# implementation is as follows, the implementation of binary search is from here, the binary search can be replaced by a simpler linear interpolation if X is linearly spaced.

namespace ConsoleApp2
{
    public class LinearInterpolator
    {
        private static int BinarySearch<T>(Span<T> list, T value)
        {
            if (list == null)
                throw new ArgumentNullException("list");
            var comp = Comparer<T>.Default;
            int lo = 0, hi = list.Length - 1;
            while (lo < hi)
            {
                int m = (hi + lo) / 2;  // this might overflow; be careful.
                if (comp.Compare(list[m], value) < 0) lo = m + 1;
                else hi = m - 1;
            }
            if (comp.Compare(list[lo], value) < 0) lo++;
            return lo;
        }

        private static int FindFirstIndexGreaterThanOrEqualTo<T>
                                  (T[] sortedList, T key)
        {
            return BinarySearch(sortedList, key);
        }

        double[] x_values;
        double[] y_values;
        public LinearInterpolator(List<double> x, List<double> y)
        {
            // quick argsort
            int[] indicies = x.AsEnumerable().Select((v, i) => new { obj = v, index = i }).OrderBy(c => c.obj).Select(c => c.index).ToArray();
            x_values = indicies.Select(i => x[i]).ToArray();
            y_values = indicies.Select(i => y[i]).ToArray();
        }

        public double Interpolate(double x)
        {
            int index = FindFirstIndexGreaterThanOrEqualTo(x_values, x);
            if (index == 0)
            {
                return y_values[0];
            }
            double y1 = y_values[index - 1];
            double y2 = y_values[index];
            double x1 = x_values[index - 1];
            double x2 = x_values[index];
            return (x - x1) / (x2 - x1) * (y2 - y1) + y1;
        }

    }

    class IntersectionFinder
    {
        public static Nullable<(double, double)> FindIntersectionOfTwoCurves(List<double> xList1, List<double> yList1,
                                                   List<double> xList2, List<double> yList2)
        {
            if (xList1 == null || yList1 == null || xList2 == null || yList2 == null ||
                xList1.Count != yList1.Count || xList2.Count != yList2.Count)
                return null;

            LinearInterpolator interpolator1 = new LinearInterpolator(xList1, yList1);
            LinearInterpolator interpolator2 = new LinearInterpolator(xList2, yList2);

            double lowerBound = Math.Max(xList1.Min(), xList2.Min());
            double upperBound = Math.Min(xList1.Max(), xList2.Max());

            if (lowerBound > upperBound) // X axes don't overlap
            {
                return null;
            }

            double diff_start = interpolator1.Interpolate(lowerBound) - interpolator2.Interpolate(lowerBound);
            double diff_end = interpolator1.Interpolate(upperBound) - interpolator2.Interpolate(upperBound);

            if ((diff_start > 0 && diff_end > 0) || (diff_start < 0 && diff_end < 0)) // intersection doesn't exist
            {
                return null;
            }

            double mid = (lowerBound + upperBound) / 2;
            double diff_mid = interpolator1.Interpolate(mid) - interpolator2.Interpolate(mid);

            int iterations = 0;
            while (Math.Abs(diff_mid) > 1e-7)
            {
                if (diff_start > diff_end) // list1 is higher
                {
                    if (diff_mid > 0) // mid is also higher, intersection in right side
                    {
                        lowerBound = mid;
                        diff_start = diff_mid;
                    }
                    else // mid is lower, intersection in left side
                    {
                        upperBound = mid;
                        diff_end = diff_mid;
                    }
                }
                else // list 2 is higher
                {
                    if (diff_mid < 0)
                    {
                        lowerBound = mid;
                        diff_start = diff_mid;
                    }
                    else
                    {
                        upperBound = mid;
                        diff_end = diff_mid;
                    }
                }
                mid = (lowerBound + upperBound) / 2;
                diff_mid = interpolator1.Interpolate(mid) - interpolator2.Interpolate(mid);
                iterations++;
                if (iterations > 10000) // prevent infinite loop if Y is discontinuous
                {
                    return null;
                }
            }

            return new(mid, interpolator1.Interpolate(mid));
        }
    }

    internal class Program
    {
        static void Main(string[] args)
        {
            List<double> xList1 = [1, 1.5, 2, 3.5, 4];
            List<double> xList2 = [1, 2, 3, 4, 5];
            List<double> yList1 = [0, 2, 4, 6, 8];
            List<double> yList2 = [8, 6, 4, 2, 0];
            Console.WriteLine(IntersectionFinder.FindIntersectionOfTwoCurves(xList1, yList1, xList2, yList2).ToString());
        }
    }
}

(2.5999999940395355, 4.799999992052714)

Answer 2

This small change can fix your problem. Don't watch for Math.Abs(y1 - y2) to be less than some epsilon (e.g. 1e-7). Instead, watch for the sign of y1 - y2 to change so that it's different from the previous iteration. That is called looking for a zero-crossing.

For example:

int? prevSign = null;
for (double x = lowerBound; x <= upperBound; x += step) {
    double y1 = interpolation1.Interpolate(x);
    double y2 = interpolation2.Interpolate(x);
    int sign = Math.Sign(y1 - y2);
    if (prevSign.HasValue) {
        if (sign != prevSign.Value) {
            return new Vec2(x, y1);
        }
    } else {
        prevSign = sign;
    }
}