diff --git a/include/bvh.hpp b/include/bvh.hpp index cd42266..5fd0204 100644 --- a/include/bvh.hpp +++ b/include/bvh.hpp @@ -4,6 +4,8 @@ #include #include +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包围盒的结点 diff --git a/include/intersection.hpp b/include/intersection.hpp new file mode 100644 index 0000000..9fd77e7 --- /dev/null +++ b/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> &IstNodePtr); + +#endif \ No newline at end of file diff --git a/include/show_libigl.hpp b/include/show_libigl.hpp index cb91271..b501672 100644 --- a/include/show_libigl.hpp +++ b/include/show_libigl.hpp @@ -4,9 +4,19 @@ #include "bvh.hpp" #include -void ShowCurve_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalCurve &curve, double sampleNum); -void ShowSurface_Igl(igl::opengl::glfw::Viewer &viewer, tinynurbs::RationalSurface &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 &curve, double sampleNum, Eigen::RowVector3d color); +void ShowSurface_Igl(tinynurbs::RationalSurface &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 \ No newline at end of file diff --git a/src/intersection.cpp b/src/intersection.cpp index 3ddb503..d6f280c 100644 --- a/src/intersection.cpp +++ b/src/intersection.cpp @@ -1,17 +1,33 @@ -#include "bvh.cpp" -#include -#include +#include "intersection.hpp" +#include "show_libigl.hpp" bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2, std::vector> &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 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 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 { diff --git a/src/main.cpp b/src/main.cpp index 6b360f1..defa5f3 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -2,40 +2,68 @@ #include #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 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 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 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 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> 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); diff --git a/src/newton.cpp b/src/newton.cpp deleted file mode 100644 index 1ab9bfc..0000000 --- a/src/newton.cpp +++ /dev/null @@ -1,155 +0,0 @@ -#include "mathdef.hpp" -/* - -// 默认id参数不用传，表示四个变量四个方程，id=1 or 2 的时候是2x2的方程组， id=3时是6x6的方程组, id=4 (x0,x1,x2... xn-1) 表示的方程组 -vector newton(const MUL_F &F, const vector &init_point, const numeric &eps, bool &isConvergence, int id) -{ - vector 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 x_val = init_point; - vector next_x; - - auto checkValue = [&F, &id, &symbols](vector &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 &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(new_f); - }; - auto GetJacobinMat = [&derivative, &symbols](const MUL_F &F, const vector &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 &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(new_f); - }; - auto GetFk = [&symbols, &fn](const MUL_F &F, const vector &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(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; -} - -*/ \ No newline at end of file diff --git a/src/show_libigl.cpp b/src/show_libigl.cpp index 1eda6f4..b9b3c0c 100644 --- a/src/show_libigl.cpp +++ b/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 &curve, double sampleNum) +void ShowCurve_Igl(tinynurbs::RationalCurve &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 &surface, double sampleNumU, double sampleNumV) +void ShowSurface_Igl(tinynurbs::RationalSurface &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); } }