Merge pull request #11 from rainman110/fix_scale_sensitivity

Fix scale sensitivity

Author: rainman110

src/internal/BSplineAlgorithms.cpp

@@ -754,8 +754,8 @@ ApproxResult BSplineAlgorithms::reparametrizeBSplineContinuouslyApprox(const Han
         approximationObj.InterpolatePoint(idx, true);
     }
 #endif
-    
-    ApproxResult result = approximationObj.FitCurveOptimal(parameters);
+
+    ApproxResult result = approximationObj.FitCurveOptimal(parameters, 1);
 
     assert(!result.curve.IsNull());
 

src/internal/BSplineApproxInterp.cpp

@@ -254,7 +254,7 @@ ApproxResult BSplineApproxInterp::FitCurveOptimal(const std::vector<double>& ini
     ApproxResult result = solve(parms, occKnots->Array1(), occMults->Array1());
     double old_error = result.error * 2.;
 
-    while(result.error > 0 && (old_error - result.error) / std::max(result.error, 1e-6) > 1e-3 && iteration < maxIter) {
+    while(result.error > 0 && (old_error - result.error) / result.error > 1e-3 && iteration < maxIter) {
         old_error = result.error;
 
         optimizeParameters(result.curve, parms);

src/internal/IntersectBSplines.cpp

@@ -232,8 +232,8 @@ 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<double>::epsilon())))
         {}
 
         virtual Standard_Integer NbVariables()  const override
@@ -310,18 +310,18 @@ namespace
             m_c2->D1(v, p2, d2);
 
             gp_Vec diff = p1.XYZ() - p2.XYZ();
-            F = diff.SquareMagnitude();
-            G(1) = 2. * diff.Dot(d1) * d_getUParam(X.Value(1));
-            G(2) = -2. * diff.Dot(d2)  * d_getVParam(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;
     };
 
-
     void CheckGradient(math_MultipleVarFunctionWithGradient& func, const math_Vector& X, double step)
     {
         int nvars = func.NbVariables();
@@ -373,6 +373,7 @@ namespace occ_gordon_internal
 
 std::vector<CurveIntersectionResult> IntersectBSplines(const Handle(Geom_BSplineCurve) curve1, const Handle(Geom_BSplineCurve) curve2, double tolerance)
 {
+    const double optimizerScale = (BSplineAlgorithms::scale(curve1) + BSplineAlgorithms::scale(curve2)) / 2.;
     auto hulls = getRangesOfIntersection(curve1, curve2, tolerance);
 
     std::list<BoundingBox> curve1_ints, curve2_ints;
@@ -425,7 +426,7 @@ std::vector<CurveIntersectionResult> IntersectBSplines(const Handle(Geom_BSpline
         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);

src/internal/IntersectBSplines.h

@@ -33,7 +33,9 @@ struct CurveIntersectionResult
  * intersect each other. If so, this process is repeated until the curve segment is almost a line segment.
  * This result is used to locally optimize into a true minimum.
  */
-std::vector<CurveIntersectionResult> IntersectBSplines(const Handle(Geom_BSplineCurve) curve1, const Handle(Geom_BSplineCurve) curve2, double tolerance=1e-5);
+std::vector<CurveIntersectionResult> IntersectBSplines(const Handle(Geom_BSplineCurve) curve1,
+                                                       const Handle(Geom_BSplineCurve) curve2,
+                                                       double absTolerance=1e-5);
 
 } // namespace occ_gordon_internal
 

tests/unittests/testctiglbsplinealgorithms.cpp

@@ -1004,7 +1004,7 @@ TEST(BSplineAlgorithms, testIntersectionFinder)
     std::vector<std::pair<double, double> > 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].first, 0.5 + std::sqrt(0.1), 1e-12);
     ASSERT_NEAR(intersection_vector[0].second, 4. / 5, 1e-13);
 }
 /*
@@ -1223,9 +1223,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);
         }
     }
 }