# Interpolate a 3D array

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:

``````# 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()))
``````

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?

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.

Good questions. I guess I am thinking more along the lines of multi-planar reformation in medical imaging link

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â€¦

Yes, thatâ€™s what I said â€” the array elements corresponds to points on a Cartesian grid and the plane slicing through that.

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?)

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

``````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
...
``````
2 Likes

Try `Images` lib: `imrotate, imresize`

1 Like

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!

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.

``````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]
``````

You need to define a boundary condition (see here) - e.g. `linear_interpolation(...; extrapolation_bc = 0.0)`.

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

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:

``````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];
``````

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)?

Can you share a full reproducer? `pts` and `dcm_heart` donâ€™t appear in your original example.

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:

``````pts = (
[203, 224, 7],
[197, 246, 123],
[229, 282, 314],
)
``````

Here is an example of an array that comes close to the `dcm_heart` array

``````function create_circle_mask(img, centroids, radius)
# initialize mask with all zeros

# 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
else
end
end
end

arr = rand(512, 512, 320)
centroids = [217, 254, 161]

``````

Sorry this is a messy thread now, but here is the full reproducer

``````function create_circle_mask(img, centroids, radius)
# initialize mask with all zeros

# 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
else
end
end
end

arr = rand(512, 512, 320)

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]
``````

`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> 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]
``````

1 Like

WOW, thank you!