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
//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);
}
//=======================================================================
//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
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)) ;
// 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)) ;
#ifndef _Geom2dAPI_Interpolate_HeaderFile
#define _Geom2dAPI_Interpolate_HeaderFile
+#include <PLib.hxx>
+
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <Standard_Handle.hxx>
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
//! 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
//! 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:
Handle(TColStd_HArray1OfReal) myParameters;
Standard_Boolean myPeriodic;
Standard_Boolean myTangentRequest;
+ Interpolate_BCType myBCType;
};
//=======================================================================
//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
//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);
}
//=======================================================================
//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
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)) ;
// 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)) ;
#ifndef _GeomAPI_Interpolate_HeaderFile
#define _GeomAPI_Interpolate_HeaderFile
+#include <PLib.hxx>
+
#include <Standard.hxx>
#include <Standard_DefineAlloc.hxx>
#include <Standard_Handle.hxx>
#include <TColgp_Array1OfVec.hxx>
#include <Geom_BSplineCurve.hxx>
-class gp_Vec;
-
//! This class is used to interpolate a BsplineCurve
//! passing through an array of points, with a C2
//! - 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
//! - 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
//! 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:
Handle(TColStd_HArray1OfReal) myParameters;
Standard_Boolean myPeriodic;
Standard_Boolean myTangentRequest;
+ Interpolate_BCType myBCType; // type of boundary conditions
};
#include <Standard_Integer.hxx>
#include <TColgp_Array1OfPnt.hxx>
#include <TColStd_Array1OfReal.hxx>
+#include <gp_Vec2d.hxx>
+#include <gp_Vec.hxx>
+#include <TColgp_Array1OfPnt2d.hxx>
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
//! 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:
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);
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);
}
}
- 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;
}
theCommands.Add("interpol",
- "interpol cname [fic]",
+ "run \"interpol\" without parameters for help",
__FILE__,
interpol, g);
theCommands.Add("tanginterpol",
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";
}
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
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
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:
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
--- /dev/null
+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
projects / occt-copy.git / commitdiff