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]
