From 2fd254e08d0645d85fddcd85d3f2d04216f8c955 Mon Sep 17 00:00:00 2001
From: Anson <anson@ansonbiggs.com>
Date: Thu, 31 Mar 2022 10:12:47 -0700
Subject: [PATCH] new volume function based on mass props

---
 src/stlProcess.jl | 67 ++++++++++++++++++++---------------------------
 test/runtests.jl  | 44 ++++++++++++++++++++++++++++---
 2 files changed, 70 insertions(+), 41 deletions(-)

diff --git a/src/stlProcess.jl b/src/stlProcess.jl
index ad736ea..0a3ef0f 100644
--- a/src/stlProcess.jl
+++ b/src/stlProcess.jl
@@ -16,43 +16,6 @@ struct Properties
     end
 end
 
-function get_volume(triangles; scale=1)
-    """
-    Reference:
-        https://people.eecs.berkeley.edu/~wkahan/VtetLang.pdf
-    """
-
-    volume = 0.0
-    for tri in triangles
-        a, b, c = tri .* scale
-
-        W = sqrt(sum((a .- b) .^ 2))
-        V = sqrt(sum((b .- c) .^ 2))
-        U = sqrt(sum((c .- a) .^ 2))
-
-        v = sqrt(sum(a .^ 2))
-        u = sqrt(sum(b .^ 2))
-        w = sqrt(sum(c .^ 2))
-
-        X = (w - U + v) * (U + v + w)
-        Y = (u - V + w) * (V + w + u)
-        Z = (v - W + u) * (W + u + v)
-
-        x = (U - v + w) * (v - w + U)
-        y = (V - w + u) * (w - u + V)
-        z = (W - u + v) * (u - v + W)
-
-        ξ = sqrt(x * Y * Z)
-        η = sqrt(y * Z * X)
-        ζ = sqrt(z * X * Y)
-        λ = sqrt(x * y * z)
-
-        volume += sqrt((ξ + η + ζ - λ) * (λ + ξ + η - ζ) * (η + ζ + λ - ξ) * (ζ + λ + ξ - η)) / (192 * u * v * w)
-    end
-
-    return volume
-end
-
 function get_mass_properties(triangles; scale=1)
     """
     Reference:
@@ -113,6 +76,34 @@ function get_mass_properties(triangles; scale=1)
     return Properties(volume, center_of_gravity, inertia)
 end
 
-export get_mass_properties, get_volume
+function fast_volume(triangles; scale=1)
+    """
+    Faster algorithm to get volume of mesh. Is just `get_mass_properties` without any calculations not relevant to volume.
+
+    Reference:
+    	https://github.com/WoLpH/numpy-stl/blob/42b6c67324e13b8712af6730f77e5e9544ef63b0/stl/base.py#L362
+    	https://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf
+    """
+    x = reduce(hcat, [[v[1] .* scale for v in tri] for tri in triangles])'
+    y = reduce(hcat, [[v[2] .* scale for v in tri] for tri in triangles])'
+    z = reduce(hcat, [[v[3] .* scale for v in tri] for tri in triangles])'
+
+    w0, w1, w2 = x[:, 1], x[:, 2], x[:, 3]
+    temp0 = w0 + w1
+    f1 = temp0 + w2
+
+    x0, x1, x2 = x[:, 1], x[:, 2], x[:, 3]
+    y0, y1, y2 = y[:, 1], y[:, 2], y[:, 3]
+    z0, z1, z2 = z[:, 1], z[:, 2], z[:, 3]
+    a1, b1, c1 = x1 - x0, y1 - y0, z1 - z0
+    a2, b2, c2 = x2 - x0, y2 - y0, z2 - z0
+    d0, d1, d2 = b1 .* c2 .- b2 .* c1, a2 .* c1 .- a1 .* c2, a1 .* b2 .- a2 .* b1
+
+    volume = sum(d0 .* f1) / 6
+
+    return volume
+end
+
+export get_mass_properties, fast_volume
 
 end # module
diff --git a/test/runtests.jl b/test/runtests.jl
index ff3c121..dedf753 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,44 @@ using MeshIO
 
 using LinearAlgebra
 
+function _check_volume(triangles; scale=1)
+    """
+    Slow algorithm just used to test the other algorithms
+    Reference:
+        https://people.eecs.berkeley.edu/~wkahan/VtetLang.pdf
+    """
+
+    volume = 0.0
+    for tri in triangles
+        a, b, c = tri .* scale
+
+        W = sqrt(sum((a .- b) .^ 2))
+        V = sqrt(sum((b .- c) .^ 2))
+        U = sqrt(sum((c .- a) .^ 2))
+
+        v = sqrt(sum(a .^ 2))
+        u = sqrt(sum(b .^ 2))
+        w = sqrt(sum(c .^ 2))
+
+        X = (w - U + v) * (U + v + w)
+        Y = (u - V + w) * (V + w + u)
+        Z = (v - W + u) * (W + u + v)
+
+        x = (U - v + w) * (v - w + U)
+        y = (V - w + u) * (w - u + V)
+        z = (W - u + v) * (u - v + W)
+
+        ξ = sqrt(x * Y * Z)
+        η = sqrt(y * Z * X)
+        ζ = sqrt(z * X * Y)
+        λ = sqrt(x * y * z)
+
+        volume += sqrt((ξ + η + ζ - λ) * (λ + ξ + η - ζ) * (η + ζ + λ - ξ) * (ζ + λ + ξ - η)) / (192 * u * v * w)
+    end
+
+    return volume
+end
+
 @testset "simple cube stl" begin
     """
     inertia math:
@@ -25,14 +63,14 @@ using LinearAlgebra
         @test all(eigvals(props.inertia) .≈ (5.0 + 1 / 3))
 
         @testset "get_volume function" begin
-            @test get_volume(stl; scale=1) ≈ props.volume rtol = 0.01
-            @test get_volume(stl; scale=4) ≈ (4 * 2)^3 rtol = 0.01
+            @test _check_volume(stl; scale=1) ≈ props.volume rtol = 0.01
+            @test _check_volume(stl; scale=4) ≈ (4 * 2)^3 rtol = 0.01
         end
     end
     @testset "Compare volumes with scaling" begin
         for scale in 1:5
             props = get_mass_properties(stl; scale=scale)
-            volume = get_volume(stl; scale=scale)
+            volume = _check_volume(stl; scale=scale)
 
             @test props.volume ≈ volume rtol = 0.01
         end