Oblique slices in Makie

Are there any way to plot an oblique slice in Makie? Something similar to Matlab’s obliqueslice function.

You can plot a heatmap and transform it such that it’s correctly positioned and angled. So the basic ingredient is there, it’s just not packaged up in a nice function so you’d have to do the transformation math yourself. And taking the slice seems a separate problem that Makie doesn’t have a function for, but that’s less related to plotting itself.

Or I guess surface makes more sense, then you don’t need transformation math either, just the correct points with values from a slicing function.

And taking the slice seems a separate problem that Makie doesn’t have a function for

This is what I actually need.
I need volumeslices function with oblique slices.

Yes but volumeslice uses heatmap because at orthogonal slices that’s fine, however for oblique you’d need surface I think, or even mesh. Can’t help you with the slice extraction though, maybe there are algorithms that are easy to port floating around somewhere, or already implemented in some Julia package.

Consider a volume, V, and a scalar field, f, defined on V. An oblique slice is defined as a plane through a point, p, in the given volume, and the normal to plane (a non-null vector).
From nullspace(normal) we extract the basis of vectors, a, b, parallel to the plane, and parameterize the plane.
On this slice is plotted a heatmap according to the values of the scalar field at the plane points.

f = (x,y,z)-> x*exp(-x^2-y^2-z^2) #scalar field
#slice defined by the plane through p and of of normal no
p  = [0, 1, 0.5]
no = [1 -0.3 1]
M = nullspace(no)
a, b = eachcol(M)
N=100
ul = vl= range(-3, 1.5, length=N)
u = ones(N, N) .* ul'
v = vl .* ones(N,N)
#slice parameterization
x = @. a[1]*u+b[1]*v+p[1]
y = @. a[2]*u+b[2]*v+p[2]
z = @. a[3]*u+b[3]*v+p[3]
surfcolor = f.(x,y,z)
vmin, vmax = extrema(surfcolor)
#with x, y, z, and surfcolor makie  can plot the slice
#---------Here is a plot with PlotlyJS--------------------#
pl =Plot(surface(x=x, y=y, z=z, surfacecolor=surfcolor, 
                 cmin=vmin, cmax=vmax, 
                 colorscale=colors.balance, 
                 colorbar_thickness=23, colorbar_len=0.7),
         Layout(width=500, height=500, font_size=11,
                scene_camera_eye=attr(x=1.8, y=1.8, z=1)))

oblique-slice-sf
Your link to a matlab example illustrates an oblique slice into a stack of RMN images. If you are interested in such a slice let us know to post an example, which differs from the above one, because the color at each slice point is calculated by interpolating the colors of the nearest two voxels on vertical direction.

2 Likes

Amazing! Thank you very much.

Your link to a matlab example illustrates an oblique slice into a stack of RMN images. If you are interested in such a slice let us know to post an example

This would be awesome. I do not have a set of images. Instead, I have three 1D arrays x,y,z which define the spatial coordinates and one 3D array (one scalar value for each of x,y,z coordinate) which I want to visualize. As a result, I also will need some sort of interpolation.

In this example is defined an oblique slice into a volume of brain images, where the slice is defined by a rotation with α degrees, of the plane z=a, about x-axis or y-axis.
a is recommended to lie between 0 and the volume height
But we can define similarly with the previous example, a slice through a point in volume and its normal.

using FileIO, Images, HTTP
using PlotlyJS
url = "https://s3.amazonaws.com/assets.datacamp.com/blog_assets/attention-mri.tif"
download(url, "attention-mri.tif")
img = load("attention-mri.tif");
vol = permutedims(channelview(img),  (1, 3, 2));
bone =  [[0.0, "rgb(255, 255, 255)"], #matplotlib colormap; I couldn't find the Colorschemes.bone
         [0.1, "rgb(220, 233, 233)"],
         [0.2, "rgb(185, 210, 210)"],
         [0.3, "rgb(156, 184, 188)"],
         [0.4, "rgb(133, 153, 165)"],
         [0.5, "rgb(112, 123, 143)"],
         [0.6, "rgb(89, 92, 121)"],
         [0.7, "rgb(66, 66, 92)"],
         [0.8, "rgb(44, 44, 62)"],
         [0.9, "rgb(21, 21, 30)"],
         [1.0, "rgb(0, 0, 0)"]];
function oblique_slice(vol::Array{T,3};  a=0, α=0,  axis="x", colorscale=bone) where T
    r, c, h = size(vol)
    xl = 0:1:c-1
    yl = 0:1:r-1
    x = ones(r, c) .* xl'
    y = yl .* ones(r,c)
    α = deg2rad(α)
    if axis == "x"
        z  = a/cos(α) .+ y*tan(α)
    elseif axis == "y" 
         z = a/cos(α) .-x*tan(α)  
    else
        error("rotate only about xaxis and yaxis")
    end
    sc = similar(z) #surfacecolor
    for idx in CartesianIndices(z) 
        s = floor(Int, z[idx])
        if 0 <= z[idx] <= h-1
            t = z[idx]-s 
            # the color at the slice point, P(x[i,j], y[i,j], z[i, j)), is the
            # interpolated color of the two nearest  voxels on the same vertical
            sc[idx] = (1-t)*vol[idx, s+1] + t*vol[idx, s+2]  
        else
            sc[idx] = NaN
        end   
    end 
    #hmax=maximum(z[(x->!isnan(x)).(sc)])
    return Plot(surface(x=x, y=y, z=z, surfacecolor=sc, 
                        colorscale=colorscale, #reversescale=true,
                        showscale=false), 
                Layout(width=500, height=500, font_size=11, 
                       scene_zaxis_range=(0, h) #or (0, hmax)
            ))
end 

fig1 = oblique_slice(vol, a=35, α=-10, axis="x")

brain-slice-x

fig2 = oblique_slice(vol, a=48, α=12, axis="y")

brain-slice-y

1 Like

Thank you. Looks like this is what I need.

This is amazing,

I try Python code to show this:
https://plotly.com/python/visualizing-mri-volume-slices/

but this one is better.

Can I ask can we create the interactive slider just like the Python, but with Julia, with the slice angle that has been chosen. So it will not be vertical up and down.

Here is the Julia version with slider: https://discourse.julialang.org/t/plotlyjs-jl-animation/70058/3

Another variation of the solution which allows the cuts of complex shape and not only the square planes. The idea is to calculate all 3D points which belong to the cut plane and then use meshscatter to plot the in-plane data pixel by pixel.

using GLMakie
using Interpolations
using LinearAlgebra


function distance_to_plane(P, point, normal)
    return dot(P - point, normal)
end


function isinplane(P, point, normal; rtol=1e-3)
    d = distance_to_plane(P, point, normal)
    return isapprox(1 + d, 1; rtol=rtol)   # more precise when compare to 1
end


# Data to plot:
xmin, xmax, Nx = -1.0, 1.0, 51
ymin, ymax, Ny = -1.0, 1.0, 51
zmin, zmax, Nz = -1.0, 1.0, 51
x = range(xmin, xmax, Nx)
y = range(ymin, ymax, Ny)
z = range(zmin, zmax, Nz)
F = zeros((Nx, Ny, Nz))
for iz=1:Nz, iy=1:Ny, ix=1:Nx
    F[ix,iy,iz] = exp(-(x[ix]^2 + y[iy]^2 + z[iz]^2))
end


# Cut plane is given by a point and normal:
point, normal = [0, 0, 0], [-1, 0, 0]

# unit vectors:
xuv, yuv, zuv = Point3f(1,0,0), Point3f(0,1,0), Point3f(0,0,1)

# transform x,y,z coordinates to points in 3D:
points = [ix * xuv + iy * yuv + iz * zuv for ix in x for iy in y for iz in z]

# find all points in the plane:
inpoints = [pt for pt in points if isinplane(pt, point, normal)]

# find values of F in all in-plane points:
itp = linear_interpolation((x, y, z), F; extrapolation_bc=Flat())
fcolors = [itp.(pt[1], pt[2], pt[3]) for pt in inpoints]

# cut F by the plane:
for iz=1:Nz, iy=1:Ny, ix=1:Nx
    P = x[ix] * xuv + y[iy] * yuv + z[iz] * zuv
    if distance_to_plane(P, point, normal) > 0
        F[ix,iy,iz] = NaN
    end
end


# Plot:
fig = Figure(resolution=(800,800), fontsize=12)
title = "point=$(point),   normal=$(normal)"
ax = Axis3(fig[1,1]; aspect=(1,1,1), perspectiveness=0, title=title)

cmap = Reverse(:Hiroshige)

mask = @. !isnan(F)
Fmin, Fmax, NF = 0.1*maximum(F[mask]), 0.9*maximum(F[mask]), 10
levels = [collect(range(Fmin, Fmax, length=NF))...]
contour!(ax, x, y, z, F; colormap=cmap, alpha=0.2, levels=levels)

# marker = Sphere(Point3f(0), 1)   # spherical marker
marker = GLMakie.Rect3D(Vec3f(-0.5), Vec3f(1))   # cubic marker
meshscatter!(
    ax, inpoints;
    color=fcolors, colormap=cmap, markersize=0.04, shading=false, transparency=true,
)

display(fig)

Though, there are still a number of issues:

  • There is a possibility that the cut plane will go through the region where there are not enough points belonging to the plane (e.g. point, normal = [0.1, 0, 0], [-1, 0, 0]). We need some interpolation routine to find the proper points, or to force the plane to stick to the nearest available cut.
  • So far, as the pixel (marker) I use a cube (it looks better than sphere at right-angle cuts). Ideally we should use some thin plate oriented parallel the cut plane.
  • I did not find the way how to control the transparency of the cut plane.
2 Likes

@fedoroff I have no working experience with Makie, so I cannot suggest how to control slice transparency. :frowning:

I made a feature request in Makie a while ago for a obliqueslice: Feature request: oblique volume slice · Issue #2119 · MakieOrg/Makie.jl · GitHub

I have a MWE in Julia, but due to my poor knowledge of Makie internals, I don’t know how to implement this directly in Makie…

Yes, I still miss this feature in Makie.