Alert! error when slicing

Specific Domains Visualization
,
1

My intension was to do Slicing an 3d object which is inside the stl file. Please help with below code

using FileIO
using GeometryBasics
using GLMakie  

rotate_x(point::GeometryBasics.Point{3, Float32}, angle_deg::Float64) = GeometryBasics.Point3f(
    point[1],
    cos(deg2rad(angle_deg)) * point[2] - sin(deg2rad(angle_deg)) * point[3],
    sin(deg2rad(angle_deg)) * point[2] + cos(deg2rad(angle_deg)) * point[3]
)

function clip_mesh(mesh::GeometryBasics.Mesh, lower_clip::Float64, upper_clip::Float64)
    vertices = GeometryBasics.coordinates(mesh)
    faces = GeometryBasics.faces(mesh)

    valid_vertex_indices = findall(v -> lower_clip < v[3] < upper_clip, vertices)
    valid_indices_set = Set(valid_vertex_indices)

    new_faces = []
    for face in faces
        face_indices = [Int(i) for i in face]
        if all(idx -> idx in valid_indices_set, face_indices)
            new_face_indices = map(idx -> findfirst(==(idx), valid_vertex_indices), face_indices)
            push!(new_faces, Face(new_face_indices...))
        end
    end

    new_vertices = vertices[valid_vertex_indices]
    return GeometryBasics.Mesh(new_vertices, new_faces)
end

function process_and_display_stl(input_file::String, rotation_angle_x::Float64, lower_clip::Float64, upper_clip::Float64)
    println("Processing STL file: $input_file")
    mesh = load(input_file)

    rotated_vertices = [rotate_x(v.position, rotation_angle_x) for v in GeometryBasics.coordinates(mesh)]
    rotated_mesh = GeometryBasics.Mesh(rotated_vertices, GeometryBasics.faces(mesh))

    clipped_mesh = clip_mesh(rotated_mesh, lower_clip, upper_clip)

    fig = Figure(resolution = (800, 800))
    ax = Axis3(fig[1, 1], title = "Processed Mesh Visualization")
    mesh!(ax, clipped_mesh, color = :gray, transparency = true)
    display(fig)
end

input_stl = "LowerJaw.stl"
rotation_angle = 90.0   
lower_clip = -20.0      
upper_clip = 10.0       

process_and_display_stl(input_stl, rotation_angle, lower_clip, upper_clip)

please find below output which i got as soon as i executed above code

image
image3469×596 63.3 KB

Please help
thanks and regards

2

The problem is that in [Int(i) for i in face] the i are of type OffsetInteger{-1, UInt32} which cannot be converted to Int via Int(...).

The underlying UInt32 value is stored as a field i, so you could use i.i. But then you of course remove the offset and need to compensate for it later. Since in the end the returned mesh should still use OffsetIntegers (I assume), it’s probably easier to change valid_vertex_indices from a Vector{Int64} to a Vector{OffsetInteger{-1, UInt32}} and then just work with face directly instead of face_indices.

(By the way, note that

julia> i = OffsetInteger{-1, Int}(1)
OffsetInteger{-1, Int64}(1)

julia> i.i  # 0-indexed equivalent of the 1-indexed 1
0

julia> i == 1  # isequal(i, Int64(1))
true

julia> s = Set(1); i in s
false

You could ‘fix’ the unintuitive Set behaviour using

julia> Base.hash(i::OffsetInteger{O, T}) where {O, T} = hash(i.i - O)

julia> i in s
true

but this is type piracy and could in principle lead to all kinds of problems.)

In clip_mesh there is also another issue: GeometryBasics.Face is a function without a unary method. You could probably use something like FT(new_face_indices...) where FT = typeof(face) (which should give NgonFace{3, OffsetInteger{-1, UInt32}}). Or better yet, extract FT out of Mesh{N, VT, FT}.

1 Like
3

Thanks for your suggestion,

i tried to modify my code:

using FileIO
using GeometryBasics
using GLMakie  

rotate_x(point::GeometryBasics.Point{3, Float32}, angle_deg::Float64) = GeometryBasics.Point3f(
    point[1],
    cos(deg2rad(angle_deg)) * point[2] - sin(deg2rad(angle_deg)) * point[3],
    sin(deg2rad(angle_deg)) * point[2] + cos(deg2rad(angle_deg)) * point[3]
)

function clip_mesh(mesh::GeometryBasics.Mesh, lower_clip::Float64, upper_clip::Float64)
    vertices = GeometryBasics.coordinates(mesh)
    faces = GeometryBasics.faces(mesh)

    valid_vertex_indices = findall(v -> lower_clip < v[3] < upper_clip, vertices)
    valid_indices_set = Set(valid_vertex_indices)

    new_faces = []
    for face in faces
        face_indices = [i.i for i in face]  
        if all(idx -> idx in valid_indices_set, face_indices)
            new_face_indices = map(idx -> findfirst(==(idx), valid_vertex_indices), face_indices)
            push!(new_faces, NgonFace{3, OffsetInteger{-1, UInt32}}(new_face_indices...))
        end
    end

    new_vertices = vertices[valid_vertex_indices]
    return GeometryBasics.Mesh(new_vertices, new_faces)
end

function process_and_display_stl(input_file::String, rotation_angle_x::Float64, lower_clip::Float64, upper_clip::Float64)
    println("Processing STL file: $input_file")
    mesh = load(input_file)

    rotated_vertices = [rotate_x(v.position, rotation_angle_x) for v in GeometryBasics.coordinates(mesh)]
    rotated_mesh = GeometryBasics.Mesh(rotated_vertices, GeometryBasics.faces(mesh))

    clipped_mesh = clip_mesh(rotated_mesh, lower_clip, upper_clip)

    fig = Figure(resolution = (800, 800))
    ax = Axis3(fig[1, 1], title = "Processed Mesh Visualization")
    mesh!(ax, clipped_mesh, color = :gray, transparency = true)
    display(fig)
end

input_stl = "LowerJawteeth.stl"
rotation_angle = 90.0        
lower_clip = -20.0           
upper_clip = 10.0            

process_and_display_stl(input_stl, rotation_angle, lower_clip, upper_clip)

after i executed the code, i get as attached output:

image
image3471×465 27.5 KB

Please help

thanks and regards

4

It would be a bit easier it you would post the output also between backticks, like your code. Also, ideally you could provide a minimal working example, where you also generate the input mesh, or provide the .stl file assuming it is quite small.

Anyway, the error shows that you need the second argument of the Mesh constructor to be something like a Vector{AbstractFace} instead of a Vector{Any}. You could change return GeometryBasics.Mesh(new_vertices, new_faces) into return GeometryBasics.Mesh(new_vertices, TriangleFace.(new_faces)). Alternatively (but probably a bit worse), directly use new_faces = AbstractFace[] instead of new_faces = []. The best approach would be new_faces = TriangleFace{OffsetInteger{-1, UInt32}}[], or simply FT[], cf. previous post.

5

Thanks for your reply

i tried to modify my code:

using FileIO
using GeometryBasics
using GLMakie  

rotate_x(point::GeometryBasics.Point{3, Float32}, angle_deg::Float64) = GeometryBasics.Point3f(
    point[1],
    cos(deg2rad(angle_deg)) * point[2] - sin(deg2rad(angle_deg)) * point[3],
    sin(deg2rad(angle_deg)) * point[2] + cos(deg2rad(angle_deg)) * point[3]
)

function clip_mesh(mesh::GeometryBasics.Mesh, lower_clip::Float64, upper_clip::Float64)
    vertices = GeometryBasics.coordinates(mesh)
    faces = GeometryBasics.faces(mesh)

    valid_vertex_indices = findall(v -> lower_clip < v[3] < upper_clip, vertices)
    valid_indices_set = Set(valid_vertex_indices)

    new_faces = TriangleFace{UInt32}[]  
    for face in faces
        face_indices = [i for i in face]  
        if all(idx -> idx in valid_indices_set, face_indices)
            new_face_indices = map(idx -> UInt32(findfirst(==(idx), valid_vertex_indices)), face_indices)
            push!(new_faces, TriangleFace(new_face_indices...))
        end
    end

    new_vertices = vertices[valid_vertex_indices]
    return GeometryBasics.Mesh(new_vertices, new_faces)
end

function process_and_save_stl(input_file::String, output_file::String, rotation_angle_x::Float64, lower_clip::Float64, upper_clip::Float64)
    println("Processing STL file: $input_file")
    mesh = load(input_file)

    rotated_vertices = [rotate_x(v.position, rotation_angle_x) for v in GeometryBasics.coordinates(mesh)]
    rotated_mesh = GeometryBasics.Mesh(rotated_vertices, GeometryBasics.faces(mesh))

    clipped_mesh = clip_mesh(rotated_mesh, lower_clip, upper_clip)

    println("Saving processed STL file to: $output_file")
    save(output_file, clipped_mesh)

    fig = Figure(size = (800, 800))
    ax = Axis3(fig[1, 1], title = "Processed Mesh Visualization")
    mesh!(ax, clipped_mesh, color = :gray, transparency = true)
    display(fig)
end

input_stl = "LowerJaw.stl"
output_stl = "LowerJaw_Processed.stl"  
rotation_angle = 90.0     
lower_clip = -20.0        
upper_clip = 10.0         

process_and_save_stl(input_stl, output_stl, rotation_angle, lower_clip, upper_clip)

And i get output as follows:

Processing STL file: LowerJaw.stl
MethodError: no method matching trailing_zeros(::OffsetInteger{-1, UInt32})

Closest candidates are:
  trailing_zeros(::BigInt)
   @ Base gmp.jl:576
  trailing_zeros(::SIMD.Vec{<:Any, <:Union{Int16, Int32, Int64, Int8, UInt16, UInt32, UInt64, UInt8}})
   @ SIMD C:\Users\OMG\.julia\packages\SIMD\cST3l\src\simdvec.jl:156
  trailing_zeros(::Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8})
   @ Base int.jl:441


Stacktrace:
  [1] hash
    @ .\float.jl:686 [inlined]
  [2] hash
    @ .\hashing.jl:30 [inlined]
  [3] hashindex
    @ .\dict.jl:156 [inlined]
  [4] ht_keyindex
    @ .\dict.jl:268 [inlined]
  [5] haskey
    @ .\dict.jl:569 [inlined]
  [6] in
    @ .\set.jl:92 [inlined]
  [7] #2
    @ .\In[1]:25 [inlined]
  [8] #214
    @ C:\Users\OMG\.julia\packages\StaticArrays\9Yt0H\src\mapreduce.jl:299 [inlined]
  [9] macro expansion
    @ C:\Users\OMG\.julia\packages\StaticArrays\9Yt0H\src\mapreduce.jl:174 [inlined]
 [10] _mapfoldl
    @ C:\Users\OMG\.julia\packages\StaticArrays\9Yt0H\src\mapreduce.jl:149 [inlined]
 [11] _mapreduce
    @ C:\Users\OMG\.julia\packages\StaticArrays\9Yt0H\src\mapreduce.jl:147 [inlined]
 [12] all
    @ C:\Users\OMG\.julia\packages\StaticArrays\9Yt0H\src\mapreduce.jl:299 [inlined]
 [13] clip_mesh(mesh::GeometryBasics.Mesh{3, Float32, GeometryBasics.Ngon{3, Float32, 3, Point{3, Float32}}, SimpleFaceView{3, Float32, 3, OffsetInteger{-1, UInt32}, Point{3, Float32}, NgonFace{3, OffsetInteger{-1, UInt32}}}}, lower_clip::Float64, upper_clip::Float64)
    @ Main .\In[1]:25
 [14] process_and_save_stl(input_file::String, output_file::String, rotation_angle_x::Float64, lower_clip::Float64, upper_clip::Float64)
    @ Main .\In[1]:47
 [15] top-level scope
    @ In[1]:68

Please help

Original stl file content
Please have the screenshots attached below :

image
image1217×706 134 KB

image
image1134×675 99.5 KB

image
image1138×648 54.9 KB

image
image951×657 105 KB

Expected final stl file content
Please have the screenshots attached below :

image
image1095×373 60.6 KB

image
image1086×506 83.1 KB

image
image1120×459 74.9 KB

Please guide me how to proceed further.

Thanks and regards

6

Note you can do mesh slicing with Comodo.jl, here is a demo for it: Comodo.jl/examples/demo_trisurfslice.jl at main · COMODO-research/Comodo.jl · GitHub

You can import an STL, slice it, and export a new STL if you like.

2 Likes
7

Using an existing library will indeed be the easiest and most likely most performant option.

I meant to upload the actual .stl file, not screenshots thereof. At this moment we cannot execute your code with your data, in order to replicate and debug your issues.

This essentially just copies face; face_indices will be a Vector{OffsetInteger{-1, UInt32}}. Equivalently you could also just use face instead of face_indices everywhere. But OffsetIntegers are not compatible with Sets, as mentioned above. If you implement Base.hash(::OffsetInteger) this should work fine.

If you’re uncomfortable with type piracy, you can use something like
to_idx(oi::OffsetInteger{O, T}) where {O, T} = oi.i - O and face_indices = to_idx.(face), which gives you a Vector{UInt32} of 1-indexed indices.

Do note that using OffsetIntegers “helps reduce copying when communicating with 0-indexed systems such as OpenGL”, which will not be the case if you use TriangleFace{UInt32}. GLMakie visualisation seems fine, though.

8

Thanks for your reply

iam unable to upload stl file here. i get message like below:

image
image1307×281 14.6 KB

please help

thanks and regards

9

You could upload it to some other platform.

But if you can live with artificial data, that’s probably easiest. E.g. mesh = GeometryBasics.mesh(Sphere(Point3f(0., 0., 0.), 1.)). The code below works for me.

Clip unit sphere to top half 
using MeshIO, FileIO, GeometryBasics, GLMakie

to_idx(oi::OffsetInteger{O, T}) where {O, T} = oi.i - O

function clip_mesh(mesh, lower_clip, upper_clip)
    vertices = coordinates(mesh)
    fcs = faces(mesh)
    valid_vertex_indices = findall(v -> lower_clip < v[3] < upper_clip, vertices)
    valid_indices_set = Set(valid_vertex_indices)

    new_faces = TriangleFace{UInt32}[]
    for face in fcs
        face_indices = to_idx.(face)
        if all(idx -> idx in valid_indices_set, face_indices)
            new_face_indices = map(idx -> UInt32(findfirst(==(idx), valid_vertex_indices)), face_indices)
            push!(new_faces, TriangleFace(new_face_indices...))
        end
    end

    new_vertices = vertices[valid_vertex_indices]
    return GeometryBasics.Mesh(new_vertices, new_faces)
end

function visualize_mesh(mesh)
    display(GLMakie.mesh(mesh, color = :gray, transparency = true))
end

function main()
    mesh = GeometryBasics.mesh(Sphere(Point3f(0., 0., 0.), 1.))
    clipped_mesh = clip_mesh(mesh, 0. - eps(), 1. + eps())
    visualize_mesh(clipped_mesh)
end

main()