Higher order local interpolation of implicit curves in Python

Question

Given a set of points describing some trajectory in the 2D plane, I would like to provide a smooth representation of this trajectory with local high order interpolation.

For instance, say we define a circle in 2D with 11 points in the figure below. I would like to add points in between each consecutive pair of points in order or produce a smooth trace. Adding points on every segment is easy enough, but it produces slope discontinuities typical for a "local linear interpolation". Of course it is not an interpolation in the classical sense, because

• the function can have multiple y values for a given x
• simply adding more points on the trajectory would be fine (no continuous representation is needed).

so I'm not sure what would be the proper vocabulary for this.

The code to produce this figure can be found below. The linear interpolation is performed with the lin_refine_implicit function. I'm looking for a higher order solution to produce a smooth trace and I was wondering if there is a way of achieving it with classical functions in Scipy? I have tried to use various 1D interpolations from scipy.interpolate without much success (again because of multiple y values for a given x).

The end goals is to use this method to provide a smooth GPS trajectory from discrete measurements, so I would think this should have a classical solution somewhere.

import numpy as np
import matplotlib.pyplot as plt

def lin_refine_implicit(x, n):
    """
    Given a 2D ndarray (npt, m) of npt coordinates in m dimension, insert 2**(n-1) additional points on each trajectory segment
    Returns an (npt*2**(n-1), m) ndarray
    """
    if n > 1:
        m = 0.5*(x[:-1] + x[1:])
        if x.ndim == 2:
            msize = (x.shape[0] + m.shape[0], x.shape[1])
        else:
            raise NotImplementedError

        x_new = np.empty(msize, dtype=x.dtype)
        x_new[0::2] = x
        x_new[1::2] = m
        return lin_refine_implicit(x_new, n-1)
    elif n == 1:
        return x
    else:
        raise ValueError
n = 11
r = np.arange(0, 2*np.pi, 2*np.pi/n)
x = 0.9*np.cos(r)
y = 0.9*np.sin(r)
xy = np.vstack((x, y)).T
xy_highres_lin = lin_refine_implicit(xy, n=3)

plt.plot(xy[:,0], xy[:,1], 'ob', ms=15.0, label='original data')
plt.plot(xy_highres_lin[:,0], xy_highres_lin[:,1], 'dr', ms=10.0, label='linear local interpolation')
plt.legend(loc='best')
plt.plot(x, y, '--k')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('GPS trajectory')
plt.show()

Kristianson

This is called parametric interpolation. scipy.interpolate.splprep provides spline approximations for such curves (provided you know the order in which the points are; if you don't know which point comes after which on the curve, the problem becomes more difficult). – Unrelenting

Answer 1

This is called parametric interpolation.

scipy.interpolate.splprep provides spline approximations for such curves. This assumes you know the order in which the points are on the curve.

If you don't know which point comes after which on the curve, the problem becomes more difficult. I think in this case, the problem is called manifold learning, and some of the algorithms in scikit-learn may be helpful in that.

Unrelenting

Answer 2

I would suggest you try to transform your cartesian coordinates into polar coordinates, that should allow you to use the standard scipy.interpolation without issues as you won't have the ambiguity of the x->y mapping anymore.

Coates

That is a valid point. Still it has a few limitations in the general case a) the choice of the origin for the polar coordinates could have an impact on the result, and it is therefore a bit ambiguous b) in the case when the trajectory is doing loops around the origin, we would have to handle the unwrapping of the polar angle,which is getting more complex c) is it specific for 2D, say lin_refine_implicit works with points in N dimensions (2D, 3D etc), while one would have to rewrite the whole code (with lots of trigonometry) to go from 2D to 3D. – Kristianson

More importantly, if possible I'd like to avoid implementing a solution myself and use some well tested and documented algorithm, if possible. I would think that one should exist, maybe in some GIS library, if not in Scipy. Thank you for your answer though, I'll follow it if I can't find any alternatives. – Kristianson