curvature.py 文件源码

def interpolate_curvature(curvatures, times, kind='cubic'):
    """ Interpolate the curvatures.

    A curvature is a parametrisation `s`, d(angle)/ ds and a parameter between [0,1].
    Return a surface f(s,t).
    """
    import numpy as np

    from scipy import interpolate
    from scipy.ndimage import measurements

    if len(curvatures) == 3 and not isinstance(curvatures[2], tuple):
        curvatures= [curvatures]

    curv_abs = [c[0] for c in curvatures]
    curves = [c[1] for c in curvatures]
    params = [c[2] for c in curvatures]

    n = len(params)
    if n==1:
        curv_abs*=2
        curves*=2
        curves[0]=np.zeros(len(curves[0]))
        params = [0.,1.]
    elif False:
        if params[0] != 0.:
            params.insert(0,0.)
            curves.insert(0,curves[0])
            curv_abs.insert(0, curv_abs[0])
        if params[-1] != 0:
            params.append(1.)
            curves.append(curves[-1])
            curv_abs.append(curv_abs[-1])

    # compute a common parametrisation s
    # We conserve the last parametrisation because the result is very sensitive
    if True:
        min_s = min(np.diff(s).min() for s in curv_abs)
        s_new = np.unique(np.array(curv_abs).flatten())
        ds = np.diff(s_new)
        k = np.cumsum(ds >= min_s)
        labels = range(k.min(), k.max()+1)
        s = np.zeros(len(labels)+1)
        s[1:] = measurements.mean(s_new[1:], k, labels)
        s[-1]=1.
        s = s_new
    else:
        s = np.array(curv_abs[-1])

    # renormalise all the curves
    curves = [np.interp(s[1:-1], old_s[1:-1], old_crv) for old_s, old_crv in zip(curv_abs,curves)]

    # interpolate correctly the curvatures
    x = s[1:-1]
    y = np.array(params)
    z = np.array(curves)

    #f = interpolate.interp2d(y, x, z, kind=kind)
    f = interpolate.RectBivariateSpline(x, y, z.T, kx=1, ky=1)
    return s, f(x,times)