2
$\begingroup$

I am using a parametrized natural spline to interpolate 2D curves. For example something like this: Airfoil

Where I parametrize both the x and the y coordinate with a parameter that increases by one between each of the given points. The interpolation looks very good, here the x coordinate as a function of the parameter t: enter image description here

The first derivative with respect to t already shows a few wiggles: enter image description here

And the second derivative is oscillating heavily: enter image description here

One interesting fact: The oscillations appear exactly where the spacing of the points changes from uniform to a cosine spacing. When I apply the spline to a circle or ellipse I obtain perfectly smooth first and second derivatives.

Anybody has an idea why the derivatives start to oscillate and how to prevent this?

Improve this question
$\endgroup$
9
  • $\begingroup$ Here the other two images: [2]: i.stack.imgur.com/5vTsj.png [3]: i.stack.imgur.com/TWSVs.png $\endgroup$
    – benno
    May 9, 2017 at 8:47
  • 1
    $\begingroup$ Where are located the kinks in the derivatives in your upper image? $\endgroup$
    – nicoguaro ♦
    May 9, 2017 at 23:04
  • $\begingroup$ The spline starts at the right (sharp) end, goes around the profile counter clockwise and ends at the right end again. So the kinks are located somewhere in between the two ends. $\endgroup$
    – benno
    May 10, 2017 at 5:19
  • $\begingroup$ Can you write down the equations for your geometry? What happens if you increase the number of points? $\endgroup$
    – nicoguaro ♦
    May 10, 2017 at 5:23
  • 1
    $\begingroup$ What is the actual parametrization you are using? Can you write it down? Does it have continuous derivatives? $\endgroup$
    – Kirill
    May 10, 2017 at 20:14

2 Answers 2

Reset to default
6
$\begingroup$

Your choice of parameterization is creating problems. Instead of spanning one in $t$ between points, span an amount proportional to the line segment between the two points in $(x,y)$ space.

I've created a Python example that demonstrates the issue. It compares a uniform parameterization $(x_t(t),y_t(t))$ (like yours, but mine is scaled so that $t \in [0,1]$) with a length-based parameterization $(x_u(u),y_u(u))$.

Here is my example airfoil:

example airfoil x,y points

Here are the $x$ and $y$ coordinates as functions of parameters $t$ and $u$.

parameterization w.r.t t and u

These are the second derivatives. Note that the uniform parameterization shows these oscillations in the second derivative, just like you've shown, but the length-based parameterization looks much smoother.

second derivatives w.r.t. t and u

Finally, I calculated the curvature using the formula

$\kappa = |x' y'' - y' x''| / (x'^2+y'^2)^{3/2}$

from here, where the derivatives are taken w.r.t. the parameter $t$ or $u$. Note that in this curve I've calculated $\kappa_t(t)$ using the formula above with derivatives w.r.t the parameter $t$, but I've plotted against $u$ so that we can see the comparison between the two calculated curvatures,which should be similar. As we can see, the uniform parameterization introduces some artifacts, but the length-based parameterization is pretty smooth.

curvature

Finally, here's the Python code I used to generate these plots

from numpy import *
from matplotlib.pyplot import *
from scipy import interpolate

# first create a list of x,y points based on an airfoil design from
# https://en.wikipedia.org/wiki/NACA_airfoil
# created so that points are not uniformly distributed
c=1.0
th=0.15
x=concatenate( (linspace(1.0,0.1,40),linspace(0.95,0.05,20)**3.0 *0.1 ) )
y=5.0*th*(0.2969*np.sqrt(x/c) - 0.1260*(x/c) - 0.3516*(x/c)**2 +0.2843*(x/c)**3 - 0.1015*(x/c)**4)

x=concatenate( (x,flipud(x)))
y=concatenate( (y,-flipud(y)))


# create parameters. t is uniformly spaced on the interval 0 to 1. 
# u is spaced proportional to edge lengths
t=arange( len(x),dtype=float) / (len(x)-1)
u=zeros( (len(x)),dtype=float)
L=0.0
u[0]=0.0
for ii in range(1,len(x)):
  L += sqrt( (x[ii]-x[ii-1])**2 + (y[ii]-y[ii-1])**2 )
  u[ii]=L
u /= L

# create two cubic splines for each x and y
csxt=interpolate.splrep(t,x,k=3)
csxu=interpolate.splrep(u,x,k=3)
csyt=interpolate.splrep(t,y,k=3)
csyu=interpolate.splrep(u,y,k=3)

# tt and uu are fine samples used for plotting
tt=linspace(0.0,1.0,10000)
uu=linspace(0.0,1.0,10000)


# evaluate x(t),y(t),x(u),y(u) and their derivatives
xt0=interpolate.splev(tt,csxt)
xt1=interpolate.splev(tt,csxt,der=1)
xt2=interpolate.splev(tt,csxt,der=2)
yt0=interpolate.splev(tt,csyt)
yt1=interpolate.splev(tt,csyt,der=1)
yt2=interpolate.splev(tt,csyt,der=2)

xu0=interpolate.splev(uu,csxu)
xu1=interpolate.splev(uu,csxu,der=1)
xu2=interpolate.splev(uu,csxu,der=2)
yu0=interpolate.splev(uu,csyu)
yu1=interpolate.splev(uu,csyu,der=1)
yu2=interpolate.splev(uu,csyu,der=2)

# calculate curvature by the formula
# |x' y'' - y' x'' | / |x'^2 + y'^2|^3/2
Kt=abs( xt1*yt2 - yt1*xt2) / sqrt(xt1*xt1 + yt1*yt1)**3
Ku=abs( xu1*yu2 - yu1*xu2) / sqrt(xu1*xu1 + yu1*yu1)**3


# interpolate between t and u, so that we can plot
# apples to apples
cstu=interpolate.splrep(t,u,k=1)
uu2=interpolate.splev(tt,cstu)

# plots
ff=figure(1)
ff.clf()
plot(x,y,'.',markersize=6)
gca().set_aspect('equal')
grid('on')
xlabel('x')
ylabel('y')
axis([-0.05,1.05,-0.2,0.2])

ff=figure(2)
ff.clf()
subplot(211)
plot(tt,xt0,'b-',label=r"$x_t(t)$")
plot(tt,yt0,'r-',label=r"$y_t(t)$")
grid('on')
xlabel('t')
legend(loc='upper right',fontsize=20)
subplot(212)
plot(uu,xu0,'b-',label=r"$x_u(u)$")
plot(uu,yu0,'r-',label=r"$y_u(u)$")
grid('on')
xlabel('u')
legend(loc='upper right',fontsize=20)



ff=figure(3)
ff.clf()
subplot(211)
plot(tt,xt2,'b-',label=r"$x_t''(t)$")
plot(tt,yt2,'r-',label=r"$y_t''(t)$")
grid('on')
xlabel('t')
legend(loc='upper right',fontsize=20)
subplot(212)
plot(uu,xu2,'b-',label=r"$x_u''(u)$")
plot(uu,yu2,'r-',label=r"$y_u''(u)$")
grid('on')
xlabel('u')
legend(loc='upper right',fontsize=20)

ff=figure(4)
ff.clf()
plot(uu2,Kt,'b-',label=r"$\kappa_t(u)$")
plot(uu,Ku,'r-',label=r"$\kappa_u(u)$")
grid('on')
xlabel('u')
gca().set_xlim([0.4,0.6])
legend(loc='upper right',fontsize=20)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much. I did not expect such a detailed answer :) I also have an arc length parametrization in my program. But i first create a spline wehre the parameter always increases by one. Then i integrate the arc length ($s = \int_{0}^{t_{max}}\sqrt{\dot{x}^2 + \dot{y}^2}dt$). I use Simpsons rule and since i already have a spline representation i can use as many integration points as i want. After the integration i know s(t) at the integration points. $\endgroup$
    – benno
    May 13, 2017 at 12:49
  • $\begingroup$ Then I create a spline for the inverse of this function: $t(s)$. (Since i have the values of t and s in two separate lists of the same size this is easy). When I now want to access x(s) i access x(t(s)). To evaluate the derivative i apply the chain rule. I think I get a more accurate arc length parametrization this way. I will now try to replace the first spline (uniform increasing parameter) with a parameter that increases by the segment length as a first guess. $\endgroup$
    – benno
    May 13, 2017 at 12:54
0
$\begingroup$

I found a more or less satisfying solution to the problem. The profile above was just something i tried after having the same problems with other "real" profiles (they were given to me in a csv file -> no analytic formula to reconstruct or redistribute the points).

The solution is as follows: I created a spline with the given points, then i redistributed the points with uniform distance (arc length). With the redistributed points I created a new spline. In the new the spline, the oscillations in the second derivative were almost gone.

I use the second derivative of the spline to determine the curvature and the normal vector.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.