Backport algorithmic changes from occ-gordon

src/geometry/BSplineAlgorithms.cpp

@@ -435,6 +435,17 @@ std::vector<double> BSplineAlgorithms::knotsFromCurveParameters(std::vector<doub
     return knots;
 }
 
+std::vector<Handle(Geom_BSplineCurve)> BSplineAlgorithms::toBSplines(const std::vector<Handle(Geom_Curve)>& curves)
+{
+    std::vector<Handle(Geom_BSplineCurve)> 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());
 

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<Handle(Geom_BSplineCurve)> toBSplines(const std::vector<Handle(Geom_Curve)>& curves);
+
     /// Concatenates a list of bspline curves
     GEOML_EXPORT static Handle(Geom_BSplineCurve) concatCurves(std::vector<Handle(Geom_BSplineCurve)> curves,
                                                               bool parByLength=true, double tolerance = 1e-6);

src/geometry/curve-networks/InterpolateCurveNetwork.cpp

@@ -29,6 +29,7 @@
 #include <algorithm>
 
 #include <math_Matrix.hxx>
+#include <Precision.hxx>
 #include <TColStd_HArray1OfReal.hxx>
 #include <GeomConvert.hxx>
 
@@ -38,28 +39,56 @@ namespace geoml
 InterpolateCurveNetwork::InterpolateCurveNetwork(const std::vector<Handle (Geom_Curve)> &profiles,
                                                            const std::vector<Handle (Geom_Curve)> &guides,
                                                            double spatialTol)
+    : InterpolateCurveNetwork(BSplineAlgorithms::toBSplines(profiles),
+                              BSplineAlgorithms::toBSplines(guides),
+                              spatialTol)
+{
+}
+
+InterpolateCurveNetwork::InterpolateCurveNetwork(const std::vector<Handle(Geom_BSplineCurve)>& profiles,
+                                                 const std::vector<Handle(Geom_BSplineCurve)>& 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<Handle(Geom_BSplineCurve)> 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<Handle(Geom_BSplineCurve)> 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<Handle (Geom_Curve)>::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<Handle (Geom_Curve)>::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<int>(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<int>(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<int>(m_guides.size());
     int nProfiles = static_cast<int>(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<double> 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<Handle (Geom_BSplineCurve)> &profiles,
+                                                  const std::vector<Handle (Geom_BSplineCurve)> &guides,
+                                                  double tol)
+{
+    return InterpolateCurveNetwork(profiles, guides, tol).Surface();
+}
+
 Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector<Handle (Geom_Curve)> &profiles, const std::vector<Handle (Geom_Curve)> &guides, double tol)
 {
     return InterpolateCurveNetwork(profiles, guides, tol).Surface();

src/geometry/curve-networks/InterpolateCurveNetwork.h

@@ -54,6 +54,16 @@ class InterpolateCurveNetwork
                                              const std::vector<Handle(Geom_Curve)>& 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<Handle(Geom_BSplineCurve)>& profiles,
+                                         const std::vector<Handle(Geom_BSplineCurve)>& 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<Handle(Geom_Curve)>& guides,
                                                               double spatialTol = 3e-4);
 
+/// Convenience function calling InterpolateCurveNetwork
+GEOML_EXPORT Handle(Geom_BSplineSurface) curveNetworkToSurface(const std::vector<Handle(Geom_BSplineCurve)>& profiles,
+                                                              const std::vector<Handle(Geom_BSplineCurve)>& guides,
+                                                              double spatialTol = 3e-4);
+
 } // namespace geoml
 
 #endif // INTERPOLATECURVENETWORK_H

src/geoml/surfaces/surfaces.cpp

@@ -38,6 +38,16 @@ Handle(Geom_BSplineSurface) interpolate_curve_network(const std::vector<Handle (
     return interpolator.Surface();
 }
 
+Handle(Geom_BSplineSurface)
+interpolate_curve_network(const std::vector<Handle(Geom_BSplineCurve)> &ucurves,
+                          const std::vector<Handle(Geom_BSplineCurve)> &vcurves,
+                          double tolerance)
+{
+    InterpolateCurveNetwork interpolator(ucurves, vcurves, tolerance);
+
+    return interpolator.Surface();
+}
+
 Handle(Geom_BSplineSurface)
 interpolate_curves(const std::vector<Handle (Geom_Curve)> &ucurves, unsigned int max_degree,
                           bool join_continuously)

src/geoml/surfaces/surfaces.h

@@ -10,6 +10,7 @@
 #include "geoml/data_structures/Array2d.h"
 
 #include <Geom_BSplineSurface.hxx>
+#include <Geom_BSplineCurve.hxx>
 #include <Geom_Curve.hxx>
 #include <gp_Vec.hxx>
 #include <gp_Pnt.hxx>
@@ -41,6 +42,16 @@ interpolate_curve_network(const std::vector<Handle(Geom_Curve)>& ucurves,
                           const std::vector<Handle(Geom_Curve)>& 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<Handle(Geom_BSplineCurve)>& ucurves,
+                          const std::vector<Handle(Geom_BSplineCurve)>& vcurves,
+                          double tolerance);
+
 /**
  * @brief Interpolates the curves by a B-spline surface in v-direction
  *