MercuryDPM  Trunk
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
NurbsUtils.cc
Go to the documentation of this file.
1 //Copyright (c) 2013-2020, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 #include <algorithm>
27 #include "NurbsUtils.h"
28 #include "Logger.h"
29 
30 namespace NurbsUtils
31 {
32 
33 bool isKnotVectorMonotonic(const std::vector<double>& knots)
34 {
35  return std::is_sorted(knots.begin(), knots.end());
36 }
37 
38 bool close(double a, double b, double eps)
39 {
40  return (std::abs(a - b) < eps) ? true : false;
41 }
42 
43 int findSpan(int degree, const std::vector<double>& knots, double u)
44 {
45  // index of last control point
46  int n = static_cast<int>(knots.size()) - degree - 2;
47 
48  // For u that is equal to last knot value
49  if (close(u, knots[n + 1]))
50  {
51  return n;
52  }
53 
54  // For values of u that lies outside the domain
55  if (u > knots[n + 1])
56  {
57  return n;
58  }
59  if (u < knots[degree])
60  {
61  return degree;
62  }
63 
64  // Binary search
65  // TODO: Replace this with std::lower_bound
66  int low = degree;
67  int high = n + 1;
68  int mid = (int) std::floor((low + high) / 2.0);
69  while (u < knots[mid] || u >= knots[mid + 1])
70  {
71  if (u < knots[mid])
72  {
73  high = mid;
74  }
75  else
76  {
77  low = mid;
78  }
79  mid = (low + high) / 2;
80  }
81  return mid;
82 }
83 
84 double bsplineOneBasis(int i, int deg, const std::vector<double>& U, double u)
85 {
86  int m = static_cast<int>(U.size()) - 1;
87  // Special case
88  if ((i == 0 && close(u, U[0])) || (i == m - deg - 1 && close(u, U[m])))
89  {
90  return 1.0;
91  }
92  // Local Property
93  if (u < U[i] || u >= U[i + deg + 1])
94  {
95  return 0.0;
96  }
97  // Initialize zeroth-degree functions
98  std::vector<double> N;
99  N.resize(deg + 1);
100  for (int j = 0; j <= deg; j++)
101  {
102  N[j] = (u >= U[i + j] && u < U[i + j + 1]) ? 1.0 : 0.0;
103  }
104  // Compute triangular table
105  for (int k = 1; k <= deg; k++)
106  {
107  double saved = (close(N[0], 0.0)) ? 0.0
108  : ((u - U[i]) * N[0]) / (U[i + k] - U[i]);
109  for (int j = 0; j < deg - k + 1; j++)
110  {
111  double Uleft = U[i + j + 1];
112  double Uright = U[i + j + k + 1];
113  if (close(N[j + 1], 0.0))
114  {
115  N[j] = saved;
116  saved = 0.0;
117  }
118  else
119  {
120  double temp = N[j + 1] / (Uright - Uleft);
121  N[j] = saved + (Uright - u) * temp;
122  saved = (u - Uleft) * temp;
123  }
124  }
125  }
126  return N[0];
127 }
128 
129 void bsplineBasis(int deg, int span, const std::vector<double>& knots, double u,
130  std::vector<double>& N)
131 {
132  N.clear();
133  N.resize(deg + 1, 0.0);
134  std::vector<double> left, right;
135  left.resize(deg + 1, 0.0);
136  right.resize(deg + 1, 0.0);
137 
138  N[0] = 1.0;
139 
140  for (int j = 1; j <= deg; j++)
141  {
142  left[j] = (u - knots[span + 1 - j]);
143  right[j] = knots[span + j] - u;
144  Mdouble saved = 0.0;
145  for (int r = 0; r < j; r++)
146  {
147  const Mdouble temp = N[r] / (right[r + 1] + left[j - r]);
148  N[r] = saved + right[r + 1] * temp;
149  saved = left[j - r] * temp;
150  }
151  N[j] = saved;
152  }
153 }
154 
155 void bsplineDerBasis(int deg, int span, const std::vector<double>& knots, double u,
156  int nDers, std::vector<std::vector<double>> &ders) {
157 
158  std::vector<double> left, right;
159  left.resize(deg + 1, 0.0);
160  right.resize(deg + 1, 0.0);
161  double saved = 0.0, temp = 0.0;
162 
163  array2<double> ndu(deg + 1, deg + 1);
164  ndu(0, 0) = 1.0;
165 
166  for (int j = 1; j <= deg; j++) {
167  left[j] = u - knots[span + 1 - j];
168  right[j] = knots[span + j] - u;
169  saved = 0.0;
170 
171  for (int r = 0; r < j; r++) {
172  // Lower triangle
173  ndu(j, r) = right[r + 1] + left[j - r];
174  temp = ndu(r, j - 1) / ndu(j, r);
175  // Upper triangle
176  ndu(r, j) = saved + right[r + 1] * temp;
177  saved = left[j - r] * temp;
178  }
179 
180  ndu(j, j) = saved;
181  }
182 
183  ders.clear();
184  ders.resize(nDers + 1);
185  for (int i = 0; i < ders.size(); i++) {
186  ders[i].resize(deg + 1, 0.0);
187  }
188 
189  for (int j = 0; j <= deg; j++) {
190  ders[0][j] = ndu(j, deg);
191  }
192 
193  array2<double> a(2, deg + 1);
194 
195  for (int r = 0; r <= deg; r++) {
196  int s1 = 0;
197  int s2 = 1;
198  a(0, 0) = 1.0;
199 
200  for (int k = 1; k <= nDers; k++) {
201  double d = 0.0;
202  int rk = r - k;
203  int pk = deg - k;
204  int j1 = 0;
205  int j2 = 0;
206 
207  if (r >= k) {
208  a(s2, 0) = a(s1, 0) / ndu(pk + 1, rk);
209  d = a(s2, 0) * ndu(rk, pk);
210  }
211 
212  if (rk >= -1) {
213  j1 = 1;
214  }
215  else {
216  j1 = -rk;
217  }
218 
219  if (r - 1 <= pk) {
220  j2 = k - 1;
221  }
222  else {
223  j2 = deg - r;
224  }
225 
226  for (int j = j1; j <= j2; j++) {
227  a(s2, j) = (a(s1, j) - a(s1, j - 1)) / ndu(pk + 1, rk + j);
228  d += a(s2, j) * ndu(rk + j, pk);
229  }
230 
231  if (r <= pk) {
232  a(s2, k) = -a(s1, k - 1) / ndu(pk + 1, r);
233  d += a(s2, k) * ndu(r, pk);
234  }
235 
236 
237  ders[k][r] = d;
238 
239  int temp = s1;
240  s1 = s2;
241  s2 = temp;
242  }
243  }
244 
245  double fac = static_cast<double>(deg);
246  for (int k = 1; k <= nDers; k++) {
247  for (int j = 0; j <= deg; j++) {
248  ders[k][j] *= fac;
249  }
250  fac *= static_cast<double>(deg - k);
251  }
252 }
253 
254 
255 }
int findSpan(int degree, const std::vector< double > &knots, double u)
Find the span of the given parameter in the knot vector.
Definition: NurbsUtils.cc:43
double Mdouble
Definition: GeneralDefine.h:34
const std::complex< Mdouble > i
Definition: ExtendedMath.h:50
double bsplineOneBasis(int i, int deg, const std::vector< double > &U, double u)
Compute a single B-spline basis function.
Definition: NurbsUtils.cc:84
A simple class for representing 2D runtime arrays.
Definition: NurbsUtils.h:41
bool isKnotVectorMonotonic(const std::vector< double > &knots)
Definition: NurbsUtils.cc:33
bool close(double a, double b, double eps)
Definition: NurbsUtils.cc:38
void bsplineDerBasis(int deg, int span, const std::vector< double > &knots, double u, int nDers, std::vector< std::vector< double >> &ders)
Compute all non-zero derivatives of B-spline basis functions.
Definition: NurbsUtils.cc:155
void bsplineBasis(int deg, int span, const std::vector< double > &knots, double u, std::vector< double > &N)
Compute all non-zero B-spline basis functions.
Definition: NurbsUtils.cc:129