JuliaGeometry / GeometryBasics.jl

Basic Geometry Types
MIT License
165 stars 54 forks source link

GeometryBasics.Mesh incorrectly converts quads to triangles #207

Open Kevin-Mattheus-Moerman opened 8 months ago

Kevin-Mattheus-Moerman commented 8 months ago

When trying to convert quadrilateral faces F and their vertices V to a mesh using:

M=GeometryBasics.Mesh(V,F)

One correctly obtains a quadrilateral mesh description:

julia> typeof(GeometryBasics.Mesh(V,F))
GeometryBasics.Mesh{3, Float64, GeometryBasics.Ngon{3, Float64, 4, Point3{Float64}}, SimpleFaceView{3, Float64, 4, Int64, Point3{Float64}, QuadFace{Int64}}}

However an error occurs if the mesh only consists of a single quad face, i.e. length(F)==1, in this case it looks like that quad is instead interpreter as a triangle:

julia> typeof(GeometryBasics.Mesh(V,F[1]))
GeometryBasics.Mesh{3, Float64, GeometryBasics.Ngon{3, Float64, 3, Point3{Float64}}, FaceView{GeometryBasics.Ngon{3, Float64, 3, Point3{Float64}}, Point3{Float64}, TriangleFace{Int64}, Vector{Point3{Float64}}, Base.ReinterpretArray{TriangleFace{Int64}, 1, Tuple{Int64, Int64, Int64}, TupleView{Tuple{Int64, Int64, Int64}, 3, 1, QuadFace{Int64}}, false}}}

A longer example and some context is given below.

Here is an example quadrilateral mesh description, in this case for a cube.

    # Create vertices       
    V=Vector{GeometryBasics.Point{3, Float64}}(undef,8)
    s = r/sqrt(3.0)
    V[1 ] = GeometryBasics.Point{3, Float64}( -s, -s, -s)
    V[2 ] = GeometryBasics.Point{3, Float64}( -s,  s, -s)
    V[3 ] = GeometryBasics.Point{3, Float64}(  s,  s, -s)
    V[4 ] = GeometryBasics.Point{3, Float64}(  s, -s, -s)
    V[5 ] = GeometryBasics.Point{3, Float64}( -s, -s,  s)
    V[6 ] = GeometryBasics.Point{3, Float64}( -s,  s,  s)
    V[7 ] = GeometryBasics.Point{3, Float64}(  s,  s,  s)
    V[8 ] = GeometryBasics.Point{3, Float64}(  s, -s,  s)

    # Create faces
    F = Vector{QuadFace{Int64}}(undef,6)
    F[1 ] = QuadFace{Int64}(1,2,3,4)
    F[2 ] = QuadFace{Int64}(8,7,6,5)
    F[3 ] = QuadFace{Int64}(5,6,2,1)
    F[4 ] = QuadFace{Int64}(6,7,3,2)    
    F[5 ] = QuadFace{Int64}(7,8,4,3)    
    F[6 ] = QuadFace{Int64}(8,5,1,4)   

I am interested in visualizing such quadrilateral meshes using GLMakie and to use per face color data. Currently GLMakie requires looping over the faces to handle such "per face colors" (default is to use vertex color data instead).

using GLMakie

fig = Figure()

ax=Axis3(fig[1, 1], aspect = :data, xlabel = "X", ylabel = "Y", zlabel = "Z", title = "Hexahedral mesh")

for i ∈ eachindex(F)    
   poly!(ax,GeometryBasics.Mesh(V,F[i]),strokewidth=5,color=C[i], transparency=false, overdraw=false,colormap = Reverse(:Spectral), colorrange=(0,6))
end

fig

The looping causes me to have to pass single faces to Geometry.Mesh which was showing the error. Note the besides being converted to triangles, which is unwanted, the quad is also not correctly converted to triangles as one can tell from this crossed/entangled look: image

This is related to a bug I've posted on GLMakie, but I think it starts with this bug here.

Thanks a lot for your help!

SimonDanisch commented 7 months ago

GeometryTypes.Mesh(V, F) should be constructed with 2 arrays. I'm not sure why it doesn't error when passing it a single quad instead of an array of quads, so I guess that's what we should fix.

You can fix this example by using Mesh(V, F[i:i]) to create a 1 element array.

But calling Poly like this is not recommended, since each poly call has a large overhead. So I'd recommend this:

V=Vector{Point3f}(undef,8)
r = 1.0
s = r/sqrt(3.0)
V[1] = Point3f( -s, -s, -s)
V[2] = Point3f( -s,  s, -s)
V[3] = Point3f(  s,  s, -s)
V[4] = Point3f(  s, -s, -s)
V[5] = Point3f( -s, -s,  s)
V[6] = Point3f( -s,  s,  s)
V[7] = Point3f(  s,  s,  s)
V[8] = Point3f(  s, -s,  s)

# Create faces
F = Vector{QuadFace{Int64}}(undef,6)
F[1] = QuadFace{Int64}(1,2,3,4)
F[2] = QuadFace{Int64}(8,7,6,5)
F[3] = QuadFace{Int64}(5,6,2,1)
F[4] = QuadFace{Int64}(6,7,3,2)    
F[5] = QuadFace{Int64}(7,8,4,3)    
F[6] = QuadFace{Int64}(8,5,1,4)   

function separate_faces(V, F)
    points = Point3f[]
    faces = QuadFace{Int64}[]
    for f in F
        fp = V[f]
        n = length(points) + 1
        push!(faces, QuadFace{Int64}((n:(n+3))...))
        append!(points, fp)
    end
    return GeometryBasics.Mesh(points, faces)
end

m = separate_faces(V, F)
face_colors = [:blue, :green, :red, :yellow, :cyan, :magenta]
vert_colors = [to_color(c) for c in face_colors for i in 1:4]
mesh(m; color=vert_colors)
wireframe!(m; color=:black, transparency=true)

This should work much better than poly.

image
Kevin-Mattheus-Moerman commented 4 months ago

As far as I'm concerned this can be closed if this point is fixed:

I'm not sure why it doesn't error when passing it a single quad instead of an array of quads, so I guess that's what we should fix.