From: aml Date: Sun, 7 Feb 2016 06:31:58 +0000 (+0300) Subject: 0027112: GeomAPI_Interpolate produces wrong result X-Git-Url: http://git.dev.opencascade.org/gitweb/?a=commitdiff_plain;h=e5976dde58caa53059d33f05f7e72e4dd9ba4908;p=occt-copy.git 0027112: GeomAPI_Interpolate produces wrong result Support of the approximation of natural boundary condition (f''= 0) --- diff --git a/src/Geom2dAPI/Geom2dAPI_Interpolate.cxx b/src/Geom2dAPI/Geom2dAPI_Interpolate.cxx index 6b62098514..d8f71b78cb 100644 --- a/src/Geom2dAPI/Geom2dAPI_Interpolate.cxx +++ b/src/Geom2dAPI/Geom2dAPI_Interpolate.cxx @@ -127,114 +127,149 @@ static void BuildParameters(const Standard_Boolean PeriodicFlag, ParametersPtr->Value(ii) + distance) ; } } + //======================================================================= //function : BuildPeriodicTangents //purpose : //======================================================================= - -static void BuildPeriodicTangent( - const TColgp_Array1OfPnt2d& PointsArray, - TColgp_Array1OfVec2d& TangentsArray, - TColStd_Array1OfBoolean& TangentFlags, - const TColStd_Array1OfReal& ParametersArray) +static void BuildPeriodicTangent(const TColgp_Array1OfPnt2d& PointsArray, + TColgp_Array1OfVec2d& TangentsArray, + TColStd_Array1OfBoolean& TangentFlags, + const TColStd_Array1OfReal& ParametersArray, + const Interpolate_BCType theBCType) { - Standard_Integer - ii, - degree ; - Standard_Real *point_array, - *parameter_array, - eval_result[2][2] ; - - gp_Vec2d a_vector ; - - if (PointsArray.Length() < 3) { - Standard_ConstructionError::Raise(); - } - - if (!TangentFlags.Value(1)) { - degree = 3 ; - if (PointsArray.Length() == 3) { - degree = 2 ; - } - point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ; - parameter_array = - (Standard_Real *) &ParametersArray.Value(1) ; - TangentFlags.SetValue(1,Standard_True) ; - PLib::EvalLagrange(ParametersArray.Value(1), - 1, - degree, - 2, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 2 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(1,a_vector) ; + gp_Vec2d aCurTangent; + Standard_Integer degree = 3; + + if (PointsArray.Length() < 3) + Standard_ConstructionError::Raise(); + if (PointsArray.Length() == 3 && + theBCType == Interpolate_CLAMPED) + { + degree = 2; + } + + if (TangentFlags.Value(1)) + return; // yet computed. + else + TangentFlags.SetValue(1,Standard_True); + + if (theBCType == Interpolate_NATURAL) + { + // Compute tangent in second point. + aCurTangent = PLib:: + ComputeLagrangeTangent2d(PointsArray, ParametersArray, PointsArray.Lower() + 1); + + // Compute tangent in the first point to satisfy + // f''(0) = 0 condition approximation: + // This is approximation due to usage of the tangent approximation in the second point + // instead of the equation below directly in solver. + // + // Formula from the Golovanov book "Geometrical modeling", 2011, pp. 18: + // + // q0 = 1.5 * (p1 - p0) / (t1 - t0) - 0.5 q1 + // + // q1=p0*(t1-t2)/((t0-t1)(t0-t2)) + p1 (2t1-t0-t2)/((t1-t0)(t1-t2))+ p2(t1-t0)/((t2-t0)(t2-t1)) + // + // Where: + // t - parameter + // p - point + // q - tangent vectors + gp_Vec2d aVec(PointsArray(1), PointsArray(2)); + aCurTangent = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent; } + else if (theBCType == Interpolate_CLAMPED) + { + Standard_Real *point_array, *parameter_array, eval_result[2][2]; + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()); + parameter_array = (Standard_Real *) &ParametersArray.Value(1); + + PLib::EvalLagrange(ParametersArray.Value(1), 1, degree, 2, point_array[0], + parameter_array[0], eval_result[0][0]); + for (Standard_Integer ii = 1 ; ii <= 2 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + } + + TangentsArray.SetValue(1,aCurTangent); } //======================================================================= //function : BuildTangents //purpose : //======================================================================= - static void BuildTangents(const TColgp_Array1OfPnt2d& PointsArray, - TColgp_Array1OfVec2d& TangentsArray, - TColStd_Array1OfBoolean& TangentFlags, - const TColStd_Array1OfReal& ParametersArray) + TColgp_Array1OfVec2d& TangentsArray, + TColStd_Array1OfBoolean& TangentFlags, + const TColStd_Array1OfReal& ParametersArray, + const Interpolate_BCType theBCType) { - Standard_Integer ii, - degree ; - Standard_Real *point_array, - *parameter_array, - - eval_result[2][2] ; - gp_Vec2d a_vector ; - - degree = 3 ; + gp_Vec2d aCurTangent; + Standard_Integer ii, degree = 3; + if (PointsArray.Length() < 3) + Standard_ConstructionError::Raise(); + if (PointsArray.Length() == 3 && + theBCType == Interpolate_CLAMPED) + { + degree = 2; + } - if ( PointsArray.Length() < 3) { - Standard_ConstructionError::Raise(); - } - if (PointsArray.Length() == 3) { - degree = 2 ; - } - if (!TangentFlags.Value(1)) { - point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ; - parameter_array = - (Standard_Real *) &ParametersArray.Value(1) ; - TangentFlags.SetValue(1,Standard_True) ; - PLib::EvalLagrange(ParametersArray.Value(1), - 1, - degree, - 2, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 2 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(1,a_vector) ; - } - if (! TangentFlags.Value(TangentFlags.Upper())) { - point_array = - (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree) ; - TangentFlags.SetValue(TangentFlags.Upper(),Standard_True) ; - parameter_array = - (Standard_Real *)&ParametersArray.Value(ParametersArray.Upper() - degree) ; - PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()), - 1, - degree, - 2, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 2 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(TangentsArray.Upper(),a_vector) ; - } -} + if (TangentFlags.Value(1) && + TangentFlags.Value(TangentFlags.Upper())) + { + return; + } + + if (theBCType == Interpolate_NATURAL) + { + // The formulas for tangent computation is stored in BuildPeriodicTangent. + if (!TangentFlags.Value(1)) + { + // Compute tangent in second point. + aCurTangent = PLib:: + ComputeLagrangeTangent2d(PointsArray, ParametersArray, PointsArray.Lower() + 1); + + gp_Vec2d aVec(PointsArray(1), PointsArray(2)); + aVec = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent; + TangentFlags.SetValue(1, Standard_True); + TangentsArray.SetValue(1, aVec); + } // if (!TangentFlags.Value(1)) + if (!TangentFlags.Value(TangentFlags.Upper())) + { + // Compute tangent in last but one point. + aCurTangent = PLib:: + ComputeLagrangeTangent2d(PointsArray, ParametersArray, PointsArray.Upper() - 1); + + Standard_Integer aLastIdx = PointsArray.Upper(); + gp_Vec2d aVec(PointsArray(aLastIdx - 1), PointsArray(aLastIdx)); + aVec = 1.5 * aVec / (ParametersArray(aLastIdx) - ParametersArray(aLastIdx - 1)) - 0.5 * aCurTangent; + TangentFlags.SetValue(TangentFlags.Upper(), Standard_True); + TangentsArray.SetValue(TangentsArray.Upper(), aVec); + } + } + else if (theBCType == Interpolate_CLAMPED) + { + Standard_Real *point_array, *parameter_array, eval_result[2][2]; + if (!TangentFlags.Value(1)) + { + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()); + parameter_array = (Standard_Real *) &ParametersArray.Value(1); + PLib::EvalLagrange(ParametersArray.Value(1), 1, degree, 2, point_array[0], parameter_array[0], eval_result[0][0]) ; + for (ii = 1 ; ii <= 2 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + TangentFlags.SetValue(1, Standard_True); + TangentsArray.SetValue(1, aCurTangent); + } // if (!TangentFlags.Value(1)) + if (!TangentFlags.Value(TangentFlags.Upper())) + { + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree); + parameter_array = (Standard_Real *) &ParametersArray.Value(ParametersArray.Upper() - degree); + PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()), 1, degree, 2,point_array[0], parameter_array[0], eval_result[0][0]); + for (ii = 1 ; ii <= 2 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + TangentFlags.SetValue(TangentFlags.Upper(), Standard_True); + TangentsArray.SetValue(TangentsArray.Upper(), aCurTangent); + } + } +} //======================================================================= //function : BuildTangents //purpose : scale the given tangent so that they have the length of @@ -310,41 +345,28 @@ static void ScaleTangents(const TColgp_Array1OfPnt2d& PointsArray, //purpose : //======================================================================= -Geom2dAPI_Interpolate::Geom2dAPI_Interpolate - (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr, - const Standard_Boolean PeriodicFlag, - const Standard_Real Tolerance) : -myTolerance(Tolerance), -myPoints(PointsPtr), -myIsDone(Standard_False), -myPeriodic(PeriodicFlag), -myTangentRequest(Standard_False) - +Geom2dAPI_Interpolate::Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d)& PointsPtr, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType) +: myTolerance(Tolerance), + myPoints(PointsPtr), + myIsDone(Standard_False), + myPeriodic(PeriodicFlag), + myTangentRequest(Standard_False), + myBCType(theBCType) { - Standard_Integer ii ; - Standard_Boolean result = - CheckPoints(PointsPtr->Array1(), - Tolerance) ; - myTangents = - new TColgp_HArray1OfVec2d(myPoints->Lower(), - myPoints->Upper()) ; - myTangentFlags = - new TColStd_HArray1OfBoolean(myPoints->Lower(), - myPoints->Upper()) ; - - if (!result) { - Standard_ConstructionError::Raise(); - } - BuildParameters(PeriodicFlag, - PointsPtr->Array1(), - myParameters) ; + Standard_Boolean result = CheckPoints(PointsPtr->Array1(), Tolerance); + myTangents = new TColgp_HArray1OfVec2d(myPoints->Lower(), myPoints->Upper()); + myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(), myPoints->Upper()); - for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) { - myTangentFlags->SetValue(ii,Standard_False) ; - } + if (!result) + Standard_ConstructionError::Raise(); - - + BuildParameters(PeriodicFlag, PointsPtr->Array1(), myParameters); + + for (Standard_Integer ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) + myTangentFlags->SetValue(ii,Standard_False); } //======================================================================= @@ -352,51 +374,40 @@ myTangentRequest(Standard_False) //purpose : //======================================================================= -Geom2dAPI_Interpolate::Geom2dAPI_Interpolate - (const Handle(TColgp_HArray1OfPnt2d)& PointsPtr, - const Handle(TColStd_HArray1OfReal)& ParametersPtr, - const Standard_Boolean PeriodicFlag, - const Standard_Real Tolerance) : -myTolerance(Tolerance), -myPoints(PointsPtr), -myIsDone(Standard_False), -myParameters(ParametersPtr), -myPeriodic(PeriodicFlag), -myTangentRequest(Standard_False) +Geom2dAPI_Interpolate::Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d) &PointsPtr, + const Handle(TColStd_HArray1OfReal) &ParametersPtr, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType) +: myTolerance(Tolerance), + myPoints(PointsPtr), + myIsDone(Standard_False), + myParameters(ParametersPtr), + myPeriodic(PeriodicFlag), + myTangentRequest(Standard_False), + myBCType(theBCType) { - Standard_Integer ii ; - - - Standard_Boolean result = - CheckPoints(PointsPtr->Array1(), - Tolerance) ; + Standard_Boolean result = CheckPoints(PointsPtr->Array1(), Tolerance); - if (PeriodicFlag) { - if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) { - Standard_ConstructionError::Raise(); - } - } - myTangents = - new TColgp_HArray1OfVec2d(myPoints->Lower(), - myPoints->Upper()) ; - myTangentFlags = - new TColStd_HArray1OfBoolean(myPoints->Lower(), - myPoints->Upper()) ; - - if (!result) { - Standard_ConstructionError::Raise(); - } - - result = - CheckParameters(ParametersPtr->Array1()) ; - if (!result) { - Standard_ConstructionError::Raise(); - } - - for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) { - myTangentFlags->SetValue(ii,Standard_False) ; - } - + if (PeriodicFlag) + { + if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) + { + Standard_ConstructionError::Raise(); + } + } + myTangents = new TColgp_HArray1OfVec2d(myPoints->Lower(), myPoints->Upper()); + myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(), myPoints->Upper()); + + if (!result) + Standard_ConstructionError::Raise(); + + result = CheckParameters(ParametersPtr->Array1()); + if (!result) + Standard_ConstructionError::Raise(); + + for (Standard_Integer ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) + myTangentFlags->SetValue(ii,Standard_False); } //======================================================================= //function : Load @@ -566,14 +577,15 @@ void Geom2dAPI_Interpolate::PerformPeriodic() mults.SetValue(num_distinct_knots ,half_order) ; if (num_points >= 3) { -// -// only enter here if there are more than 3 points otherwise -// it means we have already the tangent -// + // + // only enter here if there are more than 3 points otherwise + // it means we have already the tangent + // BuildPeriodicTangent(myPoints->Array1(), - myTangents->ChangeArray1(), - myTangentFlags->ChangeArray1(), - myParameters->Array1()) ; + myTangents->ChangeArray1(), + myTangentFlags->ChangeArray1(), + myParameters->Array1(), + myBCType); } contact_order_array.SetValue(2,1) ; parameters.SetValue(1,myParameters->Value(1)) ; @@ -787,9 +799,10 @@ void Geom2dAPI_Interpolate::PerformNonPeriodic() // if those where not given in advance // BuildTangents(myPoints->Array1(), - myTangents->ChangeArray1(), - myTangentFlags->ChangeArray1(), - myParameters->Array1()) ; + myTangents->ChangeArray1(), + myTangentFlags->ChangeArray1(), + myParameters->Array1(), + myBCType); } contact_order_array.SetValue(2,1) ; parameters.SetValue(1,myParameters->Value(1)) ; diff --git a/src/Geom2dAPI/Geom2dAPI_Interpolate.hxx b/src/Geom2dAPI/Geom2dAPI_Interpolate.hxx index 2cf1df71e8..f0178a395a 100644 --- a/src/Geom2dAPI/Geom2dAPI_Interpolate.hxx +++ b/src/Geom2dAPI/Geom2dAPI_Interpolate.hxx @@ -17,6 +17,8 @@ #ifndef _Geom2dAPI_Interpolate_HeaderFile #define _Geom2dAPI_Interpolate_HeaderFile +#include + #include #include #include @@ -31,8 +33,6 @@ class Geom2d_BSplineCurve; class StdFail_NotDone; class Standard_ConstructionError; -class gp_Vec2d; - //! This class is used to interpolate a BsplineCurve //! passing through an array of points, with a C2 @@ -54,16 +54,25 @@ public: //! Tolerance is to check if the points are not too close to one an other //! It is also used to check if the tangent vector is not too small. - //! There should be at least 2 points + //! There should be at least 2 points. //! if PeriodicFlag is True then the curve will be periodic. - Standard_EXPORT Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d)& Points, const Standard_Boolean PeriodicFlag, const Standard_Real Tolerance); + //! Last parameter sets the boundary condition type. + Standard_EXPORT Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d)& Points, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType = Interpolate_CLAMPED); //! if PeriodicFlag is True then the curve will be periodic //! Warning: - //! There should be as many parameters as there are points - //! except if PeriodicFlag is True : then there should be one more - //! parameter to close the curve - Standard_EXPORT Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d)& Points, const Handle(TColStd_HArray1OfReal)& Parameters, const Standard_Boolean PeriodicFlag, const Standard_Real Tolerance); + //! There should be as many parameters as there are points. + //! except if PeriodicFlag is True : then there should be one more. + //! parameter to close the curve. + //! Last parameter sets the boundary condition type. + Standard_EXPORT Geom2dAPI_Interpolate(const Handle(TColgp_HArray1OfPnt2d)& Points, + const Handle(TColStd_HArray1OfReal)& Parameters, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType = Interpolate_CLAMPED); //! Assigns this constrained BSpline curve to be //! tangential to vectors InitialTangent and FinalTangent @@ -100,14 +109,17 @@ Standard_EXPORT operator Handle(Geom2d_BSplineCurve)() const; //! Note: in this case, the result is given by the function Curve. Standard_EXPORT Standard_Boolean IsDone() const; + //! Get boundary condition type. + Interpolate_BCType GetBoundaryCondition() + { + return myBCType; + } - - -protected: - - - - + //! Set boundary condition type. + void SetBoundaryCondition(const Interpolate_BCType theBCType) + { + myBCType = theBCType; + } private: @@ -128,6 +140,7 @@ private: Handle(TColStd_HArray1OfReal) myParameters; Standard_Boolean myPeriodic; Standard_Boolean myTangentRequest; + Interpolate_BCType myBCType; }; diff --git a/src/GeomAPI/GeomAPI_Interpolate.cxx b/src/GeomAPI/GeomAPI_Interpolate.cxx index 0758fe7069..b691025f9d 100644 --- a/src/GeomAPI/GeomAPI_Interpolate.cxx +++ b/src/GeomAPI/GeomAPI_Interpolate.cxx @@ -94,147 +94,182 @@ static Standard_Boolean CheckParameters(const //======================================================================= //function : BuildParameters //purpose : -//======================================================================= +//======================================================================= static void BuildParameters(const Standard_Boolean PeriodicFlag, - const TColgp_Array1OfPnt& PointsArray, - Handle(TColStd_HArray1OfReal)& ParametersPtr) + const TColgp_Array1OfPnt& PointsArray, + Handle(TColStd_HArray1OfReal)& ParametersPtr) { - Standard_Integer ii, - index ; - Standard_Real distance ; - Standard_Integer - num_parameters = PointsArray.Length() ; - if (PeriodicFlag) { - num_parameters += 1 ; - } - ParametersPtr = - new TColStd_HArray1OfReal(1, - num_parameters) ; - ParametersPtr->SetValue(1,0.0e0) ; - index = 2 ; - for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) { - distance = - PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)) ; - ParametersPtr->SetValue(index, - ParametersPtr->Value(ii) + distance) ; - index += 1 ; + Standard_Integer ii, index; + Standard_Real distance; + + Standard_Integer num_parameters = PointsArray.Length(); + if (PeriodicFlag) + num_parameters += 1; + + ParametersPtr = new TColStd_HArray1OfReal(1, num_parameters); + ParametersPtr->SetValue(1, 0.0); + + index = 2; + for (ii = PointsArray.Lower() ; ii < PointsArray.Upper() ; ii++) + { + distance = PointsArray.Value(ii).Distance(PointsArray.Value(ii+1)); + ParametersPtr->SetValue(index, ParametersPtr->Value(ii) + distance); + index += 1; } - if (PeriodicFlag) { - distance = - PointsArray.Value(PointsArray.Upper()). - Distance(PointsArray.Value(PointsArray.Lower())) ; - ParametersPtr->SetValue(index, - ParametersPtr->Value(ii) + distance) ; + if (PeriodicFlag) + { + distance = PointsArray.Value(PointsArray.Upper()).Distance(PointsArray.Value(PointsArray.Lower())); + ParametersPtr->SetValue(index, ParametersPtr->Value(ii) + distance); } } + //======================================================================= //function : BuildPeriodicTangents -//purpose : +//purpose : Compute initial condition in periodic case //======================================================================= - -static void BuildPeriodicTangent( - const TColgp_Array1OfPnt& PointsArray, - TColgp_Array1OfVec& TangentsArray, - TColStd_Array1OfBoolean& TangentFlags, - const TColStd_Array1OfReal& ParametersArray) + +static void BuildPeriodicTangent(const TColgp_Array1OfPnt& PointsArray, + TColgp_Array1OfVec& TangentsArray, + TColStd_Array1OfBoolean& TangentFlags, + const TColStd_Array1OfReal& ParametersArray, + Interpolate_BCType theBCType) { - Standard_Integer - ii, - degree ; - Standard_Real *point_array, - *parameter_array, - eval_result[2][3] ; - - gp_Vec a_vector ; - - if (PointsArray.Length() < 3) { - Standard_ConstructionError::Raise(); - } - - if (!TangentFlags.Value(1)) { - degree = 3 ; - if (PointsArray.Length() == 3) { - degree = 2 ; - } - point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ; - parameter_array = - (Standard_Real *) &ParametersArray.Value(1) ; - TangentFlags.SetValue(1,Standard_True) ; - PLib::EvalLagrange(ParametersArray.Value(1), - 1, - degree, - 3, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 3 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(1,a_vector) ; + gp_Vec aCurTangent; + Standard_Integer degree = 3; + if (PointsArray.Length() < 3) + Standard_ConstructionError::Raise(); + if (PointsArray.Length() == 3 && + theBCType == Interpolate_CLAMPED) + { + degree = 2; + } + + if (TangentFlags.Value(1)) + return; // yet computed. + else + TangentFlags.SetValue(1,Standard_True); + + if (theBCType == Interpolate_NATURAL) + { + // Compute tangent in second point. + aCurTangent = PLib:: + ComputeLagrangeTangent(PointsArray, ParametersArray, PointsArray.Lower() + 1); + + // Compute tangent approximation in the first point to satisfy + // f''(0) = 0 condition approximation: + // This is approximation due to usage of the tangent approximation in the second point + // instead of the equation below directly in solver. + // + // Formula from the Golovanov book "Geometrical modeling", 2011, pp. 18: + // + // q0 = 1.5 * (p1 - p0) / (t1 - t0) - 0.5 q1 + // + // q1=p0*(t1-t2)/((t0-t1)(t0-t2)) + p1 (2t1-t0-t2)/((t1-t0)(t1-t2))+ p2(t1-t0)/((t2-t0)(t2-t1)) + // + // Where: + // t - parameter + // p - point + // q - tangent vector + gp_Vec aVec(PointsArray(1), PointsArray(2)); + aCurTangent = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent; } - } + else if (theBCType == Interpolate_CLAMPED) + { + Standard_Real *point_array, *parameter_array, eval_result[2][3]; + + if (PointsArray.Length() < 3) + Standard_ConstructionError::Raise(); + + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()); + parameter_array = (Standard_Real *) &ParametersArray.Value(1); + + PLib::EvalLagrange(ParametersArray.Value(1), 1, degree, 3, point_array[0], + parameter_array[0], eval_result[0][0]); + for (Standard_Integer ii = 1 ; ii <= 3 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + } + + TangentsArray.SetValue(1,aCurTangent); +} //======================================================================= //function : BuildTangents //purpose : //======================================================================= - + static void BuildTangents(const TColgp_Array1OfPnt& PointsArray, - TColgp_Array1OfVec& TangentsArray, - TColStd_Array1OfBoolean& TangentFlags, - const TColStd_Array1OfReal& ParametersArray) + TColgp_Array1OfVec& TangentsArray, + TColStd_Array1OfBoolean& TangentFlags, + const TColStd_Array1OfReal& ParametersArray, + const Interpolate_BCType theBCType) { - Standard_Integer ii, - degree ; - Standard_Real *point_array, - *parameter_array, - - eval_result[2][3] ; - gp_Vec a_vector ; - - degree = 3 ; + gp_Vec aCurTangent; + Standard_Integer ii, degree = 3; + if (PointsArray.Length() < 3) + Standard_ConstructionError::Raise(); + if (PointsArray.Length() == 3 && + theBCType == Interpolate_CLAMPED) + { + degree = 2; + } - if ( PointsArray.Length() < 3) { - Standard_ConstructionError::Raise(); - } - if (PointsArray.Length() == 3) { - degree = 2 ; - } - if (!TangentFlags.Value(1)) { - point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()) ; - parameter_array = - (Standard_Real *) &ParametersArray.Value(1) ; - TangentFlags.SetValue(1,Standard_True) ; - PLib::EvalLagrange(ParametersArray.Value(1), - 1, - degree, - 3, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 3 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(1,a_vector) ; - } - if (! TangentFlags.Value(TangentFlags.Upper())) { - point_array = - (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree) ; - TangentFlags.SetValue(TangentFlags.Upper(),Standard_True) ; - parameter_array = - (Standard_Real *) &ParametersArray.Value(ParametersArray.Upper() - degree) ; - PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()), - 1, - degree, - 3, - point_array[0], - parameter_array[0], - eval_result[0][0]) ; - for (ii = 1 ; ii <= 3 ; ii++) { - a_vector.SetCoord(ii,eval_result[1][ii-1]) ; - } - TangentsArray.SetValue(TangentsArray.Upper(),a_vector) ; - } -} + if (TangentFlags.Value(1) && + TangentFlags.Value(TangentFlags.Upper())) + { + return; + } + + if (theBCType == Interpolate_NATURAL) + { + // The formulas for tangent computation is stored in BuildPeriodicTangent. + if (!TangentFlags.Value(1)) + { + // Compute tangent in second point. + aCurTangent = PLib:: + ComputeLagrangeTangent(PointsArray, ParametersArray, PointsArray.Lower() + 1); + + gp_Vec aVec(PointsArray(1), PointsArray(2)); + aVec = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent; + TangentFlags.SetValue(1, Standard_True); + TangentsArray.SetValue(1, aVec); + } // if (!TangentFlags.Value(1)) + if (!TangentFlags.Value(TangentFlags.Upper())) + { + // Compute tangent in last but one point. + aCurTangent = PLib:: + ComputeLagrangeTangent(PointsArray, ParametersArray, PointsArray.Upper() - 1); + + Standard_Integer aLastIdx = PointsArray.Upper(); + gp_Vec aVec(PointsArray(aLastIdx - 1), PointsArray(aLastIdx)); + aVec = 1.5 * aVec / (ParametersArray(aLastIdx) - ParametersArray(aLastIdx - 1)) - 0.5 * aCurTangent; + TangentFlags.SetValue(TangentFlags.Upper(), Standard_True); + TangentsArray.SetValue(TangentsArray.Upper(), aVec); + } + } + else if (theBCType == Interpolate_CLAMPED) + { + Standard_Real *point_array, *parameter_array, eval_result[2][3]; + if (!TangentFlags.Value(1)) + { + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Lower()); + parameter_array = (Standard_Real *) &ParametersArray.Value(1); + PLib::EvalLagrange(ParametersArray.Value(1), 1, degree, 3, point_array[0], parameter_array[0], eval_result[0][0]) ; + for (ii = 1 ; ii <= 3 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + TangentFlags.SetValue(1, Standard_True); + TangentsArray.SetValue(1, aCurTangent); + } // if (!TangentFlags.Value(1)) + if (!TangentFlags.Value(TangentFlags.Upper())) + { + point_array = (Standard_Real *) &PointsArray.Value(PointsArray.Upper() - degree); + parameter_array = (Standard_Real *) &ParametersArray.Value(ParametersArray.Upper() - degree); + PLib::EvalLagrange(ParametersArray.Value(ParametersArray.Upper()), 1, degree, 3,point_array[0], parameter_array[0], eval_result[0][0]); + for (ii = 1 ; ii <= 3 ; ii++) + aCurTangent.SetCoord(ii,eval_result[1][ii-1]); + TangentFlags.SetValue(TangentFlags.Upper(), Standard_True); + TangentsArray.SetValue(TangentsArray.Upper(), aCurTangent); + } + } +} //======================================================================= //function : BuildTangents //purpose : scale the given tangent so that they have the length of @@ -310,40 +345,28 @@ static void ScaleTangents(const TColgp_Array1OfPnt& PointsArray, //purpose : //======================================================================= -GeomAPI_Interpolate::GeomAPI_Interpolate - (const Handle(TColgp_HArray1OfPnt)& PointsPtr, - const Standard_Boolean PeriodicFlag, - const Standard_Real Tolerance) : -myTolerance(Tolerance), -myPoints(PointsPtr), -myIsDone(Standard_False), -myPeriodic(PeriodicFlag), -myTangentRequest(Standard_False) +GeomAPI_Interpolate::GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& PointsPtr, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType) +: myTolerance(Tolerance), + myPoints(PointsPtr), + myIsDone(Standard_False), + myPeriodic(PeriodicFlag), + myTangentRequest(Standard_False), + myBCType(theBCType) { - Standard_Integer ii ; - Standard_Boolean result = - CheckPoints(PointsPtr->Array1(), - Tolerance) ; - myTangents = - new TColgp_HArray1OfVec(myPoints->Lower(), - myPoints->Upper()) ; - myTangentFlags = - new TColStd_HArray1OfBoolean(myPoints->Lower(), - myPoints->Upper()) ; - - if (!result) { - Standard_ConstructionError::Raise(); - } - BuildParameters(PeriodicFlag, - PointsPtr->Array1(), - myParameters) ; + Standard_Boolean result = CheckPoints(PointsPtr->Array1(), Tolerance); + myTangents = new TColgp_HArray1OfVec(myPoints->Lower(), myPoints->Upper()); + myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(), myPoints->Upper()); - for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) { - myTangentFlags->SetValue(ii,Standard_False) ; - } + if (!result) + Standard_ConstructionError::Raise(); - - + BuildParameters(PeriodicFlag, PointsPtr->Array1(), myParameters); + + for ( Standard_Integer ii = myPoints->Lower(); ii <= myPoints->Upper(); ii++) + myTangentFlags->SetValue(ii,Standard_False); } //======================================================================= @@ -351,51 +374,39 @@ myTangentRequest(Standard_False) //purpose : //======================================================================= -GeomAPI_Interpolate::GeomAPI_Interpolate - (const Handle(TColgp_HArray1OfPnt)& PointsPtr, - const Handle(TColStd_HArray1OfReal)& ParametersPtr, - const Standard_Boolean PeriodicFlag, - const Standard_Real Tolerance) : -myTolerance(Tolerance), -myPoints(PointsPtr), -myIsDone(Standard_False), -myParameters(ParametersPtr), -myPeriodic(PeriodicFlag), -myTangentRequest(Standard_False) +GeomAPI_Interpolate::GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& PointsPtr, + const Handle(TColStd_HArray1OfReal)& ParametersPtr, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType) +: myTolerance(Tolerance), + myPoints(PointsPtr), + myIsDone(Standard_False), + myParameters(ParametersPtr), + myPeriodic(PeriodicFlag), + myTangentRequest(Standard_False), + myBCType(theBCType) { - Standard_Integer ii ; - - - Standard_Boolean result = - CheckPoints(PointsPtr->Array1(), - Tolerance) ; + Standard_Boolean result = CheckPoints(PointsPtr->Array1(), Tolerance); + if (PeriodicFlag) + { + if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) + { + Standard_ConstructionError::Raise(); + } + } + myTangents = new TColgp_HArray1OfVec(myPoints->Lower(), myPoints->Upper()); + myTangentFlags = new TColStd_HArray1OfBoolean(myPoints->Lower(), myPoints->Upper()); - if (PeriodicFlag) { - if ((PointsPtr->Length()) + 1 != ParametersPtr->Length()) { - Standard_ConstructionError::Raise(); - } - } - myTangents = - new TColgp_HArray1OfVec(myPoints->Lower(), - myPoints->Upper()) ; - myTangentFlags = - new TColStd_HArray1OfBoolean(myPoints->Lower(), - myPoints->Upper()) ; - - if (!result) { - Standard_ConstructionError::Raise(); - } - - result = - CheckParameters(ParametersPtr->Array1()) ; - if (!result) { - Standard_ConstructionError::Raise(); - } - - for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) { - myTangentFlags->SetValue(ii,Standard_False) ; - } - + if (!result) + Standard_ConstructionError::Raise(); + + result = CheckParameters(ParametersPtr->Array1()); + if (!result) + Standard_ConstructionError::Raise(); + + for (Standard_Integer ii = myPoints->Lower(); ii <= myPoints->Upper(); ii++) + myTangentFlags->SetValue(ii,Standard_False); } //======================================================================= //function : Load @@ -573,14 +584,14 @@ void GeomAPI_Interpolate::PerformPeriodic() mults.SetValue(num_distinct_knots ,half_order) ; if (num_points >= 3) { -// -// only enter here if there are more than 3 points otherwise -// it means we have already the tangent -// + + // only enter here if there are more than 3 points otherwise + // it means we have already the tangent BuildPeriodicTangent(myPoints->Array1(), - myTangents->ChangeArray1(), - myTangentFlags->ChangeArray1(), - myParameters->Array1()) ; + myTangents->ChangeArray1(), + myTangentFlags->ChangeArray1(), + myParameters->Array1(), + myBCType); } contact_order_array.SetValue(2,1) ; parameters.SetValue(1,myParameters->Value(1)) ; @@ -789,14 +800,15 @@ void GeomAPI_Interpolate::PerformNonPeriodic() // check if the boundary conditions are set // if (num_points >= 3) { -// -// cannot build the tangents with degree 3 with only 2 points -// if those where not given in advance -// + // + // cannot build the tangents with degree 3 with only 2 points + // if those where not given in advance + // BuildTangents(myPoints->Array1(), - myTangents->ChangeArray1(), - myTangentFlags->ChangeArray1(), - myParameters->Array1()) ; + myTangents->ChangeArray1(), + myTangentFlags->ChangeArray1(), + myParameters->Array1(), + myBCType); } contact_order_array.SetValue(2,1) ; parameters.SetValue(1,myParameters->Value(1)) ; diff --git a/src/GeomAPI/GeomAPI_Interpolate.hxx b/src/GeomAPI/GeomAPI_Interpolate.hxx index de8b77b277..dd25b70406 100644 --- a/src/GeomAPI/GeomAPI_Interpolate.hxx +++ b/src/GeomAPI/GeomAPI_Interpolate.hxx @@ -17,6 +17,8 @@ #ifndef _GeomAPI_Interpolate_HeaderFile #define _GeomAPI_Interpolate_HeaderFile +#include + #include #include #include @@ -30,8 +32,6 @@ #include #include -class gp_Vec; - //! This class is used to interpolate a BsplineCurve //! passing through an array of points, with a C2 @@ -94,7 +94,11 @@ public: //! - conditions relating to the respective //! number of elements in the parallel tables //! Points and Parameters are not respected. - Standard_EXPORT GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& Points, const Standard_Boolean PeriodicFlag, const Standard_Real Tolerance); + //! Last parameter sets the boundary condition type. + Standard_EXPORT GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& Points, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType = Interpolate_CLAMPED); //! Initializes an algorithm for constructing a //! constrained BSpline curve passing through the points of the table @@ -134,7 +138,12 @@ public: //! - conditions relating to the respective //! number of elements in the parallel tables //! Points and Parameters are not respected. - Standard_EXPORT GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& Points, const Handle(TColStd_HArray1OfReal)& Parameters, const Standard_Boolean PeriodicFlag, const Standard_Real Tolerance); + //! Last parameter sets the boundary condition type. + Standard_EXPORT GeomAPI_Interpolate(const Handle(TColgp_HArray1OfPnt)& Points, + const Handle(TColStd_HArray1OfReal)& Parameters, + const Standard_Boolean PeriodicFlag, + const Standard_Real Tolerance, + const Interpolate_BCType theBCType = Interpolate_CLAMPED); //! Assigns this constrained BSpline curve to be //! tangential to vectors InitialTangent and FinalTangent @@ -173,14 +182,17 @@ Standard_EXPORT operator Handle(Geom_BSplineCurve)() const; //! Note: in this case, the result is given by the function Curve. Standard_EXPORT Standard_Boolean IsDone() const; + //! Get boundary condition type. + Interpolate_BCType GetBoundaryCondition() + { + return myBCType; + } - - -protected: - - - - + //! Set boundary condition type. + void SetBoundaryCondition(const Interpolate_BCType theBCType) + { + myBCType = theBCType; + } private: @@ -201,6 +213,7 @@ private: Handle(TColStd_HArray1OfReal) myParameters; Standard_Boolean myPeriodic; Standard_Boolean myTangentRequest; + Interpolate_BCType myBCType; // type of boundary conditions }; diff --git a/src/GeomLib/GeomLib_Interpolate.cxx b/src/GeomLib/GeomLib_Interpolate.cxx index 3b6b677262..a269923175 100644 --- a/src/GeomLib/GeomLib_Interpolate.cxx +++ b/src/GeomLib/GeomLib_Interpolate.cxx @@ -145,4 +145,3 @@ Handle(Geom_BSplineCurve) GeomLib_Interpolate::Curve() const { return myCurve ; } - diff --git a/src/GeomLib/GeomLib_Interpolate.hxx b/src/GeomLib/GeomLib_Interpolate.hxx index dcdb2ef180..a86dc159ff 100644 --- a/src/GeomLib/GeomLib_Interpolate.hxx +++ b/src/GeomLib/GeomLib_Interpolate.hxx @@ -26,11 +26,13 @@ #include #include #include +#include +#include +#include class Geom_BSplineCurve; class StdFail_NotDone; class Standard_OutOfRange; - //! this class is used to construct a BSpline curve by //! interpolation of points at given parameters The //! continuity of the curve is degree - 1 and the @@ -50,17 +52,16 @@ public: //! returns if everything went OK - Standard_Boolean IsDone() const; - + Standard_Boolean IsDone() const; + //! returns the error type if any - GeomLib_InterpolationErrors Error() const; - + GeomLib_InterpolationErrors Error() const; + //! returns the interpolated curve of the requested degree Standard_EXPORT Handle(Geom_BSplineCurve) Curve() const; - protected: diff --git a/src/GeometryTest/GeometryTest_ConstraintCommands.cxx b/src/GeometryTest/GeometryTest_ConstraintCommands.cxx index 7dd9cca968..90dbd18348 100644 --- a/src/GeometryTest/GeometryTest_ConstraintCommands.cxx +++ b/src/GeometryTest/GeometryTest_ConstraintCommands.cxx @@ -448,11 +448,22 @@ static Standard_Integer lintang (Draw_Interpretor& di,Standard_Integer n, const static Standard_Integer interpol (Draw_Interpretor& di,Standard_Integer n, const char** a) //================================================================================== { - if (n == 1) { - di <<"give a name to your curve !\n"; + if (n == 1) + { + di <<"Interpolate: \n" + << "interpol resCurveName \n" + << "interpolation via getting points from the axonometric view \n" + << "\n" + << "interpol resCurveName fileWithPoints.txt \n" + << "interpolate curve using data from file \n" + << "\n" + << "interpol {3d/2d} {Natural/Clamped} {x y [z] ...} \n" + << "interpolate dataset from input \n"; return 0; } - if (n == 2) { + if (n == 2) + { + // Read points from axonometric viewer. Standard_Integer id,XX,YY,b, i, j; di << "Pick points \n"; dout.Select(id, XX, YY, b); @@ -463,7 +474,7 @@ static Standard_Integer interpol (Draw_Interpretor& di,Standard_Integer n, const gp_Pnt2d P2d; Standard_Boolean newcurve; - if (dout.Is3D(id)) { + if (dout.Is3D(id)){ Handle(Draw_Marker3D) mark; Handle(TColgp_HArray1OfPnt) Points = new TColgp_HArray1OfPnt(1, 1); P.SetCoord((Standard_Real)XX/zoom,(Standard_Real)YY/zoom, 0.0); @@ -584,51 +595,171 @@ static Standard_Integer interpol (Draw_Interpretor& di,Standard_Integer n, const } } - else if (n == 3) { - // lecture du fichier. - // nbpoints, 2d ou 3d, puis valeurs. + else if (n == 3) + { + // Read data set from file const char* nomfic = a[2]; ifstream iFile(nomfic, ios::in); - if (!iFile) return 1; - Standard_Integer nbp, i; + if (!iFile) + return 1; + Standard_Integer nbp; Standard_Real x, y, z; iFile >> nbp; char dimen[3]; - iFile >> dimen; - if (!strcmp(dimen,"3d")) { - Handle(TColgp_HArray1OfPnt) Point = - new TColgp_HArray1OfPnt(1, nbp); - for (i = 1; i <= nbp; i++) { + std::string aBC; // boundary condition type. + iFile >> dimen >> aBC; + + Interpolate_BCType aBCType = Interpolate_CLAMPED; + if (!aBC.compare("Natural")) + aBCType = Interpolate_NATURAL; + else if (!aBC.compare("Clamped")) + aBCType = Interpolate_CLAMPED; + else + { + di << "Incorrect boundary type, \n" + << "supported boundary conditions are: \"Natural\" or \"Clamped\" \n"; + return 1; + } + + if (!strcmp(dimen,"3d")) + { + // 3D case. + Handle(TColgp_HArray1OfPnt) Point = new TColgp_HArray1OfPnt(1, nbp); + for (Standard_Integer i = 1; i <= nbp; i++) + { iFile >> x >> y >> z; Point->SetValue(i, gp_Pnt(x, y, z)); } - GeomAPI_Interpolate anInterpolator(Point, - Standard_False, - 1.0e-5) ; - anInterpolator.Perform() ; - if (anInterpolator.IsDone()) { - Handle(Geom_BSplineCurve) C = - anInterpolator.Curve(); + + GeomAPI_Interpolate anInterpolator(Point, Standard_False, 1.0e-5, aBCType); + anInterpolator.Perform(); + + if (anInterpolator.IsDone()) + { + Handle(Geom_BSplineCurve) C = anInterpolator.Curve(); DrawTrSurf::Set(a[1], C); } } - else if (!strcmp(dimen,"2d")) { - Handle(TColgp_HArray1OfPnt2d) PointPtr = - new TColgp_HArray1OfPnt2d(1, nbp); - for (i = 1; i <= nbp; i++) { + else if (!strcmp(dimen,"2d")) + { + // 2D case. + Handle(TColgp_HArray1OfPnt2d) PointPtr = new TColgp_HArray1OfPnt2d(1, nbp); + for (Standard_Integer i = 1; i <= nbp; i++) + { iFile >> x >> y; PointPtr->SetValue(i, gp_Pnt2d(x, y)); } - Geom2dAPI_Interpolate a2dInterpolator(PointPtr, - Standard_False, - 1.0e-5); - a2dInterpolator.Perform() ; - if (a2dInterpolator.IsDone()) { - Handle(Geom2d_BSplineCurve) C = a2dInterpolator.Curve() ; + + Geom2dAPI_Interpolate a2dInterpolator(PointPtr, Standard_False, 1.0e-5, aBCType); + a2dInterpolator.Perform(); + + if (a2dInterpolator.IsDone()) + { + Handle(Geom2d_BSplineCurve) C = a2dInterpolator.Curve(); DrawTrSurf::Set(a[1], C); } } + else + { + di << "Incorrect dimension type, \n" + << "supported dimensions are: \"3d\" or \"2d\" \n"; + return 1; + } } + else if (n > 4) + { + // Read points from draw command. + + // a[0] - method name. + // a[1] - result curve. + // a[2] - 3d/2d. + // a[3] - Natural / Clamped. + // a[4...] - points. + + std::string aDim(a[2]); + Standard_Integer aDimInt = 0; + if (!aDim.compare("3d")) + aDimInt = 3; + else if (!aDim.compare("2d")) + aDimInt = 2; + else + { + di << "Incorrect dimension type, \n" + << "supported dimensions are: \"3d\" or \"2d\" \n"; + return 1; + } + + std::string aBC(a[3]); + Interpolate_BCType aBCType = Interpolate_CLAMPED; + if (!aBC.compare("Natural")) + aBCType = Interpolate_NATURAL; + else if (!aBC.compare("Clamped")) + aBCType = Interpolate_CLAMPED; + else + { + di << "Incorrect boundary type, \n" + << "supported boundary conditions are: \"Natural\" or \"Clamped\" \n"; + return 1; + } + + // Check for dimension / consistence. + Standard_Integer aMod = (n - 4) % aDimInt; + if (aMod) + { + di << "Incorrect number of input points \n" + << "Number of points should correspond to the dimension \n"; + return 1; + } + + // Perform interpolation. + Standard_Integer aPntNb = (n - 4) / aDimInt, + anIdx = 4; + Standard_Real aCoords[3]; + if (aDimInt == 2) + { + // Read points. + Handle(TColgp_HArray1OfPnt2d) PointPtr = new TColgp_HArray1OfPnt2d(1, aPntNb); + for (Standard_Integer aPntIdx = 1; aPntIdx <= aPntNb; ++aPntIdx) + { + for (Standard_Integer aDimIdx = 1; aDimIdx <= aDimInt; ++aDimIdx) + aCoords[aDimIdx - 1] = Atof(a[anIdx++]); + + PointPtr->SetValue(aPntIdx, gp_Pnt2d(aCoords[0], aCoords[1])); + } + + // Compute result. + Geom2dAPI_Interpolate anInterpolator(PointPtr, Standard_False, 1.0e-5, aBCType); + anInterpolator.Perform(); + + if (anInterpolator.IsDone()) + { + Handle(Geom2d_BSplineCurve) aC2D = anInterpolator.Curve(); + DrawTrSurf::Set(a[1], aC2D); + } + } // if (aDimInt == 2) + else // 3d case + { + // Read points. + Handle(TColgp_HArray1OfPnt) PointPtr = new TColgp_HArray1OfPnt(1, aPntNb); + for (Standard_Integer aPntIdx = 1; aPntIdx <= aPntNb; ++aPntIdx) + { + for (Standard_Integer aDimIdx = 1; aDimIdx <= aDimInt; ++aDimIdx) + aCoords[aDimIdx - 1] = Atof(a[anIdx++]); + + PointPtr->SetValue(aPntIdx, gp_Pnt(aCoords[0], aCoords[1], aCoords[2])); + } + + // Compute result. + GeomAPI_Interpolate anInterpolator(PointPtr, Standard_False, 1.0e-5, aBCType); + anInterpolator.Perform(); + + if (anInterpolator.IsDone()) + { + Handle(Geom_BSplineCurve) aC3D = anInterpolator.Curve(); + DrawTrSurf::Set(a[1], aC3D); + } + } // else // 3d case + } // else if (n > 4) return 0; } @@ -816,7 +947,7 @@ void GeometryTest::ConstraintCommands(Draw_Interpretor& theCommands) theCommands.Add("interpol", - "interpol cname [fic]", + "run \"interpol\" without parameters for help", __FILE__, interpol, g); theCommands.Add("tanginterpol", diff --git a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx index 5de5e67bf7..4fd6e24452 100644 --- a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx +++ b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx @@ -408,8 +408,6 @@ void GeomliteTest::API2dCommands(Draw_Interpretor& theCommands) theCommands.Add("2dapprox", "2dapprox result nbpoint [curve] [[x] y [x] y...]",__FILE__, appro,g); - theCommands.Add("2dinterpole", "2dinterpole result nbpoint [curve] [[x] y [x] y ...]",__FILE__, - appro,g); g = "GEOMETRY curves and surfaces analysis"; diff --git a/src/PLib/PLib.cxx b/src/PLib/PLib.cxx index 665c719d5b..2bc7d202e4 100644 --- a/src/PLib/PLib.cxx +++ b/src/PLib/PLib.cxx @@ -2047,3 +2047,61 @@ void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder, } while (Error > Tol && NbIter <= MaxNbIter); } + +//======================================================================= +//function : computeLagrangeTangent2d +//purpose : compute tangent of 3-point template in middle point +// using Lagrange polynomial. +//======================================================================= +gp_Vec2d PLib::ComputeLagrangeTangent2d(const TColgp_Array1OfPnt2d& thePointsArray, + const TColStd_Array1OfReal& theParametersArray, + const Standard_Integer theIdx) +{ + gp_Vec2d aCurTangent(0.0, 0.0); + + if (theIdx == thePointsArray.Lower() || + theIdx == thePointsArray.Upper()) + return aCurTangent; + + // Compute tangent in the second point. + Standard_Real t0 = theParametersArray(theIdx - 1); + Standard_Real t1 = theParametersArray(theIdx); + Standard_Real t2 = theParametersArray(theIdx + 1); + gp_Pnt2d p0 = thePointsArray(theIdx - 1); + gp_Pnt2d p1 = thePointsArray(theIdx); + gp_Pnt2d p2 = thePointsArray(theIdx + 1); + aCurTangent.SetXY(p0.XY() * (t1-t2)/((t0-t1)*(t0-t2)) + + p1.XY() * (2*t1-t0-t2)/((t1-t0)*(t1-t2)) + + p2.XY() * (t1-t0)/((t2-t0)*(t2-t1))); + + return aCurTangent; +} + +//======================================================================= +//function : computeLagrangeTangent +//purpose : compute tangent of 3-point template in middle point +// using Lagrange polynomial. +//======================================================================= +gp_Vec PLib::ComputeLagrangeTangent(const TColgp_Array1OfPnt& thePointsArray, + const TColStd_Array1OfReal& theParametersArray, + const Standard_Integer theIdx) +{ + gp_Vec aCurTangent(0.0, 0.0, 0.0); + + if (theIdx == thePointsArray.Lower() || + theIdx == thePointsArray.Upper()) + return aCurTangent; + + // Compute tangent in the second point. + Standard_Real t0 = theParametersArray(theIdx - 1); + Standard_Real t1 = theParametersArray(theIdx); + Standard_Real t2 = theParametersArray(theIdx + 1); + gp_Pnt p0 = thePointsArray(theIdx - 1); + gp_Pnt p1 = thePointsArray(theIdx); + gp_Pnt p2 = thePointsArray(theIdx + 1); + aCurTangent.SetXYZ(p0.XYZ() * (t1-t2)/((t0-t1)*(t0-t2)) + + p1.XYZ() * (2*t1-t0-t2)/((t1-t0)*(t1-t2)) + + p2.XYZ() * (t1-t0)/((t2-t0)*(t2-t1))); + + return aCurTangent; +} \ No newline at end of file diff --git a/src/PLib/PLib.hxx b/src/PLib/PLib.hxx index 283c6f08d3..d712a6a911 100644 --- a/src/PLib/PLib.hxx +++ b/src/PLib/PLib.hxx @@ -36,6 +36,12 @@ class PLib_JacobiPolynomial; class PLib_HermitJacobi; class PLib_DoubleJacobiPolynomial; +//! This enumeration represent boundary condition type for interpolation algorithm. +enum Interpolate_BCType +{ + Interpolate_CLAMPED, + Interpolate_NATURAL +}; //! PLib means Polynomial functions library. This pk //! provides basic computation functions for @@ -337,7 +343,13 @@ public: Standard_EXPORT static void EvalLength (const Standard_Integer Degree, const Standard_Integer Dimension, Standard_Real& PolynomialCoeff, const Standard_Real U1, const Standard_Real U2, const Standard_Real Tol, Standard_Real& Length, Standard_Real& Error); + Standard_EXPORT static gp_Vec2d ComputeLagrangeTangent2d(const TColgp_Array1OfPnt2d &thePointsArray, + const TColStd_Array1OfReal &theParametersArray, + const Standard_Integer theIdx); + Standard_EXPORT static gp_Vec ComputeLagrangeTangent(const TColgp_Array1OfPnt &thePointsArray, + const TColStd_Array1OfReal &theParametersArray, + const Standard_Integer theIdx); protected: diff --git a/tests/bugs/moddata_3/bug27112_1 b/tests/bugs/moddata_3/bug27112_1 new file mode 100644 index 0000000000..8e87c5d2b1 --- /dev/null +++ b/tests/bugs/moddata_3/bug27112_1 @@ -0,0 +1,32 @@ +puts "========" +puts "OCC27112" +puts "========" +puts "" +############################################## +# GeomAPI_Interpolate produces wrong result +############################################## + +interpol result 3d Natural \ +0020.2145846565971 0043.0749267405586 -0004.42510603802374 \ +0023.9807408024132 0014.6827020376834 -0012.4125287876527 \ +0024.3667948231688 0013.8677143193155 -0012.641915904261 \ +0024.6513678260243 0013.5109991196997 -0012.7423191658556 \ +0024.7858704778464 0013.4143863517737 -0012.7695125913718 \ +0024.9455105694487 0013.3596568726249 -0012.7849172391444 \ +0025.0753806027073 0013.3631194708948 -0012.7839426245161 \ +0025.1971400830829 0013.4054685177828 -0012.7720226826757 \ +0025.3773723647216 0013.5376332070131 -0012.7348225311948 \ +0025.7223738766752 0014.0220793196683 -0012.5984677382982 \ +0026.5185459585299 0016.3036441304239 -0011.9563151411137 \ +0028.102183898912 0025.7227430797677 -0009.30579278073889 \ +0029.7854153434029 0043.0749267405586 -0004.42510603802373 + +# Result length check. +checklength result -l 63.2944 + +# Visual check. +donly result +smallview +fit +display result +checkview -screenshot -2d -path ${imagedir}/${test_image}.png \ No newline at end of file diff --git a/tests/bugs/moddata_3/bug27112_2 b/tests/bugs/moddata_3/bug27112_2 new file mode 100644 index 0000000000..d04a75cc70 --- /dev/null +++ b/tests/bugs/moddata_3/bug27112_2 @@ -0,0 +1,36 @@ +puts "========" +puts "OCC27112" +puts "========" +puts "" +############################################## +# GeomAPI_Interpolate produces wrong result +############################################## + +interpol result 3d Natural \ +1.0 1.0 0.0 \ +2.0 0.0 0.0 \ +3.0 1.0 0.0 \ +3.0 8.0 0.0 \ +4.0 10.0 0.0\ +5.0 9.0 0.0 \ +4.0 7.0 0.0 \ +0.0 5.0 0.0 \ +0.0 3.0 0.0 \ +3.0 2.0 0.0 \ +6.0 3.0 0.0 \ +6.0 5.0 0.0 \ +4.0 6.0 0.0 \ +3.0 6.0 0.0 \ +2.0 5.0 0.0 \ +2.0 4.0 0.0 \ +3.0 3.0 0.0 + +# Result length check. +checklength result -l 38.925 + +# Visual check. +donly result +top +fit +display result +checkview -screenshot -2d -path ${imagedir}/${test_image}.png \ No newline at end of file diff --git a/tests/bugs/moddata_3/bug27112_3 b/tests/bugs/moddata_3/bug27112_3 new file mode 100644 index 0000000000..c65bbe8bf4 --- /dev/null +++ b/tests/bugs/moddata_3/bug27112_3 @@ -0,0 +1,25 @@ +puts "========" +puts "OCC27112" +puts "========" +puts "" +############################################## +# GeomAPI_Interpolate produces wrong result +############################################## + +interpol result 3d Natural \ +1.0 1.0 0.0 \ +8.0 -2.0 0.0 \ +15.0 -6.0 0.0 \ +16.0 -8.0 0.0 \ +16.5 -7.0 0.0 \ +52.0 -1.0 0.0 \ + +# Result length check. +checklength result -l 59.3204 + +# Visual check. +donly result +top +fit +display result +checkview -screenshot -2d -path ${imagedir}/${test_image}.png \ No newline at end of file diff --git a/tests/bugs/moddata_3/bug27112_4 b/tests/bugs/moddata_3/bug27112_4 new file mode 100644 index 0000000000..12da204588 --- /dev/null +++ b/tests/bugs/moddata_3/bug27112_4 @@ -0,0 +1,32 @@ +puts "========" +puts "OCC27112" +puts "========" +puts "" +############################################## +# GeomAPI_Interpolate produces wrong result +# 2d case. Based on the well-known +# "titanium heat" data set +############################################## + +interpol result 2d Natural \ +0.0 1.0 \ +1.0 1.1 \ +2.0 1.1 \ +3.0 1.2 \ +4.0 1.3 \ +5.0 7.2 \ +6.0 3.1 \ +7.0 2.6 \ +8.0 1.9 \ +9.0 1.7 \ +10.0 1.6 + +# Check length +checklength result -l 19.112043076503149 + +# Visual check. +donly result +v2d +2dfit +display result +checkview -screenshot -2d -path ${imagedir}/${test_image}.png \ No newline at end of file