From 50676d343aee3ada3270ad4d7d3ddf0006bde6c8 Mon Sep 17 00:00:00 2001
From: forty-twoo <1013417276@qq.com>
Date: Sun, 4 Dec 2022 19:28:46 +0800
Subject: [PATCH] add curve surface intersection code
---
CMakeLists.txt | 2 +-
.../CMakeLists.txt | 4 +-
.../main.cpp | 209 +++++++------
.../CMakeLists.txt | 22 ++
.../main.cpp | 131 +++++++++
include/intersection.hpp | 2 +-
include/newton.hpp | 4 +-
include/show_libigl.hpp | 4 +-
src/bvh.cpp | 178 +++++------
src/intersection.cpp | 135 +++++----
src/newton.cpp | 276 +++++++++++++++++-
src/show_libigl.cpp | 146 +++++----
12 files changed, 785 insertions(+), 328 deletions(-)
create mode 100644 examples/example4_curve_surface_intersection/CMakeLists.txt
create mode 100644 examples/example4_curve_surface_intersection/main.cpp
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6ff2313..0f5398f 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -40,4 +40,4 @@ target_link_libraries(${PROJECT_NAME} PUBLIC igl::glfw)
# add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/example1_show_bvh_curve)
# add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/example2_show_bvh_surface)
add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/example3_curve_curve_intersection)
-# add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/example4_curve_surface_intersection)
+add_subdirectory(${CMAKE_CURRENT_SOURCE_DIR}/examples/example4_curve_surface_intersection)
diff --git a/examples/example3_curve_curve_intersection/CMakeLists.txt b/examples/example3_curve_curve_intersection/CMakeLists.txt
index c57785c..4c66716 100644
--- a/examples/example3_curve_curve_intersection/CMakeLists.txt
+++ b/examples/example3_curve_curve_intersection/CMakeLists.txt
@@ -7,10 +7,10 @@ set(CMAKE_CXX_STANDARD 11)
# list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
# Libigl
-# include(libigl)
+include(libigl)
# Enable the target igl::glfw
-# igl_include(glfw)
+igl_include(glfw)
# add include directories
diff --git a/examples/example3_curve_curve_intersection/main.cpp b/examples/example3_curve_curve_intersection/main.cpp
index d7b752b..6e24544 100644
--- a/examples/example3_curve_curve_intersection/main.cpp
+++ b/examples/example3_curve_curve_intersection/main.cpp
@@ -10,103 +10,124 @@
#include "bvh.hpp"
#include "show_libigl.hpp"
#include "intersection.hpp"
+#include "newton.hpp"
igl::opengl::glfw::Viewer viewer;
-int main(int argc, char *argv[])
+int main(int argc, char* argv[])
{
- viewer.data().point_size = 5;
-
- 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};
-
- double param = 0.6;
- tinynurbs::array2<double> ders = tinynurbs::bsplineDerBasis(crv1.degree, tinynurbs::findSpan(crv1.degree, crv1.knots, param), crv1.knots, param, 1);
- std::cout << "row:" << ders.rows() << " col:" << ders.cols() << std::endl;
- for (int i = 0; i < ders.rows(); i++)
- {
- for (int j = 0; j < ders.cols(); j++)
- {
- std::cout << ders(i, j) << " ";
- }
- std::cout << std::endl;
- }
-
- std::vector<double> basis = tinynurbs::bsplineBasis(crv1.degree, tinynurbs::findSpan(crv1.degree, crv1.knots, param), crv1.knots, param);
- for (auto it : basis)
- {
- std::cout << it << std::endl;
- }
-
- /***
- 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> 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);
- 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);
- */
-
- /*
- double u = 0.5;
- std::vector<double> basis = tinynurbs::bsplineBasis(crv.degree, tinynurbs::findSpan(crv.degree, crv.knots, u), crv.knots, u);
- for (auto it : basis)
- {
- std::cout << it << std::endl;
- }
- */
-
- viewer.launch();
+ viewer.data().point_size = 5;
+
+ 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),
+ glm::vec3(6, 2, 0),
+ glm::vec3(8, 4, 0)
+
+ };
+ crv1.knots = { 0, 0, 0, 1, 2, 3, 4, 5, 5, 5 }; // std::vector of floats
+ crv1.degree = 2;
+ crv1.weights = { 1, 1, 1, 1, 1, 1, 1 };
+
+
+ ShowCurve_Igl(crv1, 500, YELLOW);
+ BVH_AABB bvh_curve1(12, 2, 100, crv1);
+
+ tinynurbs::RationalCurve<double> crv2; // Planar curve using float32
+ crv2.control_points = { glm::vec3(-4, 4, 0), // std::vector of 3D points
+ glm::vec3(-5, 2, 0),
+ glm::vec3(-4, 6, 0),
+ glm::vec3(-2, -5, 0),
+ glm::vec3(1, 8, 0),
+ glm::vec3(5, -7, 0)
+ };
+ crv2.knots = { 0, 0, 0, 1, 2, 3, 4, 4, 4 }; // std::vector of floats
+ crv2.degree = 2;
+ crv2.weights = { 1, 1, 1, 1, 1, 1 };
+
+ ShowCurve_Igl(crv2, 500, YELLOW);
+ BVH_AABB bvh_curve2(12, 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);
+
+ //ShowBVHNode_Igl(bvh_curve1.bvh_aabb_node, WHITE);
+ //ShowBVHNode_Igl(bvh_curve2.bvh_aabb_node, GREEN);
+
+
+ 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;
+ std::vector<glm::vec3> ansPoints;
+ 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";
+
+ auto itsPoint = newton_curve_curve(it.first, it.second, 1e-12, true);
+ std::cout << "\n\n\n";
+ //ȥ��
+ if(ansPoints.size() == 0)
+ ansPoints.push_back(itsPoint);
+ else
+ {
+ //���̫��������ͬһ����
+ glm::dvec3 delt = itsPoint - ansPoints[ansPoints.size() - 1];
+ double dis = std::sqrt(delt.x * delt.x + delt.y * delt.y + delt.z * delt.z);
+ if(dis > 1e-5)
+ {
+ ansPoints.push_back(itsPoint);
+ std::cout << "distance between last and now: ";
+ std::cout << std::setprecision(10) << dis << "\n";
+ }
+
+ }
+
+ }
+ for(auto it : ansPoints)
+ {
+ viewer.data().add_points(Eigen::RowVector3d(it.x, it.y, it.z), BLACK);
+ std::cout << "(" << it.x << "," << it.y << "," << it.z << ")" << std::endl;
+ }
+
+
+
+
+ }
+
+ /*
+ 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);
+ */
+
+ /*
+ double u = 0.5;
+ std::vector<double> basis = tinynurbs::bsplineBasis(crv.degree, tinynurbs::findSpan(crv.degree, crv.knots, u), crv.knots, u);
+ for (auto it : basis)
+ {
+ std::cout << it << std::endl;
+ }
+ */
+
+ viewer.launch();
}
diff --git a/examples/example4_curve_surface_intersection/CMakeLists.txt b/examples/example4_curve_surface_intersection/CMakeLists.txt
new file mode 100644
index 0000000..2ebad0a
--- /dev/null
+++ b/examples/example4_curve_surface_intersection/CMakeLists.txt
@@ -0,0 +1,22 @@
+cmake_minimum_required(VERSION 3.16)
+project(example4_curve_surface_intersection)
+
+# C++ 11 is required
+set(CMAKE_CXX_STANDARD 11)
+
+# list(PREPEND CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake)
+
+# Libigl
+include(libigl)
+
+# Enable the target igl::glfw
+igl_include(glfw)
+
+
+# add include directories
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/include)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/3rdparty)
+
+add_executable(${PROJECT_NAME} main.cpp)
+
+target_link_libraries(${PROJECT_NAME} NurbsIntersection::nurbsintersect)
diff --git a/examples/example4_curve_surface_intersection/main.cpp b/examples/example4_curve_surface_intersection/main.cpp
new file mode 100644
index 0000000..87c4b06
--- /dev/null
+++ b/examples/example4_curve_surface_intersection/main.cpp
@@ -0,0 +1,131 @@
+/*
+ * @File Created: 2022-11-26, 16:52:29
+ * @Last Modified: 2022-11-26, 18:51:06
+ * @Author: forty-twoo
+ * @Copyright (c) 2022, Caiyue Li(li_caiyue@zju.edu.cn), All rights reserved.
+ */
+
+#include <igl/opengl/glfw/Viewer.h>
+#include <tinynurbs/tinynurbs.h>
+#include "bvh.hpp"
+#include "show_libigl.hpp"
+#include "intersection.hpp"
+#include "newton.hpp"
+
+igl::opengl::glfw::Viewer viewer;
+
+int main(int argc, char* argv[])
+{
+
+ viewer.data().point_size = 5;
+
+ tinynurbs::RationalCurve<double> crv; // Planar curve using float32
+ crv.control_points = { glm::vec3(-4, 0, 0), // std::vector of 3D points
+ glm::vec3(-1, -1.2, -0.5),
+ glm::vec3(2, 2.3, 1),
+ glm::vec3(-3, -1.8, -2),
+
+ };
+ crv.degree = 2;
+ crv.knots = { 0, 0, 0, 1, 2, 2, 2 }; // std::vector of floats
+ crv.weights = { 1, 1, 1, 1 };
+
+
+ tinynurbs::RationalSurface<double> srf;
+ srf.degree_u = 3;
+ srf.degree_v = 3;
+ srf.knots_u = { 0, 0, 0, 0, 1, 1, 1, 1 };
+ srf.knots_v = { 0, 0, 0, 0, 1, 1, 1, 1 };
+
+ srf.control_points = { 4, 4, {
+ glm::vec3(2, -2, -2), glm::vec3(-2.5, -2.2, -1.5), glm::vec3(-2, -2, -0.5), glm::vec3(-2, -2, 1.5),
+ glm::vec3(2, -1, -2), glm::vec3(-2.5, -1.2, -1.5), glm::vec3(-2, -1, -0.5), glm::vec3(-2, -1, 1.5),
+ glm::vec3(3, 1.2, -2), glm::vec3(-2.5, 1.2, -1.5), glm::vec3(-2, 1.5, -0.5), glm::vec3(-2, 1.5, 1.5),
+ glm::vec3(2, 2, -2), glm::vec3(-2.5, 2, -1.5), glm::vec3(-2, 2.5, -0.5), glm::vec3(-2, 2, 1.5)
+ }
+ };
+
+ srf.weights = { 4, 4,
+ {
+ 1, 1, 1, 1,
+ 1, 1, 1, 1,
+ 1, 1, 1, 1,
+ 1, 1, 1, 1,
+ }
+ };
+
+
+ BVH_AABB bvh_curve(10, 2, 100, crv);
+ BVH_AABB bvh_surface(6, 4, 100, srf);
+ ShowCurve_Igl(crv, 500, YELLOW);
+ ShowSurface_Igl(srf, 10, 10, BLUE);
+
+
+
+ double tstep = *(crv.knots.end() - 1);
+ bvh_curve.Build_NurbsCurve(crv, bvh_curve.bvh_aabb_node, 0, 0, tstep);
+
+ double ustep = *(srf.knots_u.end() - 1);
+ double vstep = *(srf.knots_v.end() - 1);
+ bvh_surface.Build_NurbsSurface(srf, bvh_surface.bvh_aabb_node, 0, 0, 0, ustep, vstep);
+
+ //ShowBVHNode_Igl(bvh_curve.bvh_aabb_node, RED);
+ //ShowBVHNode_Igl(bvh_surface.bvh_aabb_node, GREEN);
+
+
+ std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>> IstNodePtr;
+ if(BVHIntersect(bvh_curve.bvh_aabb_node, bvh_surface.bvh_aabb_node, IstNodePtr))
+ {
+ if(IstNodePtr.size() != 0)
+ std::cout << "Curve and Surface intersect!\n";
+ else
+ std::cout << "Curve and Surface not intersect!\n";
+ int idx = 0;
+ std::vector<glm::vec3> ansPoints;
+
+ 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";
+
+
+ auto itsPoint = newton_curve_surface(it.first, it.second, 1e-12, true);
+ ansPoints.push_back(itsPoint);
+ std::cout << "\n\n\n";
+ //std::cout << "(" << itsPoint.x << "," << itsPoint.y << "," << itsPoint.z << ")" << std::endl;
+ //ȥ��
+ //else
+ //{
+ // //���̫��������ͬһ����
+ // glm::dvec3 delt = itsPoint - ansPoints[ansPoints.size() - 1];
+ // double dis = std::sqrt(delt.x * delt.x + delt.y * delt.y + delt.z * delt.z);
+ // ansPoints.push_back(itsPoint);
+
+ // //if(dis > 1e-5)
+ // //{
+ // // std::cout << "distance between last and now: ";
+ // // std::cout << std::setprecision(10) << dis << "\n";
+ // //}
+
+ //}
+
+ }
+ for(auto it : ansPoints)
+ {
+ viewer.data().add_points(Eigen::RowVector3d(it.x, it.y, it.z), BLACK);
+ std::cout << "(" << it.x << "," << it.y << "," << it.z << ")" << std::endl;
+ }
+ }
+
+
+ viewer.launch();
+ return 0;
+
+
+}
diff --git a/include/intersection.hpp b/include/intersection.hpp
index f992a61..39508d9 100644
--- a/include/intersection.hpp
+++ b/include/intersection.hpp
@@ -3,6 +3,6 @@
#include "bvh.hpp"
-bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2, std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>> &IstNodePtr);
+bool BVHIntersect(BVH_AABB_NodePtr& BoxPtr1, BVH_AABB_NodePtr& BoxPtr2, std::vector<std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr>>& IstNodePtr);
#endif
\ No newline at end of file
diff --git a/include/newton.hpp b/include/newton.hpp
index 8b7b554..aabfe87 100644
--- a/include/newton.hpp
+++ b/include/newton.hpp
@@ -9,7 +9,9 @@
#define NEWTON_H_
#include "bvh.hpp"
+#include <Eigen/Core>
-std::vector<double> newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double &eps, bool &isConvergence);
+glm::vec3 newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double& eps, bool isConvergence);
+glm::vec3 newton_curve_surface(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double& eps, bool isConvergence);
#endif
\ No newline at end of file
diff --git a/include/show_libigl.hpp b/include/show_libigl.hpp
index 518c900..73aa247 100644
--- a/include/show_libigl.hpp
+++ b/include/show_libigl.hpp
@@ -21,8 +21,8 @@ extern igl::opengl::glfw::Viewer viewer;
#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 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);
diff --git a/src/bvh.cpp b/src/bvh.cpp
index f0042ea..7d14bb9 100644
--- a/src/bvh.cpp
+++ b/src/bvh.cpp
@@ -1,95 +1,107 @@
#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),
- fmin(a.Bmin.y, b.Bmin.y),
- fmin(a.Bmin.z, b.Bmin.z));
-
- ret.Bmax = glm::vec3(fmax(a.Bmax.x, b.Bmax.x),
- fmax(a.Bmax.y, b.Bmax.y),
- fmax(a.Bmax.z, b.Bmax.z));
- return ret;
+ AABB ret;
+ ret.Bmin = glm::vec3(fmin(a.Bmin.x, b.Bmin.x),
+ fmin(a.Bmin.y, b.Bmin.y),
+ fmin(a.Bmin.z, b.Bmin.z));
+
+ ret.Bmax = glm::vec3(fmax(a.Bmax.x, b.Bmax.x),
+ fmax(a.Bmax.y, b.Bmax.y),
+ fmax(a.Bmax.z, b.Bmax.z));
+ 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),
- fmin(bound.Bmin.y, p.y),
- fmin(bound.Bmin.z, p.z));
-
- ret.Bmax = glm::vec3(fmax(bound.Bmax.x, p.x),
- fmax(bound.Bmax.y, p.y),
- fmax(bound.Bmax.z, p.z));
- return ret;
+ AABB ret;
+ ret.Bmin = glm::vec3(fmin(bound.Bmin.x, p.x),
+ fmin(bound.Bmin.y, p.y),
+ fmin(bound.Bmin.z, p.z));
+
+ ret.Bmax = glm::vec3(fmax(bound.Bmax.x, p.x),
+ fmax(bound.Bmax.y, p.y),
+ fmax(bound.Bmax.z, p.z));
+ return ret;
}
-BVH_AABB_NodePtr BVH_AABB::Build_NurbsCurve(tinynurbs::RationalCurve<double> &curve, BVH_AABB_NodePtr &curNode, int curLevel, double t, double tstep)
+BVH_AABB_NodePtr BVH_AABB::Build_NurbsCurve(tinynurbs::RationalCurve<double>& curve, BVH_AABB_NodePtr& curNode, int curLevel, double t, double tstep)
{
- curNode->myBVH = this;
- if (curLevel + 1 == Hlevel)
- {
- for (double i = t; i < t + tstep; i += tstep / sampleNum)
- {
- curNode->bound = UpdateBound(curNode->bound, tinynurbs::curvePoint(curve, i));
- }
- curNode->childPtr.clear();
- curNode->param.push_back(std::make_pair(t, t + tstep));
-
- return curNode;
- }
- BVH_AABB_NodePtr childNode1 = new BVH_AABB_Node;
- BVH_AABB_NodePtr childNode2 = new BVH_AABB_Node;
-
- double tstep_ = tstep / (double)BranchNum;
- Build_NurbsCurve(curve, childNode1, curLevel + 1, t, tstep_);
- Build_NurbsCurve(curve, childNode2, curLevel + 1, t + tstep_, tstep_);
-
- curNode->bound = UnionBound(curNode->bound, childNode1->bound);
- curNode->bound = UnionBound(curNode->bound, childNode2->bound);
-
- curNode->childPtr.push_back(childNode1);
- curNode->childPtr.push_back(childNode2);
-
- curNode->param.push_back(std::make_pair(t, t + tstep));
-
- return curNode;
+ curNode->myBVH = this;
+ if(curLevel + 1 == Hlevel)
+ {
+ for(double i = t; i < t + tstep; i += tstep / sampleNum)
+ {
+ curNode->bound = UpdateBound(curNode->bound, tinynurbs::curvePoint(curve, i));
+ }
+ curNode->childPtr.clear();
+ curNode->param.push_back(std::make_pair(t, t + tstep));
+
+ return curNode;
+ }
+ BVH_AABB_NodePtr childNode1 = new BVH_AABB_Node;
+ BVH_AABB_NodePtr childNode2 = new BVH_AABB_Node;
+
+ double tstep_ = tstep / (double)BranchNum;
+ Build_NurbsCurve(curve, childNode1, curLevel + 1, t, tstep_);
+ Build_NurbsCurve(curve, childNode2, curLevel + 1, t + tstep_, tstep_);
+
+ curNode->bound = UnionBound(curNode->bound, childNode1->bound);
+ curNode->bound = UnionBound(curNode->bound, childNode2->bound);
+
+ curNode->childPtr.push_back(childNode1);
+ curNode->childPtr.push_back(childNode2);
+
+ curNode->param.push_back(std::make_pair(t, t + tstep));
+
+ return curNode;
}
-BVH_AABB_NodePtr BVH_AABB::Build_NurbsSurface(tinynurbs::RationalSurface<double> &surface, BVH_AABB_NodePtr &curNode, int curLevel, double u, double v, double ustep, double vstep)
+BVH_AABB_NodePtr BVH_AABB::Build_NurbsSurface(tinynurbs::RationalSurface<double>& surface, BVH_AABB_NodePtr& curNode, int curLevel, double u, double v, double ustep, double vstep)
{
- curNode->myBVH = this;
- if (curLevel == Hlevel)
- {
- for (double i = u; i < u + ustep; i++)
- {
- for (double j = v; j < v + vstep; j++)
- {
- curNode->bound = UpdateBound(curNode->bound, tinynurbs::surfacePoint(surface, i, j));
- }
- }
- curNode->childPtr.clear();
- curNode->param.push_back(std::make_pair(u, u + ustep));
- curNode->param.push_back(std::make_pair(v, v + vstep));
-
- return curNode;
- }
- BVH_AABB_NodePtr childNode1 = new BVH_AABB_Node;
- BVH_AABB_NodePtr childNode2 = new BVH_AABB_Node;
- BVH_AABB_NodePtr childNode3 = new BVH_AABB_Node;
- BVH_AABB_NodePtr childNode4 = new BVH_AABB_Node;
-
- double ustep_ = ustep / (double)BranchNum;
- double vstep_ = vstep / (double)BranchNum;
-
- Build_NurbsSurface(surface, childNode1, curLevel + 1, u, v, ustep_, vstep_);
- Build_NurbsSurface(surface, childNode2, curLevel + 1, u + ustep_, v, ustep_, vstep_);
- Build_NurbsSurface(surface, childNode3, curLevel + 1, u, v + vstep_, ustep_, vstep_);
- Build_NurbsSurface(surface, childNode4, curLevel + 1, u + ustep_, v + vstep_, ustep_, vstep_);
-
- curNode->param.push_back(std::make_pair(u, u + ustep));
- curNode->param.push_back(std::make_pair(v, v + vstep));
-
- return curNode;
+ curNode->myBVH = this;
+ if(curLevel + 1 == Hlevel)
+ {
+ for(double i = u; i < u + ustep; i += ustep / sampleNum)
+ {
+ for(double j = v; j < v + vstep; j += vstep / sampleNum)
+ {
+ curNode->bound = UpdateBound(curNode->bound, tinynurbs::surfacePoint(surface, i, j));
+ }
+ }
+ curNode->childPtr.clear();
+ //std::cout << curNode->bound.Bmin.x << "," << curNode->bound.Bmin.y << "," << curNode->bound.Bmin.z << "\n";
+ curNode->param.push_back(std::make_pair(u, u + ustep));
+ curNode->param.push_back(std::make_pair(v, v + vstep));
+
+ return curNode;
+ }
+ BVH_AABB_NodePtr childNode1 = new BVH_AABB_Node;
+ BVH_AABB_NodePtr childNode2 = new BVH_AABB_Node;
+ BVH_AABB_NodePtr childNode3 = new BVH_AABB_Node;
+ BVH_AABB_NodePtr childNode4 = new BVH_AABB_Node;
+
+ double ustep_ = ustep / (double)2;
+ double vstep_ = vstep / (double)2;
+
+ Build_NurbsSurface(surface, childNode1, curLevel + 1, u, v, ustep_, vstep_);
+ Build_NurbsSurface(surface, childNode2, curLevel + 1, u + ustep_, v, ustep_, vstep_);
+ Build_NurbsSurface(surface, childNode3, curLevel + 1, u, v + vstep_, ustep_, vstep_);
+ Build_NurbsSurface(surface, childNode4, curLevel + 1, u + ustep_, v + vstep_, ustep_, vstep_);
+
+ curNode->bound = UnionBound(curNode->bound, childNode1->bound);
+ curNode->bound = UnionBound(curNode->bound, childNode2->bound);
+ curNode->bound = UnionBound(curNode->bound, childNode3->bound);
+ curNode->bound = UnionBound(curNode->bound, childNode4->bound);
+
+ curNode->childPtr.push_back(childNode1);
+ curNode->childPtr.push_back(childNode2);
+ curNode->childPtr.push_back(childNode3);
+ curNode->childPtr.push_back(childNode4);
+
+ curNode->param.push_back(std::make_pair(u, u + ustep));
+ curNode->param.push_back(std::make_pair(v, v + vstep));
+
+ //std::cout << curNode->bound.Bmin.x << "," << curNode->bound.Bmin.y << "," << curNode->bound.Bmin.z << "\n";
+ return curNode;
}
\ No newline at end of file
diff --git a/src/intersection.cpp b/src/intersection.cpp
index 92845fd..a8b8469 100644
--- a/src/intersection.cpp
+++ b/src/intersection.cpp
@@ -8,77 +8,74 @@
#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)
+bool BVHIntersect(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;
+ AABB box1 = BoxPtr1->bound;
+ AABB box2 = BoxPtr2->bound;
- auto LineIntersect = [](double x1, double x2, double x1_, double x2_)
- {
- bool LIntersect = true;
- if (x1 <= x1_)
- {
- if (x1_ > x2)
- LIntersect = false;
- }
- else
- {
- if (x1 > x2_)
- LIntersect = false;
- }
- return LIntersect;
- };
- auto BoxIntersect = [&]()
- {
- double x1 = box1.Bmin.x, x1_ = box2.Bmin.x;
- double x2 = box1.Bmax.x, x2_ = box2.Bmax.x;
- double y1 = box1.Bmin.y, y1_ = box2.Bmin.y;
- double y2 = box1.Bmax.y, y2_ = box2.Bmax.y;
- double z1 = box1.Bmin.z, z1_ = box2.Bmin.z;
- double z2 = box1.Bmax.z, z2_ = box2.Bmax.z;
- if (LineIntersect(x1, x2, x1_, x2_) && LineIntersect(y1, y2, y1_, y2_) && LineIntersect(z1, z2, z1_, z2_))
- return true;
- else
- false;
- };
+ auto LineIntersect = [](double x1, double x2, double x1_, double x2_)
+ {
+ bool LIntersect = true;
+ if(x1 <= x1_)
+ {
+ if(x1_ > x2)
+ LIntersect = false;
+ } else
+ {
+ if(x1 > x2_)
+ LIntersect = false;
+ }
+ return LIntersect;
+ };
+ auto BoxIntersect = [&]()
+ {
+ double x1 = box1.Bmin.x, x1_ = box2.Bmin.x;
+ double x2 = box1.Bmax.x, x2_ = box2.Bmax.x;
+ double y1 = box1.Bmin.y, y1_ = box2.Bmin.y;
+ double y2 = box1.Bmax.y, y2_ = box2.Bmax.y;
+ double z1 = box1.Bmin.z, z1_ = box2.Bmin.z;
+ double z2 = box1.Bmax.z, z2_ = box2.Bmax.z;
+ if(LineIntersect(x1, x2, x1_, x2_) && LineIntersect(y1, y2, y1_, y2_) && LineIntersect(z1, z2, z1_, z2_))
+ return true;
+ else
+ false;
+ };
- int biggerbox = 1;
- if (glm::distance(box1.Bmax, box1.Bmin) < glm::distance(box2.Bmax, box2.Bmin))
- biggerbox = 2;
+ int biggerbox = 1;
+ if(glm::distance(box1.Bmax, box1.Bmin) < glm::distance(box2.Bmax, box2.Bmin))
+ biggerbox = 2;
- 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);
- }
- }
- else
- {
- 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;
- }
- for (auto it : BoxPtr2->childPtr)
- {
- CurveCurveBVHIntersect(BoxPtr1, it, IstNodePtr);
- }
- }
- return true;
- }
- else
- return false;
+ 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)
+ {
+ BVHIntersect(it, BoxPtr2, IstNodePtr);
+ }
+ } else
+ {
+ if(BoxPtr2->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 : BoxPtr2->childPtr)
+ {
+ BVHIntersect(BoxPtr1, it, IstNodePtr);
+ }
+ }
+ return true;
+ } else
+ return false;
}
\ No newline at end of file
diff --git a/src/newton.cpp b/src/newton.cpp
index 5b80b95..16e1871 100644
--- a/src/newton.cpp
+++ b/src/newton.cpp
@@ -6,26 +6,274 @@
*/
#include "newton.hpp"
-#include <Eigen/Core>
+#include <tinynurbs/tinynurbs.h>
+#include <iomanip>
+#include <Eigen/Dense>
-std::vector<double> newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double &eps, bool &isConvergence)
+glm::vec3 newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double& eps, bool isConvergence)
{
- double param1_st = boxPtr1->param[0].first;
- double param1_ed = boxPtr1->param[0].second;
+ double param1_st = boxPtr1->param[0].first;
+ double param1_ed = boxPtr1->param[0].second;
- double param2_st = boxPtr2->param[0].first;
- double param2_ed = boxPtr2->param[0].second;
+ double param2_st = boxPtr2->param[0].first;
+ double param2_ed = boxPtr2->param[0].second;
- double t1 = (param1_st + param1_ed) / 2.0;
- double t2 = (param2_st + param2_st) / 2.0;
+ double t1 = (param1_st + param1_ed) / 2.0;
- Eigen::MatrixXd F(2, 2);
+ double t2 = (param2_st + param2_ed) / 2.0;
- tinynurbs::RationalCurve<double> *crvPtr1 = boxPtr1->myBVH->NurbsCurvePtr;
- tinynurbs::array2<double> der_Fx = tinynurbs::bsplineDerBasis(crvPtr1->degree, tinynurbs::findSpan(crvPtr1->degree, crvPtr1->knots, t1), crvPtr1->knots, t1, 1);
- Eigen::Vector2d t(t1, t2);
+ Eigen::Matrix2d F(2, 2);
+ Eigen::Vector2d t, t_, Ft;
+ t << t1, t2;
- Eigen::MatrixXd JacobinF(2, 2);
+ tinynurbs::RationalCurve<double>* crvPtr1 = boxPtr1->myBVH->NurbsCurvePtr;
+ tinynurbs::RationalCurve<double>* crvPtr2 = boxPtr2->myBVH->NurbsCurvePtr;
- bool iterFlag = true;
+ int span_t1, span_t2, idx_1, idx_2;
+
+
+ tinynurbs::array2<double> der_ft1, der_gt2;
+
+ auto controlPts_vec1 = crvPtr1->control_points;
+ auto controlPts_vec2 = crvPtr2->control_points;
+
+ glm::vec3 derff(0), dergg(0);
+ Eigen::Matrix2d JacobinF, J_Inverse;
+
+ auto showglm = [](glm::vec3 p)
+ {
+ std::cout << "(" << p.x << "," << p.y << "," << p.z << ")";
+
+ };
+
+ bool firstItr = true, iterFlag = true;
+ glm::vec3 Pt_curve1(0), Pt_curve2(0);
+ int cnt = 0;
+ do
+ {
+ std::cout << "t : (" << t(0) << "," << t(1) << ")\n";
+ Pt_curve1 = tinynurbs::curvePoint(*crvPtr1, t(0));
+ Pt_curve2 = tinynurbs::curvePoint(*crvPtr2, t(1));
+
+ glm::dvec3 delt = Pt_curve1 - Pt_curve2;
+ double dis = std::sqrt(delt.x * delt.x + delt.y * delt.y + delt.z * delt.z);
+
+ if(dis <= eps)
+ {
+ iterFlag = false;
+ }
+ std::cout << "Curve1 point: ";
+ showglm(Pt_curve1);
+ std::cout << "\nCurve2 point: ";
+ showglm(Pt_curve2);
+
+ std::cout << std::setprecision(10) << "\nDistance on curves: " << dis << std::endl;
+
+
+ if(!firstItr)
+ {
+ t1 = t(0); t2 = t(1);
+ }
+ span_t1 = tinynurbs::findSpan(crvPtr1->degree, crvPtr1->knots, t1);
+ idx_1 = span_t1 - crvPtr1->degree - 1;
+
+ span_t2 = tinynurbs::findSpan(crvPtr2->degree, crvPtr2->knots, t2);
+ idx_2 = span_t2 - crvPtr2->degree - 1;
+
+ std::cout << "span_t1: " << span_t1 << " idx_1: " << idx_1 << " span_t2: " << span_t2 << " idx_2: " << idx_2 << '\n';
+
+ der_ft1 = tinynurbs::bsplineDerBasis(crvPtr1->degree, span_t1, crvPtr1->knots, t1, 1);
+ der_gt2 = tinynurbs::bsplineDerBasis(crvPtr2->degree, span_t2, crvPtr2->knots, t2, 1);
+
+ derff = glm::vec3(0);
+ dergg = glm::vec3(0);
+ int col1 = der_ft1.cols();
+ int col2 = der_gt2.cols();
+ for(int i = idx_1; i < idx_1 + col1; i++)
+ {
+ derff += controlPts_vec1[i] * der_ft1(1, (i - idx_1));
+ }
+ for(int i = idx_2; i < idx_2 + col2; i++)
+ {
+ dergg += controlPts_vec2[i] * der_gt2(1, (i - idx_2));
+ }
+
+ Ft << (Pt_curve1.x + Pt_curve1.z - Pt_curve2.x - Pt_curve2.z), (Pt_curve1.y + Pt_curve1.z - Pt_curve2.y - Pt_curve2.z);
+
+ JacobinF(0, 0) = derff.x + derff.z;
+ JacobinF(0, 1) = -dergg.x - dergg.z;
+ JacobinF(1, 0) = derff.y + derff.z;
+ JacobinF(1, 1) = -dergg.y - dergg.z;
+
+ std::cout << "Jacobin Mat :\n" << JacobinF << '\n';
+ J_Inverse = JacobinF.inverse();
+ t_ = t - J_Inverse * Ft;
+ t = t_;
+ firstItr = false;
+ cnt++;
+ std::cout << "The " << cnt << "th iteration done!\n------------------------------------------------------\n";
+ if(cnt > 20)
+ {
+ break;
+ }
+
+ } while(iterFlag);
+
+
+ return Pt_curve1;
+}
+
+//boxPtr1->curve_BVH boxPtr2->surface_BVH
+glm::vec3 newton_curve_surface(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double& eps, bool isConvergence)
+{
+ double param1_st = boxPtr1->param[0].first;
+ double param1_ed = boxPtr1->param[0].second;
+
+ double param2_st = boxPtr2->param[0].first;
+ double param2_ed = boxPtr2->param[0].second;
+
+ double t1 = (param1_st + param1_ed) / 2.0;
+ double t2 = (param2_st + param2_ed) / 2.0;
+ std::cout << param1_st << " " << param1_ed << '\n';
+ std::cout << param2_st << " " << param2_ed << '\n';
+
+
+ double param3_st = boxPtr2->param[1].first;
+ double param3_ed = boxPtr2->param[1].second;
+
+ double t3 = (param3_st + param3_ed) / 2.0;
+
+
+
+ Eigen::Matrix3d F(3, 3);
+ Eigen::Vector3d t, t_, Ft;
+ t << t1, t2, t3;
+
+ tinynurbs::RationalCurve<double>* crvPtr1 = boxPtr1->myBVH->NurbsCurvePtr;
+ tinynurbs::RationalSurface<double>* srfPtr2 = boxPtr2->myBVH->NurbsSurfacePtr;
+
+ int span_t1, span_t2, span_t3, idx_1, idx_2, idx_3;
+
+
+ tinynurbs::array2<double> der_ft1, der_ft2, der_ft3;
+
+ auto controlPts_vec1 = crvPtr1->control_points;
+ auto controlPts_vec2 = srfPtr2->control_points;
+
+ glm::vec3 derff1(0), derff2(0), derff3(0);
+ Eigen::Matrix3d JacobinF, J_Inverse;
+
+ auto showglm = [](glm::vec3 p)
+ {
+ std::cout << "(" << p.x << "," << p.y << "," << p.z << ")";
+
+ };
+
+ bool firstItr = true, iterFlag = true;
+ glm::vec3 Pt_curve(0), Pt_surface(0);
+ int cnt = 0;
+ do
+ {
+ Pt_curve = tinynurbs::curvePoint(*crvPtr1, t(0));
+ Pt_surface = tinynurbs::surfacePoint(*srfPtr2, t(1), t(2));
+
+ glm::dvec3 delt = Pt_curve - Pt_surface;
+ double dis = std::sqrt(delt.x * delt.x + delt.y * delt.y + delt.z * delt.z);
+
+ if(dis <= eps)
+ {
+ iterFlag = false;
+ }
+ std::cout << "Curve point: ";
+ showglm(Pt_curve);
+ std::cout << "\nSurface point: ";
+ showglm(Pt_surface);
+
+ std::cout << std::setprecision(10) << "\nDistance between them : " << dis << std::endl;
+
+
+ if(!firstItr)
+ {
+ t1 = t(0); t2 = t(1), t3 = t(2);
+ }
+ span_t1 = tinynurbs::findSpan(crvPtr1->degree, crvPtr1->knots, t1);
+ idx_1 = span_t1 - crvPtr1->degree;
+
+ span_t2 = tinynurbs::findSpan(srfPtr2->degree_u, srfPtr2->knots_u, t2);
+ idx_2 = span_t2 - srfPtr2->degree_u;
+
+ span_t3 = tinynurbs::findSpan(srfPtr2->degree_v, srfPtr2->knots_v, t3);
+ idx_3 = span_t3 - srfPtr2->degree_v;
+
+
+ //std::cout << "span_t1: " << span_t1 << " idx_1: " << idx_1 << " span_t2: " << span_t2 << " idx_2: " << idx_2 << " span_t3: " << span_t3 << " idx_3: " << idx_3 << '\n';
+
+
+ auto bas_ft1 = tinynurbs::bsplineBasis(crvPtr1->degree, span_t1, crvPtr1->knots, t1);
+ auto bas_ft2 = tinynurbs::bsplineBasis(srfPtr2->degree_u, span_t2, srfPtr2->knots_u, t2);
+ auto bas_ft3 = tinynurbs::bsplineBasis(srfPtr2->degree_v, span_t3, srfPtr2->knots_v, t3);
+
+ //std::cout << "bas_ft1.size= " << bas_ft1.size();
+ //std::cout << " bas_ft2.size= " << bas_ft2.size();
+ //std::cout << " bas_ft3.size= " << bas_ft3.size() << '\n';
+
+ der_ft1 = tinynurbs::bsplineDerBasis(crvPtr1->degree, span_t1, crvPtr1->knots, t1, 1);
+ der_ft2 = tinynurbs::bsplineDerBasis(srfPtr2->degree_u, span_t2, srfPtr2->knots_u, t2, 1);
+ der_ft3 = tinynurbs::bsplineDerBasis(srfPtr2->degree_v, span_t3, srfPtr2->knots_v, t3, 1);
+
+ derff1 = glm::vec3(0);
+ derff2 = glm::vec3(0);
+ derff3 = glm::vec3(0);
+ int col1 = der_ft1.cols();
+ int col2 = der_ft2.cols();
+ int col3 = der_ft3.cols();
+ for(int i = idx_1; i < idx_1 + col1; i++)
+ {
+ derff1 += controlPts_vec1[i] * der_ft1(1, (i - idx_1));
+ }
+ for(int i = idx_2; i < idx_2 + col2; i++)
+ {
+ for(int j = idx_3; j < idx_3 + bas_ft3.size(); j++)
+ {
+ derff2 += controlPts_vec2(i, j) * der_ft2(1, (i - idx_2)) * bas_ft3[j - idx_3];
+ }
+ }
+
+ for(int i = idx_2; i < idx_2 + bas_ft2.size(); i++)
+ {
+ for(int j = idx_3; j < idx_3 + col3; j++)
+ {
+ derff3 += controlPts_vec2(i, j) * der_ft3(1, (j - idx_3)) * bas_ft2[i - idx_2];
+ }
+ }
+
+
+ Ft << (Pt_curve.x - Pt_surface.x), (Pt_curve.y - Pt_surface.y), (Pt_curve.z - Pt_surface.z);
+
+ JacobinF(0, 0) = derff1.x;
+ JacobinF(0, 1) = -derff2.x;
+ JacobinF(0, 2) = -derff3.x;
+ JacobinF(1, 0) = derff1.y;
+ JacobinF(1, 1) = -derff2.y;
+ JacobinF(1, 2) = -derff3.y;
+ JacobinF(2, 0) = derff1.z;
+ JacobinF(2, 1) = -derff2.z;
+ JacobinF(2, 2) = -derff3.z;
+
+ //std::cout << "Jacobin Mat :\n" << JacobinF << '\n';
+ J_Inverse = JacobinF.inverse();
+ t_ = t - J_Inverse * Ft;
+ t = t_;
+ firstItr = false;
+ cnt++;
+ std::cout << "The " << cnt << "th iteration done!\n------------------------------------------------------\n";
+ if(cnt > 50)
+ {
+ break;
+ }
+
+ } while(iterFlag);
+
+
+ return Pt_curve;
}
\ No newline at end of file
diff --git a/src/show_libigl.cpp b/src/show_libigl.cpp
index b9b3c0c..861cc08 100644
--- a/src/show_libigl.cpp
+++ b/src/show_libigl.cpp
@@ -1,84 +1,108 @@
#include "bvh.hpp"
#include "show_libigl.hpp"
-void ShowCurve_Igl(tinynurbs::RationalCurve<double> &curve, double sampleNum, Eigen::RowVector3d color)
+void ShowCurve_Igl(tinynurbs::RationalCurve<double>& curve, double sampleNum, Eigen::RowVector3d color)
{
- double T = *(curve.knots.end() - 1);
+ 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, color);
- }
+ 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, color);
+ }
}
-void ShowSurface_Igl(tinynurbs::RationalSurface<double> &surface, double sampleNumU, double sampleNumV, Eigen::RowVector3d color)
+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);
+ 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, color);
- }
- }
+ /*
+ 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, color);
+ }
+ }
+ */
+ for(double u = 0; u <= U; u += U / 50)
+ {
+ for(double v = 0; v <= V; v += V / 100)
+ {
+ auto point = tinynurbs::surfacePoint(surface, u, v);
+ Eigen::RowVector3d pt;
+ pt << point.x, point.y, point.z;
+ viewer.data().add_points(pt, color);
+ }
+ }
+
+ for(double v = 0; v <= V; v += V / 50)
+ {
+ for(double u = 0; u <= U; u += U / 100)
+ {
+ auto point = tinynurbs::surfacePoint(surface, u, v);
+ Eigen::RowVector3d pt;
+ pt << point.x, point.y, point.z;
+ viewer.data().add_points(pt, color);
+ }
+ }
+ //viewer.data().point_size = 3;
}
void ShowAABB_Igl(AABB bound, Eigen::RowVector3d color)
{
- 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::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::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;
+ 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)),
- color);
- return;
+ 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)),
+ color);
+ return;
}
void ShowBVHNode_Igl(BVH_AABB_NodePtr bvhNode, Eigen::RowVector3d color)
{
- if (bvhNode == NULL)
- return;
+ if(bvhNode == NULL)
+ return;
- ShowAABB_Igl(bvhNode->bound, color);
+ ShowAABB_Igl(bvhNode->bound, color);
- for (auto it : bvhNode->childPtr)
- {
- ShowBVHNode_Igl(it, color);
- }
+ for(auto it : bvhNode->childPtr)
+ {
+ ShowBVHNode_Igl(it, color);
+ }
}