licaiyue42
/
NurbsIntersection
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

modify some files

main
forty-twoo 2 years ago
parent
16484bcfea
commit
9f0869094f
7 changed files with 283 additions and 223 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 424
      3rdparty/tinynurbs/core/basis.h
  2. 8
      README.md
  3. 19
      examples/example3_curve_curve_intersection/main.cpp
  4. 11
      include/newton.hpp
  5. 3
      include/show_libigl.hpp
  6. 9
      src/intersection.cpp
  7. 32
      src/newton.cpp

424
3rdparty/tinynurbs/core/basis.h
View File

@ -1,6 +1,6 @@
/**
* Low-level functions for evaluating B-spline basis functions and their derivatives
*
*
* Use of this source code is governed by a BSD-style license that can be found in
* the LICENSE file.
*/
@ -15,261 +15,263 @@
namespace tinynurbs
{
/**
* Find the span of the given parameter in the knot vector.
* @param[in] degree Degree of the curve.
* @param[in] knots Knot vector of the curve.
* @param[in] u Parameter value.
* @return Span index into the knot vector such that (span - 1) < u <= span
*/
template <typename T> int findSpan(unsigned int degree, const std::vector<T> &knots, T u)
{
// index of last control point
int n = static_cast<int>(knots.size()) - degree - 2;
assert(n >= 0);
/*
// For u that is equal to last knot value
if (util::close(u, knots[n + 1])) {
return n;
}
*/
// For values of u that lies outside the domain
if (u > (knots[n + 1] - std::numeric_limits<T>::epsilon()))
{
return n;
}
if (u < (knots[degree] + std::numeric_limits<T>::epsilon()))
/**
* Find the span of the given parameter in the knot vector.
* @param[in] degree Degree of the curve.
* @param[in] knots Knot vector of the curve.
* @param[in] u Parameter value.
* @return Span index into the knot vector such that (span - 1) < u <= span
*/
template <typename T>
int findSpan(unsigned int degree, const std::vector<T> &knots, T u)
{
return degree;
}
// Binary search
// TODO: Replace this with std::lower_bound
int low = degree;
int high = n + 1;
int mid = (int)std::floor((low + high) / 2.0);
while (u < knots[mid] || u >= knots[mid + 1])
{
if (u < knots[mid])
// index of last control point
int n = static_cast<int>(knots.size()) - degree - 2;
assert(n >= 0);
/*
// For u that is equal to last knot value
if (util::close(u, knots[n + 1])) {
return n;
}
*/
// For values of u that lies outside the domain
if (u > (knots[n + 1] - std::numeric_limits<T>::epsilon()))
{
high = mid;
return n;
}
else
if (u < (knots[degree] + std::numeric_limits<T>::epsilon()))
{
low = mid;
return degree;
}
mid = (int)std::floor((low + high) / 2.0);
}
return mid;
}
/**
* Compute a single B-spline basis function
* @param[in] i The ith basis function to compute.
* @param[in] deg Degree of the basis function.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @return The value of the ith basis function at u.
*/
template <typename T> T bsplineOneBasis(int i, unsigned int deg, const std::vector<T> &U, T u)
{
int m = static_cast<int>(U.size()) - 1;
// Special case
if ((i == 0 && close(u, U[0])) || (i == m - deg - 1 && close(u, U[m])))
{
return 1.0;
}
// Local property ensures that basis function is zero outside span
if (u < U[i] || u >= U[i + deg + 1])
{
return 0.0;
}
// Initialize zeroth-degree functions
std::vector<double> N;
N.resize(deg + 1);
for (int j = 0; j <= static_cast<int>(deg); j++)
{
N[j] = (u >= U[i + j] && u < U[i + j + 1]) ? 1.0 : 0.0;
}
// Compute triangular table
for (int k = 1; k <= static_cast<int>(deg); k++)
{
T saved = (util::close(N[0], 0.0)) ? 0.0 : ((u - U[i]) * N[0]) / (U[i + k] - U[i]);
for (int j = 0; j < static_cast<int>(deg) - k + 1; j++)
// Binary search
// TODO: Replace this with std::lower_bound
int low = degree;
int high = n + 1;
int mid = (int)std::floor((low + high) / 2.0);
while (u < knots[mid] || u >= knots[mid + 1])
{
T Uleft = U[i + j + 1];
T Uright = U[i + j + k + 1];
if (util::close(N[j + 1], 0.0))
if (u < knots[mid])
{
N[j] = saved;
saved = 0.0;
high = mid;
}
else
{
T temp = N[j + 1] / (Uright - Uleft);
N[j] = saved + (Uright - u) * temp;
saved = (u - Uleft) * temp;
low = mid;
}
mid = (int)std::floor((low + high) / 2.0);
}
return mid;
}
return N[0];
}
/**
* Compute all non-zero B-spline basis functions
* @param[in] deg Degree of the basis function.
* @param[in] span Index obtained from findSpan() corresponding the u and knots.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @return N Values of (deg+1) non-zero basis functions.
*/
template <typename T>
std::vector<T> bsplineBasis(unsigned int deg, int span, const std::vector<T> &knots, T u)
{
std::vector<T> N;
N.resize(deg + 1, T(0));
std::vector<T> left, right;
left.resize(deg + 1, static_cast<T>(0.0));
right.resize(deg + 1, static_cast<T>(0.0));
T saved = 0.0, temp = 0.0;
N[0] = 1.0;
for (int j = 1; j <= static_cast<int>(deg); j++)
/**
* Compute a single B-spline basis function
* @param[in] i The ith basis function to compute.
* @param[in] deg Degree of the basis function.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @return The value of the ith basis function at u.
*/
template <typename T>
T bsplineOneBasis(int i, unsigned int deg, const std::vector<T> &U, T u)
{
left[j] = (u - knots[span + 1 - j]);
right[j] = knots[span + j] - u;
saved = 0.0;
for (int r = 0; r < j; r++)
int m = static_cast<int>(U.size()) - 1;
// Special case
if ((i == 0 && close(u, U[0])) || (i == m - deg - 1 && close(u, U[m])))
{
return 1.0;
}
// Local property ensures that basis function is zero outside span
if (u < U[i] || u >= U[i + deg + 1])
{
return 0.0;
}
// Initialize zeroth-degree functions
std::vector<double> N;
N.resize(deg + 1);
for (int j = 0; j <= static_cast<int>(deg); j++)
{
temp = N[r] / (right[r + 1] + left[j - r]);
N[r] = saved + right[r + 1] * temp;
saved = left[j - r] * temp;
N[j] = (u >= U[i + j] && u < U[i + j + 1]) ? 1.0 : 0.0;
}
N[j] = saved;
// Compute triangular table
for (int k = 1; k <= static_cast<int>(deg); k++)
{
T saved = (util::close(N[0], 0.0)) ? 0.0 : ((u - U[i]) * N[0]) / (U[i + k] - U[i]);
for (int j = 0; j < static_cast<int>(deg) - k + 1; j++)
{
T Uleft = U[i + j + 1];
T Uright = U[i + j + k + 1];
if (util::close(N[j + 1], 0.0))
{
N[j] = saved;
saved = 0.0;
}
else
{
T temp = N[j + 1] / (Uright - Uleft);
N[j] = saved + (Uright - u) * temp;
saved = (u - Uleft) * temp;
}
}
}
return N[0];
}
return N;
}
/**
* Compute all non-zero derivatives of B-spline basis functions
* @param[in] deg Degree of the basis function.
* @param[in] span Index obtained from findSpan() corresponding the u and knots.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @param[in] num_ders Number of derivatives to compute (num_ders <= deg)
* @return ders Values of non-zero derivatives of basis functions.
*/
template <typename T>
array2<T> bsplineDerBasis(unsigned int deg, int span, const std::vector<T> &knots, T u,
int num_ders)
{
std::vector<T> left, right;
left.resize(deg + 1, 0.0);
right.resize(deg + 1, 0.0);
T saved = 0.0, temp = 0.0;
array2<T> ndu(deg + 1, deg + 1);
ndu(0, 0) = 1.0;
for (int j = 1; j <= static_cast<int>(deg); j++)
/**
* Compute all non-zero B-spline basis functions
* @param[in] deg Degree of the basis function.
* @param[in] span Index obtained from findSpan() corresponding the u and knots.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @return N Values of (deg+1) non-zero basis functions.
*/
template <typename T>
std::vector<T> bsplineBasis(unsigned int deg, int span, const std::vector<T> &knots, T u)
{
left[j] = u - knots[span + 1 - j];
right[j] = knots[span + j] - u;
saved = 0.0;
std::vector<T> N;
N.resize(deg + 1, T(0));
std::vector<T> left, right;
left.resize(deg + 1, static_cast<T>(0.0));
right.resize(deg + 1, static_cast<T>(0.0));
T saved = 0.0, temp = 0.0;
N[0] = 1.0;
for (int r = 0; r < j; r++)
for (int j = 1; j <= static_cast<int>(deg); j++)
{
// Lower triangle
ndu(j, r) = right[r + 1] + left[j - r];
temp = ndu(r, j - 1) / ndu(j, r);
// Upper triangle
ndu(r, j) = saved + right[r + 1] * temp;
saved = left[j - r] * temp;
left[j] = (u - knots[span + 1 - j]);
right[j] = knots[span + j] - u;
saved = 0.0;
for (int r = 0; r < j; r++)
{
temp = N[r] / (right[r + 1] + left[j - r]);
N[r] = saved + right[r + 1] * temp;
saved = left[j - r] * temp;
}
N[j] = saved;
}
ndu(j, j) = saved;
return N;
}
array2<T> ders(num_ders + 1, deg + 1, T(0));
for (int j = 0; j <= static_cast<int>(deg); j++)
/**
* Compute all non-zero derivatives of B-spline basis functions
* @param[in] deg Degree of the basis function.
* @param[in] span Index obtained from findSpan() corresponding the u and knots.
* @param[in] knots Knot vector corresponding to the basis functions.
* @param[in] u Parameter to evaluate the basis functions at.
* @param[in] num_ders Number of derivatives to compute (num_ders <= deg)
* @return ders Values of non-zero derivatives of basis functions.
*/
template <typename T>
array2<T> bsplineDerBasis(unsigned int deg, int span, const std::vector<T> &knots, T u,
int num_ders)
{
ders(0, j) = ndu(j, deg);
}
std::vector<T> left, right;
left.resize(deg + 1, 0.0);
right.resize(deg + 1, 0.0);
T saved = 0.0, temp = 0.0;
array2<T> a(2, deg + 1);
array2<T> ndu(deg + 1, deg + 1);
ndu(0, 0) = 1.0;
for (int r = 0; r <= static_cast<int>(deg); r++)
{
int s1 = 0;
int s2 = 1;
a(0, 0) = 1.0;
for (int k = 1; k <= num_ders; k++)
for (int j = 1; j <= static_cast<int>(deg); j++)
{
T d = 0.0;
int rk = r - k;
int pk = deg - k;
int j1 = 0;
int j2 = 0;
left[j] = u - knots[span + 1 - j];
right[j] = knots[span + j] - u;
saved = 0.0;
if (r >= k)
for (int r = 0; r < j; r++)
{
a(s2, 0) = a(s1, 0) / ndu(pk + 1, rk);
d = a(s2, 0) * ndu(rk, pk);
// Lower triangle
ndu(j, r) = right[r + 1] + left[j - r];
temp = ndu(r, j - 1) / ndu(j, r);
// Upper triangle
ndu(r, j) = saved + right[r + 1] * temp;
saved = left[j - r] * temp;
}
if (rk >= -1)
{
j1 = 1;
}
else
{
j1 = -rk;
}
ndu(j, j) = saved;
}
if (r - 1 <= pk)
{
j2 = k - 1;
}
else
{
j2 = deg - r;
}
array2<T> ders(num_ders + 1, deg + 1, T(0));
for (int j = j1; j <= j2; j++)
{
a(s2, j) = (a(s1, j) - a(s1, j - 1)) / ndu(pk + 1, rk + j);
d += a(s2, j) * ndu(rk + j, pk);
}
for (int j = 0; j <= static_cast<int>(deg); j++)
{
ders(0, j) = ndu(j, deg);
}
array2<T> a(2, deg + 1);
for (int r = 0; r <= static_cast<int>(deg); r++)
{
int s1 = 0;
int s2 = 1;
a(0, 0) = 1.0;
if (r <= pk)
for (int k = 1; k <= num_ders; k++)
{
a(s2, k) = -a(s1, k - 1) / ndu(pk + 1, r);
d += a(s2, k) * ndu(r, pk);
}
T d = 0.0;
int rk = r - k;
int pk = deg - k;
int j1 = 0;
int j2 = 0;
if (r >= k)
{
a(s2, 0) = a(s1, 0) / ndu(pk + 1, rk);
d = a(s2, 0) * ndu(rk, pk);
}
if (rk >= -1)
{
j1 = 1;
}
else
{
j1 = -rk;
}
if (r - 1 <= pk)
{
j2 = k - 1;
}
else
{
j2 = deg - r;
}
for (int j = j1; j <= j2; j++)
{
a(s2, j) = (a(s1, j) - a(s1, j - 1)) / ndu(pk + 1, rk + j);
d += a(s2, j) * ndu(rk + j, pk);
}
ders(k, r) = d;
if (r <= pk)
{
a(s2, k) = -a(s1, k - 1) / ndu(pk + 1, r);
d += a(s2, k) * ndu(r, pk);
}
int temp = s1;
s1 = s2;
s2 = temp;
ders(k, r) = d;
int temp = s1;
s1 = s2;
s2 = temp;
}
}
}
T fac = static_cast<T>(deg);
for (int k = 1; k <= num_ders; k++)
{
for (int j = 0; j <= static_cast<int>(deg); j++)
T fac = static_cast<T>(deg);
for (int k = 1; k <= num_ders; k++)
{
ders(k, j) *= fac;
for (int j = 0; j <= static_cast<int>(deg); j++)
{
ders(k, j) *= fac;
}
fac *= static_cast<T>(deg - k);
}
fac *= static_cast<T>(deg - k);
}
return ders;
}
return ders;
}
} // namespace tinynurbs

8
README.md
View File

@ -4,12 +4,18 @@
**STL**, **Libigl**, **Eigen** and the dependencies of the **igl::opengl::glfw::Viewer (OpenGL, glad and GLFW)**. The CMake build system will automatically download libigl and its dependencies using CMake FetchContent, thus requiring no setup on your part.
## Platform
The project could run both on Linux and Windows.
## Compile
```
mkdir build
cd build
cmake ..
make
make #if you are using visual studio, find the .sln project and build it.
```
## Run

19
examples/example3_curve_curve_intersection/main.cpp
View File

@ -28,6 +28,25 @@ int main(int argc, char *argv[])
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);

11
include/newton.hpp
View File

@ -1,4 +1,15 @@
/*
* @File Created: 2022-11-26, 17:43:27
* @Last Modified: 2022-11-26, 18:50:24
* @Author: forty-twoo
* @Copyright (c) 2022, Caiyue Li(li_caiyue@zju.edu.cn), All rights reserved.
*/
#ifndef NEWTON_H_
#define NEWTON_H_
#include "bvh.hpp"
std::vector<double> newton_curve_curve(BVH_AABB_NodePtr boxPtr1, BVH_AABB_NodePtr boxPtr2, const double &eps, bool &isConvergence);
#endif

3
include/show_libigl.hpp
View File

@ -21,7 +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);

9
src/intersection.cpp
View File

@ -13,15 +13,6 @@ bool CurveCurveBVHIntersect(BVH_AABB_NodePtr &BoxPtr1, BVH_AABB_NodePtr &BoxPtr2
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 = true;

32
src/newton.cpp
View File

@ -1 +1,31 @@
#include "newton.hpp"
/*
* @File Created: 2022-11-26, 17:43:34
* @Last Modified: 2022-11-26, 18:50:31
* @Author: forty-twoo
* @Copyright (c) 2022, Caiyue Li(li_caiyue@zju.edu.cn), All rights reserved.
*/
#include "newton.hpp"
#include <Eigen/Core>
std::vector<double> 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 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;
Eigen::MatrixXd F(2, 2);
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::MatrixXd JacobinF(2, 2);
bool iterFlag = true;
}
Write Preview
Loading…
Cancel
Save