Stream: helpdesk (published)

Topic: Matlab interp3/slice


view this post on Zulip Dale Black (Apr 02 2022 at 21:40):

Does anyone know a seamless way to reimplement this interp3() functionality in julia, specifically with the syntax as the first example below?

Vq = interp3(X,Y,Z,V,Xq,Yq,Zq)

view this post on Zulip Dale Black (Apr 03 2022 at 21:41):

Or if this slice() functionality has a julia equivalent, that would be better/easier slice(X,Y,Z,V,xslice,yslice,zslice)

view this post on Zulip Benoit Pasquier (Apr 04 2022 at 03:45):

Doesn't Interpolations.jl do the trick?

view this post on Zulip Dale Black (Apr 04 2022 at 21:07):

Ahh okay, it does look like that's what I am after. I am just now realizing that it's not the missing functionality but that I don't quite understand my problem as well as I should. Are there any examples of multi-planar reformation in Julia that you know of because I haven't found any from the simple google searches and that's what I need

view this post on Zulip Benoit Pasquier (Apr 05 2022 at 02:12):

Not sure what that is but if I understand it's just an interpolation? If you have your X ,Y, Z, V, Xq, Yq, and Zq you should be able to just make the interpolation with the Interpolations.jl package, right? Maybe make a small example of these and try to create the slices you're looking for?

view this post on Zulip Dale Black (Apr 05 2022 at 16:03):

I guess I am confused about how Interpolations.jl takes in all of the arguments X, Y, Z, V, Xq, Yq, Zq? For example

itp = interpolate(Xs, Ys, Zs, V, BSpline(Linear()))
itp(Xq, Yq, Zq)

Doesn't work and I can't figure out what order of arguments I should be using if X Y``Z are the points and Xq, Yq, Zq are the query points

view this post on Zulip Dale Black (Apr 05 2022 at 16:04):

This is the error I get if this is helpful btw

MethodError: no method matching interpolate(::Array{Float64, 3}, ::Array{Float64, 3}, ::Array{Float64, 3}, ::Array{Int16, 3}, ::Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}})

Closest candidates are:

interpolate(::AbstractArray, !Matched::Interpolations.NoInterp) at ~/.julia/packages/Interpolations/Glp9h/src/nointerp/nointerp.jl:1

interpolate(!Matched::Type{TWeights}, !Matched::Type{TC}, ::Any, !Matched::IT) where {TWeights, TC, IT<:Union{Interpolations.NoInterp, Tuple{Vararg{Union{Interpolations.NoInterp, Interpolations.BSpline}}}, Interpolations.BSpline}} at ~/.julia/packages/Interpolations/Glp9h/src/b-splines/b-splines.jl:159

interpolate(!Matched::Type{TWeights}, !Matched::Type{TC}, ::Any, !Matched::IT, !Matched::Real, !Matched::Int64) where {TWeights, TC, IT<:Union{Interpolations.NoInterp, Tuple{Vararg{Union{Interpolations.NoInterp, Interpolations.BSpline}}}, Interpolations.BSpline}} at ~/.julia/packages/Interpolations/Glp9h/src/b-splines/b-splines.jl:164

...

PlanarReformation(::Array{Int16, 3}, ::Vector{Float64}, ::Matrix{Float64}, ::Matrix{Float64}, ::Float64, ::Float64)@Other: 15
top-level scope@Local: 1[inlined]

view this post on Zulip Benoit Pasquier (Apr 06 2022 at 00:40):

@Mark Kittisopikul maybe the docs need a 2D and 3D examples?

view this post on Zulip Benoit Pasquier (Apr 06 2022 at 00:41):

Otherwise, does this work for you?

f(x, y, z) = sin(x + 1) * (1 + cos(y)) * exp(-z^2)
xs = range(0, 2π, length=11)
ys = range(-π/2, π/2, length=6)
zs = range(-1, 1, length=15)
V = [f(x, y, z) for x in xs, y in ys, z in zs]
itp = interpolate(V, BSpline(Linear())) # interpolate linearly between the data points
stp = scale(itp, xs, ys, zs) # re-scale to the actual domain
[itp(x, y, z) for x in Xq, y in Yq, z in Zq] # <- make your own Xq, Yq, Zq here

view this post on Zulip Dale Black (Apr 18 2022 at 18:59):

Sorry for such a late response, I had some important deadlines. This example is super helpful and should be enough to help me get started I think!

view this post on Zulip Dale Black (Apr 19 2022 at 00:12):

Okay, so in the MatLab function, I can use 3 different 3D arrays as the sample points (Xs Ys and Zs) along with a 3D array V. Interpolations.jl doesn't like this

Xs = rand(512, 512, 40)
Ys = rand(512, 512, 40)
Zs = rand(512, 512, 40)

V = rand(512, 512, 40)

So if I have those arrays, is there a way to pass them into the interpolate function that makes sense? In matlab its just interp3(Xs, Ys, Zs, V, ...) which is what I am hoping to recreate in Julia

view this post on Zulip Dale Black (Apr 19 2022 at 00:25):

It looks like the interp3 in MatLab takes the first 4 arguments Xs, Ys, Zs and V and then uses griddedInterpolant to create a new interpolant F so I will look into that step first using Interpolations.jl

view this post on Zulip Benoit Pasquier (Apr 19 2022 at 05:56):

Are you using a regular grid? This surely error in MATLAB too:

>> Xs = rand(3, 4, 5);
>> Ys = rand(3, 4, 5);
>> Zs = rand(3, 4, 5);
>> Vs = rand(3, 4, 5);
>> Xq = rand(2, 3, 4);
>> Yq = rand(2, 3, 4);
>> Zq = rand(2, 3, 4);
>> Vq = interp3(Xs, Ys, Zs, Vs, Xq, Yq, Zq);
Error using griddedInterpolant
Grid arrays must have NDGRID structure.

Error in interp3 (line 144)
            F = griddedInterpolant(X, Y, Z, V, method,extrap);

view this post on Zulip Benoit Pasquier (Apr 19 2022 at 05:57):

If your grid is regular, then no need to have full 3D arrays for x, y, and y: A vector suffices for each.

view this post on Zulip Benoit Pasquier (Apr 19 2022 at 05:59):

Can you make a copy-pastable minimal working example (MWE) of MATLAB code that does not error? And then we can see if we can do the same with Interpolations.jl :shrug:

view this post on Zulip Dale Black (Apr 19 2022 at 16:19):

Here is a MWE in MatLab that I would like to reimplement in Julia

>> [Xs, Ys, Zs] = meshgrid(1:512, 1:512, 1:56);
>> V = rand(512, 512, 56);
>> plane_list = rand(4225, 3);
>> plane = interp3(Xs, Ys, Zs, V, plane_list(:,1), plane_list(:,2), plane_list(:,3));

% Result: size(plane) = 4225 x 1

view this post on Zulip Dale Black (Apr 19 2022 at 16:21):

But it seems like the problem is coming from my misunderstanding of the regular grid usage in Julia?

view this post on Zulip Dale Black (Apr 19 2022 at 16:30):

What do you mean a vector will suffice for 3D regular grids? Is it something like this?

V = rand(512, 512, 40)
x = 1:size(V, 1)
y = 1:size(V, 2)
z = 1:size(V, 3)
itp = interpolate(x, y, z, V, Gridded(Linear()))

If so, that still isn't working (if I am even doing it right)

view this post on Zulip Dale Black (Apr 19 2022 at 16:31):

Nor is this

V = rand(512, 512, 40)
x = collect(1:size(V, 1))
y = collect(1:size(V, 2))
z = collect(1:size(V, 3))
itp = interpolate(x, y, z, V, Gridded(Linear()))

view this post on Zulip Nils (Apr 20 2022 at 10:37):

Have you looked at the Interpolations.jl docs? Specifically those on gridded interpolation here? Quoting from there:

The general syntax is

itp = interpolate(nodes, A, options...)

where nodes = (xnodes, ynodes, ...) specifies the positions along each axis at which the array A is sampled for arbitrary ("rectangular") samplings.

view this post on Zulip Nils (Apr 20 2022 at 10:38):

In your case:

julia> itp = interpolate((x, y, z), V, Gridded(Linear()));

julia> itp(2.3, 3.5, 11.7)
0.5487285843274945

view this post on Zulip Dale Black (Apr 20 2022 at 12:00):

Ahhhhhh how did I miss that. I will give it a shot later today but that looks like my answer haha. Thanks!


Last updated: Nov 06 2024 at 04:40 UTC