marcomusy / vedo

A python module for scientific analysis of 3D data based on VTK and Numpy
https://vedo.embl.es
MIT License
2.06k stars 267 forks source link

Mesh split and disassemble/assemble animation #498

Open ttsesm opened 3 years ago

ttsesm commented 3 years ago

Hi @marcomusy,

I have a mesh which I would like to split in a random number of sub-meshes as you do in this example but without discarding any parts. If you see in the image below at the current state some of the parts after the split are removed and only 40 pieces are kept.

import vedo as vd

em = vd.Volume(vd.dataurl + "embryo.tif").isosurface(80)
splitem = em.splitByConnectivity(maxdepth=40)#[0:9]

vd.show(em, "Before", at=0, N=2, axes=1)
vd.show(splitem, "After", at=1, interactive=True).close()

image

Is there a flag for the maxdepth= parameter to consider the complete initial mesh? Ideally I would like the complete mesh to be split in 40 (or whichever number is given each time) pieces.

Thereafter, I would like to create a kind of animation similar to this one where I move each part aside in the space (or as an inventory list on a side pane) and then put them again together giving this assemply/disassemply outlook. Is there any example which I could follow? If not how easy would that be to do it with vedo?

Thanks.

marcomusy commented 3 years ago

Is there a flag for the maxdepth= parameter to consider the complete initial mesh?

just choose a large number.

Long time ago I did this (inspired by the famous manim package):


from vedo import dataurl, Volume, Axes, show
from vedo.applications import Animation

iso = Volume(dataurl+"embryo.tif").isosurface(80)

meshes = iso.splitByConnectivity(maxdepth=40)

axes = Axes(meshes[0]) # build axes manually

for m in meshes:
    cm = m.centerOfMass()
    m.pos(cm*7) # explode meshes
# show(meshes, axes=1).close()

anim = Animation()     # a vedo.Plotter object
anim.timeResolution = 0.01  # secs

anim.fadeIn(meshes[1:], t=0, duration=0.2)
anim.fadeIn(meshes[0],  t=1, duration=1)
# anim.fadeIn(axes, t=0, duration=0.2) # will not work right now!

for m in meshes:
    anim.move(m, (0,0,0), style="quadratic")
    anim.rotate(m, axis="y", angle=360)

anim.totalDuration = 5 # can now shrink/expand total duration
anim.play()

# edit with https://ezgif.com/

embrio_splt

I thought nobody would ever care about.. so it's a bit outdated and not well tested, but it should work for simple things, or you can create your own animation by moving objects in a loop. Check out examples/other/animation1.py.

ttsesm commented 3 years ago

Thanks Marco,

Appreciated, the animation looks quite cool I will play with it and see what kind of adaptations I might be able to do.

just choose a large number.

It wouldn't make more sense the mesh to be split to the exact given number instead of just giving a large number. This might also be not desired, since someone might just want to split the mesh in 3-4 big pieces instead and keeping all the initial structure.

ttsesm commented 3 years ago

Actually to make it a bit more explicit of what I what I want to do check on the following. Imagine that I have some meshes as one .obj files which each has been reassembled together from individual single meshes and of which I would like obtain separately. For example: image image

The above mesh as you can see has 6 submeshes, which I would like to individualize.

DoraColumn_Reassembled.zip Reassembled_Tombstone.zip

marcomusy commented 3 years ago

Uhm since they have a "GroupIds" array you can also use mesh.threshold() to extract any parts, or split them normally with splitByConnectivity() as discussed above... not sure if answering your question.

ttsesm commented 3 years ago

Uhm since they have a "GroupIds" array you can also use mesh.threshold() to extract any parts, or split them normally with splitByConnectivity() as discussed above... not sure if answering your question.

Hmm, actually the splitByConnectivity() should be working but for some reason the .obj loader does not load the mesh correctly:

m = vd.Mesh("DoraColumn_Reassembled.obj")

# return the list of the largest 6 connected meshes:
splitem = m.splitByConnectivity(maxdepth=6)

vd.show(m, "Before", at=0, N=2, axes=1)
vd.show(splitem, "After", at=1, interactive=True).close()

image

You can see one of the segmented parts is loaded "wrongly" which leads to a faulty split.

marcomusy commented 3 years ago

I don't quite understand why that happens... maybe the pieces are slightly overlapping so that they cannot be split (?)

The thresholding works though:

import vedo as vd

path = "data/vedo/DoraColumn_Reassembled/"
m = vd.Mesh(path+"DoraColumn_Reassembled.obj").print()

m0 = m.clone().threshold("GroupIds",-0.5,0.5, on="cells").texture(path+"DoraColumn1_med.jpg")
m1 = m.clone().threshold("GroupIds", 0.5,1.5, on="cells").texture(path+"DoraColumn2_med.jpg")
m2 = m.clone().threshold("GroupIds", 1.5,2.5, on="cells").texture(path+"DoraColumn3_med.jpg")
m3 = m.clone().threshold("GroupIds", 2.5,3.5, on="cells").texture(path+"DoraColumn4_low.jpg")
m4 = m.clone().threshold("GroupIds", 3.5,4.5, on="cells").texture(path+"DoraColumn5_low.jpg")
m5 = m.clone().threshold("GroupIds", 4.5,5.5, on="cells").texture(path+"DoraColumn6_low.jpg")

vd.show(m, at=0, shape="6/1", axes=1, sharecam=0)
vd.show(m0, "piece 0", at=1)
vd.show(m1, "piece 1", at=2)
vd.show(m2, "piece 2", at=3)
vd.show(m3, "piece 3", at=4)
vd.show(m4, "piece 4", at=5)
vd.show(m5, "piece 5", at=6)
vd.interactive()

Screenshot from 2021-11-04 17-05-13

ttsesm commented 3 years ago

Thanks a lot. Yes, it is really strange. Actually every obj file that I load it has one of the pieces always multicolored, I do not think it is overlapping or something. If you try the Tombstone attachment you will see the same issue.

marcomusy commented 3 years ago

It is multicolored because vedo tries to colorize the texture coords of each piece (try press 4 to cycle through the existing point arrays).

I dont have a good explanation for the split method..

ttsesm commented 3 years ago

Using the 4 to cycle through the existing point arrays finds all the 6 ones sub-meshes. I do not understand why the split method does not work either.

The threshold alternative is ok but it considers too much manual work :-(. Also any idea why the texture was not loaded correctly?

marcomusy commented 3 years ago

I still think that the meshes might be overlapping ..

Also any idea why the texture was not loaded correctly?

why do you say that?

ttsesm commented 3 years ago

They should not, because in Meshlab when I use the select by faces tool (which selects all the connected faces) and pick faces on each one of the sub-meshes I get the individual pieces. Thus this means that they should not, otherwise if two pieces were overlapping I would be getting both selected.

Because in your example with the threshold alternative while you are loading the texture file still the texture is not showing on the objects.

ttsesm commented 2 years ago

Hi @marcomusy,

Why the position of the object is not set o (0,0,0): image

I am loading the file with mesh = vd.Mesh(mesh_file, c=color).pos(x=0, y=0, z=0) but you see that at the z axis the object is not centered and from the print() I see that the center of mash in not at (0, 0, 0). It should be there, shouldn't be?

Mesh/Points 
           file: ../data/tombstone/Tombstone1_low.obj
          color: teal1, rgb=(0.0824, 0.345, 0.247), alpha=1.0
         points: 150000
          cells: 50000
       polygons: 50000
       position: (0.0, 0.0, 0.0)
          scale: (1.00, 1.00, 1.00)
 center of mass: (-31.2, 30.3, -151)
   average size: 184.192
  diagonal size: 724.139
         bounds: x=(-235, 241) y=(-76.8, 132) z=(-392, 112)
  clicked point: [ 219.021061      9.90228517 -198.20774671]
    scalar mode: Default   coloring = Default
   active array: material_0 (point data)  MaterialIds (cell data)
     point data: name=material_0 (2 FLOAT, np.float32), range=(8.80e-4,0.999)
     point data: name=Normals (3 FLOAT, np.float32), range=(-1.00,0.998)
      cell data: name=MaterialIds (1 INT), range=(0,0)
      cell data: name=GroupIds (1 FLOAT), range=(0,0)
marcomusy commented 2 years ago

It looks normal. A Mesh has vertices whose coordinates can have any arbitrary offset! If you wish to center the mesh so that its center of mass coincides with the cartesian origin do: mesh.shift(-mesh.centerOfMass())

ttsesm commented 2 years ago

Is it possible to retrieve the center points of the cells in a grid?

grid = Grid(sx=sx, sy=sy, resx=gridRes[0], resy=gridRes[1])

I have created a grid based on the amount of objects that I have and I would like to create an animation where each object moves to each own grid cell: image

marcomusy commented 2 years ago

yes! grid.cellCenters()

ttsesm commented 2 years ago

Hhmm, I think I am doing something wrong or missing something: Peek 2022-01-14 11-37

    anim = Animation()  # a vedo.Plotter object
    anim.timeResolution = 0.01  # secs
    axes = Axes(mesh_models[0])  # build axes manually

    anim.fadeIn(mesh_models, t=0, duration=0.2)
    anim.fadeIn(grid, t=0, duration=0.2)
    # anim.fadeIn(axes, t=0, duration=0.2) # will not work right now!

    for i, mm in enumerate(mesh_models):
        anim.move(mm, grid.cellCenters()[i])

    anim.totalDuration = 5  # can now shrink/expand total duration
    anim.play()

I do not understand why the objects move across the grid instead of along :thinking:

I tried also with a plotter (since I would like to have the axes) but I am also getting something strange: Peek 2022-01-14 11-49

    # Setup the scene
    vd.show(mesh_models, grid, axes=1, viewup="z", interactive=0)

    for t in np.arange(0, 1, 0.005):
        for i, mesh_model in enumerate(mesh_models):
            mesh_model.pos(grid.cellCenters()[i])
            plotter = vd.show(mesh_model, grid, axes=1)
        if plotter.escaped:
            break  # if ESC button is hit during the loop

    interactive().close()

I guess here I would need to extract the intermediate points from (0, 0, 0) to each cell center and update in each loop, but again it seems that the given positions are somehow wrong.

Any idea?

marcomusy commented 2 years ago

Please whenever possible post the full working code. This seems to work ok:

from vedo import *
from vedo.applications import Animation

mesh_models = load("tombstone/Tombstone*_low.obj")

grid = Grid(sx=1500, sy=2000, resx=3, resy=2)

gpts = Points(grid.cellCenters())

anim = Animation()  # a vedo.Plotter object
anim.timeResolution = 0.01  # secs
anim.totalDuration = 5      # can shrink/expand duration

anim.fadeIn(mesh_models, t=0, duration=0.2)
anim.fadeIn(grid, t=0, duration=0.2)
anim.fadeIn(gpts, t=0, duration=0.2)

for i, mm in enumerate(mesh_models):
    mm.rotateX(70).color(i)
    anim.move(mm, grid.cellCenters()[i])

anim.show() # needed on mac OSX, press q
anim.play()
Screenshot 2022-01-14 at 20 41 20
ttsesm commented 2 years ago

Ok, I the problem was on the transformation that I am applying and not reseting it:

def transform_mesh(m):
    cm = m.centerOfMass()
    m.shift(-cm)
    elli = pcaEllipsoid(m, pvalue=0.5)

    ax1 = versor(elli.axis1)
    ax2 = versor(elli.axis2)
    ax3 = versor(elli.axis3)

    T = np.array([ax1, ax2, ax3])  # the transposed matrix is already the inverse
    # print(T)
    # print(T@ax1)

    return m.applyTransform(T, reset=True) # <-- I had to enable reset

and the code for doing the same with a plotter is the following:

def update_pos(points, n_points):
    p1, p2 = points
    diff = p2 - p1
    t = np.linspace(0, 1, n_points+20)[1:-19]
    return (p1[np.newaxis, :] + t[:, np.newaxis] * diff[np.newaxis, :]).squeeze()

vd.show(mesh_models, grid, axes=1, viewup="y", interactive=0)

for t in np.arange(0, 1, 0.005):
    for i, mesh_model in enumerate(mesh_models):
        mesh_model.pos(update_pos([mesh_model.pos(), grid.cellCenters()[i]], 1))
        plotter = vd.show(mesh_models, grid)

    if plotter.escaped:
        break  # if ESC button is hit during the loop

interactive().close()

anim

Btw, why is the camera viewAngle setting doesn't seem to work? I've tried to change the scene setting to the following vd.show(mesh_models, grid, axes=1, camera={'viewAngle': (45)}, viewup="z", interactive=0) in order to get a view similar to this: image but it doesn't seem to have any affect to the angle view since I am getting the default one: image

marcomusy commented 2 years ago

either you specify viewup or the camera ... what I would do is to choose one orientation that you like manually in a scene then press C and copy-paste the parameters that are printed out.

Also you can remove the z-axis just by setting axes=dict(ztitle="").

ttsesm commented 2 years ago

Hi @marcomusy,

Aside of the disassemble/assemble animation is there a way to fracture a mesh into multiple pieces. Something similar to this addon from Blender https://www.youtube.com/watch?v=ZG_ZMnKzVTQ.

Imagine you have a solid (or non solid, which I guess it would need to be solidified somehow maybe by using the volume???) sphere (this could be any 3D mesh) and then I want to fracture this in random pieces where each piece is fractured differently from the other. Then when you put the pieces back together you should get the original sphere.

The fracturing though I would prefer not to be in a sharp cut like a slice, though it might contain such fractures but not only such. I would like to have something similar like when you break an object and you can get any possible fracture as kind of exemplified in the blender addon in the video.

marcomusy commented 2 years ago

yes. try:

from vedo import *
import numpy as np
np.random.seed(0)

s = Sphere(quads=True, res=15).clean()
res = 0.02  #control the tetras resolution

# fill the space w/ points
pts = (np.random.rand(10000, 3)-0.5)*2
fillpts = s.insidePoints(pts).subsample(res).scale(0.9)  # make pts uniform
seeds = s.clone().subsample(0.3).ps(12).c('black') # pick uniform pts on sphere
printc("# of pieces, #points in sphere:", seeds.N(), fillpts.N())

tmesh = delaunay3D(merge(fillpts,s))
# assign a closest point to each tetrahedron
cids = []
for p in tmesh.cellCenters():
    cid = seeds.closestPoint(p, returnPointId=True)
    cids.append(cid)

tmesh.celldata["fragment"] = np.array(cids)

pieces = []
for i in range(seeds.NPoints()):
    tc = tmesh.clone().threshold(above=i-0.1, below=i+0.1)
    mc = tc.tomesh(fill=False).color(i)
    pieces.append(mc)

############### animate
plt = Plotter(interactive=1)
plt.show(pieces, seeds, "press q to make it explode")
for i in range(15):
    for pc in pieces:
        cm = pc.centerOfMass()
        pc.shift(cm/15)
    plt.render()
plt.interactive().close()
Screenshot 2022-03-07 at 19 11 41

press q:

Screenshot 2022-03-07 at 19 12 02
marcomusy commented 2 years ago

PS: because it relies in delauney3D() the above example only works for convex shapes, otherwise you need to provide a tetrahedral mesh (vtk doesn't do that).

ttsesm commented 2 years ago

This is nice, thanks. Let me understand it though.

You create the mesh (in this case the sphere), and then you fill the space with random points: image

Then you keep the points falling inside the mesh, and from these points you subsample the number of pieces.

From there you use the original mesh and the inside points to create a volume tetrahedral mesh based on the delauney3D triangulation: image

Thereafter you create the pieces by labeling each tetrahedral face as to the closest seed point. image image

Some questions though:

  1. Can you elaborate a bit on the command tc = tmesh.clone().threshold(above=i-0.1, below=i+0.1) I did not understand how the threshold() method works.
  2. Can I get a smoother resolution on the fragmented sides, I tried to increase the resolution but still the cuts look too sharp. I guess this might be related to what you mentioned about using the delauney method and the convex shapes or is it how the fillpts points are distributed inside the shape?

    otherwise you need to provide a tetrahedral mesh (vtk doesn't do that)

  3. Can you also explain a bit about that and what do you mean?
marcomusy commented 2 years ago
  1. that makes a copy then extracts all (and only) tets that have scalar between i-0.1 and i+0.1, so it picks the individual pieces.
  2. if you subset more you get less internal points and less surface points (seed), so in a way that controls resolution
  3. deluaney3d only fill with tets convex surfaces (or it makes them convex by "wrapping" them into the tightest convex shape). If you start from a polygonal concave surface you need some external package to generate the tetrahedral mesh. The you can feed it to the above pipeline, eg:
"""Add a custom scalar to a TetMesh to segment it.
Press q to make it explode"""
from vedo import TetMesh, Plotter, dataurl, printc

n = 20000
f1 = 0.005  # control the tetras resolution
f2 = 0.15   # control the nr of seeds

tmesh = TetMesh(dataurl+'limb_ugrid.vtk')
surf = tmesh.tomesh(fill=False)

# pick uniform pts on the surface
seeds = surf.clone().subsample(f2).ps(10).c('black')
printc("#pieces:", seeds.N())

# assign to each tetrahedron the id of the closest seed point
cids = []
for p in tmesh.cellCenters():
    cid = seeds.closestPoint(p, returnPointId=True)
    cids.append(cid)
tmesh.celldata["fragment"] = cids

pieces = []
for i in range(seeds.NPoints()):
    tc = tmesh.clone().threshold(name="fragment", above=i-0.1, below=i+0.1)
    mc = tc.tomesh(fill=False).color(i)
    pieces.append(mc)

############### animate
plt = Plotter(size=(1200,800), axes=1)
plt.show(__doc__, pieces)
for i in range(20):
    for pc in pieces:
        cm = pc.centerOfMass()
        pc.shift(cm/25)
    plt.render()
plt.interactive().close()

Screenshot from 2022-03-08 14-40-24

marcomusy commented 2 years ago

PS: One can use pymeshfix + tetgen to generalize, eg:

"""Segment a TetMesh with a custom scalar.
Press q to make it explode"""
from vedo import Mesh, TetMesh, Plotter, Text2D, dataurl
import tetgen
import pymeshfix

n = 20000
f1 = 0.005  # control the tetras resolution
f2 = 0.15   # control the nr of seeds

# repair and tetralize the closed surface
amesh = Mesh(dataurl+'bunny.obj')
meshfix = pymeshfix.MeshFix(amesh.points(), amesh.faces())
meshfix.repair() # will make it manifold
repaired = Mesh(meshfix.mesh)
tet = tetgen.TetGen(repaired.points(), repaired.faces())
tet.tetrahedralize(order=1, mindihedral=20, minratio=1.5)
tmesh = TetMesh(tet.grid)

surf = tmesh.tomesh(fill=False)
txt = Text2D(__doc__, font="Brachium")

# pick points on the surface and use subsample to make them uniform
seeds = surf.clone().subsample(f2).ps(10).c('black')

# assign to each tetrahedron the id of the closest seed point
cids = []
for p in tmesh.cellCenters():
    cid = seeds.closestPoint(p, returnPointId=True)
    cids.append(cid)
tmesh.celldata["fragment"] = cids
#tmesh.celldata.select("fragment")# bug, has no effect, needs name=...

pieces = []
for i in range(seeds.NPoints()):
    tc = tmesh.clone().threshold(name="fragment", above=i-0.1, below=i+0.1)
    mc = tc.tomesh(fill=False).color(i)
    pieces.append(mc)

############### animate
plt = Plotter(size=(1200,800), axes=1)
plt.show(txt, pieces)
for i in range(20):
    for pc in pieces:
        cm = pc.centerOfMass()
        pc.shift(cm/25)
    txt.text(f"{__doc__}\n\nNr. of pieces = {seeds.N()}")
    plt.render()
plt.interactive().close()

Screenshot from 2022-03-08 15-31-56 Screenshot from 2022-03-08 15-32-14

ttsesm commented 2 years ago

This is cool.

Yes except of tetgen and pymeshfix I also found that pygalmesh can be used as well. I do not know if Octree from vtk can also be used instead, I saw some people mentioning about it for such a task.

Now I need to figure out how to get smoother fragmentations, I will play a bit with the parameters to see if I can optimize the cutting part.

Thanks a lot :+1:

ttsesm commented 2 years ago

Hi @marcomusy, is it possible to load .msh tetrahedra files with vedo? Usually used with the gmsh lib.

marcomusy commented 2 years ago

It should - if it doesnt load in vedo you may use meshio package and then pass it to vedo, check out examples/other/pygmsh_cut.py

ttsesm commented 2 years ago

The tmesh = vd.TetMesh("test.msh") directly to the file doesn't work. It gives 2022-05-12 10:43:46.451 ( 189.979s) [ B1CBF740] vtkDataReader.cxx:543 ERR| vtkUnstructuredGridReader (0x5559b098fd90): Unrecognized file type: $MeshFormat for file: test.msh. Then through meshio seems to work fine:

tmesh_meshio = meshio.read("test.msh")
msh = vd.TetMesh([tmesh_meshio.points, np.vstack((cell.data for cell in tmesh_meshio.cells if cell.type=="tetra"))]).tomesh()
ttsesm commented 2 years ago

@marcomusy from the previous example (with the split) I am trying to subdivide and smooth the extracted pieces:

pieces = []
for i in range(seeds.NPoints()):
    tc = tmesh.clone().threshold(name="fragment", above=i-0.1, below=i+0.1)
    mc = tc.tomesh(fill=False).color(i).subdivide(3).smooth()
    pieces.append(mc)

but it gives me the following error:

ERROR:root:Algorithm vtkLoopSubdivisionFilter(0x55b2082dd180) returned failure for request: vtkInformation (0x55b1ff34e980)
2022-05-12 13:47:19.954 (1501.291s) [        16727740]vtkWindowedSincPolyData:108    ERR| vtkWindowedSincPolyDataFilter (0x7f067400c860): No data to smooth!
ERROR:root:No data to smooth!
2022-05-12 13:47:27.410 (1508.747s) [        16727740]vtkSubdivisionFilter.cx:60     ERR| vtkLoopSubdivisionFilter (0x55b2082dd0e0): No data to subdivide
ERROR:root:No data to subdivide
2022-05-12 13:47:27.410 (1508.747s) [        16727740]       vtkExecutive.cxx:752    ERR| vtkCompositeDataPipeline (0x55b2082bffa0): Algorithm vtkLoopSubdivisionFilter(0x55b2082dd0e0) returned failure for request: vtkInformation (0x55b1ff4412a0)
  Debug: Off
  Modified Time: 365532050
  Reference Count: 1
  Registered Events: (none)
  Request: REQUEST_DATA
  FROM_OUTPUT_PORT: 0
  ALGORITHM_AFTER_FORWARD: 1
  FORWARD_DIRECTION: 0
ERROR:root:Algorithm vtkLoopSubdivisionFilter(0x55b2082dd0e0) returned failure for request: vtkInformation (0x55b1ff4412a0)
2022-05-12 13:47:27.410 (1508.747s) [        16727740]vtkWindowedSincPolyData:108    ERR| vtkWindowedSincPolyDataFilter (0x7f067400c860): No data to smooth!
ERROR:root:No data to smooth!

apparently my mc.faces() and mc.cells() lists are empty after calling .tomesh(fill=False): image image However, this shouldn't be happening, right? There should be at least the boundary faces there.

ttsesm commented 2 years ago

Also to speed up the split into pieces instead of creating seeds and iterating over the face centers in order to label them (which as I notice in a mesh with dense tetrahedra is quite time consuming, thus quite slow), do you think it would make sense to do the split as a combination of the intersect and minus boolean operations?

This means though that I would need a really random shaped 3d model to apply the boolean operations. I was thinking something like a 3d dendrogram (or neural tree, I do not really know how these are called) image

which would grow random branches from the center of mass of the 3d mesh until its surface and then we could use this for the boolean operations between the two models? Or it could be something simpler.

marcomusy commented 2 years ago
ttsesm commented 2 years ago
  • uhm i dont understand.. is the rabbit example working for you?

The rabbit example working fine.

  • it doesnt make much sense to subdivide the mesh version of the tetmesh (that is just meant for visualization!), but you can subdivide the tetmesh itself.

Why not? In principle at the end I want to keep the individual mesh version of the pieces. However, because the fragmented sides are kind of too rough I would like to smooth them. One way to do that is to create more tetrahedras or the second that I thought is to just smooth the surface of the meshes that correspond to the meshes. I've tried the first one but once you start increasing the amount of tetrahedras the post-processing becomes to much time consuming. Thus, I thought the solution of smoothing and to apply some smoothness in order to keep the amount of tetrahedras in logical numbers and then fix the roughness by smoothing. At least I wanted to check how it looks like. However, I am getting the error I've mentioned above. Test the code bellow:


"""Add a custom scalar to a TetMesh to segment it.
Press q to make it explode"""
from vedo import show, Mesh, Sphere, TetMesh, Plotter, dataurl, printc, Text2D
import tetgen
import pymeshfix

n = 20000 f1 = 0.005 # control the tetras resolution f2 = 0.15 # control the nr of seeds

repair and tetralize the closed surface

amesh = Mesh(dataurl+'bunny.obj') meshfix = pymeshfix.MeshFix(amesh.points(), amesh.faces()) meshfix.repair() # will make it manifold repaired = Mesh(meshfix.mesh) tet = tetgen.TetGen(repaired.points(), repaired.faces()) tet.tetrahedralize(order=1, mindihedral=210, minratio=1.2, maxvolume=0.01) tmesh = TetMesh(tet.grid)

surf = tmesh.tomesh(fill=False)

pick uniform pts on the surface

seeds = surf.clone().subsample(f2).ps(10).c('black') printc("#pieces:", seeds.N())

assign to each tetrahedron the id of the closest seed point

cids = [] for p in tmesh.cellCenters(): cid = seeds.closestPoint(p, returnPointId=True) cids.append(cid) tmesh.celldata["fragment"] = cids

pieces = [] for i in range(seeds.NPoints()): tc = tmesh.clone().threshold(name="fragment", above=i-0.1, below=i+0.1) mc = tc.tomesh(fill=False).color(i) mc = mc.subdivide(3).smooth() <-------------------------------------- this one give me the error!!!!!!! pieces.append(mc)

############### animate plt = Plotter(size=(1200,800), axes=1) plt.show(doc, pieces) for i in range(20): for pc in pieces: cm = pc.centerOfMass() pc.shift(cm/25) plt.render() plt.interactive().close()



The `mc = mc.subdivide(3).smooth()` command should work since I am just applying some operations on a `Mesh` object, shouldn't be?

> * do not use boolean with complex shapes , it will most probably fail and will be even slower!
>
Ok
marcomusy commented 2 years ago

The mc = mc.subdivide(3).smooth() command should work since I am just applying some operations on a Mesh object, shouldn't be?

Yes, the problem is much more serious .. infact tetgen doesn't seem to work well:

Screenshot 2022-05-14 at 01 10 47

you may notice that the wireframe part is very irregular and non-manifold, that is why the subdivision and smoothing fails.

The vedo tetralizer seems a bit better but still there are a few degenerate tets...:

from vedo import *

amesh = Mesh(dataurl+'bunny.obj')
amesh.fillHoles().cap().smooth()

tmesh = amesh.tetralize(side=0.015, debug=False)
#tmesh.write('mytetmesh.vtk')  # save to disk!

surf = tmesh.tomesh(fill=False)
surf.subdivide(method=1).smooth()

show(surf, axes=1).close()
Screenshot 2022-05-14 at 01 20 33

sorry in this moment it is not clear to me why that happens.

ttsesm commented 2 years ago

Hhhmm I see, I've also used another tetrahedrealizer but the problem remains (though I forced the output to be manifold). Check on the attached file. You can load it as:

tmesh_meshio = meshio.read("bunny.msh")
tmesh = TetMesh([tmesh_meshio.points, np.vstack((cell.data for cell in tmesh_meshio.cells if cell.type=="tetra"))])

bunny.zip