You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
96 lines
3.4 KiB
96 lines
3.4 KiB
#include <cmath>
|
|
#include <cstdio>
|
|
#include <solve.h> // polymesh_t
|
|
|
|
// Vector subtraction: a - b
|
|
inline vector3d vec3_sub(const vector3d& a, const vector3d& b) { return {a.x - b.x, a.y - b.y, a.z - b.z}; }
|
|
|
|
// Vector cross product: a × b
|
|
inline vector3d vec3_cross(const vector3d& a, const vector3d& b)
|
|
{
|
|
return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x};
|
|
}
|
|
|
|
// Vector dot product: a · b
|
|
inline double vec3_dot(const vector3d& a, const vector3d& b) { return a.x * b.x + a.y * b.y + a.z * b.z; }
|
|
|
|
// Vector length: ||a||
|
|
inline double vec3_norm(const vector3d& a) { return std::sqrt(a.x * a.x + a.y * a.y + a.z * a.z); }
|
|
|
|
// Compute triangle area: 0.5 * || (b-a) × (c-a) ||
|
|
inline double triangle_area(const vector3d& a, const vector3d& b, const vector3d& c)
|
|
{
|
|
vector3d ab = vec3_sub(b, a);
|
|
vector3d ac = vec3_sub(c, a);
|
|
vector3d cross_product = vec3_cross(ab, ac);
|
|
return 0.5 * vec3_norm(cross_product);
|
|
}
|
|
|
|
// Compute triangle signed volume contribution (w.r.t origin): (a · (b × c)) / 6
|
|
inline double triangle_signed_volume(const vector3d& a, const vector3d& b, const vector3d& c)
|
|
{
|
|
vector3d cross_product = vec3_cross(b, c);
|
|
return vec3_dot(a, cross_product) / 6.0;
|
|
}
|
|
|
|
// Compute polygon face area (triangulate as fan)
|
|
double polygon_area(const vector3d* vertices, const uint32_t* face_indices, uint32_t n)
|
|
{
|
|
if (n < 3) return 0.0;
|
|
double area = 0.0;
|
|
const vector3d& v0 = vertices[face_indices[0]];
|
|
for (uint32_t i = 1; i < n - 1; ++i) {
|
|
const vector3d& v1 = vertices[face_indices[i]];
|
|
const vector3d& v2 = vertices[face_indices[i + 1]];
|
|
area += triangle_area(v0, v1, v2);
|
|
}
|
|
return area;
|
|
}
|
|
|
|
// Compute polygon signed volume contribution (sum of triangle contributions)
|
|
double polygon_signed_volume(const vector3d* vertices, const uint32_t* face_indices, uint32_t n)
|
|
{
|
|
if (n < 3) return 0.0;
|
|
double volume = 0.0;
|
|
const vector3d& v0 = vertices[face_indices[0]];
|
|
for (uint32_t i = 1; i < n - 1; ++i) {
|
|
const vector3d& v1 = vertices[face_indices[i]];
|
|
const vector3d& v2 = vertices[face_indices[i + 1]];
|
|
volume += triangle_signed_volume(v0, v1, v2);
|
|
}
|
|
return volume;
|
|
}
|
|
|
|
// 🟩 Main: compute total mesh surface area
|
|
double compute_surface_area(const polymesh_t* mesh)
|
|
{
|
|
if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { return 0.0; }
|
|
|
|
double total_area = 0.0;
|
|
const uint32_t* face_ptr = mesh->faces; // current index into faces array
|
|
|
|
for (uint32_t i = 0; i < mesh->num_faces; ++i) {
|
|
uint32_t n = mesh->vertex_counts[i];
|
|
total_area += polygon_area(mesh->vertices, face_ptr, n);
|
|
face_ptr += n; // move to the next face's start index
|
|
}
|
|
|
|
return total_area;
|
|
}
|
|
|
|
// 🟦 Main: compute total mesh volume (requires closed mesh with outward normals)
|
|
double compute_volume(const polymesh_t* mesh)
|
|
{
|
|
if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) { return 0.0; }
|
|
|
|
double total_volume = 0.0;
|
|
const uint32_t* face_ptr = mesh->faces;
|
|
|
|
for (uint32_t i = 0; i < mesh->num_faces; ++i) {
|
|
uint32_t n = mesh->vertex_counts[i];
|
|
total_volume += polygon_signed_volume(mesh->vertices, face_ptr, n);
|
|
face_ptr += n;
|
|
}
|
|
|
|
return std::abs(total_volume); // ensure positive volume
|
|
}
|