Visual Servoing Platform  version 3.6.1 under development (2024-05-21)
vpNurbs.cpp
1 /*
2  * ViSP, open source Visual Servoing Platform software.
3  * Copyright (C) 2005 - 2023 by Inria. All rights reserved.
4  *
5  * This software is free software; you can redistribute it and/or modify
6  * it under the terms of the GNU General Public License as published by
7  * the Free Software Foundation; either version 2 of the License, or
8  * (at your option) any later version.
9  * See the file LICENSE.txt at the root directory of this source
10  * distribution for additional information about the GNU GPL.
11  *
12  * For using ViSP with software that can not be combined with the GNU
13  * GPL, please contact Inria about acquiring a ViSP Professional
14  * Edition License.
15  *
16  * See https://visp.inria.fr for more information.
17  *
18  * This software was developed at:
19  * Inria Rennes - Bretagne Atlantique
20  * Campus Universitaire de Beaulieu
21  * 35042 Rennes Cedex
22  * France
23  *
24  * If you have questions regarding the use of this file, please contact
25  * Inria at visp@inria.fr
26  *
27  * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
28  * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
29  *
30  * Description:
31  * This class implements the Non Uniform Rational B-Spline (NURBS)
32  */
33 
34 #include <cmath> // std::fabs
35 #include <limits> // numeric_limits
36 #include <visp3/core/vpColVector.h>
37 #include <visp3/me/vpNurbs.h>
38 /*
39  Compute the distance d = |Pw1-Pw2|
40 */
41 inline double distance(const vpImagePoint &iP1, double w1, const vpImagePoint &iP2, double w2)
42 {
43  double distancei = iP1.get_i() - iP2.get_i();
44  double distancej = iP1.get_j() - iP2.get_j();
45  double distancew = w1 - w2;
46  return sqrt(vpMath::sqr(distancei) + vpMath::sqr(distancej) + vpMath::sqr(distancew));
47 }
48 
49 vpNurbs::vpNurbs() : weights() { p = 3; }
50 
51 vpNurbs::vpNurbs(const vpNurbs &nurbs) : vpBSpline(nurbs), weights(nurbs.weights) { }
52 
53 vpImagePoint vpNurbs::computeCurvePoint(double l_u, unsigned int l_i, unsigned int l_p, std::vector<double> &l_knots,
54  std::vector<vpImagePoint> &l_controlPoints, std::vector<double> &l_weights)
55 {
56  vpBasisFunction *N = nullptr;
57  N = computeBasisFuns(l_u, l_i, l_p, l_knots);
58  vpImagePoint pt;
59 
60  double ic = 0;
61  double jc = 0;
62  double wc = 0;
63  for (unsigned int j = 0; j <= l_p; j++) {
64  ic = ic + N[j].value * (l_controlPoints[l_i - l_p + j]).get_i() * l_weights[l_i - l_p + j];
65  jc = jc + N[j].value * (l_controlPoints[l_i - l_p + j]).get_j() * l_weights[l_i - l_p + j];
66  wc = wc + N[j].value * l_weights[l_i - l_p + j];
67  }
68 
69  pt.set_i(ic / wc);
70  pt.set_j(jc / wc);
71 
72  if (N != nullptr)
73  delete[] N;
74 
75  return pt;
76 }
77 
79 {
80  vpBasisFunction *N = nullptr;
81  N = computeBasisFuns(u);
82  vpImagePoint pt;
83 
84  double ic = 0;
85  double jc = 0;
86  double wc = 0;
87  for (unsigned int j = 0; j <= p; j++) {
88  ic = ic + N[j].value * (controlPoints[N[0].i + j]).get_i() * weights[N[0].i + j]; // N[0].i = findSpan(u)-p
89  jc = jc + N[j].value * (controlPoints[N[0].i + j]).get_j() * weights[N[0].i + j];
90  wc = wc + N[j].value * weights[N[0].i + j];
91  }
92 
93  pt.set_i(ic / wc);
94  pt.set_j(jc / wc);
95 
96  if (N != nullptr)
97  delete[] N;
98 
99  return pt;
100 }
101 
102 vpMatrix vpNurbs::computeCurveDers(double l_u, unsigned int l_i, unsigned int l_p, unsigned int l_der,
103  std::vector<double> &l_knots, std::vector<vpImagePoint> &l_controlPoints,
104  std::vector<double> &l_weights)
105 {
106  vpMatrix derivate(l_der + 1, 3);
107  vpBasisFunction **N = nullptr;
108  N = computeDersBasisFuns(l_u, l_i, l_p, l_der, l_knots);
109 
110  for (unsigned int k = 0; k <= l_der; k++) {
111  derivate[k][0] = 0.0;
112  derivate[k][1] = 0.0;
113  derivate[k][2] = 0.0;
114 
115  for (unsigned int j = 0; j <= l_p; j++) {
116  derivate[k][0] = derivate[k][0] + N[k][j].value * (l_controlPoints[l_i - l_p + j]).get_i();
117  derivate[k][1] = derivate[k][1] + N[k][j].value * (l_controlPoints[l_i - l_p + j]).get_j();
118  derivate[k][2] = derivate[k][2] + N[k][j].value * (l_weights[l_i - l_p + j]);
119  }
120  }
121 
122  if (N != nullptr) {
123  for (unsigned int i = 0; i <= l_der; i++)
124  delete[] N[i];
125  delete[] N;
126  }
127  return derivate;
128 }
129 
130 vpMatrix vpNurbs::computeCurveDers(double u, unsigned int der)
131 {
132  vpMatrix derivate(der + 1, 3);
133  vpBasisFunction **N = nullptr;
134  N = computeDersBasisFuns(u, der);
135 
136  for (unsigned int k = 0; k <= der; k++) {
137  derivate[k][0] = 0.0;
138  derivate[k][1] = 0.0;
139  derivate[k][2] = 0.0;
140  for (unsigned int j = 0; j <= p; j++) {
141  derivate[k][0] = derivate[k][0] + N[k][j].value * (controlPoints[N[0][0].i - p + j]).get_i();
142  derivate[k][1] = derivate[k][1] + N[k][j].value * (controlPoints[N[0][0].i - p + j]).get_j();
143  derivate[k][2] = derivate[k][2] + N[k][j].value * (weights[N[0][0].i - p + j]);
144  }
145  }
146 
147  if (N != nullptr)
148  delete[] N;
149 
150  return derivate;
151 }
152 
153 vpImagePoint *vpNurbs::computeCurveDersPoint(double l_u, unsigned int l_i, unsigned int l_p, unsigned int l_der,
154  std::vector<double> &l_knots, std::vector<vpImagePoint> &l_controlPoints,
155  std::vector<double> &l_weights)
156 {
157  std::vector<vpImagePoint> A;
158  vpImagePoint pt;
159  for (unsigned int j = 0; j < l_controlPoints.size(); j++) {
160  pt = l_controlPoints[j];
161  pt.set_i(pt.get_i() * l_weights[j]);
162  pt.set_j(pt.get_j() * l_weights[j]);
163  A.push_back(pt);
164  }
165 
166  vpMatrix Awders = computeCurveDers(l_u, l_i, l_p, l_der, l_knots, A, l_weights);
167 
168  vpImagePoint *CK = new vpImagePoint[l_der + 1];
169 
170  for (unsigned int k = 0; k <= l_der; k++) {
171  double ic = Awders[k][0];
172  double jc = Awders[k][1];
173  for (unsigned int j = 1; j <= k; j++) {
174  double tmpComb = static_cast<double>(vpMath::comb(k, j));
175  ic = ic - tmpComb * Awders[k][2] * (CK[k - j].get_i());
176  jc = jc - tmpComb * Awders[j][2] * (CK[k - j].get_j());
177  }
178  CK[k].set_ij(ic / Awders[0][2], jc / Awders[0][2]);
179  }
180  return CK;
181 }
182 
183 
184 vpImagePoint *vpNurbs::computeCurveDersPoint(double u, unsigned int der)
185 {
186  unsigned int i = findSpan(u);
187  return computeCurveDersPoint(u, i, p, der, knots, controlPoints, weights);
188 }
189 
190 
191 void vpNurbs::curveKnotIns(double l_u, unsigned int l_k, unsigned int l_s, unsigned int l_r, unsigned int l_p,
192  std::vector<double> &l_knots, std::vector<vpImagePoint> &l_controlPoints,
193  std::vector<double> &l_weights)
194 {
195  vpMatrix Rw(l_p + 1, 3);
196  std::vector<vpImagePoint>::iterator it1;
197  std::vector<double>::iterator it2;
198  vpImagePoint pt;
199  double w = 0;
200 
201  for (unsigned int j = 0; j <= l_p - l_s; j++) {
202  Rw[j][0] = (l_controlPoints[l_k - l_p + j]).get_i() * l_weights[l_k - l_p + j];
203  Rw[j][1] = (l_controlPoints[l_k - l_p + j]).get_j() * l_weights[l_k - l_p + j];
204  Rw[j][2] = l_weights[l_k - l_p + j];
205  }
206 
207  it1 = l_controlPoints.begin();
208  l_controlPoints.insert(it1 + (int)l_k - (int)l_s, l_r, pt);
209  it2 = l_weights.begin();
210  l_weights.insert(it2 + (int)l_k - (int)l_s, l_r, w);
211 
212  unsigned int L = 0;
213  double alpha;
214  for (unsigned int j = 1; j <= l_r; j++) {
215  L = l_k - l_p + j;
216 
217  for (unsigned int i = 0; i <= l_p - j - l_s; i++) {
218  alpha = (l_u - l_knots[L + i]) / (l_knots[i + l_k + 1] - l_knots[L + i]);
219  Rw[i][0] = alpha * Rw[i + 1][0] + (1.0 - alpha) * Rw[i][0];
220  Rw[i][1] = alpha * Rw[i + 1][1] + (1.0 - alpha) * Rw[i][1];
221  Rw[i][2] = alpha * Rw[i + 1][2] + (1.0 - alpha) * Rw[i][2];
222  }
223 
224  pt.set_ij(Rw[0][0] / Rw[0][2], Rw[0][1] / Rw[0][2]);
225  l_controlPoints[L] = pt;
226  l_weights[L] = Rw[0][2];
227 
228  pt.set_ij(Rw[l_p - j - l_s][0] / Rw[l_p - j - l_s][2], Rw[l_p - j - l_s][1] / Rw[l_p - j - l_s][2]);
229  l_controlPoints[l_k + l_r - j - l_s] = pt;
230  l_weights[l_k + l_r - j - l_s] = Rw[l_p - j - l_s][2];
231  }
232 
233  for (unsigned int j = L + 1; j < l_k - l_s; j++) {
234  pt.set_ij(Rw[j - L][0] / Rw[j - L][2], Rw[j - L][1] / Rw[j - L][2]);
235  l_controlPoints[j] = pt;
236  l_weights[j] = Rw[j - L][2];
237  }
238 
239  it2 = l_knots.begin();
240  l_knots.insert(it2 + (int)l_k, l_r, l_u);
241 }
242 
243 void vpNurbs::curveKnotIns(double u, unsigned int s, unsigned int r)
244 {
245  unsigned int i = findSpan(u);
246  curveKnotIns(u, i, s, r, p, knots, controlPoints, weights);
247 }
248 
249 
250 void vpNurbs::refineKnotVectCurve(double *l_x, unsigned int l_r, unsigned int l_p, std::vector<double> &l_knots,
251  std::vector<vpImagePoint> &l_controlPoints, std::vector<double> &l_weights)
252 {
253  unsigned int a = findSpan(l_x[0], l_p, l_knots);
254  unsigned int b = findSpan(l_x[l_r], l_p, l_knots);
255  b++;
256 
257  unsigned int n = (unsigned int)l_controlPoints.size();
258  unsigned int m = (unsigned int)l_knots.size();
259 
260  for (unsigned int j = 0; j < n; j++) {
261  l_controlPoints[j].set_ij(l_controlPoints[j].get_i() * l_weights[j], l_controlPoints[j].get_j() * l_weights[j]);
262  }
263 
264  std::vector<double> l_knots_tmp(l_knots);
265  std::vector<vpImagePoint> l_controlPoints_tmp(l_controlPoints);
266  std::vector<double> l_weights_tmp(l_weights);
267 
268  vpImagePoint pt;
269  double w = 0;
270 
271  for (unsigned int j = 0; j <= l_r; j++) {
272  l_controlPoints.push_back(pt);
273  l_weights.push_back(w);
274  l_knots.push_back(w);
275  }
276 
277  for (unsigned int j = b + l_p; j <= m - 1; j++)
278  l_knots[j + l_r + 1] = l_knots_tmp[j];
279 
280  for (unsigned int j = b - 1; j <= n - 1; j++) {
281  l_controlPoints[j + l_r + 1] = l_controlPoints_tmp[j];
282  l_weights[j + l_r + 1] = l_weights_tmp[j];
283  }
284 
285  unsigned int i = b + l_p - 1;
286  unsigned int k = b + l_p + l_r;
287 
288  {
289  unsigned int j = l_r + 1;
290  do {
291  j--;
292  while (l_x[j] <= l_knots[i] && i > a) {
293  l_controlPoints[k - l_p - 1] = l_controlPoints_tmp[i - l_p - 1];
294  l_weights[k - l_p - 1] = l_weights_tmp[i - l_p - 1];
295  l_knots[k] = l_knots_tmp[i];
296  k--;
297  i--;
298  }
299 
300  l_controlPoints[k - l_p - 1] = l_controlPoints[k - l_p];
301  l_weights[k - l_p - 1] = l_weights[k - l_p];
302 
303  for (unsigned int l = 1; l <= l_p; l++) {
304  unsigned int ind = k - l_p + l;
305  double alpha = l_knots[k + l] - l_x[j];
306  // if (vpMath::abs(alpha) == 0.0)
307  if (std::fabs(alpha) <= std::numeric_limits<double>::epsilon()) {
308  l_controlPoints[ind - 1] = l_controlPoints[ind];
309  l_weights[ind - 1] = l_weights[ind];
310  }
311  else {
312  alpha = alpha / (l_knots[k + l] - l_knots_tmp[i - l_p + l]);
313  l_controlPoints[ind - 1].set_i(alpha * l_controlPoints[ind - 1].get_i() +
314  (1.0 - alpha) * l_controlPoints[ind].get_i());
315  l_controlPoints[ind - 1].set_j(alpha * l_controlPoints[ind - 1].get_j() +
316  (1.0 - alpha) * l_controlPoints[ind].get_j());
317  l_weights[ind - 1] = alpha * l_weights[ind - 1] + (1.0 - alpha) * l_weights[ind];
318  }
319  }
320  l_knots[k] = l_x[j];
321  k--;
322  } while (j != 0);
323  }
324 
325  for (unsigned int j = 0; j < n; j++) {
326  l_controlPoints[j].set_ij(l_controlPoints[j].get_i() / l_weights[j], l_controlPoints[j].get_j() / l_weights[j]);
327  }
328 }
329 
330 void vpNurbs::refineKnotVectCurve(double *x, unsigned int r)
331 {
332  refineKnotVectCurve(x, r, p, knots, controlPoints, weights);
333 }
334 
335 unsigned int vpNurbs::removeCurveKnot(double l_u, unsigned int l_r, unsigned int l_num, double l_TOL, unsigned int l_s,
336  unsigned int l_p, std::vector<double> &l_knots,
337  std::vector<vpImagePoint> &l_controlPoints, std::vector<double> &l_weights)
338 {
339  unsigned int n = (unsigned int)l_controlPoints.size();
340  unsigned int m = n + l_p + 1;
341 
342  for (unsigned int j = 0; j < n; j++) {
343  l_controlPoints[j].set_ij(l_controlPoints[j].get_i() * l_weights[j], l_controlPoints[j].get_j() * l_weights[j]);
344  }
345 
346  unsigned int ord = l_p + 1;
347  double fout = (2 * l_r - l_s - l_p) / 2.;
348  unsigned int last = l_r - l_s;
349  unsigned int first = l_r - l_p;
350  unsigned int tblSize = 2 * l_p + 1;
351  vpImagePoint *tempP = new vpImagePoint[tblSize];
352  double *tempW = new double[tblSize];
353  vpImagePoint pt;
354  unsigned int t = 0;
355  double alfi = 0;
356  double alfj = 0;
357  unsigned int i, j;
358 
359  for (t = 0; t < l_num; t++) {
360  unsigned int off = first - 1;
361  tempP[0] = l_controlPoints[off];
362  tempW[0] = l_weights[off];
363  tempP[last + 1 - off] = l_controlPoints[last + 1];
364  tempW[last + 1 - off] = l_weights[last + 1];
365  i = first;
366  j = last;
367  unsigned int ii = 1;
368  unsigned int jj = last - off;
369  int remflag = 0;
370  while (j - i > t) {
371  alfi = (l_u - l_knots[i]) / (l_knots[i + ord + t] - l_knots[i]);
372  alfj = (l_u - l_knots[j - t]) / (l_knots[j + ord] - l_knots[j - t]);
373  pt.set_i((l_controlPoints[i].get_i() - (1.0 - alfi) * tempP[ii - 1].get_i()) / alfi);
374  tempP[ii].set_i((l_controlPoints[i].get_i() - (1.0 - alfi) * tempP[ii - 1].get_i()) / alfi);
375  tempP[ii].set_j((l_controlPoints[i].get_j() - (1.0 - alfi) * tempP[ii - 1].get_j()) / alfi);
376  tempW[ii] = ((l_weights[i] - (1.0 - alfi) * tempW[ii - 1]) / alfi);
377  tempP[jj].set_i((l_controlPoints[j].get_i() - alfj * tempP[jj + 1].get_i()) / (1.0 - alfj));
378  tempP[jj].set_j((l_controlPoints[j].get_j() - alfj * tempP[jj + 1].get_j()) / (1.0 - alfj));
379  tempW[jj] = ((l_weights[j] - alfj * tempW[jj + 1]) / (1.0 - alfj));
380  i++;
381  j--;
382  ii++;
383  jj--;
384  }
385 
386  if (j - i < t) {
387  double distancei = tempP[ii - 1].get_i() - tempP[jj + 1].get_i();
388  double distancej = tempP[ii - 1].get_j() - tempP[jj + 1].get_j();
389  double distancew = tempW[ii - 1] - tempW[jj + 1];
390  double distance = sqrt(vpMath::sqr(distancei) + vpMath::sqr(distancej) + vpMath::sqr(distancew));
391  if (distance <= l_TOL)
392  remflag = 1;
393  }
394  else {
395  alfi = (l_u - l_knots[i]) / (l_knots[i + ord + t] - l_knots[i]);
396  double distancei =
397  l_controlPoints[i].get_i() - (alfi * tempP[ii + t + 1].get_i() + (1.0 - alfi) * tempP[ii - 1].get_i());
398  double distancej =
399  l_controlPoints[i].get_j() - (alfi * tempP[ii + t + 1].get_j() + (1.0 - alfi) * tempP[ii - 1].get_j());
400  double distancew = l_weights[i] - (alfi * tempW[ii + t + 1] + (1.0 - alfi) * tempW[ii - 1]);
401  double distance = sqrt(vpMath::sqr(distancei) + vpMath::sqr(distancej) + vpMath::sqr(distancew));
402  if (distance <= l_TOL)
403  remflag = 1;
404  }
405  if (remflag == 0)
406  break;
407  else {
408  i = first;
409  j = last;
410  while (j - i > t) {
411  l_controlPoints[i].set_i(tempP[i - off].get_i());
412  l_controlPoints[i].set_j(tempP[i - off].get_j());
413  l_weights[i] = tempW[i - off];
414  l_controlPoints[j].set_i(tempP[j - off].get_i());
415  l_controlPoints[j].set_j(tempP[j - off].get_j());
416  l_weights[j] = tempW[j - off];
417  i++;
418  j--;
419  }
420  }
421  first--;
422  last++;
423  }
424  if (t == 0) {
425  delete[] tempP;
426  delete[] tempW;
427  return t;
428  }
429  for (unsigned int k = l_r + 1; k <= m; k++)
430  l_knots[k - t] = l_knots[k];
431  j = (unsigned int)fout;
432  i = j;
433  for (unsigned int k = 1; k < t; k++) {
434  if (k % 2)
435  i++;
436  else
437  j--;
438  }
439  for (unsigned int k = i + 1; k <= n; k++) {
440  l_controlPoints[j].set_i(l_controlPoints[k].get_i());
441  l_controlPoints[j].set_j(l_controlPoints[k].get_j());
442  l_weights[j] = l_weights[k];
443  j++;
444  }
445  for (unsigned int k = 0; k < t; k++) {
446  l_knots.erase(l_knots.end() - 1);
447  l_controlPoints.erase(l_controlPoints.end() - 1);
448  }
449 
450  for (unsigned int k = 0; k < l_controlPoints.size(); k++)
451  l_controlPoints[k].set_ij(l_controlPoints[k].get_i() / l_weights[k], l_controlPoints[k].get_j() / l_weights[k]);
452 
453  delete[] tempP;
454  delete[] tempW;
455  return t;
456 }
457 
458 unsigned int vpNurbs::removeCurveKnot(double l_u, unsigned int l_r, unsigned int l_num, double l_TOL)
459 {
460  return removeCurveKnot(l_u, l_r, l_num, l_TOL, 0, p, knots, controlPoints, weights);
461 }
462 
463 void vpNurbs::globalCurveInterp(std::vector<vpImagePoint> &l_crossingPoints, unsigned int l_p,
464  std::vector<double> &l_knots, std::vector<vpImagePoint> &l_controlPoints,
465  std::vector<double> &l_weights)
466 {
467  if (l_p == 0) {
468  // vpERROR_TRACE("Bad degree of the NURBS basis functions");
469  throw(vpException(vpException::badValue, "Bad degree of the NURBS basis functions"));
470  }
471 
472  l_knots.clear();
473  l_controlPoints.clear();
474  l_weights.clear();
475  unsigned int n = (unsigned int)l_crossingPoints.size() - 1;
476  unsigned int m = n + l_p + 1;
477 
478  double d = 0;
479  for (unsigned int k = 1; k <= n; k++)
480  d = d + distance(l_crossingPoints[k], 1, l_crossingPoints[k - 1], 1);
481 
482  // Compute ubar
483  std::vector<double> ubar;
484  ubar.push_back(0.0);
485  for (unsigned int k = 1; k < n; k++) {
486  ubar.push_back(ubar[k - 1] + distance(l_crossingPoints[k], 1, l_crossingPoints[k - 1], 1) / d);
487  }
488  ubar.push_back(1.0);
489 
490  // Compute the knot vector
491  for (unsigned int k = 0; k <= l_p; k++)
492  l_knots.push_back(0.0);
493 
494  double sum = 0;
495  for (unsigned int k = 1; k <= l_p; k++)
496  sum = sum + ubar[k];
497 
498  // Centripetal method
499  for (unsigned int k = 1; k <= n - l_p; k++) {
500  l_knots.push_back(sum / l_p);
501  sum = sum - ubar[k - 1] + ubar[l_p + k - 1];
502  }
503 
504  for (unsigned int k = m - l_p; k <= m; k++)
505  l_knots.push_back(1.0);
506 
507  vpMatrix A(n + 1, n + 1);
508  vpBasisFunction *N;
509 
510  for (unsigned int i = 0; i <= n; i++) {
511  unsigned int span = findSpan(ubar[i], l_p, l_knots);
512  N = computeBasisFuns(ubar[i], span, l_p, l_knots);
513  for (unsigned int k = 0; k <= l_p; k++)
514  A[i][span - l_p + k] = N[k].value;
515  delete[] N;
516  }
517  // vpMatrix Ainv = A.inverseByLU();
518  vpMatrix Ainv;
519  A.pseudoInverse(Ainv);
520  vpColVector Qi(n + 1);
521  vpColVector Qj(n + 1);
522  vpColVector Qw(n + 1);
523  for (unsigned int k = 0; k <= n; k++) {
524  Qi[k] = l_crossingPoints[k].get_i();
525  Qj[k] = l_crossingPoints[k].get_j();
526  }
527  Qw = 1;
528  vpColVector Pi = Ainv * Qi;
529  vpColVector Pj = Ainv * Qj;
530  vpColVector Pw = Ainv * Qw;
531 
532  vpImagePoint pt;
533  for (unsigned int k = 0; k <= n; k++) {
534  pt.set_ij(Pi[k], Pj[k]);
535  l_controlPoints.push_back(pt);
536  l_weights.push_back(Pw[k]);
537  }
538 }
539 
541 {
542  std::vector<vpImagePoint> v_crossingPoints;
543  l_crossingPoints.front();
544  vpMeSite s = l_crossingPoints.value();
545  vpImagePoint pt(s.get_ifloat(), s.get_jfloat());
546  vpImagePoint pt_1 = pt;
547  v_crossingPoints.push_back(pt);
548  l_crossingPoints.next();
549  while (!l_crossingPoints.outside()) {
550  s = l_crossingPoints.value();
551  pt.set_ij(s.get_ifloat(), s.get_jfloat());
552  if (vpImagePoint::distance(pt_1, pt) >= 10) {
553  v_crossingPoints.push_back(pt);
554  pt_1 = pt;
555  }
556  l_crossingPoints.next();
557  }
558  globalCurveInterp(v_crossingPoints, p, knots, controlPoints, weights);
559 }
560 
561 void vpNurbs::globalCurveInterp(const std::list<vpImagePoint> &l_crossingPoints)
562 {
563  std::vector<vpImagePoint> v_crossingPoints;
564  for (std::list<vpImagePoint>::const_iterator it = l_crossingPoints.begin(); it != l_crossingPoints.end(); ++it) {
565  v_crossingPoints.push_back(*it);
566  }
567  globalCurveInterp(v_crossingPoints, p, knots, controlPoints, weights);
568 }
569 
570 
571 void vpNurbs::globalCurveInterp(const std::list<vpMeSite> &l_crossingPoints)
572 {
573  std::vector<vpImagePoint> v_crossingPoints;
574  vpMeSite s = l_crossingPoints.front();
575  vpImagePoint pt(s.get_ifloat(), s.get_jfloat());
576  vpImagePoint pt_1 = pt;
577  v_crossingPoints.push_back(pt);
578  std::list<vpMeSite>::const_iterator it = l_crossingPoints.begin();
579  ++it;
580  for (; it != l_crossingPoints.end(); ++it) {
581  vpImagePoint pt_tmp(it->get_ifloat(), it->get_jfloat());
582  if (vpImagePoint::distance(pt_1, pt_tmp) >= 10) {
583  v_crossingPoints.push_back(pt_tmp);
584  pt_1 = pt_tmp;
585  }
586  }
587  globalCurveInterp(v_crossingPoints, p, knots, controlPoints, weights);
588 }
589 
590 void vpNurbs::globalCurveInterp() { globalCurveInterp(crossingPoints, p, knots, controlPoints, weights); }
591 
592 void vpNurbs::globalCurveApprox(std::vector<vpImagePoint> &l_crossingPoints, unsigned int l_p, unsigned int l_n,
593  std::vector<double> &l_knots, std::vector<vpImagePoint> &l_controlPoints,
594  std::vector<double> &l_weights)
595 {
596  l_knots.clear();
597  l_controlPoints.clear();
598  l_weights.clear();
599  unsigned int m = (unsigned int)l_crossingPoints.size() - 1;
600 
601  double d = 0;
602  for (unsigned int k = 1; k <= m; k++)
603  d = d + distance(l_crossingPoints[k], 1, l_crossingPoints[k - 1], 1);
604 
605  // Compute ubar
606  std::vector<double> ubar;
607  ubar.push_back(0.0);
608  for (unsigned int k = 1; k < m; k++)
609  ubar.push_back(ubar[k - 1] + distance(l_crossingPoints[k], 1, l_crossingPoints[k - 1], 1) / d);
610  ubar.push_back(1.0);
611 
612  // Compute the knot vector
613  for (unsigned int k = 0; k <= l_p; k++)
614  l_knots.push_back(0.0);
615 
616  d = (double)(m + 1) / (double)(l_n - l_p + 1);
617 
618  for (unsigned int j = 1; j <= l_n - l_p; j++) {
619  double i = floor(j * d);
620  double alpha = j * d - i;
621  l_knots.push_back((1.0 - alpha) * ubar[(unsigned int)i - 1] + alpha * ubar[(unsigned int)i]);
622  }
623 
624  for (unsigned int k = 0; k <= l_p; k++)
625  l_knots.push_back(1.0);
626 
627  // Compute Rk
628  std::vector<vpImagePoint> Rk;
629  vpBasisFunction *N;
630  for (unsigned int k = 1; k <= m - 1; k++) {
631  unsigned int span = findSpan(ubar[k], l_p, l_knots);
632  if (span == l_p && span == l_n) {
633  N = computeBasisFuns(ubar[k], span, l_p, l_knots);
634  vpImagePoint pt(l_crossingPoints[k].get_i() - N[0].value * l_crossingPoints[0].get_i() -
635  N[l_p].value * l_crossingPoints[m].get_i(),
636  l_crossingPoints[k].get_j() - N[0].value * l_crossingPoints[0].get_j() -
637  N[l_p].value * l_crossingPoints[m].get_j());
638  Rk.push_back(pt);
639  delete[] N;
640  }
641  else if (span == l_p) {
642  N = computeBasisFuns(ubar[k], span, l_p, l_knots);
643  vpImagePoint pt(l_crossingPoints[k].get_i() - N[0].value * l_crossingPoints[0].get_i(),
644  l_crossingPoints[k].get_j() - N[0].value * l_crossingPoints[0].get_j());
645  Rk.push_back(pt);
646  delete[] N;
647  }
648  else if (span == l_n) {
649  N = computeBasisFuns(ubar[k], span, l_p, l_knots);
650  vpImagePoint pt(l_crossingPoints[k].get_i() - N[l_p].value * l_crossingPoints[m].get_i(),
651  l_crossingPoints[k].get_j() - N[l_p].value * l_crossingPoints[m].get_j());
652  Rk.push_back(pt);
653  delete[] N;
654  }
655  else {
656  Rk.push_back(l_crossingPoints[k]);
657  }
658  }
659 
660  vpMatrix A(m - 1, l_n - 1);
661  // Compute A
662  for (unsigned int i = 1; i <= m - 1; i++) {
663  unsigned int span = findSpan(ubar[i], l_p, l_knots);
664  N = computeBasisFuns(ubar[i], span, l_p, l_knots);
665  for (unsigned int k = 0; k <= l_p; k++) {
666  if (N[k].i > 0 && N[k].i < l_n)
667  A[i - 1][N[k].i - 1] = N[k].value;
668  }
669  delete[] N;
670  }
671 
672  vpColVector Ri(l_n - 1);
673  vpColVector Rj(l_n - 1);
674  vpColVector Rw(l_n - 1);
675  for (unsigned int i = 0; i < l_n - 1; i++) {
676  double sum = 0;
677  for (unsigned int k = 0; k < m - 1; k++)
678  sum = sum + A[k][i] * Rk[k].get_i();
679  Ri[i] = sum;
680  sum = 0;
681  for (unsigned int k = 0; k < m - 1; k++)
682  sum = sum + A[k][i] * Rk[k].get_j();
683  Rj[i] = sum;
684  sum = 0;
685  for (unsigned int k = 0; k < m - 1; k++)
686  sum = sum + A[k][i]; // The crossing points weigths are equal to 1.
687  Rw[i] = sum;
688  }
689 
690  vpMatrix AtA = A.AtA();
691  vpMatrix AtAinv;
692  AtA.pseudoInverse(AtAinv);
693 
694  vpColVector Pi = AtAinv * Ri;
695  vpColVector Pj = AtAinv * Rj;
696  vpColVector Pw = AtAinv * Rw;
697 
698  vpImagePoint pt;
699  l_controlPoints.push_back(l_crossingPoints[0]);
700  l_weights.push_back(1.0);
701  for (unsigned int k = 0; k < l_n - 1; k++) {
702  pt.set_ij(Pi[k], Pj[k]);
703  l_controlPoints.push_back(pt);
704  l_weights.push_back(Pw[k]);
705  }
706  l_controlPoints.push_back(l_crossingPoints[m]);
707  l_weights.push_back(1.0);
708 }
709 
710 
711 void vpNurbs::globalCurveApprox(vpList<vpMeSite> &l_crossingPoints, unsigned int n)
712 {
713  std::vector<vpImagePoint> v_crossingPoints;
714  l_crossingPoints.front();
715  while (!l_crossingPoints.outside()) {
716  vpMeSite s = l_crossingPoints.value();
717  vpImagePoint pt(s.get_ifloat(), s.get_jfloat());
718  v_crossingPoints.push_back(pt);
719  l_crossingPoints.next();
720  }
721  globalCurveApprox(v_crossingPoints, p, n, knots, controlPoints, weights);
722 }
723 
724 
725 void vpNurbs::globalCurveApprox(const std::list<vpImagePoint> &l_crossingPoints, unsigned int n)
726 {
727  std::vector<vpImagePoint> v_crossingPoints;
728  for (std::list<vpImagePoint>::const_iterator it = l_crossingPoints.begin(); it != l_crossingPoints.end(); ++it) {
729  v_crossingPoints.push_back(*it);
730  }
731  globalCurveApprox(v_crossingPoints, p, n, knots, controlPoints, weights);
732 }
733 
734 void vpNurbs::globalCurveApprox(const std::list<vpMeSite> &l_crossingPoints, unsigned int n)
735 {
736  std::vector<vpImagePoint> v_crossingPoints;
737  for (std::list<vpMeSite>::const_iterator it = l_crossingPoints.begin(); it != l_crossingPoints.end(); ++it) {
738  vpImagePoint pt(it->get_ifloat(), it->get_jfloat());
739  v_crossingPoints.push_back(pt);
740  }
741  globalCurveApprox(v_crossingPoints, p, n, knots, controlPoints, weights);
742 }
743 
744 void vpNurbs::globalCurveApprox(unsigned int n)
745 {
746  globalCurveApprox(crossingPoints, p, n, knots, controlPoints, weights);
747 }
Class that provides tools to compute and manipulate a B-Spline curve.
Definition: vpBSpline.h:107
static unsigned int findSpan(double l_u, unsigned int l_p, const std::vector< double > &l_knots)
Definition: vpBSpline.cpp:79
static vpBasisFunction ** computeDersBasisFuns(double l_u, unsigned int l_i, unsigned int l_p, unsigned int l_der, const std::vector< double > &l_knots)
Definition: vpBSpline.cpp:228
static vpBasisFunction * computeBasisFuns(double l_u, unsigned int l_i, unsigned int l_p, const std::vector< double > &l_knots)
Definition: vpBSpline.cpp:142
Implementation of column vector and the associated operations.
Definition: vpColVector.h:163
error that can be emitted by ViSP classes.
Definition: vpException.h:59
@ badValue
Used to indicate that a value is not in the allowed range.
Definition: vpException.h:85
Class that defines a 2D point in an image. This class is useful for image processing and stores only ...
Definition: vpImagePoint.h:82
void set_j(double jj)
Definition: vpImagePoint.h:305
double get_j() const
Definition: vpImagePoint.h:125
static double distance(const vpImagePoint &iP1, const vpImagePoint &iP2)
void set_ij(double ii, double jj)
Definition: vpImagePoint.h:316
void set_i(double ii)
Definition: vpImagePoint.h:294
double get_i() const
Definition: vpImagePoint.h:114
Provide simple list management.
Definition: vpList.h:109
void next(void)
position the current element on the next one
Definition: vpList.h:244
void front(void)
Position the current element on the first element of the list.
Definition: vpList.h:317
bool outside(void) const
Test if the current element is outside the list (on the virtual element)
Definition: vpList.h:350
type & value(void)
return the value of the current element
Definition: vpList.h:263
static double sqr(double x)
Definition: vpMath.h:201
static long double comb(unsigned int n, unsigned int p)
Definition: vpMath.h:390
Implementation of a matrix and operations on matrices.
Definition: vpMatrix.h:146
vpMatrix AtA() const
Definition: vpMatrix.cpp:645
vpMatrix pseudoInverse(double svThreshold=1e-6) const
Definition: vpMatrix.cpp:2343
void insert(const vpMatrix &A, unsigned int r, unsigned int c)
Definition: vpMatrix.cpp:5781
Performs search in a given direction(normal) for a given distance(pixels) for a given 'site'....
Definition: vpMeSite.h:65
double get_ifloat() const
Definition: vpMeSite.h:192
double get_jfloat() const
Definition: vpMeSite.h:198
Class that provides tools to compute and manipulate a Non Uniform Rational B-Spline curve.
Definition: vpNurbs.h:92
static void refineKnotVectCurve(double *l_x, unsigned int l_r, unsigned int l_p, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:250
std::vector< double > weights
Vector which contains the weights associated to each control Points.
Definition: vpNurbs.h:95
static vpImagePoint computeCurvePoint(double l_u, unsigned int l_i, unsigned int l_p, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:53
static void curveKnotIns(double l_u, unsigned int l_k, unsigned int l_s, unsigned int l_r, unsigned int l_p, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:191
vpNurbs()
Definition: vpNurbs.cpp:49
static vpMatrix computeCurveDers(double l_u, unsigned int l_i, unsigned int l_p, unsigned int l_der, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:102
static vpImagePoint * computeCurveDersPoint(double l_u, unsigned int l_i, unsigned int l_p, unsigned int l_der, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:153
static void globalCurveApprox(std::vector< vpImagePoint > &l_crossingPoints, unsigned int l_p, unsigned int l_n, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:592
static unsigned int removeCurveKnot(double l_u, unsigned int l_r, unsigned int l_num, double l_TOL, unsigned int l_s, unsigned int l_p, std::vector< double > &l_knots, std::vector< vpImagePoint > &l_controlPoints, std::vector< double > &l_weights)
Definition: vpNurbs.cpp:335
void globalCurveInterp()
Definition: vpNurbs.cpp:590