I am doing 3D interpolation with the nodes being [x,y,z], where x is a 1*nx vector, y is 1*ny,1*nx, the corresponding values to be interpolated is V which is a nx*ny*nz matrix. I use Dierckx for 1d/2d interpolation, but it does not allow for 3d version. My nodes are not uniform. For example, x = [1.2, 1.5, 1.9, 2.9], y = [1.4, 1.6],z = [3.,3.6,4.]. GridInterpolations seem to work, but it cannot evaluate multiple points at the same time after the interpolation step. I prefer to use piecewise cubic spline. I wonder if anyone knows what package suits my need.
What do you mean by “at the same time”? Can’t you just broadcast the interpolation call or use a comprehension if you want to apply it to a vector of points?
I meant something like broadcasting. My understanding is that it’s not straightforward to broadcast the interpolation functions in GridInterpolations. I think I’ll write a loop. Thanks!
You could try something like this:
using GridInterpolations
grid = RectangleGrid(0.:2.0, 0.:2.0, 0.:2.0) # rectangular grid
gridData = rand(3, 3, 3)
x = [0.1 0.2 0.3; 0. 1.5 2.0; 1.3 1.6 0.0; 0.2 0.5 1.0]
interpolate.(Ref(grid), Ref(gridData), eachrow(x))
2 Likes
Does this or any other package provide the interpolation for the other direction, that is, given x and fx and a rectangular grid, estimate the values in griddata?
Parallelization wont be as easy because it will require binning the x coordinates and then staging the writes into sub grids to avoid conflicting updates on the same griddata location.
1 Like
Packages like GMT.jl or DIVAnd.jl can grid irregular data…
1 Like
My package RadialBasisFunctions.jl can do this. The function is called regrid (see docs). You can map any set of points where you have the function values to any other set of points you wish, for any dimension as well. The input must be a vector of Tuple or SVector - something you can infer the length of. This is useful if your data is changing and you need to performthe mapping a lot. If you only need it once, it won’t be as fast as other methods, but it might be convenient to use.
using RadialBasisFunctions
using StaticArrays
using Statistics
using Random
f(x) = 1 + sin(4 * x[1]) + cos(3 * x[1]) + sin(2 * x[2])
N = 1_000
x = map(x -> SVector{2}(rand(MersenneTwister(x), 2)), 1:N)
y = f.(x)
x2 = map(x -> SVector{2}(rand(2)), 1:100)
r = regrid(x, x2)
f_at_x2 = r(y)
1 Like
Thank you. If you only need it once...: Do you mean if y is a short vector, then it won’t be as fast?
Yes, it’s not as fast as other packages if you need to query for different x2 everytime. You can do that, but you have to rebuild everytime. This function is optimal for when you have a known x and x2 and then you just want to get y2 where y is an input. It isn’t clear to me what your exact scenario is. What are your knowns and unknowns/queries?
For my case, my x2 is relatively large: 1042 points. So I think this package may be a good fit. Thanks!
It’s still a little rough around the edges. It came out of my graduate work. Now that I’ve finished and have more time, I’m going to be dedicating some time to it and tidying it up. Please make issues int he repo for any comments or questions, etc.
There is also GitHub - eljungsk/ScatteredInterpolation.jl: Interpolation of scattered data which works pretty well !
2 Likes
Thanks. I wonder which method you use in your package. Is it possible to incorporate cubic spline? I’m writing a code doing this on my own. But I’m not an expert in Julia and coding, so if someone could provide a high performance package with high speed that would be great! If cubic spline is not the method you want to incorporate, that is totally fine!
Given how things are implemented in ScatteredInterpolation, it should be fairly easy to add your own RBF kernel if not already available.
As for ScatteredInterpolation.jl, you can add your own kernels as well. You just need to make a new type which is subtype of AbstractRadialBasis. Is there a reason you need to do cubic splines specifically? ScatteredInterpolations.jl is a good choice.
Okay, I’ll take a look at it. Thank you. I want to use cubic spline because for the specific problems I’m solving, most previously related research uses cubic spline which proves to have a better performance.
Over RBFs with polynomials? In terms of accuracy, these packages perform better in general than vanilla cubic splines
I would get something working first, and only worry about finding better algorithms when you see where the problems and performance bottlenecks are.
(Also, the relative performance of different algorithms can be very implementation-dependent, so I wouldn’t assume that someone else’s experience with another language and another library will be identical to yours.)
1 Like
I also have KernelInterpolation.jl, which is similar to ScatteredInterpolations.jl and RadialBasisFunctions.jl.
1 Like
@JoshuaLampert, FYI, I tried your excellent package for interpolating the scattered points of the sin(x*y) function as in the example in this post, and got what appears to be excellent results for that type of input data:
KernelInterpolation code
using KernelInterpolation, Plots
# 1 - INPUT DATA:
n = 200
xs, ys = 2Ď€*(rand(n) .- 0.5), 2Ď€*(rand(n) .- 0.5)
zs = 100*sin.(xs .* ys)
# 2 - KERNEL INTERPOLATION:
nodes = NodeSet([xs ys])
values = zs
kernel = GaussKernel{dim(nodes)}(shape_parameter = 1.0)
itp = KernelInterpolation.interpolate(nodes, values, kernel)
N = 100
nodeset = homogeneous_hypercube(N, (-π, -π), (π, π))
itp_values = itp.(nodeset)
# 3 - PLOT RESULTS:
x = unique(values_along_dim(nodeset, 1))
y = unique(values_along_dim(nodeset, 2))
heatmap(x, y, reshape(itp_values, (N,N))', ratio=1, lims=(-Ď€,Ď€))
scatter!(xs, ys, marker_z = zs, ms=3, msw=0.5)
2 Likes