diff --git a/src/geometry/BSplineAlgorithms.cpp b/src/geometry/BSplineAlgorithms.cpp index 2be5afa..75ca93c 100644 --- a/src/geometry/BSplineAlgorithms.cpp +++ b/src/geometry/BSplineAlgorithms.cpp @@ -435,6 +435,17 @@ std::vector BSplineAlgorithms::knotsFromCurveParameters(std::vector BSplineAlgorithms::toBSplines(const std::vector& curves) +{ + std::vector result; + + std::transform(curves.begin(), curves.end(), std::back_inserter(result), [](const auto& curve) { + return GeomConvert::CurveToBSplineCurve(curve); + }); + + return result; +} + double BSplineAlgorithms::scale(const TColgp_Array2OfPnt& points) { double theScale = 0.; @@ -704,7 +715,7 @@ ApproxResult BSplineAlgorithms::reparametrizeBSplineContinuouslyApprox(const Han } #endif - ApproxResult result = approximationObj.FitCurveOptimal(parameters); + ApproxResult result = approximationObj.FitCurveOptimal(parameters, 1); assert(!result.curve.IsNull()); diff --git a/src/geometry/BSplineAlgorithms.h b/src/geometry/BSplineAlgorithms.h index fab26d9..7e41daf 100644 --- a/src/geometry/BSplineAlgorithms.h +++ b/src/geometry/BSplineAlgorithms.h @@ -279,6 +279,9 @@ class BSplineAlgorithms /// Trims a bspline curve GEOML_EXPORT static Handle(Geom_BSplineCurve) trimCurve(const Handle(Geom_BSplineCurve)& curve, double umin, double umax); + /// Converts a curve array into a B-spline array + GEOML_EXPORT static std::vector toBSplines(const std::vector& curves); + /// Concatenates a list of bspline curves GEOML_EXPORT static Handle(Geom_BSplineCurve) concatCurves(std::vector curves, bool parByLength=true, double tolerance = 1e-6); diff --git a/src/geometry/curve-networks/InterpolateCurveNetwork.cpp b/src/geometry/curve-networks/InterpolateCurveNetwork.cpp index bb9dc41..72d9011 100644 --- a/src/geometry/curve-networks/InterpolateCurveNetwork.cpp +++ b/src/geometry/curve-networks/InterpolateCurveNetwork.cpp @@ -29,6 +29,7 @@ #include #include +#include #include #include @@ -38,28 +39,56 @@ namespace geoml InterpolateCurveNetwork::InterpolateCurveNetwork(const std::vector &profiles, const std::vector &guides, double spatialTol) + : InterpolateCurveNetwork(BSplineAlgorithms::toBSplines(profiles), + BSplineAlgorithms::toBSplines(guides), + spatialTol) +{ +} + +InterpolateCurveNetwork::InterpolateCurveNetwork(const std::vector& profiles, + const std::vector& guides, + double spatialTol) : m_hasPerformed(false) , m_spatialTol(spatialTol) { - // check whether there are any u-directional and v-directional B-splines in the vectors if (profiles.size() < 2) { - throw Error("There must be at least two profiles for the curve network interpolation.", geoml::MATH_ERROR); + throw Error("There must be at least two profiles for the curve network interpolation.", MATH_ERROR); } - if (guides.size() < 2) { - throw Error("There must be at least two guides for the curve network interpolation.", geoml::MATH_ERROR); + if (guides.size() < 2) { + throw Error("There must be at least two guides for the curve network interpolation.", MATH_ERROR); + } + + std::vector uniqueProfiles; + for (const auto& profile : profiles) { + const bool isUnique = std::none_of(uniqueProfiles.begin(), uniqueProfiles.end(), [&](const Handle(Geom_BSplineCurve)& curve) { + return profile->IsEqual(curve, Precision::Confusion()); + }); + if (isUnique) { + uniqueProfiles.push_back(profile); + } + } + + std::vector uniqueGuides; + for (const auto& guide : guides) { + const bool isUnique = std::none_of(uniqueGuides.begin(), uniqueGuides.end(), [&](const Handle(Geom_BSplineCurve)& curve) { + return guide->IsEqual(curve, Precision::Confusion()); + }); + if (isUnique) { + uniqueGuides.push_back(guide); + } } - - m_profiles.reserve(profiles.size()); - m_guides.reserve(guides.size()); - // Copy the curves - for (std::vector::const_iterator it = profiles.begin(); it != profiles.end(); ++it) { - m_profiles.push_back(GeomConvert::CurveToBSplineCurve(*it)); + if (uniqueProfiles.size() < 2) { + throw Error("There must be at least two unique profiles for the curve network interpolation.", MATH_ERROR); } - for (std::vector::const_iterator it = guides.begin(); it != guides.end(); ++it) { - m_guides.push_back(GeomConvert::CurveToBSplineCurve(*it)); + + if (uniqueGuides.size() < 2) { + throw Error("There must be at least two unique guides for the curve network interpolation.", MATH_ERROR); } + + m_profiles = uniqueProfiles; + m_guides = uniqueGuides; } @@ -90,38 +119,13 @@ void InterpolateCurveNetwork::ComputeIntersections(math_Matrix& intersection_par } // for closed curves else if (currentIntersections.size() == 2) { - - // only the u-directional B-spline curves are closed - if (profiles[0]->IsClosed()) { - - if (spline_v_idx == 0) { - intersection_params_u(spline_u_idx, spline_v_idx) = std::min(currentIntersections[0].first, currentIntersections[1].first); - } - else if (spline_v_idx == static_cast(guides.size() - 1)) { - intersection_params_u(spline_u_idx, spline_v_idx) = std::max(currentIntersections[0].first, currentIntersections[1].first); - } - - // intersection_params_vector[0].second == intersection_params_vector[1].second - intersection_params_v(spline_u_idx, spline_v_idx) = currentIntersections[0].second; - } - - // only the v-directional B-spline curves are closed - if (guides[0]->IsClosed()) { - - if (spline_u_idx == 0) { - intersection_params_v(spline_u_idx, spline_v_idx) = std::min(currentIntersections[0].second, currentIntersections[1].second); - } - else if (spline_u_idx == static_cast(profiles.size() - 1)) { - intersection_params_v(spline_u_idx, spline_v_idx) = std::max(currentIntersections[0].second, currentIntersections[1].second); - } - // intersection_params_vector[0].first == intersection_params_vector[1].first - intersection_params_u(spline_u_idx, spline_v_idx) = currentIntersections[0].first; - } - -// // TODO: both u-directional splines and v-directional splines are closed -// else if (intersection_params_vector.size() == 4) { - -// } + // Closed curves produce the same intersection at both ends of + // their parameter range. Use the lower parameter here; the + // duplicated boundary curve is added after sorting. + intersection_params_u(spline_u_idx, spline_v_idx) = + std::min(currentIntersections[0].first, currentIntersections[1].first); + intersection_params_v(spline_u_idx, spline_v_idx) = + std::min(currentIntersections[0].second, currentIntersections[1].second); } else if (currentIntersections.size() > 2) { @@ -160,7 +164,6 @@ void InterpolateCurveNetwork::SortCurves(math_Matrix& intersection_params_u, mat void InterpolateCurveNetwork::MakeCurvesCompatible() { - // reparametrize into [0,1] for (CurveArray::iterator it = m_profiles.begin(); it != m_profiles.end(); ++it) { BSplineAlgorithms::reparametrizeBSpline(*(*it), 0., 1., 1e-15); @@ -173,18 +176,63 @@ void InterpolateCurveNetwork::MakeCurvesCompatible() int nGuides = static_cast(m_guides.size()); int nProfiles = static_cast(m_profiles.size()); - // now find all intersections of all B-splines with each other - math_Matrix intersection_params_u(0, nProfiles - 1, - 0, nGuides - 1); - math_Matrix intersection_params_v(0, nProfiles - 1, - 0, nGuides - 1); + math_Matrix tmp_intersection_params_u(0, nProfiles - 1, 0, nGuides - 1); + math_Matrix tmp_intersection_params_v(0, nProfiles - 1, 0, nGuides - 1); - ComputeIntersections(intersection_params_u, intersection_params_v); + ComputeIntersections(tmp_intersection_params_u, tmp_intersection_params_v); + SortCurves(tmp_intersection_params_u, tmp_intersection_params_v); - // sort intersection_params_u and intersection_params_v and u-directional and v-directional B-spline curves - SortCurves(intersection_params_u, intersection_params_v); + const bool isClosedProfile = m_profiles.front()->IsClosed() || m_profiles.front()->IsPeriodic(); + const bool isClosedGuides = m_guides.front()->IsClosed() || m_guides.front()->IsPeriodic(); + + if (isClosedProfile && isClosedGuides) { + throw Error("Closed in both U and V directions surface isn't supported at this time"); + } + + math_Matrix intersection_params_u(0, isClosedGuides ? nProfiles : nProfiles - 1, + 0, isClosedProfile ? nGuides : nGuides - 1); + math_Matrix intersection_params_v(0, isClosedGuides ? nProfiles : nProfiles - 1, + 0, isClosedProfile ? nGuides : nGuides - 1); + + if (isClosedProfile) { + m_guides.push_back(m_guides.front()); + ++nGuides; + + for (int spline_u_idx = 0; spline_u_idx < nProfiles; ++spline_u_idx) { + for (int spline_v_idx = 0; spline_v_idx < nGuides - 1; ++spline_v_idx) { + intersection_params_u(spline_u_idx, spline_v_idx) = tmp_intersection_params_u(spline_u_idx, spline_v_idx); + intersection_params_v(spline_u_idx, spline_v_idx) = tmp_intersection_params_v(spline_u_idx, spline_v_idx); + } + + intersection_params_u(spline_u_idx, nGuides - 1) = + tmp_intersection_params_u(spline_u_idx, 0) < 1e-5 + ? 1.0 + : tmp_intersection_params_u(spline_u_idx, 0); + intersection_params_v(spline_u_idx, nGuides - 1) = tmp_intersection_params_v(spline_u_idx, 0); + } + } + else if (isClosedGuides) { + m_profiles.push_back(m_profiles.front()); + ++nProfiles; + + for (int spline_v_idx = 0; spline_v_idx < nGuides; ++spline_v_idx) { + for (int spline_u_idx = 0; spline_u_idx < nProfiles - 1; ++spline_u_idx) { + intersection_params_u(spline_u_idx, spline_v_idx) = tmp_intersection_params_u(spline_u_idx, spline_v_idx); + intersection_params_v(spline_u_idx, spline_v_idx) = tmp_intersection_params_v(spline_u_idx, spline_v_idx); + } + + intersection_params_u(nProfiles - 1, spline_v_idx) = tmp_intersection_params_u(0, spline_v_idx); + intersection_params_v(nProfiles - 1, spline_v_idx) = + tmp_intersection_params_v(0, spline_v_idx) < 1e-5 + ? 1.0 + : tmp_intersection_params_v(0, spline_v_idx); + } + } + else { + intersection_params_u = tmp_intersection_params_u; + intersection_params_v = tmp_intersection_params_v; + } - // eliminate small inaccuracies of the intersection parameters: EliminateInaccuraciesNetworkIntersections(m_profiles, m_guides, intersection_params_u, intersection_params_v); std::vector newParametersProfiles; @@ -344,6 +392,27 @@ void InterpolateCurveNetwork::EliminateInaccuraciesNetworkIntersections(const st } } +void InterpolateCurveNetwork::EnsureC2() +{ + if (m_gordonSurf.IsNull()) { + return; + } + + const int minUMult = std::max(1, m_gordonSurf->UDegree() - 2); + for (int iu = 2; iu <= m_gordonSurf->NbUKnots() - 1; ++iu) { + if (m_gordonSurf->UMultiplicity(iu) > minUMult) { + m_gordonSurf->RemoveUKnot(iu, minUMult, m_spatialTol); + } + } + + const int minVMult = std::max(1, m_gordonSurf->VDegree() - 2); + for (int iv = 2; iv <= m_gordonSurf->NbVKnots() - 1; ++iv) { + if (m_gordonSurf->VMultiplicity(iv) > minVMult) { + m_gordonSurf->RemoveVKnot(iv, minVMult, m_spatialTol); + } + } +} + Handle(Geom_BSplineSurface) InterpolateCurveNetwork::Surface() { @@ -408,9 +477,18 @@ void InterpolateCurveNetwork::Perform() m_skinningSurfGuides = builder.SurfaceGuides(); m_tensorProdSurf = builder.SurfaceIntersections(); + EnsureC2(); + m_hasPerformed = true; } +Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector &profiles, + const std::vector &guides, + double tol) +{ + return InterpolateCurveNetwork(profiles, guides, tol).Surface(); +} + Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector &profiles, const std::vector &guides, double tol) { return InterpolateCurveNetwork(profiles, guides, tol).Surface(); diff --git a/src/geometry/curve-networks/InterpolateCurveNetwork.h b/src/geometry/curve-networks/InterpolateCurveNetwork.h index c2475f0..21e2f6c 100644 --- a/src/geometry/curve-networks/InterpolateCurveNetwork.h +++ b/src/geometry/curve-networks/InterpolateCurveNetwork.h @@ -54,6 +54,16 @@ class InterpolateCurveNetwork const std::vector& guides, double spatialTolerance); + /** + * @brief InterpolateCurveNetwork interpolated a curve network of already converted B-spline curves + * @param profiles The profiles to be interpolated + * @param guides The guides curves to be interpolated + * @param spatialTolerance Maximum allowed distance between each guide and profile (in theory they must intersect) + */ + GEOML_EXPORT InterpolateCurveNetwork(const std::vector& profiles, + const std::vector& guides, + double spatialTolerance); + GEOML_EXPORT operator Handle(Geom_BSplineSurface) (); /// Returns the interpolation surface @@ -91,6 +101,8 @@ class InterpolateCurveNetwork math_Matrix & intersection_params_u, math_Matrix & intersection_params_v) const; + void EnsureC2(); + bool m_hasPerformed; double m_spatialTol; @@ -106,6 +118,11 @@ GEOML_EXPORT Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector const std::vector& guides, double spatialTol = 3e-4); +/// Convenience function calling InterpolateCurveNetwork +GEOML_EXPORT Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector& profiles, + const std::vector& guides, + double spatialTol = 3e-4); + } // namespace geoml #endif // INTERPOLATECURVENETWORK_H diff --git a/src/geoml/surfaces/surfaces.cpp b/src/geoml/surfaces/surfaces.cpp index 1547c06..3ef4ecf 100644 --- a/src/geoml/surfaces/surfaces.cpp +++ b/src/geoml/surfaces/surfaces.cpp @@ -38,6 +38,16 @@ Handle(Geom_BSplineSurface) interpolate_curve_network(const std::vector &ucurves, + const std::vector &vcurves, + double tolerance) +{ + InterpolateCurveNetwork interpolator(ucurves, vcurves, tolerance); + + return interpolator.Surface(); +} + Handle(Geom_BSplineSurface) interpolate_curves(const std::vector &ucurves, unsigned int max_degree, bool join_continuously) diff --git a/src/geoml/surfaces/surfaces.h b/src/geoml/surfaces/surfaces.h index 7f5fc4a..50b9ba8 100644 --- a/src/geoml/surfaces/surfaces.h +++ b/src/geoml/surfaces/surfaces.h @@ -10,6 +10,7 @@ #include "geoml/data_structures/Array2d.h" #include +#include #include #include #include @@ -41,6 +42,16 @@ interpolate_curve_network(const std::vector& ucurves, const std::vector& vcurves, double tolerance); +/** + * @brief Interpolates the curve network by a B-spline surface + * + * Variant for already converted B-spline curves. + */ +GEOML_API_EXPORT Handle(Geom_BSplineSurface) +interpolate_curve_network(const std::vector& ucurves, + const std::vector& vcurves, + double tolerance); + /** * @brief Interpolates the curves by a B-spline surface in v-direction * diff --git a/src/math/IntersectBSplines.cpp b/src/math/IntersectBSplines.cpp index 7d5f2cb..6b3d4b9 100644 --- a/src/math/IntersectBSplines.cpp +++ b/src/math/IntersectBSplines.cpp @@ -222,8 +222,10 @@ namespace class CurveCurveDistanceObjective : public math_MultipleVarFunctionWithGradient { public: - CurveCurveDistanceObjective(const Handle(Geom_Curve)& c1, const Handle(Geom_Curve)& c2) - : m_c1(c1), m_c2(c2) + CurveCurveDistanceObjective(const Handle(Geom_Curve)& c1, const Handle(Geom_Curve)& c2, double scale) + : m_c1(c1) + , m_c2(c2) + , m_invScaleSqr(1. / sqr(std::max(scale, std::numeric_limits::epsilon()))) {} virtual Standard_Integer NbVariables() const override @@ -300,15 +302,16 @@ namespace m_c2->D1(v, p2, d2); gp_Vec diff = p1.XYZ() - p2.XYZ(); - F = diff.SquareMagnitude(); - G(1) = 2. * diff.Dot(d1)* (m_c1->LastParameter()-m_c1->FirstParameter()) * d_getUParam(X.Value(1)); - G(2) = -2. * diff.Dot(d2) * (m_c2->LastParameter()-m_c2->FirstParameter()) * d_getUParam(X.Value(2)); + F = diff.SquareMagnitude() * m_invScaleSqr; + G(1) = 2. * diff.Dot(d1) * d_getUParam(X.Value(1)) * m_invScaleSqr; + G(2) = -2. * diff.Dot(d2) * d_getVParam(X.Value(2)) * m_invScaleSqr; return true; } private: const Handle(Geom_Curve) m_c1, m_c2; + const double m_invScaleSqr; }; @@ -320,6 +323,7 @@ namespace geoml std::vector IntersectBSplines(const Handle(Geom_BSplineCurve) curve1, const Handle(Geom_BSplineCurve) curve2, double tolerance) { + const double optimizerScale = (geoml::BSplineAlgorithms::scale(curve1) + geoml::BSplineAlgorithms::scale(curve2)) / 2.; auto hulls = getRangesOfIntersection(curve1, curve2, tolerance); std::list curve1_ints, curve2_ints; @@ -372,7 +376,7 @@ std::vector IntersectBSplines(const Handle(Geom_ auto c1 = BSplineAlgorithms::trimCurve(curve1, boxes.b1.range.min,boxes.b1.range.max); auto c2 = BSplineAlgorithms::trimCurve(curve2, boxes.b2.range.min,boxes.b2.range.max); - CurveCurveDistanceObjective obj(c1, c2); + CurveCurveDistanceObjective obj(c1, c2, optimizerScale); // The objective is designed such that x=[0, 0] is in the middle of the parameter space of both curves math_Vector guess(1, 2); diff --git a/tests/unittests/testctiglbsplinealgorithms.cpp b/tests/unittests/testctiglbsplinealgorithms.cpp index 0354218..f077dba 100644 --- a/tests/unittests/testctiglbsplinealgorithms.cpp +++ b/tests/unittests/testctiglbsplinealgorithms.cpp @@ -1013,8 +1013,8 @@ TEST(BSplineAlgorithms, testIntersectionFinder) std::vector > intersection_vector = BSplineAlgorithms::intersections(spline_u, spline_v); // splines should intersect at u = 0.5 + std::sqrt(0.1) and v = 4. / 5 - ASSERT_NEAR(intersection_vector[0].first, 0.5 + std::sqrt(0.1), 1e-13); - ASSERT_NEAR(intersection_vector[0].second, 4. / 5, 1e-13); + ASSERT_NEAR(intersection_vector[0].first, 0.5 + std::sqrt(0.1), 1e-12); + ASSERT_NEAR(intersection_vector[0].second, 4. / 5, 1e-12); } /* TEST(BSplineAlgorithms, testSortBSpline) @@ -1232,9 +1232,9 @@ TEST(BSplineAlgorithms, testCreateGordonSurfaceGeneral) gp_Pnt point_curve2 = spline_u4->Value(u_value); // represents y(z) = (z - 0.5)^2 with offset 2 in x-direction gp_Pnt right_point(point_curve1.X() * (1. - v_value) + point_curve2.X() * v_value, point_curve1.Y() * (1. - v_value) + point_curve2.Y() * v_value, point_curve1.Z() * (1. - v_value) + point_curve2.Z() * v_value); - ASSERT_NEAR(surface_point.X(), right_point.X(), 1e-13); - ASSERT_NEAR(surface_point.Y(), right_point.Y(), 1e-13); - ASSERT_NEAR(surface_point.Z(), right_point.Z(), 1e-13); + ASSERT_NEAR(surface_point.X(), right_point.X(), 1e-12); + ASSERT_NEAR(surface_point.Y(), right_point.Y(), 1e-12); + ASSERT_NEAR(surface_point.Z(), right_point.Z(), 1e-12); } } }