mirror of
https://gitlab.com/MisterBiggs/stl-process.git
synced 2025-08-03 12:01:33 +00:00
203 lines
7.4 KiB
Julia
203 lines
7.4 KiB
Julia
module stlProcess
|
||
|
||
using MeshIO
|
||
using FileIO
|
||
using Roots
|
||
using LinearAlgebra
|
||
|
||
struct Properties
|
||
volume::Float64
|
||
center_of_gravity::Vector{Float64}
|
||
inertia::Matrix{Float64}
|
||
surface_area::Float64
|
||
characteristic_length::Float64
|
||
sb_values::Vector{Float64}
|
||
|
||
function Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
|
||
@assert size(center_of_gravity) == (3,)
|
||
@assert size(sb_values) == (3,)
|
||
@assert size(inertia) == (3, 3)
|
||
|
||
return new(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
|
||
end
|
||
end
|
||
|
||
function get_mass_properties(triangles; scale=1.0)
|
||
"""
|
||
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])'
|
||
points = collect(Set(reduce(vcat, [[Array(Float64.(v)) .* scale for v in tri] for tri in triangles])))
|
||
|
||
function subexpression(x)
|
||
w0, w1, w2 = x[:, 1], x[:, 2], x[:, 3]
|
||
temp0 = w0 + w1
|
||
f1 = temp0 + w2
|
||
temp1 = w0 .* w0
|
||
temp2 = temp1 .+ w1 .* temp0
|
||
f2 = temp2 .+ w2 .* f1
|
||
f3 = w0 .* temp1 .+ w1 .* temp2 .+ w2 .* f2
|
||
g0 = f2 .+ w0 .* (f1 .+ w0)
|
||
g1 = f2 .+ w1 .* (f1 .+ w1)
|
||
g2 = f2 .+ w2 .* (f1 .+ w2)
|
||
return f1, f2, f3, g0, g1, g2
|
||
end
|
||
|
||
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
|
||
|
||
f1x, f2x, f3x, g0x, g1x, g2x = subexpression(x)
|
||
f1y, f2y, f3y, g0y, g1y, g2y = subexpression(y)
|
||
f1z, f2z, f3z, g0z, g1z, g2z = subexpression(z)
|
||
|
||
intg = zeros(10)
|
||
intg[1] = sum(d0 .* f1x)
|
||
intg[2:4] = [sum(d0 .* f2x), sum(d1 .* f2y), sum(d2 .* f2z)]
|
||
intg[5:7] = [sum(d0 .* f3x), sum(d1 .* f3y), sum(d2 .* f3z)]
|
||
intg[8] = sum(d0 .* (y0 .* g0x + y1 .* g1x + y2 .* g2x))
|
||
intg[9] = sum(d1 .* (z0 .* g0y + z1 .* g1y + z2 .* g2y))
|
||
intg[10] = sum(d2 .* (x0 .* g0z + x1 .* g1z + x2 .* g2z))
|
||
intg = intg ./ [6, 24, 24, 24, 60, 60, 60, 120, 120, 120]
|
||
|
||
volume = intg[1]
|
||
|
||
center_of_gravity = intg[2:4] / volume
|
||
cogsq = center_of_gravity .^ 2
|
||
|
||
inertia = zeros((3, 3))
|
||
inertia[1, 1] = intg[6] + intg[7] - volume .* (cogsq[2] + cogsq[3])
|
||
inertia[2, 2] = intg[5] + intg[7] - volume .* (cogsq[3] + cogsq[1])
|
||
inertia[3, 3] = intg[5] + intg[6] - volume .* (cogsq[1] + cogsq[2])
|
||
inertia[1, 2] = inertia[2, 1] = -(intg[8] - volume .* center_of_gravity[1] .* center_of_gravity[2])
|
||
inertia[2, 3] = inertia[3, 2] = -(intg[9] - volume .* center_of_gravity[2] .* center_of_gravity[3])
|
||
inertia[1, 3] = inertia[3, 1] = -(intg[10] - volume .* center_of_gravity[3] .* center_of_gravity[1])
|
||
inertia = inertia / volume
|
||
|
||
# https://math.stackexchange.com/questions/128991/how-to-calculate-the-area-of-a-3d-triangle
|
||
surface_area = sum(norm.(eachrow([x0 y0 z0] - [x1 y1 z1]) .× eachrow([x1 y1 z1] - [x2 y2 z2])) / 2)
|
||
|
||
characteristic_length = calc_characteristic_length(points, inertia, center_of_gravity)
|
||
|
||
sb_values = solid_body(points)
|
||
|
||
return Properties(volume, center_of_gravity, inertia, surface_area, characteristic_length, sb_values)
|
||
end
|
||
|
||
function solid_body(points)
|
||
pts = reduce(hcat, points)
|
||
|
||
Xsb = abs(reduce(-, extrema(@view pts[1, :])))
|
||
Ysb = abs(reduce(-, extrema(@view pts[2, :])))
|
||
Zsb = abs(reduce(-, extrema(@view pts[3, :])))
|
||
|
||
return [Xsb, Ysb, Zsb]
|
||
end
|
||
|
||
function calc_characteristic_length(points, inertia, center_of_gravity)
|
||
"""
|
||
Calculate the Characteristic Length using eigenvectors.
|
||
|
||
More information on characteristic length: https://ntrs.nasa.gov/api/citations/20090019123/downloads/20090019123.pdf
|
||
More information on the geometry used: https://math.stackexchange.com/q/100447
|
||
"""
|
||
|
||
characteristic_points = []
|
||
# find the characteristic points for each eigenvector
|
||
for eig in eachrow(eigvecs(inertia))
|
||
|
||
# Find 3 points for each direction of the eigenvector
|
||
θs = map(
|
||
point -> acos(
|
||
clamp(dot(eig, point .- center_of_gravity) / (norm(point .- center_of_gravity) * norm(eig)), -1.0, 1.0),
|
||
),
|
||
points,
|
||
)
|
||
sort_index = sortperm(θs)
|
||
min_points = points[sort_index[1:3]]
|
||
max_points = points[sort_index[(end - 2):end]]
|
||
min_point = 0
|
||
max_point = 0
|
||
|
||
# Create a parameterized function using the eigenvector and center_of_gravity
|
||
line_func(t) = center_of_gravity .+ eig .* t
|
||
|
||
try
|
||
# create a plane using 3 points, then find where the parameterized function intercepts it
|
||
max_normal = (max_points[1] - max_points[2]) × (max_points[2] - max_points[3])
|
||
max_t = find_zero(t -> let
|
||
l = line_func(t)
|
||
sum(max_normal .* (l .- max_points[1]))
|
||
end, norm(eig))
|
||
|
||
max_point = line_func(max_t)
|
||
|
||
# Build a plane in the opposite direction, and find the intercept again
|
||
min_normal = (min_points[1] - min_points[2]) × (min_points[2] - min_points[3])
|
||
min_t = find_zero(t -> let
|
||
l = line_func(t)
|
||
sum(min_normal .* (l .- min_points[1]))
|
||
end, -norm(eig))
|
||
|
||
min_point = line_func(min_t)
|
||
|
||
catch e
|
||
if isa(e, Roots.ConvergenceFailed)
|
||
@warn "Characteristic Length Algorithm failed to converge, this usually means stl is flat. Setting length in dir to 0."
|
||
else
|
||
@warn "Unknown error when calculating Characteristic Length."
|
||
end
|
||
min_point = 0
|
||
max_point = 0
|
||
end
|
||
|
||
push!(characteristic_points, (min_point, max_point))
|
||
end
|
||
|
||
# Calculate the characteristic length.
|
||
return sum([sqrt(sum((max - min) .^ 2)) for (min, max) in characteristic_points]) / 3.0
|
||
end
|
||
|
||
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
|
||
|
||
function find_scale(trianges; desired_volume=1)
|
||
return find_zero(scale -> fast_volume(trianges; scale=scale) - desired_volume, 2.0)
|
||
end
|
||
|
||
export get_mass_properties, fast_volume, find_scale
|
||
|
||
end # module
|