Skip to content

Instantly share code, notes, and snippets.

@mr-glt

mr-glt/mesh_exporter.cpp Secret

Created May 4, 2024 14:54
Show Gist options
  • Save mr-glt/714163f561d032460a6a9ea903684c2a to your computer and use it in GitHub Desktop.
Save mr-glt/714163f561d032460a6a9ea903684c2a to your computer and use it in GitHub Desktop.
Download ZIP
mesh_exporter.cpp
Raw
mesh_exporter.cpp
#include "libvcad/writers/mesh_exporter.h"
// STD
#include <iostream>
#include <vector>
#include "libvcad/utils/marching_cubes_defs.h"
// OMP
#include <omp.h>
using namespace std;
using namespace glm;
struct Cube
{
vec3 p[8];
double dist_vals[8];
uint8_t cube_index;
};
struct Triangle
{
Triangle() {}
Triangle(vec3 p0, vec3 p1, vec3 p2, vec3 n)
{
p[0] = p0;
p[1] = p1;
p[2] = p2;
normal = n;
}
vec3 p[3];
vec3 normal;
};
int GetTriangleValIndex(int i, int j)
{
return TriangleTable[j + (i * 16)];
}
vec3 Interpolate(float iso_value, vec3 p1, vec3 p2, float v1, float v2)
{
float mu = (iso_value - v1) / (v2 - v1);
float x = p1.x + mu * (p2.x - p1.x);
float y = p1.y + mu * (p2.y - p1.y);
float z = p1.z + mu * (p2.z - p1.z);
return vec3(x, y, z);
}
vector<Triangle> TriangulateCubeMaterial(const Cube& cube, double lower_material_iso, double upper_material_iso)
{
std::vector<Triangle> triangles;
// Calculate cube index
int cube_index = 0;
if (cube.dist_vals[0] >= lower_material_iso && cube.dist_vals[0] <= upper_material_iso) cube_index |= 1;
if (cube.dist_vals[1] >= lower_material_iso && cube.dist_vals[1] <= upper_material_iso) cube_index |= 2;
if (cube.dist_vals[2] >= lower_material_iso && cube.dist_vals[2] <= upper_material_iso) cube_index |= 4;
if (cube.dist_vals[3] >= lower_material_iso && cube.dist_vals[3] <= upper_material_iso) cube_index |= 8;
if (cube.dist_vals[4] >= lower_material_iso && cube.dist_vals[4] <= upper_material_iso) cube_index |= 16;
if (cube.dist_vals[5] >= lower_material_iso && cube.dist_vals[5] <= upper_material_iso) cube_index |= 32;
if (cube.dist_vals[6] >= lower_material_iso && cube.dist_vals[6] <= upper_material_iso) cube_index |= 64;
if (cube.dist_vals[7] >= lower_material_iso && cube.dist_vals[7] <= upper_material_iso) cube_index |= 128;
if (EdgeTable[cube_index] == 0)
return triangles;
/* Find the vertices where the surface intersects the cube */
vec3 intersections[12] = { vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0),
vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0),
vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0) };
// Instead of interpolating, just set the coordinates to middle of the edge
// NOTE: this will make things looks less smooth
if ((EdgeTable[cube_index] & 1) != 0)
intersections[0] = (cube.p[0] + cube.p[1]) / 2.0f;
if ((EdgeTable[cube_index] & 2) != 0)
intersections[1] = (cube.p[1] + cube.p[2]) / 2.0f;
if ((EdgeTable[cube_index] & 4) != 0)
intersections[2] = (cube.p[2] + cube.p[3]) / 2.0f;
if ((EdgeTable[cube_index] & 8) != 0)
intersections[3] = (cube.p[3] + cube.p[0]) / 2.0f;
if ((EdgeTable[cube_index] & 16) != 0)
intersections[4] = (cube.p[4] + cube.p[5]) / 2.0f;
if ((EdgeTable[cube_index] & 32) != 0)
intersections[5] = (cube.p[5] + cube.p[6]) / 2.0f;
if ((EdgeTable[cube_index] & 64) != 0)
intersections[6] = (cube.p[6] + cube.p[7]) / 2.0f;
if ((EdgeTable[cube_index] & 128) != 0)
intersections[7] = (cube.p[7] + cube.p[4]) / 2.0f;
if ((EdgeTable[cube_index] & 256) != 0)
intersections[8] = (cube.p[0] + cube.p[4]) / 2.0f;
if ((EdgeTable[cube_index] & 512) != 0)
intersections[9] = (cube.p[1] + cube.p[5]) / 2.0f;
if ((EdgeTable[cube_index] & 1024) != 0)
intersections[10] = (cube.p[2] + cube.p[6]) / 2.0f;
if ((EdgeTable[cube_index] & 2048) != 0)
intersections[11] = (cube.p[3] + cube.p[7]) / 2.0f;
for (int i = 0; GetTriangleValIndex(cube_index, i) != -1; i += 3)
{
vec3 p0 = intersections[GetTriangleValIndex(cube_index, i + 0)];
vec3 p1 = intersections[GetTriangleValIndex(cube_index, i + 1)];
vec3 p2 = intersections[GetTriangleValIndex(cube_index, i + 2)];
auto normal = cross(p1 - p0, p2 - p0);
normal = normalize(normal);
triangles.emplace_back(p0, p1, p2, normal);
}
return triangles;
}
vector<Triangle> TriangulateCubeSurface(const Cube& cube, double iso_value)
{
std::vector<Triangle> triangles;
// Calulate cube index
int cube_index = 0;
if (cube.dist_vals[0] <= iso_value) cube_index |= 1;
if (cube.dist_vals[1] <= iso_value) cube_index |= 2;
if (cube.dist_vals[2] <= iso_value) cube_index |= 4;
if (cube.dist_vals[3] <= iso_value) cube_index |= 8;
if (cube.dist_vals[4] <= iso_value) cube_index |= 16;
if (cube.dist_vals[5] <= iso_value) cube_index |= 32;
if (cube.dist_vals[6] <= iso_value) cube_index |= 64;
if (cube.dist_vals[7] <= iso_value) cube_index |= 128;
if (EdgeTable[cube_index] == 0)
return triangles;
/* Find the vertices where the surface intersects the cube */
vec3 intersections[12] = { vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0),
vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0),
vec3(0,0,0), vec3(0,0,0), vec3(0,0,0), vec3(0,0,0) };
uint8_t colors[12] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
if ((EdgeTable[cube_index] & 1) != 0)
intersections[0] = Interpolate(iso_value, cube.p[0], cube.p[1], cube.dist_vals[0], cube.dist_vals[1]);
if ((EdgeTable[cube_index] & 2) != 0)
intersections[1] = Interpolate(iso_value, cube.p[1], cube.p[2], cube.dist_vals[1], cube.dist_vals[2]);
if ((EdgeTable[cube_index] & 4) != 0)
intersections[2] = Interpolate(iso_value, cube.p[2], cube.p[3], cube.dist_vals[2], cube.dist_vals[3]);
if ((EdgeTable[cube_index] & 8) != 0)
intersections[3] = Interpolate(iso_value, cube.p[3], cube.p[0], cube.dist_vals[3], cube.dist_vals[0]);
if ((EdgeTable[cube_index] & 16) != 0)
intersections[4] = Interpolate(iso_value, cube.p[4], cube.p[5], cube.dist_vals[4], cube.dist_vals[5]);
if ((EdgeTable[cube_index] & 32) != 0)
intersections[5] = Interpolate(iso_value, cube.p[5], cube.p[6], cube.dist_vals[5], cube.dist_vals[6]);
if ((EdgeTable[cube_index] & 64) != 0)
intersections[6] = Interpolate(iso_value, cube.p[6], cube.p[7], cube.dist_vals[6], cube.dist_vals[7]);
if ((EdgeTable[cube_index] & 128) != 0)
intersections[7] = Interpolate(iso_value, cube.p[7], cube.p[4], cube.dist_vals[7], cube.dist_vals[4]);
if ((EdgeTable[cube_index] & 256) != 0)
intersections[8] = Interpolate(iso_value, cube.p[0], cube.p[4], cube.dist_vals[0], cube.dist_vals[4]);
if ((EdgeTable[cube_index] & 512) != 0)
intersections[9] = Interpolate(iso_value, cube.p[1], cube.p[5], cube.dist_vals[1], cube.dist_vals[5]);
if ((EdgeTable[cube_index] & 1024) != 0)
intersections[10] = Interpolate(iso_value, cube.p[2], cube.p[6], cube.dist_vals[2], cube.dist_vals[6]);
if ((EdgeTable[cube_index] & 2048) != 0)
intersections[11] = Interpolate(iso_value, cube.p[3], cube.p[7], cube.dist_vals[3], cube.dist_vals[7]);
for (int i = 0; GetTriangleValIndex(cube_index, i) != -1; i += 3)
{
vec3 p0 = intersections[GetTriangleValIndex(cube_index, i + 0)];
vec3 p1 = intersections[GetTriangleValIndex(cube_index, i + 1)];
vec3 p2 = intersections[GetTriangleValIndex(cube_index, i + 2)];
auto normal = cross(p1 - p0, p2 - p0);
normal = normalize(normal);
triangles.emplace_back(p0, p1, p2, normal);
}
return triangles;
}
void WriteSTLFile(const std::vector<Triangle>& triangles, const std::string& outputFileName) {
std::ofstream stlFile(outputFileName, std::ios::out | std::ios::binary);
if (!stlFile.is_open()) {
std::cerr << "Failed to open STL file for writing." << std::endl;
return;
}
// Write a 80-byte header
std::string header("Custom STL Export");
if (header.size() > 80) {
header = header.substr(0, 80);
} else {
header += std::string(80 - header.size(), ' ');
}
stlFile.write(header.c_str(), 80);
// Write the number of triangles (4 bytes)
uint32_t numTriangles = static_cast<uint32_t>(triangles.size());
stlFile.write(reinterpret_cast<char*>(&numTriangles), sizeof(uint32_t));
// Write each triangle
for (const Triangle& triangle : triangles)
{
// Calculate and write the normal vector (12 bytes, floats)
vec3 normal = normalize(cross(triangle.p[1] - triangle.p[0], triangle.p[2] - triangle.p[0]));
stlFile.write(reinterpret_cast<const char*>(&normal.x), 4);
stlFile.write(reinterpret_cast<const char*>(&normal.y), 4);
stlFile.write(reinterpret_cast<const char*>(&normal.z), 4);
// Write the vertex coordinates (36 bytes, floats)
for (int i = 0; i < 3; i++)
{
stlFile.write(reinterpret_cast<const char*>(&triangle.p[i].x), 4);
stlFile.write(reinterpret_cast<const char*>(&triangle.p[i].y), 4);
stlFile.write(reinterpret_cast<const char*>(&triangle.p[i].z), 4);
}
// Write a short attribute byte count (2 bytes, not used)
uint16_t attributeByteCount = 0;
stlFile.write(reinterpret_cast<char*>(&attributeByteCount), sizeof(uint16_t));
}
stlFile.close();
}
void MeshExporter::exportTreeToMaterialMesh(std::shared_ptr<Node> root, vec3 min, vec3 max, vec3 voxel_size, const string& path,
uint8_t material_chan, double lower_material_iso, double upper_material_iso)
{
float geometry_iso = 0.0;
vec3 space_size(abs(max.x - min.x), abs(max.y - min.y), abs(max.z - min.z));
long x_dim = ceil(space_size.x / voxel_size.x);
long y_dim = ceil(space_size.y / voxel_size.y);
long z_dim = ceil(space_size.z / voxel_size.z);
double half_vox_size_x = voxel_size.x / 2.0;
double half_vox_size_y = voxel_size.y / 2.0;
double half_vox_size_z = voxel_size.z / 2.0;
vector<Triangle> triangles;
// Get total number of voxels
long total_voxels = x_dim * y_dim * z_dim;
// Get start time
auto start = std::chrono::high_resolution_clock::now();
std::cout << "Total voxels: " << total_voxels << std::endl;
// Create a vector of triangles for each omp thread
vector<vector<Triangle>> thread_triangles(omp_get_max_threads());
// Reserve space on the vectors
for (int i = 0; i < thread_triangles.size(); i++)
thread_triangles[i].reserve(total_voxels / 100);
// Make a clone of the tree for each thread
// This is necessary because the tree is not thread safe
vector<shared_ptr<Node>> thread_roots(omp_get_max_threads());
for (int i = 0; i < thread_roots.size(); i++)
thread_roots[i] = root->clone();
#pragma omp parallel for
for (long x = 0; x < x_dim; x++)
{
for (long y = 0; y < y_dim; y++)
{
for (long z = 0; z < z_dim; z++)
{
auto thread_root = thread_roots[omp_get_thread_num()];
vec3 center(min.x + (x * voxel_size.x), min.y + (y * voxel_size.y), min.z + (z * voxel_size.z));
Cube c;
// Extract cube vertices
c.p[0] = vec3(center.x - half_vox_size_x, center.y - half_vox_size_y, center.z - half_vox_size_z);
c.p[1] = vec3(center.x + half_vox_size_x, center.y - half_vox_size_y, center.z - half_vox_size_z);
c.p[2] = vec3(center.x + half_vox_size_x, center.y - half_vox_size_y, center.z + half_vox_size_z);
c.p[3] = vec3(center.x - half_vox_size_x, center.y - half_vox_size_y, center.z + half_vox_size_z);
c.p[4] = vec3(center.x - half_vox_size_x, center.y + half_vox_size_y, center.z - half_vox_size_z);
c.p[5] = vec3(center.x + half_vox_size_x, center.y + half_vox_size_y, center.z - half_vox_size_z);
c.p[6] = vec3(center.x + half_vox_size_x, center.y + half_vox_size_y, center.z + half_vox_size_z);
c.p[7] = vec3(center.x - half_vox_size_x, center.y + half_vox_size_y, center.z + half_vox_size_z);
c.dist_vals[0] = thread_root->evaluate(c.p[0].x, c.p[0].y, c.p[0].z) <= geometry_iso ?
thread_root->distribution(c.p[0].x, c.p[0].y, c.p[0].z).at(material_chan) : -1.0;
c.dist_vals[1] = thread_root->evaluate(c.p[1].x, c.p[1].y, c.p[1].z) <= geometry_iso ?
thread_root->distribution(c.p[1].x, c.p[1].y, c.p[1].z).at(material_chan) : -1.0;
c.dist_vals[2] = thread_root->evaluate(c.p[2].x, c.p[2].y, c.p[2].z) <= geometry_iso ?
thread_root->distribution(c.p[2].x, c.p[2].y, c.p[2].z).at(material_chan) : -1.0;
c.dist_vals[3] = thread_root->evaluate(c.p[3].x, c.p[3].y, c.p[3].z) <= geometry_iso ?
thread_root->distribution(c.p[3].x, c.p[3].y, c.p[3].z).at(material_chan) : -1.0;
c.dist_vals[4] = thread_root->evaluate(c.p[4].x, c.p[4].y, c.p[4].z) <= geometry_iso ?
thread_root->distribution(c.p[4].x, c.p[4].y, c.p[4].z).at(material_chan) : -1.0;
c.dist_vals[5] = thread_root->evaluate(c.p[5].x, c.p[5].y, c.p[5].z) <= geometry_iso ?
thread_root->distribution(c.p[5].x, c.p[5].y, c.p[5].z).at(material_chan) : -1.0;
c.dist_vals[6] = thread_root->evaluate(c.p[6].x, c.p[6].y, c.p[6].z) <= geometry_iso ?
thread_root->distribution(c.p[6].x, c.p[6].y, c.p[6].z).at(material_chan) : -1.0;
c.dist_vals[7] = thread_root->evaluate(c.p[7].x, c.p[7].y, c.p[7].z) <= geometry_iso ?
thread_root->distribution(c.p[7].x, c.p[7].y, c.p[7].z).at(material_chan) : -1.0;
vector<Triangle> local_triangles = TriangulateCubeMaterial(c, lower_material_iso, upper_material_iso);
// Add triangles to the global list in the correct thread slot
std::copy(local_triangles.begin(), local_triangles.end(), std::back_inserter(thread_triangles[omp_get_thread_num()]));
}
}
}
for(auto t : thread_triangles)
std::copy(t.begin(), t.end(), std::back_inserter(triangles));
WriteSTLFile(triangles, path);
}
void MeshExporter::exportTreeToSurfaceMesh(std::shared_ptr<Node> root, vec3 min, vec3 max, vec3 voxel_size,
const std::string& path, double iso_value)
{
vec3 space_size(abs(max.x - min.x), abs(max.y - min.y), abs(max.z - min.z));
long x_dim = ceil(space_size.x / voxel_size.x);
long y_dim = ceil(space_size.y / voxel_size.y);
long z_dim = ceil(space_size.z / voxel_size.z);
double half_vox_size_x = voxel_size.x / 2.0;
double half_vox_size_y = voxel_size.y / 2.0;
double half_vox_size_z = voxel_size.z / 2.0;
vector<Triangle> triangles;
// Get total number of voxels
long total_voxels = x_dim * y_dim * z_dim;
// Get start time
auto start = std::chrono::high_resolution_clock::now();
std::cout << "Total voxels: " << total_voxels << std::endl;
// Create a vector of triangles for each omp thread
vector<vector<Triangle>> thread_triangles(omp_get_max_threads());
// Reserve space on the vectors
for (int i = 0; i < thread_triangles.size(); i++)
thread_triangles[i].reserve(total_voxels / 100);
// Make a clone of the tree for each thread
// This is necessary because the tree is not thread safe
vector<shared_ptr<Node>> thread_roots(omp_get_max_threads());
for (int i = 0; i < thread_roots.size(); i++)
thread_roots[i] = root->clone();
#pragma omp parallel for
for (long x = 0; x < x_dim; x++)
{
for (long y = 0; y < y_dim; y++)
{
for (long z = 0; z < z_dim; z++)
{
auto thread_root = thread_roots[omp_get_thread_num()];
vec3 center(min.x + (x * voxel_size.x), min.y + (y * voxel_size.y), min.z + (z * voxel_size.z));
Cube c;
// Extract cube vertices
c.p[0] = vec3(center.x - half_vox_size_x, center.y - half_vox_size_y, center.z - half_vox_size_z);
c.p[1] = vec3(center.x + half_vox_size_x, center.y - half_vox_size_y, center.z - half_vox_size_z);
c.p[2] = vec3(center.x + half_vox_size_x, center.y - half_vox_size_y, center.z + half_vox_size_z);
c.p[3] = vec3(center.x - half_vox_size_x, center.y - half_vox_size_y, center.z + half_vox_size_z);
c.p[4] = vec3(center.x - half_vox_size_x, center.y + half_vox_size_y, center.z - half_vox_size_z);
c.p[5] = vec3(center.x + half_vox_size_x, center.y + half_vox_size_y, center.z - half_vox_size_z);
c.p[6] = vec3(center.x + half_vox_size_x, center.y + half_vox_size_y, center.z + half_vox_size_z);
c.p[7] = vec3(center.x - half_vox_size_x, center.y + half_vox_size_y, center.z + half_vox_size_z);
c.dist_vals[0] = thread_root->evaluate(c.p[0].x, c.p[0].y, c.p[0].z);
c.dist_vals[1] = thread_root->evaluate(c.p[1].x, c.p[1].y, c.p[1].z);
c.dist_vals[2] = thread_root->evaluate(c.p[2].x, c.p[2].y, c.p[2].z);
c.dist_vals[3] = thread_root->evaluate(c.p[3].x, c.p[3].y, c.p[3].z);
c.dist_vals[4] = thread_root->evaluate(c.p[4].x, c.p[4].y, c.p[4].z);
c.dist_vals[5] = thread_root->evaluate(c.p[5].x, c.p[5].y, c.p[5].z);
c.dist_vals[6] = thread_root->evaluate(c.p[6].x, c.p[6].y, c.p[6].z);
c.dist_vals[7] = thread_root->evaluate(c.p[7].x, c.p[7].y, c.p[7].z);
vector<Triangle> local_triangles = TriangulateCubeSurface(c, iso_value);
// Add triangles to the global list in the correct thread slot
std::copy(local_triangles.begin(), local_triangles.end(), std::back_inserter(thread_triangles[omp_get_thread_num()]));
}
}
}
for(auto t : thread_triangles)
std::copy(t.begin(), t.end(), std::back_inserter(triangles));
// Convert to indexed triangle format
std::vector<Point_3> vertices;
std::vector<Indexed_Triangle> index_triangles;
for (const auto& triangle : triangles)
for (auto i : triangle.p) vertices.push_back(Point_3(i.x, i.y, i.z));
for (size_t i = 0; i < triangles.size(); i++) index_triangles.push_back(Indexed_Triangle{i * 3, i * 3 + 1, i * 3 + 2});
// Create a surface mesh
SurfaceMesh mesh(vertices, index_triangles);
mesh.write_stl(path);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment