scipy - Distance between a point and a curve in python

Question

I have several points and a curve discribed as two lists including positions. I try to get the list of differences between the points and the curves. I tried to follow this web, but I do not understand this command:

X = fmin_cobyla(objective, x0=[0.5,0.5], cons=[c1])

What are the correct arguments in my case, please?

import numpy as np  
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla

data = np.loadtxt('O_Cout.dat', unpack=True, usecols=[0, 2])

z_1v1 = np.polyfit(data[0], data[1], 2)
f_1v1 = np.poly1d(z_1v1)
# Creating more points on the streamline - defining new time with more steps
x_new = list(np.arange(0,100000,1))
y_new = f_1v1(x_new)

# Plot figure with size
fig, ax = plt.subplots()

ax.scatter(data[0], data[1])

ax.plot(x_new, y_new)


def objective(X):
    x,y = X
    return np.sqrt((x - P[0])**2 + (y - P[1])**2)

def c1(X):
    x,y = X
    return f(x) - y


for i in range(len(data[1])-1):
    P = (data[0][i], data[1][i])
    print(P)

    X = fmin_cobyla(objective, x0=[0.5,0.5], cons=[c1])

    print ('The minimum distance is', objective(X))


# Save the figure
plt.tight_layout()
plt.savefig('OC_parabola.png')

Answer 1

The script you found is meant for a known function f(x) but IIUC you don't know f(x): your curve is only defined by the coordinates (x,y) and you don't know f(x) so that y=f(x).

In this case, you can use the same basics.

Given a point P

and a curve defined by the coordinates (x,y), the distance between the point P and a point of the curve can be simply defined by

that we wish to minimize, i.e. find the minimum/a in the defined domain.

For example

import numpy as np  
import matplotlib.pyplot as plt

# Here I define a function f(x) to
# generate y coordinates, but let us
# suppose we don't know it and that
# we got only x and y
def f(x):
    return np.cbrt( np.exp(2*x) -1 )

# This is what we really got
x = np.linspace(-2, 2, 1000)
y = f(x)

# The point P
P = (.5, .5)

fig, ax = plt.subplots(figsize=(7, 7))
ax.plot(x, y, lw=4)
ax.plot(*P, 'or')
ax.text(
    P[0], P[1], 
    f"  P ({P[0]}, {P[1]})", 
    ha='left', va='center',
    fontsize=15
)
ax.set(
    xlim=(-2, 2),
    ylim=(-2, 2),
)
plt.show()

Let us define the function d, the distance between a point P and the curve

def distance(x, y, x0, y0):
    d_x = x - x0
    d_y = y - y0
    dis = np.sqrt( d_x**2 + d_y**2 )
    return dis

and now compute d between the given P and (x,y) and find the minimum

from scipy.signal import argrelmin

# compute distance
dis = distance(x, y, P[0], P[1])
# find the minima
min_idxs = argrelmin(dis)[0]
# take the minimum
glob_min_idx = min_idxs[np.argmin(dis[min_idxs])]
# coordinates and distance
min_x = x[glob_min_idx]
min_y = y[glob_min_idx]
min_d = dis[glob_min_idx]

and plot results

fig, ax = plt.subplots(figsize=(7, 7))

ax.plot(x, y, lw=4)
ax.plot(
    [P[0], min_x],
    [P[1], min_y],
    'k--', lw=1,
    label=f'distance {min_d:.2f}'
)
ax.plot(*P, 'or')
ax.text(
    P[0], P[1], 
    f"  P ({P[0]}, {P[1]})", 
    ha='left', va='center',
    fontsize=15
)
ax.set(
    xlim=(-2, 2),
    ylim=(-2, 2),
)
ax.legend()
plt.show()

EDIT

Improving, it can be defined a simple function to return all minimum distances, for example

import numpy as np  
import matplotlib.pyplot as plt

def distance(x, y, x0, y0):
    """
    Return distance between point
    P[x0,y0] and a curve (x,y)
    """
    d_x = x - x0
    d_y = y - y0
    dis = np.sqrt( d_x**2 + d_y**2 )
    return dis

def min_distance(x, y, P, precision=5):
    """
    Compute minimum/a distance/s between
    a point P[x0,y0] and a curve (x,y)
    rounded at `precision`.
    
    ARGS:
        x, y      (array)
        P         (tuple)
        precision (int)
        
    Returns min indexes and distances array.
    """
    # compute distance
    d = distance(x, y, P[0], P[1])
    d = np.round(d, precision)
    # find the minima
    glob_min_idxs = np.argwhere(d==np.min(d)).ravel()
    return glob_min_idxs, d

that works even if there is more than one minimum

def f(x):
    return x**2

x = np.linspace(-2, 2, 1000)
y = f(x)

P = (0, 1)

min_idxs, dis = min_distance(x, y, P)

fig, ax = plt.subplots(figsize=(7, 7))

ax.plot(x, y, lw=4)
for idx in min_idxs:
    ax.plot(
        [P[0], x[idx]],
        [P[1], y[idx]],
        '--', lw=1,
        label=f'distance {dis[idx]:.2f}'
    )
ax.plot(*P, 'or')
ax.text(
    P[0], P[1], 
    f"  P ({P[0]}, {P[1]})", 
    ha='left', va='center',
    fontsize=15
)
ax.set(
    xlim=(-2, 2),
    ylim=(-1, 3),
)
ax.legend()
plt.show()

def f(x):
    return np.sqrt(4 - x**2)

x = np.linspace(-2, 2, 21)
y = f(x)

P = (0, 0)

min_idxs, dis = min_distance(x, y, P)

fig, ax = plt.subplots(figsize=(7, 7))

ax.plot(x, y, lw=4)
for idx in min_idxs:
    ax.plot(
        [P[0], x[idx]],
        [P[1], y[idx]],
        '--', lw=1,
        label=f'distance {dis[idx]:.2f}'
    )
ax.plot(*P, 'or')
ax.text(
    P[0], P[1], 
    f"  P ({P[0]}, {P[1]})", 
    ha='left', va='center',
    fontsize=15
)
ax.set(
    xlim=(-2, 2),
    ylim=(-1, 3),
)
ax.legend(loc='upper left', bbox_to_anchor=(1,1))
plt.show()