add voxel func

CMakeLists.txt

@@ -35,6 +35,7 @@ endif()
 
 if(Octree_BUILD_TEST)
     add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/octree_test)
+    add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/solid_octree_test)
     add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/sdf_test)
     add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/udf_test)
     add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/search_assign_test)
@@ -45,6 +46,7 @@ endif()
 enable_testing()
 if(Octree_BUILD_TEST)
     add_test(NAME octree_test COMMAND $<TARGET_FILE:octree_test>)
+    add_test(NAME solid_octree_test COMMAND $<TARGET_FILE:solid_octree_test>)
     add_test(NAME sdf_test COMMAND $<TARGET_FILE:sdf_test>)
     add_test(NAME udf_test COMMAND $<TARGET_FILE:udf_test>)
     add_test(NAME search_assign_test COMMAND $<TARGET_FILE:search_assign_test>)

include/Octree/AABB.h

@@ -15,10 +15,10 @@
 #include <array>
 
 namespace Octree {
-    class AABB {
+class AABB {
     Eigen::Vector3d _min, _max;
 
-    public:
+public:
     explicit AABB(const Eigen::Vector3d &min, const Eigen::Vector3d &max);
 
     explicit AABB(double minx, double maxx, double miny, double maxy, double minz, double maxz);
@@ -27,7 +27,7 @@ namespace Octree {
 
     explicit AABB() = default;
 
-    public:
+public:
     Eigen::Vector3d size() const;
 
     Eigen::Vector3d min() const;
@@ -63,20 +63,34 @@ namespace Octree {
                                  const Eigen::Vector3d &b,
                                  const Eigen::Vector3d &c) const;
 
+    /**
+     * If the AABB intersects with another one
+     * @param aabb
+     * @return
+     */
+    bool intersect_aabb(const AABB &aabb) const;
+
     /**
      * Divide this AABB into 8 parts
      * @return
      */
     std::array<AABB, 8> half_divide() const;
 
-    public:
+public:
 
     /**
      * Extend boundary with @param step
      * @param step
      */
     void extend(double step);
-    };
+
+public:
+    static AABB ToCubeAABBWithSize(const AABB &aabb, const double res);
+
+    static AABB ToCubeAABBWithSize(const AABB &aabb, const double res, int &num_on_axis);
+
+    static AABB ToIntegerSizeAABB(const AABB &aabb, const double res, int &num_on_x, int &num_on_y, int &num_on_z);
+};
 }
 
 #endif //OCTREE_AABB_H

include/Octree/BaseOctree.h

@@ -19,40 +19,39 @@
 #include "AABB.h"
 
 namespace Octree {
-    class BaseOctree;
+class BaseOctree;
 
-    class OctreeNode;
+class OctreeNode;
 
-    using std::shared_ptr;
+using std::shared_ptr;
-    using std::weak_ptr;
+using std::weak_ptr;
-    using std::vector;
+using std::vector;
-    using std::set;
+using std::set;
-    using std::array;
+using std::array;
 
 #if 0
-    using octree_sptr = shared_ptr<SDFOctree>;
+using octree_sptr = shared_ptr<SDFOctree>;
-    using octree_wptr = weak_ptr<SDFOctree>;
+using octree_wptr = weak_ptr<SDFOctree>;
 #endif
 
-    using ocnode_sptr = shared_ptr<OctreeNode>;
+using ocnode_sptr = shared_ptr<OctreeNode>;
-    using ocnode_wptr = weak_ptr<OctreeNode>;
+using ocnode_wptr = weak_ptr<OctreeNode>;
-    using ocnode_ptr = OctreeNode *;
+using ocnode_ptr = OctreeNode *;
 
-
+struct OctreeNode : public std::enable_shared_from_this<OctreeNode> {
-    struct OctreeNode : public std::enable_shared_from_this<OctreeNode> {
     friend class BaseOctree;
     friend class OctreeBuilder;
 
-    protected:
+protected:
     uint32_t level;
     array<ocnode_sptr, 8> children;
     ocnode_ptr father;
     AABB aabb;
 
-    public:
+public:
     set<uint32_t> tri_ids;
 
-    public:
+public:
     explicit OctreeNode(ocnode_ptr father, const AABB &aabb, uint32_t level);
 
     /**
@@ -93,24 +92,26 @@ namespace Octree {
      */
     ocnode_sptr map_node(const Eigen::Vector3d &pos);
 
-    public:  //Getter
+    void map_voxel(const AABB &aabb, std::set<ocnode_sptr> &nodes);
-        uint32_t get_level(){ return level; }
+
+public:  //Getter
+    uint32_t get_level() { return level; }
 
     AABB get_aabb() { return aabb; }
-    };
+};
 
-    class BaseOctree {
+class BaseOctree {
-    public:
+public:
     uint32_t max_level;
 
     ocnode_sptr root;
 
-    public:
+public:
     explicit BaseOctree(uint32_t max_level, const AABB &aabb);
 
     explicit BaseOctree(uint32_t max_level, ocnode_sptr root);
 
-    public:
+public:
     /**
      * Find the leaf that contains the @param pos
      * @param pos
@@ -118,6 +119,13 @@ namespace Octree {
      */
     ocnode_sptr map_node(const Eigen::Vector3d &pos);
 
+    /**
+     * Find the leaves that intersect with @param aabb
+     * @param aabb
+     * @return
+     */
+    std::set<ocnode_sptr> map_voxel(const AABB &aabb);
+
     /**
      * Return the root node
      * @return
@@ -135,8 +143,7 @@ namespace Octree {
      * @return
      */
     AABB aabb() const;
-    };
+};
 }
 
-
 #endif //OCTREE_BASEOCTREE_H

include/Octree/OctreeBuilder.h

@@ -15,7 +15,7 @@
 #include <Eigen/Core>
 #include <pMesh/mesh/TriangleMesh.h>
 #include "AABB.h"
-
+#include "BaseOctree.h"
 
 namespace Octree {
 

include/Octree/voxel/SolidVoxelBuilder.h

@@ -0,0 +1,21 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#ifndef OCTREE_SRC_SOLIDVOXELBUILDER_H_
+#define OCTREE_SRC_SOLIDVOXELBUILDER_H_
+
+#include "Octree/voxel/Voxel.h"
+
+namespace Octree {
+class SolidVoxelBuilder {
+
+};
+}
+
+#endif //OCTREE_SRC_SOLIDVOXELBUILDER_H_

include/Octree/voxel/SurfaceVoxelBuilder.h

@@ -0,0 +1,35 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#ifndef OCTREE_INCLUDE_OCTREE_VOXEL_SURFACEVOXELBUILDER_H_
+#define OCTREE_INCLUDE_OCTREE_VOXEL_SURFACEVOXELBUILDER_H_
+
+
+#include "Octree/BaseOctree.h"
+#include "Octree/OctreeBuilder.h"
+
+namespace Octree {
+class Voxel;
+
+class SurfaceVoxelBuilder {
+private:
+    AABB aabb;
+    int _x, _y, _z;
+    double block_size;
+    BaseOctree &octree;
+    const pMesh::Triangle3dMesh<> &mesh;
+
+public:
+    explicit SurfaceVoxelBuilder(BaseOctree &octree, const pMesh::Triangle3dMesh<> &mesh, double block_size);
+
+    void build(Voxel &voxel);
+};
+}
+
+#endif //OCTREE_INCLUDE_OCTREE_VOXEL_SURFACEVOXELBUILDER_H_

include/Octree/voxel/Voxel.h

@@ -0,0 +1,42 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#ifndef OCTREE_INCLUDE_OCTREE_VOXEL_VOXEL_H_
+#define OCTREE_INCLUDE_OCTREE_VOXEL_VOXEL_H_
+
+#include <vector>
+#include <Octree/AABB.h>
+
+namespace Octree {
+class Voxel {
+public:
+    AABB aabb;
+    int _x, _y, _z;
+    std::vector<bool> voxels;
+
+public:
+    Voxel() = default;
+    Voxel(int _x, int _y, int _z);
+    Voxel(const std::vector<bool> &value, int _x, int _y, int _z);
+
+public:
+    /**
+     * Tool function: calc 1D index from 3D coordinate
+     * @param x
+     * @param y
+     * @param z
+     * @return
+     */
+    int id_at(int x, int y, int z) const;
+
+    bool at(int x, int y, int z) const;
+};
+}
+
+#endif //OCTREE_INCLUDE_OCTREE_VOXEL_VOXEL_H_

src/AABB.cpp

@@ -15,23 +15,23 @@
 /********************* TOOLS ************************************/
 
 namespace AABBTools {
-    static bool AXISTEST(double rad, double p0, double p1) {
+static bool AXISTEST(double rad, double p0, double p1) {
     return (std::min(p0, p1) > rad || std::max(p0, p1) < -rad);
-    }
+}
 
-    inline double vec_min_component(double a, double b, double c) {
+inline double vec_min_component(double a, double b, double c) {
     return std::min(std::min(a, b), c);
-    }
+}
 
-    inline double vec_max_component(double a, double b, double c) {
+inline double vec_max_component(double a, double b, double c) {
     return std::max(std::max(a, b), c);
-    }
+}
 
-    bool directionTest(double a, double b, double c, double hs) {
+bool directionTest(double a, double b, double c, double hs) {
     return (vec_min_component(a, b, c) > hs || vec_max_component(a, b, c) < -hs);
-    }
+}
 
-    bool planeBoxOverlap(const Eigen::Vector3d &normal, double d, const Eigen::Vector3d &hs) {
+bool planeBoxOverlap(const Eigen::Vector3d &normal, double d, const Eigen::Vector3d &hs) {
     Eigen::Vector3d vmin = {(normal[0] > 0.0) ? -hs[0] : +hs[0],
                             (normal[1] > 0.0) ? -hs[1] : +hs[1],
                             (normal[2] > 0.0) ? -hs[2] : +hs[2]};
@@ -40,7 +40,7 @@ namespace AABBTools {
     Eigen::Vector3d vmax = -vmin;
     if (normal.dot(vmax) + d >= 0.0) return true;
     return false;
-    }
+}
 }
 
 /*********************************************************/
@@ -55,17 +55,14 @@ Octree::AABB::AABB(double minx, double maxx, double miny, double maxy, double mi
     _min(minx, miny, minz), _max(maxx, maxy, maxz) {
 }
 
-
 Octree::AABB::AABB(const std::array<double, 6> &aabb) :
     AABB(aabb[0], aabb[1], aabb[2], aabb[3], aabb[4], aabb[5]) {
 }
 
-
 Eigen::Vector3d Octree::AABB::size() const {
     return _max - _min;
 }
 
-
 Eigen::Vector3d Octree::AABB::center() const {
     return (_max + _min) / 2.0;
 }
@@ -141,7 +138,6 @@ bool Octree::AABB::intersect_triangle(const Eigen::Vector3d &a,
     return true; // box and triangle overlaps
 }
 
-
 bool Octree::AABB::contains(const Eigen::Vector3d &v) const {
     return (v.x() >= _min.x() && v.y() >= _min.y() && v.z() >= _min.z()) &&
         (v.x() <= _max.x() && v.y() <= _max.y() && v.z() <= _max.z());
@@ -170,3 +166,79 @@ void Octree::AABB::extend(double step) {
     this->_min -= Eigen::Vector3d(step, step, step);
     this->_max += Eigen::Vector3d(step, step, step);
 }
+
+Octree::AABB Octree::AABB::ToCubeAABBWithSize(const Octree::AABB &aabb, const double res) {
+    int _;
+    return ToCubeAABBWithSize(aabb, res, _);
+}
+
+Octree::AABB Octree::AABB::ToCubeAABBWithSize(const Octree::AABB &aabb, const double res, int &num_on_axis) {
+    AABB cube = aabb;
+    Eigen::Vector3i num;
+
+    auto fit_to_integer_size = [](double &m, double &M, double res) -> int {
+        double S = M - m;
+        int num = std::ceil(S / res);
+
+        // padding num to pow of 2
+        if (num & (num - 1)) {
+            num |= num >> 1;
+            num |= num >> 2;
+            num |= num >> 4;
+            num |= num >> 8;
+            num |= num >> 16;
+            num = num + 1;
+        } else {
+            num = num == 0 ? 1 : num;
+        }
+
+        double error = (num * res - S) / 2.0;
+        m -= error;
+        M += error;
+        return num;
+    };
+    num.x() = fit_to_integer_size(cube._min.x(), cube._max.x(), res);
+    num.y() = fit_to_integer_size(cube._min.y(), cube._max.y(), res);
+    num.z() = fit_to_integer_size(cube._min.z(), cube._max.z(), res);
+    int max_id = std::max_element(num.data(), num.data() + 3) - num.data();
+
+    Eigen::Vector3d half_size = cube.size() / 2.0;
+    Eigen::Vector3d center = cube.center();
+
+    for (int i = 0; i < 3; ++i) {
+        if (max_id == i) continue;
+        cube._min.data()[i] = center.data()[i] - half_size.data()[max_id];
+        cube._max.data()[i] = center.data()[i] + half_size.data()[max_id];
+    }
+    num_on_axis = num.data()[max_id];
+    return cube;
+}
+
+bool Octree::AABB::intersect_aabb(const Octree::AABB &other) const {
+    return (_min.x() <= other._max.x() && _max.x() >= other._min.x()) &&
+        (_min.y() <= other._max.y() && _max.y() >= other._min.y()) &&
+        (_min.z() <= other._max.z() && _max.z() >= other._min.z());
+}
+
+
+Octree::AABB Octree::AABB::ToIntegerSizeAABB(const Octree::AABB &aabb,
+                                             const double res,
+                                             int &num_on_x,
+                                             int &num_on_y,
+                                             int &num_on_z) {
+    AABB cube = aabb;
+
+    auto fit_to_integer_size = [](double &m, double &M, double res) -> int {
+        double S = M - m;
+        int num = std::ceil(S / res);
+        double error = (num * res - S) / 2.0;
+        m -= error;
+        M += error;
+        return num;
+    };
+    num_on_x = fit_to_integer_size(cube._min.x(), cube._max.x(), res);
+    num_on_y = fit_to_integer_size(cube._min.y(), cube._max.y(), res);
+    num_on_z = fit_to_integer_size(cube._min.z(), cube._max.z(), res);
+
+    return cube;
+}

src/BaseOctree.cpp

@@ -10,13 +10,25 @@
 
 #include "Octree/BaseOctree.h"
 
-
 Octree::ocnode_sptr Octree::BaseOctree::map_node(const Eigen::Vector3d &pos) {
     auto fd = (root == nullptr) ? nullptr : root->map_node(pos);
     return fd;
     return fd == nullptr ? root : fd;
 }
 
+void Octree::OctreeNode::map_voxel(const AABB &aabb, set<ocnode_sptr> &nodes) {
+    if (aabb.intersect_aabb(this->aabb)) {
+        if (this->is_leaf()) {
+            nodes.insert(shared_from_this());
+        } else {
+            for (auto &child : this->children) {
+                if(child != nullptr)
+                    child->map_voxel(aabb, nodes);
+            }
+        }
+    }
+}
+
 Octree::BaseOctree::BaseOctree(uint32_t max_level, const Octree::AABB &aabb) :
     BaseOctree(max_level, std::make_shared<OctreeNode>(nullptr, aabb, 0)) {
 }
@@ -26,14 +38,13 @@ Octree::BaseOctree::BaseOctree(uint32_t max_level, Octree::ocnode_sptr root) : m
 }
 
 Octree::AABB Octree::BaseOctree::aabb() const {
-    if(root == nullptr){
+    if (root == nullptr) {
         return AABB();
-    }else{
+    } else {
         return root->aabb;
     }
 }
 
-
 Octree::OctreeNode::OctreeNode(ocnode_ptr father, const Octree::AABB &aabb, uint32_t level) : father(father),
                                                                                               level(level),
                                                                                               aabb(aabb) {
@@ -53,12 +64,10 @@ void Octree::OctreeNode::split_if_no_children() {
     }
 }
 
-
 Octree::ocnode_sptr Octree::OctreeNode::make_child(const Octree::AABB &aabb) {
-    return std::make_shared<OctreeNode>(this, aabb, level+1);
+    return std::make_shared<OctreeNode>(this, aabb, level + 1);
 }
 
-
 bool Octree::OctreeNode::is_leaf() const {
     return std::all_of(children.begin(), children.end(), [](ocnode_sptr child) { return child == nullptr; });
 }
@@ -72,7 +81,7 @@ Octree::ocnode_sptr Octree::OctreeNode::map_node(const Eigen::Vector3d &pos) {
     bool ct = aabb.contains(pos);
     if (ct) {
         if (is_leaf()) return shared_from_this();
-        for (const auto &child: this->children) {
+        for (const auto &child : this->children) {
             if (child == nullptr) continue;
             auto rt = child->map_node(pos);
             if (rt != nullptr) return rt;
@@ -81,3 +90,9 @@ Octree::ocnode_sptr Octree::OctreeNode::map_node(const Eigen::Vector3d &pos) {
     }
     return nullptr;
 }
+
+std::set<Octree::ocnode_sptr> Octree::BaseOctree::map_voxel(const Octree::AABB &aabb) {
+    std::set<ocnode_sptr> s;
+    this->root->map_voxel(aabb, s);
+    return std::move(s);
+}

src/voxel/SolidVoxelBuilder.cpp

@@ -0,0 +1,10 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#include "Octree/voxel/SolidVoxelBuilder.h"

src/voxel/SurfaceVoxelBuilder.cpp

@@ -0,0 +1,60 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#include "Octree/voxel/SurfaceVoxelBuilder.h"
+#include "Octree/voxel/Voxel.h"
+#include <queue>
+
+Octree::SurfaceVoxelBuilder::SurfaceVoxelBuilder(BaseOctree &octree,
+                                                 const pMesh::Triangle3dMesh<> &mesh,
+                                                 double block_size) :
+    block_size(block_size), octree(octree), mesh(mesh) {
+    aabb = AABB::ToIntegerSizeAABB(octree.aabb(), block_size, _x, _y, _z);
+}
+
+void Octree::SurfaceVoxelBuilder::build(Octree::Voxel &voxel) {
+    voxel._x = _x;
+    voxel._y = _y;
+    voxel._z = _z;
+    voxel.aabb = aabb;
+    voxel.voxels.clear();
+    voxel.voxels = std::vector<bool>(_x * _y * _z, false);
+
+    BOOST_LOG_TRIVIAL(debug) << "Voxel size: " << _x << " " << _y << " " << _z;
+
+    for (int i = 0; i < _x; ++i) {
+        double base_x = this->aabb.min().x() + i * block_size;
+        for (int j = 0; j < _y; ++j) {
+            double base_y = this->aabb.min().y() + j * block_size;
+            for (int k = 0; k < _z; ++k) {
+                double base_z = this->aabb.min().z() + k * block_size;
+
+                AABB vox({base_x, base_y, base_z}, {base_x + block_size, base_y + block_size, base_z + block_size});
+                auto candidates_vox = octree.map_voxel(vox);
+                bool exit_flag = false;
+                for (auto candidate : candidates_vox) {
+                    for (uint32_t tid : candidate->tri_ids) {
+                        const pMesh::Triangle3dMesh<>::Face &face = mesh.faces[tid].attr;
+                        if (vox.intersect_triangle(mesh.vertices[face.vertices[0].id()].attr.coordinate,
+                                                   mesh.vertices[face.vertices[1].id()].attr.coordinate,
+                                                   mesh.vertices[face.vertices[2].id()].attr.coordinate)) {
+                            voxel.voxels[voxel.id_at(i, j, k)] = true;
+                            exit_flag = true;
+                            break;
+                        }
+                    }
+                    if(exit_flag) {
+                        break;
+                    }
+                }
+            }
+        }
+    }
+}
+

src/voxel/Voxel.cpp

@@ -0,0 +1,24 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2022/2/23
+ * @email: yjxkwp\@foxmail.com
+ * @description: 
+ * ------------------------------------
+**/
+
+#include "Octree/voxel/Voxel.h"
+
+Octree::Voxel::Voxel(const std::vector<bool> &value, int _x, int _y, int _z): voxels(value), _x(_x), _y(_y), _z(_z) {
+
+}
+
+Octree::Voxel::Voxel(int _x, int _y, int _z): Voxel(std::vector<bool>(_x * _y * _z, false), _x, _y, _z) {
+
+}
+int Octree::Voxel::id_at(int x, int y, int z) const {
+    return ((z * _y) + y) * _x + x;
+}
+bool Octree::Voxel::at(int x, int y, int z) const {
+    return voxels[id_at(x,y,z)];
+}

tests/solid_octree_test/.gitignore

@@ -0,0 +1 @@
+test-path.h

tests/solid_octree_test/CMakeLists.txt

@@ -0,0 +1,6 @@
+configure_file(
+        "${CMAKE_CURRENT_SOURCE_DIR}/test-path.h.in"
+        "${CMAKE_CURRENT_SOURCE_DIR}/test-path.h"
+)
+add_executable(solid_octree_test main.cpp)
+target_link_libraries(solid_octree_test Octree pMesh)

tests/solid_octree_test/main.cpp

@@ -0,0 +1,178 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2021/11/17
+ * @email: yjxkwp@foxmail.com
+ * @site: https://donot.fit
+ * @description: 
+ * ------------------------------------
+**/
+
+#include <iostream>
+#include "Octree/sdf/SDFOctree.h"
+#include "Octree/OctreeBuilder.h"
+#include "Octree/OctreeTraverser.h"
+#include "Octree/sdf/SDFTraversalSampler.h"
+#include <igl/marching_cubes.h>
+#include <igl/voxel_grid.h>
+#include <igl/writeOBJ.h>
+#include <pMesh/mesh/TriangleMesh.h>
+#include <pMesh/mesh/HexahedronMesh.h>
+#include <pMesh/io/reader/OBJReader.h>
+#include <pMesh/io/reader/VTKReader.h>
+#include <pMesh/io/writer/VTKWriter.h>
+#include <pMesh/io/adapter/DefaultReadAdapter.h>
+#include <pMesh/io/adapter/DefaultWriteAdapter.h>
+#include <boost/timer.hpp>
+#include <boost/log/trivial.hpp>
+#include <pMesh/io/reader/BaseReader.h>
+#include <pMesh/io/adapter/DefaultReadAdapter.h>
+#include "test-path.h"
+
+#include "Octree/voxel/SurfaceVoxelBuilder.h"
+#include "Octree/voxel/Voxel.h"
+
+class VTKTraverser : public pMesh::io::BaseReader, public Octree::OctreeTraverser {
+    int node_cnt = 0;
+public:
+    explicit VTKTraverser(const Octree::BaseOctree &octree, const TraverseStrategy strategy = DFS);
+
+    void visit_node(Octree::OctreeNode &node) override;
+
+    pMesh::io::ReadAdapter *reader_adapter;
+
+    bool operator>>(pMesh::io::ReadAdapter &adapter) override;
+};
+
+
+int main() {
+    pMesh::io::fs_path data_base_path = TEST_DATA_BASE_PATH;
+    BOOST_LOG_TRIVIAL(debug) << "base data path is " << boost::filesystem::absolute(data_base_path);
+
+    auto mesh_path = data_base_path / "stanford-bunny.obj";
+    auto out_octree_path = data_base_path / "stanford-bunny-octree.vtk";
+    auto out_voxel_path = data_base_path / "stanford-bunny-voxel.vtk";
+    auto out_error_path = data_base_path / "stanford-bunny-error_points.vtk";
+
+    pMesh::Triangle3dMesh<> mesh;
+    pMesh::io::OBJReader(mesh_path)
+        >> pMesh::io::DefaultSurfaceReadAdapter<>(mesh, false)();
+
+    BOOST_LOG_TRIVIAL(info) << "Load Face #" << mesh.f_size();
+
+    const int level = 4;
+    const double block_size = 0.005;
+
+    boost::timer t;
+    Octree::AABB aabb(mesh.aabb());
+    Octree::BaseOctree octree(level, aabb);
+
+    Octree::OctreeBuilder builder(mesh, octree);
+    builder.build();
+    BOOST_LOG_TRIVIAL(debug) << "time elapse " << t.elapsed();
+
+    auto aabb_size = octree.aabb().size();
+    BOOST_LOG_TRIVIAL(debug) << "Global AABB " << aabb_size.transpose();
+
+    auto final_aabb_size = aabb_size / pow(2, level);
+    BOOST_LOG_TRIVIAL(debug) << "Local AABB " << final_aabb_size.transpose();
+
+    pMesh::HexahedronMesh<> hexahedron_mesh;
+    VTKTraverser(octree) >> pMesh::io::DefaultVolumeReadAdapter<>(hexahedron_mesh)();
+
+    pMesh::io::VTKWriter(12, out_octree_path)
+        << pMesh::io::DefaultVolumeWriteAdapter<>(hexahedron_mesh)();
+
+    t.restart();
+    Octree::Voxel voxel;
+    Octree::SurfaceVoxelBuilder voxel_builder(octree, mesh, block_size);
+    voxel_builder.build(voxel);
+    BOOST_LOG_TRIVIAL(debug) << "SurfaceVoxelBuilder time elapse " << t.elapsed();
+
+    // output voxel
+    pMesh::HexahedronMesh<> voxelMesh;
+    {
+        auto aabb = voxel.aabb;
+        int _x = voxel._x;
+        int _y = voxel._y;
+        int _z = voxel._z;
+
+        for (int i = 0; i < _x; ++i) {
+            double base_x = aabb.min().x() + i * block_size;
+            for (int j = 0; j < _y; ++j) {
+                double base_y = aabb.min().y() + j * block_size;
+                for (int k = 0; k < _z; ++k) {
+                    double base_z = aabb.min().z() + k * block_size;
+
+                    if(voxel.voxels[voxel.id_at(i, j, k)]){
+                        Octree::AABB v({base_x, base_y, base_z}, {base_x + block_size, base_y + block_size, base_z + block_size});
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(0, {v.min().x(),v.min().y(),v.min().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(1, {v.max().x(),v.min().y(),v.min().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(2, {v.max().x(),v.max().y(),v.min().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(3, {v.min().x(),v.max().y(),v.min().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(4, {v.min().x(),v.min().y(),v.max().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(5, {v.max().x(),v.min().y(),v.max().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(6, {v.max().x(),v.max().y(),v.max().z()}));
+                        voxelMesh.vertices.emplace_back(pMesh::Volume::Vertex(7, {v.min().x(),v.max().y(),v.max().z()}));
+
+                        int base_id = voxelMesh.v_size();
+                        pMesh::Volume::Cell cell;
+                        using PPP = pMesh::Volume::VertexHandle;
+                        cell.vertices = {
+                            PPP(base_id + 0), PPP(base_id + 1), PPP(base_id + 2), PPP(base_id + 3),
+                            PPP(base_id + 4), PPP(base_id + 5), PPP(base_id + 6), PPP(base_id + 7)
+                        };
+                        voxelMesh.cells.emplace_back(cell);
+                    }
+                }
+            }
+        }
+    }
+
+    pMesh::io::VTKWriter(12, out_voxel_path)
+        << pMesh::io::DefaultVolumeWriteAdapter<>(voxelMesh)();
+
+    return 0;
+}
+
+void VTKTraverser::visit_node(Octree::OctreeNode &node) {
+    if (node.get_level() != octree.get_max_level()) {
+//        assert(!node.is_leaf());
+        assert(node.tri_ids.empty());
+        return;
+    }
+
+    auto min = node.get_aabb().min();
+    auto max = node.get_aabb().max();
+    reader_adapter->feed_vertex({min.x(), min.y(), min.z()});
+    reader_adapter->feed_vertex({max.x(), min.y(), min.z()});
+    reader_adapter->feed_vertex({max.x(), max.y(), min.z()});
+    reader_adapter->feed_vertex({min.x(), max.y(), min.z()});
+    reader_adapter->feed_vertex({min.x(), min.y(), max.z()});
+    reader_adapter->feed_vertex({max.x(), min.y(), max.z()});
+    reader_adapter->feed_vertex({max.x(), max.y(), max.z()});
+    reader_adapter->feed_vertex({min.x(), max.y(), max.z()});
+
+    reader_adapter->feed_collection(
+        {
+            0 + node_cnt, 1 + node_cnt, 2 + node_cnt, 3 + node_cnt, 4 + node_cnt,
+            5 + node_cnt, 6 + node_cnt, 7 + node_cnt
+        });
+    node_cnt += 8;
+}
+
+bool VTKTraverser::operator>>(pMesh::io::ReadAdapter &adapter) {
+    reader_adapter = &adapter;
+    node_cnt = 0;
+    adapter.start();
+
+    this->traverse();
+
+    adapter.end();
+    return true;
+}
+
+VTKTraverser::VTKTraverser(const Octree::BaseOctree &octree, const OctreeTraverser::TraverseStrategy strategy)
+    : OctreeTraverser(
+    octree, strategy) {
+}

tests/solid_octree_test/test-path.h.in

@@ -0,0 +1,16 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2021/12/1
+ * @email: yjxkwp@foxmail.com
+ * @site: https://donot.fit
+ * @description: 
+ * ------------------------------------
+**/
+
+#ifndef OCTREE_SDF_TEST_PATH_H
+#define OCTREE_SDF_TEST_PATH_H
+
+#define TEST_DATA_BASE_PATH "@TEST_DATA_BASE_PATH@"
+
+#endif //OCTREE_SDF_TEST_PATH_H