POV-Ray : Newsgroups : povray.text.scene-files : Even spacing of points along natural cubic splines Server Time19 May 2024 15:43:37 EDT (-0400)
 Even spacing of points along natural cubic splines (Message 1 to 2 of 2)
 From: Tor Olav Kristensen Subject: Even spacing of points along natural cubic splines Date: 12 Dec 2010 09:30:00 Message:
```
{
"@context": "https://schema.org",
"@type": "DiscussionForumPosting",
"headline": "Even spacing of points along natural cubic splines",
"dateCreated": "2010-12-12T14:30:00+00:00",
"datePublished": "2010-12-12T14:30:00+00:00",
"author": {
"@type": "Person",
"name": "Tor Olav Kristensen"
}
}
I wrote the code below in order to solve problems similar to what is described

From: yoyodunno
Subject: help with splines
Date: 9th Dec 2010

--
Tor Olav

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10
// POV-Ray SDL file for even spacing of points along natural cubic splines
// Copyright 2010 by Tor Olav Kristensen
// Email: t o r . o l a v . k [_A_T_] g m a i l . c o m
// http://subcube.com
// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#include "colors.inc"
#include "functions.inc"

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#local vDirection = vnormalize(pEnd - pStart);

union {
difference {
cone { <0, 0, 0>, 0, -(A + B)*vDirection, Radius + C }
translate pEnd
}
cylinder { pStart, pEnd - A*vDirection, Radius }
}

#end // macro Arrow

#macro Size1A(Array1A)

dimension_size(Array1A, 1)

#end // macro Size1A

#macro ExtractComponent1A(Vectors1A, v0)

#local SizeU = Size1A(Vectors1A);
#local Component1A = array[SizeU];

#local I = 0;
#while (I < SizeU)
#local Component1A[I] = vdot(Vectors1A[I], v0);
#local I = I + 1;
#end // while

Component1A

#end // macro ExtractComponent1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#macro TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z) {
d1_X/
f_r(d1_X, d1_Y, d1_Z)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_)
)
}

#end // macro TangentFunctionX_1A

#macro TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z) {
d1_Y/
f_r(d1_X, d1_Y, d1_Z)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_)
)
}

#end // macro TangentFunctionY_1A

#macro TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z) {
d1_Z/
f_r(d1_X, d1_Y, d1_Z)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_)
)
}

#end // macro TangentFunctionZ_1A

#macro BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn,
d2_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
(d1_Y*d2_Z - d1_Z*d2_Y)/
f_r(
d1_Y*d2_Z - d1_Z*d2_Y,
d1_Z*d2_X - d1_X*d2_Z,
d1_X*d2_Y - d1_Y*d2_X
)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_),
d2_X_Fn(u_),
d2_Y_Fn(u_),
d2_Z_Fn(u_)
)
}

#end // macro BinormalFunctionX_1A

#macro BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn,
d2_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
(d1_Z*d2_X - d1_X*d2_Z)/
f_r(
d1_Y*d2_Z - d1_Z*d2_Y,
d1_Z*d2_X - d1_X*d2_Z,
d1_X*d2_Y - d1_Y*d2_X
)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_),
d2_X_Fn(u_),
d2_Y_Fn(u_),
d2_Z_Fn(u_)
)
}

#end // macro BinormalFunctionY_1A

#macro BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn,
d2_Z_Fn)

#local TempFn =
function(d1_X, d1_Y, d1_Z, d2_X, d2_Y, d2_Z) {
(d1_X*d2_Y - d1_Y*d2_X)/
f_r(
d1_Y*d2_Z - d1_Z*d2_Y,
d1_Z*d2_X - d1_X*d2_Z,
d1_X*d2_Y - d1_Y*d2_X
)
}

function(u_) {
TempFn(
d1_X_Fn(u_),
d1_Y_Fn(u_),
d1_Z_Fn(u_),
d2_X_Fn(u_),
d2_Y_Fn(u_),
d2_Z_Fn(u_)
)
}

#end // macro BinormalFunctionZ_1A

#macro NormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

#local B_FnY_1A = BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local B_FnZ_1A = BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local T_FnY_1A = TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
#local T_FnZ_1A = TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

function(u_) {
B_FnY_1A(u_)*T_FnZ_1A(u_) - B_FnZ_1A(u_)*T_FnY_1A(u_)
}

#end // macro NormalFunctionX_1A

#macro NormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

#local B_FnX_1A = BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local B_FnZ_1A = BinormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local T_FnX_1A = TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
#local T_FnZ_1A = TangentFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

function(u_) {
B_FnZ_1A(u_)*T_FnX_1A(u_) - B_FnX_1A(u_)*T_FnZ_1A(u_)
}

#end // macro NormalFunctionY_1A

#macro NormalFunctionZ_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn, d2_Y_Fn, d2_Z_Fn)

#local B_FnX_1A = BinormalFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local B_FnY_1A = BinormalFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn, d2_X_Fn,
d2_Y_Fn, d2_Z_Fn);
#local T_FnX_1A = TangentFunctionX_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);
#local T_FnY_1A = TangentFunctionY_1A(d1_X_Fn, d1_Y_Fn, d1_Z_Fn);

function(u_) {
B_FnX_1A(u_)*T_FnY_1A(u_) - B_FnY_1A(u_)*T_FnX_1A(u_)
}

#end // macro NormalFunctionZ_1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#macro IntegrateTrapezoidFunction(Fn, A, B, M)

#local D = (B - A)/M;

function(n) {
select(n - 1, 0, ((Fn(A) + Fn(A + n*D))/2 + sum(i, 1, n - 1, Fn(A +
i*D)))*D)
}

#end // macro IntegrateTrapezoidFunction

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#macro NaturalCubicSplineCoefficients1A3A(tt, yy)

#local Size = Size1A(tt);
#local N = Size - 1;
#local hh = array[Size];
#local uu = array[Size];
#local bb = array[Size];
#local vv = array[Size];
#local zz = array[Size];
#local kk = array[Size][3];

#local I = 0;
#while (I <= N - 1)
#local hh[I] = tt[I+1] - tt[I];
#local kk[I][1] = 6*(yy[I+1] - yy[I])/hh[I];
#local I = I + 1;
#end // while

#local uu[1] = 2*(hh[0] + hh[1]);
#local vv[1] = kk[1][1] - kk[0][1];
#local I = 2;
#while (I <= N - 1)
#local uu[I] = 2*(hh[I] + hh[I-1]) - pow(hh[I-1], 2)/uu[I-1];
#local vv[I] = kk[I][1] - kk[I-1][1] - hh[I-1]*vv[I-1]/uu[I-1];
#local I = I + 1;
#end // while

#local zz[N] = 0;
#local I = N - 1;
#while (I >= 1)
#local zz[I] = (vv[I] - hh[I]*zz[I+1])/uu[I];
#local I = I - 1;
#end // while
#local zz[0] = 0;

#local I = 0;
#while (I <= N - 1)
#local kk[I][0] = (zz[I+1] - zz[I])/(6*hh[I]);
#local kk[I][1] = zz[I]/2;
#local kk[I][2] = -hh[I]/6*(zz[I+1] + 2*zz[I]) + (yy[I+1] - yy [I])/hh[I];
#local I = I + 1;
#end // while

kk

#end // macro NaturalCubicSplineCoefficients1A3A

#macro NaturalCubicSplineFunction1A(tt, yy, kk)

#local Size = Size1A(tt);
#local N = Size - 1;
#local ff = array[Size];

#local I = 0;
#while (I <= N - 1)
#local ff[I] =
function(t_) {
yy[I] + (t_ - tt[I])*(kk[I][2] + (t_ - tt[I])*(kk[I][1] + (t_ -
tt[I])*kk[I][0]))
}
#local I = I + 1;
#end // while

function(t_) {
select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
#local I = 0;
#while (I <= N - 1)
select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
#local I = I + 1;
#end // while
select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
}

#end // macro NaturalCubicSplineFunction1A

#macro Der1_NaturalCubicSplineFunction1A(tt, kk)

#local Size = Size1A(tt);
#local N = Size - 1;
#local ff = array[Size];

#local I = 0;
#while (I <= N - 1)
#local ff[I] =
function(t_) {
kk[I][2] + (t_ - tt[I])*(2*kk[I][1] + (t_ - tt[I])*3*kk[I][0])
}
#local I = I + 1;
#end // while

function(t_) {
select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
#local I = 0;
#while (I <= N - 1)
select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
#local I = I + 1;
#end // while
select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
}

#end // macro Der1_NaturalCubicSplineFunction1A

#macro Der2_NaturalCubicSplineFunction1A(tt, kk)

#local Size = Size1A(tt);
#local N = Size - 1;
#local ff = array[Size];

#local I = 0;
#while (I <= N - 1)
#local ff[I] =
function(t_) {
2*kk[I][1] + (t_ - tt[I])*6*kk[I][0]
}
#local I = I + 1;
#end // while

function(t_) {
select(t_ - tt[0], 1, 0, 0)*ff[0](t_) +
#local I = 0;
#while (I <= N - 1)
select(t_ - tt[I], 0, select(t_ - tt[I+1], 1, 0))*ff[I](t_) +
#local I = I + 1;
#end // while
select(t_ - tt[N], 0, 1, 1)*ff[N-1](t_)
}

#end // macro Der2_NaturalCubicSplineFunction1A

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#declare Times =
array[5] {
-2,
-1,
0,
0.5,
1
}

#declare Points =
array[5] {
<-2, 0, 40>,
<0, 0, -2>,
<1, 0, 2>,
<0, 6, 2>,
<8, 2, 2>
}

#declare NaturalSpline =
spline {
natural_spline
Times[0], Points[0]
Times[1], Points[1]
Times[2], Points[2]
Times[3], Points[3]
Times[4], Points[4]
}

#declare Start = -2;
#declare End = 1;
#declare Span = End - Start;

#declare Intervals = 200;

#declare I = 0;
#while (I < Intervals)
#declare T = Start + I/Intervals*Span;
sphere {
NaturalSpline(T), 0.04
pigment { color White }
}
#declare I = I + 1;
#end // while

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

#declare CompX = ExtractComponent1A(Points, x);
#declare CompY = ExtractComponent1A(Points, y);
#declare CompZ = ExtractComponent1A(Points, z);

#declare CoeffsX = NaturalCubicSplineCoefficients1A3A(Times, CompX)
#declare CoeffsY = NaturalCubicSplineCoefficients1A3A(Times, CompY)
#declare CoeffsZ = NaturalCubicSplineCoefficients1A3A(Times, CompZ)

#declare FnX = NaturalCubicSplineFunction1A(Times, CompX, CoeffsX);
#declare FnY = NaturalCubicSplineFunction1A(Times, CompY, CoeffsY);
#declare FnZ = NaturalCubicSplineFunction1A(Times, CompZ, CoeffsZ);

#declare Der1_FnX = Der1_NaturalCubicSplineFunction1A(Times, CoeffsX);
#declare Der1_FnY = Der1_NaturalCubicSplineFunction1A(Times, CoeffsY);
#declare Der1_FnZ = Der1_NaturalCubicSplineFunction1A(Times, CoeffsZ);

#declare Der2_FnX = Der2_NaturalCubicSplineFunction1A(Times, CoeffsX);
#declare Der2_FnY = Der2_NaturalCubicSplineFunction1A(Times, CoeffsY);
#declare Der2_FnZ = Der2_NaturalCubicSplineFunction1A(Times, CoeffsZ);

#declare T_FnX = TangentFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ);
#declare T_FnY = TangentFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ);
#declare T_FnZ = TangentFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ);

#declare B_FnX = BinormalFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);
#declare B_FnY = BinormalFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);
#declare B_FnZ = BinormalFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);

#declare N_FnX = NormalFunctionX_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);
#declare N_FnY = NormalFunctionY_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);
#declare N_FnZ = NormalFunctionZ_1A(Der1_FnX, Der1_FnY, Der1_FnZ, Der2_FnX,
Der2_FnY, Der2_FnZ);

#declare Arrows =
union {
object { Arrow(0*x, 4*x, 0.05) pigment { color Red } }
object { Arrow(0*y, 4*y, 0.05) pigment { color Green } }
object { Arrow(0*z, 4*z, 0.05) pigment { color Blue } }
}

Arrows

#declare Intervals = 300;

#declare LengthFn =
IntegrateTrapezoidFunction(
function(t_) { f_r(Der1_FnX(t_), Der1_FnY(t_), Der1_FnZ(t_)) }
Start, End,
Intervals
)

#declare TotalLength = LengthFn(Intervals);

#declare Knots = Intervals + 1;
#declare tt_ = array[Knots]
#declare yy_ = array[Knots]

#declare I = 0;
#while (I < Knots)
#declare tt_[I] = Start + I/Intervals*Span;
#declare yy_[I] = Start + LengthFn(I)/TotalLength*Span;
#declare I = I + 1;
#end // while

#declare kk_ = NaturalCubicSplineCoefficients1A3A(yy_, tt_)
#declare InverseSplineFn = NaturalCubicSplineFunction1A(yy_, tt_, kk_)

#declare SmallArrows =
object {
Arrows
scale <1, 1, 1>/6
}

#declare Intervals = 400;
#declare Knots = Intervals + 1;

#declare I = 0;
#while (I < Knots)
#declare T_ = Start + I/Intervals*Span;
#declare T = InverseSplineFn(T_);
//  #declare T = T_; // Uncomment this line to see it without correction
object {
SmallArrows
matrix <
T_FnX(T), T_FnY(T), T_FnZ(T),
B_FnX(T), B_FnY(T), B_FnZ(T),
N_FnX(T), N_FnY(T), N_FnZ(T),
FnX(T),   FnY(T),   FnZ(T)
>
}
#declare I = I + 1;
#end // while

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10

light_source { <-1, 6, -2>*5 color White }

camera {
location <1, 1, -1>*8
look_at <2, 1, 2>
}

// ===== 1 ======= 2 ======= 3 ======= 4 ======= 5 ======= 6 ======= 7 ======= 8
======= 9 ======= 10
```
 From: Tor Olav Kristensen Subject: Re: Even spacing of points along natural cubic splines Date: 12 Dec 2010 10:25:00 Message:
```
{
"@context": "https://schema.org",
"@type": "DiscussionForumPosting",
"@id": "#web.4d04e86cfc787f1ac734aecd0%40news.povray.org",
"headline": "Re: Even spacing of points along natural cubic splines",
"dateCreated": "2010-12-12T15:25:00+00:00",
"datePublished": "2010-12-12T15:25:00+00:00",
"author": {
"@type": "Person",
"name": "Tor Olav Kristensen"
}
}
"Tor Olav Kristensen" <tor### [at] TOBEREMOVEDgmailcom> wrote:
> I wrote the code below in order to solve problems similar to what is described
>
>
> From: yoyodunno
> Subject: help with splines
> Date: 9th Dec 2010

Sorry for the line breaks.

See the attached file.

--
Tor Olav
```

Attachments: