测试udf翻转

CMakeLists.txt

@@ -37,6 +37,7 @@ if(Octree_BUILD_TEST)
     add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/tests/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)
 endif()
 
 enable_testing()
@@ -44,4 +45,5 @@ if(Octree_BUILD_TEST)
     add_test(NAME octree_test COMMAND $<TARGET_FILE: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>)
 endif()

include/Octree/VoxelDenseData.h

@@ -14,13 +14,14 @@
 #include <vector>
 #include <Eigen/Core>
 #include <Eigen/Dense>
+#include <tuple>
 
 namespace Octree {
     /**
      * All Dense value will be sampled use coordinate [0, 1]^3
      */
     class VoxelDenseData {
-    protected:
+    public:
         int _x, _y, _z;  // 3D range [start with 0]
         std::vector<double> value;
 
@@ -45,6 +46,16 @@ namespace Octree {
          * @return
          */
         int get_array_index(int x, int y, int z) const;
+
+    public:
+        // TODO: 测试用，会删除
+        void neg_at_pos(int x, int y, int z) {
+            value[get_array_index(x, y, z)] = -abs(value[get_array_index(x, y, z)]);
+        }
+
+        void neg_at_index(uint32_t ind){
+            value[ind] = -abs(value[ind]);
+        }
     };
 }
 

tests/search_assign_test/.gitignore

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

tests/search_assign_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(search_assign_test main.cpp search_assign.cpp search_assign.h)
+target_link_libraries(search_assign_test Octree pMesh)

tests/search_assign_test/main.cpp

@@ -0,0 +1,114 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2021/11/17
+ * @email: yjxkwp@foxmail.com
+ * @site: https://donot.fit
+ * @description: 
+ * ------------------------------------
+**/
+
+#include <iostream>
+#include <igl/marching_cubes.h>
+#include <igl/voxel_grid.h>
+#include <igl/writeOBJ.h>
+#include <igl/in_element.h>
+#include <pMesh/mesh/TriangleMesh.h>
+#include <pMesh/mesh/HexahedronMesh.h>
+#include <pMesh/io/reader/OBJReader.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 "Octree/udf/UDFOctree.h"
+#include "Octree/OctreeBuilder.h"
+#include "Octree/OctreeTraverser.h"
+#include "Octree/udf/UDFTraversalSampler.h"
+#include "test-path.h"
+#include "search_assign.h"
+
+
+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);
+    bool use_udf = true;
+    auto out_mc_path = data_base_path / (use_udf? "sphere-udf-mc.obj": "sphere-sdf-mc.obj");
+
+    Eigen::AlignedBox<double, 3> aabb(Eigen::Vector3d(0,0,0), Eigen::Vector3d(1, 1, 1));
+
+    int aabb_res_x = 20;
+    int aabb_res_y = 20;
+    int aabb_res_z = 20;
+    int aabb_len = aabb_res_x * aabb_res_y * aabb_res_z;
+
+    double x0 = 0;
+    double y0 = 0;
+    double z0 = 0;
+    double R = 1.2;
+
+    auto sphere_shape = [=](double x, double y, double z) {
+        return sqrt(pow(x - x0, 2) + pow(y - y0, 2) + pow(z - z0, 2)) - R;
+    };
+    Octree::UDFOctree octree(0, Octree::AABB(aabb.min(), aabb.max()));
+    octree.get_root()->udf = std::make_shared<Octree::UDFData>(std::vector<double>(aabb_len), aabb_res_x, aabb_res_y,
+                                                               aabb_res_z);
+    auto &udf_data = *octree.get_root()->udf;
+    BOOST_LOG_TRIVIAL(debug) << "Calculating UDF";
+
+    for (int x = 0; x < aabb_res_x; ++x) {
+        double pos_x = ((aabb_res_x - x) * aabb.min().x() + x * aabb.max().x()) / (1.0 * aabb_res_x);
+        for (int y = 0; y < aabb_res_y; ++y) {
+            double pos_y = ((aabb_res_y - y) * aabb.min().y() + y * aabb.max().y()) / (1.0 * aabb_res_y);
+            for (int z = 0; z < aabb_res_z; ++z) {
+                double pos_z = ((aabb_res_z - z) * aabb.min().z() + z * aabb.max().z()) / (1.0 * aabb_res_z);
+
+                int ind = udf_data.get_array_index(x, y, z);
+                if(use_udf)
+                    udf_data.value[ind] = abs(sphere_shape(pos_x, pos_y, pos_z));
+                else
+                    udf_data.value[ind] = sphere_shape(pos_x, pos_y, pos_z);
+            }
+        }
+    }
+    BOOST_LOG_TRIVIAL(debug) << "End computing UDF";
+
+    if(use_udf) {
+        int start[3] = {aabb_res_x/2, aabb_res_y/2, aabb_res_z/2};
+        double block_size[3];
+        // TODO: 
+        double threshold = 0;
+        search_assign(udf_data, start, threshold);
+    }
+
+    if (1) {
+        Eigen::MatrixXd GV;
+        Eigen::RowVector3i res;
+        const int s = 200;
+        igl::voxel_grid(aabb, s, 1, GV, res);
+
+        // compute values
+        std::cout << "Computing distances..." << std::endl;
+        Eigen::VectorXd S = Eigen::VectorXd(GV.rows()), B;
+
+        {
+            for (int i = 0; i < GV.rows(); ++i) {
+                auto node = octree.map_node(GV.row(i));
+                if (node == nullptr) {
+                    S[i] = 1000;
+                } else {
+                    S[i] = node->get_udf(GV.row(i));
+                }
+            }
+        }
+        std::cout << "Marching cubes..." << std::endl;
+        Eigen::MatrixXd SV, BV;
+        Eigen::MatrixXi SF, BF;
+
+        igl::marching_cubes(S, GV, res(0), res(1), res(2), 0, SV, SF);
+        igl::writeOBJ(out_mc_path.string(), SV, SF);
+    }
+    return 0;
+}

tests/search_assign_test/search_assign.cpp

@@ -0,0 +1,79 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2021/12/6
+ * @email: yjxkwp@foxmail.com
+ * @site: https://donot.fit
+ * @description: 
+ * ------------------------------------
+**/
+
+#include "search_assign.h"
+#include <pMesh/mesh/TriangleMesh.h>
+#include <iostream>
+
+bool *visited;
+
+void
+dfs_for_assigning(Octree::UDFData &udf, int current_x, int current_y, int current_z, int ind, const double threshold,
+                  std::function<void(const double *)> recorder) {
+    if (visited[ind] || ind == -1) return;
+    visited[ind] = true;
+
+    if (abs(udf.value[ind]) < threshold) {
+        return;
+//        double pos[] = {current_x * 1.0, current_y * 1.0, current_z * 1.0, udf.value[ind]};
+//        recorder(pos);
+    }
+
+    udf.value[ind] = -abs(udf.value[ind]);
+
+    // sub routine
+    for (int i = -1; i <= 1; ++i) {
+        int x = current_x + i;
+        if (x < 0 || x >= udf._x) continue;
+
+        for (int j = -1; j <= 1; ++j) {
+            int y = current_y + j;
+            if (y < 0 || y >= udf._y) continue;
+
+            for (int k = -1; k <= 1; ++k) {
+                int z = current_z + k;
+                if (z < 0 || z >= udf._z) continue;
+
+                int next_ind = udf.get_array_index(x, y, z);
+                dfs_for_assigning(udf, x, y, z, next_ind, threshold, recorder);
+            }
+        }
+    }
+
+    // sub routine
+//    for (auto nb: neighbors) {
+//        if (nb[0] < 0 || nb[0] >= udf._x ||
+//            nb[1] < 0 || nb[1] >= udf._y ||
+//            nb[2] < 0 || nb[2] >= udf._z) {
+//            continue;
+//        } else {
+//            int next_ind = udf.get_array_index(nb[0], nb[1], nb[2]);
+//            dfs_for_assigning(udf, nb[0], nb[1], nb[2], next_ind, threshold, recorder);
+//        }
+//    }
+}
+
+void search_assign(Octree::UDFData &udf, int *start, double threshold) {
+    visited = new bool[udf._x * udf._y * udf._z];
+
+    pMesh::Triangle3dMesh point_cloud;
+    auto recoder = [&](const double *pos) {
+//        point_cloud.vertices.emplace_back(pMesh::Triangle3dMesh::VertexField())
+//        std::cout << pos[0] << " " << pos[1] << " " << pos[2] << std::endl;
+    };
+
+    // TODO: 选取开始点
+    int next_ind = udf.get_array_index(start[0], start[1], start[2]);
+    dfs_for_assigning(udf, start[0], start[1], start[2], next_ind, threshold, recoder);
+    // TODO: 爆栈问题
+
+    delete visited;
+    visited = nullptr;
+};

tests/search_assign_test/search_assign.h

@@ -0,0 +1,34 @@
+/**
+ * ------------------------------------
+ * @author: Weipeng Kong
+ * @date: 2021/12/6
+ * @email: yjxkwp@foxmail.com
+ * @site: https://donot.fit
+ * @description: 测试DFS赋值
+ * ------------------------------------
+**/
+
+#ifndef OCTREE_SEARCH_ASSIGN_H
+#define OCTREE_SEARCH_ASSIGN_H
+
+#include "Octree/udf/UDFData.h"
+
+
+/**
+ *   _________________________
+ *   |       |       |       |
+ *   |       |       |       |
+ *   |_______|_______|_______|
+ *   |       |   ↑   |       |
+ *   |       | ← d → |       |
+ *   |_______|___↓___|_______|
+ *   |       |       |       |
+ *   |       |       |       |
+ *   |_______|_______|_______|
+ */
+
+void search_assign(Octree::UDFData &udf, int *start, double threshold);
+
+
+
+#endif //OCTREE_SEARCH_ASSIGN_H

tests/search_assign_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

tests/udf_test/main.cpp

@@ -12,6 +12,7 @@
 #include <igl/marching_cubes.h>
 #include <igl/voxel_grid.h>
 #include <igl/writeOBJ.h>
+#include <igl/in_element.h>
 #include <pMesh/mesh/TriangleMesh.h>
 #include <pMesh/mesh/HexahedronMesh.h>
 #include <pMesh/io/reader/OBJReader.h>