29 gsVector3d<T> resultNormal = initFrom3DPlaneFit(points3D, margin);
30 if(outNormal != NULL) *outNormal = resultNormal;
34 gsCurveLoop<T>::gsCurveLoop(
const std::vector<T> angles3D,
const std::vector<T> lengths3D,
const std::vector<bool> isConvex, T margin,
const bool unitSquare)
36 this->initFrom3DByAngles(angles3D, lengths3D, isConvex, margin, unitSquare);
49 slv.allRoots(*c, roots, 0, u(0,0) ) ;
55 m_curves[i]->eval_into( xx, e );
57 for(
index_t r=0; r!=e.cols(); r++)
70 gsWarn <<
"isOnCurve only works for B-splines.\n";
76 bool gsCurveLoop<T>::isOn(gsMatrix<T>
const &u, T & paramResult, T tol )
78 for(
int i=0; i<this->numCurves(); i++)
79 if( m_curves[i]->isOn(u, tol) )
81 this->parameterOf(u,i,paramResult,tol);
88 bool gsCurveLoop<T>::is_ccw()
92 for (
typename std::vector< gsCurve<T> *>::const_iterator it=
93 m_curves.begin(); it!= m_curves.end() ; it++ )
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);
100 return ( result > 0 );
104 void gsCurveLoop<T>::reverse()
106 for (
typename std::vector< gsCurve<T> *>::const_iterator it=
107 m_curves.begin(); it!= m_curves.end() ; it++ )
110 size_t size=m_curves.size();
111 for (
size_t i=0;i<size/2;i++)
114 m_curves[i]=m_curves[size-i-1];
115 m_curves[size-i-1]=swap;
120 void gsCurveLoop<T>::translate(gsVector<T>
const & v)
122 for ( iterator it = m_curves.begin(); it != m_curves.end(); ++it)
129 typename std::vector< gsCurve<T> *>::const_iterator it;
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");
136 for ( it= m_curves.begin()+1; it!= m_curves.end() ; it++ )
149 for (
int i=0; i<n; ++i )
164 int n = m_curves.size();
166 for(
int i = startIndex; i != endIndex; i = (i + 1) % n)
168 result->insertCurve(m_curves[i]);
170 if(startIndex < endIndex)
172 m_curves.erase(m_curves.begin() + startIndex, m_curves.begin() + endIndex);
175 m_curves.push_back(newCurveThisFace);
179 m_curves.insert(m_curves.begin() + startIndex, newCurveThisFace);
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);
201 std::vector<T> result;
202 slv.allRoots(*c, result, direction, abscissa) ;
206 gsWarn<<
"Could not get intersection for this type of curve!\n";
207 return std::vector<T>();
216 switch (numEndPoints)
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++ )
244 interval = (*it)->parameterRange();
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;
259 int numCurves = m_curves.size();
261 assert(numCurves != 0);
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;
268 for(
int i = 1; i < numCurves; i++)
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 );
282 size_t n = signedAngles.size();
283 assert(lengths.size() == n);
285 std::vector<T> scaledAngles(signedAngles);
287 for(
size_t i = 0; i < n; i++)
289 totalAngle += scaledAngles[i];
293 gsWarn <<
"Treatment of non-positive total turning angle is not implemented.\n";
298 T angleScale = (T)(2.0 * EIGEN_PI) / totalAngle;
299 for(
size_t i = 0; i < n; i++)
301 scaledAngles[i] *= angleScale;
302 if(
math::abs(scaledAngles[i]) >= (T)(EIGEN_PI))
304 gsWarn <<
"Scaled turning angle exceeded pi, treatment of this has not been implemented.\n";
311 std::vector<T> cumulativeAngles(n);
312 cumulativeAngles[0] = 0;
313 for(
size_t i = 1; i < n; i++)
315 cumulativeAngles[i] = cumulativeAngles[i - 1] + scaledAngles[i];
323 for(
size_t j = 0; j < n; j++)
325 constraintCoefs(0, j) = math::cos(cumulativeAngles[j]);
326 constraintCoefs(1, j) = math::sin(cumulativeAngles[j]);
333 for(
size_t i = 0; i < n; i++)
335 invL(i, i) = (T)(1) / lengths[i];
340 constraintRhs.setZero();
344 for(
size_t j = 0; j < n; j++)
346 if(resultLengths(j, 0) <= 0)
348 gsWarn <<
"Could not compute satisfactory lengths for approximating polygon.\n";
359 for(
size_t i = 1; i < n; i++)
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]);
364 adjustPolygonToUnitSquare(result, margin);
371 std::vector<T> result;
372 size_t n = m_curves.size();
373 for(
size_t i = 0; i < n; i++)
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));
386 assert(corners.cols() == 2);
387 size_t n = corners.rows();
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();
395 scaleFactor += margin * (T)(2);
396 scaleFactor = (T)(1) / scaleFactor / std::max(maxx - minx, maxy - miny);
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++)
401 corners.row(i) *= scaleFactor;
402 corners.row(i) += shiftVector;
411 size_t n = points3D.size();
412 GISMO_ASSERT(n >= 3,
"Must have at least 3 points to construct a polygon");
420 for(
size_t i = 0; i < n; i++)
422 com += *(points3D[i]);
428 for(
size_t i = 0; i < n; i++)
430 shiftedPoints.col(i) = *(points3D[i]) - com;
434 gsEigen::JacobiSVD< gsEigen::Matrix<T, Dynamic, Dynamic> > svd(
435 shiftedPoints * shiftedPoints.transpose(), gsEigen::ComputeFullU);
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();
443 std::vector< gsVector<T> > projPts;
445 for(
size_t i = 0; i < n; i++)
448 projPts.push_back(x);
449 for(
size_t d = 0; d < 2; d++)
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);
457 T scale = (T(1) - (T)(2) * margin) / std::max((bbox(0, 1) - bbox(0, 0)), bbox(1, 1) - bbox(1, 0));
459 shift << margin - scale * bbox(0, 0), margin - scale * bbox(1, 0);
460 for(
size_t i = 0; i < n; i++)
462 projPts[i] = projPts[i] * scale + shift;
466 for(
size_t i = 0; i < n; i++)
468 size_t i1 = (i + 1) % n;
470 tcp << projPts[i](0), projPts[i](1), projPts[i1](0), projPts[i1](1);
478 return p1.cross(p2).normalized().transpose();
485 size_t n = isConvex.size();
486 assert(n == points3D.size());
488 std::vector<T> angles;
489 std::vector<T> lengths;
491 for(
size_t i = 0; i < n; i++)
499 T length0 = change0.norm();
500 T length1 = change1.norm();
501 T acosvalue= change0.dot(change1) / length0 / length1;
504 else if(acosvalue<-1)
506 angles.push_back( math::acos( acosvalue));
507 lengths.push_back(length1);
510 return initFrom3DByAngles(angles, lengths, isConvex, margin);
517 size_t np = isConvex.size();
519 for(
size_t i = 0; i < np; i++)
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);
527 for(
size_t i = 0; i < np; i++)
529 size_t iprev = (i + np - 1) % np;
533 cps.row(iprev * 4 + 3) = u;
537 cps.row(iprev * 4 + 2) = u + (prevPt - u) / 3;
538 cps.row(i * 4 + 1) = u + (nextPt - u) / 3;
542 cps.row(iprev * 4 + 2) = u - (nextPt - u) / 3;
543 cps.row(i * 4 + 1) = u - (prevPt - u) / 3;
547 adjustPolygonToUnitSquare(cps, margin);
550 for(
size_t i = 0; i < np; i++)
559 bool gsCurveLoop<T>::initFrom3DByAngles(
const std::vector<T>& angles3D,
const std::vector<T>& lengths3D,
const std::vector<bool>& isConvex, T margin,
bool unitSquare)
561 size_t n = isConvex.size();
562 assert(angles3D.size() == n);
563 assert(lengths3D.size() == n);
569 if (unitSquare==
true && n==4)
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)
575 corners4<<0,0,1,0,1,1,0,1;
576 for(
size_t i = 0; i < n; i++)
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) ));
587 std::vector<T> signedAngles(n);
588 for(
size_t i = 0; i < n; i++)
590 signedAngles[i] = (T)(isConvex[i]? 1: -1) * angles3D[i];
595 if(!success)
return false;
598 for(
size_t i = 0; i < n; i++)
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) ));
610 size_t nCur = m_curves.size();
611 T offset = minu + maxu;
612 for(
size_t i = 0; i < nCur; i++)
615 size_t ncp = coefs.rows();
616 for(
size_t j = 0; j < ncp; j++)
618 coefs(j, 0) = offset - coefs(j, 0);
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");
632 cp = curve(curveId).coefs();
634 int deg = curve(curveId).basis().degree(0);
637 gsWarn<<
"Pls extend to code to curves of degree more than 1, degree 1 is temporarily used.\n";
639 b = cp.row(cp.rows() - 1);
642 GISMO_ASSERT(supp.rows() == 1 && supp.cols() == 2,
"Unexpected support matrix");
644 u << (supp(0, 0) + supp(0, 1)) / (T)(2);
645 curve(curveId).eval_into (u, n);
646 n.transposeInPlace();
648 tcp0 << a(0,0),a(0,1),n(0,0),n(0,1);
650 tcp1 << n(0,0),n(0,1),b(0,0),b(0,1);
653 removeCurve(m_curves.begin()+curveId);
655 m_curves.insert(m_curves.begin()+curveId,tcurve0);
656 m_curves.insert(m_curves.begin()+curveId+1,tcurve1);
666 m_curves.erase(begin, end);
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