marcomusy / vedo

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

Fill between two non-overlapping boundaries #1043

Closed JeffreyWardman closed 4 months ago

JeffreyWardman commented 4 months ago

I'd like to enclose the bottom of in-between the two meshes seen here:

image

I can get the bottom planes and fill the meshes both entirely but I cannot figure out how to subtract one flat plane from the other.

Is there a simple way to subtract the inner filled boundary from the outer such that only the salmon pink component remains below? image

Note:

image

marcomusy commented 4 months ago

have you tried any of the cut_with_*() methods?

JeffreyWardman commented 4 months ago

Ah yes, I tried cut_with_mesh and cut_with_cookiecutter but without any success. The issue is probably around there not being any points inside although I've tried merging the internal with the triangulated external before cut_with_* without success.

boundaries.zip

boundaries_reordered_pts.zip

marcomusy commented 4 months ago

do you need to keep the original points fixed or can resample? if so you may use Spline() or KSpline() through the original points then join_with_strips()

JeffreyWardman commented 4 months ago

For my use case I need to keep the original points. What I really want to find is the internal area on the flat plane inside the walls of the mesh.

marcomusy commented 4 months ago

i must say i'm a bit confused there seems to be something strange in the dataset

from vedo import *

b0 = Mesh("boundaries_reordered_pts/external_boundaries.vtk")
b1 = Mesh("boundaries_reordered_pts/internal_boundaries.vtk")

# a0 = b0.clone().triangulate().area()
# a1 = b1.clone().triangulate().area()
# print(f"area={a0-a1}")

b0.join()
b1.join()

print(b0)
print(b1)

bm = merge(b0, b1)
bm.triangulate()
bm.lw(1).alpha(0.5)
print(bm)

show(bm, axes=8).close()

a0 = b0.clone().triangulate().area()
a1 = b1.clone().triangulate().area()
print(f"area={a0-a1}")
Screenshot 2024-02-07 at 12 59 09

and i get a triangulation error from vtk

2024-02-07 13:02:28.156 (   0.756s)
 [          2884C4]vtkContourTriangulator.:80     ERR| vtkContourTriangulator (0x6000038816c0): 
Triangulation failed, output might have holes.

but not clear to me why.. plus the vtkContourTriangulator filter is normally very robust..

JeffreyWardman commented 4 months ago

I'm not too sure why. The distribution of points are heavily concentrated in areas. I'm using meshlib's thickenMesh (with smoothing) to create the internal wall (see https://github.com/MeshInspector/MeshLib/issues/2162#issuecomment-1926822536). I can't see anything irregular. I reordered the vertices myself in the example I shared with you based off of nearest distance in an approach very similar to #1019.

Using the originally unordered vertices results in the below. So it might be related to the ordering?

Screenshot 2024-02-08 at 7 54 40 am
JeffreyWardman commented 4 months ago

By the way, if this is out of scope of vedo then please let me know.

j042 commented 4 months ago

444.0692336040022

I think ordering vertices based on nearest distance is a good approach, but I don't think it guarantees an untangled polygon. However, it didn't seem to be an issue here. Vertices in your files weren't coplanar and my guess is that vtkContourTriangulator is not a big fan of nonplanar vertices.

import vedo

e = vedo.load("boundaries_reordered_pts/external_boundaries.vtk")
i = vedo.load("boundaries_reordered_pts/internal_boundaries.vtk")

# enforce coplanarity
e.vertices[:, -1] = 0.
i.vertices[:, -1] = 0.

print(e.triangulate().area() - i.triangulate().area())

As external and internal lines have different winding orientations (outer - ccw, inner - cw), merging those coplanar objects should give you in-between mesh.

e = vedo.load("boundaries_reordered_pts/external_boundaries.vtk")
i = vedo.load("boundaries_reordered_pts/internal_boundaries.vtk")
e.vertices[:, -1] = 0.
i.vertices[:, -1] = 0.

ei = vedo.merge(e, i)

print(ei.triangulate().area())
# 444.0694090021098

This value is close to in-between mesh created with triangles - 444.0692335663805.

marcomusy commented 4 months ago

By the way, if this is out of scope of vedo then please let me know.

it's on its boundaries :)) With Jaewook's suggestion this seems to work:

from vedo import *

e = Mesh("/Users/mmusy/Downloads/boundaries_reordered_pts/external_boundaries.vtk")
i = Mesh("/Users/mmusy/Downloads/boundaries_reordered_pts/internal_boundaries.vtk")

# enforce coplanarity
e.vertices[:, -1] = 0.
i.vertices[:, -1] = 0.

m = merge(e.join_segments(), i.join_segments())
m.triangulate().lw(1).c("blue9")

show(m, f"Area={m.area()}", axes=8).close()
Screenshot 2024-02-08 at 13 07 54
JeffreyWardman commented 4 months ago

@j042 Wow. Thank you so, so much! Such a simple and elegant solution. Thanks both you and @marcomusy for your help. That worked perfectly!

I'll rename the issue to be a bit clearer and close the issue.

marcomusy commented 4 months ago

..don't forget to put back the z values before computing the areas if that is what you need!