param_at_point gives wrong answer because it uses a local minimizer

[bdlucas1]

This issue concerns the code currently in the dev branch meant to fix #708 and #570, so this is not a duplicate of those issues.

The problem is that the minimizer used finds some local minimum, not necessarily a global minimum, so it can return the wrong result. A simple example is finding points on a circle:

curve = bd.Edge.make_circle(10)

n = 101
ts = np.linspace(0, 1, n)
points = [curve@t for t in ts]
computed_ts = [curve.param_at_point(p) for p in points]

for t, ct in zip(ts, computed_ts):
    print(f"t {t:.2f}, computed t {ct:.6f}, difference {abs(ct-t):.6f}")

Many cases are correct, but in many cases it has found a very wrong answer - the point on the other end of the circle - as a local minimum, for example:

t 0.68, computed t 0.680000, difference 0.000000
t 0.69, computed t 0.690000, difference 0.000000
t 0.70, computed t 0.700000, difference 0.000000
t 0.71, computed t 0.710000, difference 0.000000
t 0.72, computed t 0.720000, difference 0.000000
t 0.73, computed t 0.730000, difference 0.000000
t 0.74, computed t 0.000000, difference 0.740000
t 0.75, computed t 0.000000, difference 0.750000
t 0.76, computed t 0.000000, difference 0.760000
t 0.77, computed t 0.000000, difference 0.770000
t 0.78, computed t 0.000000, difference 0.780000
t 0.79, computed t 0.000000, difference 0.790000

This image shows the points at the computed parameter values. It should show 100 points evenly spaced around the circle, but some have been "moved" to an endpoint of the circle.

A more complex example:

curve = bd.Edge.make_spline([(-2,0,0), (-10,1,0), (0,0,0), (1,-2,0), (2,0,0), (1,1,0),])

Image shows points at the computed parameter values - should be 100 points evenly spaced along the curve. Some find a local minimum at the other end of the curve, but in some cases it finds a local minimum at a different point within the curve.

For your consideration here's an approach to finding a global minimum that seems to work pretty well. In my tests performance matches the performance of the current code in cases where the current code gives the right answer.

def param_at_point(curve, point, tol=1e-6, logmaxit=10):

    point = bd.Vector(point)

    # Function that we seek to find a zero of; this will also be a local minimum.
    f = lambda x: (curve.position_at(x[0]) - point).length

    # Sweep the parameter range from 0 to 1 repeatedly, on each sweep examining successively narrower
    # intervals, starting with the interval [0, 1], until we find an interval where a local minimum
    # of f is also a zero.
    #
    # For most simple curves we will succeed on the first iteration after examining the interval [0,1].
    # For somewhat more complex curves only a few iterations are required.
    #
    # logitmax=10 means we end after sweeping the range from 0 to 1 with intervals of length 1/1024.
    # This should be adequate for any reasonable curve, but could still fail on highly convoluted curves;
    # I guess it would take a spline with something on the order of 1024 knots
    
    n = 1
    for sweep in range(logmaxit):

        # Sweep the parameter range from 0 to 1 with n intervals of length 1/n
        for i in range(n):

            # Find a local minimum within the current interval [i/n, (i+1)/n].
            lo, hi = i/n, (i+1)/n
            p_lo, p_hi = curve.position_at(lo), curve.position_at(hi)
            d_lo, d_hi = (p_lo - point).length, (p_hi - point).length
            x0 = lo + d_lo / (d_hi + d_lo) * (hi - lo)
            result = scipy.optimize.minimize(f, x0=[x0], bounds=[(lo,hi)], method="SLSQP", tol=tol)

            # We've found a local minimum. If it is also a zero of f, to within
            # a tolerance on f(x) determined by the x tolerance we requested from the minimizer,
            # return the result; otherwise keep searching.
            if result.fun < (curve@(result.x[0]+tol) - curve@(result.x[0]-tol)).length:
                return result.x[0]

        # halve our search intervals for next sweep
        n *= 2

    # If the point is not on the curve we won't find a zero of f, so we get here.
    # We save a bit of computation by not checking this up front, making the common case faster,
    # at the expense of making it slow to determine that the point is not on the curve.
    # With logitmax=10 I measured about 1 sec to determine that a point is not on the curve.
    # This could be made faster by decreasing logitmax, so also decreasing complexity of curves
    # we can handle.
    raise Exception("point not on curve")