curvature calculation incorrect at path edges

[nberliner]

The current implementation of curvature calculation uses the local curvature expression. Calculation of the derivatives is done exploiting the fact that (f*g)'= f'*g = f*g' with * the convolution (see wikipedia and stackoverflow). This however implicates that we are not calculating the derivative of our function but of the function convolved with a gaussian.

In general this might work fine for many data points, whereas smoothing might become problematic for low sampling or few datapoints with significant edge contributions. As an example I generate points on a half circle with radius 5 and plotted the calculated curvature at each point. Towards the edges of the path the radii values are deviating significantly from the true value!

This plot can be generated using the following code:

import numpy as np
import matplotlib.pylab as plt
from scipy.ndimage import gaussian_filter1d

class Color:
    """
    Helper to assign colors to float or integer values mapped to a given range.
    """
    def __init__(self, scaleMin=None, scaleMax=None):
        self.Nglobal = dict()
        self.cmap = plt.get_cmap('seismic_r')

        self.scaleMin = scaleMin
        self.scaleMax = scaleMax

        if scaleMin == scaleMax and scaleMin is not None:
            print('Warning: trying to set zero scaling range!')
            self.scaleMin = scaleMin
            self.scaleMax = scaleMin * 1.1

    def __call__(self, N):

        if self.scaleMin is None and self.scaleMax is not None:
            c = float(N) / self.scaleMax
        elif self.scaleMin is not None and self.scaleMax is None:
            c = float(N)  - self.scaleMin
        elif self.scaleMin is None and self.scaleMax is None:
            c = float(N)
        else:
            c = (float(N) - self.scaleMin) / ( self.scaleMax - self.scaleMin )

        return self.cmap(c)

def curvature(x, y, sigma=1, absolute=False):
    # Credit goes here: http://stackoverflow.com/a/28310758/1922650
    #first and second derivative
    x1 = gaussian_filter1d(x,  sigma=sigma, order=1, mode='wrap')
    x2 = gaussian_filter1d(x1, sigma=sigma, order=1, mode='wrap')
    y1 = gaussian_filter1d(y,  sigma=sigma, order=1, mode='wrap')
    y2 = gaussian_filter1d(y1, sigma=sigma, order=1, mode='wrap')
    if absolute:
        return np.abs(x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)
    else:
        return (x1*y2 - y1*x2) / np.power(x1**2 + y1**2, 3./2)

# The user defined parameters
radius   = 5.0
sampling = 30
sigma    = 1

# Generate the test data
alpha = np.linspace(-np.pi/2,np.pi/2, sampling)
X = radius*np.cos(alpha)
Y = radius*np.sin(alpha)

# Calculate the radius from the curvature
C = 1/curvature(X, Y, sigma=sigma)
curvatureMax   = np.max(C)
curvatureMin   = np.min(C)

color  = Color(scaleMin=curvatureMin, scaleMax=curvatureMax)
cColor = [ color(i) for i in C ]

# Plot the test data color coded for curvature
fig = plt.figure(figsize=(14,7), dpi=120)
ax1  = fig.add_subplot(121)
ax1.set_title("Generated data")
ax1.set_xlim( [-radius-1,radius+1] )
ax1.set_ylim( [-radius-1,radius+1] )
ax1.scatter(x=X, y=Y, c=cColor, alpha=1, edgecolor='none', s=10, cmap=plt.cm.seismic_r)

ax2 = fig.add_subplot(122)
ax2.set_title("Calculated curvature")
ax2.plot(C)
ax2.axhline(y=radius, color='red')

[nberliner]

One way around this issue would be to directly calculate the curvature from the triangle that is spanned by the points along the path (see here for some math refresher). This would look something like this: and should give the exact values for each point. I am unsure however if we will need to add some smoothing after curvature calculation (some form is already implemented and should be easy to adapt).

[nberliner]

I looked into it a bit more and with the proposed method above the edge still is not well defined. In addition, and more problematic, we would loose the information about the sign of the curvature.

Hence, instead of using the above method I switched to np.gradient() to calculate the derivatives (previously this was done using gaussian_filter1d() and continue to use the local curvature expression . The edges will still be problematic but much less so compared to the gaussian_filter1d() method. See below for a test done on simulated data on a half-circle with radius 5.

[nberliner]

I fixed another problem with displaying the curvature values and the curvature looks now fine to me. I will close this issue for now.

Note: so far the curvature color coding was recalculated for every contour segment, i.e. the maximum and minimum values defining the color spectrum were reset every time. This makes the plots not comparable between each other since the displayed colors do not correspond to the same curvature values. This is now fixed and the same color coding is used for all plots.