Add functionality for reparameterizing B-splines without changing geometry (Fix issue #890)

src/fuselage/CCPACSFuselageProfile.cpp

@@ -203,7 +203,9 @@ void CCPACSFuselageProfile::BuildWires(WireCache& cache) const
         occPoints->SetValue(occPoints->Upper(), pMiddle);
     }
 
-    CTiglInterpolatePointsWithKinks interp(occPoints, kinks, params, 0.5, 3);
+    // Here, the B-spline is reparameterized based on the CPACS profile after setting it up at first for accuracy reasons
+    // Also, a tolerance is passed used for knot insertion and removal during the reparameterization algorithm
+    CTiglInterpolatePointsWithKinks interp(occPoints, kinks, params, 0.5, 3, CTiglInterpolatePointsWithKinks::Algo::InterpolateFirstThenReparametrize, 1e-8);
     auto spline = interp.Curve();
 
     if (mirrorSymmetry) {
@@ -213,9 +215,14 @@ void CCPACSFuselageProfile::BuildWires(WireCache& cache) const
         CTiglBSplineAlgorithms::reparametrizeBSpline(*spline, umin, umax);
     }
 
-    // we reparametrize the spline to get better performing lofts.
-    // there might be a small accuracy loss though.
-    spline = CTiglBSplineAlgorithms::reparametrizeBSplineNiceKnots(spline).curve;
+    // Reparamaterization based on a ParamMap defined in the CPACS file within CTiglInterpolatePointsWithKinks does not get along with reparametrizeBSplineNiceKnots.
+    // The geometry is changed in a way that is not acceptable anymore. However out of performance reasons, the feature should not be rejected in the most cases.
+    // Due to those conflicts, the function is only called, when there are no parameters defined in the CPACS file:
+    if (params.empty()) {
+        // we reparametrize the spline to get better performing lofts.
+        // there might be a small accuracy loss though.
+        spline = CTiglBSplineAlgorithms::reparametrizeBSplineNiceKnots(spline).curve;
+     }
 
     // Create wires
     TopoDS_Edge edge = BRepBuilderAPI_MakeEdge(spline).Edge();

src/geometry/CTiglBSplineAlgorithms.cpp

@@ -314,9 +314,147 @@ namespace
     {
         return array2GetRow<TColgp_Array2OfPnt, TColgp_HArray1OfPnt, Handle_TColgp_HArray1OfPnt>(matrix, rowIndex);
     }
-
+
+    // Some helper functions for CTiglBSplineAlgorithms::reparameterizePiecewiseLinear ...
+    void knotRefinementForReparam(const Handle(Geom_BSplineCurve) curve, std::vector<double> const& paramsOld, double tolerance)
+    {
+        // This function matches the first step of the reparameterization algorithm in The NURBS Book (2nd edition), p. 251
+
+        Standard_Integer nbKnots = curve->NbKnots();
+
+        //Start with increasing the multiplicity up to p
+        Standard_Integer degree = curve->Degree();
+
+        // Pass first and last index of knots (to increase all multiplicities until mult >= degree)
+        curve->IncreaseMultiplicity(1, nbKnots, degree);
+
+        // Add all knots u_i=f(s_i)
+        // Since we assume a piecewise linear reparameterization function f, the knots of f are just the new parameters
+        // Hence, evaluating f at one knot gives the corresponding old parameter
+        TColStd_Array1OfReal knotsParams(1, paramsOld.size());
+        TColStd_Array1OfInteger multsParams(1, paramsOld.size());
+        for (int paramIdx = 0; paramIdx < paramsOld.size(); paramIdx++) {
+            knotsParams.SetValue(paramIdx+1, paramsOld[paramIdx]);
+            multsParams.SetValue(paramIdx+1, degree);
+        }
+        curve->InsertKnots(knotsParams, multsParams, tolerance);
+    }
+
+    TColStd_Array1OfReal calcNewKnotVectorS(const Handle(Geom_BSplineCurve) curve, std::vector<double> const& paramsOld, std::vector<double> const& paramsNew, double tolerance)
+    {
+        // This function matches the second step of the reparameterization algorithm in The NURBS Book (2nd edition), p. 251
+
+        double newKnot;
+
+        // Set up knot vector S' consisting of s_i=g(u_i) [inverse reparametrization fct]
+        Standard_Integer nbKnots = curve->NbKnots(); // Size of distinct knots
+        TColStd_Array1OfReal newKnots(1, nbKnots);
+        for (int knotIdx = 1; knotIdx <= nbKnots; knotIdx++) {
+            // Calc s_i=g(u_i) of picewise linear fct
+            // defined as inverse of interpolation (x_i,y_i)=(paramsNew[i], paramsOld[i])
+            newKnot = details::calcReparamfctInv(paramsOld, paramsNew, curve->Knot(knotIdx), tolerance);
+            newKnots.SetValue(knotIdx, newKnot);
+        }
+        return newKnots;
+    }
+
+    Standard_Integer findKnotIndex(TColStd_Array1OfReal knots, double value)
+    {
+        for(int idx = 1; idx <= knots.Size(); idx++) {
+            if(fabs(knots[idx] - value) <= 1e-12)
+                return idx;
+        }
+        return -1;
+    }
+
+    void removeKnotsAfterReparam(const Handle(Geom_BSplineCurve) &curve,
+                                 const Handle(Geom_BSplineCurve) curveOrigin,
+                                 std::vector<double> const& paramsOld, std::vector<double> const& paramsNew,
+                                 double tolerance)
+    {
+        // This function matches the fourth step of the reparameterization algorithm in The NURBS Book (2nd edition), p. 251
+
+        TColStd_Array1OfReal knotsSBar = curve->Knots();
+
+        Standard_Integer nbKnots = curve->NbKnots();
+        double knotS, knotU;
+        Standard_Integer multS, multU, idxU, multNew;
+        // Exclude first and last knot
+        // Iterate backwards to account for shrinking size of knots vector by deletion if multNew is 0
+        for (int knotIdx = nbKnots-1; knotIdx >= 2; knotIdx--) {
+            knotS = knotsSBar[knotIdx];
+            multS = 0;
+            // Define multiplicity of knotS in original knot vector S of reparam fct
+            // By construction, the new parameters are exactely these knots and they all have mult=1
+            // => If knotS is found in these, its mult in this knot vector is exactely 1, 0 else
+            if(std::find_if(paramsNew.begin(), paramsNew.end(), [&knotS](double v){ return fabs(knotS-v) <= 1e-12; }) != paramsNew.end())
+                multS = 1;
+            // Used as original function, not inverse [swap interpolation dimensions]
+            knotU = details::calcReparamfctInv(paramsNew, paramsOld, knotS, tolerance);
+            idxU = findKnotIndex(curveOrigin->Knots(), knotU);
+            if(idxU != -1)
+                multU = curveOrigin->Multiplicity(idxU);
+            else
+                multU = 0;
+
+            //multNew = max(pq - p + multU, pq - q + multS) = max(multU, multS) [since q=1 and pq=p]
+            multNew = std::max(multU, multS);
+            // Decrease mult of knot to multNew
+            curve->RemoveKnot(knotIdx, multNew, tolerance);
+        }
+    }
+
 } // namespace
 
+namespace details
+{
+    double calcReparamfctInv(std::vector<double> const& paramsOld, std::vector<double> const& paramsNew, double u, double tolerance)
+    {
+        // Calc inverse of reparametrization function s=g(u)=f⁻¹(s)
+        // We set up BSpline of degree 1 -> g is piecewise linear
+        // In the interval [p_i,p_{i+1}], g(u)=m*u+n
+        // p_i and p_{i+1} are old parameters that need to be determined in the first place
+
+        // u needs to be in range of paramsOld
+        if (u < paramsOld[0] || paramsOld[paramsOld.size()-1] < u) {
+            if (paramsOld[0] - tolerance <= u && u <= paramsOld[paramsOld.size()-1] + tolerance) {
+                u = Clamp(u, paramsOld[0], paramsOld[paramsOld.size()-1]);
+            }
+            else {
+                throw tigl::CTiglError("Argument of reparameterization function out of range in calcReparamfctInv");
+            }
+        }
+
+        // First: Find interval [p_i,p_{i+1}] to determine boundaries
+        // a_x = lowerBound, b_x = upperBound
+
+        double a_x, a_y, b_x, b_y;
+        auto it = std::find_if(
+            paramsOld.begin(),
+            paramsOld.end(),
+            [&u](double v){ return u <= v; }
+        );
+
+        int paramsIdx = std::distance(paramsOld.begin(), it);
+        // Catch case u=0, since distance would be 0 which causes illegal vector access by subtraction of 1
+        if (paramsIdx  == 0) {
+            ++paramsIdx;
+        }
+
+        a_x = paramsOld[paramsIdx - 1];
+        b_x = paramsOld[paramsIdx];
+        a_y = paramsNew[paramsIdx - 1];
+        b_y = paramsNew[paramsIdx];
+
+        // Calc slope m
+        double m = (b_y-a_y) / (b_x-a_x);
+        // Calc intercept n
+        double n = a_y - a_x*m;
+
+        // Calc inverse g(u) [linear function value]
+        return m*u + n;
+    }
+} // namespace details
 
 namespace tigl
 {
@@ -836,6 +974,46 @@ CTiglApproxResult CTiglBSplineAlgorithms::reparametrizeBSplineNiceKnots(Handle(G
     return CTiglBSplineAlgorithms::reparametrizeBSplineContinuouslyApprox(spline, {umin, umax}, {umin, umax}, nSegmentsNew + target_degree);
 }
 
+Handle(Geom_BSplineCurve) CTiglBSplineAlgorithms::reparameterizePiecewiseLinear(Handle(Geom_BSplineCurve) curve,
+                                                                                std::vector<double> const& paramsOld,
+                                                                                std::vector<double> const& paramsNew,
+                                                                                double tolerance)
+{
+    // Apply reparameterization on a given B-Spline curve defined by old and new parameters
+    // Assumption:  - Reparameterization function given as 1D-interpolation of (x=paramsNew, y=paramsOld)
+    //              - Picewise linear function -> B-Spline of degree 1 (q=1 => pq=1)
+
+    // Check inputs
+    if (paramsOld.size() != paramsNew.size())
+        throw CTiglError("Parameter sizes for reparameterization do not match in reparameterizePiecewiseLinear");
+
+    if (tolerance < 0)
+        throw CTiglError("Tolerance for reparameterization has to be non-negative in reparameterizePiecewiseLinear");
+
+    const Standard_Integer degree = curve->Degree();
+
+    // Create copy as its original information (e.g. poles) is needed later
+    Handle(Geom_BSplineCurve) curveCopy = Handle(Geom_BSplineCurve)::DownCast(curve->Copy());
+
+    // Step 1
+    knotRefinementForReparam(curve, paramsOld, tolerance);
+
+    // Step 2
+    TColStd_Array1OfReal newKnots = calcNewKnotVectorS(curve, paramsOld, paramsNew, tolerance);
+
+    // Step 3 is not necessary as all poles remain the same (NURBS book, p. 250)
+    // -> Set up B-Spline from poles resulting from first knot refinement (step 1) and the new knots (and mults)
+
+    // Internal knots should appear with multiplicity p*q according to the NURBS book, boundary knots have multiplicity p+1 (step 2)
+    // For that reason, the multiplicities of the curve after knot refinement are used when setting up the curve later
+    Handle(Geom_BSplineCurve) curveReparameterized = new Geom_BSplineCurve(curve->Poles(), newKnots, curve->Multiplicities(), degree);
+
+    // Step 4
+    removeKnotsAfterReparam(curveReparameterized, curveCopy, paramsOld, paramsNew, tolerance);
+
+    return curveReparameterized;
+}
+
 math_Matrix CTiglBSplineAlgorithms::bsplineBasisMat(int degree, const TColStd_Array1OfReal& knots, const TColStd_Array1OfReal& params, unsigned int derivOrder)
 {
     Standard_Integer ncp = knots.Length() - degree - 1;

src/geometry/CTiglBSplineAlgorithms.h

@@ -158,6 +158,24 @@ class CTiglBSplineAlgorithms
      */
     TIGL_EXPORT static CTiglApproxResult reparametrizeBSplineNiceKnots(Handle(Geom_BSplineCurve) spline);
 
+    /**
+     * @brief reparameterizePiecewiseLinear:
+     *          Apply reparameterization on a given B-Spline curve defined by old and new parameters
+     *          Based on algorithm found in The NURBS book (2nd edition), p. 251, and the explanations
+     *          Here, we use a piecewise linear reparameterization function (degree q=1) interpolating the wanted parameters
+     *          As a result, the degree stays the same after reparameterization
+     * @param paramsOld:
+     *          Array of the old parameters
+     * @param paramsNew:
+     *          Array of the new parameters
+     * @param tolerance:
+     *          Define the tolerance used for OpenCascade knot removal funtion [Geom_BSplineCurve::RemoveKnot(Index, M, Tolerance)]
+     */
+    TIGL_EXPORT static Handle(Geom_BSplineCurve) reparameterizePiecewiseLinear(Handle(Geom_BSplineCurve) curve,
+                                                                               std::vector<double> const& paramsOld, std::vector<double> const& paramsNew,
+                                                                               double tolerance);
+
+
     /**
      * @brief reparametrizeBSplineContinuouslyApprox:
      *          Reparametrizes a given B-spline by giving an array of its old parameters that should have the values of the given array of new parameters after this function call.
@@ -318,4 +336,21 @@ class CTiglBSplineAlgorithms
 };
 } // namespace tigl
 
+namespace details
+{
+    /// Helper function for reparameterization to calculate inverse of reparameterization function (paramsOld -> paramsNew)
+    /**
+     * @brief calcReparamfctInv
+     * @param paramsOld:
+     *          Parameters of original B-Spline curve (domain of inverse reparameterization function)
+     * @param paramsNew:
+     *          New parameters of B-Spline curve (range of inverse reparameterization function)
+     * @param u:
+     *          Function argument, has to be inside range of paramsOld
+     * @param tolerance
+     * @return the wanted inverse of u (new parameter corresponding to old parameter)
+     */
+    TIGL_EXPORT double calcReparamfctInv(std::vector<double> const& paramsOld, std::vector<double> const& paramsNew, double u, double tolerance);
+} // namespace details
+
 #endif // CTIGLBSPLINEALGORITHMS_H

src/geometry/CTiglInterpolatePointsWithKinks.cpp

@@ -127,15 +127,19 @@ namespace tigl
 {
 
 CTiglInterpolatePointsWithKinks::CTiglInterpolatePointsWithKinks(const Handle(TColgp_HArray1OfPnt) & points,
-                                                                const std::vector<unsigned int>& kinkIndices,
-                                                                const ParamMap& parameters,
-                                                                double alpha,
-                                                                unsigned int maxDegree)
+                                                                 const std::vector<unsigned int>& kinkIndices,
+                                                                 const ParamMap& parameters,
+                                                                 double alpha,
+                                                                 unsigned int maxDegree,
+                                                                 Algo algo,
+                                                                 const double tolerance)
     : m_pnts(points)
     , m_kinks(kinkIndices)
     , m_params(parameters)
     , m_alpha(alpha)
     , m_maxDegree(maxDegree)
+    , m_algo(algo)
+    , m_tolerance(tolerance)
     , m_result(*this, &CTiglInterpolatePointsWithKinks::ComputeResult)
 {
 }
@@ -152,7 +156,14 @@ std::vector<double> CTiglInterpolatePointsWithKinks::Parameters() const
 
 void CTiglInterpolatePointsWithKinks::ComputeResult(Result& result) const
 {
-    auto params = m_params;
+
+    ParamMap params;
+    // m_params should only be used for Algo::InterpolateBasedOnParameters
+    // When Algo::InterpolateFirstThenReparametrize is used, the parameterization should be applied after interpolating
+    // So, params stay empty to ensure TiGL uses its default parameters
+    if (m_algo == Algo::InterpolateBasedOnParameters) {
+        params = m_params;
+    }
 
     // assuming iterating over parameters is sorted in keys
     auto new_params = computeParams(m_pnts, params, m_alpha);
@@ -229,6 +240,19 @@ void CTiglInterpolatePointsWithKinks::ComputeResult(Result& result) const
 
     result.curve = CTiglBSplineAlgorithms::concatCurves(curve_segments, false);
     result.parameters = new_params;
+
+    if (m_algo == Algo::InterpolateFirstThenReparametrize) {
+        // We reparametrize the spline to get better performing lofts
+        const double tolerance = m_tolerance; // Define tolerance for knot removal in reparameterizePiecewiseLinear
+        // Get the above calculated params (Default params calculated without setting them explicitely) => paramsOld
+        // Wanted parameters [m_params] appear as map -> transform and interpolate to all m_pnts => paramsNew
+        std::vector<double> paramsOld = new_params;
+        params = m_params;
+        std::vector<double> paramsNew = computeParams(m_pnts, params, 0.5);
+        result.curve = CTiglBSplineAlgorithms::reparameterizePiecewiseLinear(result.curve, paramsOld, paramsNew, tolerance);
+        // Overwrite the object's variable with the parameters that are set after applying the reparameterization
+        result.parameters = paramsNew;
+    }
 }
 
 

src/geometry/CTiglInterpolatePointsWithKinks.h

@@ -41,11 +41,18 @@ using ParamMap = std::map<unsigned int, double>;
 class CTiglInterpolatePointsWithKinks
 {
 public:
+    enum class Algo {
+        InterpolateBasedOnParameters,       /**< generates an interpolating curve that matches the parameter map directly. The geometry of the curve is directly influenced by the parameter map */
+        InterpolateFirstThenReparametrize,  /**< interpolates first and then reparametrizes to match the parameter map. The geometry of the curve is independent of the parameter map*/
+    };
+
     TIGL_EXPORT CTiglInterpolatePointsWithKinks(const Handle(TColgp_HArray1OfPnt) & points,
                                                 const std::vector<unsigned int>& kinkIndices,
                                                 const ParamMap& parameters,
                                                 double alpha = 0.5,
-                                                unsigned int maxDegree=3);
+                                                unsigned int maxDegree=3,
+                                                Algo algo = Algo::InterpolateBasedOnParameters,
+                                                const double tolerance = 1e-8);
 
 
     TIGL_EXPORT Handle(Geom_BSplineCurve) Curve() const;
@@ -58,6 +65,8 @@ class CTiglInterpolatePointsWithKinks
     ParamMap m_params;
     double m_alpha;
     unsigned int m_maxDegree;
+    Algo m_algo;
+    const double m_tolerance;
 
     struct Result {
         Handle(Geom_BSplineCurve) curve;