Verified the existence with method I described in Grassmann.jl A\b 3x faster than Julia's StaticArrays.jl
First, you need to convert the mesh into Grassmann format:
using Grassmann
function mesh_init(m::GeometryBasics.Mesh)
c,f = coordinates(m),faces(m)
s = size(eltype(f))[1]; V = SubManifold(ℝ^s)
p = ChainBundle([Chain{V,1}(1.0,k...) for k ∈ c])
t = ChainBundle([Chain{p,1}(SVector(k)) for k ∈ f])
return (p,t)
end
Then define the following checkpoint for determining the tetrahedron containing a point:
function checkpoint(P::Chain{V,1},t) where V
p = Manifold(t)
for i ∈ 1:length(t)
P ∈ Chain{V,1}(p[value(t[i])]) && (return i)
end
return 0
end
Then you can determine whether there exists a tetrahedron with point P or not (returns 0 if not).
julia> p,t = mesh_init(result) # result from TetGen
(Λ¹⟨++++⟩×211, Λ¹Λ¹⟨++++⟩×211×566)
julia> checkpoint(Chain{SubManifold(ℝ^4),1}(1,2,2,2),t)
122
Note that I am using homogeneous coordinates with (1,x,y,z) coordinates for the points, so you will have to use 4 dimensional points with an additional 1 coordinate.
You will find that indeed the mesh data contains simplices having points in the problematic regions.
PS, it would be great if Makie could natively support Grassmann so I don’t need GeometryBasics as a compatibility layer. Instead of GeometryBasics.Point I prefer using Grassmann.Chain{V,1}, because I need that type in order to do my calculations. As you can see, I am using homogeneous coordinates with Chain{V,1} types in Grassmann to represent mesh data. Currently, I need to convert this mesh data back into GeometryBasics.Point for plotting (using cache), which works fine, but I would much prefer being able to work with Grassmann.Chain{V,1} directly without conversions.
GLMakie.mesh(t, color=(:blue, 0.1), transparency=true)
wireframe!(ans[end][1])
As you can see, it’s possible to use mesh function on Grassmann data as well, but it requires the conversion process with caching I mentioned above. I’d prefer to have Makie work with Grassmann point data directly instead of converting to GeometryTypes.Point.