0027112: GeomAPI_Interpolate produces wrong result

Support of the approximation of natural boundary condition (f''= 0)

src/Geom/FILES

@@ -41,6 +41,8 @@ Geom_Geometry.hxx
 Geom_HSequenceOfBSplineSurface.hxx
 Geom_Hyperbola.cxx
 Geom_Hyperbola.hxx
+Geom_Interpolate.cxx
+Geom_Interpolate.hxx
 Geom_Line.cxx
 Geom_Line.hxx
 Geom_OffsetCurve.cxx

src/Geom/Geom_Interpolate.cxx

@@ -0,0 +1,73 @@
+// Copyright (c) 1994-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#include <Geom_Interpolate.hxx>
+
+//=======================================================================
+//function : computeLagrangeTangent2d
+//purpose  : compute tangent of 3-point template in middle point
+//           using Lagrange polynomial.
+//=======================================================================
+gp_Vec2d Geom_Interpolate::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 Geom_Interpolate::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;
+}

src/Geom/Geom_Interpolate.hxx

@@ -0,0 +1,48 @@
+// Copyright (c) 1994-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _Geom_Interpolate_HeaderFile
+#define _Geom_Interpolate_HeaderFile
+
+#include <gp_Vec2d.hxx>
+#include <gp_Vec.hxx>
+#include <TColgp_Array1OfPnt2d.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <TColStd_Array1OfReal.hxx>
+
+
+//! This enumeration represent boundary condition type for interpolation algorithm.
+enum GeomInterpolate_BCType
+{
+  GeomInterpolate_CLAMPED,
+  GeomInterpolate_NATURAL
+};
+
+//!  This class contains static method which are used in interpolation algorithms:
+//! GeomAPI_Interpolate
+//! Geom2dAPI_Interpolate
+class Geom_Interpolate
+{
+public:
+
+  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);
+};
+
+#endif

src/Geom2dAPI/Geom2dAPI_Interpolate.cxx

@@ -127,114 +127,150 @@ 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 GeomInterpolate_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 (TangentFlags.Value(1))
+    return; // yet computed.
+  else
+    TangentFlags.SetValue(1,Standard_True);
+
+  if (theBCType == GeomInterpolate_NATURAL)
+  {
+    // Compute tangent in second point.
+    aCurTangent = Geom_Interpolate::
+      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 == GeomInterpolate_CLAMPED)
+  {
+    Standard_Real *point_array, *parameter_array, eval_result[2][2];
+
+    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, 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 GeomInterpolate_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,
+           aFirstTangent, // Resulting tangents.
+           aLastTangent; // Resulting tangents.
+  Standard_Integer ii, degree = 3;
+  if (PointsArray.Length() == 3)
+    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 == GeomInterpolate_NATURAL)
+  {
+    // The formulas for tangent computation is stored in BuildPeriodicTangent.
+    if (!TangentFlags.Value(1))
+    {
+      // Compute tangent in second point.
+      aCurTangent = Geom_Interpolate::
+        ComputeLagrangeTangent2d(PointsArray, ParametersArray, PointsArray.Lower() + 1);
+
+      gp_Vec2d aVec(PointsArray(1), PointsArray(2));
+      aFirstTangent = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent;
+    } // if (!TangentFlags.Value(1))
+    if (!TangentFlags.Value(TangentFlags.Upper()))
+    {
+      // Compute tangent in last but one point.
+      aCurTangent = Geom_Interpolate::
+        ComputeLagrangeTangent2d(PointsArray, ParametersArray, PointsArray.Upper() - 1);
+
+      Standard_Integer aLastIdx = PointsArray.Upper();
+      gp_Vec2d aVec(PointsArray(aLastIdx - 1), PointsArray(aLastIdx));
+      aLastTangent = 1.5 * aVec / (ParametersArray(aLastIdx) - ParametersArray(aLastIdx - 1)) - 0.5 * aCurTangent;
+    }
+  }
+  else if (theBCType == GeomInterpolate_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);
+      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++)
+        aFirstTangent.SetCoord(ii,eval_result[1][ii-1]);
+    } // if (!TangentFlags.Value(1))
+    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++)
+        aLastTangent.SetCoord(ii,eval_result[1][ii-1]);
+    }
+  }
+
+  if (!TangentFlags.Value(1))
+  {
+    TangentFlags.SetValue(1, Standard_True);
+    TangentsArray.SetValue(1, aFirstTangent);
+  }
+  if (!TangentFlags.Value(TangentFlags.Upper()))
+  {
+    TangentFlags.SetValue(TangentFlags.Upper(), Standard_True);
+    TangentsArray.SetValue(TangentsArray.Upper(), aLastTangent);
+  }
+}
 //=======================================================================
 //function : BuildTangents
 //purpose  : scale the given tangent so that they have the length of
@@ -310,41 +346,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 GeomInterpolate_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()) ;
+  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) ;
+  if (!result)
+    Standard_ConstructionError::Raise();
 
-  for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
-    myTangentFlags->SetValue(ii,Standard_False) ;
-  }
+  BuildParameters(PeriodicFlag, PointsPtr->Array1(), myParameters);
 
- 
-		 
+  for (Standard_Integer ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++)
+    myTangentFlags->SetValue(ii,Standard_False);
 }
 
 //=======================================================================
@@ -352,51 +375,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 GeomInterpolate_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 +578,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 +800,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)) ;

src/Geom2dAPI/Geom2dAPI_Interpolate.hxx

@@ -17,6 +17,8 @@
 #ifndef _Geom2dAPI_Interpolate_HeaderFile
 #define _Geom2dAPI_Interpolate_HeaderFile
 
+#include <Geom_Interpolate.hxx>
+
 #include <Standard.hxx>
 #include <Standard_DefineAlloc.hxx>
 #include <Standard_Handle.hxx>
@@ -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 GeomInterpolate_BCType theBCType = GeomInterpolate_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 GeomInterpolate_BCType theBCType = GeomInterpolate_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.
+  GeomInterpolate_BCType GetBoundaryCondition()
+  {
+    return myBCType;
+  }
 
-
-
-protected:
-
-
-
-
+  //! Set boundary condition type.
+  void GetBoundaryCondition(const GeomInterpolate_BCType theBCType)
+  {
+    myBCType = theBCType;
+  }
 
 private:
 
@@ -128,6 +140,7 @@ private:
   Handle(TColStd_HArray1OfReal) myParameters;
   Standard_Boolean myPeriodic;
   Standard_Boolean myTangentRequest;
+  GeomInterpolate_BCType myBCType;
 
 
 };

src/GeomAPI/GeomAPI_Interpolate.cxx

@@ -94,147 +94,178 @@ 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 ;
+  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;
   }
-  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 ;
-  }
-  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,
+                                 GeomInterpolate_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)
+    degree = 2;
+
+  if (TangentFlags.Value(1))
+    return; // yet computed.
+  else
+    TangentFlags.SetValue(1,Standard_True);
+
+  if (theBCType == GeomInterpolate_NATURAL)
+  {
+    // Compute tangent in second point.
+    aCurTangent = Geom_Interpolate::
+      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 == GeomInterpolate_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)
-{
- Standard_Integer ii,
- degree ;
- Standard_Real *point_array,
- *parameter_array,
- 
- eval_result[2][3] ;
- gp_Vec a_vector ;
- 
- degree = 3 ;
 
- 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) ;
- }
-} 
+static void BuildTangents(const TColgp_Array1OfPnt&      PointsArray,
+                          TColgp_Array1OfVec&            TangentsArray,
+                          TColStd_Array1OfBoolean&       TangentFlags,
+                          const TColStd_Array1OfReal&    ParametersArray,
+                          const GeomInterpolate_BCType       theBCType)
+{
+  gp_Vec aFirstTangent, // Resulting tangents.
+         aLastTangent; // Resulting tangents.
+  Standard_Integer ii, degree = 3;
+  if (PointsArray.Length() == 3)
+    degree = 2;
+
+  gp_Vec aCurTangent;
+
+  if (TangentFlags.Value(1) &&
+      TangentFlags.Value(TangentFlags.Upper()))
+  {
+    return;
+  }
+
+  if (theBCType == GeomInterpolate_NATURAL)
+  {
+    // The formulas for tangent computation is stored in BuildPeriodicTangent.
+    if (!TangentFlags.Value(1))
+    {
+      // Compute tangent in second point.
+      aCurTangent = Geom_Interpolate::
+        ComputeLagrangeTangent(PointsArray, ParametersArray, PointsArray.Lower() + 1);
+
+      gp_Vec aVec(PointsArray(1), PointsArray(2));
+      aFirstTangent = 1.5 * aVec / (ParametersArray(2) - ParametersArray(1)) - 0.5 * aCurTangent;
+    } // if (!TangentFlags.Value(1))
+    if (!TangentFlags.Value(TangentFlags.Upper()))
+    {
+      // Compute tangent in last but one point.
+      aCurTangent = Geom_Interpolate::
+        ComputeLagrangeTangent(PointsArray, ParametersArray, PointsArray.Upper() - 1);
+
+      Standard_Integer aLastIdx = PointsArray.Upper();
+      gp_Vec aVec(PointsArray(aLastIdx - 1), PointsArray(aLastIdx));
+      aLastTangent = 1.5 * aVec / (ParametersArray(aLastIdx) - ParametersArray(aLastIdx - 1)) - 0.5 * aCurTangent;
+    }
+  }
+  else if (theBCType == GeomInterpolate_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++)
+        aFirstTangent.SetCoord(ii,eval_result[1][ii-1]);
+    } // 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++)
+        aLastTangent.SetCoord(ii,eval_result[1][ii-1]);
+    }
+  }
+
+  if (!TangentFlags.Value(1))
+  {
+    TangentFlags.SetValue(1, Standard_True);
+    TangentsArray.SetValue(1, aFirstTangent);
+  }
+  if (!TangentFlags.Value(TangentFlags.Upper()))
+  {
+    TangentFlags.SetValue(TangentFlags.Upper(), Standard_True);
+    TangentsArray.SetValue(TangentsArray.Upper(), aLastTangent);
+  }
+}
 //=======================================================================
 //function : BuildTangents
 //purpose  : scale the given tangent so that they have the length of
@@ -310,40 +341,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 GeomInterpolate_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()) ;
+  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) ;
+  if (!result)
+    Standard_ConstructionError::Raise();
 
-  for (ii = myPoints->Lower() ; ii <= myPoints->Upper() ; ii++) {
-    myTangentFlags->SetValue(ii,Standard_False) ;
-  }
+  BuildParameters(PeriodicFlag, PointsPtr->Array1(), myParameters);
 
- 
-		 
+  for ( Standard_Integer ii = myPoints->Lower(); ii <= myPoints->Upper(); ii++)
+    myTangentFlags->SetValue(ii,Standard_False);
 }
 
 //=======================================================================
@@ -351,51 +370,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 GeomInterpolate_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 +580,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 +796,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)) ;

src/GeomAPI/GeomAPI_Interpolate.hxx

@@ -17,6 +17,8 @@
 #ifndef _GeomAPI_Interpolate_HeaderFile
 #define _GeomAPI_Interpolate_HeaderFile
 
+#include <Geom_Interpolate.hxx>
+
 #include <Standard.hxx>
 #include <Standard_DefineAlloc.hxx>
 #include <Standard_Handle.hxx>
@@ -30,8 +32,6 @@
 #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
@@ -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 GeomInterpolate_BCType theBCType = GeomInterpolate_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 GeomInterpolate_BCType theBCType = GeomInterpolate_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.
+  GeomInterpolate_BCType GetBoundaryCondition()
+  {
+    return myBCType;
+  }
 
-
-
-protected:
-
-
-
-
+  //! Set boundary condition type.
+  void GetBoundaryCondition(const GeomInterpolate_BCType theBCType)
+  {
+    myBCType = theBCType;
+  }
 
 private:
 
@@ -201,6 +213,7 @@ private:
   Handle(TColStd_HArray1OfReal) myParameters;
   Standard_Boolean myPeriodic;
   Standard_Boolean myTangentRequest;
+  GeomInterpolate_BCType myBCType; // type of boundary conditions
 
 
 };

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;
+
+    GeomInterpolate_BCType aBCType = GeomInterpolate_CLAMPED;
+    if (!aBC.compare("Natural"))
+      aBCType = GeomInterpolate_NATURAL;
+    else if (!aBC.compare("Clamped"))
+      aBCType = GeomInterpolate_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]);
+    GeomInterpolate_BCType aBCType = GeomInterpolate_CLAMPED;
+    if (!aBC.compare("Natural"))
+      aBCType = GeomInterpolate_NATURAL;
+    else if (!aBC.compare("Clamped"))
+      aBCType = GeomInterpolate_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",

src/GeomliteTest/GeomliteTest_API2dCommands.cxx

@@ -401,8 +401,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";
 

tests/bugs/moddata_3/bug27112_1

@@ -0,0 +1,33 @@
+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.
+mkedge anEdge result
+checkprops anEdge -l 63.2944
+
+# Visual check.
+donly result
+smallview
+fit
+display result
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/moddata_3/bug27112_2

@@ -0,0 +1,37 @@
+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.
+mkedge anEdge result
+checkprops anEdge -l 38.925
+
+# Visual check.
+donly result
+top
+fit
+display result
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png

tests/bugs/moddata_3/bug27112_3

@@ -0,0 +1,26 @@
+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.
+mkedge anEdge result
+checkprops anEdge -l 59.3204
+
+# Visual check.
+donly result
+top
+fit
+display result
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png