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

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]

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

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

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.

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

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

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]