G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsCurveLoop.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 # include <gsModeling/gsCurveLoop.h>
17 # include <gsNurbs/gsKnotVector.h>
18 # include <gsNurbs/gsBSpline.h>
19 # include <gsNurbs/gsBSplineSolver.h>
21 
22 namespace gismo
23 {
24 
25 template<class T>
26 gsCurveLoop<T>::gsCurveLoop(const std::vector<gsVector3d<T> *> points3D, const std::vector<bool> , T margin, gsVector3d<T> *outNormal)
27 {
28  // attempt to choose a curve loop by using the plane of best fit
29  gsVector3d<T> resultNormal = initFrom3DPlaneFit(points3D, margin);
30  if(outNormal != NULL) *outNormal = resultNormal;
31 }
32 
33 template<class T>
34 gsCurveLoop<T>::gsCurveLoop(const std::vector<T> angles3D, const std::vector<T> lengths3D, const std::vector<bool> isConvex, T margin, const bool unitSquare)
35 {
36  this->initFrom3DByAngles(angles3D, lengths3D, isConvex, margin, unitSquare);
37 }
38 
39 template<class T>
40 bool gsCurveLoop<T>::parameterOf(gsMatrix<T> const &u, int i, T & result, T tol)
41 {
43  gsMatrix<T> e;
44  gsBSpline<T> * c = dynamic_cast<gsBSpline<T> *>(m_curves[i]);
45 
46  if ( c !=0 )
47  {
48  std::vector<T> roots;
49  slv.allRoots(*c, roots, 0, u(0,0) ) ;
50 
51  if(roots.size()==0)
52  return false;
53 
54  gsAsMatrix<T> xx(roots);
55  m_curves[i]->eval_into( xx, e );
56 
57  for(index_t r=0; r!=e.cols(); r++)
58  {
59  //gsDebug<<"matrix e (1,r): "<< e(1,r) <<"\n";
60  if( math::abs(e( 1, r )-u(1,0))< tol )
61  {
62  result = roots[r];
63  return true;
64  }
65  }
66  return false;
67  }
68  else
69  {
70  gsWarn <<"isOnCurve only works for B-splines.\n";
71  return false;
72  }
73 }
74 
75 template<class T>
76 bool gsCurveLoop<T>::isOn(gsMatrix<T> const &u, T & paramResult, T tol )
77 {
78  for(int i=0; i<this->numCurves(); i++)
79  if( m_curves[i]->isOn(u, tol) )
80  {
81  this->parameterOf(u,i,paramResult,tol);
82  return true;
83  }
84  return false;
85 }
86 
87 template<class T>
88 bool gsCurveLoop<T>::is_ccw()
89 {
90  T result(0);
91 
92  for ( typename std::vector< gsCurve<T> *>::const_iterator it=
93  m_curves.begin(); it!= m_curves.end() ; it++ )
94  {
95  index_t nrows = (*it)->coefs().rows();
96  for(index_t i=0; i!=nrows-1;++i)
97  result +=(*it)->coefs()(i,0)* (*it)->coefs()(i+1,1)-
98  (*it)->coefs()(i,1)* (*it)->coefs()(i+1,0);
99  }
100  return ( result > 0 );
101 }
102 
103 template<class T>
104 void gsCurveLoop<T>::reverse()
105 {
106  for ( typename std::vector< gsCurve<T> *>::const_iterator it=
107  m_curves.begin(); it!= m_curves.end() ; it++ )
108  (*it)->reverse();
109  gsCurve<T>* swap;
110  size_t size=m_curves.size();
111  for ( size_t i=0;i<size/2;i++)
112  {
113  swap = m_curves[i];
114  m_curves[i]=m_curves[size-i-1];
115  m_curves[size-i-1]=swap;
116  }
117 }
118 
119 template<class T>
120 void gsCurveLoop<T>::translate(gsVector<T> const & v)
121 {
122  for ( iterator it = m_curves.begin(); it != m_curves.end(); ++it)
123  (*it)->translate(v);
124 }
125 
126 template<class T>
128 {
129  typename std::vector< gsCurve<T> *>::const_iterator it;
130 
131  GISMO_ASSERT( m_curves.size(), "CurveLoop is empty.\n");
132  GISMO_ASSERT( m_curves.front(), "Something went wrong. Invalid pointer in gsCurveLoop member.\n");
133 
134  typename gsCurve<T>::uPtr loop = m_curves.front()->clone();
135 
136  for ( it= m_curves.begin()+1; it!= m_curves.end() ; it++ )
137  {
138  loop->merge( *it ) ;
139  }
140  return loop;
141 }
142 
143 template<class T>
144 gsMatrix<T> gsCurveLoop<T>::normal( int const & c, gsMatrix<T> const & u )
145 {
146  int n = u.cols();
147  gsMatrix<T> result(2,n);
148 
149  for ( int i=0; i<n; ++i )
150  {
151  gsMatrix<T> b = m_curves[c]->deriv( u.col(i) ) ;
152  b.normalize(); // unit tangent vector
153  result(0,i) = b(1);
154  result(1,i) = -b(0);
155  }
156  return result;
157 }
158 
159 template<class T>
160 typename gsCurveLoop<T>::uPtr gsCurveLoop<T>::split(int startIndex, int endIndex,
161  gsCurve<T> * newCurveThisFace,
162  gsCurve<T> * newCurveNewFace)
163 {
164  int n = m_curves.size();
165  uPtr result(new gsCurveLoop<T>(newCurveNewFace));
166  for(int i = startIndex; i != endIndex; i = (i + 1) % n)
167  {
168  result->insertCurve(m_curves[i]);
169  }
170  if(startIndex < endIndex)
171  {
172  m_curves.erase(m_curves.begin() + startIndex, m_curves.begin() + endIndex);
173  if(startIndex == 0)
174  {
175  m_curves.push_back(newCurveThisFace);
176  }
177  else
178  {
179  m_curves.insert(m_curves.begin() + startIndex, newCurveThisFace);
180  }
181  }
182  else
183  {
184  m_curves.erase(m_curves.begin() + startIndex, m_curves.end());
185  m_curves.erase(m_curves.begin(), m_curves.begin() + endIndex);
186  m_curves.push_back(newCurveThisFace);
187  }
188  return result;
189 }
190 
191 template <class T>
192 std::vector<T> gsCurveLoop<T>::lineIntersections(int const & direction , T const & abscissa)
193 {
194  // For now only Curve loops with ONE curve are supported
195  gsBSplineSolver<T> slv;
196  typename gsBSpline<T>::uPtr c = memory::convert_ptr<gsBSpline<T> >(singleCurve());
197  // = dynamic_cast<gsBSpline<T> *>(m_curves[0]))
198 
199  if ( c )
200  {
201  std::vector<T> result;
202  slv.allRoots(*c, result, direction, abscissa) ;
203  /*slv.allRoots(*c, result, direction, abscissa,1e-7,1000) ;*/
204  return result;
205  }
206  gsWarn<<"Could not get intersection for this type of curve!\n";
207  return std::vector<T>();
208 }
209 
210 template <class T>
211 gsMatrix<T> gsCurveLoop<T>::sample(int npoints, int numEndPoints) const
212 {
213  GISMO_ASSERT(npoints>=2, "");
214  GISMO_ASSERT(numEndPoints>=0 && numEndPoints<=2, "");
215  int np; // new number of points
216  switch (numEndPoints)
217  {
218  case (0):
219  np = npoints-2;
220  break;
221  case (1):
222  np = npoints-1;
223  break;
224  case (2):
225  np = npoints;
226  break;
227  default:
228  np = 0;
229  break;
230  }
231 
232  gsMatrix<T> u(2, m_curves.size() * np);
233  int i=0;
234  gsMatrix<T> interval;
235  gsMatrix<T> pts(1,np);
236  gsMatrix<T> uCols(2,np);
237  T a,b;
238  int firstInd=0;
239  int secondInd=np-1;
240  if (numEndPoints==0) {firstInd=1;secondInd=npoints-2;}
241  typename std::vector< gsCurve<T> *>::const_iterator it;
242  for ( it= m_curves.begin(); it!= m_curves.end() ; it++ )
243  {
244  interval = (*it)->parameterRange();
245  a = interval(0,0);
246  b = interval(0,1);
247  for (int ii=firstInd;ii<=secondInd;ii++) pts(0,ii-firstInd)= a + (T)(ii)*(b-a)/(T)(npoints-1);
248  (*it)->eval_into( pts, uCols );
249  u.middleCols( i * np,np ) = uCols;
250  i++ ;
251  };
252  return u;
253 }
254 
255 template <class T>
257 {
258  gsMatrix<T, 2, 2> result;
259  int numCurves = m_curves.size();
260  T offset = 0.0001;
261  assert(numCurves != 0); // bounding box does not exist if there are no curves
262  gsMatrix<T> *coefs = &(m_curves[0]->coefs());
263  result(0, 0) = coefs->col(0).minCoeff() - offset;
264  result(1, 0) = coefs->col(1).minCoeff() - offset;
265  result(0, 1) = coefs->col(0).maxCoeff() + offset;
266  result(1, 1) = coefs->col(1).maxCoeff() + offset;
267 
268  for(int i = 1; i < numCurves; i++)
269  {
270  coefs = &(m_curves[i]->coefs());
271  result(0, 0) = math::min(result(0, 0), coefs->col(0).minCoeff()-offset );
272  result(1, 0) = math::min(result(1, 0), coefs->col(1).minCoeff()-offset );
273  result(0, 1) = math::max(result(0, 1), coefs->col(0).maxCoeff()+offset );
274  result(1, 1) = math::max(result(1, 1), coefs->col(1).maxCoeff()+offset );
275  }
276  return result;
277 }
278 
279 template<class T>
280 bool gsCurveLoop<T>::approximatingPolygon(const std::vector<T> &signedAngles, const std::vector<T> &lengths, T margin, gsMatrix<T> &result)
281 {
282  size_t n = signedAngles.size();
283  assert(lengths.size() == n);
284 
285  std::vector<T> scaledAngles(signedAngles); // copy
286  T totalAngle(0);
287  for(size_t i = 0; i < n; i++)
288  {
289  totalAngle += scaledAngles[i];
290  }
291  if(totalAngle <= 0)
292  {
293  gsWarn << "Treatment of non-positive total turning angle is not implemented.\n";
294  return false;
295  }
296  // Scale the turning angles so they add to 2 * pi. These will be the
297  // angles we use at each vertex in the domain.
298  T angleScale = (T)(2.0 * EIGEN_PI) / totalAngle;
299  for(size_t i = 0; i < n; i++)
300  {
301  scaledAngles[i] *= angleScale;
302  if(math::abs(scaledAngles[i]) >= (T)(EIGEN_PI))
303  {
304  gsWarn << "Scaled turning angle exceeded pi, treatment of this has not been implemented.\n";
305  return false;
306  }
307  }
308 
309  // Compute the cumulative turning angle (i.e. the total angle turned
310  // by the time we reach a given vertex).
311  std::vector<T> cumulativeAngles(n);
312  cumulativeAngles[0] = 0;
313  for(size_t i = 1; i < n; i++)
314  {
315  cumulativeAngles[i] = cumulativeAngles[i - 1] + scaledAngles[i];
316  }
317 
318  // Build a matrix to provide the coefficients of the following
319  // constraints (RHS of equations is built later):
320  // total change of x coordinate (as we go round the polygon) == 0;
321  // total change of y coordinate (as we go round the polygon) == 0;
322  gsMatrix<T> constraintCoefs(2, n);
323  for(size_t j = 0; j < n; j++)
324  {
325  constraintCoefs(0, j) = math::cos(cumulativeAngles[j]);
326  constraintCoefs(1, j) = math::sin(cumulativeAngles[j]);
327  }
328 
329  // Find a critical point of the following cost function subject to the above constraints:
330  // (1 / 2) * ( sum of ( (length i) / (original length i) ) ^ 2 - 1 ).
331  gsMatrix<T> invL(n, n);
332  invL.setZero();
333  for(size_t i = 0; i < n; i++)
334  {
335  invL(i, i) = (T)(1) / lengths[i];
336  }
337  gsMatrix<T> ones(n, 1);
338  ones.setOnes();
339  gsMatrix<T> constraintRhs(2, 1);
340  constraintRhs.setZero();
341 
342  gsMatrix<T> resultLengths = optQuadratic(invL, ones, constraintCoefs, constraintRhs);
343 
344  for(size_t j = 0; j < n; j++)
345  {
346  if(resultLengths(j, 0) <= 0)
347  {
348  gsWarn << "Could not compute satisfactory lengths for approximating polygon.\n";
349  return false;
350  }
351  }
352 
353  // Work out the positions of the corners of this polygon.
354  // (Corner n == corner 0). Initially, we put the first vertex at (0,0). Then
355  // we apply a scale & shift to fit it all inside the box [0,1]x[0,1].
356  result.resize(n, 2);
357  result(0, 0) = 0; // set up the first corner
358  result(0, 1) = 0;
359  for(size_t i = 1; i < n; i++) // set up corners other than the first
360  {
361  result(i, 0) = result(i - 1, 0) + resultLengths(i - 1) * math::cos(cumulativeAngles[i - 1]);
362  result(i, 1) = result(i - 1, 1) + resultLengths(i - 1) * math::sin(cumulativeAngles[i - 1]);
363  }
364  adjustPolygonToUnitSquare(result, margin);
365  return true;
366 }
367 
368 template<class T>
369 std::vector<T> gsCurveLoop<T>::domainSizes() const
370 {
371  std::vector<T> result;
372  size_t n = m_curves.size();
373  for(size_t i = 0; i < n; i++)
374  {
375  gsMatrix<T> range = m_curves[i]->parameterRange();
376  GISMO_ASSERT(range.rows() == 1 && range.cols() == 2, "Expected 1x2 matrix for parameter range of curve");
377  result.push_back(range(0, 1) - range(0, 0));
378  }
379 
380  return result;
381 }
382 
383 template <class T>
385 {
386  assert(corners.cols() == 2);
387  size_t n = corners.rows();
388 
389  T minx = corners.col(0).minCoeff();
390  T maxx = corners.col(0).maxCoeff();
391  T miny = corners.col(1).minCoeff();
392  T maxy = corners.col(1).maxCoeff();
393 
394  T scaleFactor(1);
395  scaleFactor += margin * (T)(2);
396  scaleFactor = (T)(1) / scaleFactor / std::max(maxx - minx, maxy - miny);
397  gsMatrix<T> shiftVector(1, 2);
398  shiftVector << 0.5 - 0.5 * scaleFactor * (maxx + minx), 0.5 - 0.5 * scaleFactor * (maxy + miny);
399  for(size_t i = 0; i < n; i++)
400  {
401  corners.row(i) *= scaleFactor;
402  corners.row(i) += shiftVector;
403  }
404 }
405 
406 template<class T>
407 gsVector3d<T> gsCurveLoop<T>::initFrom3DPlaneFit(const std::vector<gsVector3d<T> *> points3D, T margin)
408 {
409  // attempt to construct a curve loop using the plane of best fit
410 
411  size_t n = points3D.size(); // number of corners
412  GISMO_ASSERT(n >= 3, "Must have at least 3 points to construct a polygon");
413 
414  freeAll(m_curves);
415  m_curves.clear();
416 
417  // find the center of mass of the points
418  gsVector3d<T> com;
419  com << 0, 0, 0;
420  for(size_t i = 0; i < n; i++)
421  {
422  com += *(points3D[i]);
423  }
424  com /= n;
425 
426  // construct a matrix with the shifted points
427  gsMatrix<T> shiftedPoints(3, n);
428  for(size_t i = 0; i < n; i++)
429  {
430  shiftedPoints.col(i) = *(points3D[i]) - com;
431  }
432 
433  // compute a singular value decomposition
434  gsEigen::JacobiSVD< gsEigen::Matrix<T, Dynamic, Dynamic> > svd(
435  shiftedPoints * shiftedPoints.transpose(), gsEigen::ComputeFullU);
436 
437  // extract plane projection matrix from the svd
438  gsMatrix<T> svd_u(svd.matrixU());
439  GISMO_ASSERT(svd_u.rows() == 3 && svd_u.cols() == 3, "Unexpected svd matrix result");
440  gsMatrix<T> projMat = svd_u.block(0, 0, 3, 2).transpose();
441 
442  // project all the points down onto the plane, keep track of min, max of coords
443  std::vector< gsVector<T> > projPts;
444  gsMatrix<T> bbox(2, 2); // first column is min of each dimension, second column is max
445  for(size_t i = 0; i < n; i++)
446  {
447  gsVector<T> x = projMat * shiftedPoints.col(i);
448  projPts.push_back(x);
449  for(size_t d = 0; d < 2; d++)
450  {
451  if(i == 0 || x(d) < bbox(d, 0)) bbox(d, 0) = x(d);
452  if(i == 0 || x(d) > bbox(d, 1)) bbox(d, 1) = x(d);
453  }
454  }
455 
456  // scale and shift so everything fits inside the bounding box minus margin.
457  T scale = (T(1) - (T)(2) * margin) / std::max((bbox(0, 1) - bbox(0, 0)), bbox(1, 1) - bbox(1, 0));
458  gsVector<T> shift(2);
459  shift << margin - scale * bbox(0, 0), margin - scale * bbox(1, 0);
460  for(size_t i = 0; i < n; i++)
461  {
462  projPts[i] = projPts[i] * scale + shift;
463  }
464 
465  // interpolate linearly
466  for(size_t i = 0; i < n; i++)
467  {
468  size_t i1 = (i + 1) % n;
469  gsMatrix<T> tcp(2, 2);
470  tcp << projPts[i](0), projPts[i](1), projPts[i1](0), projPts[i1](1);
471  this->insertCurve(new gsBSpline<T>( 0, 1, 0, 1, give(tcp) ));
472  }
473 
474  // return the normal to the plane of best fit. we could just grab the last
475  // column from svd_u but then we'd have to worry about the orientation.
476  gsVector3d<T> p1 = projMat.row(0);
477  gsVector3d<T> p2 = projMat.row(1);
478  return p1.cross(p2).normalized().transpose();
479 
480 }
481 
482 template<class T>
483 bool gsCurveLoop<T>::initFrom3DByAngles(const std::vector<gsVector3d<T> *> points3D, const std::vector<bool> isConvex, T margin)
484 {
485  size_t n = isConvex.size();
486  assert(n == points3D.size());
487 
488  std::vector<T> angles;
489  std::vector<T> lengths;
490 
491  for(size_t i = 0; i < n; i++)
492  {
493  gsVector3d<T> *prev = points3D[(i + n - 1) % n];
494  gsVector3d<T> *cur = points3D[i];
495  gsVector3d<T> *next = points3D[(i + 1) % n];
496 
497  gsVector3d<T> change0 = *cur - *prev;
498  gsVector3d<T> change1 = *next - *cur;
499  T length0 = change0.norm();
500  T length1 = change1.norm();
501  T acosvalue= change0.dot(change1) / length0 / length1;
502  if(acosvalue>1)
503  acosvalue=1;
504  else if(acosvalue<-1)
505  acosvalue=-1;
506  angles.push_back( math::acos( acosvalue));
507  lengths.push_back(length1);
508  }
509 
510  return initFrom3DByAngles(angles, lengths, isConvex, margin);
511 }
512 
513 template<class T>
514 void gsCurveLoop<T>::initFromIsConvex(const std::vector<bool> isConvex, T margin)
515 {
516  // find the corners of a regular polygon
517  size_t np = isConvex.size();
518  gsMatrix<T> corners(np, 2);
519  for(size_t i = 0; i < np; i++)
520  {
521  T angle = (T)i * (T)EIGEN_PI * (T)(2) / (T)(np);
522  corners(i, 0) = math::cos(angle);
523  corners(i, 1) = math::sin(angle);
524  }
525  // choose control points of cubic splines which ensure the correct angle signs
526  gsMatrix<T> cps(np * 4, 2);
527  for(size_t i = 0; i < np; i++)
528  {
529  size_t iprev = (i + np - 1) % np;
530  gsMatrix<T> prevPt = corners.row(iprev);
531  gsMatrix<T> u = corners.row(i);
532  gsMatrix<T> nextPt = corners.row((i + 1) % np);
533  cps.row(iprev * 4 + 3) = u;
534  cps.row(i * 4) = u;
535  if(isConvex[i])
536  {
537  cps.row(iprev * 4 + 2) = u + (prevPt - u) / 3;
538  cps.row(i * 4 + 1) = u + (nextPt - u) / 3;
539  }
540  else
541  {
542  cps.row(iprev * 4 + 2) = u - (nextPt - u) / 3;
543  cps.row(i * 4 + 1) = u - (prevPt - u) / 3;
544  }
545  }
546  // put all the control points inside a box
547  adjustPolygonToUnitSquare(cps, margin);
548  // construct the splines
549  curves().clear();
550  for(size_t i = 0; i < np; i++)
551  {
552  gsMatrix<T> tcp = cps.block(i * 4, 0, 4, 2);
553  this->insertCurve(new gsBSpline<T>( 0, 1, 0, 3, give(tcp) ));
554  }
555 }
556 
557 
558 template<class T>
559 bool gsCurveLoop<T>::initFrom3DByAngles(const std::vector<T>& angles3D, const std::vector<T>& lengths3D, const std::vector<bool>& isConvex, T margin, bool unitSquare)
560 {
561  size_t n = isConvex.size();
562  assert(angles3D.size() == n);
563  assert(lengths3D.size() == n);
564 
565  freeAll(m_curves);
566  m_curves.clear();
567 
568  gsMatrix<T> corners4(4, 2);
569  if (unitSquare==true && n==4)
570  {
571  bool isAllConvex=true;
572  for (size_t i=0;i<4;i++) {if (isConvex.at(i)==false) {isAllConvex=false;break;}}
573  if (isAllConvex==true)
574  {
575  corners4<<0,0,1,0,1,1,0,1;
576  for(size_t i = 0; i < n; i++)
577  {
578  gsMatrix<T> tcp(2, 2);
579  tcp << corners4(i, 0), corners4(i, 1), corners4((i + 1) % n, 0), corners4((i + 1) % n, 1);
580  this->insertCurve(new gsBSpline<T>( 0, 1, 0, 1, give(tcp) ));
581  }
582  return true;
583  }
584  }
585 
586  // compute the 2D angles corresponding to the 3D ones
587  std::vector<T> signedAngles(n);
588  for(size_t i = 0; i < n; i++)
589  {
590  signedAngles[i] = (T)(isConvex[i]? 1: -1) * angles3D[i];
591  }
592 
593  gsMatrix<T> corners;
594  bool success = gsCurveLoop<T>::approximatingPolygon(signedAngles, lengths3D, margin, corners);
595  if(!success) return false;
596 
597  // create a loop of B-splines that are all straight lines.
598  for(size_t i = 0; i < n; i++)
599  {
600  gsMatrix<T> tcp(2, 2);
601  tcp << corners(i, 0), corners(i, 1), corners((i + 1) % n, 0), corners((i + 1) % n, 1);
602  this->insertCurve(new gsBSpline<T>( 0, 1, 0, 1, give(tcp) ));
603  }
604  return true;
605 }
606 
607 template<class T>
608 void gsCurveLoop<T>::flip1(T minu, T maxu)
609 {
610  size_t nCur = m_curves.size();
611  T offset = minu + maxu;
612  for(size_t i = 0; i < nCur; i++)
613  {
614  gsMatrix<T> &coefs = m_curves[i]->coefs();
615  size_t ncp = coefs.rows();
616  for(size_t j = 0; j < ncp; j++)
617  {
618  coefs(j, 0) = offset - coefs(j, 0);
619  }
620  }
621 }
622 
623 template<class T>
624 gsMatrix<T> gsCurveLoop<T>::splitCurve(size_t curveId, T lengthRatio)
625 {
626  GISMO_UNUSED(lengthRatio);
627  GISMO_ASSERT(lengthRatio>0 && lengthRatio<1, "the second parameter *lengthRatio* must be between 0 and 1");
628  GISMO_ASSERT(curveId < m_curves.size(), "the first parameter *curveID* exceeds the number of curves of the loop");
629  // the two vertices a and b of the curve curveId, counting counterclockwise from z+
630  gsMatrix<T> a(2,1),b(2,1),n(2,1);
631  gsMatrix<T> cp,tcp0(2,2),tcp1(2,2);
632  cp = curve(curveId).coefs();
633  //int deg = curve(curveId).degree();
634  int deg = curve(curveId).basis().degree(0);
635  //GISMO_ASSERT(deg==1,"Pls extend to code to curves of degree more than 1");
636  if (deg!=1)
637  gsWarn<<"Pls extend to code to curves of degree more than 1, degree 1 is temporarily used.\n";
638  a = cp.row(0);
639  b = cp.row(cp.rows() - 1);
640  //n = (1-lengthRatio)*a + lengthRatio*b; // new vertex
641  gsMatrix<T> supp = curve(curveId).support();
642  GISMO_ASSERT(supp.rows() == 1 && supp.cols() == 2, "Unexpected support matrix");
643  gsMatrix<T> u(1, 1);
644  u << (supp(0, 0) + supp(0, 1)) / (T)(2);
645  curve(curveId).eval_into (u, n);
646  n.transposeInPlace();
647  // create the two new curves
648  tcp0 << a(0,0),a(0,1),n(0,0),n(0,1);
649  gsBSpline<T> * tcurve0 = new gsBSpline<T>( 0, 1, 0, 1, give(tcp0) );
650  tcp1 << n(0,0),n(0,1),b(0,0),b(0,1);
651  gsBSpline<T> * tcurve1 = new gsBSpline<T>( 0, 1, 0, 1, give(tcp1) );
652  // remove the original curve
653  removeCurve(m_curves.begin()+curveId);
654  // add the two new curves
655  m_curves.insert(m_curves.begin()+curveId,tcurve0);
656  m_curves.insert(m_curves.begin()+curveId+1,tcurve1);
657 
658  return n;
659 }
660 
661 template <class T>
662 void gsCurveLoop<T>::removeCurves(iterator begin, iterator end)
663 {
664  // free the curves before removing them
665  freeAll(begin, end);
666  m_curves.erase(begin, end);
667 }
668 
669 }
Knot vector for B-splines.
A fixed-size, statically allocated 3D vector.
Definition: gsVector.h:218
gsMatrix< T > normal(int const &c, gsMatrix< T > const &u)
Definition: gsCurveLoop.hpp:144
Abstract base class representing a curve.
Definition: gsCurve.h:30
Creates a mapped object or data pointer to a matrix without copying data.
Definition: gsLinearAlgebra.h:126
S give(S &x)
Definition: gsMemory.h:266
memory::unique_ptr< gsBSpline > uPtr
Unique pointer for gsBSpline.
Definition: gsBSpline.h:63
#define index_t
Definition: gsConfig.h:32
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
A B-spline function of one argument, with arbitrary target dimension.
Definition: gsBSpline.h:50
uPtr split(int startIndex, int endIndex, gsCurve< T > *newCurveThisFace, gsCurve< T > *newCurveNewFace)
Definition: gsCurveLoop.hpp:160
static bool approximatingPolygon(const std::vector< T > &signedAngles, const std::vector< T > &lengths, T margin, gsMatrix< T > &result)
Definition: gsCurveLoop.hpp:280
#define gsWarn
Definition: gsDebug.h:50
gsMatrix< T > sample(int npoints=50, int numEndPoints=2) const
Sample npoints uniformly distributed (in parameter domain) points on the loop.
Definition: gsCurveLoop.hpp:211
gsVector3d< T > initFrom3DPlaneFit(const std::vector< gsVector3d< T > * > points3D, T margin)
Initialize a curve loop from some 3d vertices by projecting them onto the plane of best fit and const...
Definition: gsCurveLoop.hpp:407
memory::unique_ptr< gsCurveLoop > uPtr
Unique pointer for gsCurveLoop.
Definition: gsCurveLoop.h:47
GISMO_DEPRECATED index_t direction(index_t s)
Returns the parametric direction that corresponds to side s.
Definition: gsBoundary.h:1048
void flip1(T minu=0, T maxu=1)
flip a curve loop in the u direction
Definition: gsCurveLoop.hpp:608
void freeAll(It begin, It end)
Frees all pointers in the range [begin end)
Definition: gsMemory.h:312
Represents a B-spline curve/function with one parameter.
static void adjustPolygonToUnitSquare(gsMatrix< T > &corners, T const margin)
Definition: gsCurveLoop.hpp:384
gsMatrix< T > splitCurve(size_t curveId, T lengthRatio=.5)
Definition: gsCurveLoop.hpp:624
bool initFrom3DByAngles(const std::vector< gsVector3d< T > * > points3D, const std::vector< bool > isConvex, T margin)
Definition: gsCurveLoop.hpp:483
Utility functions required by gsModeling classes.
gsCurveLoop()
Default empty constructor.
Definition: gsCurveLoop.h:51
#define GISMO_UNUSED(x)
Definition: gsDebug.h:112
memory::unique_ptr< gsCurve > uPtr
Unique pointer for gsCurve.
Definition: gsCurve.h:38
void initFromIsConvex(const std::vector< bool > isConvex, T margin)
Initialize a curve loop from a list of bools indicating whether the corners are convex or not...
Definition: gsCurveLoop.hpp:514
gsCurve< T >::uPtr singleCurve() const
Return the curve-loop as a single new B-Spline curve.
Definition: gsCurveLoop.hpp:127
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4488
Definition: gsBSplineSolver.h:113
std::vector< T > domainSizes() const
return a vector containing whose nth entry is the length of the domin of the nth curve.
Definition: gsCurveLoop.hpp:369
gsMatrix< T > optQuadratic(gsMatrix< T > const &A, gsMatrix< T > const &b, gsMatrix< T > const &C, gsMatrix< T > const &d)
Find X which solves: min (AX-b)^T (AX-b), s.t. CX=d.
Definition: gsModelingUtils.hpp:253
Provides classes and functions to solve equations involving B-splines.
A closed loop given by a collection of curves.
Definition: gsCurveLoop.h:36
Interface for gsCurveLoop class, representing a loop of curves, in anticlockwise order.
std::vector< T > lineIntersections(int const &direction, T const &abscissa)
Definition: gsCurveLoop.hpp:192