isSelfIntersecting function gives apparently incorrect result

I'm trying to use the isSelfIntersecting function from the Lidar Toolbox to check whether or not a surfaceMesh object is self intersecting. In most cases it works fine, but in certain situations it returns true when the surface is not self intersecting.

The problem seems to arise in situations where all vertices in the mesh are contained, or nearly contained, by one plane (as in the example below). I wondered if it might be an issue with leaving the face normals undefined, but I get the same issue even if I define them when calling surfaceMesh (not included in example for readability).

Is this an issue with the way I'm defining the mesh and using the function? Or is this a known issue with isSelfIntersecting?

Many thanks for any help!

% define the vertices for the mesh
vertices = [0	0.0785926524819716	0.0150000000000000; ...
    0.00223394141607314	0.0787352254127876	0.0140425952680485; ...
    0.00428764473102951	0.0788687035475108	0.0131624353567936; ...
    0	0.0822253462122803	0.0150000000000000; ...
    0.00224235572430482	0.0822979726488426	0.0140389888997977; ...
    0.00430286261362254	0.0823659662099516	0.0131559129739490; ...
    0	0.0860000000000000	0.0150000000000000; ...
    0.00225109241982507	0.0860000000000000	0.0140352443565181; ...
    0.00431866355685131	0.0860000000000000	0.0131491406913786];

%define the faces for the mesh
faces = [1	2	5;
    1	4	5;
    2	3	6;
    2	5	6;
    4	5	8;
    4	7	8;
    5	6	9;
    5	8	9];

%create the surface mesh object
mesh = surfaceMesh(vertices,faces);

%plot the mesh and verify visually that it does not self intersect
%(it should be an inclined plane sloping downward along the x-axis)
patch('vertices',vertices,'faces',faces,'EdgeColor',[1 .5 0],'FaceColor',[.5 1 .5]);
view(3)

%call isSelfintersecting, returns true despite mesh appearing to be
%non-self intersecting
isSelfIntersecting(mesh)

Answer 1

[Deep]

Hi Ed,

No, you have defined the mesh correctly and are using isSelfIntersection just fine. This is likely due to floating point errors. On inputting your vertices data, it appears slightly differently in MATLAB workspace variables due to the precision limitations of double.

If you look closely, you can see some differences from your original input. For example, in the 2nd column 1st row, there's a diff of 4e-16. Unfortunately, this seems to be a limitation with floating point computation. Hope this helps!

[Ed Wheatcroft]

Hi Deep,

Thanks for your answer, and for finding this discrepancy.

I tried reducing the precision of vertices using:

decPlaces = 10
vertices = round(vertices,decPlaces);

That fixed the issue so long as I kept decPlaces at or below 10, so this would seem to confirm your theory.

I'm a bit confused as to how this came about in the first place - the mesh in the example is a very small piece of a much larger mesh which exhibited this same issue, which was in turn generated by another MATLAB script. I pasted this piece directly out of the workspace to give a minumum working example, so I'm wondering where the extra decimal place of precison appeared from in vertices(1,2) for example.

[Deep]

I am not aware why that would happen. This might require more information about your system setup and these scripts. Perhaps the output was copied from a computer with a different hardware platform? You may create a new question with more information about these scripts so it can be identified why this is happening.

[Ed Wheatcroft]

Ok, no worries :) I actually found another configuration which gives an apparently erroneous result while I was playing around with this. This time, the mesh is clearly self intersecting, but isSelfIntersecting returns false.

In this new example the entries in vertices should be exactly representable using double (or indeed single) precision, so I wonder if this is even a separate issue.

Either way, it seems like isSelfIntersecting isn't totally reliable.

vertices = [0 0 0;      
    1 0 0;              
    0 1 0;              
    0 0 1;              
    0.5 0.5 1;          
    0.5 0.5 -1];        

faces = [1 2 3;
    1 2 4;
    1 4 5;
    1 5 6];

mesh = surfaceMesh(vertices,faces);

patch('vertices',vertices,'faces',faces,'EdgeColor',[1 .5 0],'FaceColor',[.5 1 .5]);
view(3)

isSelfIntersecting(mesh)