Rewrite of is_closed check.

[jb-leger]

This PR contains to (very linked) changes:

• The rewriting of the is_closed check. (tl;dr; the computation is done using float64 to avoid errors on floating point aritmetics, and the maximum error is now computed relatively to the absolute value (since 1e-4 was not appropriate in all cases)
• Adding a exact is_closed check. (tl;dr; assert each edge is present exactly two times, one time in each direction (an edge can be two time in the same direction).

Rewriting of is_closed check

I had problems, false negative, with a mesh (which is closed, I am sure of that property):

In [3]: m.update_normals()

In [4]: m.is_closed()
Your mesh is not closed, the mass methods will not function
            correctly on this mesh.  For more info:
            https://github.com/WoLpH/numpy-stl/issues/69
Out[4]: False

In [5]: m.normals.sum(axis=0)
Out[5]: array([-7.1411133e-03,  1.3427734e-03, -1.7878906e+02], dtype=float32)

In [6]: m.normals[:,2].sum()
Out[6]: 0.0234375

In [7]: m.data.size
Out[7]: 577544

The problem is floating point arithmetic on 32 bits with a large number of triangles. If we do the sum on float64:

In [8]: np.asarray(m.normals, dtype=np.float64).sum(axis=0)
Out[8]: array([ 5.96046448e-08, -6.82390237e-07, -8.84249806e-04])

However, we can quantify the maximum allowed error by using a pessimistic approach. Each initial value (in float32) have a error of abs(value)*eps where eps is the relative error of float32, then the maximum error is sum(abs(value))*eps:

In [9]: np.abs(np.asarray(m.normals, dtype=np.float64)).sum(axis=0) * np.finfo(np.float32).eps
Out[9]: array([0.00313204, 0.0031321 , 0.01907349])

Then we obtain a is_closed` test without (I hope) false negative.

Remark: the problem appear, even with a low number (like 1e4) of triangles.

Exact version of is_closed

As shown in #198, the is_closed can lead false positive. Because the sum to zero of the normals is a necessary condition, but not a sufficient one.

The #198 ask to check "if each edge is defined exactly twice", but again, this is a necessary condition, but not a sufficient one (even with adding the previous check with the normal).

The necessary and sufficient condition is the following: each undirected edge must be defined twice, each directed edge must be defined exactly once. The direction used is counterclockwise w.r.t. the normal (triangle which have a opposite normal w.r.t. the cross product are reversed).

This check is implemented with exact=True. And success to detect example of #198 as unclosed.

However, this check is quite longer than the not exact check on normal, therefore the not exact test is conserved by default, and this check is added as a option.

Remark: fixes #198.

[wolph]

This looks like a great addition, thank you! I'm testing it now and plan to create a new release this week :)

[jb-leger]

FYI, I am working on another check I need, the closeness is not sufficient for my usage and I am curently implementing a self intersection test. If you want to include that in a new release, it should be a good idea to wait 2-3 days.

[wolph]

I'm planning to do a little maintenance anyhow so you've got time :)

[wolph]

Just wondering, are you close to finishing it yet? :)

[jb-leger]

Apologies for the delay. No, I don't have a functional optimized self-intersecting test. (I have one which is very time consuming). I will open another PR in the future if I success. Do not wait for this test for a new release (2-3 months is more likely than 2-3 days). Sorry for the hope.

[wolph]

No worries, all help is always appreciated!

I wouldn't mind adding an extra slow but accurate method either as an option. If you want to include it that is :)

[wolph]

Github doesn't appear to understand the merge, but it's been merged and released.

Thank you for the help!