Skip to content

Commit

Permalink
FIX: CurveIntersection: Type fixes (#713)
Browse files Browse the repository at this point in the history 
* fix compiler warnings and code styling for curve intersection and B-spline
* bugs fix: uv gsMatrix->gsVector, return result in getPotentialIntersectionRanges

---------

Co-authored-by: jiyess <jiye@mail.dlut.edu.cn>
  • Loading branch information
@hverhelst @jiyess
hverhelst and jiyess committed Jun 7, 2024
1 parent 00ffb4f commit 67bdd38
Show file tree
Hide file tree
Showing 3 changed files with 376 additions and 296 deletions.
9 changes: 6 additions & 3 deletions examples/bSplineCurve_example.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,7 +89,8 @@ int main(int argc, char *argv[])
gsInfo << "Done. Re-run with --trim to learn basic trim/merge operations\n";

// Basic intersection operations between two BSpline curves - @Ye
if (intersect){
if (intersect)
{
gsMatrix<real_t> ctrPts1(4, 2);
ctrPts1 << 0,0, 1,1, 2,1, 3,1;
gsBSpline<real_t> bsp1(0, 1, 0, 3, ctrPts1);
Expand All @@ -101,10 +102,12 @@ int main(int argc, char *argv[])

gsInfo << intersectPts.size() << " intersections are found!" << "\n";
gsMatrix<> iPts(bsp1.geoDim(), intersectPts.size());
for (size_t j = 0; j < intersectPts.size(); ++j) {
for (size_t j = 0; j < intersectPts.size(); ++j)
{
iPts.col(j) = intersectPts[j].getPoint();
}
if (!intersectPts.empty()) {
if (!intersectPts.empty())
{
gsWriteParaviewPoints(iPts, "intersect");
}

Expand Down
111 changes: 63 additions & 48 deletions src/gsNurbs/gsBSpline.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,9 @@ void gsBSpline<T>::merge( gsGeometry<T> * otherG )
"Cannot merge a closed curve with anything." );

// check degree
const int mDeg = this ->basis().degree();
const int oDeg = other->basis().degree();
const int deg = math::max(mDeg,oDeg);
const short_t mDeg = this ->basis().degree();
const short_t oDeg = other->basis().degree();
const short_t deg = math::max(mDeg,oDeg);

other->gsBSpline::degreeElevate( deg - oDeg ); // degreeElevate(0) does nothing (and very quickly)
this ->gsBSpline::degreeElevate( deg - mDeg );
Expand All @@ -73,8 +73,8 @@ void gsBSpline<T>::merge( gsGeometry<T> * otherG )
mKnots.append( oKnots.begin()+deg+1, oKnots.end());

// merge coefficients
int n= this->coefsSize();
int skip = continuous ? 1 : 0;
unsigned n= this->coefsSize();
index_t skip = continuous ? 1 : 0;
this->m_coefs.conservativeResize( n + other->coefsSize() -skip, gsEigen::NoChange ) ;

this->m_coefs.block( n,0,other->coefsSize()-skip,other->geoDim() ) =
Expand All @@ -99,13 +99,15 @@ gsBSpline<T> gsBSpline<T>::segmentFromTo(T u0, T u1, T tolerance) const
KnotVectorType & knots = copy.basis().knots();

// some constants
const int p = basis().degree(); // degree
const short_t p = basis().degree(); // degree
const index_t multStart = p + 1 - knots.multiplicity(u0); // multiplicity
const index_t multEnd = p + 1 - knots.multiplicity(u1); // multiplicity
const index_t multEnd = p + 1 - knots.multiplicity(u1); // multiplicity

// insert the knot, such that its multiplicity is p+1
if (multStart>0) { copy.insertKnot(u0, 0, multStart); }
if (multEnd>0) { copy.insertKnot(u1, 0, multEnd ); }
if (multStart>0)
copy.insertKnot(u0, 0, multStart);
if (multEnd>0)
copy.insertKnot(u1, 0, multEnd );

gsMatrix<T>& coefs = copy.coefs();
const index_t tDim = coefs.cols();
Expand All @@ -117,7 +119,8 @@ gsBSpline<T> gsBSpline<T>::segmentFromTo(T u0, T u1, T tolerance) const
index_t nL2 = knots.uFind(u1).firstAppearance();
bool isEnd = math::abs(u1 - this->domainEnd()) < tolerance;
// if ( isEnd ) { nL2 += 1; } // Adjust for end parameter
if ( isEnd ) { nL2 = copy.numCoefs(); } // Adjust for end parameter
if ( isEnd )
nL2 = copy.numCoefs(); // Adjust for end parameter

// Prepare control points for new geometry
gsMatrix<T> coefRes = coefs.block(nL, 0, nL2-nL, tDim);
Expand All @@ -132,57 +135,69 @@ gsBSpline<T> gsBSpline<T>::segmentFromTo(T u0, T u1, T tolerance) const
}

template<class T>
gsMultiPatch<T> gsBSpline<T>::toBezier(T tolerance) const {
gsMultiPatch<T> bezierSegments;
gsMultiPatch<T> gsBSpline<T>::toBezier(T tolerance) const
{
gsMultiPatch<T> bezierSegments;

gsBSpline<T> currentSegment(*this);
gsBSpline<T> leftPart;
gsBSpline<T> currentSegment(*this);
gsBSpline<T> leftPart;

for (auto iter = this->knots().ubegin() + 1; iter != this->knots().uend() - 1; ++iter) {
currentSegment.splitAt(*iter, leftPart, currentSegment, tolerance);
bezierSegments.addPatch(leftPart);
}
for (auto iter = this->knots().ubegin() + 1; iter != this->knots().uend() - 1; ++iter)
{
currentSegment.splitAt(*iter, leftPart, currentSegment, tolerance);
bezierSegments.addPatch(leftPart);
}

bezierSegments.addPatch(currentSegment); // Add the last segment
return bezierSegments;
bezierSegments.addPatch(currentSegment); // Add the last segment
return bezierSegments;
}

template<class T>
std::vector<internal::gsCurveIntersectionResult<T>> gsBSpline<T>::intersect(const gsBSpline<T>& other,
T tolerance, T curvatureTolerance) const {
std::vector<internal::gsBoundingBoxPair<T>> hulls = internal::getPotentialIntersectionRanges<T>(*this, other, curvatureTolerance);

std::vector<internal::gsCurveIntersectionResult<T>> results;
for (const auto &hull : hulls) {
gsBSpline<T> crv1 = this->segmentFromTo(hull.b1.getRange().getMin(), hull.b1.getRange().getMax());
gsBSpline<T> crv2 = other.segmentFromTo(hull.b2.getRange().getMin(), hull.b2.getRange().getMax());

internal::gsCurveCurveDistanceSystem<T> obj(crv1, crv2);
gsMatrix<T, 2, 1> uv;
uv(0, 0) = 0.5 * (crv1.domainStart() + crv1.domainEnd());
uv(1, 0) = 0.5 * (crv2.domainStart() + crv2.domainEnd());
T distance = obj.compute(uv, tolerance);

if (distance < math::max((T)1e-10, tolerance)) {
internal::gsCurveIntersectionResult<T> result(uv(0), uv(1), 0.5 * (crv1.eval(uv.row(0)) + crv2.eval(uv.row(1))));
results.push_back(result);
T tolerance,
T curvatureTolerance) const
{
std::vector<internal::gsBoundingBoxPair<T>> hulls = internal::getPotentialIntersectionRanges<T>(*this, other, curvatureTolerance);

std::vector<internal::gsCurveIntersectionResult<T>> results;
for (typename std::vector<internal::gsBoundingBoxPair<T>>::const_iterator hull = hulls.begin(); hull!=hulls.end(); hull++)
{
gsBSpline<T> crv1 = this->segmentFromTo(hull->b1.getRange().getMin(), hull->b1.getRange().getMax());
gsBSpline<T> crv2 = other.segmentFromTo(hull->b2.getRange().getMin(), hull->b2.getRange().getMax());

internal::gsCurveCurveDistanceSystem<T> obj(crv1, crv2);
gsVector<T,2> uv;
uv(0) = 0.5 * (crv1.domainStart() + crv1.domainEnd());
uv(1) = 0.5 * (crv2.domainStart() + crv2.domainEnd());
T distance = obj.compute(uv, tolerance);

if (distance < math::max(1e-10, tolerance))
{
gsMatrix<T> uCrv1(1,1), uCrv2(1,1);
uCrv1(0,0) = uv(0);
uCrv2(0,0) = uv(1);
gsMatrix<T> pt = 0.5 * (crv1.eval(uCrv1) + crv2.eval(uCrv2));
internal::gsCurveIntersectionResult<T> result(uv(0), uv(1), pt);
results.push_back(result);
}
}
}

return results;
}

template<class T>
T gsBSpline<T>::pseudoCurvature() const {
int coefsSize = m_coefs.rows();

T len = (m_coefs.row(0)-m_coefs.row(coefsSize-1)).norm();
T total = 0.0;
for (int ipt = 0; ipt != coefsSize - 1; ++ipt) {
T dist = (m_coefs.row(ipt) - m_coefs.row(ipt + 1)).norm();
total += dist;
}
return total / len;
T gsBSpline<T>::pseudoCurvature() const
{
index_t coefsSize = m_coefs.rows();

T len = (m_coefs.row(0)-m_coefs.row(coefsSize-1)).norm();
T total = 0.0;
for (index_t ipt = 0; ipt != coefsSize - 1; ++ipt)
{
T dist = (m_coefs.row(ipt) - m_coefs.row(ipt + 1)).norm();
total += dist;
}
return total / len;
}

template<class T>
Expand Down
Loading

0 comments on commit 67bdd38

Please sign in to comment.