vedo icon indicating copy to clipboard operation
vedo copied to clipboard

How to load/get only the surface of a solidified mesh

Open ttsesm opened this issue 2 years ago • 22 comments

Hi Marco,

I have a mesh which also has internal fragmentation but as far as I know it is not tetrahedralized and it is a normal .obj file (see attached). image

I want to load it or after load to get only the surface structure and ignore the internal points/faces. I've tried to use the .tomesh(fill=False) method but apparently there isn't anything relevant for the Mesh class. Thus, is there any other efficient way (build in function, I would prefer it from looping over the points and checking whether the point is inside the mesh) to do that.

Thanks piece_1.zip .

ttsesm avatar Jul 04 '22 17:07 ttsesm

Hi, It's not a trivial problem unfortunately the mesh seems non-manifold, so i'm not sure how you can do it, i tried

from vedo import *

settings.useDepthPeeling=1

m = Mesh('/home/musy/Downloads/piece_1.obj')
m.computeNormals(consistency=False)
m.color('k5').alpha(0.3).flat()

pts = m.cellCenters()
n = m.celldata["Normals"]

pp = Points(pts+n/10000)
ids = m.insidePoints(pp, returnIds=True)
print(ids)
# m.deleteCells(ids) # core dumps
pp_in = Points(pts[ids], c='r')

show(m, pp_in, axes=1)

Screenshot from 2022-07-04 22-13-53

which identifies most of the internal cells but then it crashes when trying to remove them. I would investigate with pymeshlab if there is anything that's ready to use for such a case.

marcomusy avatar Jul 04 '22 20:07 marcomusy

Thanks, I will have a look as well.

ttsesm avatar Jul 05 '22 08:07 ttsesm

I found this article http://meshlabstuff.blogspot.com/2009/04/how-to-remove-internal-faces-with.html but it doesn't seem to be consistent :thinking:

Also the attached file is a bit more interesting case to be tested (this is by using your code from above): image

small_dino.zip

ttsesm avatar Jul 05 '22 10:07 ttsesm

Hi Marco,

A small update, I hope it can be useful. I've tried to play a bit with the boundary() function and it seems to work quite nicely at least for the given use cases. Considering what you mentioned that the meshes seem to be non-manifold I've tried to extract the boundaries with the nonManifoldEdges= parameter set to True. Thus, as you can see bellow it successfully extracted these internal geometry: ezgif com-gif-maker(2) Then I've tried to activate the returnCellIds= or returnPointIds= in order to get the ids of these faces/points to be removed but I was getting a strange array structure which didn't make sense plus the following message:

Error evaluating: thread_id: pid_1216403_id_140314406714048
frame_id: 94146529420160
scope: FRAME
attrs: b
Traceback (most recent call last):
  File "/usr/share/pycharm/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_vars.py", line 314, in resolve_compound_variable_fields
    return _typeName, _resolve_default_variable_fields(var, resolver, offset)
  File "/usr/share/pycharm/plugins/python-ce/helpers/pydev/_pydevd_bundle/pydevd_vars.py", line 264, in _resolve_default_variable_fields
    return resolver.get_dictionary(VariableWithOffset(var, offset) if offset else var)
  File "/usr/share/pycharm/plugins/python-ce/helpers/pydev/pydevd_plugins/extensions/types/pydevd_plugin_numpy_types.py", line 93, in get_dictionary
    ret['min'] = obj.min()
  File "/home/ttsesm/Development/breaking_bad/venv_bb/lib/python3.9/site-packages/numpy/core/_methods.py", line 44, in _amin
    return umr_minimum(a, axis, None, out, keepdims, initial, where)
ValueError: zero-size array to reduction operation minimum which has no identity

I do not know whether this could be of any help but it should be possible to remove this internal geometry based on the above discovery. The snippet I used:

m = vd.Mesh(mesh)#.tetralize().tomesh(fill=False)
b = m.boundaries(nonManifoldEdges=True) # activate returnCellIds=True for getting the error message
m.color('k5').alpha(0.3).flat()

vd.show(m, b, "Whole mesh", at=0, N=2, axes=1)
vd.show(b, "Only nonManifoldEdges boundaries", at=1, interactive=True).close()

ttsesm avatar Jul 05 '22 12:07 ttsesm

Yes ! We got the same idea:

pip install -U git+https://github.com/marcomusy/vedo.git

Then

from vedo import *

settings.useDepthPeeling=1

# m = Mesh('/home/musy/Downloads/piece_1.obj')
m = Mesh('clipped.vtk')
m.computeNormals(consistency=False)
m.color('k5').alpha(0.3).flat()

bb = m.boundaries(nonManifoldEdges=True, boundaryEdges=False).lw(2)
bb_ids = m.boundaries(nonManifoldEdges=True, boundaryEdges=False, returnCellIds=True)
bb_ids = np.unique(bb_ids).astype(int)
print(bb_ids)
# m.deleteCells(bb_ids).clean()

labs = m.labels('id', cells=True)
show(m, bb, labs, axes=1)

Screenshot from 2022-07-05 14-41-52

we are getting close to it but it looks like the problem is now that the cells ID are not the ones that we would like to kill. Maybe you can dig a bit more into it ? clipped.zip

marcomusy avatar Jul 05 '22 12:07 marcomusy

Yes, I am looking at it. For some reason it seems that the returning face ids are wrong... :thinking:

how can I create a separate mesh from the returning face ids to visualize only these faces to confirm that they correspond to what the boundary output is? using m.faces()[bb_ids] doesn't work. I guess I would need to get the point ids as well and then do:

points = [m.points()[i] for in in pp_ids]
faces = [m.faces()[i] for i in bb_ids]

mm = Mesh(points, faces)

ttsesm avatar Jul 05 '22 14:07 ttsesm

ok, investigating a bit more I've noticed two things:

  1. The returned cell ids from the boundaries(returnCellIds=True) method are mostly correct except that there are also some small faces which are not on the boundary but are returned as well. Check on the screenshots bellow: image image image You can clearly see these small faces which my guess is that shouldn't be returned.

  2. Then when you apply the deleteCells() function for some reason it deletes these small faces but keeps the actual boundary faces that should have been deleted instead. image

That's quite a strange behavior. I will try to see if I can find something inside the boundaries() source code.

ttsesm avatar Jul 06 '22 10:07 ttsesm

The returnPointIds=True parameter seems to return the correct boundary points (points in red dots): image

but the returnCellIds=True returns these extra non-boundary faces (or better partially boundary since one edge is usually as I see it always on a boundary) :thinking:

If I use the deleteCellsByPointIndex() instead of the deleteCells() it removes the internal geometry but also the external one which is associated with these points.

ttsesm avatar Jul 06 '22 13:07 ttsesm

interesting findings, you may need to unpack the two methods to see why the cell ids in the first case seem wrong. The deleteCellsByPointIndex also deletes those extra triangles because they are attached to at least one red point, so that is the (undesired but) correct behavior.

..or just wait one week and I will look into it :)

marcomusy avatar Jul 06 '22 16:07 marcomusy

interesting findings, you may need to unpack the two methods to see why the cell ids in the first case seem wrong.

Yes, I was trying that. Didn't manage to find much though :-(.

The deleteCellsByPointIndex also deletes those extra triangles because they are attached to at least one red point, so that is the (undesired but) correct behavior.

Yes, that's fine. I understood that this is the correct behavior.

..or just wait one week and I will look into it :)

Well sure, if I do not manage to get anything until then most likely I will need your support ;-).

ttsesm avatar Jul 07 '22 10:07 ttsesm

@marcomusy I didn't manage to have any progress on this, thus if you could have a look on it when you find some time I would appreciated it ;-).

ttsesm avatar Jul 11 '22 07:07 ttsesm

Hi I just pushed a new version. It is not perfect in your case because the non-manifold regions are very close to each other so there are faces that get deleted only because all their vertices belong to some non-manifold face even if an edge does not belong to a non-manifold face. Check out:

from vedo import *

settings.useDepthPeeling=1

m = Mesh('piece_1.obj').clean()
m.color('k5').flat()

# bb = m.boundaries(nonManifoldEdges=True, boundaryEdges=False)
ids = m.boundaries(nonManifoldEdges=True, boundaryEdges=False, returnCellIds=True)

m.deleteCells(ids).clean()

show(m, axes=1)

marcomusy avatar Jul 18 '22 18:07 marcomusy

Hi Marco, thanks for the update. Interesting, yes it still doesn't work perfectly. At least now the inner cells are removed. However, now the detected edges also look quite different. Are you sure that they are correct?

Because for example now I get these isolated edges as well as the edges in the feet of the dino which I was not getting previously (https://user-images.githubusercontent.com/10189018/177328460-7fb10a56-5de0-40f8-a478-b63152bb2222.gif): image

ttsesm avatar Jul 19 '22 15:07 ttsesm

Hi - yes you are right - I just fixed it in the master!

marcomusy avatar Jul 19 '22 15:07 marcomusy

Thanks, indeed now is better. I still get the holes on the surface as you mentioned, thus I've tried to apply some hole filling with fillholes() but it doesn't seem to do anything: image

I tried to play with different sizes nothing seems to change :thinking:

bb_ids = m.boundaries(nonManifoldEdges=True, boundaryEdges=True, returnCellIds=True)
m.deleteCells(bb_ids).clean()
m.fillHoles(size=50)

ttsesm avatar Jul 20 '22 11:07 ttsesm

A small update on this, considering that I cannot avoid the holes in the output. I've managed to resolve the holes issue with the meshfix lib as described here. Any other approach to fill the holes, either with vedo, or other libs trimesh, pymeshlab, etc. failed.

Output with meshfix.repair(): image

ttsesm avatar Oct 07 '22 15:10 ttsesm

nice there is already an example for meshfix in the vedo library too: examples/other/remesh_meshfix.py

marcomusy avatar Oct 07 '22 16:10 marcomusy

you may want to have a look at meshlab as well: https://www.meshlab.net/

there are python bindings (pymeshlab) which are really easy to combine with vedo as the data-structure is identical:

# from vedo to pymeshlab
ms = pymeshlab.MeshSet()
mesh = pymeshlab.Mesh(vertices, faces)
ms.add_mesh(mesh)

# and back
vertices = ms.current_mesh().vertex_matrix()
faces = ms.current_mesh().face_matrix()

m2 = vedo.Mesh([vertices, faces])
p.add(m2, render=False)

I use it myself for boolean operations as they are more stable than the ones in vtk (https://github.com/RubendeBruin/pymeshup)

RubendeBruin avatar Oct 08 '22 12:10 RubendeBruin

We already have examples for interfacing pymeshlab in examples/other/pymeshlab1.py which works for pymeshlab-2021.10, but now with the latest release it seems that the their api changed so i'll need to update it..

marcomusy avatar Oct 08 '22 15:10 marcomusy

@marcomusy if I remember correctly you also had an interface for porting from open3d to vedo, right? but I cannot find it anymore :thinking:

ttsesm avatar Oct 10 '22 09:10 ttsesm

Actually I never looked into open3d!

marcomusy avatar Oct 10 '22 09:10 marcomusy

Actually I never looked into open3d!

Ok, then maybe I've seen it somewhere else... :thinking:

Anyways, actually I've created the following functions in case you are interested to add them. It is as simple as that:

import vedo as vd
import open3d as o3d

def o3d2vedo(o3d_mesh):
    m = vd.Mesh([np.array(o3d_mesh.vertices), np.array(o3d_mesh.triangles)])

    # you could also check whether normals and color are present in order to port with the above vertices/faces
    return m

and the counter function:

import vedo as vd
import open3d as o3d

def vedo2open3d(vd_mesh):
    """
    Return an `open3d.geometry.TriangleMesh` version of
    the current mesh.

    Returns
    ---------
    open3d : open3d.geometry.TriangleMesh
      Current mesh as an open3d object.
    """
    # create from numpy arrays
    o3d_mesh = o3d.geometry.TriangleMesh(
        vertices=o3d.utility.Vector3dVector(vd_mesh.points()),
        triangles=o3d.utility.Vector3iVector(vd_mesh.faces()))

    # I need to add some if check here in case color and normals info are not existing
    # o3d_mesh.vertex_colors = o3d.utility.Vector3dVector(vd_mesh.pointdata["RGB"]/255)
    # o3d_mesh.vertex_normals = o3d.utility.Vector3dVector(vd_mesh.pointdata["Normals"])

    return o3d_mesh

ttsesm avatar Oct 10 '22 09:10 ttsesm