Getting interpolated values at fractions between knots

[jkrumbiegel]

Hi Jeffrey,

I thought about adding your Catmull-Rom splines to my Animations.jl package, for things like camera trajectories through space for example. What I would need is a function to get the point that lies on the spline between two knots at a certain fraction. So f(0) would be knot 1, f(1) would be knot 2 and f(0.5) would be halfway between the two. If I understand correctly, the two functions that your package offers generate lots of points at equal spacing, but is my request also possible? Also, knowing the arc length between two points beforehand would be useful, in order to control animation speed separately.

Thanks for your help,
Julius

[JeffreySarnoff]

Happy to hear this. In src/centripetal/arcbased.jl there is this function:

 approx_catmullrom_arclength(pt0, pt1, pt2, pt3)

Given 4 ND points along a centripetal Catmull-Rom span,
roughly approximate the arclength of the curvilinear segment
that would be determined by the two bounding points
and the tangents they determine [the arc between p1 and p2].
This well-behaved approximation was developed by Jens Gravesen

I could use the same automagical generation of p0 and p3 if you really want to feed two points at a time .. for animation and most other continuous-like processes, providing the leading and trailing points as well as the two edge points gives cleaner results. Otherwise, an approximation is used.

"So f(0) would be knot 1, f(1) would be knot 2 and f(0.5) would be halfway between the two."
If you have available the preceding and following points p(a) and p(z) where p(0)=f(0) and p(1)=f(1) where p(a), a=0-delta1 p(z) z=1+delta2 where the deltas are small (e.g. 0.1) or match (e.g. 1.0) then your value f(0.5) will be more accurate.

I am happy to export / add + export what will work well for you. Could you show me an example of some successive data/function points [xs and f(x)s]?

[jkrumbiegel]

If I understand correctly, you always need at least four points to have a spline, that is fine with me per se, but I would still like to have the option to interpolate with different timings between pairs of knots, given that the overall spline has at least four points.

Let's say I have four camera positions and fit a spline through them, then maybe I want to first animate from point 1 to 2 quite fast, but then go slower to 3, etc. That's why I need to know the arc lengths between every pair of knots, as well as the ability to interpolate the positions at arbitrary fractions.

Does approx_catmullrom_arclength give the length of the full four-point spline? I don't quite understand the docstring, yet. For four control points I would expect 3 different arclengths.

[JeffreySarnoff]

approx_catmullrom_arclength gives the approximate length of the splined curve through the second and third points. It is a good approximation. The first and fourth points are used to better determine the curve, and so obtain a more accurate estimate.

[JeffreySarnoff]

In the README, there is a description of using interpoint spacing to govern the number of intermediating points [arclength relative allocations] (https://github.com/JeffreySarnoff/CatmullRom.jl#arclength-relative-allocation). Is this the sort of speed control that you are seeking? If not, what do you want to do differently?

[JeffreySarnoff]

I can give you the ability to interpolate for whatever % of the interknot distances you choose. That requires some coding. If I understand correctly, you want to give a sequence of knots and for each adjacent pair obtain the arc distance between them and the splined coordinate halfway between them. It is easy to obtain a halfway-x by linear interpolation, it is less easy to move halfway along the spline and then return the corresponding x.

[jkrumbiegel]

Halfway was just a random example, I need any percentage between two knots. Depending on the number of frames for an animation and the chosen easing function, the interpolated points will be spaced apart in arbitrary ways, but always be specified by their fraction along the spline. That's because the speed of a dot moving along the spline depends on the distance it travels, so depending on the easing function the animation uses I might need to calculate positions for each entry in let's say a vector of [0, 0.1, 0.25, 0.6, 0.9, 1] spline length fractions. That could be one animation.

[JeffreySarnoff]

That is no problem. How about linear interpolation for the xs?

[JeffreySarnoff]

Actually I get a better approximation using the spline length fractions to generate new knots and use those knots to find a better than linear approximation to the xs.

[jkrumbiegel]

which xs do you mean? the x coordinates of my knots? I don't know much about the mathematics behind the splines if you're referring to something there

[JeffreySarnoff]

I do mean the x coordinates of your knots and the x coordinates of you spline length fractions.
Each spline length fraction needs to become coordinates for the corresponding interpolated knot.

SplineLengthCoords_xy = splineat(0.1, spline) ...

[JeffreySarnoff]

Please show me a short set of coordinates (your initial points, in sequence) .. 6 is enough. You can estimate them, or make them up if that easier.

[JeffreySarnoff]

note to myself:
https://www.habrador.com/tutorials/interpolation/3-move-along-curve/
https://www.geometrictools.com/Documentation/MovingAlongCurveSpecifiedSpeed.pdf
https://gamedev.stackexchange.com/questions/14985/determine-arc-length-of-a-catmull-rom-spline-to-move-at-a-constant-speed
https://gamedev.stackexchange.com/questions/47354/catmull-rom-spline-constant-speed
https://easings.net/en

using LinearAlgebra: norm
ccoords=[[Float64.(x)...,] for x in cameracoords];
xyzchange=norm.(diff(ccoords))[1:end]
relxyzchange = xyzchange ./ sum(xyzchange)
cumrelxyzchange = cumsum(relxyzchange)
pushfirst!(cumrelxyzchange, 0.0)
m = vcat(cumrelxyzchange',reduce(hcat,ccoords))
m = permutedims(m)
wxyz = [m[i,:] for i=1:size(m)[1]]

m
7×4 Array{Float64,2}:
 0.0       0.0   0.0  0.0
 0.381936  1.0  10.0  0.0
 0.507982  2.0  11.0  3.0
 0.561727  1.0  12.0  3.0
 0.67574   3.0  13.0  1.0
 0.76883   1.0  14.0  0.0
 1.0       0.0  20.0  0.0

wxyz
 [0.0, 0.0, 0.0, 0.0]
 [0.38193624872274273, 1.0, 10.0, 0.0]
 [0.5079815132220372, 2.0, 11.0, 3.0]
 [0.5617273945910012, 1.0, 12.0, 3.0]
 [0.6757396261215276, 3.0, 13.0, 1.0]
 [0.7688302233501427, 1.0, 14.0, 0.0]
 [1.0, 0.0, 20.0, 0.0]


julia> catmullrom(wxyz,1)
4-element Array{Array{Float64,1},1}:
 [0.0, 0.13902702210457504, 0.27195110540413064, 0.38193624872274273, 0.43095816546243493, 0.4724098534875167, 0.5079815132220372, 0.5268285052860687, 0.5431682448626223, 0.5617273945910012, 0.5970204225643049, 0.6370215304433786, 0.6757396261215276, 0.7059266641271809, 0.7350661270574904, 0.7688302233501427, 0.8404011484418478, 0.923456321065393, 1.0]
 [0.0, 0.3041489336458593, 0.5998005980466854, 1.0, 1.4237692481311095, 1.8820409748215985, 2.0, 1.6918442707242118, 1.2385787179538434, 1.0, 1.5270065101988246, 2.5094325705883724, 3.0, 2.561801606073483, 1.7220023216537015, 1.0, 0.45857841019141643, 0.20967935631405243, 0.0]           
 [0.0, 3.6260658318269923, 7.395893156899611, 10.0, 10.565991214326147, 10.747275433770806, 11.0, 11.322475943300562, 11.675414152823436, 12.0, 12.36927835945959, 12.678412275885982, 13.0, 13.245818739923497, 13.50180597124347, 14.0, 16.03602251031095, 18.673511687861996, 20.0]          
 [0.0, -0.24568865083301272, -0.49137730166602545, 0.0, 0.97227355948038, 2.187791045785854, 3.0, 3.1560552178212045, 3.1401390997340077, 3.0, 2.4742498433323306, 1.6849829485596453, 1.0, 0.5811444467981155, 0.24743637912482885, 0.0, -0.21745181858467416, -0.2206671023668562, 0.0]       

julia> floatvecs(seq_of_coords, ::Type{T}=Float64) where {T} = map(x->[T.(x)...,], seq_of_coords)
floatvecs (generic function with 2 methods)

julia> cameracoords
7-element Array{Tuple{Int64,Int64,Int64},1}:
 (0, 0, 0)
 (0, 10, 0)
 (0, 11, 3)
 (1, 12, 3)
 (3, 13, 1)
 (0, 14, 0)
 (0, 20, 0)

julia> cameraxyz = floatvecs(cameracoords)
7-element Array{Array{Float64,1},1}:
 [0.0, 0.0, 0.0]
 [0.0, 10.0, 0.0]
 [0.0, 11.0, 3.0]
 [1.0, 12.0, 3.0]
 [3.0, 13.0, 1.0]
 [0.0, 14.0, 0.0]
 [0.0, 20.0, 0.0]

julia> catmullrom_polys = CatmullRom.catmullrom_polys
catmullrom_polys (generic function with 1 method)

julia> polys = catmullrom_polys(cameraxyz)
4×3 Array{Polynomials.Poly{Float64},2}:
 Poly(-0.896094*x^2 + 0.896094*x^3)                  …  Poly(1.92019*x + 3.95737*x^2 - 2.87757*x^3)
 Poly(0.599254*x + 0.649575*x^2 - 0.24883*x^3)          Poly(3.0 + 0.803984*x - 1.04896*x^2 + 0.244981*x^3)
 Poly(1.0 + 1.67774*x + 3.01027*x^2 - 2.68801*x^3)      Poly(3.0 - 0.814175*x - 2.88296*x^2 + 1.69714*x^3)
 Poly(3.0 - 0.384566*x - 6.51018*x^2 + 3.89474*x^3)     Poly(1.0 - 1.56528*x + 0.704117*x^2 - 0.138841*x^3)

julia> using Polynomials

julia> """
           catmullrom_polys(points)
       Traverse a sequence of points (coordinates in 2,3,..n dims),
       obtaining interpolatory centripetal Catmull-Rom polynomials.
       Yields separate sequences for each coordinate dimension,
           each sequence is ordered from first to last interpoint span.
       Each poly within a sequence fits an interpoint span with respect to
           the sequence's coordinate axis.
       _polys_ is an `(n_points - 3) x (n_coords) 2D array`
          columns are of one coordinate axis
          each row holds polys that interpolate one interpoint span across all coordinates.
       """
"    catmullrom_polys(points)\nTraverse a sequence of points (coordinates in 2,3,..n dims),\nobtaining interpolatory centripetal Catmull-Rom polynomials.\nYields separate sequences for each coordinate dimension, \n    each sequence is ordered from first to last interpoint span.\nEach poly within a sequence fits an interpoint span with respect to\n    the sequence's coordinate axis.\n_polys_ is an `(n_points - 3) x (n_coords) 2D array`\n   columns are of one coordinate axis \n   each row holds polys that interpolate one interpoint span across all coordinates.\n"

julia> first_xyz_polys = polys[1,:]
3-element Array{Poly{Float64},1}:
 Poly(-0.896094363372106*x^2 + 0.896094363372106*x^3)
 Poly(10.0 + 2.664128249735226*x - 3.625096444993443*x^2 + 1.9609681952582172*x^3)
 Poly(1.9201949994086551*x + 3.957373254730034*x^2 - 2.877568254138689*x^3)

julia> first_xyz_polys_at_half = map(polyval(ply, 0.5), first_xyz_polys)
ERROR: UndefVarError: ply not defined
Stacktrace:
 [1] top-level scope at REPL[93]:1

julia> first_xyz_polys_at_half = map(ply->polyval(ply, 0.5), first_xyz_polys)
3-element Array{Float64,1}:
 -0.11201179542151325
 10.67091103802653
  1.5897447816194998

julia> second_xyz_polys = polys[2,:]
3-element Array{Poly{Float64},1}:
 Poly(0.5992544178491149*x + 0.6495750861548379*x^2 - 0.24882950400395276*x^3)
 Poly(11.0 + 0.8672491406746515*x + 0.3930873492090886*x^2 - 0.26033648988374003*x^3)
 Poly(3.0 + 0.8039841684766105*x - 1.0489649195425734*x^2 + 0.24498075106596284*x^3)

julia> second_xyz_polys_at_quarter = map(ply->polyval(ply, 0.25), second_xyz_polys)
3-element Array{Float64,1}:
  0.18652408634689432
 11.237312486839798
  3.1392635588831475

julia> xs, ys, zs = zip(first_xyz_polys_at_half, second_xyz_polys_at_quarter)
Base.Iterators.Zip{Tuple{Array{Float64,1},Array{Float64,1}}}(([-0.11201179542151325, 10.67091103802653, 1.5897447816194998], [0.18652408634689432, 11.237312486839798, 3.1392635588831475]))

julia> collect(ans)
3-element Array{Tuple{Float64,Float64},1}:
 (-0.11201179542151325, 0.18652408634689432)
 (10.67091103802653, 11.237312486839798)
 (1.5897447816194998, 3.1392635588831475)

julia> xs, ys, zs = zip(first_xyz_polys_at_half, second_xyz_polys_at_quarter)
Base.Iterators.Zip{Tuple{Array{Float64,1},Array{Float64,1}}}(([-0.11201179542151325, 10.67091103802653, 1.5897447816194998], [0.18652408634689432, 11.237312486839798, 3.1392635588831475]))

julia> floatvecs(ans)
3-element Array{Array{Float64,1},1}:
 [-0.11201179542151325, 0.18652408634689432]
 [10.67091103802653, 11.237312486839798]
 [1.5897447816194998, 3.1392635588831475]

julia> xs
(-0.11201179542151325, 0.18652408634689432)

julia> ys
(10.67091103802653, 11.237312486839798)

julia> zs
(1.5897447816194998, 3.1392635588831475)

[jkrumbiegel]

I would do something like this:

cameracoords = [
    (0, 0, 0),
    (0, 10, 0),
    (0, 11, 3),
    (3, 12, 3),
    (3, 13, 0),
    (0, 14, 0),
    (0, 20, 0),
]

anim = CatmullAnimation(cameracoords, timestamps_for_knots, easingfunction)

dt = 1/30
timestamps = 0:dt:5
positions = at.(anim, timestamps)

This would be going up, then something like an upward spiral and then up again. So I would expect one spline to go through all these points, but the spline shape doesn't yet determine the speed of the animation. This I would control with easings functions and time stamps.

[JeffreySarnoff]

In your example, timestamps_for_knots are one timestamp for each cameracoord and successive cameracoords have increasing timestamps?

[jkrumbiegel]

yes I think that's how I would implement it

[JeffreySarnoff]

There is a bug that appears using this with Julia v1.2+. I have to fix that first. Standby :).

[JeffreySarnoff]

Fixed. (Actually your input brought something to light, so thanks for that).

[jkrumbiegel]

that's good :) thank you, so does that mean there's now a function that I can use or was this more of a preliminary step?

[JeffreySarnoff]

preliminary, necessary. I need to play around a little and see what is easiest to use, then export that.

[JeffreySarnoff]

You could start playing with the following. It will become cleaner for use.

pkg> up
pkg> add Polynomials
pkg> add CatmullRom

you will need this function (it will be in the next revision of the package)
coordvecs(coordseq, ::Type{T}=Float64) where {T<} = map(x->[T.(x)...,], coordseq)
It does this

vec_of_coords = [(1, 2), (2, 4)]
coordvecs(vec_of_coords) = [[1.0, 2.0], [2.0, 4.0]]

tup_of_coords = ((1, 2), (2, 4))
coordvecs(tup_of_coords, Float32) = ([1.0f0, 2.0f0], [2.0f0, 4.0f0])

You want to work with the polynomial expressions of the spline segments directly.
catmullrom_polys = CatmullRom.catmullrom_polys

catmullrom_polys takes a sequence of coords (as prepared with coordvecs)
and returns the polynomials that spline between neighboring coords (neighbors with respect to the order in which the coords are given, so there will be a polynomial spline connecting (coords[i], coords[i+1]) for i=2:length(coords)-1. There will be size(coords)[1] - 3 because one point is used twice and the endpoints do not have connections on both sides.

If you have 12 [x,y,z] coordinates, catmullrom_polys(coords) will return (12-3==9, length([x,y,z])==3) a 9x3 array of polynomials. Each column maps to a coordinate, (here x, y, or z), and each row maps to an adjacent pair of coords. An example is given below.

using Polynomials
using CatmullRom

coordvecs(coordseq, ::Type{T}=Float64) where {T<:AbstractFloat} = map(x->[T.(x)...,], coordseq)

catmullrom_polys = CatmullRom.catmullrom_polys

camera_a_coords = [
           (0, 0, 0),
           (0, 10, 0),
           (0, 11, 3),
           (1, 12, 3),
           (3, 13, 1),
           (0, 14, 0),
           (0, 20, 0),
       ];
a_coords = coordvecs(camera_a_coords);

camera_b_coords = [
           (0, 0, 0),
           (1, 5, 2),
           (2, 11, 3),
           (1,  8, 3),
           (3, 11, 1),
           (1, 14, 3),
           (0, 10, 2),
       ];
b_coords = coordvecs(camera_b_coords);

const X = 1
const Y = 2
const Z = 3

# get polynomials for each coordinate dimension
# a_polys_x is a vector of length `size(a_polys)[1]` that holds the x-coordinate splines
a_polys_x, a_polys_y, a_polys_z = [a_polys[:, i] for i in (X,Y,Z)];

# get polynomials for each dimension of a 3D spline
# a_polys_2to3 is a vector of length `size(a_polys)[2]` that holds the x_spline, y_spine, z_spline
# that joins the second coordinate with the third coordinate
a_polys_2to3 = a_polys[1, :]
a_polys_2to3_x, a_polys_2to3_y, a_polys_2to3_z = a_polys_2to3

# now you can check the interpolation (don't worry about the indexing calcs)
function Polynomials.polyval(plys, value)
    map(s->polyval.(s, value), map(i->plys[i,:], 1:size(plys)[2]))
end

# we can do something to provide the omitted penultimate spline
polyval(a_polys, 0.0) # omits spline from nowhere to first
polyval(a_polys, 1.0) # omits spline from last to nowhere and the prior spline 
a_coords

# to evaluate at a given value 0.0..1.0
# use polyval with whichever subset of the polys is appropriate

[jkrumbiegel]

Ok I've tried this out a little bit in 2D so I could plot what's going on. I don't yet know how to access the parts of the splines connecting the 1st and 2nd as well as n-1th and nth spline. Is there a trick to getting these polys?

Now I see that the spacing of uniformly spaced values fed into the polys is not uniform anymore, I think I would want uniform so that the speed on the spline doesn't depend on its curvature. I'll have to fiddle with this a bit.

This is what I tried:

twod_coords = [
    (0, 0),
    (2, 0),
    (2, 1),
    (1, 0.5),
    (0.5, 2)]

coordvec = coordvecs(twod_coords)

polys = catmullrom_polys(coordvec)

PyPlot.close_figs()
figure()
scatter(first.(twod_coords), (x->x[2]).(twod_coords))
x_poly1 = polyval(polys[1, 1], 0:0.1:1)
y_poly1 = polyval(polys[1, 2], 0:0.1:1)
scatter(x_poly1, y_poly1, alpha=0.5)

x_poly2 = polyval(polys[2, 1], 0:0.1:1)
y_poly2 = polyval(polys[2, 2], 0:0.1:1)
scatter(x_poly2, y_poly2, alpha=0.5)

gcf()

[JeffreySarnoff]

good progress -- there is a function to help you with the extremal splines, but first: what you want with respect to uniform speed .. is that uniform distances along the x axis, uniform distances within each of splines, uniform distances over all of the splines (along the larger constructed curve through the data), or something else?

[jkrumbiegel]

I mean uniform movement along the splines, like a car travelling with constant speed on a winding road. If I have uniform movement on each partial spline I also have uniform movement along the full spline. I would preferably want both versions, depending on the use case.

[jkrumbiegel]

That is, I have the full spline if I can generate uniformly spaced values on each partial spline and know the length of each partial spline as well

[JeffreySarnoff]

I have been traveling and will travelling through next week. I am looking into supporting your need, although I do not expect to do much until I return.