construct bvh

CMakeLists.txt

@@ -32,4 +32,4 @@ source_group("3rdparty" FILES ${THRDPARTY_FILES})
 add_executable(${PROJECT_NAME} ${SRC_FILES} ${HEADER_FILES} ${THRDPARTY_FILES} )
 
 
-target_link_libraries(${PROJECT_NAME} PUBLIC igl::glfw)
+target_link_libraries(${PROJECT_NAME} PUBLIC igl::glfw cln ginac)

include/bvh.hpp

@@ -45,8 +45,8 @@ public:
     };
     BVH_AABB_NodePtr Build_NurbsCurve(tinynurbs::RationalCurve<double> &curve, BVH_AABB_NodePtr &curNode, int curLevel, double t, double tstep);
     BVH_AABB_NodePtr Build_NurbsSurface(tinynurbs::RationalSurface<double> &surface, BVH_AABB_NodePtr &curNode, int curLevel, double u, double v, double ustep, double vstep);
-    AABB &UnionBound(const AABB &a, const AABB &b);
-    AABB &UpdateBound(const AABB &bound, const glm::vec3 &p);
+    AABB UnionBound(const AABB &a, const AABB &b);
+    AABB UpdateBound(const AABB &bound, const glm::vec3 &p);
 };
 
 #endif // !BVH_H_

include/show_in_libigl.hpp

@@ -1,5 +0,0 @@
-#include "bvh.hpp"
-#include <igl/opengl/glfw/Viewer.h>
-
-void ShowCurve_Igl(igl::opengl::glfw::Viewer viewer, tinynurbs::RationalCurve<double> &curve, double sampleNum);
-void ShowSurface_Igl(igl::opengl::glfw::Viewer viewer, tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV);

include/show_libigl.hpp

@@ -0,0 +1,12 @@
+#ifndef SHOW_LIBIGL_H_
+#define SHOW_LIBIGL_H_
+
+#include "bvh.hpp"
+#include <igl/opengl/glfw/Viewer.h>
+
+void ShowCurve_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalCurve<double> &curve, double sampleNum);
+void ShowSurface_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV);
+void ShowBVHNode(igl::opengl::glfw::Viewer &viewer, BVH_AABB_NodePtr bvhNode);
+void ShowAABB(igl::opengl::glfw::Viewer &viewer, AABB bound);
+
+#endif

src/bvh.cpp

@@ -1,6 +1,6 @@
 #include "bvh.hpp"
 
-AABB &BVH_AABB::UnionBound(const AABB &a, const AABB &b)
+AABB BVH_AABB::UnionBound(const AABB &a, const AABB &b)
 {
     AABB ret;
     ret.Bmin = glm::vec3(fmin(a.Bmin.x, b.Bmin.x),
@@ -13,7 +13,7 @@ AABB &BVH_AABB::UnionBound(const AABB &a, const AABB &b)
     return ret;
 }
 
-AABB &BVH_AABB::UpdateBound(const AABB &bound, const glm::vec3 &p)
+AABB BVH_AABB::UpdateBound(const AABB &bound, const glm::vec3 &p)
 {
     AABB ret;
     ret.Bmin = glm::vec3(fmin(bound.Bmin.x, p.x),
@@ -28,7 +28,7 @@ AABB &BVH_AABB::UpdateBound(const AABB &bound, const glm::vec3 &p)
 
 BVH_AABB_NodePtr BVH_AABB::Build_NurbsCurve(tinynurbs::RationalCurve<double> &curve, BVH_AABB_NodePtr &curNode, int curLevel, double t, double tstep)
 {
-    if (curLevel == Hlevel)
+    if (curLevel + 1 == Hlevel)
     {
         for (double i = t; i < t + tstep; i += tstep / sampleNum)
         {

src/main.cpp

@@ -1,34 +1,46 @@
 #include <igl/opengl/glfw/Viewer.h>
+#include <tinynurbs/tinynurbs.h>
+#include "bvh.hpp"
+#include "show_libigl.hpp"
 
 int main(int argc, char *argv[])
 {
-    // Inline mesh of a cube
-    const Eigen::MatrixXd V = (Eigen::MatrixXd(8, 3) << 0.0, 0.0, 0.0,
-                               0.0, 0.0, 1.0,
-                               0.0, 1.0, 0.0,
-                               0.0, 1.0, 1.0,
-                               1.0, 0.0, 0.0,
-                               1.0, 0.0, 1.0,
-                               1.0, 1.0, 0.0,
-                               1.0, 1.0, 1.0)
-                                  .finished();
-    const Eigen::MatrixXi F = (Eigen::MatrixXi(12, 3) << 0, 6, 4,
-                               0, 2, 6,
-                               0, 3, 2,
-                               0, 1, 3,
-                               2, 7, 6,
-                               2, 3, 7,
-                               4, 6, 7,
-                               4, 7, 5,
-                               0, 4, 5,
-                               0, 5, 1,
-                               1, 5, 7,
-                               1, 7, 3)
-                                  .finished();
-
-    // Plot the mesh
     igl::opengl::glfw::Viewer viewer;
-    viewer.data().set_mesh(V, F);
-    viewer.data().set_face_based(true);
+
+    viewer.data().point_size = 5;
+
+    /*
+    tinynurbs::RationalCurve<double> crv;      // Planar curve using float32
+    crv.control_points = {glm::vec3(-3, 0, 0), // std::vector of 3D points
+                          glm::vec3(-2, 2, 3),
+                          glm::vec3(-1, 5, 0),
+                          glm::vec3(-3, -1, -1),
+                          glm::vec3(7, 6, 4)};
+    crv.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+    crv.degree = 2;
+    crv.weights = {1, 1, 1, 1, 1};
+
+    ShowCurve_Igl(viewer, crv, 100);
+    BVH_AABB bvh_curve(5, 2, 100);
+    double tstep = *(crv.knots.end() - 1);
+    bvh_curve.Build_NurbsCurve(crv, bvh_curve.bvh_aabb_node, 0, 0, tstep);
+    ShowBVHNode(viewer, bvh_curve.bvh_aabb_node);
+    */
+    tinynurbs::RationalCurve<double> crv;        // Planar curve using float32
+    crv.control_points = {glm::vec3(-3, -3, -3), // std::vector of 3D points
+                          glm::vec3(-2, -2, -2),
+                          glm::vec3(-1, -1, -1),
+                          glm::vec3(0, 0, 0),
+                          glm::vec3(1, 1, 1)};
+    crv.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+    crv.degree = 2;
+    crv.weights = {1, 1, 1, 1, 1};
+
+    ShowCurve_Igl(viewer, crv, 100);
+    BVH_AABB bvh_curve(5, 2, 100);
+    double tstep = *(crv.knots.end() - 1);
+    bvh_curve.Build_NurbsCurve(crv, bvh_curve.bvh_aabb_node, 0, 0, tstep);
+    ShowBVHNode(viewer, bvh_curve.bvh_aabb_node);
+
     viewer.launch();
 }

src/newton.cpp

@@ -1,6 +1,8 @@
 #include "newton.hpp"
 
 // 默认id参数不用传，表示四个变量四个方程，id=1 or 2 的时候是2x2的方程组， id=3时是6x6的方程组, id=4 (x0,x1,x2... xn-1) 表示的方程组
+
+/*
 vector<numeric> newton(const MUL_F &F, const vector<numeric> &init_point, const numeric &eps, bool &isConvergence, int id)
 {
     vector<symbol> symbols = {get_symbol("u"), get_symbol("v"), get_symbol("s"), get_symbol("t")};
@@ -150,3 +152,5 @@ vector<numeric> newton(const MUL_F &F, const vector<numeric> &init_point, const
     //    }
     return x_val;
 }
+
+*/

src/show_in_libigl.cpp

@@ -1,25 +0,0 @@
-#include "show_in_libigl.hpp"
-
-void ShowCurve_Igl(igl::opengl::glfw::Viewer viewer, tinynurbs::RationalCurve<double> &curve, double sampleNum)
-{
-    double T = *(curve.knots.end() - 1);
-
-    for (double i = 0; i < T; i += T / sampleNum)
-    {
-        auto point = tinynurbs::curvePoint(curve, i);
-        pt << point.x, point.y, point.z;
-        viewer.data().add_points(pt, Eigen::RowVector3d(0, 1, 0));
-    }
-}
-
-void ShowSurface_Igl(igl::opengl::glfw::Viewer viewer, tinynurbs::RationalCurve<double> &curve, double sampleNumU, double sampleNumV)
-{
-    double T = *(curve.knots.end() - 1);
-
-    for (double i = 0; i < T; i += T / sampleNum)
-    {
-        auto point = tinynurbs::curvePoint(curve, i);
-        pt << point.x, point.y, point.z;
-        viewer.data().add_points(pt, Eigen::RowVector3d(0, 1, 0));
-    }
-}

src/show_libigl.cpp

@@ -0,0 +1,85 @@
+#include "bvh.hpp"
+#include "show_libigl.hpp"
+
+void ShowCurve_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalCurve<double> &curve, double sampleNum)
+{
+    double T = *(curve.knots.end() - 1);
+
+    for (double i = 0; i < T; i += T / sampleNum)
+    {
+        auto point = tinynurbs::curvePoint(curve, i);
+        Eigen::RowVector3d pt;
+        pt << point.x, point.y, point.z;
+        viewer.data().add_points(pt, Eigen::RowVector3d(1, 1, 0));
+    }
+}
+
+void ShowSurface_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV)
+{
+    double U = *(surface.knots_u.end() - 1);
+    double V = *(surface.knots_v.end() - 1);
+
+    for (double u = 0; u <= U; u += U / sampleNumU)
+    {
+        for (double v = 0; v <= V; v += V / sampleNumV)
+        {
+            auto point = tinynurbs::surfacePoint(surface, u, v);
+            Eigen::RowVector3d pt;
+            pt << point.x, point.y, point.z;
+            viewer.data().add_points(pt, Eigen::RowVector3d(1, 0.7, 0.5));
+        }
+    }
+}
+
+void ShowAABB(igl::opengl::glfw::Viewer &viewer, AABB bound)
+{
+
+    Eigen::Vector3d m;
+    m << bound.Bmin.x, bound.Bmin.y, bound.Bmin.z;
+    Eigen::Vector3d M;
+    M << bound.Bmax.x, bound.Bmax.y, bound.Bmax.z;
+
+    Eigen::MatrixXd V_box(8, 3);
+    V_box << m(0), m(1), m(2),
+        M(0), m(1), m(2),
+        M(0), M(1), m(2),
+        m(0), M(1), m(2),
+        m(0), m(1), M(2),
+        M(0), m(1), M(2),
+        M(0), M(1), M(2),
+        m(0), M(1), M(2);
+
+    Eigen::MatrixXi E_box(12, 2);
+    E_box << 0, 1,
+        1, 2,
+        2, 3,
+        3, 0,
+        4, 5,
+        5, 6,
+        6, 7,
+        7, 4,
+        0, 4,
+        1, 5,
+        2, 6,
+        7, 3;
+
+    for (unsigned i = 0; i < E_box.rows(); ++i)
+        viewer.data().add_edges(
+            V_box.row(E_box(i, 0)),
+            V_box.row(E_box(i, 1)),
+            Eigen::RowVector3d(1, 1, 1));
+    return;
+}
+
+void ShowBVHNode(igl::opengl::glfw::Viewer &viewer, BVH_AABB_NodePtr bvhNode)
+{
+    if (bvhNode == NULL)
+        return;
+
+    ShowAABB(viewer, bvhNode->bound);
+
+    for (auto it : bvhNode->childPtr)
+    {
+        ShowBVHNode(viewer, it);
+    }
+}