add curve surface intersection code

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)

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

examples/example3_curve_curve_intersection/main.cpp

@@ -10,6 +10,7 @@
 #include "bvh.hpp"
 #include "show_libigl.hpp"
 #include "intersection.hpp"
+#include "newton.hpp"
 
 igl::opengl::glfw::Viewer viewer;
 
@@ -23,53 +24,43 @@ int main(int argc, char *argv[])
 		glm::vec3(-1, 4, 0),
 		glm::vec3(1, -5, 0),
 		glm::vec3(2, 1, 0),
-                           glm::vec3(4, -1, 0)};
+		glm::vec3(4, -1, 0),
-    crv1.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+		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};
+	crv1.weights = { 1, 1, 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);
+	ShowCurve_Igl(crv1, 500, YELLOW);
-    for (auto it : basis)
+	BVH_AABB bvh_curve1(12, 2, 100, crv1);
-    {
-        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(-5, 2, 0),
-                           glm::vec3(1, -1, 0),
+		glm::vec3(-4, 6, 0),
-                           glm::vec3(3, -3, 0),
+		glm::vec3(-2, -5, 0),
-                           glm::vec3(4, -4, 0)};
+		glm::vec3(1, 8, 0),
-    crv2.knots = {0, 0, 0, 1, 2, 3, 3, 3}; // std::vector of floats
+		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};
+	crv2.weights = { 1, 1, 1, 1, 1, 1 };
 
-    ShowCurve_Igl(crv2, 5000, YELLOW);
+	ShowCurve_Igl(crv2, 500, YELLOW);
-    BVH_AABB bvh_curve2(11, 2, 100, crv2);
+	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))
 	{
@@ -78,6 +69,7 @@ int main(int argc, char *argv[])
 		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);
@@ -88,7 +80,36 @@ int main(int argc, char *argv[])
 				<< 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;
+		}
+
+
+
+
 	}
 
 	/*

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)

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;
+
+
+}

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

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

src/bvh.cpp

@@ -60,16 +60,17 @@ BVH_AABB_NodePtr BVH_AABB::Build_NurbsCurve(tinynurbs::RationalCurve<double> &cu
 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)
+	if(curLevel + 1 == Hlevel)
 	{
-        for (double i = u; i < u + ustep; i++)
+		for(double i = u; i < u + ustep; i += ustep / sampleNum)
 		{
-            for (double j = v; j < v + vstep; j++)
+			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));
 
@@ -80,16 +81,27 @@ BVH_AABB_NodePtr BVH_AABB::Build_NurbsSurface(tinynurbs::RationalSurface<double>
 	BVH_AABB_NodePtr childNode3 = new BVH_AABB_Node;
 	BVH_AABB_NodePtr childNode4 = new BVH_AABB_Node;
 
-    double ustep_ = ustep / (double)BranchNum;
+	double ustep_ = ustep / (double)2;
-    double vstep_ = vstep / (double)BranchNum;
+	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;
 }

src/intersection.cpp

@@ -8,7 +8,7 @@
 #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;
@@ -20,8 +20,7 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
 		{
 			if(x1_ > x2)
 				LIntersect = false;
-        }
+		} else
-        else
 		{
 			if(x1 > x2_)
 				LIntersect = false;
@@ -60,25 +59,23 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
 			}
 			for(auto it : BoxPtr1->childPtr)
 			{
-                CurveCurveBVHIntersect(it, BoxPtr2, IstNodePtr);
+				BVHIntersect(it, BoxPtr2, IstNodePtr);
 			}
-        }
+		} else
-        else
 		{
 			if(BoxPtr2->childPtr.size() == 0)
 			{
-                std::pair<BVH_AABB_NodePtr, BVH_AABB_NodePtr> pi = std::make_pair(BoxPtr2, BoxPtr1);
+				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)
 			{
-                CurveCurveBVHIntersect(BoxPtr1, it, IstNodePtr);
+				BVHIntersect(BoxPtr1, it, IstNodePtr);
 			}
 		}
 		return true;
-    }
+	} else
-    else
 		return false;
 }

src/newton.cpp

@@ -6,9 +6,11 @@
  */
 
 #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;
@@ -17,15 +19,261 @@ std::vector<double> newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePt
 	double param2_ed = boxPtr2->param[0].second;
 
 	double t1 = (param1_st + param1_ed) / 2.0;
-    double t2 = (param2_st + param2_st) / 2.0;
 
-    Eigen::MatrixXd F(2, 2);
+	double t2 = (param2_st + param2_ed) / 2.0;
+
+	Eigen::Matrix2d F(2, 2);
+	Eigen::Vector2d t, t_, Ft;
+	t << t1, t2;
+
+	tinynurbs::RationalCurve<double>* crvPtr1 = boxPtr1->myBVH->NurbsCurvePtr;
+	tinynurbs::RationalCurve<double>* crvPtr2 = boxPtr2->myBVH->NurbsCurvePtr;
+
+	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::array2<double> der_Fx = tinynurbs::bsplineDerBasis(crvPtr1->degree, tinynurbs::findSpan(crvPtr1->degree, crvPtr1->knots, t1), crvPtr1->knots, t1, 1);
+	tinynurbs::RationalSurface<double>* srfPtr2 = boxPtr2->myBVH->NurbsSurfacePtr;
-    Eigen::Vector2d t(t1, t2);
+
+	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);
 
-    Eigen::MatrixXd JacobinF(2, 2);
 
-    bool iterFlag = true;
+	return Pt_curve;
 }

src/show_libigl.cpp

@@ -18,6 +18,7 @@ void ShowSurface_Igl(tinynurbs::RationalSurface<double> &surface, double sampleN
 	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)
@@ -28,6 +29,29 @@ void ShowSurface_Igl(tinynurbs::RationalSurface<double> &surface, double sampleN
 			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)