Repository: mikedh/trimesh

Issue:
trimesh.load improperly reads multi-body ASCII STLs
Hello!

Really great tool you have written! I have found a rather simple issue reading ASCII STL which include more than one solid body. Attached is an ASCII STL file which includes two solids, denoted as the lines between "solid bodyA" and "endsolid bodyA".

Below is an image from paraview, which shows the two separate STL triangle groups (solids) in the file.

![image](https://user-images.githubusercontent.com/34963757/99145238-5720f300-266d-11eb-9679-f4393ae46516.png)

trimesh however seems to only read the first one:
```

trimesh.load('multiBodyTest.stl', file_type='stl')
for facet in mesh.facets:
    mesh.visual.face_colors[facet] = trimesh.visual.random_color()
mesh.show()
```
![image](https://user-images.githubusercontent.com/34963757/99145292-de6e6680-266d-11eb-97d1-361f83b71768.png)

I'm not sure what the intended behavior is (i.e. feature vs. bug). Would be nice to read the file in and somehow access the individual solids, or at minimum, ignore the distinction between solids and simply read the whole file in ignoring those lines.
[multiBodyTest.stl.gz](https://github.com/mikedh/trimesh/files/5540506/multiBodyTest.stl.gz)



Discussion:

Patch:
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 0267b20f..1dca390d 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -73,7 +73,9 @@ jobs:
         sudo apt-get install -qq -y libspatialindex-dev libgeos-dev
     - name: Install Brew On Mac
       if: matrix.os == 'macos-latest'
-      run: brew install spatialindex geos
+      run: |
+        brew install spatialindex geos
+        brew cask install openscad
     - name: Install Pytest
       run: |
         python -m pip install --upgrade pip
diff --git a/README.md b/README.md
index f766dd06..87e4b695 100644
--- a/README.md
+++ b/README.md
@@ -57,6 +57,13 @@ mesh = trimesh.Trimesh(vertices=[[0, 0, 0], [0, 0, 1], [0, 1, 0]],
                        faces=[[0, 1, 2]],
                        process=False)
 
+# some formats represent multiple meshes with multiple instances
+# the loader tries to return the datatype which makes the most sense
+# which will for scene-like files will return a `trimesh.Scene` object.
+# if you *always* want a straight `trimesh.Trimesh` you can ask the
+# loader to "force" the result into a mesh through concatenation
+mesh = trimesh.load('models/CesiumMilkTruck.glb', force='mesh')
+
 # mesh objects can be loaded from a file name or from a buffer
 # you can pass any of the kwargs for the `Trimesh` constructor
 # to `trimesh.load`, including `process=False` if you would like
diff --git a/models/cubevt.obj b/models/cubevt.obj
new file mode 100644
index 00000000..283013ee
--- /dev/null
+++ b/models/cubevt.obj
@@ -0,0 +1,61 @@
+o Box
+v 0.5 0.5 0.5
+v 0.5 0.5 -0.5
+v 0.5 -0.5 0.5
+v 0.5 -0.5 -0.5
+v -0.5 0.5 -0.5
+v -0.5 0.5 0.5
+v -0.5 -0.5 -0.5
+v -0.5 -0.5 0.5
+v -0.5 0.5 -0.5
+v 0.5 0.5 -0.5
+v -0.5 0.5 0.5
+v 0.5 0.5 0.5
+v -0.5 -0.5 0.5
+v 0.5 -0.5 0.5
+v -0.5 -0.5 -0.5
+v 0.5 -0.5 -0.5
+v -0.5 0.5 0.5
+v 0.5 0.5 0.5
+v -0.5 -0.5 0.5
+v 0.5 -0.5 0.5
+v 0.5 0.5 -0.5
+v -0.5 0.5 -0.5
+v 0.5 -0.5 -0.5
+v -0.5 -0.5 -0.5
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+f 1/1/1 3/3/3 2/2/2
+f 3/3/3 4/4/4 2/2/2
+f 5/5/5 7/7/7 6/6/6
+f 7/7/7 8/8/8 6/6/6
+f 9/9/9 11/11/11 10/10/10
+f 11/11/11 12/12/12 10/10/10
+f 13/13/13 15/15/15 14/14/14
+f 15/15/15 16/16/16 14/14/14
+f 17/17/17 19/19/19 18/18/18
+f 19/19/19 20/20/20 18/18/18
+f 21/21/21 23/23/23 22/22/22
+f 23/23/23 24/24/24 22/22/22
diff --git a/models/multibody.stl b/models/multibody.stl
new file mode 100644
index 00000000..0e9b435c
--- /dev/null
+++ b/models/multibody.stl
@@ -0,0 +1,228 @@
+solid bodyB
+  facet normal -1 -1.34230289E-014 1.65623108E-015
+    outer loop
+      vertex -0.440789347 -0.528135696 0.0626728809
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.440789347 -0.713575677 0.0680076416
+    endloop
+  endfacet
+  facet normal -1 -2.72317242E-014 2.69443431E-014
+    outer loop
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.440789347 -0.528135696 0.0626728809
+      vertex -0.440789347 -0.529415204 -0.0519321219
+    endloop
+  endfacet
+  facet normal -1 2.26640913E-007 8.91020051E-007
+    outer loop
+      vertex -0.510789397 -0.528032778 0.0626652343
+      vertex -0.510789533 -0.718809786 -0.0414456238
+      vertex -0.510789435 -0.713575676 0.068007641
+    endloop
+  endfacet
+  facet normal -1 -1.56375913E-007 1.5928758E-006
+    outer loop
+      vertex -0.510789533 -0.718809786 -0.0414456238
+      vertex -0.510789397 -0.528032778 0.0626652343
+      vertex -0.51078958 -0.529730896 -0.0519249498
+    endloop
+  endfacet
+  facet normal 9.28622855E-009 -0.0287562293 -0.999586454
+    outer loop
+      vertex -0.510789435 -0.713575676 0.068007641
+      vertex -0.440789347 -0.528135696 0.0626728809
+      vertex -0.440789347 -0.713575677 0.0680076416
+    endloop
+  endfacet
+  facet normal 6.68754162E-005 -0.028781449 -0.999585726
+    outer loop
+      vertex -0.440789347 -0.528135696 0.0626728809
+      vertex -0.510789435 -0.713575676 0.068007641
+      vertex -0.510789397 -0.528032778 0.0626652343
+    endloop
+  endfacet
+  facet normal 6.82196909E-013 0.0553381089 0.998467673
+    outer loop
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.51078958 -0.529730896 -0.0519249498
+      vertex -0.510789533 -0.718809786 -0.0414456238
+    endloop
+  endfacet
+  facet normal -0.00014702111 0.0552838452 0.998470668
+    outer loop
+      vertex -0.51078958 -0.529730896 -0.0519249498
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.440789347 -0.529415204 -0.0519321219
+    endloop
+  endfacet
+  facet normal 2.73756404E-012 0.998858557 -0.0477659113
+    outer loop
+      vertex -0.510789435 -0.713575676 0.068007641
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.510789533 -0.718809786 -0.0414456238
+    endloop
+  endfacet
+  facet normal 3.94218571E-009 0.998858558 -0.0477659088
+    outer loop
+      vertex -0.440789347 -0.718809786 -0.0414456238
+      vertex -0.510789435 -0.713575676 0.068007641
+      vertex -0.440789347 -0.713575677 0.0680076416
+    endloop
+  endfacet
+  facet normal 0.00451085653 -0.999880043 0.0148172661
+    outer loop
+      vertex -0.440789347 -0.529415204 -0.0519321219
+      vertex -0.510789397 -0.528032778 0.0626652343
+      vertex -0.51078958 -0.529730896 -0.0519249498
+    endloop
+  endfacet
+  facet normal -0.00147137708 -0.9999366 0.0111637921
+    outer loop
+      vertex -0.510789397 -0.528032778 0.0626652343
+      vertex -0.440789347 -0.529415204 -0.0519321219
+      vertex -0.440789347 -0.528135696 0.0626728809
+    endloop
+  endfacet
+endsolid bodyB
+solid bodyA
+  facet normal 0.290553109 -3.06075535E-006 0.956858867
+    outer loop
+      vertex 0.0999005454 -0.369621646 0.16772795
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.0232856607 -0.369621589 0.190992293
+    endloop
+  endfacet
+  facet normal 0.290553241 -3.03091855E-006 0.956858827
+    outer loop
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.0999005454 -0.369621646 0.16772795
+      vertex 0.0999012214 -9.06014682E-007 0.167728916
+    endloop
+  endfacet
+  facet normal -0.279158286 1.12280344E-006 -0.9602451
+    outer loop
+      vertex 0.0485044684 -0.369621893 0.287966019
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.125241892 -0.369621955 0.265657248
+    endloop
+  endfacet
+  facet normal -0.279102736 -1.13846098E-005 -0.960261247
+    outer loop
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.0485044684 -0.369621893 0.287966019
+      vertex 0.048387693 1.05677683E-006 0.287995578
+    endloop
+  endfacet
+  facet normal -1.58162397E-006 -1 -2.7532619E-006
+    outer loop
+      vertex 0.125241892 -0.369621955 0.265657248
+      vertex 0.0232856607 -0.369621589 0.190992293
+      vertex 0.0999005454 -0.369621646 0.16772795
+    endloop
+  endfacet
+  facet normal -1.60317747E-006 -1 -2.72383025E-006
+    outer loop
+      vertex 0.0232856607 -0.369621589 0.190992293
+      vertex 0.125241892 -0.369621955 0.265657248
+      vertex 0.0485044684 -0.369621893 0.287966019
+    endloop
+  endfacet
+  facet normal 0.290553241 3.03090369E-006 0.956858827
+    outer loop
+      vertex 0.0999005454 0.369621646 0.16772795
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.0999012214 -9.06014682E-007 0.167728916
+    endloop
+  endfacet
+  facet normal 0.290553109 3.06074129E-006 0.956858867
+    outer loop
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.0999005454 0.369621646 0.16772795
+      vertex 0.0232856607 0.369621589 0.190992293
+    endloop
+  endfacet
+  facet normal -0.279102736 1.13846749E-005 -0.960261247
+    outer loop
+      vertex 0.0485044684 0.369621893 0.287966019
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.048387693 1.05677683E-006 0.287995578
+    endloop
+  endfacet
+  facet normal -0.279158286 -1.12281024E-006 -0.9602451
+    outer loop
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.0485044684 0.369621893 0.287966019
+      vertex 0.125241892 0.369621955 0.265657248
+    endloop
+  endfacet
+  facet normal -1.58162397E-006 1 -2.7532619E-006
+    outer loop
+      vertex 0.0232856607 0.369621589 0.190992293
+      vertex 0.125241892 0.369621955 0.265657248
+      vertex 0.0999005454 0.369621646 0.16772795
+    endloop
+  endfacet
+  facet normal -1.60317747E-006 1 -2.72383025E-006
+    outer loop
+      vertex 0.125241892 0.369621955 0.265657248
+      vertex 0.0232856607 0.369621589 0.190992293
+      vertex 0.0485044684 0.369621893 0.287966019
+    endloop
+  endfacet
+  facet normal -0.968111427 1.11605422E-006 0.25051999
+    outer loop
+      vertex 0.125241892 -0.369621955 0.265657248
+      vertex 0.0999012214 -9.06014682E-007 0.167728916
+      vertex 0.0999005454 -0.369621646 0.16772795
+    endloop
+  endfacet
+  facet normal -0.968111495 1.04182999E-006 0.250519728
+    outer loop
+      vertex 0.0999012214 -9.06014682E-007 0.167728916
+      vertex 0.125241892 -0.369621955 0.265657248
+      vertex 0.125242366 1.11914757E-006 0.265657543
+    endloop
+  endfacet
+  facet normal -0.968111495 -1.11604956E-006 0.250519728
+    outer loop
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.0999005454 0.369621646 0.16772795
+      vertex 0.0999012214 -9.06014682E-007 0.167728916
+    endloop
+  endfacet
+  facet normal -0.968111427 -1.04183601E-006 0.25051999
+    outer loop
+      vertex 0.0999005454 0.369621646 0.16772795
+      vertex 0.125242366 1.11914757E-006 0.265657543
+      vertex 0.125241892 0.369621955 0.265657248
+    endloop
+  endfacet
+  facet normal 0.9678088 -8.71112912E-007 -0.251686566
+    outer loop
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.0485044684 -0.369621893 0.287966019
+      vertex 0.0232856607 -0.369621589 0.190992293
+    endloop
+  endfacet
+  facet normal 0.968111426 0.000325890961 -0.25051978
+    outer loop
+      vertex 0.0485044684 -0.369621893 0.287966019
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+      vertex 0.048387693 1.05677683E-006 0.287995578
+    endloop
+  endfacet
+  facet normal 0.968111479 8.74758278E-007 -0.250519787
+    outer loop
+      vertex 0.0232856607 0.369621589 0.190992293
+      vertex 0.048387693 1.05677683E-006 0.287995578
+      vertex 0.023286254 -8.48891568E-007 0.190993295
+    endloop
+  endfacet
+  facet normal 0.967808748 -0.000325890506 -0.251686551
+    outer loop
+      vertex 0.048387693 1.05677683E-006 0.287995578
+      vertex 0.0232856607 0.369621589 0.190992293
+      vertex 0.0485044684 0.369621893 0.287966019
+    endloop
+  endfacet
+endsolid bodyA
diff --git a/models/nancolor.obj b/models/nancolor.obj
new file mode 100644
index 00000000..1f8a3a0f
--- /dev/null
+++ b/models/nancolor.obj
@@ -0,0 +1,21 @@
+v -0.500000 -0.500000 -0.500000 0.5365 0.9971 0.5398
+v -0.500000 -0.500000 0.500000 0.6006 0.9889 0.5296
+v -0.500000 0.500000 -0.500000 0.5224 0.9981 0.5371
+v -0.500000 0.500000 0.500000 0.5783 0.9922 0.5400
+v 0.500000 -0.500000 -0.500000 0.6500 0.9760 0.5295
+v 0.500000 -0.500000 0.500000 0.6701 0.9693 0.5287
+v 0.500000 0.500000 -0.500000 nan nan nan
+v 0.500000 0.500000 0.500000 0.6701 0.9693 0.5287
+f 2 4 1
+f 5 2 1
+f 1 4 3
+f 3 5 1
+f 2 8 4
+f 6 2 5
+f 6 8 2
+f 4 8 3
+f 7 5 3
+f 3 8 7
+f 7 6 5
+f 8 6 7
+
diff --git a/trimesh/exchange/obj.py b/trimesh/exchange/obj.py
index 59da6fe5..821bd960 100644
--- a/trimesh/exchange/obj.py
+++ b/trimesh/exchange/obj.py
@@ -180,9 +180,9 @@ def load_obj(file_obj,
                 # survive index errors as sometimes we
                 # want materials without UV coordinates
                 uv = vt[mask_vt]
-            except BaseException as E:
+            except BaseException:
+                log.warning('index failed on UV coordinates, skipping!')
                 uv = None
-                raise E
 
             # mask vertices and use new faces
             mesh.update({'vertices': v[mask_v].copy(),
diff --git a/trimesh/exchange/stl.py b/trimesh/exchange/stl.py
index 78432a62..46f42293 100644
--- a/trimesh/exchange/stl.py
+++ b/trimesh/exchange/stl.py
@@ -151,50 +151,57 @@ def load_stl_ascii(file_obj):
       face_normals: (m,3) float, normal vector of each face
     """
 
-    # the first line is the header
-    header = file_obj.readline()
-    # make sure header is a string, not bytes
-    if hasattr(header, 'decode'):
-        try:
-            header = util.decode_text(header)
-        except BaseException:
-            header = ''
-    # save header to metadata
-    metadata = {'header': header}
-
     # read all text into one string
-    text = file_obj.read()
-    # try to convert bytes to string
-    text = util.decode_text(text)
-    # split by endsolid keyword
-    text = text.lower().split('endsolid')[0]
-    # create array of splits
-    blob = np.array(text.strip().split())
-
-    # there are 21 'words' in each face
-    face_len = 21
-
-    # length of blob should be multiple of face_len
-    if (len(blob) % face_len) != 0:
-        raise HeaderError('Incorrect length STL file!')
-    face_count = int(len(blob) / face_len)
-
-    # this offset is to be added to a fixed set of tiled indices
-    offset = face_len * np.arange(face_count).reshape((-1, 1))
-    normal_index = np.tile([2, 3, 4], (face_count, 1)) + offset
-    vertex_index = np.tile([8, 9, 10,
-                            12, 13, 14,
-                            16, 17, 18], (face_count, 1)) + offset
-
-    # faces are groups of three sequential vertices
-    faces = np.arange(face_count * 3).reshape((-1, 3))
-    face_normals = blob[normal_index].astype('<f8')
-    vertices = blob[vertex_index.reshape((-1, 3))].astype('<f8')
-
-    return {'vertices': vertices,
-            'faces': faces,
-            'metadata': metadata,
-            'face_normals': face_normals}
+    raw = util.decode_text(file_obj.read()).strip().lower()
+
+    # split into solid body
+    solids = raw.split('endsolid')
+
+    kwargs = {}
+
+    for solid in solids:
+
+        stripped = solid.split('solid', 1)
+        if len(stripped) != 2:
+            continue
+        header, text = stripped[1].split('\n', 1)
+
+        name = header.strip()
+        if name in kwargs or len(name) == 0:
+            name = '{}_{}'.format(name, util.unique_id())
+
+        # create array of splits
+        blob = np.array(text.strip().split())
+
+        # there are 21 'words' in each face
+        face_len = 21
+
+        # length of blob should be multiple of face_len
+        if (len(blob) % face_len) != 0:
+            util.log.warning('skipping solid!')
+            continue
+        face_count = int(len(blob) / face_len)
+
+        # this offset is to be added to a fixed set of tiled indices
+        offset = face_len * np.arange(face_count).reshape((-1, 1))
+        normal_index = np.tile([2, 3, 4], (face_count, 1)) + offset
+        vertex_index = np.tile([8, 9, 10,
+                                12, 13, 14,
+                                16, 17, 18], (face_count, 1)) + offset
+
+        # faces are groups of three sequential vertices
+        faces = np.arange(face_count * 3).reshape((-1, 3))
+        face_normals = blob[normal_index].astype('<f8')
+        vertices = blob[vertex_index.reshape((-1, 3))].astype('<f8')
+
+        kwargs[name] = {'vertices': vertices,
+                        'faces': faces,
+                        'face_normals': face_normals}
+
+    if len(kwargs) == 1:
+        return next(iter(kwargs.values()))
+
+    return {'geometry': kwargs}
 
 
 def export_stl(mesh):
diff --git a/trimesh/inertia.py b/trimesh/inertia.py
index d4e8a8f9..dfdb2f94 100644
--- a/trimesh/inertia.py
+++ b/trimesh/inertia.py
@@ -12,10 +12,6 @@ import numpy as np
 
 from trimesh import util
 
-# a matrix where all non- diagonal terms are -1.0
-# and all diagonal terms are 1.0
-negate_nondiagonal = (np.eye(3, dtype=np.float64) * 2) - 1
-
 
 def cylinder_inertia(mass, radius, height, transform=None):
     """
@@ -95,7 +91,7 @@ def principal_axis(inertia):
     # np.linalg.svd, np.linalg.eig, np.linalg.eigh
     # moment of inertia is square symmetric matrix
     # eigh has the best numeric precision in tests
-    components, vectors = np.linalg.eigh(inertia * negate_nondiagonal)
+    components, vectors = np.linalg.eigh(inertia)
 
     # eigh returns them as column vectors, change them to row vectors
     vectors = vectors.T
@@ -136,9 +132,8 @@ def transform_inertia(transform, inertia_tensor):
         raise ValueError('inertia_tensor must be (3, 3)!')
 
     transformed = util.multi_dot([rotation,
-                                  inertia_tensor * negate_nondiagonal,
+                                  inertia_tensor,
                                   rotation.T])
-    transformed *= negate_nondiagonal
     return transformed
 
 
diff --git a/trimesh/triangles.py b/trimesh/triangles.py
index b6d621df..6b1fdf28 100644
--- a/trimesh/triangles.py
+++ b/trimesh/triangles.py
@@ -266,11 +266,11 @@ def mass_properties(triangles,
         (volume * (center_mass[[0, 2]]**2).sum())
     inertia[2, 2] = integrated[4] + integrated[5] - \
         (volume * (center_mass[[0, 1]]**2).sum())
-    inertia[0, 1] = (
+    inertia[0, 1] = - (
         integrated[7] - (volume * np.product(center_mass[[0, 1]])))
-    inertia[1, 2] = (
+    inertia[1, 2] = - (
         integrated[8] - (volume * np.product(center_mass[[1, 2]])))
-    inertia[0, 2] = (
+    inertia[0, 2] = - (
         integrated[9] - (volume * np.product(center_mass[[0, 2]])))
     inertia[2, 0] = inertia[0, 2]
     inertia[2, 1] = inertia[1, 2]
diff --git a/trimesh/version.py b/trimesh/version.py
index cc301490..ada695aa 100644
--- a/trimesh/version.py
+++ b/trimesh/version.py
@@ -1,1 +1,1 @@
-__version__ = '3.8.15'
+__version__ = '3.8.17'
diff --git a/trimesh/visual/color.py b/trimesh/visual/color.py
index 961f6044..d3d75b93 100644
--- a/trimesh/visual/color.py
+++ b/trimesh/visual/color.py
@@ -69,11 +69,14 @@ class ColorVisuals(Visuals):
                                           dtype=np.uint8),
             'material_shine': 77.0}
 
-        if face_colors is not None:
-            self.face_colors = face_colors
+        try:
+            if face_colors is not None:
+                self.face_colors = face_colors
 
-        if vertex_colors is not None:
-            self.vertex_colors = vertex_colors
+            if vertex_colors is not None:
+                self.vertex_colors = vertex_colors
+        except ValueError:
+            util.log.warning('unable to convert colors!')
 
     @caching.cache_decorator
     def transparency(self):
@@ -570,12 +573,16 @@ def to_rgba(colors, dtype=np.uint8):
     # integer value for opaque alpha given our datatype
     opaque = np.iinfo(dtype).max
 
+    if colors.dtype.kind == 'f':
+        # replace any `nan` or `inf` values with zero
+        colors[~np.isfinite(colors)] = 0.0
+
     if (colors.dtype.kind == 'f' and colors.max() < (1.0 + 1e-8)):
         colors = (colors * opaque).round().astype(dtype)
     elif (colors.max() <= opaque):
         colors = colors.astype(dtype)
     else:
-        raise ValueError('colors non- convertible!')
+        raise ValueError('colors non-convertible!')
 
     if util.is_shape(colors, (-1, 3)):
         # add an opaque alpha for RGB colors


Test patch:
diff --git a/tests/data/mass_properties.json b/tests/data/mass_properties.json
index 4c6a1f1b..97b16697 100644
--- a/tests/data/mass_properties.json
+++ b/tests/data/mass_properties.json
@@ -1,58 +1,58 @@
 [
-    {
-        "area": 9.388587, 
-        "center_mass": [
-            0, 
-            0.436105, 
-            -0.599006
-        ], 
-        "density": 0.036127, 
-        "filename": "angle_block.STL", 
-        "inertia": [
-            [
-                0.008312, 
-                0.0, 
-                0.0
-            ], 
-            [
-                0.0, 
-                0.012823, 
-                0.000807
-            ], 
-            [
-                0.0, 
-                0.000807, 
-                0.009005
-            ]
-        ], 
-        "mass": 0.041358
-    }, 
-    {
-        "area": 18.13860637, 
-        "center_mass": [
-            1.25, 
-            1.21830489, 
-            0.20084261
-        ], 
-        "density": 0.09754366, 
-        "filename": "idler_riser.STL", 
-        "inertia": [
-            [
-                0.12108243, 
-                0.0, 
-                0.0
-            ], 
-            [
-                0.0, 
-                0.09131981, 
-                0.00013893
-            ], 
-            [
-                0.0, 
-                0.00013893, 
-                0.2057651
-            ]
-        ], 
-        "mass": 0.14522859
-    }
+  {
+    "area": 9.388587,
+    "center_mass": [
+      0,
+      0.436105,
+      -0.599006
+    ],
+    "density": 0.036127,
+    "filename": "angle_block.STL",
+    "inertia": [
+      [
+        0.008312,
+        -0.0,
+        -0.0
+      ],
+      [
+        -0.0,
+        0.012823,
+        -0.000807
+      ],
+      [
+        -0.0,
+        -0.000807,
+        0.009005
+      ]
+    ],
+    "mass": 0.041358
+  },
+  {
+    "area": 18.13860637,
+    "center_mass": [
+      1.25,
+      1.21830489,
+      0.20084261
+    ],
+    "density": 0.09754366,
+    "filename": "idler_riser.STL",
+    "inertia": [
+      [
+        0.12108243,
+        -0.0,
+        -0.0
+      ],
+      [
+        -0.0,
+        0.09131981,
+        -0.00013893
+      ],
+      [
+        -0.0,
+        -0.00013893,
+        0.2057651
+      ]
+    ],
+    "mass": 0.14522859
+  }
 ]
\ No newline at end of file
diff --git a/tests/test_inertia.py b/tests/test_inertia.py
index 269e8703..ebdc8970 100644
--- a/tests/test_inertia.py
+++ b/tests/test_inertia.py
@@ -4,6 +4,21 @@ except BaseException:
     import generic as g
 
 
+# tetToTris extracts outward facing triangles
+# of a tetrahedron given in this order
+#    1
+#  / | \
+# 0-----2
+#  \ | /
+#    3
+
+triIdxs = g.np.ravel([[2, 1, 0], [3, 0, 1], [3, 2, 0], [3, 1, 2]])
+
+
+def tetToTris(tet):
+    return g.np.reshape(tet[triIdxs], (-1, 3))
+
+
 class InertiaTest(g.unittest.TestCase):
 
     def test_inertia(self):
@@ -99,6 +114,161 @@ class InertiaTest(g.unittest.TestCase):
                 elif hasattr(p.primitive, 'center'):
                     p.primitive.center = g.np.random.random(3)
 
+    def test_tetrahedron(self):
+        # Based on the 'numerical example' of the paper:
+
+        # Explicit Exact Formulas for the 3-D Tetrahedron Inertia Tensor
+        # in Terms of its Vertex Coordinates, [F. Tonon, 2004]
+        # http://thescipub.com/pdf/jmssp.2005.8.11.pdf
+
+        # set up given vertices
+        vertices = g.np.float32([[8.3322, -11.86875, 0.93355],
+                                 [0.75523, 5., 16.37072],
+                                 [52.61236, 5., -5.3858],
+                                 [2., 5., 3.]])
+
+        # set up a simple trimesh tetrahedron
+        tris = tetToTris(g.np.int32([0, 1, 2, 3]))
+        tm_tet = g.trimesh.Trimesh(vertices, tris)
+
+        # 'ground truth' values from the paper
+        # however, there are some minor flaws
+
+        # a small deviation in the last decimal of the mass-centers' x:
+        CM_gt = g.np.float32([15.92492, 0.78281, 3.72962])
+
+        # moment of inertia values
+        a_mi, b_mi, c_mi = [43520.33257, 194711.28938, 191168.76173]
+        # principle inertia values
+        a_pi, b_pi, c_pi = [4417.6615, -46343.16662, 11996.20119]
+
+        # NOTE: I guess there is a mistake in the paper
+        # b' (Eq. 9e) computes from x and z values
+        # c' (Eq. 9f) computes from x and y values
+        # therefore the given matrix E_Q (Eq. 2) is not correct
+        # b' and c' need to be swapped like this:
+        MI_gt = g.np.float32([[a_mi, -c_pi, -b_pi],
+                              [-c_pi, b_mi, -a_pi],
+                              [-b_pi, -a_pi, c_mi]])
+
+        # check center of mass
+        assert g.np.allclose(CM_gt, tm_tet.center_mass)
+
+        # check moment of inertia tensor
+        assert g.np.allclose(MI_gt, tm_tet.moment_inertia)
+
+    def test_cube_with_tetras(self):
+
+        # set up a unit cube, vertices in this order:
+        #   1-----2
+        #  /|    /|
+        # 0-+---3 |
+        # | 5---+-6
+        # |/    |/
+        # 4-----7
+
+        vertices = g.np.float32([[-1, -1, 1],
+                                 [-1, 1, 1],
+                                 [1, 1, 1],
+                                 [1, -1, 1],
+                                 [-1, -1, -1],
+                                 [-1, 1, -1],
+                                 [1, 1, -1],
+                                 [1, -1, -1]]) * 0.5
+
+        # 6 quad faces for the cube
+        quads = g.np.int32([[3, 2, 1, 0],
+                            [0, 1, 5, 4],
+                            [1, 2, 6, 5],
+                            [2, 3, 7, 6],
+                            [3, 0, 4, 7],
+                            [4, 5, 6, 7]])
+
+        # set up two different tetraherdalizations of a cube
+        # using 5 tets (1 in the middle and 4 around)
+        tetsA = g.np.int32([[0, 1, 3, 4],
+                            [1, 2, 3, 6],
+                            [1, 4, 5, 6],
+                            [3, 4, 6, 7],
+                            [1, 3, 4, 6]])
+        # using 6 sets (around the cube diagonal)
+        tetsB = g.np.int32([[0, 1, 2, 6],
+                            [0, 1, 6, 5],
+                            [0, 4, 5, 6],
+                            [0, 4, 6, 7],
+                            [0, 3, 7, 6],
+                            [0, 2, 3, 6]])
+
+        # create a trimesh cube from vertices and faces
+        tm_cube = g.trimesh.Trimesh(vertices, quads)
+
+        # https://en.wikipedia.org/wiki/List_of_moments_of_inertia
+        # ground truth for a unit cube with side length s = mass m = 1
+        # I = 1/6 * m * s^2
+        MI_gt = g.np.eye(3) / 6
+        assert g.np.allclose(MI_gt, tm_cube.moment_inertia)
+
+        # compare both tetrahedralizations
+        # should be equivalent to each other and to the cube
+        for tets in [tetsA, tetsB]:
+            # create trimesh tets from vertices and triangles
+            tm_tets = [g.trimesh.Trimesh(vertices, tetToTris(t)) for t in tets]
+
+            # get mass, center of mass, and moment of inertia for each tet
+            Ms = [tm_tet.mass for tm_tet in tm_tets]
+            CMs = [tm_tet.center_mass for tm_tet in tm_tets]
+            MIs = [tm_tet.moment_inertia for tm_tet in tm_tets]
+
+            # compute total mass and center of mass
+            mass = g.np.sum(Ms)
+            center_mass = g.np.dot(Ms, CMs)
+
+            # compare with unit cube
+            assert g.np.abs(tm_cube.mass - mass) < 1e-6
+            assert g.np.allclose(tm_cube.center_mass, center_mass)
+
+            # the moment of inertia tensors for the individual tetrahedra
+            # have to be re-assembled, using the parallel-axis-theorem
+            # https://en.wikipedia.org/wiki/Parallel_axis_theorem#Tensor_generalization
+
+            E = g.np.eye(3)
+            MI = g.np.zeros((3, 3), g.np.float32)
+            for i in range(len(tm_tets)):
+                R = CMs[i] - center_mass
+                MI += MIs[i] + Ms[i] * (g.np.dot(R, R) * E - g.np.outer(R, R))
+
+            assert g.np.allclose(MI_gt, MI)
+
+
+class MassTests(g.unittest.TestCase):
+
+    def setUp(self):
+        # inertia numbers pulled from solidworks
+        self.truth = g.data['mass_properties']
+        self.meshes = dict()
+        for data in self.truth:
+            filename = data['filename']
+            self.meshes[filename] = g.get_mesh(filename)
+
+    def test_mass(self):
+
+        for truth in self.truth:
+            mesh = self.meshes[truth['filename']]
+            calc = g.trimesh.triangles.mass_properties(
+                triangles=mesh.triangles,
+                density=truth['density'],
+                skip_inertia=False)
+
+            for key, value in calc.items():
+                if key not in truth:
+                    continue
+                if not g.np.allclose(calc[key], truth[key], atol=1e-2):
+                    raise ValueError('{}({}):\n{}\n!=\n{}'.format(
+                        truth['filename'],
+                        key,
+                        calc[key],
+                        g.np.array(truth[key])))
+
 
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()
diff --git a/tests/test_obj.py b/tests/test_obj.py
index 2875b7ac..ce4533b7 100644
--- a/tests/test_obj.py
+++ b/tests/test_obj.py
@@ -200,6 +200,14 @@ class OBJTest(g.unittest.TestCase):
 
         assert g.np.isclose(e.area, s.area, rtol=.01)
 
+    def test_edge_cases(self):
+        # a mesh with some NaN colors
+        n = g.get_mesh('nancolor.obj')
+        assert n.faces.shape == (12, 3)
+
+        v = g.get_mesh('cubevt.obj')
+        assert v.faces.shape == (12, 3)
+
 
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()
diff --git a/tests/test_stl.py b/tests/test_stl.py
index f3f7c865..aa2af99e 100644
--- a/tests/test_stl.py
+++ b/tests/test_stl.py
@@ -67,6 +67,10 @@ class STLTests(g.unittest.TestCase):
         assert g.np.allclose(m.vertex_attributes['random'], random[v_mask])
         assert len(m.vertices) == len(v_mask)
 
+    def test_ascii_multibody(self):
+        s = g.get_mesh('multibody.stl')
+        assert len(s.geometry) == 2
+
 
 if __name__ == '__main__':
     g.trimesh.util.attach_to_log()
diff --git a/tests/test_util.py b/tests/test_util.py
index 05241301..77be5746 100644
--- a/tests/test_util.py
+++ b/tests/test_util.py
@@ -189,42 +189,6 @@ class ContainsTest(unittest.TestCase):
             raise ValueError('contains test does not match truth!')
 
 
-class MassTests(unittest.TestCase):
-
-    def setUp(self):
-        # inertia numbers pulled from solidworks
-        self.truth = g.data['mass_properties']
-        self.meshes = dict()
-        for data in self.truth:
-            filename = data['filename']
-            self.meshes[filename] = g.get_mesh(filename)
-
-    def test_mass(self):
-        def check_parameter(a, b):
-            diff = np.abs(np.array(a) - np.array(b))
-            check = (diff < TOL_CHECK).all()
-            return check
-
-        for truth in self.truth:
-            mesh = self.meshes[truth['filename']]
-            calculated = trimesh.triangles.mass_properties(
-                triangles=mesh.triangles, density=truth['density'], skip_inertia=False)
-
-            parameter_count = 0
-            for parameter in calculated.keys():
-                if not (parameter in truth):
-                    continue
-                parameter_ok = check_parameter(calculated[parameter],
-                                               truth[parameter])
-                if not parameter_ok:
-                    log.error('Parameter %s failed on file %s!',
-                              parameter, truth['filename'])
-                self.assertTrue(parameter_ok)
-                parameter_count += 1
-            log.info('%i mass parameters confirmed for %s',
-                     parameter_count, truth['filename'])
-
-
 class IOWrapTests(unittest.TestCase):
 
     def test_io_wrap(self):