search

wolph / numpy-stl

Simple library to make working with STL files (and 3D objects in general) fast and easy.
http://numpy-stl.readthedocs.org/
BSD 3-Clause "New" or "Revised" License
628 stars 107 forks source link

Possible error in volume calculations #210

Closed seastate closed 1 year ago

seastate commented 1 year ago

@wolph, thanks for this very useful module.

As far as I can see, the volume calculation is giving incorrect answers. I've pushed a demonstration of the error to the vol_test branch of my fork of your repo (@seastate/numpy-stl).

If you import the test_vol script, it imports the file 'plankter.stl'. This shape is a very small ellipsoid, which appears to import and plot correctly. The script prints the volume and CoG calculated by numpy-stl, the ideal geometrical properties, and the volume and CoG calculated by a matlab/octave script I'm recoding into python.

The descretization is coarse, and the idealized and octave volume and CoG are not exactly equal but close. The numpy-stl volume is very far off (the CoG calculation is harder to evaluate in this simple test, but it seem the volume test is simpler to debug anyway).

It's possible that I'm making a mistake. It's also possible that there is truncation error because the shape is small*, but running on a 64-bit machine I wouldn't expect that.

Thanks again.

System info: Ubuntu Mate 22.04 Python 3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0] on linux

pip show output: Name: numpy-stl Version: 3.0.1 Summary: Library to make reading, writing and modifying both binary and ascii STL files easy. Home-page: https://github.com/WoLpH/numpy-stl/ Author: Rick van Hattem Author-email: Wolph@Wol.ph License: BSD Location: /home/dg/.local/lib/python3.10/site-packages Requires: numpy, python-utils Required-by:

*It's so small because it is input for a hydrodynamic model of microorganism swimming, not for 3D printing :)

seastate commented 1 year ago

P.S. Total areas agree between the two numerical methods, both get 4.911966e-09

seastate commented 1 year ago

This is an issue with inconsistent ordering of the vectors in the stl file I'm using, so that the normals are not consistently pointing outwards. I used a simple algorithm suitable for this simple shape, using sign of the dot product of the vector from the origin to the centroid with the normal vector, to correct vector directions. The numpy-stl calculations are then correct.

I'm working on a more general algorithm that counts intersections of segments with triangles to work on more general shapes. I'll post when I have a candidate fix.

seastate commented 1 year ago

I've pushed a version of test_vol.py with several tests, the last of which (get_mass_props6) uses an intersection-based correction for normal vector direction (enforcing they all point outwards, regardless of vertex order). In this test the corrected algorithm appears to get the correct answer, consistent with the idealized geometry, the octave code, and the results using the simple vector-based scheme for correcting normals direction.

BTW it also uses a simple vector-based algorithm for calculating CoG, volume etc. which to me seems more intuitive and concise than the algebraic algorithm currently in numpy-stl. I'd be interested to hear your perspectives.

wolph commented 1 year ago

@wolph, thanks for this very useful module.

As far as I can see, the volume calculation is giving incorrect answers. I've pushed a demonstration of the error to the vol_test branch of my fork of your repo (@seastate/numpy-stl).

Quite possible to be honest... I'm no linear algebra expert and I honestly don't know in which scenarios the current algorithm might be right or wrong. The only thing I do know is that the results in the tests correspond with the output that meshlab produces. But that assumes that meshlab is correct to begin with ;)

It's possible that I'm making a mistake. It's also possible that there is truncation error because the shape is small*, but running on a 64-bit machine I wouldn't expect that.

Even though your machine is 64-bit, the STL meshes use 32 bit floating point numbers to be compatible with the binary STL format. If you're not using the binary STL format and don't need to have the stl writing speedups available, you could increase that to 64 bit.

I should note however... that is completely not tested and I wouldn't be surprised if it breaks things :P

Your method looks good to me so I wouldn't mind merging it. I'm guessing that each method has some upsides/downsides in specific cases (i.e. performance vs accuracy) so it might be a good idea to support multiple methods anyhow