Splitting connected components and getting the perimeter

[ManuGraiph]

Hi Marco,

I wanted to ask you a couple questions..

1. If i have a planar mesh like the one in the picture:

is it possible to "split" it into those two separated, connected, regions, basically like splitting them by a 'connected components' logic?

2. Let's say i have one of those two connected areas (or any other mesh), can i get the perimeter of it somehow (the external perimeter, the hull)?

3. Random bonus question: Would it be possible to add native support for torch in parallel with numpy in vedo? The torch.tensor subclass is a 99.9% drop-in replacement for numpy, with the plus that you can allocate it on GPU and speed up operations by a lot.

Thanks in advance!

[marcomusy]

Hi @ManuGraiph consider the following example:

from vedo import *

c1 = Circle()
c2 = Circle().scale([2,1,1]).rotate_x(20).shift(4,1,0)

l1 = Line(c1.boundaries()).c("blue").lw(4)
l2 = Line(c2.boundaries()).c("purple").lw(4)

print("lengths", l1.length(), l2.length())

# merge two circles into one mesh and split them back
circles = merge(c1, c2, flag=False)
circles.split(flag=True)
print(circles)

# merge two circles into one mesh and split them back
circles = merge(c1, c2, flag=True)
c3, c4 = circles.split(flag=False)
print(circles)

show(c3, c4, l1, l2, axes=1, viewup="z")

check out: https://vedo.embl.es/autodocs/content/vedo/vedo/mesh.html#Mesh.split https://vedo.embl.es/autodocs/content/vedo/vedo/pointcloud.html#merge

[ManuGraiph]

Every day i'm more and more convinced this is the best python library out there.

Thanks a real lot!!

[marcomusy]

Thanks Manu for your kind words, your user feedback and contribution helps constantly improving it!

[ManuGraiph]

Thanks to you for such an excellent work! If i had more knowledge on vtk i'd love to contrib in this project.

[ManuGraiph]

Sorry to reopen, a couple more questions arised and didn't wanna open a new post.

1. Is there any way to "filter" a mesh and keep only the cells that have a certain property in one of their celldata entries? Let's say i have a numerical value in a celldata property called 'x', can i, somehow, keep only the cells that have a value of.. 5 in that attrib?

2. Is there any special way to use the boolean method? I'm trying to get the intersection & union of 2 meshes (basically to get an iou) and for some reason they keep giving me an area of 0 (-> mesh1.boolean("plus", mesh2).area()).

3. Regarding the boundaries question, how can i get the POINTS from the boundaries but in order? If i get the boundaries.points or Line.points (the same) they are not in "order", like, they are not in the order of the contour.

Thanks a lot in advance!

[marcomusy]

Hi @ManuGraiph

1. Yes check out example vedo -r examples/basic/mesh_threshold.py and this variant

from vedo import *

man = Mesh(dataurl+"man.vtk")

# a dummy scalar field is added to the mesh
cellcents = man.cell_centers().coordinates  # get the cell centers
man.celldata["myscalars"] = cellcents[:, 0] + 37

man.cmap("coolwarm", "myscalars").add_scalarbar(horizontal=1)

# make a copy and threshold the mesh
cutman = man.clone().threshold("myscalars", 37, 37.5, on="cells").lw(1)

# distribute the meshes on the 2 renderers
show([man, cutman], N=2, elevation=-30, axes=11).close()

2. Is the output as exapected? Are you computing the normals? Have you tried to use triangulate()?

3. Have a look at this. Hope it helps.

[ManuGraiph]

I'll check it tomorrow, thx a lot again!

[ManuGraiph]

Ok, threshold worked perfect, thanks, literally what i needed. (2) and (3) point to the same goal, so, i'll stick to (2).

I have 2 "horizontal" meshes that cover the same area (let's call em vi and vf), but they are different, the first is the one shown above, the other one is a "reconstruction" based on some technique we are testing. Now, my goal is to calculate the iou between those two.

I set the z coordinate of both meshes points to 0, to get sure they are on the same plane, so, i have:

Now, when i do:

vi.boolean("intersect",vf).area() vi.boolean("plus",vf).area()

Both are 0, they says "there is no intersection". I'm not sure if this can be done for planar meshes or if must be done for "closed-volume" meshes.

[marcomusy]

Try

from vedo import *

square1 = Rectangle().scale(2).c('r5').alpha(0.5)
square1.triangulate().subdivide(5, method=1).lw(1)

square2 = Rectangle().scale(1).rotate_z(25).shift(.1,.2).c('b5').alpha(0.5)
square2.triangulate().subdivide(5, method=1).lw(1)

square1.distance_to(square2)
print(square1)
cut_square = square1.clone().cut_with_scalar(0.001, "Distance")

show([[square1, square2], cut_square], N=2, axes=1)

Other possibilities: https://vedo.embl.es/autodocs/content/vedo/vedo/pointcloud.html#Points.cut_with_line https://vedo.embl.es/autodocs/content/vedo/vedo/pointcloud.html#Points.cut_with_cookiecutter https://vedo.embl.es/autodocs/content/vedo/vedo/pointcloud.html#Points.cut_with_point_loop

Boolean indeed only works for closed surfaces in 3d!

[ManuGraiph]

Ok, worked with that approach, there is an obvious trade off between the cpu time and the subdivisions, but for this specific purpose it's acceptable, thanks!.

I'm gonna benchmark it with the other approach, using the join_segments fxn and getting the intersection and union manually.

Thanks a real lot for the help Marco!!

[ManuGraiph]

I'm struggling a bit with the join_segments method, my understanding is that i should be able to do something like this (correct me if i'm wrong):

a = mesh.boundaries() b = Line(a) c = b.join_segments()

In order to get the border points of a mesh, in the correct order, however, my python gets stuck when i use the join_segments method, am i doing something wrong?

[marcomusy]

Hi, how many points has got the line? if you want you can attach a file so i can have a look. mesh.boundaries().write("test.ply")

[ManuGraiph]

Sure, it has +-9300.

bounds.tar.gz

[marcomusy]

Sorry, my mistake: I need the vtk file as ply does not save lines.. mesh.boundaries().write("test.vtk")

[marcomusy]

Oh maybe it's just this:

from vedo import *

s = Rectangle().subdivide(5, method=1)
s.cut_with_sphere([0.5,0.5,0], 0.5).scale([1,1,1])
s.subdivide(5, method=1)
b = s.boundaries()

line = Line(b)
print(line)

segs = b.join().join_segments()

print(segs[0].length())
show(s,b, axes=8)

[ManuGraiph]

There it is in vtk format, thanks!

bounds_vtk.tar.gz

[marcomusy]

Ok, try:

from vedo import *

b = load("test.vtk")
line = Line(b)
print(line)

segs = b.join().join_segments()
print(segs[0])

show(segs, b, axes=8)

[ManuGraiph]

Ok, i found the issue, i was not using the .join() method, which gets it stuck, with the .join() in-between works perfect!!!

I'm gonna try to calculate the iou using shapely instead of the cut we talked about last week to see if i can improve the speed.

Thanks a lot for the time!!

[marcomusy]

What about this solution:

from vedo import *
import time

square1 = Mesh(dataurl+"dolfin_fine.vtk")
square1.subdivide(method=1).lw(1)

square2 = Rectangle().scale(0.9).rotate_z(25).shift(0.1,0.3)
square2.triangulate().subdivide(5, method=1)
square2.c('blue5').alpha(0.5).lw(1)

t0 = time.time()
a1 = square1.area()
a2 = square2.area()

esq = square2.clone().z(-1).extrude(2)
square1.cut_with_mesh(esq)
a3 = square1.area()

# compute intersection over union area
iou = a3 / (a1 + a2 - a3)

t1 = time.time()
print(f"IOU  : {iou:.4f}")
print(f"Area1: {a1:.4f}, Area2: {a2:.4f}, Area3: {a3:.4f}")
print(f"Time : {t1-t0:.4f} seconds")
show(square1, square2, axes=8)

IOU  : 0.4035
Area1: 0.9027, Area2: 0.8100, Area3: 0.4924
Time : 0.0855 seconds

[ManuGraiph]

I tried some approach with extrude a couple days ago, basically to use boolean. I'm gonna test this one!

The one i'm currently using is:

    from shapely.geometry import Polygon

    vip = vi.boundaries().join().join_segments()[0].points[:, 0:2]
    vfp = vf.boundaries().join().join_segments()[0].points[:, 0:2]

    p1 = Polygon(vip)
    p2 = Polygon(vfp)

    iou = p1.intersection(p2).area / p1.union(p2).area

Vi and vf are the two meshes i'm trying to cross.

EDIT: Tested, i'm getting the same time with both methods, but the last one you posted is giving me incorrect results for the ious

[ManuGraiph]

Hi,

I ran into a little issue :(.. I guess where it's coming from, but i'm not sure if it can be fixed.

I have a mesh, which, after getting their boundaries and doing a .join() on them, i get this (which is correct):

But.. when i do a join_segments() on it, i get only the lower element:

My guess is that it's happening cause the element is connected to the "larger" part in a single point, so the segments somehow get dizzy.

Is there any way to modify this behaviour and to get the full connected region?

Thanks in advance!

[marcomusy]

Hi , sorry I was away! Have you already solved this?

[ManuGraiph]

Hi Marco! No worries :)!

Nop, i haven't been able to solve it.. When i have situations like the above (or an intersection composed of 2 separated bodies), the perimeter gets dizzy.

Actually, i don't know how to move forward, cause with 2 separate shapes, the concept of IoU doesn't really makes sense anymore :/.

[marcomusy]

I see. If I understand it correctly the problem is underdetermined in general. One way of "injecting more information" in the system could be to subdivide the line - which is currently comprised of only corner points - by adding a sufficient nr on intermediate points, so that at crossings it becomes better defined what the direction of the loop is?

[ManuGraiph]

Well, yes and no, the scope of what i'm doing is to see "how well" the shape/area/etc of a planar body is reconstructed/recovered by probing info in some specific points of a fixed grid, hence, i can't "interpolate" between them cause it will pretty much defy the point of it.

I think i'll just flag those cases and analyze them manually.