diff --git a/src/stlProcess.jl b/src/stlProcess.jl index 250ed1b..ad736ea 100644 --- a/src/stlProcess.jl +++ b/src/stlProcess.jl @@ -16,15 +16,52 @@ struct Properties end end -function get_mass_properties(triangles) +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: https://github.com/WoLpH/numpy-stl/blob/42b6c67324e13b8712af6730f77e5e9544ef63b0/stl/base.py#L362 https://www.geometrictools.com/Documentation/PolyhedralMassProperties.pdf """ - x = reduce(hcat, [[v[1] for v in tri] for tri in triangles])' - y = reduce(hcat, [[v[2] for v in tri] for tri in triangles])' - z = reduce(hcat, [[v[3] for v in tri] for tri in triangles])' + 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])' function subexpression(x) w0, w1, w2 = x[:, 1], x[:, 2], x[:, 3] @@ -76,6 +113,6 @@ function get_mass_properties(triangles) return Properties(volume, center_of_gravity, inertia) end -export get_mass_properties +export get_mass_properties, get_volume end # module diff --git a/test/runtests.jl b/test/runtests.jl index 6b69f92..e97b8c8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,12 +8,25 @@ using MeshIO using LinearAlgebra @testset "simple cube stl" begin - # https://hepweb.ucsd.edu/ph110b/110b_notes/node26.html - stl = load(raw"test\test_assets\cubeblender.stl") + """ + inertia math: + https://hepweb.ucsd.edu/ph110b/110b_notes/node26.html - props = get_mass_properties(stl) + Cube is a standard cube with a length of 2 for each side. + """ + cube_path = raw"test_assets\cubeblender.stl" + stl = load(cube_path) - @test props.volume == 8.0 - @test all(props.center_of_gravity .== 0.0) - @test all(eigvals(props.inertia) .≈ (5.0 + 1 / 3)) + @testset "get_mass_properties" begin + props = get_mass_properties(stl) + + @test props.volume == 8.0 + @test all(props.center_of_gravity .== 0.0) + @test all(eigvals(props.inertia) .≈ (5.0 + 1 / 3)) + + @testset "get_volume" begin + @test get_volume(stl; scale=1) ≈ props.volume rtol = 0.01 + @test get_volume(stl; scale=4) ≈ (4 * 2)^3 rtol = 0.01 + end + end end