CODEWITHSUNDEEP
pythoncurve
newstudent
newstudent

Reputation: 444

Curvature of set of points

data file 

import matplotlib.pylab as plt
import numpy as np

#initial data
data=np.loadtxt('profile_nonoisebigd02.txt')
x=data[:,0]
y=data[:,1]

initial profile 

#first derivatives 
dx= np.gradient(data[:,0])
dy = np.gradient(data[:,1])

#second derivatives 
d2x = np.gradient(dx)
d2y = np.gradient(dy)

#calculation of curvature from the typical formula
curvature = np.abs(dx * d2y - d2x * dy) / (dx * dx + dy * dy)**1.5

curvature

Can anyone help me out as to where am I going wrong with the curvature ? The set of points give me a parabola, but curvature is not what I expect.

Upvotes: 1

Views: 7278

Answers (1)

hilberts_drinking_problem
hilberts_drinking_problem

Reputation: 11602

It seems that your data just isn't smooth enough; I used pandas to replace x, y, dx, dy, d2x, d2y and curvature by rolling means for different values window sizes. As the window size increases, the curvature starts to look more and more like what you'd expect to see for a smooth parabola (legend gives window size):

enter image description here

For the reference, here is the plot of your original data:

enter image description here

The code used to create the smoothed frames:

def get_smooth(smoothing=10, return_df=False):
    data=np.loadtxt('profile_nonoisebigd02.txt')

    if return_df:
        return pd.DataFrame(data)

    df = pd.DataFrame(data).sort_values(by=0).reset_index(drop=True).rolling(smoothing).mean().dropna()

    # first derivatives
    df['dx'] = np.gradient(df[0])
    df['dy'] = np.gradient(df[1])

    df['dx'] = df.dx.rolling(smoothing, center=True).mean()
    df['dy'] = df.dy.rolling(smoothing, center=True).mean()

    # second derivatives
    df['d2x'] = np.gradient(df.dx)
    df['d2y'] = np.gradient(df.dy)

    df['d2x'] = df.d2x.rolling(smoothing, center=True).mean()
    df['d2y'] = df.d2y.rolling(smoothing, center=True).mean()


    # calculation of curvature from the typical formula
    df['curvature'] = df.eval('abs(dx * d2y - d2x * dy) / (dx * dx + dy * dy) ** 1.5')
    # mask = curvature < 100

    df['curvature'] = df.curvature.rolling(smoothing, center=True).mean()

    df.dropna(inplace=True)
    return df[0], df.curvature

Upvotes: 5

Related Questions