From 5defc9b468e6df08e8b5aa078a4fa23be4dd8061 Mon Sep 17 00:00:00 2001 From: Dtouch Date: Mon, 4 Nov 2024 20:13:02 +0800 Subject: [PATCH] minimal implementation --- .gitignore | 9 ++ CMakeLists.txt | 8 ++ Geometry.hpp | 177 +++++++++++++++++++++++++++++ Intersection.hpp | 281 +++++++++++++++++++++++++++++++++++++++++++++++ STLReader.hpp | 66 +++++++++++ Voxel.hpp | 85 ++++++++++++++ main.cpp | 24 ++++ 7 files changed, 650 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 Geometry.hpp create mode 100644 Intersection.hpp create mode 100644 STLReader.hpp create mode 100644 Voxel.hpp create mode 100644 main.cpp diff --git a/.gitignore b/.gitignore index e257658..c4c0fb5 100644 --- a/.gitignore +++ b/.gitignore @@ -32,3 +32,12 @@ *.out *.app +*.stl + +*-build-*/ +build/ +Build/ + +.idea/ +.vs/ +.vscode/ diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..4a25dc9 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,8 @@ +cmake_minimum_required(VERSION 3.27) +project(STL2Voxel) + +set(CMAKE_CXX_STANDARD 20) + +include_directories(.) + +add_executable(STL2Voxel main.cpp) diff --git a/Geometry.hpp b/Geometry.hpp new file mode 100644 index 0000000..115b90f --- /dev/null +++ b/Geometry.hpp @@ -0,0 +1,177 @@ +#pragma once +#include +#include + +template +class Vec3 { +public: + Vec3(T a_, T b_, T c_) : a(a_), b(b_), c(c_) {} + Vec3() : a(0), b(0), c(0) {} + + union { T x, a; }; + union { T y, b; }; + union { T z, c; }; + + T operator[](int i) const { + if (i == 0) { + return a; + } else if (i == 1) { + return b; + } else if (i == 2) { + return c; + } + return a; + } + + T &operator[](int i) { + if (i == 0) { + return a; + } else if (i == 1) { + return b; + } else if (i == 2) { + return c; + } + return a; + } + + Vec3 operator+(const Vec3 &other) const { + return {a + other.a, b + other.b, c + other.c}; + } + + Vec3 operator/(const float w) const { + return {a / w, b / w, c / w}; + } + + Vec3 operator-(const Vec3 &other) const { + return {a - other.a, b - other.b, c - other.c}; + } + + Vec3 operator*(const float &w) const { + return {a * w, b * w, c * w}; + } + + float length() const { + return sqrt(a * a + b * b + c * c); + } + + Vec3 norm() const { + return *this / length(); + } + + Vec3 cross(const Vec3 &other) const { + return {b * other.c - c * other.b, c * other.a - a * other.c, a * other.b - b * other.a}; + } + + float dot(const Vec3 &other) const { + return a * other.a + b * other.b + c * other.c; + } + + bool operator==(const Vec3& other) const { + return x == other.x && y == other.y && z == other.z; + } +}; + +// 为 Vec3 创建哈希函数 +namespace std { + template + struct hash> { + size_t operator()(const Vec3& vec) const { + return hash()(vec.x) ^ hash()(vec.y) ^ hash()(vec.z); + } +}; +} + +typedef Vec3 Vec3f; +typedef Vec3 Vec3u; + +class AABB { +public: + Vec3f min = {std::numeric_limits::max(), std::numeric_limits::max(), + std::numeric_limits::max()}; + Vec3f max = {std::numeric_limits::lowest(), std::numeric_limits::lowest(), + std::numeric_limits::lowest()}; + Vec3f center() const { + return (min + max) / 2; + } + Vec3f size() const { + return max - min; + } + AABB operator+(const Vec3f& v) const { + return {min + v, max + v}; + } +}; + +class Mesh { +public: + std::vector vertices; + std::vector faces; + // 获得指定面的 AABB + AABB getAABB(const std::vector& faceIndices) const { + AABB aabb; + for (const size_t &idx: faceIndices) { + const Vec3u &face = faces[idx]; + for (int i = 0; i < 3; ++i) { + const Vec3f &v = vertices[face[i]]; + for (int j = 0; j < 3; ++j) { + aabb.min[j] = std::min(aabb.min[j], v[j]); + aabb.max[j] = std::max(aabb.max[j], v[j]); + } + } + } + return aabb; + } + AABB getAABB() const { + AABB aabb; + for (const Vec3f &v: vertices) { + for (int j = 0; j < 3; ++j) { + aabb.min[j] = std::min(aabb.min[j], v[j]); + aabb.max[j] = std::max(aabb.max[j], v[j]); + } + } + return aabb; + } +}; + + +//class LineSegment { +//public: +// Vec3f start; +// Vec3f end; +// +// /** +// * @brief Construct a new Line Segment object +// * @param s start point +// * @param e end point +// * @param r radius +// */ +// LineSegment(const Vec3f &s, const Vec3f &e, const float r = 0) : start(s), end(e) { +// length = (end - start).length(); +// dir = (end - start).norm(); +// } +// +// float getLength() const { +// return length; +// } +// +// Vec3f getDir() const { +// return dir; +// } +// +//private: +// float length; +// Vec3f dir; +//}; + +class Ray { +public: + Vec3f start; + Vec3f dir; + + Ray(const Vec3f &s, const Vec3f &d) : start(s), dir(d) { + dir = dir.norm(); + } +}; + + + + diff --git a/Intersection.hpp b/Intersection.hpp new file mode 100644 index 0000000..8eb0e4d --- /dev/null +++ b/Intersection.hpp @@ -0,0 +1,281 @@ +#pragma once + +#include +#include "cmath" +#include "algorithm" +#include "Geometry.hpp" + +struct IntersectRes { + bool hit; + float t; +}; + +IntersectRes +triangleRayIntersection(const Ray &ray, const Vec3f &a, const Vec3f &b, const Vec3f &c) { + Vec3f e1 = b - a; + Vec3f e2 = c - a; + Vec3f s = ray.start - a; + Vec3f s1 = ray.dir.cross(e2); + Vec3f s2 = s.cross(e1); + float s1_dot_e1 = s1.dot(e1); + float b1 = s1.dot(s) / s1_dot_e1; + float b2 = s2.dot(ray.dir) / s1_dot_e1; + if (b1 >= 0. && b2 >= 0. && b1 + b2 <= 1.) { + // hit + return {true, s2.dot(e2) / s1_dot_e1}; + } + return {false, 0.}; +} + +bool isClose(double a, double b, double rel_tol = 1e-9, double abs_tol = 0.0) { + return std::fabs(a - b) <= std::max(rel_tol * std::max(std::fabs(a), std::fabs(b)), abs_tol); +} + +int sign(double x) { + return (x > 0) - (x < 0); +} + +float dot2(Vec3f a) { + return a.dot(a); +} + +float pointTriangleDistance(const Vec3f &p, const Vec3f &a, const Vec3f &b, const Vec3f &c) { + Vec3f ab = b - a, ap = p - a; + Vec3f bc = c - b, bp = p - b; + Vec3f ca = a - c, cp = p - c; + Vec3f nor = ab.cross(ca); + + return sqrt( + (sign(ab.cross(nor).dot(ap)) + + sign(bc.cross(nor).dot(bp)) + + sign(ca.cross(nor).dot(cp)) < 2.0) + ? + fmin(fmin( + dot2(ab * std::clamp(ab.dot(ap) / dot2(ab), 0.0f, 1.0f) - ap), + dot2(bc * std::clamp(bc.dot(bp) / dot2(bc), 0.0f, 1.0f) - bp)), + dot2(ca * std::clamp(ca.dot(cp) / dot2(ca), 0.0f, 1.0f) - cp)) + : + nor.dot(ap) * nor.dot(ap) / dot2(nor)); +} + +class BVHNode { +public: + size_t left; + size_t right; + size_t parent; + AABB boundingBox; + + BVHNode(size_t l, size_t r, size_t p, const AABB &aabb) : left(l), right(r), parent(p), boundingBox(aabb) {} + + BVHNode() : left(0), right(0), parent(0), boundingBox() {} +}; + + +class BVH { +public: + const Mesh &mesh; + std::vector faceCenters; + std::vector nodes; + + BVH(const Mesh &mesh_) : mesh(mesh_) { + faceCenters.reserve(mesh.faces.size()); + for (const Vec3u &face: mesh.faces) { + Vec3f center = (mesh.vertices[face[0]] + mesh.vertices[face[1]] + mesh.vertices[face[2]]) / 3.0f; + faceCenters.push_back(center); + } + std::vector indicesList(mesh.faces.size()); + std::iota(indicesList.begin(), indicesList.end(), 0); + nodes.resize(2 * mesh.faces.size() - 1); + size_t nowIdx = 0; + dfsBuild(indicesList, computeAABB(indicesList), nowIdx); + } + + size_t + dfsBuild(std::vector &indicesList, AABB aabb, size_t &nowIdx) { + const size_t nodeIdx = nowIdx; + nowIdx++; + if (indicesList.size() == 1) { + // leaf + nodes[nodeIdx] = {0, indicesList[0], 0, aabb}; + return nodeIdx; + } + // longest axis + int longAxis = -1; + float longAxisLen = -1; + for (int i = 0; i < 3; ++i) { + const float axisLen = aabb.max[i] - aabb.min[i]; + if (axisLen > longAxisLen) { + longAxisLen = axisLen; + longAxis = i; + } + } + + // split indices list + const size_t k = indicesList.size() / 2; +// std::nth_element(indicesList.begin(), indicesList.begin() + k - 1, indicesList.end(), +// [&](const size_t &a, const size_t &b) { +// return faceCenters[a][longAxis] < faceCenters[b][longAxis]; +// }); + nth(indicesList, k - 1, longAxis); + std::vector leftIndices(indicesList.begin(), indicesList.begin() + k); + std::vector rightIndices(indicesList.begin() + k, indicesList.end()); + + const AABB leftAABB = computeAABB(leftIndices); + const AABB rightAABB = computeAABB(rightIndices); + const size_t leftIdx = dfsBuild(leftIndices, leftAABB, nowIdx); + const size_t rightIdx = dfsBuild(rightIndices, rightAABB, nowIdx); + nodes[nodeIdx] = {leftIdx, rightIdx, 0, aabb}; + nodes[leftIdx].parent = nodeIdx; + nodes[rightIdx].parent = nodeIdx; + return nodeIdx; + } + + unsigned int intersectWithRay(const Ray &ray) const { + return recursiveRayIntersection(ray, 0); + } + + void writeBVHAsObj(const std::string &path) { + std::ofstream file(path); + for (const BVHNode &node: nodes) { + Vec3f vertex; + for (int i = 0; i < 2; ++i) { + vertex[0] = i == 0 ? node.boundingBox.min[0] : node.boundingBox.max[0]; + for (int j = 0; j < 2; ++j) { + vertex[1] = j == 0 ? node.boundingBox.min[1] : node.boundingBox.max[1]; + for (int k = 0; k < 2; ++k) { + vertex[2] = k == 0 ? node.boundingBox.min[2] : node.boundingBox.max[2]; + file << "v " << vertex[0] << " " << vertex[1] << " " << vertex[2] << std::endl; + } + } + } + } + for (size_t i = 0; i < nodes.size(); ++i) { + const BVHNode &node = nodes[i]; + + file << "l " << i * 8 + 1 << " " << i * 8 + 2 << std::endl; + file << "l " << i * 8 + 3 << " " << i * 8 + 4 << std::endl; + file << "l " << i * 8 + 5 << " " << i * 8 + 6 << std::endl; + file << "l " << i * 8 + 7 << " " << i * 8 + 8 << std::endl; + file << "l " << i * 8 + 1 << " " << i * 8 + 3 << std::endl; + file << "l " << i * 8 + 2 << " " << i * 8 + 4 << std::endl; + file << "l " << i * 8 + 5 << " " << i * 8 + 7 << std::endl; + file << "l " << i * 8 + 6 << " " << i * 8 + 8 << std::endl; + file << "l " << i * 8 + 1 << " " << i * 8 + 5 << std::endl; + file << "l " << i * 8 + 2 << " " << i * 8 + 6 << std::endl; + file << "l " << i * 8 + 3 << " " << i * 8 + 7 << std::endl; + file << "l " << i * 8 + 4 << " " << i * 8 + 8 << std::endl; + + } + file.close(); + } + + void writeLeavesAsObj(const std::string &path) { + std::ofstream file(path); + for (const BVHNode &node: nodes) { + if (node.left == 0) { + Vec3f vertex; + for (int i = 0; i < 2; ++i) { + vertex[0] = i == 0 ? node.boundingBox.min[0] : node.boundingBox.max[0]; + for (int j = 0; j < 2; ++j) { + vertex[1] = j == 0 ? node.boundingBox.min[1] : node.boundingBox.max[1]; + for (int k = 0; k < 2; ++k) { + vertex[2] = k == 0 ? node.boundingBox.min[2] : node.boundingBox.max[2]; + file << "v " << vertex[0] << " " << vertex[1] << " " << vertex[2] << std::endl; + } + } + } + } + } + for (size_t i = 0; i < nodes.size(); ++i) { + const BVHNode &node = nodes[i]; + if (node.left == 0) { + file << "l " << i * 8 + 1 << " " << i * 8 + 2 << std::endl; + file << "l " << i * 8 + 3 << " " << i * 8 + 4 << std::endl; + file << "l " << i * 8 + 5 << " " << i * 8 + 6 << std::endl; + file << "l " << i * 8 + 7 << " " << i * 8 + 8 << std::endl; + file << "l " << i * 8 + 1 << " " << i * 8 + 3 << std::endl; + file << "l " << i * 8 + 2 << " " << i * 8 + 4 << std::endl; + file << "l " << i * 8 + 5 << " " << i * 8 + 7 << std::endl; + file << "l " << i * 8 + 6 << " " << i * 8 + 8 << std::endl; + file << "l " << i * 8 + 1 << " " << i * 8 + 5 << std::endl; + file << "l " << i * 8 + 2 << " " << i * 8 + 6 << std::endl; + file << "l " << i * 8 + 3 << " " << i * 8 + 7 << std::endl; + file << "l " << i * 8 + 4 << " " << i * 8 + 8 << std::endl; + } + } + file.close(); + } + +private: + AABB computeAABB(const std::vector &indices) { + AABB aabb; + for (const size_t &idx: indices) { + const Vec3u &face = mesh.faces[idx]; + for (int i = 0; i < 3; ++i) { + const Vec3f &vertex = mesh.vertices[face[i]]; + for (int j = 0; j < 3; ++j) { + aabb.min[j] = std::min(aabb.min[j], vertex[j]); + aabb.max[j] = std::max(aabb.max[j], vertex[j]); + } + } + } + return {aabb.min, aabb.max}; + } + + unsigned int recursiveRayIntersection(const Ray &ray, size_t nodeIdx) const { + // segment-box intersection test + const AABB &aabb = nodes[nodeIdx].boundingBox; + bool hit = false; + for (int i = 0; !hit && i < 3; ++i) { + float t_min = (aabb.min[i] - ray.start[i]) / ray.dir[i]; + float t_max = (aabb.max[i] - ray.start[i]) / ray.dir[i]; + if (t_min > t_max) { + std::swap(t_min, t_max); + } + if (t_max < 0) return 0; + Vec3f hitPt = ray.start + ray.dir * t_min; + int otherPlane1 = (i + 1) % 3, otherPlane2 = (i + 2) % 3; + if (hitPt[otherPlane1] >= aabb.min[otherPlane1] && + hitPt[otherPlane1] <= aabb.max[otherPlane1] && + hitPt[otherPlane2] >= aabb.min[otherPlane2] && + hitPt[otherPlane2] <= aabb.max[otherPlane2]) { + hit = true; + } + } + if (hit) { + if (nodes[nodeIdx].left == 0) { + // leaf + Vec3u face = mesh.faces[nodes[nodeIdx].right]; + IntersectRes res = triangleRayIntersection(ray, mesh.vertices[face[0]], mesh.vertices[face[1]], + mesh.vertices[face[2]]); + if (!res.hit || res.t < 0) return 0; + return 1; + } else { + // check children + return recursiveRayIntersection(ray, nodes[nodeIdx].left) + + recursiveRayIntersection(ray, nodes[nodeIdx].right); + } + } + return 0; + } + + // implement nth, without std::nth_element + // kth is 0-based + void nth(std::vector &indicesList, size_t kth, int longAxis) { + recursiveChoose(indicesList, 0, indicesList.size() - 1, kth, longAxis); + } + + void recursiveChoose(std::vector &indicesList, size_t begin, size_t end, size_t kth, int longAxis) { + if (begin >= end) return; + int i = begin, j = end; + while (i < j) { + while (i < j && faceCenters[indicesList[j]][longAxis] >= faceCenters[indicesList[begin]][longAxis]) j--; + while (i < j && faceCenters[indicesList[i]][longAxis] <= faceCenters[indicesList[begin]][longAxis]) i++; + std::swap(indicesList[i], indicesList[j]); + } + std::swap(indicesList[begin], indicesList[i]); + if (i == kth) return; + if (i < kth) recursiveChoose(indicesList, i + 1, end, kth, longAxis); + else recursiveChoose(indicesList, begin, i - 1, kth, longAxis); + } +}; \ No newline at end of file diff --git a/STLReader.hpp b/STLReader.hpp new file mode 100644 index 0000000..105d601 --- /dev/null +++ b/STLReader.hpp @@ -0,0 +1,66 @@ +#pragma once +#include "Geometry.hpp" +#include +#include +#include +#include +#include +#include + +bool readSTL(const std::string& filename, Mesh& mesh) { + std::ifstream file(filename); + if (!file.is_open()) { + std::cerr << "Failed to open file: " << filename << std::endl; + return false; + } + + std::string line; + std::unordered_map, unsigned int> vertexMap; // 用于去重 + unsigned int vertexIndex = 0; + + while (std::getline(file, line)) { + std::istringstream stream(line); + std::string keyword; + + if (line.find("facet normal") != std::string::npos) { + std::getline(file, line); // Read the outer loop + + Vec3 vertices[3]; + for (int i = 0; i < 3; ++i) { + std::getline(file, line); + stream.clear(); + stream.str(line); + stream >> keyword >> vertices[i].x >> vertices[i].y >> vertices[i].z; + + // 检查是否已存在该顶点 + if (vertexMap.find(vertices[i]) == vertexMap.end()) { + vertexMap[vertices[i]] = vertexIndex; + mesh.vertices.push_back(vertices[i]); + vertexIndex++; + } + } + + // 添加面索引 + mesh.faces.emplace_back(vertexMap[vertices[0]], vertexMap[vertices[1]], vertexMap[vertices[2]]); + + std::getline(file, line); // Read endloop + std::getline(file, line); // Read endfacet + } + } + + file.close(); + return true; +} + +//int main() { +// Mesh mesh; +// if (readSTL("../plate.stl", mesh)) { +// std::cout << "Successfully read STL file." << std::endl; +// std::cout << "Vertices count: " << mesh.vertices.size() << std::endl; +// std::cout << "Faces count: " << mesh.faces.size() << std::endl; +// } else { +// std::cerr << "Failed to read STL file." << std::endl; +// } +// return 0; +//} + diff --git a/Voxel.hpp b/Voxel.hpp new file mode 100644 index 0000000..3f3fbc3 --- /dev/null +++ b/Voxel.hpp @@ -0,0 +1,85 @@ +#pragma once + +#include +#include "Intersection.hpp" + +class VoxelGrid { +private: + int width, height, depth, xySize; + Vec3 cellSize; + AABB aabb; + std::vector data; // 每个体素只占用1 bit + +public: + VoxelGrid(const AABB &aabb, const Vec3 &cellSize) + : cellSize(cellSize) { + // 根据 AABB 和 cellSize 计算宽、高、深 + width = static_cast(std::lround((aabb.max.x - aabb.min.x) / cellSize.x)); + height = static_cast(std::lround((aabb.max.y - aabb.min.y) / cellSize.y)); + depth = static_cast(std::lround((aabb.max.z - aabb.min.z) / cellSize.z)); + + xySize = width * height; + int totalSize = (xySize * depth + 7) / 8; // 每8个体素占用一个字节 + data.resize(totalSize, 0); + this->aabb = aabb; // 使用传入的 AABB + } + + VoxelGrid() = default; + + // 根据索引获取单元格的 AABB + AABB getCellAABB(int x, int y, int z) const { + Vec3 cellMin(x * cellSize.x, y * cellSize.y, z * cellSize.z); + Vec3 cellMax((x + 1) * cellSize.x, (y + 1) * cellSize.y, (z + 1) * cellSize.z); + return AABB(cellMin, cellMax) + aabb.min; + } + + // 根据索引写入 + void setVoxel(int x, int y, int z, bool value) { + assert(x >= 0 && x < width && y >= 0 && y < height && z >= 0 && z < depth); + int index = x + y * width + z * xySize; + int byteIndex = index / 8; + int bitIndex = index % 8; + + if (value) { + data[byteIndex] |= (1 << bitIndex); // 设置为1 + } else { + data[byteIndex] &= ~(1 << bitIndex); // 设置为0 + } + } + + // 根据索引读取 + bool getVoxel(int x, int y, int z) const { + assert(x >= 0 && x < width && y >= 0 && y < height && z >= 0 && z < depth); + int index = x + y * width + z * xySize; + int byteIndex = index / 8; + int bitIndex = index % 8; + + return (data[byteIndex] >> bitIndex) & 1; // 返回位值 + } + + static VoxelGrid generateFromMesh(const Mesh &mesh, const Vec3 &cellSize) { + auto sceneAABB = mesh.getAABB(); + auto voxelGrid = VoxelGrid(sceneAABB, cellSize); + BVH bvh(mesh); + int solvedVoxel = 0; + for (int z = 0; z < voxelGrid.depth; ++z) { + for (int y = 0; y < voxelGrid.height; ++y) { + for (int x = 0; x < voxelGrid.width; ++x) { + auto start = voxelGrid.getCellAABB(x, y, z).center(); + auto dir = Vec3f(1, 2, 3).norm(); + if (bvh.intersectWithRay({start, dir}) % 2 == 1) { + voxelGrid.setVoxel(x, y, z, true); + } else { + voxelGrid.setVoxel(x, y, z, false); + } + + solvedVoxel++; + if (solvedVoxel % int(1e5) == 0) { + std::cout << "Solved " << solvedVoxel << " voxels." << std::endl; + } + } + } + } + return voxelGrid; + } +}; \ No newline at end of file diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..7385201 --- /dev/null +++ b/main.cpp @@ -0,0 +1,24 @@ +#include +#include "STLReader.hpp" +#include "Voxel.hpp" +#include + +int main() { + Mesh mesh; + if (readSTL("../plate.stl", mesh)) { + std::cout << "Successfully read STL file." << std::endl; + std::cout << "Vertices count: " << mesh.vertices.size() << std::endl; + std::cout << "Faces count: " << mesh.faces.size() << std::endl; + } else { + std::cerr << "Failed to read STL file." << std::endl; + } + Vec3f cellSize = {0.1, 0.1, 0.1}; +// Vec3 cellSize = {0.01, 0.01, 0.01}; + // timer + auto begin = std::chrono::steady_clock::now(); + auto voxelGrid = VoxelGrid::generateFromMesh(mesh, cellSize); + auto end = std::chrono::steady_clock::now(); + std::cout << "Time elapsed: " << std::chrono::duration_cast(end - begin).count() << "ms" + << std::endl; + return 0; +}