# 3D interpolation in Julia

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