|
 |
I wrote the code below in order to solve problems similar to what is described
in this thread:
http://news.povray.org/povray.general/thread/%3Cweb.4d01143d847db539663718460%40news.povray.org%3E/
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
#macro Arrow(pStart, pEnd, Radius)
#local A = Radius*3;
#local B = Radius*2;
#local C = Radius;
#local vDirection = vnormalize(pEnd - pStart);
union {
difference {
cone { <0, 0, 0>, 0, -(A + B)*vDirection, Radius + C }
cone { -A*vDirection, Radius, -(A + 2*B)*vDirection, Radius + 2*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
Post a reply to this message
|
|
com> wrote: