Weird issues when checking for intersection between meshes

[ManuGraiph]

Hi!

I have the following situation/error and i haven't been able to figure out why it happens. I can't post the data, since it's work product.

I have a lot of meshes, all of them being just a planar "square" made by two triangles joined at one side, rotated in axis z and x, and centered in different positions in the space (x,y,z). Now, i'm checking for intersections between those meshes, using two nested loops that iterate pivoting over each one of them and checking for intersections with the rest using the intersectWith() method.

Now, the weird behaviour is the following: Sometimes when checking two meshes for intersections, python crashes, without any error or warning, it just crashes. I tried to debug/guess what was happening and i was able to find out that when two meshes have the EXACT same midpoint position (the x,y,z of the center) the intersectWith() method crashed (no idea why). I tried to quick-fix it by adding noise to the planes center coordinates, which so far fixed it for those cases, but i'm still getting the same error in other planes where apparently there is nothing related to that.

I'm asking here in case it's a known issue; i suspect that something like.. idk, when 2 vertices are exactly the same or something similar happens, the function crashes, but i don't have a better clue about it.

This is an example of two planes that were crashing python when trying to calculate the intersection:

Any help or suggestion is more than welcome!

[ManuGraiph]

I wrote a reproductible example:

m1v = np.array([[17903.93947473, 11418.12995419,  2572.68075251],[17889.43531656, 11423.09300945,  2546.89320749],[17913.65204344, 11446.51421055,  2572.68075251],[17899.14788527, 11451.47726581,  2546.89320749]])
m2v = np.array([[17897.68174692, 11454.30365475,  2573.22453367],[17917.8546173 , 11442.45883457,  2554.44274633],[17882.4917027 , 11428.43354543,  2573.22453367],[17902.66457308, 11416.58872525,  2554.44274633]])

m1f = [[0,1,2],[1,2,3]]

m1 = vd.Mesh([m1v,m1f])
m2 = vd.Mesh([m2v,m1f])

m1.intersectWith(m2,tol=1e-06)

That gives me a python crash.

[marcomusy]

Hi the first thing i notice is that one face you create is reversed. but that doesn't seem to be the origin of the problem, in fact your code above doesn't crash on my mac... although i get a warning message

2022-09-20 21:29:04.675 (   0.783s) [           73938]      vtkDelaunay2D.cxx:1409  WARN| vtkDelaunay2D (0x600003b64000): Edge not recovered, polygon fill suspect
2022-09-20 21:29:04.675 (   0.783s) [           73938]       vtkExecutive.cxx:753    ERR| vtkCompositeDataPipeline (0x600003b61340): Algorithm vtkIntersectionPolyDataFilter(0x600003764540) returned failure for request: vtkInformation (0x600001779fc0)
  Debug: Off
  Modified Time: 1258
  Reference Count: 1
  Registered Events: (none)
  Request: REQUEST_DATA
  FROM_OUTPUT_PORT: 0
  ALGORITHM_AFTER_FORWARD: 1
  FORWARD_DIRECTION: 0

which to be honest i don't quite understand...

try out this version:

import numpy as np
import vedo as vd 

m1v = np.array([[17903.93947473, 11418.12995419,  2572.68075251],[17889.43531656, 11423.09300945,  2546.89320749],[17913.65204344, 11446.51421055,  2572.68075251],[17899.14788527, 11451.47726581,  2546.89320749]])
m2v = np.array([[17897.68174692, 11454.30365475,  2573.22453367],[17917.8546173 , 11442.45883457,  2554.44274633],[17882.4917027 , 11428.43354543,  2573.22453367],[17902.66457308, 11416.58872525,  2554.44274633]])

m1f = [[0,1,2],[3,2,1]]

m1 = vd.Mesh([m1v,m1f]).subdivide(method=1).bc("green5").lw(1)
m2 = vd.Mesh([m2v,m1f]).subdivide(method=1).bc("blue5").lw(1)

inters = m1.intersectWith(m2)
vd.show(m1,m2, inters, axes=1)

[ManuGraiph]

Adding that when creating the mesh did the trick, now it didn't crash, but i think it's just that with that subdivide the verts just don't fall into whatever is causing the issue.

I think it's something with the vertices positions. As i said in my first message, i was having the crash when the mid points were the same and i fixed it by adding noise, and when i wasn't using a seed, the crash happened in random pairs, so i'm guessing it's somithng really psecific related to the verts positions.

[marcomusy]

I also guess it is the subdivide().. checking intersections in a nested loop might become slow (?), you may consider writing a custom intersectWith() function if needed..

[ManuGraiph]

If instead of the subdivide i just apply a m1.shift(dx=0.01,dy=0,dz=0) it also works, it's something 100% related to the node specific positions. The thing is, i can't fix it programatically, cause it makes python crash.

What version of vtk are you using? The issue is actually in the bf.Update() of the intersectWith() function.

--

About the function, of course, the speed is also a big issue in this case. I basically have thousands of those square meshes and i have to check all the intersections between them and get the intersection line when they do.

[ManuGraiph]

What version of vtk and python are you using @marcomusy ?

[marcomusy]

vedo --info

vedo version      : 2022.3.0.dev0   https://vedo.embl.es
vtk version       : 9.0.3
python version    : 3.9.7 (default, Sep 16 2021, 08:50:36) [Clang 10.0.0 ]
python interpreter: /Users/mmusy/Software/miniconda3/bin/python
vedo installation : /Users/mmusy/Software/miniconda3/lib/python3.9/site-packages/vedo
system            : Darwin 21.6.0 posix x86_64
monitor (primary) : None, resolution=(1280, 800), x=0, y=0

[ManuGraiph]

Can't find the 9.0.3 version of vtk anywhere, not on pip, conda, or even the wheel. I'm sure the issue is in vtk.

[marcomusy]

I just tested it on vtk9.2.0 and it crashed python without the subdivide(), and it works fine with it.

if your planes are just 4 coplanar points you may consider solving the problem of the intersection analytically which would increase the speed by a lot.

[ManuGraiph]

Yeah, basically anything that you do to alter those specific points (subdivide, a shift, rounding the verts, etc) make it work.

Could you please elaborate a bit on the equations for that? :O I need to check for intersections and to get the line formed by the intersection when they do. I assumed using meshes would be faster than an analytical method.

[marcomusy]

Yeah, basically anything that you do to alter those specific points (subdivide, a shift, rounding the verts, etc) make it work.

yes... it really looks like a numerical instability in vtk... I'll try to submit an issue there..

Could you please elaborate a bit on the equations for that? :O I need to check for intersections and to get the line formed by the intersection when they do. I assumed using meshes would be faster than an analytical method.

it should be an exercise in linear algebra check out

https://math.stackexchange.com/questions/3328849/intersection-between-two-finite-planes

https://www.kristakingmath.com/blog/parametric-equations-intersection-of-planes#:~:text=The%20intersection%20of%20two%20planes,will%20always%20be%20a%20line.&text=where%20r%200%20r_0%20r,vectors%20of%20the%20two%20planes.

And/or

import numpy as np
import vedo

# crashing
# m1v = np.array([[17903.93947473, 11418.12995419, 2572.68075251],[17889.4353166, 11423.09300945, 2546.89320749],[17913.6520434, 11446.51421055, 2572.68075251],[17899.14788527, 11451.47726581, 2546.89320749]])
# m2v = np.array([[17897.68174692, 11454.30365475, 2573.22453367],[17917.8546173, 11442.45883457, 2554.44274633],[17882.4917027, 11428.43354543, 2573.22453367],[17902.66457308, 11416.58872525, 2554.44274633]])

# not crashing
m1v = np.array([[17903.9, 11418.1, 2572.7],[17889.4, 11423.1, 2546.9],[17913.6, 11446.5, 2572.7],[17899.1, 11451.5, 2546.9]])
m2v = np.array([[17897.7, 11454.3, 2573.2],[17917.8, 11442.5, 2554.4],[17882.5 ,11428.4, 2573.2],[17902.7, 11416.6, 2554.4]])

m1f = [[0,1,2], [3,2,1]]

a = np.cross(m1v[2]-m1v[0], m1v[1]-m1v[0])
b = np.cross(m2v[2]-m2v[0], m2v[1]-m2v[0])

c = np.cross(a,b)
c /= np.linalg.norm(c)

m1 = vedo.Mesh([m1v,m1f]).lw(1).bc("green5")
m2 = vedo.Mesh([m2v,m1f]).lw(1).bc("blue5")

arr1 = vedo.Arrow(m1v[0], m1v[0]+a/50)
arr2 = vedo.Arrow(m2v[0], m2v[0]+b/50)
arr3 = vedo.Arrow(m1v[0], m1v[0]+c*20)

inters = m1.intersectWith(m2)
vedo.show(m1,m2, inters, arr1, arr2, arr3, axes=1)

[ManuGraiph]

Well, yeah, the intersection between two infinite planes is an easy to implement solution, but i'm not so sure if it will be so easy with finite planes, that's why i thought it would be easier and faster with meshes.

I'm gonna give it a try making each "plane" just a triangle (a single face) to see if the problem persists, i'll report back in a bit.

Thanks a lot!

[ManuGraiph]

Well, as a mini-update, i tested using one triangular element per plane, so far i've ran a couple hundred thousand intersection tests and the issue hasn't popped.

Looking at the last picture i posted, i think i figured it out, i think the function crashes when two arists (the lines connecting the vertices) from both meshes intersect. Since the "intersection line" is passing through that point, i think it crashes when it tries to find the intersection points from each face with each other and it turns out two are "the same".

[ManuGraiph]

Ok, not sure of how wise what i'm doing is, but lowering the tolerance to 1*10^-3 or higher makes it work.

What is the tolerance for?

[marcomusy]

OK, thanks for investigating... the tolerance should be to define the how close the algorithm should check i guess, I submitted an issue to the vtk developers as it's definitely a bug in the class.

[ManuGraiph]

It still crashes on some points, but making the tol higher seems to reduce it a lot. When using one-element meshes i haven't had any crash, so i think it's related to the internal arist lines of the mesh..

[ManuGraiph]

Oh, btw, could you please tell me if there is any update from the vtk team about this bug?

[marcomusy]

yes, sure, this is the thread: https://discourse.vtk.org/t/vtkintersectionpolydatafilter-crashing/9428

[ManuGraiph]

Thanks! Just out of curiosity, did you test it with vtk 9.1.0?

[marcomusy]

I tested on vtk9.2.0rc2 but my guess is that the bug is there in all versions.

[ManuGraiph]

Yeah. I'm still convinced it's related to the internal vertex lines intersecting or 'numerically intersecting' (passing too close to each other). For a 2-elements mesh, each one can intersect the other mesh at 2 points, resulting in 4 intersection points, if the mid-arist intersect each other, the 2 'central' intersection points are gonna be too close, and my guess is, if that distance is lower than the tolerance, it crashes. Maybe it assumes both are the same point and then the output doesn't have the expected length or similar.