0031668: Visualization - WebGL sample doesn't work on Emscripten 1.39
[occt.git] / src / ShapeConstruct / ShapeConstruct_ProjectCurveOnSurface.cxx
index 4cb1f72..0568000 100644 (file)
 //szv#4 S4163
 //:s5 abv 22.04.99  Adding debug printouts in catch {} blocks
 //#1  svv 11.01.00  Porting on DEC
-#include <ShapeConstruct_ProjectCurveOnSurface.ixx>
 
-#include <Standard_ErrorHandler.hxx>
-#include <Standard_Failure.hxx>
-
-#include <Precision.hxx>
-#include <TColStd_Array1OfInteger.hxx>
-#include <TColgp_Array1OfPnt.hxx>
-
-#include <GeomAPI_PointsToBSpline.hxx>
-#include <Geom2dAPI_Interpolate.hxx>
-#include <GeomAPI_Interpolate.hxx>
-#include <Geom2dAdaptor.hxx>
+#include <Approx_CurveOnSurface.hxx>
+#include <Geom2d_BSplineCurve.hxx>
+#include <Geom2d_Circle.hxx>
+#include <Geom2d_Curve.hxx>
+#include <Geom2d_Ellipse.hxx>
+#include <Geom2d_Hyperbola.hxx>
 #include <Geom2d_Line.hxx>
+#include <Geom2d_Parabola.hxx>
 #include <Geom2d_TrimmedCurve.hxx>
-#include <Geom2d_BSplineCurve.hxx>
+#include <Geom2dAdaptor.hxx>
+#include <Geom2dAPI_Interpolate.hxx>
+#include <Geom_BezierSurface.hxx>
+#include <Geom_BoundedCurve.hxx>
 #include <Geom_BSplineCurve.hxx>
-#include <Geom_TrimmedCurve.hxx>
-#include <Geom_RectangularTrimmedSurface.hxx>
-#include <Geom_SurfaceOfLinearExtrusion.hxx>
-#include <Geom_SphericalSurface.hxx>
+#include <Geom_BSplineSurface.hxx>
+#include <Geom_Curve.hxx>
 #include <Geom_OffsetSurface.hxx>
 #include <Geom_Plane.hxx>
-#include <GeomProjLib.hxx>
-#include <GeomAdaptor_HSurface.hxx>
+#include <Geom_RectangularTrimmedSurface.hxx>
+#include <Geom_SphericalSurface.hxx>
+#include <Geom_Surface.hxx>
+#include <Geom_SurfaceOfLinearExtrusion.hxx>
+#include <Geom_TrimmedCurve.hxx>
 #include <GeomAdaptor_HCurve.hxx>
-
+#include <GeomAdaptor_HSurface.hxx>
+#include <GeomAPI_Interpolate.hxx>
+#include <GeomAPI_PointsToBSpline.hxx>
+#include <GeomProjLib.hxx>
+#include <gp_Pnt2d.hxx>
+#include <ElCLib.hxx>
+#include <NCollection_Sequence.hxx>
+#include <Precision.hxx>
+#include <ProjLib_CompProjectedCurve.hxx>
+#include <ProjLib_HCompProjectedCurve.hxx>
+#include <ProjLib_ProjectedCurve.hxx>
+#include <ShapeAnalysis.hxx>
 #include <ShapeAnalysis_Curve.hxx>
 #include <ShapeAnalysis_Surface.hxx>
+#include <ShapeConstruct_ProjectCurveOnSurface.hxx>
 #include <ShapeExtend.hxx>
+#include <Standard_ErrorHandler.hxx>
+#include <Standard_Failure.hxx>
+#include <Standard_Type.hxx>
+#include <TColgp_Array1OfPnt.hxx>
+#include <TColStd_Array1OfInteger.hxx>
+#include <IntRes2d_Domain.hxx>
+#include <IntCurve_IntConicConic.hxx>
 
-#include <ProjLib_ProjectedCurve.hxx>
-#include <ProjLib_CompProjectedCurve.hxx>
-#include <ProjLib_HCompProjectedCurve.hxx>
-#include <Approx_CurveOnSurface.hxx>
+#include <algorithm>
+IMPLEMENT_STANDARD_RTTIEXT(ShapeConstruct_ProjectCurveOnSurface,Standard_Transient)
 
 #define NCONTROL 23
 
+
+static void AdjustSecondPointToFirstPoint(const gp_Pnt2d& theFirstPoint,
+                                          gp_Pnt2d& theSecondPoint,
+                                          const Handle(Geom_Surface)& theSurf)
+{
+  if (theSurf->IsUPeriodic())
+  {
+    Standard_Real UPeriod = theSurf->UPeriod();
+    Standard_Real NewU = ElCLib::InPeriod(theSecondPoint.X(),
+                                          theFirstPoint.X() - UPeriod/2,
+                                          theFirstPoint.X() + UPeriod/2);
+    theSecondPoint.SetX(NewU);
+  }
+  if (theSurf->IsVPeriodic())
+  {
+    Standard_Real VPeriod = theSurf->VPeriod();
+    Standard_Real NewV = ElCLib::InPeriod(theSecondPoint.Y(),
+                                          theFirstPoint.Y() - VPeriod/2,
+                                          theFirstPoint.Y() + VPeriod/2);
+    theSecondPoint.SetY(NewV);
+  }
+}
+
+
 //=======================================================================
 //function : ShapeConstruct_ProjectCurveOnSurface
 //purpose  : 
@@ -154,51 +194,26 @@ ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
    return myAdjustOverDegen;
  }
 
-
-//=======================================================================
-//function : NbSurfIntervals
-//purpose  : work-around of bug in standard method 
-//           GeomAdaptor_Surface->NbIntervals() (PRO16346)
-//=======================================================================
-
-static Standard_Integer NbSurfIntervals(const Handle(GeomAdaptor_HSurface)& GAS, const GeomAbs_Shape cont)
-{
-  Standard_Integer NbU = 0;
-  if (GAS->GetType() == GeomAbs_SurfaceOfExtrusion) {
-    // extract the surface
-    Handle(Geom_SurfaceOfLinearExtrusion) surf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(GAS->ChangeSurface().Surface());
-    // build a 3d adaptor curve
-    GeomAdaptor_Curve Adaptor3dCurve(surf->BasisCurve(), GAS->FirstUParameter(), GAS->LastUParameter());
-    if (Adaptor3dCurve.GetType() == GeomAbs_BSplineCurve)
-      NbU = Adaptor3dCurve.NbIntervals(cont);
-  }
-  if (NbU == 0)
-    NbU = GAS->NbUIntervals(cont);
-  return NbU * (GAS->NbVIntervals(cont));
-}
-
 //=======================================================================
 //function : Status
 //purpose  : 
 //=======================================================================
 
- Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status Status) const
+ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status theStatus) const
 {
-  return ShapeExtend::DecodeStatus (myStatus, Status);
+  return ShapeExtend::DecodeStatus (myStatus, theStatus);
 }
 
 //=======================================================================
 //function : Perform
 //purpose  : 
 //=======================================================================
-
 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
                                                                 const Standard_Real First,
                                                                 const Standard_Real Last,
                                                                 Handle(Geom2d_Curve)& c2d,
-                                                                const GeomAbs_Shape,
-                                                                const Standard_Integer,
-                                                                const Standard_Integer)
+                                                                const Standard_Real TolFirst,
+                                                                const Standard_Real TolLast)
 {
   myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
@@ -208,7 +223,6 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curv
     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
     return Standard_False;
   }
-
 //  Projection Analytique
   Handle(Geom_Curve) crv3dtrim = c3d;
   if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
@@ -240,40 +254,130 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curv
   if ( ! bspl.IsNull() ) {
     Standard_Integer nint = 0;
     for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
-      if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last ) nint++;
+    {
+      if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last )
+        nint++;
+    }
     Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
     while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
 #ifdef OCCT_DEBUG
     if ( nbPini > NCONTROL ) 
-      cout << "Warning: number of points for projecting is " << nbPini << endl;
+      std::cout << "Warning: number of points for projecting is " << nbPini << std::endl;
 #endif
   }
 
 //    $$$$    end :92 (big BSplineCurve C0)
   
   // this number should be "parametric dependent"
-  TColgp_Array1OfPnt   points(1, nbPini);
-  TColStd_Array1OfReal params(1, nbPini);
-
-  Standard_Integer iPnt;
+  TColgp_SequenceOfPnt points;
+  TColStd_SequenceOfReal params;
+  NCollection_Sequence<Standard_Real> aKnotCoeffs;
   gp_Pnt p3d;
+  Standard_Integer iPnt;
+
+  // In case of bspline compute parametrization speed on each 
+  // knot interval inside [aFirstParam, aLastParam].
+  // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
+  // use PerformByProjLib algorithm.
+  if(!bspl.IsNull())
+  {
+    Standard_Real aFirstParam = First; // First parameter of current interval.
+    Standard_Real aLastParam = Last; // Last parameter of current interval.
+
+    // First index computation.
+    Standard_Integer anIdx = 1;
+    for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
+    {
+      if(bspl->Knot(anIdx) > First)
+      {
+        break;
+      }
+    }
+
+    GeomAdaptor_Curve aC3DAdaptor(c3d);
+    Standard_Real aMinParSpeed = Precision::Infinite(); // Minimal parameterization speed.
+    for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
+    {
+      // Fill current knot interval.
+      aLastParam = Min(Last, bspl->Knot(anIdx));
+      Standard_Integer aNbIntPnts = NCONTROL;
+      // Number of inner points is adapted according to the length of the interval
+      // to avoid a lot of calculations on small range of parameters.
+      if (anIdx > 1)
+      {
+        const Standard_Real aLenThres = 1.e-2;
+        const Standard_Real aLenRatio =
+            (aLastParam - aFirstParam) / (bspl->Knot(anIdx) - bspl->Knot(anIdx - 1));
+        if (aLenRatio < aLenThres)
+        {
+          aNbIntPnts = Standard_Integer(aLenRatio / aLenThres * aNbIntPnts);
+          if (aNbIntPnts < 2)
+            aNbIntPnts = 2;
+        }
+      }
+      Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
+      Standard_Integer anIntIdx;
+      gp_Pnt p3d1, p3d2;
+      // Start filling from first point.
+      aC3DAdaptor.D0(aFirstParam, p3d1);
+
+      Standard_Real aLength3d = 0.0;
+      for(anIntIdx = 1; anIntIdx < aNbIntPnts; anIntIdx++)
+      {
+        Standard_Real aParam = aFirstParam + aStep * anIntIdx;
+        aC3DAdaptor.D0 (aParam, p3d2);
+        const Standard_Real aDist = p3d2.Distance(p3d1);
+
+        aLength3d += aDist;
+        p3d1 = p3d2;
+
+        aMinParSpeed = Min(aMinParSpeed, aDist / aStep);
+      }
+      const Standard_Real aCoeff = aLength3d / (aLastParam - aFirstParam);
+      if (Abs(aCoeff) > gp::Resolution())
+        aKnotCoeffs.Append(aCoeff);
+      aFirstParam = aLastParam;
+    }
+
+    Standard_Real anEvenlyCoeff = 0;
+    if (aKnotCoeffs.Size() > 0)
+    {
+      anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) / 
+                      *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
+    }
+
+    const Standard_Real aMaxQuotientCoeff = 1500.0;
+    if (anEvenlyCoeff > aMaxQuotientCoeff &&
+        aMinParSpeed > Precision::Confusion() )
+    {
+      PerformByProjLib(c3d, First, Last, c2d);
+      // PerformByProjLib fail detection:
+      if (!c2d.IsNull())
+      {
+        return Status (ShapeExtend_DONE);
+      }
+    }
+  }
+
   Standard_Real deltaT, t;
   deltaT = (Last - First) / (nbPini-1);
   nbrPnt = nbPini;
-  for (iPnt = 1; iPnt <= nbPini; iPnt ++) {
+  for (iPnt = 1; iPnt <= nbPini; iPnt ++)
+  {
     if      (iPnt == 1)      t = First;
     else if (iPnt == nbPini) t = Last;
     else                     t = First + (iPnt - 1) * deltaT;
 
     c3d->D0 (t, p3d);
-    points(iPnt) = p3d;
-    params(iPnt) = t;
+    points.Append(p3d);
+    params.Append(t);
   }
 
-//  CALCUL par approximation
-
-  TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
-  ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
+  //  CALCUL par approximation
+  TColgp_SequenceOfPnt2d pnt2d;
+  ApproxPCurve (nbrPnt,c3d,TolFirst,TolLast,
+                points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
+  nbPini = points.Length();
   if (!c2d.IsNull()) {
     myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
     return Standard_True;
@@ -307,7 +411,6 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curv
   c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 //  c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
 //  Faut-il aussi reprendre la C3D ?
-
   myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
   return Status (ShapeExtend_DONE);
 }
@@ -321,9 +424,9 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(G
                                                                        const Standard_Real First,
                                                                        const Standard_Real Last,
                                                                        Handle(Geom2d_Curve)& c2d,
-                                                                       const GeomAbs_Shape continuity,
-                                                                       const Standard_Integer maxdeg,
-                                                                       const Standard_Integer nbinterval)
+                                                                       const GeomAbs_Shape /*continuity*/,
+                                                                       const Standard_Integer /*maxdeg */,
+                                                                       const Standard_Integer /*nbinterval */)
 {
   //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
   c2d.Nullify();
@@ -332,91 +435,62 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(G
     return Standard_False;
   }
    
-  try {
+  try
+  {
     OCC_CATCH_SIGNALS
     Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
-    Standard_Real URes = GAS->ChangeSurface().UResolution ( myPreci );
-    Standard_Real VRes = GAS->ChangeSurface().VResolution ( myPreci );
     Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
-    ProjLib_CompProjectedCurve Projector ( GAS, GAC, URes, VRes );
-     
-    Standard_Real ubeg, ufin;
-    Standard_Integer nbSol = Projector.NbCurves();
-    if (nbSol==1) {
-      Projector.Bounds(1, ubeg, ufin);
-      if((ubeg<=First)&&(ufin>=Last)) {
-       Standard_Integer nbintervals = ( nbinterval < 1 ? 
-                                       NbSurfIntervals(GAS, GeomAbs_C3)+GAC->NbIntervals(GeomAbs_C3)+2:
-                                       nbinterval);
-       Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve();
-       HProjector->Set(Projector);
-       Handle(Adaptor2d_HCurve2d) HPCur = HProjector;
-       Approx_CurveOnSurface appr(HPCur, GAS, First, Last, myPreci,
-                                  continuity, maxdeg,
-                                  nbintervals,
-                                  Standard_False, Standard_True);
-       if ( appr.IsDone() )
-         c2d = appr.Curve2d();
-      }
-#ifdef OCCT_DEBUG
-      else
-       cout<<"Warning: ProjLib cutting pcurve "<< First << " -> " << ubeg <<" ; "<< Last << " -> " << ufin << endl;
-#endif     
+    ProjLib_ProjectedCurve Projector(GAS, GAC);
+
+    switch (Projector.GetType())
+    {
+    case GeomAbs_Line : 
+      c2d = new Geom2d_Line(Projector.Line()); 
+      break;
+    case GeomAbs_Circle : 
+      c2d = new Geom2d_Circle(Projector.Circle());
+      break;
+    case GeomAbs_Ellipse :
+      c2d = new Geom2d_Ellipse(Projector.Ellipse());
+      break;
+    case GeomAbs_Parabola : 
+      c2d = new Geom2d_Parabola(Projector.Parabola()); 
+      break;
+    case GeomAbs_Hyperbola : 
+      c2d = new Geom2d_Hyperbola(Projector.Hyperbola()); 
+      break;
+    case GeomAbs_BSplineCurve :
+      c2d = Projector.BSpline();
+      break;
+    default:
+      // Not possible, handling added to avoid gcc warning.
+      break;
     }
-#ifdef OCCT_DEBUG
-    else cout<<"Warning: ProjLib "<< nbSol << " curves in ProjLib"<<endl;
-#endif
-    if(c2d.IsNull()) {
+
+    if(c2d.IsNull())
+    {
       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
       return Standard_False;
     }
-    else {
+    else 
+    {
       myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
       return Standard_True;
     }
     
   }
-  catch(Standard_Failure) {
+  catch(Standard_Failure const& anException) {
 #ifdef OCCT_DEBUG
-    cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
-    Standard_Failure::Caught()->Print(cout); cout << endl;
+    std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
+    anException.Print(std::cout); std::cout << std::endl;
 #endif
+    (void)anException;
     myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
     c2d.Nullify();
   }
   return Standard_False;
 }
 
-//=======================================================================
-//function : PerformAdvanced
-//purpose  : 
-//=======================================================================
-
-Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
-                                                                       const Standard_Real First,
-                                                                       const Standard_Real Last,
-                                                                       Handle(Geom2d_Curve)& c2d)
-{
-  Standard_Boolean hasResult = Standard_False;
-  Standard_Integer nbintervals;
-  
-  Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
-// &&   (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
-
-  if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
-  if (isStandard) {
-    Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
-    Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
-    nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
-    isStandard = (nbintervals < 2);
-  }
-  if (isStandard) {
-    hasResult = PerformByProjLib(c3d, First, Last, c2d);
-  }
-  if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
-  return hasResult;
-}
-
 //=======================================================================
 //function : ProjectAnalytic
 //purpose  : 
@@ -463,19 +537,337 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   return result;
 }
 
+ //! Fix possible period jump and handle walking period parameter.
+ static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
+                                               Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
+                                               Standard_Real thePeriod, // Period on objective coord
+                                               Standard_Integer theSavedPoint, // Point number to choose period
+                                               Standard_Real theSavedParam) // Param from cashe to choose period
+{
+  Standard_Real aSavedParam;
+  Standard_Integer aSavedPoint;
+  Standard_Real aMinParam = 0.0, aMaxParam = thePeriod;
+  if (theSavedPoint < 0) {
+    // normalize to first period by default
+    aSavedParam = 0.5 * thePeriod;
+    aSavedPoint = 0;
+  }
+  else {
+    aSavedParam =  theSavedParam;
+    aSavedPoint = theSavedPoint;
+    while (aMinParam > aSavedParam) {
+      aMinParam -= thePeriod;
+      aMaxParam -= thePeriod;
+    }
+    while (aMaxParam < aSavedParam) {
+      aMinParam += thePeriod;
+      aMaxParam += thePeriod;
+    }
+  }
+
+  Standard_Real aFixIsoParam = aMinParam;
+  Standard_Boolean isIsoLine = Standard_False;
+  if (aMaxParam - aSavedParam < Precision::PConfusion() ||
+      aSavedParam - aMinParam < Precision::PConfusion()) {
+    aFixIsoParam = aSavedParam;
+    isIsoLine = Standard_True;
+  }
+  // normalize all coordinates to [aMinParam, aMaxParam)
+  for (Standard_Integer i = 0; i < 4; i++) {
+    Standard_Real aParam = thePnt[i].Coord(theIdx);
+    Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(aParam, aMinParam, aMaxParam);
+    aParam += aShift;
+    // Walk over period coord -> not walking on another isoline in parameter space.
+    if (isIsoLine) {
+      if (aMaxParam - aParam < Precision::PConfusion() || aParam - aMinParam < Precision::PConfusion())
+        aParam = aFixIsoParam;
+    }
+    else {
+      if (aMaxParam - aParam < Precision::PConfusion())
+        aParam = aMaxParam;
+      if (aParam - aMinParam < Precision::PConfusion())
+          aParam = aMinParam;
+    }
+
+    thePnt[i].SetCoord(theIdx, aParam);
+  }
+
+  // find possible period jump and increasing or decreasing coordinates vector
+  Standard_Boolean isJump = Standard_False;
+  Standard_Real aPrevDiff = 0.0;
+  Standard_Real aSumDiff = 1.0;
+  for (Standard_Integer i = 0; i < 3; i++) {
+    Standard_Real aDiff = thePnt[i + 1].Coord(theIdx) - thePnt[i].Coord(theIdx);
+    if (aDiff < -Precision::PConfusion()) {
+      aSumDiff *= -1.0;
+    }
+    //if first derivative changes its sign then period jump may exists in this place
+    if (aDiff * aPrevDiff < -Precision::PConfusion()) {
+      isJump = Standard_True;
+    }
+    aPrevDiff = aDiff;
+  }
+
+  if (!isJump)
+    return Standard_False;
+
+  if (aSumDiff > 0) { // decreasing sequence (parameters decrease twice(--) and one period jump(+))
+    for (Standard_Integer i = aSavedPoint; i > 0; i--)
+      if (thePnt[i].Coord(theIdx) > thePnt[i - 1].Coord(theIdx)) {
+        thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) + thePeriod);
+      }
+    for (Standard_Integer i = aSavedPoint; i < 3; i++)
+      if (thePnt[i].Coord(theIdx) < thePnt[i + 1].Coord(theIdx)) {
+        thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) - thePeriod);
+      }
+  }
+  else {// increasing sequence (parameters increase twice(++) and one period jump(-))
+    for (Standard_Integer i = aSavedPoint; i > 0; i--)
+      if (thePnt[i].Coord(theIdx) < thePnt[i - 1].Coord(theIdx)) {
+        thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) - thePeriod);
+      }
+    for (Standard_Integer i = aSavedPoint; i < 3; i++)
+      if (thePnt[i].Coord(theIdx) > thePnt[i + 1].Coord(theIdx)) {
+        thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) + thePeriod);
+      }
+  }
+
+  // Do not return false, because for nonlinear 2d curves vector of parameters
+  // may change its first derivative and shifted parameters will be broken for this case.
+  return Standard_True;
+ }
+
+//=======================================================================
+//function : getLine
+//purpose  : 
+//=======================================================================
+
+ Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
+   const TColgp_SequenceOfPnt& thepoints,
+   const TColStd_SequenceOfReal& theparams,
+   TColgp_SequenceOfPnt2d& thePnt2ds,
+   Standard_Real theTol,
+   Standard_Boolean &isRecompute,
+   Standard_Boolean &isFromCashe) const 
+ {
+   Standard_Integer nb =  thepoints.Length();
+   gp_Pnt aP[4];
+   aP[0] = thepoints(1);
+   aP[1] = thepoints(2);
+   aP[2] = thepoints(nb - 1);
+   aP[3] = thepoints(nb);
+   gp_Pnt2d aP2d[4];
+   Standard_Integer i = 0;
+
+   Standard_Real aTol2 = theTol * theTol;
+   Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
+   Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
+
+   // Workaround:
+   // Protection against bad "tolerance" shapes.
+   if (aTol2 > 1.0)
+   {
+     theTol = Precision::Confusion();
+     aTol2 = theTol * theTol;
+   }
+   if (aTol2 < Precision::SquareConfusion())
+     aTol2 = Precision::SquareConfusion();
+   Standard_Real anOldTol2 = aTol2;
+   // auxiliary variables to choose period for connection with previous 2dcurve (if exist)
+   Standard_Integer aSavedPointNum = -1;
+   gp_Pnt2d aSavedPoint;
+
+   // project first and last points
+   for( ; i < 4; i +=3)
+   {
+     Standard_Integer j;
+     for (j = 0; j < myNbCashe; j++)
+     {
+       if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
+       {
+         aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
+         aSavedPointNum = i;
+         aSavedPoint = myCashe2d[j];
+         if (i == 0)
+           isFromCashe = Standard_True;
+         break;
+       }
+     }
+     
+     if (j >= myNbCashe)
+       aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
+     
+     Standard_Real aDist = mySurf->Gap();
+     Standard_Real aCurDist = aDist * aDist;
+     if( aTol2 < aDist * aDist)
+       aTol2 = aCurDist;
+   }
+   
+   if ( isPeriodicU || isPeriodicV )
+   {
+     // Compute second and last but one c2d points.
+     for(i = 1; i < 3; i++)
+     {
+       Standard_Integer j;
+       for (j = 0; j < myNbCashe; j++)
+       {
+         if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
+         {
+           aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
+           aSavedPointNum = i;
+           aSavedPoint = myCashe2d[j];
+           break;
+         }
+       }
+       
+       if (j >= myNbCashe)
+         aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
+       
+       Standard_Real aDist = mySurf->Gap();
+       Standard_Real aCurDist = aDist * aDist;
+       if( aTol2 < aDist * aDist)
+         aTol2 = aCurDist;
+     }
+
+     if (isPeriodicU)
+     {
+       isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
+     }
+
+     if (isPeriodicV)
+     {
+       isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
+     }
+   }
+
+   if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
+     // Do not try to make line, because in this case may be very special case when 3d curve
+     // go over the pole of, e.g., sphere, and partly lies along seam
+     // (see ApproxPCurve() for more information).
+     return 0;
+   }
+
+   thePnt2ds.SetValue(1, aP2d[0]);
+   thePnt2ds.SetValue(nb, aP2d[3]);
+
+   // Restore old tolerance in 2d space to avoid big gap cases.
+   aTol2 = anOldTol2;
+   // Check that straight line in 2d with parameterisation as in 3d will fit 
+   // fit 3d curve at all points.
+   Standard_Real dPar = theparams(nb) - theparams(1);
+   if ( Abs(dPar) < Precision::PConfusion() )
+     return 0;
+   gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
+   gp_Vec2d aVec = aVec0 / dPar;
+   Handle(Geom_Surface) aSurf = mySurf->Surface();
+   Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
+   if (isNormalCheck) {
+     for(i = 1; i <= nb; i++)
+     {
+       gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
+       gp_Pnt aCurP;
+       gp_Vec aNormalVec, aDu, aDv;
+       aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
+       aNormalVec = aDu ^ aDv;
+       if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
+         isNormalCheck = Standard_False;
+         break;
+       }
+       gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
+       Standard_Real aDist = aNormalLine.Distance(thepoints(i));
+       if (aDist > theTol)
+         return 0;
+     }
+   }
+   if (!isNormalCheck) {
+     Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()). 
+                                     SquareDistance(thepoints(1));
+     aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
+     for(i = 2; i <  nb; i++)
+     {
+       gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
+       gp_Pnt aCurP;
+       aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
+       Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
+     
+       if(Abs (aFirstPointDist - aDist1) > aTol2)
+         return 0;
+     }
+   }
+
+   // check if pcurve can be represented by Geom2d_Line (parameterised by length)
+   Standard_Real aLLength = aVec0.Magnitude();
+   if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
+   {
+     gp_XY aDirL = aVec0.XY() / aLLength;
+     gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
+     return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
+   }
+
+   // create straight bspline
+   TColgp_Array1OfPnt2d aPoles(1, 2);
+   aPoles(1) = aP2d[0];
+   aPoles(2) = aP2d[3];
+
+   TColStd_Array1OfReal aKnots(1,2);
+   aKnots(1) = theparams(1);
+   aKnots(2) = theparams(theparams.Length());
+
+   TColStd_Array1OfInteger aMults(1,2);
+   aMults(1) = 2;
+   aMults(2) = 2;
+   Standard_Integer aDegree = 1;
+   Handle(Geom2d_BSplineCurve) abspl2d =
+     new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
+   return abspl2d;
+ }
+
 //=======================================================================
 //function : ApproxPCurve
 //purpose  : 
 //=======================================================================
 
- Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
-                                                                    const TColgp_Array1OfPnt& points,
-                                                                    const TColStd_Array1OfReal& params,
-                                                                    TColgp_Array1OfPnt2d& pnt2d,
-                                                                    Handle(Geom2d_Curve)& c2d) 
+  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
+                                                                      const Handle(Geom_Curve)& c3d,
+                                                                      const Standard_Real TolFirst,
+                                                                      const Standard_Real TolLast,
+                                                                      TColgp_SequenceOfPnt& points,
+                                                                      TColStd_SequenceOfReal& params,
+                                                                      TColgp_SequenceOfPnt2d& pnt2d,
+                                                                      Handle(Geom2d_Curve)& c2d) 
 {
-  Standard_Boolean isDone = Standard_True;
-  
+  // for performance, first try to handle typical case when pcurve is straight
+  Standard_Boolean isRecompute = Standard_False;
+  Standard_Boolean isFromCasheLine = Standard_False;
+  for (Standard_Integer iseq = 1; iseq <= nbrPnt; iseq++)
+  {
+    gp_Pnt2d aP2d(0.,0.);
+    pnt2d.Append(aP2d);
+  }
+  c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
+  if(!c2d.IsNull())
+  {
+    // fill cashe
+    Standard_Boolean ChangeCycle = Standard_False;
+    if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
+       myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
+      ChangeCycle = Standard_True;
+    myNbCashe = 2;
+    if(ChangeCycle) {
+      myCashe3d[0] = points(1);
+      myCashe3d[1] = points(nbrPnt);
+      myCashe2d[0] = pnt2d(1);
+      myCashe2d[1] = pnt2d(nbrPnt);
+    }
+    else {
+      myCashe3d[1] = points(1);
+      myCashe3d[0] = points(nbrPnt);
+      myCashe2d[1] = pnt2d(1);
+      myCashe2d[0] = pnt2d(nbrPnt);
+    }
+    return Standard_True;
+  }
+    Standard_Boolean isDone = Standard_True;
   // test if the curve 3d is a boundary of the surface 
   // (only for Bezier or BSpline surface)
   
@@ -500,7 +892,6 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   
   gp_Pnt p3d;
   gp_Pnt2d p2d;
-  Standard_Integer i;
   Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
   Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
   
@@ -584,35 +975,35 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   
   Standard_Real gap = myPreci; //:q1
   Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
+  // auxiliaruy variables to shift 2dcurve, according to previous
+  Standard_Boolean isFromCashe = Standard_False;
+  gp_Pnt2d aSavedPoint;
   if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
     //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
     if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
       ChangeCycle = Standard_True;
   //for( i = 1; i <= nbrPnt; i ++) {
   for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
-    if(ChangeCycle) //skl for OCC3430
-      i=nbrPnt-ii+1;
-    else
-      i=ii;
-    p3d = points(i);
+    const Standard_Integer aPntIndex = ChangeCycle ? (nbrPnt - ii + 1) : ii;
+    p3d = points (aPntIndex);
     if (isoParam) {
       
       if (isoPar2d3d) {
-       if (isoPar2 > isoPar1) tPar = params(i);
-       else                   tPar = t1 + t2 - params(i);
+       if (isoPar2 > isoPar1) tPar = params (aPntIndex);
+       else                   tPar = t1 + t2 - params(aPntIndex);
       } else if (!isAnalytic) {
        // projection to iso
-       if      (i==1)      tPar = isoPar1;
-       else if (i==nbrPnt) tPar = isoPar2;
+       if      (aPntIndex == 1)      tPar = isoPar1;
+       else if (aPntIndex == nbrPnt) tPar = isoPar2;
        else {
-         tPar = pout(i);
+         tPar = pout(aPntIndex);
          //:S4030  ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
        }
       }
       
       if (!isoPar2d3d && isAnalytic) {
-       if      (i == 1)      p2d = valueP1;
-       else if (i == nbrPnt) p2d = valueP2;
+       if      (aPntIndex == 1)      p2d = valueP1;
+       else if (aPntIndex == nbrPnt) p2d = valueP2;
         else {
          p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
                                      Precision::Confusion()+1000*gap); //:q1
@@ -625,20 +1016,45 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
     }
     
     else {
-      if     ( (i == 1)      && p1OnIso)  p2d = valueP1;
-      else if( (i == nbrPnt) && p2OnIso)  p2d = valueP2;
+      if     ( (aPntIndex == 1)      && p1OnIso)  p2d = valueP1;
+      else if( (aPntIndex == nbrPnt) && p2OnIso)  p2d = valueP2;
       else  {// general case (not an iso)  mais attention aux singularites !
-       if ( ii==1 ) {
-         //:q9 abv 23 Mar 99: use cashe as 1st approach
-         Standard_Integer j; // svv #1
-         for ( j=0; j < myNbCashe; j++ ) 
-           if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci ) {
-             p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci, 
-                                          Precision::Confusion()+gap);
-             break;
-           }
-         if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci); 
-       }
+        // first and last points are already computed by getLine()
+        if (aPntIndex == 1 || aPntIndex == nbrPnt)
+        {
+          if (!isRecompute)
+          {
+            p2d = pnt2d (aPntIndex);
+            gap = mySurf->Gap();
+            if (aPntIndex == 1) {
+              isFromCashe = isFromCasheLine;
+              aSavedPoint = p2d;
+            }
+            continue;
+          }
+          else
+          {
+            //:q9 abv 23 Mar 99: use cashe as 1st approach
+            Standard_Integer j; // svv #1
+            for (j = 0; j < myNbCashe; ++j)
+            {
+              if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
+              {
+                p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci, Precision::Confusion()+gap);
+                if (aPntIndex == 1)
+                {
+                  isFromCashe = Standard_True;
+                  aSavedPoint = myCashe2d[j];
+                }
+                break;
+              }
+            }
+            if (j >= myNbCashe)
+            {
+              p2d = mySurf->ValueOfUV(p3d, myPreci);
+            }
+          }
+        }
         else {
           p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
                                        Precision::Confusion()+1000*gap); //:q1
@@ -646,12 +1062,12 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
        gap = mySurf->Gap();
       }
     }
-    pnt2d (i) = p2d;
+    pnt2d (aPntIndex) = p2d;
     if ( ii > 1 ) {
       if(ChangeCycle)
-        p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
+        p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex + 1).XY() );
       else
-        p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
+        p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex - 1).XY() );
     }
   }
 
@@ -660,6 +1076,37 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
     mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
   }
+
+  //Check the extremities of 3d curve for coinciding with singularities of surf
+  //Standard_Integer NbSing = mySurf->NbSingularities(Precision::Confusion());
+  gp_Pnt PointFirst = points.First(), PointLast = points.Last();
+  Standard_Real aTolFirst = (TolFirst == -1)? Precision::Confusion() : TolFirst;
+  Standard_Real aTolLast  = (TolLast == -1)?  Precision::Confusion() : TolLast;
+  for (Standard_Integer i = 1; ; i++)
+  {
+    Standard_Real aPreci, aFirstPar, aLastPar;
+    gp_Pnt aP3d;
+    gp_Pnt2d aFirstP2d, aLastP2d;
+    Standard_Boolean IsUiso;
+    if (!mySurf->Singularity(i, aPreci, aP3d, aFirstP2d, aLastP2d, aFirstPar, aLastPar, IsUiso))
+      break;
+    if (aPreci <= Precision::Confusion() &&
+        PointFirst.Distance(aP3d) <= aTolFirst)
+    {
+      CorrectExtremity(c3d, params, pnt2d,
+                       Standard_True, //first point
+                       aFirstP2d,
+                       IsUiso);
+    }
+    if (aPreci <= Precision::Confusion() &&
+        PointLast.Distance(aP3d) <= aTolLast)
+    {
+      CorrectExtremity(c3d, params, pnt2d,
+                       Standard_False, //last point
+                       aFirstP2d,
+                       IsUiso);
+    }
+  }
   
   //  attention aux singularites ... (hors cas iso qui les traite deja)
   //  if (!isoParam) {
@@ -676,34 +1123,56 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   Standard_Real Up = ul - uf;
   Standard_Real Vp = vl - vf;
   Standard_Real dist2d;
+  const Standard_Real TolOnUPeriod = Precision::Confusion() * Up;
+  const Standard_Real TolOnVPeriod = Precision::Confusion() * Vp;
 #ifdef OCCT_DEBUG
   if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
-    cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
+    std::cout << "WARNING : Recadrage incertain sur U & VClosed" << std::endl;
   }
 #endif
   // Si la surface est UCLosed, on recadre les points
   if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
     // Premier point dans le domain [uf, ul]
     Standard_Real prevX, firstX = pnt2d (1).X();
-    while (firstX < uf)  {  firstX += Up;   pnt2d (1).SetX(firstX);  }
-    while (firstX > ul)  {  firstX -= Up;   pnt2d (1).SetX(firstX);  }
+    if (!isFromCashe) {
+      // do not shift 2dcurve, if it connects to previous
+      while (firstX < uf)  {  firstX += Up;   pnt2d (1).SetX(firstX);  }
+      while (firstX > ul)  {  firstX -= Up;   pnt2d (1).SetX(firstX);  }
+    }
+    // shift first point, according to cashe
+    if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
+      Standard_Real aMinParam = uf, aMaxParam = ul;
+      while (aMinParam > aSavedPoint.X()) {
+        aMinParam -= Up;
+        aMaxParam -= Up;
+      }
+      while (aMaxParam < aSavedPoint.X()) {
+        aMinParam += Up;
+        aMaxParam += Up;
+      }
+      Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
+      firstX += aShift;
+      pnt2d(1).SetX(firstX);
+    }
     prevX = firstX;
     
     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
     Standard_Real minX = firstX, maxX = firstX;
+    Standard_Boolean ToAdjust = Standard_False;
     
     // On decalle toujours le suivant
-    for (i = 2; i <= nbrPnt; i++) {
-      //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
-      Standard_Real CurX = pnt2d (i).X();
+    for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
+    {
+      //      dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
+      Standard_Real CurX = pnt2d (aPntIter).X();
       dist2d = Abs (CurX - prevX);
-      if (dist2d > ( Up / 2) ) {
-       if        (CurX > prevX + Up/2) {
-         while (CurX > prevX + Up/2) {  CurX -= Up;  pnt2d (i).SetX (CurX);  }
-       } else if (CurX < prevX - Up/2) {
-         while (CurX < prevX - Up/2) {  CurX += Up;  pnt2d (i).SetX (CurX);  }
-       }
-       
+      if (dist2d > ( Up / 2) )
+      {
+        InsertAdditionalPointOrAdjust(ToAdjust, 1, Up, TolOnUPeriod,
+                                      CurX, prevX,
+                                      c3d,
+                                      aPntIter,
+                                      points, params, pnt2d);
       }
       prevX = CurX;
       if ( minX > CurX ) minX = CurX;      //:97
@@ -711,12 +1180,16 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
     }
     
     //:97
-    Standard_Real midX = 0.5 * ( minX + maxX );
-    Standard_Real shiftX=0.;
-    if ( midX > ul ) shiftX = -Up;
-    else if ( midX < uf ) shiftX = Up;
-    if ( shiftX != 0. ) 
-      for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
+    if (!isFromCashe) {
+      // do not shift 2dcurve, if it connects to previous
+      Standard_Real midX = 0.5 * ( minX + maxX );
+      Standard_Real shiftX=0.;
+      if ( midX > ul ) shiftX = -Up;
+      else if ( midX < uf ) shiftX = Up;
+      if ( shiftX != 0. ) 
+        for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
+          pnt2d (aPntIter).SetX ( pnt2d (aPntIter).X() + shiftX );
+      }
   }
   // Si la surface est VCLosed, on recadre les points
   // Same code as UClosed : optimisation souhaitable !!
@@ -728,24 +1201,45 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
     // Premier point dans le domain [vf, vl]
     Standard_Real prevY, firstY = pnt2d (1).Y();
-    while (firstY < vf)  {  firstY += Vp;  pnt2d (1).SetY(firstY);  }
-    while (firstY > vl)  {  firstY -= Vp;  pnt2d (1).SetY(firstY);  }
+    if (!isFromCashe) {
+      // do not shift 2dcurve, if it connects to previous
+      while (firstY < vf)  {  firstY += Vp;  pnt2d (1).SetY(firstY);  }
+      while (firstY > vl)  {  firstY -= Vp;  pnt2d (1).SetY(firstY);  }
+    }
+    // shift first point, according to cashe
+    if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
+      Standard_Real aMinParam = vf, aMaxParam = vl;
+      while (aMinParam > aSavedPoint.Y()) {
+        aMinParam -= Vp;
+        aMaxParam -= Vp;
+      }
+      while (aMaxParam < aSavedPoint.Y()) {
+        aMinParam += Vp;
+        aMaxParam += Vp;
+      }
+      Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
+      firstY += aShift;
+      pnt2d(1).SetY(firstY);
+    }
     prevY = firstY;
     
     //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
     Standard_Real minY = firstY, maxY = firstY;
+    Standard_Boolean ToAdjust = Standard_False;
     
     // On decalle toujours le suivant
-    for (i = 2; i <= nbrPnt; i ++) {
+    for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
+    {
       //      dist2d = pnt2d (i-1).Distance(pnt2d (i));
-      Standard_Real CurY = pnt2d (i).Y();
+      Standard_Real CurY = pnt2d (aPntIter).Y();
       dist2d = Abs (CurY - prevY);
-      if (dist2d > ( Vp / 2) ) {
-       if        (CurY > prevY + Vp/2) {
-         while (CurY > prevY + Vp/2) {  CurY -= Vp;  pnt2d (i).SetY (CurY);  }
-       } else if (CurY < prevY - Vp/2) {
-         while (CurY < prevY - Vp/2) {  CurY += Vp;  pnt2d (i).SetY (CurY);  }
-       }
+      if (dist2d > ( Vp / 2) )
+      {
+        InsertAdditionalPointOrAdjust(ToAdjust, 2, Vp, TolOnVPeriod,
+                                      CurY, prevY,
+                                      c3d,
+                                      aPntIter,
+                                      points, params, pnt2d);
       }
       prevY = CurY;
       if ( minY > CurY ) minY = CurY;      //:97
@@ -753,20 +1247,24 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
     }
     
     //:97
-    Standard_Real midY = 0.5 * ( minY + maxY );
-    Standard_Real shiftY=0.;
-    if ( midY > vl ) shiftY = -Vp;
-    else if ( midY < vf ) shiftY = Vp;
-    if ( shiftY != 0. ) 
-      for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
+    if (!isFromCashe) {
+      // do not shift 2dcurve, if it connects to previous
+      Standard_Real midY = 0.5 * ( minY + maxY );
+      Standard_Real shiftY=0.;
+      if ( midY > vl ) shiftY = -Vp;
+      else if ( midY < vf ) shiftY = Vp;
+      if ( shiftY != 0. ) 
+        for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
+          pnt2d(aPntIter).SetY ( pnt2d(aPntIter).Y() + shiftY );
+      }
   }
   
   //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
   //#78 rln 12.03.99 S4135
   if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
-    for (i = 2; i <= nbrPnt; i++) {
+    for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter) {
       //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
-      dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
+      dist2d = Abs (pnt2d (aPntIter).Y() - pnt2d (aPntIter - 1).Y());
       if (dist2d > ( Vp / 2) ) {
        // ATTENTION : il faut regarder ou le decalage se fait.
        // si plusieurs points sont decalles, il faut plusieurs passes
@@ -782,10 +1280,10 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
        Standard_Boolean currOnLast  = Standard_False;
        
        //  .X ?  plutot .Y ,  non ?
-       Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
-       Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
-       Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
-       Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
+       Standard_Real distPrevVF = Abs(pnt2d (aPntIter - 1).Y() - vf);
+       Standard_Real distPrevVL = Abs(pnt2d (aPntIter - 1).Y() - vl);
+       Standard_Real distCurrVF = Abs(pnt2d (aPntIter).Y() - vf);
+       Standard_Real distCurrVL = Abs(pnt2d (aPntIter).Y() - vl);
        
        Standard_Real theMin = distPrevVF;
        prevOnFirst = Standard_True;
@@ -809,30 +1307,30 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
        }
        //  Modifs RLN/Nijni  3-DEC-1997
        if (prevOnFirst) {
-         // on decalle le point (i-1) en V Last
-         gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of  vl RLN/Nijni
-         pnt2d (i-1) = newPrev;
+         // on decalle le point (aPntIter - 1) en V Last
+         gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vf); // instead of  vl RLN/Nijni
+         pnt2d (aPntIter - 1) = newPrev;
        }
        else if (prevOnLast) {
-         // on decalle le point (i-1) en V first
-         gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of  vf RLN/Nijni
-         pnt2d (i-1) = newPrev;
+         // on decalle le point (aPntIter - 1) en V first
+         gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vl); // instead of  vf RLN/Nijni
+         pnt2d (aPntIter - 1) = newPrev;
        }
        else if (currOnFirst) {
-         // on decalle le point (i) en V Last
-         gp_Pnt2d newCurr(pnt2d (i).X(),vf);  // instead of vl  RLN/Nijni
-         pnt2d (i) = newCurr;
+         // on decalle le point (aPntIter) en V Last
+         gp_Pnt2d newCurr(pnt2d (aPntIter).X(),vf);  // instead of vl  RLN/Nijni
+         pnt2d (aPntIter) = newCurr;
        }
        else if (currOnLast) {
-         // on decalle le point (i) en V First
-         gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf  RLN/Nijni
-         pnt2d (i) = newCurr;
+         // on decalle le point (aPntIter) en V First
+         gp_Pnt2d newCurr(pnt2d (aPntIter).X(), vl); // instead of vf  RLN/Nijni
+         pnt2d (aPntIter) = newCurr;
        }
        // on verifie
 #ifdef OCCT_DEBUG
-       dist2d = pnt2d (i-1).Distance(pnt2d (i));
+       dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
        if (dist2d > ( Vp / 2) ) {
-         cout << "Echec dans le recadrage" << endl;
+         std::cout << "Echec dans le recadrage" << std::endl;
        }
 #endif
       }
@@ -864,7 +1362,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
        Standard_Integer OnBound=0, PrevOnBound=0;
        Standard_Integer ind; // svv #1
        Standard_Boolean start = Standard_True;
-       for ( ind=1; ind <= nbrPnt; ind++ ) {
+       for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
          Standard_Real CurX = pnt2d(ind).X();
          // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
          if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
@@ -878,7 +1376,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
          PrevOnBound = OnBound;
        }
         // if found, adjust seam part
-        if ( ind <= nbrPnt ) {
+        if ( ind <= pnt2d.Length() ) {
          PrevX = ( myAdjustOverDegen ? uf : ul );
          Standard_Real dU = Up/2 + Precision::PConfusion();
          if ( PrevOnBound ) {
@@ -891,7 +1389,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
          }
          else if ( OnBound ) {
            pnt2d(ind).SetX ( PrevX );
-           for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
+           for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
              Standard_Real CurX = pnt2d(j).X();
              while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
              while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
@@ -909,7 +1407,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
        Standard_Integer OnBound=0, PrevOnBound=0;
        Standard_Integer ind; // svv #1
        Standard_Boolean start = Standard_True;
-       for ( ind=1; ind <= nbrPnt; ind++ ) {
+       for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
          Standard_Real CurY = pnt2d(ind).Y();
          // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
          if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
@@ -923,7 +1421,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
          PrevOnBound = OnBound;
        }
         // if found, adjust seam part
-        if ( ind <= nbrPnt ) {
+        if ( ind <= pnt2d.Length() ) {
          PrevY = ( myAdjustOverDegen ? vf : vl );
          Standard_Real dV = Vp/2 + Precision::PConfusion();
          if ( PrevOnBound ) {
@@ -936,7 +1434,7 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
          }
          else if ( OnBound ) {
            pnt2d(ind).SetY ( PrevY );
-           for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
+           for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
              Standard_Real CurY = pnt2d(j).Y();
              while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
              while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
@@ -954,31 +1452,16 @@ Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(G
   //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
   //   myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
     myCashe3d[0] = points(1);
-    myCashe3d[1] = points(nbrPnt);
+    myCashe3d[1] = points.Last();
     myCashe2d[0] = pnt2d(1);
-    myCashe2d[1] = pnt2d(nbrPnt);
+    myCashe2d[1] = pnt2d.Last();
   }
   else {
     myCashe3d[1] = points(1);
-    myCashe3d[0] = points(nbrPnt);
+    myCashe3d[0] = points.Last();
     myCashe2d[1] = pnt2d(1);
-    myCashe2d[0] = pnt2d(nbrPnt);
+    myCashe2d[0] = pnt2d.Last();
   }
-  
-  if (isoParam && isoPar2d3d) {
-    
-    // create directly isoparametrics (PCurve)
-    
-    gp_Vec2d aVec2d(valueP1, valueP2);
-    gp_Dir2d aDir2d(aVec2d);
-    gp_Pnt2d P0;
-    
-    if (isoTypeU) P0.SetCoord(valueP1.X(), valueP1.Y() - params(1)*aDir2d.Y());
-    else          P0.SetCoord(valueP1.X() - params(1)*aDir2d.X(), valueP1.Y());
-    
-    c2d = new Geom2d_Line(P0, aDir2d);
-  }
-
   return isDone;
 }
 
@@ -1029,14 +1512,15 @@ Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(con
     C2d = new Geom2d_BSplineCurve  ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
     return C2d;
   }
-  catch(Standard_Failure) {
-#ifdef OCCT_DEBUG //:s5
+  catch(Standard_Failure const& anException) {
+#ifdef OCCT_DEBUG
+ //:s5
     //    debug ...
     Standard_Integer nbp = params->Length();
     Standard_Integer nb2 = points2d->Length();
-    cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
-    Standard_Failure::Caught()->Print(cout); 
-    cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
+    std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
+    anException.Print(std::cout); 
+    std::cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
 //     if (nb2 > nbp) nb2 = nbp;
 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
 //     //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
@@ -1046,6 +1530,7 @@ Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(con
 //       dbl.AddXYZ (quoi);
 //     }
 #endif
+     (void)anException;
      C2d.Nullify();
   }
   return C2d;
@@ -1073,14 +1558,15 @@ Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(con
     myInterPol2d.Perform();
     if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
   }
-  catch(Standard_Failure) {
-#ifdef OCCT_DEBUG //:s5
+  catch(Standard_Failure const& anException) {
+#ifdef OCCT_DEBUG
+//:s5
 // //    debug ...
     Standard_Integer nbp = params->Length();
     Standard_Integer nb2 = points2d->Length();
-    cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
-    Standard_Failure::Caught()->Print(cout); 
-    cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
+    std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
+    anException.Print(std::cout); 
+    std::cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
 //     if (nb2 > nbp) nb2 = nbp;
 //     Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
 // //    dbl.AddString ("NbP2d/NbParams puis  X Y Param -> mini");
@@ -1090,6 +1576,7 @@ Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(con
 //       dbl.AddXYZ (quoi);
 //     }
 #endif
+    (void)anException;
     C2d.Nullify();
   }
   return C2d;
@@ -1114,16 +1601,219 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
     myInterPol.Perform();
     if (myInterPol.IsDone()) C3d = myInterPol.Curve();
   }
-  catch(Standard_Failure) {
-    C3d.Nullify();
-#ifdef OCCT_DEBUG //:s5
-    cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
-    Standard_Failure::Caught()->Print(cout); cout << endl;
+  catch(Standard_Failure const& anException) {
+#ifdef OCCT_DEBUG
+ //:s5
+    std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
+    anException.Print(std::cout); std::cout << std::endl;
 #endif
+    (void)anException;
+    C3d.Nullify();
   }
   return C3d;
 }
 
+//============================================================================================
+//function : CorrectExtremity
+//purpose  : corrects first or last 2d point of future curve 
+//           in the case when it coincids with a singularity of surface
+//============================================================================================
+
+void ShapeConstruct_ProjectCurveOnSurface::CorrectExtremity(const Handle(Geom_Curve)& theC3d,
+                                                            const TColStd_SequenceOfReal& theParams,
+                                                            TColgp_SequenceOfPnt2d& thePnt2d,
+                                                            const Standard_Boolean theIsFirstPoint,
+                                                            const gp_Pnt2d& thePointOnIsoLine,
+                                                            const Standard_Boolean theIsUiso)
+{
+  Standard_Integer NbPnt = thePnt2d.Length();
+  Standard_Integer IndCoord = (theIsUiso)? 2 : 1;
+  Standard_Real SingularityCoord = thePointOnIsoLine.Coord(3-IndCoord);
+  gp_Pnt2d EndPoint = (theIsFirstPoint)? thePnt2d(1) : thePnt2d(NbPnt);
+  Standard_Real FinishCoord = EndPoint.Coord(3-IndCoord); //the constant coord of isoline
+  
+  gp_Dir2d aDir = (theIsUiso)? gp::DY2d() : gp::DX2d();
+  gp_Lin2d anIsoLine(EndPoint, aDir);
+  IntRes2d_Domain Dom1, Dom2;
+
+  Standard_Boolean IsPeriodic = (theIsUiso)?
+    mySurf->Surface()->IsVPeriodic() : mySurf->Surface()->IsUPeriodic();
+
+  gp_Pnt2d FirstPointOfLine, SecondPointOfLine;
+  Standard_Real FinishParam, FirstParam, SecondParam;
+
+  if (theIsFirstPoint)
+  {
+    FirstPointOfLine = thePnt2d(3);
+    SecondPointOfLine = thePnt2d(2);
+    FinishParam = theParams(1);
+    FirstParam = theParams(3);
+    SecondParam = theParams(2);
+  }
+  else //last point
+  {
+    FirstPointOfLine = thePnt2d(NbPnt-2);
+    SecondPointOfLine = thePnt2d(NbPnt-1);
+    FinishParam = theParams(NbPnt);
+    FirstParam = theParams(NbPnt-2);
+    SecondParam = theParams(NbPnt-1);
+  }
+
+  if (SingularityCoord > FinishCoord &&
+      SecondPointOfLine.Coord(3-IndCoord) > FinishCoord)
+    return; //the curve passes through the singularity, do nothing
+  if (SingularityCoord < FinishCoord &&
+      SecondPointOfLine.Coord(3-IndCoord) < FinishCoord)
+    return; //the curve passes through the singularity, do nothing
+  //Check correctness of <EndPoint>
+  {
+    const Standard_Real aPrevDist = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
+    const Standard_Real aCurDist  = Abs(EndPoint.Coord(IndCoord) - SecondPointOfLine.Coord(IndCoord));
+    if (aCurDist <= 2 * aPrevDist)
+      return;
+  }
+  
+  gp_Pnt2d FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
+    gp_Pnt2d(SecondPointOfLine.X(), FinishCoord); //first approximation of <FinishPoint>
+
+  for (;;)
+  {
+    if (Abs(SecondPointOfLine.Coord(3-IndCoord) - FinishCoord) <= 2*Precision::PConfusion())
+      break;
+
+    gp_Vec2d aVec(FirstPointOfLine, SecondPointOfLine);
+    Standard_Real aSqMagnitude = aVec.SquareMagnitude();
+    if (aSqMagnitude <= 1.e-32)
+      break;
+    aDir.SetCoord(aVec.X(), aVec.Y());
+    
+    gp_Lin2d aLine(FirstPointOfLine, aDir);
+    IntCurve_IntConicConic Intersector(anIsoLine, Dom1,
+                                       aLine,     Dom2,
+                                       1.e-10, 1.e-10);
+    if (Intersector.IsDone() && !Intersector.IsEmpty())
+    {
+      IntRes2d_IntersectionPoint IntPoint = Intersector.Point(1);
+      FinishPoint = IntPoint.Value();
+    }
+    else
+      FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
+        gp_Pnt2d(SecondPointOfLine.X(), FinishCoord);
+    
+    gp_Pnt2d PrevPoint = FirstPointOfLine;
+    FirstPointOfLine = SecondPointOfLine;
+    FirstParam = SecondParam;
+    SecondParam = (FirstParam + FinishParam)/2;
+    if (Abs(SecondParam - FirstParam) <= 2*Precision::PConfusion())
+      break;
+    gp_Pnt aP3d;
+    theC3d->D0(SecondParam, aP3d);
+    SecondPointOfLine = mySurf->NextValueOfUV(FirstPointOfLine, aP3d,
+                                              myPreci, Precision::Confusion());
+    if (IsPeriodic)
+      AdjustSecondPointToFirstPoint(FirstPointOfLine, SecondPointOfLine, mySurf->Surface());
+    
+    //Check <SecondPointOfLine> to be enough close to <FirstPointOfLine>
+    //because when a projected point is too close to singularity,
+    //the non-constant coordinate becomes random.
+    const Standard_Real aPrevDist = Abs(FirstPointOfLine.Coord(IndCoord) - PrevPoint.Coord(IndCoord));
+    const Standard_Real aCurDist  = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
+    if (aCurDist > 2 * aPrevDist)
+      break;
+  }
+
+  if (theIsFirstPoint)
+    thePnt2d(1) = FinishPoint;
+  else
+    thePnt2d(NbPnt) = FinishPoint;
+}
+
+//============================================================================================
+//function : InsertAdditionalPointOrAdjust
+//purpose  : If the current point is too far from the previous point 
+//           (more than half-period of surface), it can happen in two cases:
+//           1. Real current step on corresponding coordinate is small, all we need is adjust;
+//           2. Current step on corresponding coordinate is really bigger than half-period of
+//           surface in this parametric direction, so we must add additional point to exclude
+//           such big intervals between points in 2d space.
+//============================================================================================
+
+void ShapeConstruct_ProjectCurveOnSurface::
+InsertAdditionalPointOrAdjust(Standard_Boolean& ToAdjust,
+                              const Standard_Integer theIndCoord,
+                              const Standard_Real Period,
+                              const Standard_Real TolOnPeriod,
+                              Standard_Real& CurCoord,
+                              const Standard_Real prevCoord,
+                              const Handle(Geom_Curve)& c3d,
+                              Standard_Integer& theIndex,
+                              TColgp_SequenceOfPnt& points,
+                              TColStd_SequenceOfReal& params,
+                              TColgp_SequenceOfPnt2d& pnt2d)
+{
+  Standard_Real CorrectedCurCoord = ElCLib::InPeriod(CurCoord,
+                                                     prevCoord - Period/2,
+                                                     prevCoord + Period/2);
+  if (!ToAdjust)
+  {
+    Standard_Real CurPar = params(theIndex);
+    Standard_Real PrevPar = params(theIndex-1);
+    Standard_Real MidPar = (PrevPar + CurPar)/2;
+    gp_Pnt MidP3d;
+    c3d->D0(MidPar, MidP3d);
+    gp_Pnt2d MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
+    Standard_Real MidCoord = MidP2d.Coord(theIndCoord);
+    MidCoord = ElCLib::InPeriod(MidCoord, prevCoord - Period/2, prevCoord + Period/2);
+    Standard_Real FirstCoord = prevCoord, LastCoord = CorrectedCurCoord;
+    if (LastCoord < FirstCoord)
+    {Standard_Real tmp = FirstCoord; FirstCoord = LastCoord; LastCoord = tmp;}
+    if (LastCoord - FirstCoord <= TolOnPeriod)
+      ToAdjust = Standard_True;
+    else if (FirstCoord <= MidCoord && MidCoord <= LastCoord)
+      ToAdjust = Standard_True;
+    else //add mid point
+    {
+      //Standard_Real RefU = prevX;
+      Standard_Boolean Success = Standard_True;
+      Standard_Real FirstT = PrevPar; //params(i-1)
+      Standard_Real LastT  = CurPar;  //params(i)
+      MidCoord = MidP2d.Coord(theIndCoord);
+      while (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod ||
+             Abs(CurCoord - MidCoord) >= Period/2 - TolOnPeriod)
+      {
+        if (MidPar - FirstT <= Precision::PConfusion() ||
+            LastT - MidPar <= Precision::PConfusion())
+        {
+          Success = Standard_False;
+          break; //wrong choice
+        }
+        if (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod)
+          LastT = (FirstT + LastT)/2;
+        else
+          FirstT = (FirstT + LastT)/2;
+        MidPar = (FirstT + LastT)/2;
+        c3d->D0(MidPar, MidP3d);
+        MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
+        MidCoord = MidP2d.Coord(theIndCoord);
+      }
+      if (Success)
+      {
+        points.InsertBefore(theIndex, MidP3d);
+        params.InsertBefore(theIndex, MidPar);
+        pnt2d.InsertBefore(theIndex, MidP2d);
+        theIndex++;
+      }
+      else
+        ToAdjust = Standard_True; 
+    } //add mid point
+  } //if (!ToAdjust)
+  if (ToAdjust)
+  {
+    CurCoord = CorrectedCurCoord;
+    pnt2d(theIndex).SetCoord (theIndCoord, CurCoord);
+  }
+}
+
 //=======================================================================
 //function : CheckPoints
 //purpose  : 
@@ -1139,7 +1829,8 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
 
   // will store 0 when the point is to be removed, 1 otherwise
   TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
-  for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
+  for (i = firstElem; i<=lastElem ; i++)
+    tmpParam.SetValue(i,1);
   Standard_Real DistMin2 = RealLast();
   gp_Pnt Prev = points->Value (lastValid);
   gp_Pnt Curr;
@@ -1164,13 +1855,13 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
     return;
 
 #ifdef OCCT_DEBUG
-  cout << "Warning : removing 3d points for interpolation" << endl;
+  std::cout << "Warning : removing 3d points for interpolation" << std::endl;
 #endif
   // Build new HArrays
   Standard_Integer newLast = lastElem - nbPntDropped;
   if ((newLast - firstElem + 1)  < 2) {
 #ifdef OCCT_DEBUG
-    cout << "Too many degenerated points for 3D interpolation" << endl;
+    std::cout << "Too many degenerated points for 3D interpolation" << std::endl;
 #endif
     return;
   }
@@ -1235,13 +1926,13 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
     return;
 
 #ifdef OCCT_DEBUG
-  cout << "Warning : removing 2d points for interpolation" << endl;
+  std::cout << "Warning : removing 2d points for interpolation" << std::endl;
 #endif
   // Build new HArrays
   Standard_Integer newLast = lastElem - nbPntDropped;
   if ((newLast - firstElem + 1)  < 2) {
 #ifdef OCCT_DEBUG
-    cout << "Too many degenerated points for 2D interpolation" << endl;
+    std::cout << "Too many degenerated points for 2D interpolation" << std::endl;
 #endif
     //pdn 12.02.99 S4135 Creating pcurve with minimal length.
     tmpParam.SetValue(firstElem,1);
@@ -1260,7 +1951,7 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
   for (i = firstElem; i <= lastElem ; i++) {
     if (tmpParam.Value(i) == 1) { 
 #ifdef OCCT_DEBUG
-      cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
+      std::cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << std::endl;
 #endif
       newPnts->SetValue(newCurr, points->Value(i));
       newParams->SetValue(newCurr, params->Value(i));
@@ -1268,7 +1959,7 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
     }
     else {
 #ifdef OCCT_DEBUG
-      cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << endl;
+      std::cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " <<  params->Value(i) << std::endl;
 #endif
     }
   }
@@ -1284,8 +1975,8 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
 
  Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
-                                                                         const TColgp_Array1OfPnt& points,
-                                                                         const TColStd_Array1OfReal& params,
+                                                                         const TColgp_SequenceOfPnt& points,
+                                                                         const TColStd_SequenceOfReal& params,
                                                                          Standard_Boolean& isoTypeU,
                                                                          Standard_Boolean& p1OnIso,
                                                                          gp_Pnt2d& valueP1,
@@ -1432,7 +2123,7 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
     if (mp[0] > 0 && 
        ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
       p1OnIso = Standard_True;
-      mind2[0] = currd2[0];
+      mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file 
       if (isoU) valueP1.SetCoord(isoVal, tp[0]);
       else      valueP1.SetCoord(tp[0], isoVal);
     }
@@ -1514,18 +2205,17 @@ Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(cons
   }  */
   return isoParam;
   }  // RAJOUT
-  catch(Standard_Failure) {
-//  pb : on affiche ce qu on peut
+  catch(Standard_Failure const& anException) {
 #ifdef OCCT_DEBUG
+//  pb : on affiche ce qu on peut
     for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
-      cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
-       points(numpnt).X()<<"  "<<points(numpnt).Y()<<"  "<<points(numpnt).Z()<<")"<<endl;
+      std::cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
+       points(numpnt).X()<<"  "<<points(numpnt).Y()<<"  "<<points(numpnt).Z()<<")"<<std::endl;
     }
+    std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
+    anException.Print(std::cout); std::cout << std::endl;
 #endif
-#ifdef OCCT_DEBUG //:s5
-    cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
-    Standard_Failure::Caught()->Print(cout); cout << endl;
-#endif
+    (void)anException;
     return Standard_False;
   }
 }