(help) 3D Spline Interpolation

[malu05]

My math sucks and i have to gather code from around the web to get my stuff working.

Here is my 3D spline interpolation script. And it works just fine.

The problem is that with this method i have to define the tangent for each point.

expandcollapse popup

#include <GUIConstantsEx.au3>
#include <Misc.au3>
#Include <WinAPI.au3>
#include <Math.au3>

dim $angles[3]
$pi = 4 * ATan(1) ;equals 3.14159265358979
$radToDeg = 180 / $pi

dim $_p0x = -8901.369140625
dim $_p0y = -113.60546875
dim $_p0z = 82.1371917724609

dim $_p1x = -8903.7734375
dim $_p1y = -113.341896057129
dim $_p1z = 83.5231399536133

dim $_p2x = -8902.8359375
dim $_p2y = -112.653060913086
dim $_p2z = 83.6386947631836

dim $_p3x = -8901.927734375
dim $_p3y = -112.125366210938
dim $_p3z = 83.6458740234375

dim $_p4x = "na"
dim $_p4y = "na"
dim $_p4z = "na"

dim $_p5x = -8901.908203125
dim $_p5y = -119.156959533691
dim $_p5z = 83.7604446411133

dim $_p6x = -8899.4599609375
dim $_p6y = -120.456336975098
dim $_p6z = 83.8341217041016



local $P_ONE[3] = [$_p0x,$_p0y,$_p0z]
local $P_TWO[3] = [$_p1x,$_p1y,$_p1z]
local $P_TREE[3] = [$_p2x,$_p2y,$_p2z]
local $P_FOUR[3] = [$_p3x,$_p3y,$_p3z]
CubicPathing($P_ONE,$P_TWO,$P_TREE,$P_FOUR)

local $P_Pre[3] = [$_p2x,$_p2y,$_p2z]
local $P_ONE[3] = [$_p3x,$_p3y,$_p3z]
local $P_TWO[3] = [$_p4x,$_p4y,$_p4z]
local $P_TREE[3] = [$_p5x,$_p5y,$_p5z]
local $P_FOUR[3] = [$_p6x,$_p6y,$_p6z]
CubicPathing($P_ONE,$P_TWO,$P_TREE,$P_FOUR,$P_Pre)

func CubicPathing($P_ONE,$P_TWO,$P_TREE,$P_FOUR,$P_PRE="na")
local $p0x = $P_ONE[0]
local $p0y = $P_ONE[1]
local $p0z = $P_ONE[2]

local $p1x = $P_TWO[0]
local $p1y = $P_TWO[1]
local $p1z = $P_TWO[2]

local $p2x = $P_TREE[0]
local $p2y = $P_TREE[1]
local $p2z = $P_TREE[2]

local $p3x = $P_FOUR[0]
local $p3y = $P_FOUR[1]
local $p3z = $P_FOUR[2]


if $P_PRE = "na" Then
    ;nothing
Else
    local $pPx = $P_PRE[0]
local $pPy = $P_PRE[1]
local $pPz = $P_PRE[2]
    $p1x = $p0x-($pPx-$p0x)
    $p1y = $p0y-($pPy-$p0y)
    $p1z = $p0z-($pPz-$p0z)
EndIf
for $i = 0 to 1.001 step +0.001
    $_struc = "float posX;float posY;float posZ"
    $_strucT  = DllStructCreate($_struc)
    DllStructSetData($_strucT,"posX",(1-$i)^3 * $p0x + 3*(1-$i)^2 * $i * $p1x + 3*(1-$i)* $i^2 * $p2x + $i^3 * $p3x)
    DllStructSetData($_strucT,"posY",(1-$i)^3 * $p0y + 3*(1-$i)^2 * $i * $p1y + 3*(1-$i)* $i^2 * $p2y + $i^3 * $p3y)
    DllStructSetData($_strucT,"posZ",(1-$i)^3 * $p0z + 3*(1-$i)^2 * $i * $p1z + 3*(1-$i)* $i^2 * $p2z + $i^3 * $p3z)    
    ConsoleWrite(" X: " &DllStructGetData($_strucT,"posX")&" Y: " &DllStructGetData($_strucT,"posY")&" Z: " &DllStructGetData($_strucT,"posZ")&@crlf)
Next
EndFunc

What i wanted to do is just having to define the points that the spline will go through

Like so:

This brings me back to my math problem. I can read though articles like this http://paulbourke.net/miscellaneous/interpolation/ but i don't really understand it so i don't really know how i can turn it into autoit code.

Is there a easy way to calculate tangents for a smooth 3d interpolation or would anyone help me understand the article?

[jchd]

Welcome to interpolation!

Without control points, a spline isn't a spline!

You're doing a piecewise cubic interpolation.

There is, of course, an infinity of curves that will pass through your points so you have to force the choice of one by defining endpoints tangents.

There a huge number of algorithms available, each of them claiming to find the "best" interpolation, but each of them defining how they judge what "best" is.

Then you have splines, B-splines, beta-splines, Bezier, Rabenki, De Casteljau curves or surfaces, Lagrange, Hermite, Bessel, trigonometric polynomials, etc.

All of them rely on minimizing "something (generally a quadratic quantity) to determine the "best" curve or surface.

The litterature on these subjects is soo voluminous that you couldn't read it in your lifetime.

I'd advise you to Google for Paul Faget De Casteljau to find exposition of his "Formes à Pôles" (in french, probably "Polar Forms" in english) and relating algorithms.

I believe that you can derive a workable implementation for your "ideal" 3D curve, without ressorting to a concept of extremum nor requiring you to provide control points yourself (or tangents).

[malu05]

Welcome to interpolation!

Without control points, a spline isn't a spline!

You're doing a piecewise cubic interpolation.

There is, of course, an infinity of curves that will pass through your points so you have to force the choice of one by defining endpoints tangents.

There a huge number of algorithms available, each of them claiming to find the "best" interpolation, but each of them defining how they judge what "best" is.

Then you have splines, B-splines, beta-splines, Bezier, Rabenki, De Casteljau curves or surfaces, Lagrange, Hermite, Bessel, trigonometric polynomials, etc.

All of them rely on minimizing "something (generally a quadratic quantity) to determine the "best" curve or surface.

The litterature on these subjects is soo voluminous that you couldn't read it in your lifetime.

I'd advise you to Google for Paul Faget De Casteljau to find exposition of his "Formes à Pôles" (in french, probably "Polar Forms" in english) and relating algorithms.

I believe that you can derive a workable implementation for your "ideal" 3D curve, without ressorting to a concept of extremum nor requiring you to provide control points yourself (or tangents).

Thanks for the reply.

However as noted in my post my math skills is my problem here.

When i google stuff like Paul Faget De Casteljau's Formes à Pôles i get a bunch of math mumbo jumbo that i really can't translate to anything. What i needed was a bit help getting the formula down in code format so i could understand it.

[jchd]

What I don't get is why people insist trying to perform heavy numeric computations in AutoIt, often using vectors or N-dimensional matrices instead of using the tools which have been carefully developped and polished for performing such task very efficiently and using much more expressive languages when it comes to mathematical work.

If you want a quick result, why not use existing implementations wich use tools fitted to the task?

If you want to have your own hand at it and learn how this kind of job is done, you can use the built-in Lagrange function in SciLab:

function[P]=lagrange(X,Y)//X nodes,Y values;P is the numerical Lagrange polynomial interpolation
n=length(X);// n is the number of nodes. (n-1) is the degree
x=poly(0,"x");P=0;
for i=1:n, L=1;
  for j=[1:i-1,i+1:n] L=L*(x-X(j))/(X(i)-X(j));end
  P=P+L*Y(i); 
end
endfunction

First feed it the X and Y vectors to obtain a polynomial for X-->Y interpolation, P_y. Then repeat with X and Z vectors and obtain P_z. After that you just have to evaluate P_y(x) at any you need to find the corresponding y. Evaluation of P_z(x) at the same x gives you the z a this point.

And you're done.

Do it piecewise (no more than a handful of points at the time). Don't forget that n points give polynomials of degree n-1, so if you feed 300 points, you'll be dealing with x .. x²⁹⁹ with potentially large coefficients.

Using SciLab will give you "some" guarantee of stability (anyway much better than any naïve implementation in AutoIt) but De Casteljau algorithm is one of the most stable for evaluating large polynomials.

[malu05]

What I don't get is why people insist trying to perform heavy numeric computations in AutoIt, often using vectors or N-dimensional matrices instead of using the tools which have been carefully developped and polished for performing such task very efficiently and using much more expressive languages when it comes to mathematical work.

If you want a quick result, why not use existing implementations wich use tools fitted to the task?

If you want to have your own hand at it and learn how this kind of job is done, you can use the built-in Lagrange function in SciLab:

function[P]=lagrange(X,Y)//X nodes,Y values;P is the numerical Lagrange polynomial interpolation
n=length(X);// n is the number of nodes. (n-1) is the degree
x=poly(0,"x");P=0;
for i=1:n, L=1;
  for j=[1:i-1,i+1:n] L=L*(x-X(j))/(X(i)-X(j));end
  P=P+L*Y(i); 
end
endfunction

First feed it the X and Y vectors to obtain a polynomial for X-->Y interpolation, P_y. Then repeat with X and Z vectors and obtain P_z. After that you just have to evaluate P_y(x) at any you need to find the corresponding y. Evaluation of P_z(x) at the same x gives you the z a this point.

And you're done.

Do it piecewise (no more than a handful of points at the time). Don't forget that n points give polynomials of degree n-1, so if you feed 300 points, you'll be dealing with x .. x²⁹⁹ with potentially large coefficients.

Using SciLab will give you "some" guarantee of stability (anyway much better than any naïve implementation in AutoIt) but De Casteljau algorithm is one of the most stable for evaluating large polynomials.

Yea, i do agree that autoIt is not fit for this, but since its a part of my application i need it to be in the script and not use a external app to do the calculation.

The method i posted in the first post works as intended and the interpolation is flawless for its purpose. But some of my users find the tangents a bit intimidating. So its not directly a matter of finding another "better" algorithm but rather to avoid having to use tangents if that would be using a algorithm that don't use them or simply pre-calculate the best suited tangents.

[jchd]

Then "invent" extra endpoints and apply your cubic thing. You could just use some point in the segment P0 to P1 (before P0 and another in segment Pn-1 to Pn (after Pn). It's a bit crude but it could be enough for your needs.