add some files

include/bvh.hpp

@@ -4,6 +4,8 @@
 #include <iostream>
 #include <tinynurbs/tinynurbs.h>
 
+class BVH_AABB;
+
 //包围盒结构，存储盒子对角的两个点
 typedef struct AABB
 {
@@ -15,6 +17,12 @@ typedef struct AABB
         Bmax = glm::vec3(minNum, minNum, minNum);
         Bmin = glm::vec3(maxNum, maxNum, maxNum);
     }
+    void ShowAABB()
+    {
+
+        std::cout << "min_point (" << Bmin.x << ", " << Bmin.y << ", " << Bmin.z << ") , max_point (" << Bmax.x << ", "
+                  << Bmax.y << ", " << Bmax.z << ")\n";
+    }
 } AABB;
 
 // AABB包围盒的结点

include/intersection.hpp

@@ -0,0 +1,7 @@
+#ifndef INTERSECTION_H_
+#define INTERSECTION_H_
+
+#include "bvh.hpp"
+bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2, std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>> &IstNodePtr);
+
+#endif

include/show_libigl.hpp

@@ -4,9 +4,19 @@
 #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);
+extern igl::opengl::glfw::Viewer viewer;
+
+#define RED Eigen ::RowVector3d(1, 0, 0)
+#define GREEN Eigen ::RowVector3d(0, 1, 0)
+#define BLUE Eigen ::RowVector3d(0, 0, 1)
+#define WHITE Eigen ::RowVector3d(1, 1, 1)
+#define BLACK Eigen ::RowVector3d(0, 0, 0)
+#define YELLOW Eigen ::RowVector3d(1, 1, 0)
+#define PINK Eigen ::RowVector3d(0.8, 0.2, 0.6)
+
+void ShowCurve_Igl(tinynurbs::RationalCurve<double> &curve, double sampleNum, Eigen::RowVector3d color);
+void ShowSurface_Igl(tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV, Eigen::RowVector3d color);
+void ShowBVHNode_Igl(BVH_AABB_NodePtr bvhNode, Eigen::RowVector3d color);
+void ShowAABB_Igl(AABB bound, Eigen::RowVector3d color);
 
 #endif

src/intersection.cpp

@@ -1,17 +1,33 @@
-#include "bvh.cpp"
-#include <iostream>
-#include <algorithm>
+#include "intersection.hpp"
+#include "show_libigl.hpp"
 
 bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2, std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>> &IstNodePtr)
 {
     AABB box1 = BoxPtr1->bound;
     AABB box2 = BoxPtr2->bound;
 
+    /*
+    std::cout << "Box1 : ";
+    box1.ShowAABB();
+    ShowAABB_Igl(box1, Eigen::RowVector3d(1, 0, 0));
+    std::cout << "Box2 : ";
+    box2.ShowAABB();
+    ShowAABB_Igl(box2, Eigen::RowVector3d(0, 0, 1));
+    */
+
     auto LineIntersect = [](double x1, double x2, double x1_, double x2_)
     {
-        bool LIntersect = false;
-        if ((x2 <= x2_ && x2 >= x1_) || (x1_ <= x2 && x1_ >= x1))
-            LIntersect = true;
+        bool LIntersect = true;
+        if (x1 <= x1_)
+        {
+            if (x1_ > x2)
+                LIntersect = false;
+        }
+        else
+        {
+            if (x1 > x2_)
+                LIntersect = false;
+        }
         return LIntersect;
     };
     auto BoxIntersect = [&]()
@@ -34,17 +50,19 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
 
     if (BoxIntersect())
     {
+        // std::cout << "Above boxes intersect! The bigger box is " << biggerbox << '\n';
         if (biggerbox == 1)
         {
             if (BoxPtr1->childPtr.size() == 0)
             {
                 std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr> pi = std::make_pair(BoxPtr1, BoxPtr2);
+                std::cout << "We find a intersection box pair\n";
                 IstNodePtr.push_back(pi);
                 return true;
             }
             for (auto it : BoxPtr1->childPtr)
             {
-                CurveCurveBVHIntersect(it, BoxPtr2, IstNodePtr)
+                CurveCurveBVHIntersect(it, BoxPtr2, IstNodePtr);
             }
         }
         else
@@ -52,6 +70,7 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
             if (BoxPtr2->childPtr.size() == 0)
             {
                 std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr> pi = std::make_pair(BoxPtr2, BoxPtr1);
+                std::cout << "We find a intersection box pair\n";
                 IstNodePtr.push_back(pi);
                 return true;
             }
@@ -60,6 +79,7 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
                 CurveCurveBVHIntersect(BoxPtr1, it, IstNodePtr);
             }
         }
+        return true;
     }
     else
     {

src/main.cpp

@@ -2,40 +2,68 @@
 #include <tinynurbs/tinynurbs.h>
 #include "bvh.hpp"
 #include "show_libigl.hpp"
+#include "intersection.hpp"
+
+igl::opengl::glfw::Viewer viewer;
 
 int main(int argc, char *argv[])
 {
-    igl::opengl::glfw::Viewer viewer;
 
     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};
+    tinynurbs::RationalCurve<double> crv1;      // Planar curve using float32
+    crv1.control_points = {glm::vec3(-3, 0, 0), // std::vector of 3D points
+                           glm::vec3(-1, 4, 0),
+                           glm::vec3(1, -5, 0),
+                           glm::vec3(2, 1, 0),
+                           glm::vec3(4, -1, 0)};
+    crv1.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+    crv1.degree = 2;
+    crv1.weights = {1, 1, 1, 1, 1};
 
-    ShowCurve_Igl(viewer, crv, 100);
-    BVH_AABB bvh_curve(8, 2, 10, crv);
-    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);
+    ShowCurve_Igl(crv1, 5000, YELLOW);
+    BVH_AABB bvh_curve1(11, 2, 100, crv1);
+    // 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};
-     */
+    tinynurbs::RationalCurve<double> crv2;      // Planar curve using float32
+    crv2.control_points = {glm::vec3(-4, 4, 0), // std::vector of 3D points
+                           glm::vec3(-2, 2, 0),
+                           glm::vec3(1, -1, 0),
+                           glm::vec3(3, -3, 0),
+                           glm::vec3(4, -4, 0)};
+    crv2.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+    crv2.degree = 2;
+    crv2.weights = {1, 1, 1, 1, 1};
+
+    ShowCurve_Igl(crv2, 5000, YELLOW);
+    BVH_AABB bvh_curve2(11, 2, 100, crv2);
+
+    double tstep1 = *(crv1.knots.end() - 1);
+    double tstep2 = *(crv2.knots.end() - 1);
+    bvh_curve1.Build_NurbsCurve(crv1, bvh_curve1.bvh_aabb_node, 0, 0, tstep1);
+    bvh_curve2.Build_NurbsCurve(crv2, bvh_curve2.bvh_aabb_node, 0, 0, tstep2);
+
+    std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>> IstNodePtr;
+    if (CurveCurveBVHIntersect(bvh_curve1.bvh_aabb_node, bvh_curve2.bvh_aabb_node, IstNodePtr))
+    {
+        if (IstNodePtr.size() != 0)
+            std::cout << "Curve1 and Curve2 intersect!\n";
+        else
+            std::cout << "Curve1 and Curve2 not intersect!\n";
+        int idx = 0;
+        for (auto it : IstNodePtr)
+        {
+            ShowAABB_Igl(it.first->bound, RED);
+            ShowAABB_Igl(it.second->bound, GREEN);
+
+            std::cout << "The " << idx++ << "th boxes pair lie on\n"
+                      << "     min_point ("
+                      << it.first->bound.Bmin.x << ", " << it.first->bound.Bmin.y << ", " << it.first->bound.Bmin.z << ") , max_point (" << it.first->bound.Bmax.x << ", " << it.first->bound.Bmax.y << ", " << it.first->bound.Bmax.z << ")\n";
+            std::cout << "     min_point ("
+                      << it.second->bound.Bmin.x << ", " << it.second->bound.Bmin.y << ", " << it.second->bound.Bmin.z << ") , max_point (" << it.second->bound.Bmax.x << ", " << it.second->bound.Bmax.y << ", " << it.second->bound.Bmax.z << ")\n";
+        }
+    }
 
     /*
     ShowCurve_Igl(viewer, crv, 100);

src/newton.cpp

@@ -1,155 +0,0 @@
-#include "mathdef.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")};
-    if (id == 1)
-        symbols = {get_symbol("u"), get_symbol("v")};
-    if (id == 2)
-        symbols = {get_symbol("s"), get_symbol("t")};
-    if (id == 3)
-        symbols = {get_symbol("u"), get_symbol("v"), get_symbol("s"), get_symbol("t"), get_symbol("p"), get_symbol("q")};
-    if (id == 4)
-    {
-        symbols = {};
-        for (int i = 0; i < init_point.size(); i++)
-        {
-            symbols.push_back(get_symbol("x" + to_string(i)));
-        }
-    }
-    vector<numeric> x_val = init_point;
-    vector<numeric> next_x;
-
-    auto checkValue = [&F, &id, &symbols](vector<numeric> &point)
-    {
-        ex res = 0;
-        lst l = {};
-        for (int i = 0; i < F.f.size(); i++)
-        {
-            l.append(symbols[i] == point[i]);
-        }
-        for (int i = 0; i < F.f.size(); i++)
-        {
-            ex f = F.f[i].func;
-            res += abs(evalf(f.subs(l)));
-        }
-        return res;
-    };
-    auto derivative = [&symbols](const ex &f, const symbol x, const vector<numeric> &val)
-    {
-        ex diff_f = f.diff(x);
-        lst L = {};
-        for (int i = 0; i < val.size(); i++)
-        {
-            L.append(symbols[i] == val[i]);
-        }
-
-        ex new_f = evalf(diff_f.subs(L));
-        return ex_to<numeric>(new_f);
-    };
-    auto GetJacobinMat = [&derivative, &symbols](const MUL_F &F, const vector<numeric> &xk)
-    {
-        lst l = {};
-        int n = xk.size();
-        for (int i = 0; i < F.f.size(); i++)
-        {
-            for (int j = 0; j < n; j++)
-            {
-                auto tmp = derivative(F.f[i].func, symbols[j], xk);
-                l.append(tmp);
-            }
-        }
-        matrix result = matrix(F.f.size(), n, l);
-        return result;
-    };
-    auto fn = [&symbols](const ex &f, const vector<numeric> &val)
-    {
-        lst L = {};
-        for (int i = 0; i < val.size(); i++)
-        {
-            L.append(symbols[i] == val[i]);
-        }
-        ex new_f = evalf(f.subs(L));
-        return ex_to<numeric>(new_f);
-    };
-    auto GetFk = [&symbols, &fn](const MUL_F &F, const vector<numeric> &xk)
-    {
-        lst l = {};
-        int n = F.f.size();
-        for (int i = 0; i < n; i++)
-        {
-            auto tmp = fn(F.f[i].func, xk);
-            l.append(tmp);
-        }
-        matrix result = matrix(n, 1, l);
-        return result;
-    };
-    int n = F.f.size();
-    int limit = 300;
-    isConvergence = true;
-    numeric totalDelta;
-    do
-    {
-        next_x.clear();
-        matrix grad = GetJacobinMat(F, x_val);
-
-        matrix grad_t = grad.transpose();
-        matrix f_xk = GetFk(F, x_val);
-
-        matrix JF_1;
-        try
-        {
-            JF_1 = grad.inverse();
-        }
-        catch (const std::exception &e)
-        {
-            cout << "grad = " << grad << endl;
-            JF_1 = grad_t.mul(grad);
-            isConvergence = false;
-            for (int i = 0; i < F.f.size(); i++)
-            {
-                cout << "f = " << F.f[i].func << endl;
-            }
-            cout << __FILE__ << ":" << __LINE__ << " " << e.what() << endl;
-            throw runtime_error("need to diff");
-            return x_val;
-        }
-
-        auto result = JF_1.mul(f_xk);
-        for (int i = 0; i < n; i++)
-        {
-            numeric step = (-ex_to<numeric>(result(i, 0)));
-            next_x.push_back(x_val[i] + step);
-        }
-        totalDelta = 0;
-        for (int i = 0; i < n; i++)
-        {
-            totalDelta = max(totalDelta, abs(x_val[i] - next_x[i]));
-        }
-
-        if (totalDelta < eps)
-            break;
-        x_val = next_x;
-        if (checkValue(x_val) < eps)
-            break;
-        limit--;
-        if (!limit)
-            break;
-    } while (true);
-    // cout << "limit = " << 300-limit << endl;
-    if (checkValue(x_val) > eps)
-    {
-        isConvergence = false;
-        //        cout << "totaldelta = " << totalDelta << endl;
-        //        cout << "checkValue = " << checkValue(x_val) << endl;
-    }
-    //    if (totalDelta > eps) {
-    //
-    //        isConvergence = false;
-    //    }
-    return x_val;
-}
-
-*/

src/show_libigl.cpp

@@ -1,7 +1,7 @@
 #include "bvh.hpp"
 #include "show_libigl.hpp"
 
-void ShowCurve_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalCurve<double> &curve, double sampleNum)
+void ShowCurve_Igl(tinynurbs::RationalCurve<double> &curve, double sampleNum, Eigen::RowVector3d color)
 {
     double T = *(curve.knots.end() - 1);
 
@@ -10,11 +10,10 @@ void ShowCurve_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalCurve<d
         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));
+        viewer.data().add_points(pt, color);
     }
 }
-
-void ShowSurface_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV)
+void ShowSurface_Igl(tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV, Eigen::RowVector3d color)
 {
     double U = *(surface.knots_u.end() - 1);
     double V = *(surface.knots_v.end() - 1);
@@ -26,12 +25,12 @@ void ShowSurface_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalSurfa
             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));
+            viewer.data().add_points(pt, color);
         }
     }
 }
 
-void ShowAABB(igl::opengl::glfw::Viewer &viewer, AABB bound)
+void ShowAABB_Igl(AABB bound, Eigen::RowVector3d color)
 {
 
     Eigen::Vector3d m;
@@ -67,19 +66,19 @@ void ShowAABB(igl::opengl::glfw::Viewer &viewer, AABB bound)
         viewer.data().add_edges(
             V_box.row(E_box(i, 0)),
             V_box.row(E_box(i, 1)),
-            Eigen::RowVector3d(1, 1, 1));
+            color);
     return;
 }
 
-void ShowBVHNode(igl::opengl::glfw::Viewer &viewer, BVH_AABB_NodePtr bvhNode)
+void ShowBVHNode_Igl(BVH_AABB_NodePtr bvhNode, Eigen::RowVector3d color)
 {
     if (bvhNode == NULL)
         return;
 
-    ShowAABB(viewer, bvhNode->bound);
+    ShowAABB_Igl(bvhNode->bound, color);
 
     for (auto it : bvhNode->childPtr)
     {
-        ShowBVHNode(viewer, it);
+        ShowBVHNode_Igl(it, color);
     }
 }