[Code challenge] Point at euclidian length

Let’s take back the challenges we were doing in the old GH forum. The one who offers the most efficient solution and meets all the criteria wins. The prize is to be the best no camper sniper :cowboy_hat_face:

The problem is given a curve and a position on that curve, find the point on that curve whose distance between the two positions is a given euclidean distance.

Conditions:

  1. Inputs: curve, curve parameter, desired length.
  2. Outputs: desired point, parameter at desired point.
  3. Efficiency: it should be computed 10,000 times in less than a second or two.
  4. Accuracy: the error must be below 1e-12.
  5. Generality. Supporting unsolvable cases (such as the length of the curve being less than the desired length) gives an advantage against other shooters with similar quality proposals.

Details:

  • For cases with more than one solution, you should opt for the one of the nearest parameter.
  • If the desired point is outside the domain of the curve, you can return an empty result.
  • It should work for any curve.
  • (Optional) If the input length is negative, the desired point must be in the previous section domain of the input position.

To have a place to start and understand the difficulty of the problem quickly, here are 3 proposals of mine that don’t really work. But I’ll keep drawing my weapons.
PointAtEuclidianLength-CodeChallenge.gh (15.8 KB)

If you have thought that intersecting a sphere with the curve will compute so fast, you will be the last one to shoot!

I’m waiting for your head shots.

Do you mean something like this? (I was experimenting with the “DivideCurveEquidistant” problem discussed recently but didn’t get to finish the full division part, only finding one Chord of the specified distance):

Fig. Approved by Montserrat :smiley: (turn down sound). The blue and yellow curves has equal distance:

// Rolf

1 Like

I guess, but does it meet the efficiency criteria (3) and the others?

Hm, the component does some other irrelevant stuff as well so… but a quick test gives 7.4 sec for
0.0001 and 13 sec for 0.0000000000001 (1e–13), so no, not under 2 sec, at least not without some tweaking… :thinking:

// Rolf

Not bad! Can you share the approach you’ve used, at least explain it?

The challenge is to make it efficient. The solution can be approached in many ways, or make it precise with intersections, but it is a challenge to make it fast and robust. And if a solution can be extended to surfaces, it will be assigned the veteran’s medal.

To be pedantic, I think the phrasing of the challenge definition needs to be refined a little - the point implies it always exists and is always unique.
Perhaps ‘find the first point in the increasing parameter space of the curve which is at a given Euclidean distance (if it exists)’.
Also, do the input curves have to be smooth? If they can be kinky I think the problem is harder.

1 Like

:+1:. The problem will be very different if it’s to find any/nearest/all point(s);

I updated the description a while ago to include in the requirement 5 that if the solution has several points it should be returned to the nearest parameter.

If the length is negative, the point in the domain section before the input parameter should be returned.

It should work for any curve.

I will update the description in a details section.

1 Like

I’d have to clean it up a bit before sharing. Too much other (irrelevant) stuff still in it.

I’ll see when I have the time, in just a minute I have a customer at the door buying some device I just sold.

I’ll be back later on this. My plan was to make it ready to a full-fledged DivideEquidistantCurve component. That will happen, but not today. :slight_smile:

// Rolf

1 Like

about 100 times faster.

image

for the 2nd case

image

Share if you want to shoot :cowboy_hat_face:

:wink: I’m still testing its robustness

1 Like

@Dani_Abalde, I think the test curve should have an interval in which the curve is “coming closer” to the start point, just to prove that the algorithm doesn’t get stuck when the length shrinks instead of expands during the algorithm…

Something like this; The Red lines shows where the chord (or whatever name) gets shorter as the end point advances but despite the negative trend the algorithm still advances along the curve… and eventually finds the blue length:


Edit: My test curve:
test_curve.gh (1.3 KB)

// Rolf

If no conditions are placed on the curve in the problem definition, a solution should even work on a curve like this:

2 Likes

Can you post the curve?

// Rolf

Please upload the curves with which to check the algorithms.

If your algorithm only finds one point, it’s better than nothing. If it finds several, it’s not difficult or slow to return the one with the closest to and highest than the input parameter.

trickycurve.3dm (10.7 KB)
(it even doubles back exactly overlapping itself in some places, but it is still a valid single joined curve)

I don’t actually think it always makes sense to try and make an algorithm which is guaranteed to work on absolutely all cases including even the most degenerate ones. Something that works in the vast majority of cases likely to occur in practical use is usually enough, ideally with some knowledge of what these limitations are.

I completely agree. Don’t be shy about upload proposals, the best, though imperfect, wins.

It seems my algorithm handles the trickycurve all the way. The “squary” part is tricky indeed… but I think it handles all the tricks? I stepped 1000 x 0.001 steps (t) with a tolerance of 0.0000000000001 (1e-13) and the algorithm iterated 8900 loops in total for the entire curve.

That’s only 8.9 iterations per step (t), but long intervals of the curve is quite simple for this specific length of the “chord”: (Edit: Hm, the clip doesn’t show… uploading it as a s zip then)
TrickyCurve_04.zip (2.8 MB)

// Rolf

1 Like

~50x faster. The algorithm is a binary search that automatically extracts monotone intervals from a curve.

gives the correct result for Daniel’s curve. The code doesn’t handle when the point doesn’t exist or the distance is negative, while they are easy to implement.

using System;
using System.Collections.Generic;
using System.Linq;
using Grasshopper.Kernel;
using Rhino;
using Rhino.Geometry;

namespace PancakeAlgo.Geometry
{
    internal static class GeometryExtensions
    {
        internal static IEnumerable<Interval> DivideIntoSubDomains(this Interval range, int division)
        {
            var segment = range.Length / division;
            for (var i = 0; i < division; i++)
            {
                yield return new Interval(range.T0 + segment * i, range.T0 + segment * (i + 1));
            }
        }

        internal static IEnumerable<double> GetDiscontinuities(this Curve crv, double start, Continuity type)
        {
            var max = crv.Domain.Max;
            yield return start;

            for (; ; )
            {
                if (!crv.GetNextDiscontinuity(type, start, max, out var t))
                {
                    if (Math.Abs(start - max) > RhinoMath.ZeroTolerance)
                        yield return start;
                    yield break;
                }

                yield return t;
                start = t;
            }
        }

        internal static IEnumerable<Interval> DivideDomainAtKnots(this Curve crv, double start)
        {
            var max = crv.Domain.Max;

            var kinks = GetDiscontinuities(crv, start, Continuity.Cinfinity_continuous)
                .Concat(GetDiscontinuities(crv, start, Continuity.C2_locus_continuous))
                .Distinct(DoubleEquality).ToArray();
            Array.Sort(kinks);

            for (var i = 0; i < kinks.Length - 1; i++)
                yield return new Interval(kinks[i], kinks[i + 1]);
        }

        private static readonly DoubleEqualityWithError DoubleEquality = new DoubleEqualityWithError();
        private class DoubleEqualityWithError : IEqualityComparer<double>
        {
            public bool Equals(double x, double y) => Math.Abs(x - y) < RhinoMath.ZeroTolerance;

            public int GetHashCode(double obj) => obj.GetHashCode();
        }
    }
    public class pcgPointLengthAway : GH_Component
    {
        /// <summary>
        /// Initializes a new instance of the pcgPointLengthAway class.
        /// </summary>
        public pcgPointLengthAway()
          : base("pcgPointLengthAway", "pcgPointLenAway",
              "Description",
              "Pancake", "Geometry")
        {
        }

        protected override void RegisterInputParams(GH_Component.GH_InputParamManager pManager)
        {
            pManager.AddCurveParameter("Curve", "C", "", GH_ParamAccess.item);
            pManager.AddNumberParameter("Parameter", "t", "", GH_ParamAccess.item);
            pManager.AddNumberParameter("Distance", "L", "", GH_ParamAccess.list);
            pManager.AddNumberParameter("Tolerance", "T", "", GH_ParamAccess.item, 1e-12);
            pManager.AddIntegerParameter("Division", "D", "", GH_ParamAccess.item, 1);
            pManager.AddNumberParameter("Bias", "B", "", GH_ParamAccess.item, 0.5);
        }

        protected override void RegisterOutputParams(GH_Component.GH_OutputParamManager pManager)
        {
            pManager.AddNumberParameter("Parameter", "t", "", GH_ParamAccess.list);
        }

        protected override void SolveInstance(IGH_DataAccess DA)
        {
            Curve crv = null;
            double t = 0;
            var listOfD = new List<double>();

            var tolerance = 1e-12;
            var div = 1;

            DA.GetData(0, ref crv);
            DA.GetData(1, ref t);
            DA.GetDataList(2, listOfD);
            DA.GetData(3, ref tolerance);
            DA.GetData(4, ref div);
            DA.GetData(5, ref BinaryBias);

            var origin = crv.PointAt(t);
            var domains = crv.DivideDomainAtKnots(t)
                .Select(it => it.DivideIntoSubDomains(div)).SelectMany(it => it).ToArray();
            var paramTolerance = EstimateParameterTolerance(crv, tolerance);

            var result = listOfD.AsParallel().AsOrdered().Select(d =>
            {
                var squareD = d * d;
                for (var i = 0; i < domains.Length; i++)
                {
                    if (LookForOnSegment(crv, domains[i].T0, domains[i].T1, origin, squareD, paramTolerance, tolerance, out var p))
                    {
                        return p;
                    }
                }
                return double.PositiveInfinity;
            });

            DA.SetDataList(0, result);
        }

        private static double EstimateParameterTolerance(Curve crv, double distTolerance)
        {
            var length = crv.GetLength();
            var domainLength = crv.Domain.Length;

            return distTolerance / length * domainLength * 1e-2;
        }

        private static double BinaryBias = 0.5;
        private static bool LookForOnSegment(Curve crv, double start, double end, Point3d origin,
            double squareDist, double paramTolerance, double distTolerance, out double param)
        {
            param = 0;
            var lowerBound = start;
            var upperBound = end;

            var ptrAtStart = crv.PointAt(lowerBound);
            if (Math.Abs((ptrAtStart - origin).SquareLength - squareDist) <= distTolerance)
            {
                param = lowerBound;
                return true;
            }

            var ptrAtEnd = crv.PointAt(upperBound);
            if (Math.Abs((ptrAtEnd - origin).SquareLength - squareDist) <= distTolerance)
            {
                param = upperBound;
                return true;
            }

            var iteration = 0;
            var bias = BinaryBias;

            for (; ; ++iteration)
            {
                var middleGround = (upperBound - lowerBound) * bias + lowerBound;

                var p = crv.PointAt(middleGround);

                var delta = (p - origin).SquareLength - squareDist;

                if (Math.Abs(delta) <= distTolerance)
                {
                    param = middleGround;
                    return true;
                }

                if (upperBound - lowerBound < paramTolerance)
                    return false;

                if (delta > 0)
                {
                    upperBound = middleGround;
                    continue;
                }
                else if (delta < 0)
                {
                    lowerBound = middleGround;
                    continue;
                }

                throw new ArithmeticException();
            }
        }

        public override Guid ComponentGuid => new Guid("d1b7aae7-7280-417f-9bf6-0845f5394eac");
    }
}
6 Likes