Clipping a mesh with a boundary function

[stla]

Clipping a mesh involves three ingredients: the mesh to be clipped, a function associating a scalar to each vertex of the mesh, called the boundary function, and a threshold value. Then the clipping consists in removing the vertices whose corresponding values by the boundary function are greater (or smaller) than the the threshold. This is not enough: the boundary constructed in this way must be treated in order to be regular (smooth).

For example, in this SO post, the Togliatti surface, here constructed in a box:

has been clipped to a sphere.

The first method simply consists in removing the triangles which go outside this sphere, but this gives an irregular boundary:

For this particular case of a sphere, one can achieve a nice, smooth result by using spherical coordinates. This is possible also because the mesh is a mesh of an isosurface here.

But the SO member @user255430 (who is nobody but the author of the R package rgl) gave a more general solution by introducing the clipping with a boundary function. Here the function associates to each vertex its Euclidean norm and the threshold is the sphere radius. This method is not restricted to isosurfaces, it works for any mesh.

Remark. In the C++ library CGAL, one can clip a mesh by another mesh instead of clipping a mesh by a boundary function fn. Then the clipping by the boundary function fn is identical to the CGAL clipping by a mesh of the isosurface fn < threshold.

Clipping is closely similar to intersection. Here is the intersection of a ball and a tetrahedron:

One can achieve the same result by clipping the tetrahedron to the sphere, in other words clipping with the boundary function fn(vertex) = norm(vertex).

The method described in the aforementionned SO post is public so we can borrow its ideas. Note that this is restricted to triangle meshes.

Where should we put the clipping function in Meshes.jl ? With the transforms?

[juliohm]

Thank you for detailed explanation @stla , it helps a lot 💯

Can you confirm the condition to retain vertices? Is it fn(v) == threshold or fn(v) <= threshold?

[stla]

It is fn(v) <= threshold or fn(v) >= threshold. If you use both separately you split the mesh into two parts.

I already implemented it and I noticed it is slow. This is due to these two lines:

topo = convert(HalfEdgeTopology, topology(tmesh)) # this is slow
∂₂₀ = Boundary{2,0}(topo) # this is slow as well

We need the triangles only. Isn't there a more efficient way to get them?

[juliohm]

As I explained on Discourse, check the sort=false option in the HalfEdgeTopology constructor.

[juliohm]

Let's assume that the condition to retain vertices is fn(v) <= t. In this case, it looks like we are indeed performing an intersection between the interior of the surface fn(v) == t and the other argument. The main difference is the post-processing step that takes the resulting set of disconnected vertices and fixes the gaps. Would that be the correct interpretation?

If we could split this operation into filtering + post-processing I think we would have more flexibility combining these two later? For example, assume we have the filtering step implemented as a transform that takes a mesh and returns a new mesh where fn(v) <= t. Can we then introduce a new Repair{K} that takes this mesh and fills the gaps to achieve the desired "clip" result?

[juliohm]

Regarding the slowness of the topology conversion, we need to double check if we really need to perform the conversion. We usually perform the conversion to HalfEdgeTopology when we are using a specific topological relation that is efficient for that topology. It seems that you are only using Boundary{2,0} and that is implemented for all kinds of topologies, so no need to convert?

[stla]

As I explained on Discourse, check the sort=false option in the HalfEdgeTopology constructor.

I tried topo = HalfEdgeTopology(topology(tmesh).connec; sort = false) but I didn't notice any difference of speed. Still very slow.

[stla]

It seems that you are only using Boundary{2,0} and that is implemented for all kinds of topologies, so no need to convert?

Aahh, will try that. Yes we only need the triangles

[juliohm]

Check the latest master branch the Repair{1}. I just refactored it to avoid the Boundary relation altogether. If you are only looping over the triangles, you can do for elem in elements(topology(mesh))

[stla]

I don't understand why ∂₂₀ = Boundary{2,0}(topo) is slow. When I type it in the Julia console, it is immediate.

[juliohm]

It should be immediate, why you are saying it is slow? Should we move this conversation to Zulip for quicker interaction?

[stla]

I'm tired now, I need a pause. I give you the code if you want to try.

using Meshes
using MeshViz
using GLMakie


function clipMesh(mesh, fn, bound = 0, greater = true, minvertices = 0)
  # We triangulate the mesh if necessary.
  tmesh = eltype(mesh) <: Triangle ? mesh : simplexify(mesh)

  # We add vertices as desired using the Loop subdivision.
  while nvertices(tmesh) < minvertices
    tmesh = refine(tmesh, TriSubdivision())
  end

  # We calculate the values associated to the vertices.
  verts = coordinates.(Meshes.vertices(tmesh))
  values = [fn(vert) for vert in verts]

  # We treat the particular cases.
  minvalue = minimum(values)
  maxvalue = maximum(values)
  if greater
    if bound < minvalue
      return tmesh
    end
    if bound > maxvalue
      return nothing # empty mesh
    end
  else
    if bound < minvalue
      return nothing # empty mesh
    end
    if bound > maxvalue
      return tmesh
    end
  end

  # We transform `values` in order that the vertices to keep 
  #   (hereafter called the "good" vertices) are those having 
  #   a positive value.
  values = values .- bound
  if !greater
    values = -values
  end

  # We will have to construct new vertices on some edges; this dictionary 
  #   will map these edges to the new vertex indices.
  new_vertices_map = Dict{Tuple{Int,Int},Int}()

  # This function will construct the new vertices, it takes as arguments
  #   the two indices of the modified edge and it returns the index of the 
  #   new vertex; it also updates the number of vertices.
  nverts = length(verts)
  function make_vertex(keep, discard)
    edge = (keep, discard)
    if !haskey(new_vertices_map, edge)
      positive_value = values[keep]
      negative_value = values[discard]
      alpha = positive_value / (positive_value - negative_value)
      newvertex = (1 - alpha) * verts[keep] + alpha * verts[discard]
      push!(verts, newvertex)
      nverts = nverts + 1
      new_vertices_map[edge] = nverts
    end
    return new_vertices_map[edge]
  end

  # Recall that vertices to keep are those with a positive value.
  # We'll call them the "good" vertices.
  # `tokeep` is a Boolean vector whose true values correspond 
  #   to the good vertices. 
  tokeep = values .>= 0

  # Get the initial triangles.
  ##topo = convert(HalfEdgeTopology, topology(tmesh))
  ##topo = HalfEdgeTopology(topology(tmesh).connec; sort = false)
  ##trgls = topology(tmesh).connec
  ##ntriangles = length(trgls)
  topo = topology(tmesh)
  ∂₂₀ = Boundary{2,0}(topo)
  ntriangles = nelements(topo)

  # Each triangle has a count of good vertices: 0, 1, 2 or 3.
  # Triangles with a 0 count are discarded, and those with a 3 
  #   count are left unchanged.
  counts = Vector{Int}(undef, 0)

  # We'll store the new triangles in a matrix; 
  #   it will initially contain the initial triangles.
  # We initially use a vector instead of a matrix, and we will 
  #   reshape it later, this will be more efficient than successively 
  #   applying some `hcat`.
  triangles = Vector{Int}()

  # Vector to store the connectivities corresponding to the triangles.
  connec = Vector{Connectivity}(undef, 0)

  # Let's initialize all this stuff now.
  for i in 1:ntriangles
    a, b, c = triangle = ∂₂₀(i)
    nkeep = tokeep[a] + tokeep[b] + tokeep[c]
    if nkeep > 0
      push!(counts, nkeep)
      push!(triangles, triangle...)
      push!(connec, connect((a, b, c)))
    end
  end
  # Reshape `triangles` to a matrix.
  triangles = reshape(triangles, 3, :)

  # We have to treat the triangles with one or two good vertices.
  # We start by those with one good vertex.
  singles = findall(==(1), counts)
  for i in singles
    triangle = triangles[:, i]
    positive = tokeep[triangle]
    index_to_keep = findfirst(positive)
    keep = triangle[index_to_keep]
    discard1, discard2 = deleteat!(triangle, positive)
    newindices = [make_vertex(keep, discard1), make_vertex(keep, discard2)]
    if index_to_keep != 2
      reverse!(newindices) # to deal with the orientation
    end
    connec[i] = connect((newindices[1], keep, newindices[2]))
  end

  # Now we treat the triangles with two good vertices.
  doubles = findall(==(2), counts)
  for i in doubles
    triangle = triangles[:, i]
    negative = .!tokeep[triangle]
    index_to_discard = findfirst(negative)
    if index_to_discard == 2
      reverse!(triangle) # deal with orientation
    end
    discard = triangle[index_to_discard]
    keep1, keep2 = deleteat!(triangle, negative)
    newindex = make_vertex(keep1, discard)
    connec[i] = connect((keep1, keep2, newindex))
    newconnectivity = connect((newindex, keep2, make_vertex(keep2, discard)))
    push!(connec, newconnectivity)
  end
 
  # Finally, the clipped mesh.
  cmesh = Meshes.SimpleMesh(Meshes.Point.(verts), connec)

  # There are unused vertices, we remove them.
  cmesh |> Repair{1}()
end

Example:

using MarchingCubes
using Meshes
using MeshViz
using GLMakie
using LinearAlgebra

# isosurface function 
phi = (1 + sqrt(5)) / 2
function f(x, y, z)
  - 4 * (phi^2*x^2 - y^2) * (phi^2*y^2 - z^2) * (phi^2*z^2 - x^2) + 
    (1 + 2*phi) * (x^2 + y^2 + z^2 - 1)^2
end

# make the isosurface
n = 200
voxel = zeros(Float32, n, n, n)
rg = LinRange{Float32,Int32}(-1.8, 1.8, n)
coords = convert(Vector{Float32}, collect(rg))
for k in 1:n, j in 1:n, i in 1:n
  voxel[i, j, k] = f(rg[i], rg[j], rg[k])
end
mc = MC(voxel, Int32; x = coords, y = coords, z = coords)
march(mc, Float32(0))
mesh = MarchingCubes.makemesh(Meshes, mc) # it is in Float64 mode

println("marching cubes done")

bound = convert(Float32, 3)
fn(v) = v[1] * v[1] + v[2] * v[2] + v[3] * v[3]
cmesh = clipMesh(mesh, fn, bound, false)

With this example there's no modification of some faces, there's only some discarded faces. So the code after the comment # We have to treat the triangles with one or two good vertices. does nothing.

[stla]

Have you tried ?

I tried ProfileView and I get this graphic:

All the red bars come from halfedge.jl.

[juliohm]

@stla can you please clarify which code you ran to get the profile output above? The construction of the Boundary relation object?

[stla]

using ProfileView
@profview clipMesh(mesh, fn, bound, false)

To be run outside VS code.

[juliohm]

As I mentioned we need to profile specific functions from halfedge.jl. The clipMesh function above does different calls and we can't track if it is an issue with our functions or their usage.

[stla]

It's possible to trace back with ProfileView. I didn't try yet.

[stla]

Strange, it is faster now. I just upgraded Julia to 1.8.5.

[stla]

I confirm it is fast now. And I don't get the same graph. I have some Any in red, how to solve that?

[juliohm]

These Any types you are seeing are due to some type instability. You are probably returning different types in a given function depending on different conditions and the compiler is struggling to figure out the return type.

I suggest that we convert this code into a PR with small steps as we've been doing so far. But first it would be nice to finish #412 , which is almost ready.