# Interpolate a 3D array
**URL:** https://discourse.julialang.org/t/interpolate-a-3d-array/93885
**Category:** Signal and Image Processing
**Tags:** interpolations, meshes
**Created:** [February 1, 2023, 8:27pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885 "2023-02-01T20:27:14Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 1, 2023, 8:27pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/1 "2023-02-01T20:27:15Z")
If I have a 3D array along with a unit normal vector and point for a plane that intersects that array, how would I use `Interpolations.jl` to create a new 3D array that is the same size as the original array but where all the voxels are interpolated to be perpendicular to that plane? So far I have an interpolation but I am not sure how to use it:
```julia
# Define original array
A = rand(512, 512, 320)
# Define the normal vector of the plane
normal = [0.980998, -0.0695175, -0.181134]
# Define a point on the plane
point = [190, 195, 7]
A_x, A_y, A_z = 1:512, 1:512, 1:320
nodes = (A_x, A_y, A_z)
itp = interpolate(nodes, A, Gridded(Linear()))
```
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 1, 2023, 9:20pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/2 "2023-02-01T21:20:03Z")
> [@Dale\_James\_Black](#):
>
> a plane that intersects that array,
What does it mean for a plane to intersect the array? Do you mean that each entry `A[i,j,k]` in the array corresponds to a point `(x[i],y[j],z[k])` on a grid for a box in space that intersects the plane?
> [@Dale\_James\_Black](#):
>
> where all the voxels are interpolated to be perpendicular to that plane
What does it mean for a voxel to be interpolated perpendicular to a plane? The entries `A[i,j,k]` in your example are scalars.
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 1, 2023, 9:31pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/3 "2023-02-01T21:31:47Z")
Good questions. I guess I am thinking more along the lines of multi-planar reformation in medical imaging [link](https://www.google.com/url?sa=i&url=https%3A%2F%2Fwww.radiantviewer.com%2Fdicom-viewer-manual%2F3d-multiplanar-reconstructions.html&psig=AOvVaw0jTJv4FeMzrXbnjQwEfdzZ&ust=1675373194433000&source=images&cd=vfe&ved=0CA4QjRxqFwoTCOCbqNKh9fwCFQAAAAAdAAAAABAI)
> Blockquote
> What does it mean for a plane to intersect the array? Do you mean that each entry `A[i,j,k]` in the array corresponds to a point `(x[i],y[j],z[k])` on a grid for a box in space that intersects the plane?
I don’t think each point in the array corresponds to a point on the plane. I am thinking about some plane at an arbitrary angle that “slices” through the array. I am then trying to use this plane to reformat the original array. Hopefully that makes it a little more clear, but it’s still not super clear in my own head…
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 1, 2023, 9:42pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/4 "2023-02-01T21:42:01Z")
> [@Dale\_James\_Black](#):
>
> I don’t think each point in the array corresponds to a point on the plane. I am thinking about some plane at an arbitrary angle that “slices” through the array.
Yes, that’s what I said — the array elements corresponds to points on a Cartesian grid and the plane slicing through that.
> [@Dale\_James\_Black](#):
>
> I am then trying to use this plane to reformat the original array.
That’s the part that is unclear to me. What are the new values that you want?
(My first thought is that you wanted the values in the array interpolated onto a grid in the plane, i.e. for visualizing an arbitrary slice of a 3d array. This seems to be what is happening in the “multiplanar reconstruction” that you linked. But you said that you want a 3d array out, not a 2d array, so I’m still confused. Perhaps you want a 3d array consisting of a sequence of slices parallel to the given plane?)
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 1, 2023, 10:01pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/5 "2023-02-01T22:01:46Z")
> Blockquote
> That’s the part that is unclear to me. What are the new values that you want?
> (My first thought is that you wanted the values in the array interpolated onto a grid in the plane, i.e. for visualizing an arbitrary slice of a 3d array. This seems to be what is happening in the “multiplanar reconstruction” that you linked. But you said that you want a 3d array out, not a 2d array, so I’m still confused. Perhaps you want a 3d array consisting of a sequence of slices parallel to the given plane?)
Ahh, yes that is a better way to describe it. I would like a 3D array of slices parallel to the given plane. So yes, a sequence of 2D slices is more accurate
---
### Author:  [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)
#### Post date: [February 1, 2023, 10:25pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/6 "2023-02-01T22:25:42Z")
```julia
julia> using Interpolations, Meshes
julia> A_itp = linear_interpolation(axes(A), A);
julia> plane = Plane(Point(point), Vec(normal))
Plane{Float64}(Point(190.0, 195.0, 7.0), [0.06951752486064443, 0.9975604794114311, -0.006356386842015778], [0.18113406477660973, -0.006356386842015778, 0.9834378714101962])
julia> locations = [plane(i, j) + Vec(normal)*k for i in -5:5, j in -5:5, k in -5:5]
11×11×11 Array{Point3, 3}:
[:, :, 1] =
Point(183.84175205181373, 190.39156703715292, 3.020262577159097) … Point(185.65309269957982, 190.32800316873278, 12.85464129126106)
...
julia> [A_itp(SVector(loc.coords...)...) for loc in locations]
11×11×11 Array{Float64, 3}:
[:, :, 1] =
0.682766 0.455143 0.344485 0.499102 0.495368 0.283909 0.714595 0.153907 0.417154 0.502021 0.602994
0.699649 0.775464 0.496942 0.452601 0.499833 0.184824 0.512752 0.422986 0.268267 0.387782 0.29775
0.307134 0.611958 0.27292 0.666716 0.38154 0.426064 0.412537 0.697101 0.292527 0.506895 0.500114
...
```
---
### Author:  [@pxshen](https://discourse.julialang.org/u/pxshen)
#### Post date: [February 1, 2023, 10:28pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/7 "2023-02-01T22:28:20Z")
Try `Images` lib: `imrotate, imresize`
---
### Author:  [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)
#### Post date: [February 1, 2023, 11:01pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/8 "2023-02-01T23:01:28Z")
Linking [related thread](https://discourse.julialang.org/t/oblique-slices-in-makie/83879).
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 1, 2023, 11:01pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/9 "2023-02-01T23:01:54Z")
Thank you, that looks like it’s exactly what I am looking for. I will investigate a bit more but for now I will go with this!
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 1, 2023, 11:51pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/10 "2023-02-01T23:51:07Z")
Do I need to involve an `extrapolate` function of some sort? One problem I am running into is bad indexing when I make the `locations` range larger than `-5:5`. Or is this just because of the linear interpolation?
e.g.
```julia
r = 6
locations = [plane(i, j) + Meshes.Vec(unit_normal)*k for i in -r:r, j in -r:r, k in -r:r]
test_arr = [A_itp(loc.coords...) for loc in locations]
BoundsError: attempt to access 512×512×320 extrapolate(scale(interpolate(::Array{Float64, 3}, BSpline(Linear())), (Base.OneTo(512), Base.OneTo(512), Base.OneTo(320))), Throw()) with element type Float64 at index [155.6834977746412, 198.78313469043155, -0.9594702586293882]
```
---
### Author:  [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)
#### Post date: [February 2, 2023, 12:03am UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/11 "2023-02-02T00:03:57Z")
You need to define a boundary condition (see [here](http://juliamath.github.io/Interpolations.jl/latest/extrapolation/)) - e.g. `linear_interpolation(...; extrapolation_bc = 0.0)`.
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 2, 2023, 12:04am UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/12 "2023-02-02T00:04:51Z")
~~Oh gotcha. It also seems like using `extrapolate()` solves this issue. Is there a preferred approach?~~
`linear_interpolation(; extrapolation_bc=0)` is much faster so that makes sense to use
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 12:21am UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/13 "2023-02-03T00:21:25Z")
Okay, everything works based on your suggestion! One question though - how would you make it so that the interpolation places the object of interest in the center? If you can point me in the right direction, I will publish a simple MultiPlanarReformation.jl package based on your implementation, since I think that would be very useful to have in this ecosystem.
For example, this image below is the 3D array that I have and the plane that I used for interpolation comes from the line of air through the middle
When I perform interpolation based on your code, the reformatted array looks exactly like I would expect, just not centered:
Here is the code snippet:
```julia
plane = Plane(Meshes.Point(pts[1]), Meshes.Point(pts[2]), Meshes.Point(pts[3]))
A_itp = linear_interpolation(axes(dcm_heart), dcm_heart; extrapolation_bc = 0.0);
r, r2 = 512/2, 320/2
locations = [plane(i, j) + normal(plane)*k for i in -r:r, j in -r:r, k in -r2:r2]
mpr = [A_itp(loc.coords...) for loc in locations];
```
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 12:32am UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/14 "2023-02-03T00:32:34Z")
Does the plane need to be initialized in a way that matches the origin of the plane to the center point of the original 3D array (or something along those lines)?
---
### Author:  [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)
#### Post date: [February 3, 2023, 12:54am UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/15 "2023-02-03T00:54:06Z")
Can you share a full reproducer? `pts` and `dcm_heart` don’t appear in your original example.
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 5:43pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/16 "2023-02-03T17:43:24Z")
The `dcm_heart` array is kinda hard to reproduce since it come’s from a 3D DICOM image. But it’s `size=(512, 512, 320)` with a cylindrical shape near the center, and zeros everywhere else.
`pts` is simpler:
```julia
pts = (
[203, 224, 7],
[197, 246, 123],
[229, 282, 314],
)
```
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 5:52pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/17 "2023-02-03T17:52:13Z")
Here is an example of an array that comes close to the `dcm_heart` array
```julia
function create_circle_mask(img, centroids, radius)
# initialize mask with all zeros
mask = zeros(size(img))
# define the center of the circle
x0, y0 = centroids[1], centroids[2]
# set all pixels inside the circle to 1 and all pixels outside the circle to 0
for x in axes(img, 1), y in axes(img, 2)
if ((x - x0)^2 + (y - y0)^2) <= radius^2
mask[x, y] = 1
else
mask[x, y] = 0
end
end
return Bool.(mask)
end
arr = rand(512, 512, 320)
centroids = [217, 254, 161]
mask = create_circle_mask(arr[:, :, 161], centroids, 75)
cylinder_arr = arr .* mask
```
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 6:04pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/18 "2023-02-03T18:04:31Z")
Sorry this is a messy thread now, but here is the full reproducer
```julia
function create_circle_mask(img, centroids, radius)
# initialize mask with all zeros
mask = zeros(size(img))
# define the center of the circle
x0, y0 = centroids[1], centroids[2]
# set all pixels inside the circle to 1 and all pixels outside the circle to 0
for x in axes(img, 1), y in axes(img, 2)
if ((x - x0)^2 + (y - y0)^2) <= radius^2
mask[x, y] = 1
else
mask[x, y] = 0
end
end
return Bool.(mask)
end
arr = rand(512, 512, 320)
mask = create_circle_mask(arr[:, :, 162], centroids, 75)
cylinder_arr = arr .* mask
points = (
[203, 224, 7],
[197, 246, 123],
[229, 282, 314],
)
plane = Plane(Meshes.Point(points[1]), Meshes.Point(points[2]), Meshes.Point(points[3]))
mpr_itp = linear_interpolation(axes(cylinder_arr), cylinder_arr; extrapolation_bc = 0.0)
r, r2 = 512/2, 320/2
locations = [plane(i, j) + normal(plane)*k for i in -r:r, j in -r:r, k in -r2:r2]
mpr2 = [mpr_itp(loc.coords...) for loc in locations]
```
---
### Author:  [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)
#### Post date: [February 3, 2023, 6:27pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/19 "2023-02-03T18:27:11Z")
`Meshes.plane` uses the first point it’s given as the origin of the resulting plane. If you want a plane centered around a different point, you can apply an offset to each of your points. It’ll be faster to compute the locations on the fly:
```julia
julia> offset = Point(centroids) - Point(first(pts))
3-element Vec3 with indices SOneTo(3):
14.0
30.0
154.0
julia> plane = Plane((Point.(pts) .+ Ref(offset))...)
Plane{Float64}(Point(217.0, 254.0, 161.0), [0.0004935295610762698, 0.18606852260760806, 0.9825366463003782], [0.9999860524254478, -0.0052584582476063145, 0.0004935295610762698])
julia> offset_interpolate(itp, plane, i, j, k) = itp((plane(i, j) + Meshes.normal(plane)*k).coords...)
offset_interpolate (generic function with 1 method)
mpr2 = [offset_interpolate(mpr_itp, plane, i, j, k) for i in -r:r, j in -r:r, k in -r2:r2]
```
---
### Author:  [@Dale\_James\_Black](https://discourse.julialang.org/u/Dale_James_Black)
#### Post date: [February 3, 2023, 6:58pm UTC](https://discourse.julialang.org/t/interpolate-a-3d-array/93885/20 "2023-02-03T18:58:04Z")
WOW, thank you!