marcomusy / vedo

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

Pointdata and celldata lost when decimating mesh #992

Open JeffreyWardman opened 9 months ago

JeffreyWardman commented 9 months ago
import numpy as np
import vedo

mesh = vedo.load("")
mesh.pointdata["counter"] = np.arange(mesh.npoints)
mesh.decimate(n=20000)

assert len(mesh.pointdata.keys()) != 0

The above loses all cell and point data information. The trick mentioned in #571 and #620 therefore does not work. I'd like to have the point data resolved at least for my use case.

marcomusy commented 9 months ago

you can use

mesh.decimate(n=20000, method="pro")

Another option is to interpolate back the data with

decimated_mesh.interpolate_data_from(original_mesh, n=3)

ps: i edited the examples in the links.

JeffreyWardman commented 9 months ago

Unfortunately neither solution is suitable for me. The former really struggles with curves and creates a lot of triangles that look like dents. The latter takes far too long (to the point where I haven't seen it complete after a couple of minutes). Is there a way of including the pointdata labels with the default method?

marcomusy commented 9 months ago

How big is the mesh? the interpolate_data seems to do it fast on 500k points mesh. Anyways I have updated the decimate method and split ii in 3 different methods:

from vedo import *

m1 = Mesh(dataurl + "dragon.ply")
m1.lw(1).lighting('off')
m1.generate_random_data()
# print(m1)

print("..decimating")
m2 = m1.clone()

m2.decimate(0.25)
# m2.decimate_pro(0.25)
# m2.decimate_binned([200,200,200], use_clustering=False)
# m2.decimate_binned([200,200,200], use_clustering=True)

# print("..interpolate_data_from")
# m2.interpolate_data_from(m1, n=3)
print(m2)
printc("is_manifold:", m2.is_manifold(), invert=1)

show(m1, m2, N=2, axes=1, zoom=4)

Screenshot from 2023-12-19 20-56-55

unfortunately i just realized that - at least in this specific mesh - all the filters seem to produce a non-manifold output.. but I cannot fix that, it's a problem of the upstream vtk apparently..

JeffreyWardman commented 9 months ago

Do you have a link for the issue in VTK's GItLab repository?

Is there a way to allow decimate to return the removed/remaining point indices? That would resolve the issue. Not sure why interpolate_data_from wasn't working before but using it produces poor results (expect it to be green, the other colour is another label; it might be treating the labels as floats then rounding rather then as integers):

image (Ignore the border between black and green having a bunch of colours. That's a different issue). The main issue are the orange and aqua being in the green).

I found a bug where the vertex positions change somewhat due to converting from float64 to float32. Can be verified via:

indices = np.where(np.all(original_vertices[:, None, :] == new_vertices, axis=-1))
# vs
float32_indices = np.where(np.all(original_vertices.astype(np.float32)[:, None, :] == new_vertices, axis=-1))
raphaelsulzer commented 8 months ago

Hi,

I have a similar issue:

I have a mesh with custom data per face, stored in mesh.celldata["mydata"] I now iteratively collapse edges of the mesh and call mesh.clean() after each iteration to update topology (a slightly modified version of vedo.Mesh.collapse_edges()). While mesh.clean() correctly removes faces, e.g. when a triangle edge is collapsed, it does not update the celldata array. Is there any way to update the celldata, or e.g. get a map from input to output cells when calling mesh.clean()?

NB: the trick to move "mydata" to pointdata before mesh.clean() and then back does not work, because I do not want to interpolate "mydata".

Louis-Pujol commented 8 months ago

Hi all !

I worked a few time ago on celldata preservation with quadric decimation. Regarding the algorithm, it seemed feasible to keep track of successive edge collapses. Unfortunately, vtk implementation did not keep track of this information, so after calling vtkQuadricDecimation, you can't find a correspondence between vertices in fine and coarse representation.

I uploaded my work in the package fast-simplifcation, it is a python wrapper around the C++ implementation fast-quadric-mesh_simplification (fqmr), an approximate and faster version of vtk quadric decimation. An advantage of fqmr over vtkQuadricDecimation is that the implementation is contained in one file, much more easy to customize than the vtk implementation with all the dependencies. The modification I did roughly consisted in storing the indices of vertices at each collapse, so it came with almost no time overhead compared to just call the decimation function.

Preservation of CellData is not accessible in one line, but you have access to an array called indice_mapping that gives you the map from the vertices of the fine mesh to the vertices of the coarse mesh. Then, using a scatter operation (mean in the following example), you can transfer signal from fine vertices to coarse vertices. I think you can do something similar for faces with the information contained in indice_mapping.

If there is a shared interest in this feature, it could be interesting to write it cleanly somewhere to make it accessible to the broader community.

Here is an example:

import fast_simplification
import vedo
import numpy as np

# Load a mesh and define a signal on vertices
mesh = vedo.Sphere()
points = mesh.vertices
faces = mesh.cells
signal = points[:, 2]
mesh.pointdata["signal"] = signal

# Decimate the mesh and compute the mapping between the original vertices
# and the decimated ones with fast-simplification
points_decim, faces_decim, collapses = fast_simplification.simplify(
    points,
    faces,
    0.9,
    return_collapses=True
    )

points_decim, faces_decim, indice_mapping = fast_simplification.replay_simplification(
    points,
    faces,
    collapses
)

# Compute the average of the signal on the decimated vertices (scatter operation)
import numpy as np
unique_values, counts = np.unique(indice_mapping, return_counts=True)
a = np.zeros(len(unique_values), dtype=signal.dtype)
np.add.at(a, indice_mapping, signal) # scatter addition
a /= counts # divide by the counts of each vertex indice to get the average

# Create a new mesh with the decimated vertices and the averaged signal
decimated_mesh = vedo.Mesh([points_decim, faces_decim])
decimated_mesh.pointdata["signal"] = a

vedo.show(mesh, decimated_mesh, N=2, axes=1)

And a screenshot of the output: Capture d’écran du 2024-01-13 15-59-39

marcomusy commented 8 months ago

This looks fantastic, thanks a lot Louis for the contribution! I added your script as an example here https://github.com/marcomusy/vedo/blob/master/examples/other/fast_simpl.py

related: #1003 #1007

Louis-Pujol commented 8 months ago

Thanks @marcomusy. Just a quick remark: in the docstring of the example it is indicated that data is transferred from faces to faces, it is actually transferred from vertices to vertices. It could be possible to find the faces/faces correspondence, I didn't do that the time I worked on fast simplification as it wasn't my use case. I keep it in mind for future updates.

For the moment, if you want to transfer data from faces to faces, I would recommend to add a small tweak to transfer data from faces to vertices and vice-versa.

marcomusy commented 8 months ago

@Louis-Pujol Would you mind opening PR to correct the doc string?

@raphaelsulzer I hope this can help addressing your case!

raphaelsulzer commented 8 months ago

@raphaelsulzer I hope this can help addressing your case!

As you mentioned in the previous issue #1007, the newest version of vedo is able to preserve/update point and celldata when doing an edge collapse. This is exactly what I needed.