~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");
}
}
(turn down sound). The blue and yellow curves has equal distance:
. The problem will be very different if it’s to find any/nearest/all point(s);
I’m still testing its robustness
