Derivatives of Polylines

This demonstration shows some examples of computing curvature and torsion for polylines using standard formulas and the formulae of Langer, Belyaev and Seidel.

First, we import some data (note: you'll probably have to change the path for trefoil_data.dat to reproduce this on your system).

TrefoilTable = Import["/Users/cantarel/courses/math4500/trefoil_data.dat","Table"];

ListPointPlot3D[TrefoilTable,BoxRatios->Automatic]

We now compute the curvature function for this curve using some of the various estimators that we've covered so far in class. To do this, we need to set up a little code:

FindCDEF[i_,P_] := Module[{c,d,e,f,CC,DD, EE,FF},

CC = P[[i-1]] - P[[i-2]];
DD = P[[i]] - P[[i-1]];
EE = P[[i+1]] - P[[i]];
FF = P[[i+2]] - P[[i+1]];

c = Sqrt[CC.CC];
d = Sqrt[DD.DD];
e = Sqrt[EE.EE];
f = Sqrt[FF.FF];

{c,d,e,f,CC,DD,EE,FF}

];

kThree[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

phi = VectorAngle[DD,EE];

2*phi/(d + e)

];

kFive[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

(1/(c + 2(d+e) + f))((2e+f)kThree[i-1,P] + (c+2d)kThree[i+1,P])

];

kFivePrime[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

(3/(c + 2(d+e)+f))(kThree[i+1,P] - kThree[i-1,P])

];

We now try out the three and five point curvature estimators on our sample data.

Plotting the data shows us which one is less sensitive to the noise in the computation:

kThreeData = Table[kThree[i,TrefoilTable],{i,5,1195}];
kFiveData  = Table[kFive[i,TrefoilTable],{i,5,1195}];

Show[kThreePlot = ListLinePlot[{kThreeData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["Three Point Curvature Data",16],ImageSize->{640,480}]]

Show[kFivePlot = ListLinePlot[{kFiveData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["Five Point Curvature Data",16],ImageSize->{640,480}]]

What other methods should we compare this to? Rawdon defined the curvature of a polyline to be the radius of curvature of a circle arc spliced into the polyline. This is given by the formula:

kMinRad[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

phi = VectorAngle[DD,EE];

(2*Tan[phi/2])/Min[d,e]

];

kMinRadData  = Table[kMinRad[i,TrefoilTable],{i,5,1195}];

Show[kMinRadPlot = ListLinePlot[{kMinRadData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["MinRad Curvature Data",16],ImageSize->{640,480}]]

The differences between these datasets are rather subtle, so we view them in slideshow format:

kPlots = {kMinRadPlot,kThreePlot,kFivePlot};

SlideView[kPlots,ImageSize->{640,480}]

What about the torsion of this dataset? Well, torsion requires us to do a little bit of coding to implement the various elements of the LBS five-point estimator for torsion:

b[i_,P_] := Module[{Incoming,Outgoing},

Incoming = P[[i]] - P[[i-1]];
Outgoing = P[[i+1]] - P[[i]];
Cross[Incoming,Outgoing]/Norm[Cross[Incoming,Outgoing]]

];

t[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

(d e/(d+e))*((1/d^2)*DD + (1/e^2)*EE)

];

etaD[i_,P_] := Dot[Cross[b[i,P],b[i-1,P]],t[i,P]];

etaE[i_,P_] := Dot[Cross[b[i+1,P],b[i,P]],
t[i,P]];

tauD[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

3*etaD[i,P]/(c+d + e)

];

tauE[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

3*etaE[i,P]/(d + e + f)

];

tauHat[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

(1/(c+d+e+f))*((f + 2e - d)*tauD[i,P] + (c + 2d-e)*tauE[i,P])

];

tauFive[i_,P_] := Module[{c,d,e,f,CC,DD,EE,FF,phi},

{c,d,e,f,CC,DD,EE,FF} = FindCDEF[i,P];

tauHat[i,P] + ((c e - d f + d^2 - e^2)/(3*(c+d+e+f)))*(kFivePrime[i,P]/kFive[i,P]) * tauHat[i,P]

];

Now the functions tauD and tauE are four point (linear) approximations to the torsion at i, while tauFive is a five point (quadratic) approximation to the torsion at i. We go ahead and try to see the difference in the data.

tauDData  = Table[tauD[i,TrefoilTable],{i,5,1195}];

tauDPlot = ListLinePlot[{tauDData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauD torsion Data",16],ImageSize->{640,480}];

tauEData  = Table[tauE[i,TrefoilTable],{i,5,1195}];

tauEPlot = ListLinePlot[{tauEData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauE torsion Data",16],ImageSize->{640,480}];

tauFiveData  = Table[tauFive[i,TrefoilTable],{i,5,1195}];

tauFivePlot = ListLinePlot[{tauFiveData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauFive torsion Data",16],ImageSize->{640,480}];

SlideView[{tauDPlot,tauEPlot,tauFivePlot},ImageSize->{640,480}]

What we discover is that while the five point data certainly looks a little better than the three point data, basically all of the torsion data looks like soup. We conclude that the underlying polyline is not resolving the torsion data well enough to figure out what the function is supposed to look like.

We now want to test the curvature and torsion functions against a curve that we know the curvature and torsion of to get a sense of the speed of convergence (and to make sure that the above soup is not merely a function of a programming error).

r[t_] := {(2 + Cos[3*t/2])*Cos[t],(2 + Cos[3*t/2])*Sin[t], Sin[3*t/2]};

ParametricPlot3D[r[t],{t,0,6*Pi},BoxRatios-> Automatic]

We now define the curvature and torsion of the curve as a function using Mathematica's D operator for differentiation.

kappa[t_] := Module[{rp,rpp,rppp},

rp := D[r[x],x] /. x -> t;
rpp := D[r[x],x,x] /. x -> t;
rppp := D[r[x],x,x,x] /. x -> t;

Norm[Cross[rp,rpp]]/(Norm[rp]^3)

];


tau[t_] := Module[{rp,rpp,rppp},

rp := D[r[x],x] /. x -> t;
rpp := D[r[x],x,x] /. x -> t;
rppp := D[r[x],x,x,x] /. x -> t;

-Dot[Cross[rp,rpp],rppp]/(Norm[Cross[rp,rpp]]^2)

];

We can now generate some data for this curve by inscribing some points in the curve.

n := 1000;

R = Table[r[6.0*Pi*i/n],{i,0,n}];
newkThreeData = Table[{6.0*Pi*i/n,kThree[i,R]},{i,5,n-5}];

SlideView[{ListLinePlot[{newkThreeData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["kThree Data",16],ImageSize->{640,480}],Plot[kappa[t],{t,0,6*Pi},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["kappa Data",16],ImageSize->{640,480}]}]

We see from this that the curvature estimator seems to work as advertised. What about the torsion formula?

n := 1000;

R = Table[r[6.0*Pi*i/n],{i,0,n}];
newtauFiveData = Table[{6.0*Pi*i/n,tauFive[i,R]},{i,5,n-5}];

SlideView[{ListLinePlot[{newtauFiveData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauFive Data",16],ImageSize->{640,480}],Plot[tau[t],{t,6.0*Pi*5/n,6*Pi - 6.0*Pi*5/n},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tau Data",16],ImageSize->{640,480}]}]

   errorTable := Table[{6.0*Pi*i/n,tauFive[i,R] - tau[6.0*Pi*i/n]},{i,5,n-5}];
	
	ListLinePlot[{errorTable},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["errorTable",16],ImageSize->{640,480}]

We conclude that the estimators work as advertised, at least on noiseless curves. It's worth doing a final experiment to see whether we can see a clearer difference between the kThree and kFive estimators for noisier data.

n := 1000;
Noise := 10^(-6);

R = Table[r[6.0*Pi*i/n] + {RandomReal[{-Noise,Noise}],RandomReal[{-Noise,Noise}],RandomReal[{-Noise,Noise}]},{i,0,n}];

newkThreeData = Table[{6.0*Pi*i/n,kThree[i,R]},{i,5,n-5}];
newkFiveData = Table[{6.0*Pi*i/n,kFive[i,R]},{i,5,n-5}];

SlideView[{ListLinePlot[{newkThreeData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["kThree Data",16],ImageSize->{640,480}],

ListLinePlot[{newkFiveData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["kFive Data",16],ImageSize->{640,480}],

Plot[kappa[t],{t,0,6*Pi},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["kappa Data",16],ImageSize->{640,480}]}]

newtauFiveData = Table[{6.0*Pi*i/n,tauFive[i,R]},{i,5,n-5}];
newtauDData = Table[{6.0*Pi*i/n,tauD[i,R]},{i,5,n-5}];


SlideView[{ListLinePlot[{newtauDData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauD Data",16],ImageSize->{640,480}],

ListLinePlot[{newtauFiveData},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tauFive Data",16],ImageSize->{640,480}],

Plot[tau[t],{t,0,6*Pi},Filling->Axis,PlotRange->Full,AxesOrigin->{0,0},PlotLabel->Style["tau Data",16],ImageSize->{640,480}]}]

Experimenting with this shows how much more sensitive the torsion estimator is to noise in the data than the curvature estimator, even though both are in principle O(d,e)^2 in error.