Visualize 3D polygon of brain regions given points representing exterior border

vandenman · July 6, 2023, 9:15am

I’m trying to create a 3d plot of some brain regions. However, the data I have is not a mesh but a series of 3d points that indicate the exterior border of the regions. How can I plot this with GLMakie?

Below is the code I’m using. The two data files, ROI_MNI_V4_List.mat and ROI_MNI_V4_Border.mat can be obtained from EBRAINS - Knowledge Graph Search and then going to Get data.

using GLMakie
import Colors, GeometryBasics
import MAT # for reading matlab data files
# [5ae59095] Colors v0.12.10
# [e9467ef8] GLMakie v0.8.6
# [5c1252a2] GeometryBasics v0.4.7
# [23992714] MAT v0.10.5

# wherever you downloaded the aal atlas
root = joinpath(pwd(), "aal_for_SPM12")
roi_list   = MAT.matread(joinpath(root, "ROI_MNI_V4_List.mat"))
roi_border = MAT.matread(joinpath(root, "ROI_MNI_V4_Border.mat"))

distinct_ids = unique(vec(Int.(roi_list["ROI"]["ID"])))
ids = eachindex(distinct_ids)
ids_indices = [findall(==(distinct_ids[id]), vec(roi_border["BORDER_V"])) for id in ids]

# a vector with the boundary regions for each id
border_points = [
    GeometryBasics.Point3f.(eachcol(view(roi_border["BORDER_XYZ"], :, indices)))
    for indices in ids_indices
]

mesh(border_points[1]) # visualize only the first region for now

The code above sort of plots something okay, but the mesh is not a volume but a disjoint set of triangles.

I tried a bunch of things, but I also noticed that e.g., GeometryBasics.mesh fails completely for this type of data and returns an empty Mesh:

julia> GeometryBasics.mesh(border_points[1])
Mesh{3, Float32, Triangle}

I figured that the border points essentially represent a 3d polygon, but that doesn’t help much either:

julia> polygon = GeometryBasics.Polygon(border_points[1])
julia> GeometryBasics.mesh(polygon)
ERROR: MethodError: no method matching earcut_triangulate(::Vector{Vector{Point{3, Float32}}})

Closest candidates are:
  earcut_triangulate(::Vector{Vector{Point2{Float64}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/triangulation.jl:181
  earcut_triangulate(::Vector{Vector{Point{2, Float32}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/triangulation.jl:189
  earcut_triangulate(::Vector{Vector{Point2{Int64}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/triangulation.jl:197
  ...

Stacktrace:
 [1] faces(polygon::GeometryBasics.Polygon{3, Float32, Point{3, Float32}, GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}, Vector{GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/triangulation.jl:224
 [2] decompose(#unused#::Type{GeometryBasics.NgonFace{3, GeometryBasics.OffsetInteger{-1, UInt32}}}, primitive::GeometryBasics.Polygon{3, Float32, Point{3, Float32}, GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}, Vector{GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/interfaces.jl:103
 [3] decompose_triangulate_fallback(primitive::GeometryBasics.Polygon{3, Float32, Point{3, Float32}, GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}, Vector{GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}}}; pointtype::Type, facetype::Type)
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/meshes.jl:91
 [4] decompose_triangulate_fallback
   @ ~/.julia/packages/GeometryBasics/Du43H/src/meshes.jl:89 [inlined]
 [5] mesh(polygon::GeometryBasics.Polygon{3, Float32, Point{3, Float32}, GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}, Vector{GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}}}; pointtype::Type, facetype::Type, normaltype::Nothing)
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/meshes.jl:161
 [6] mesh(polygon::GeometryBasics.Polygon{3, Float32, Point{3, Float32}, GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}, Vector{GeometryBasics.LineString{3, Float32, Point{3, Float32}, Base.ReinterpretArray{GeometryBasics.Line{3, Float32}, 1, Tuple{Point{3, Float32}, Point{3, Float32}}, GeometryBasics.TupleView{Tuple{Point{3, Float32}, Point{3, Float32}}, 2, 1, Vector{Point{3, Float32}}}, false}}}})
   @ GeometryBasics ~/.julia/packages/GeometryBasics/Du43H/src/meshes.jl:158
 [7] top-level scope

For what it’s worth, if I visualize all regions then this does produce some rather artsy figures:

but as soon as you look at a different angle it’s clear that everything is messed up

Does anybody have any ideas on how to fix this?

johnh · July 6, 2023, 10:10am

A quick search returns this

sdanisch · July 6, 2023, 10:38am

I don’t think surface in plotly will do this, since it’s not a surface which is described by the points.
The simplest I could find is to compute the convex hull (which is probably not optimal, not sure what other algorithms there are to compute the triangles out of an unordered pointcloud).
This works:

function compute_hull(points3d::AbstractMatrix)
    hull = chull(collect(points3d))
    tris = [GLTriangleFace(hull.simplices[i, :]...) for i in 1:size(hull.simplices, 1)]
    return normal_mesh(Makie.to_vertices(points3d), tris)
end
# a vector with the boundary regions for each id
border_meshes = [
    compute_hull(view(roi_border["BORDER_XYZ"], :, indices)')
    for indices in ids_indices
]

mesh(border_meshes; color=1:length(border_meshes), colormap=(:viridis, 0.5), transparency=true) # visualize only the first region for now

vandenman · July 6, 2023, 11:25am

Thanks that works perfectly!

To save anybody trying to replicate this some time, the function chull comes from the library QHull.