0026196: Wrong result obtained by projection algorithm.
[occt.git] / src / ProjLib / ProjLib_ComputeApproxOnPolarSurface.cxx
index 347fd4a..07d2aef 100644 (file)
@@ -14,7 +14,7 @@
 // commercial license or contractual agreement.
 
 #include <ProjLib_ComputeApproxOnPolarSurface.hxx>
-#include <AppCont_Function2d.hxx>
+#include <AppCont_Function.hxx>
 #include <ElSLib.hxx>
 #include <ElCLib.hxx>
 #include <BSplCLib.hxx>
 //#include <GeomLib_IsIso.hxx>
 //#include <GeomLib_CheckSameParameter.hxx>
 
-#ifdef DEB
+#ifdef OCCT_DEBUG
 #ifdef DRAW
 #include <DrawTrSurf.hxx>
+#include <Geom2d_Curve.hxx>
 #endif
 //static Standard_Integer compteur = 0;
 #endif
 
+struct aFuncStruct
+{
+  Handle(Adaptor3d_HSurface) mySurf; // Surface where to project.
+  Handle(Adaptor3d_HCurve)   myCurve; // Curve to project.
+  Handle(Adaptor2d_HCurve2d) myInitCurve2d; // Initial 2dcurve projection.
+  Standard_Real mySqProjOrtTol; // Used to filter non-orthogonal projected point.
+  Standard_Real myTolU;
+  Standard_Real myTolV;
+  Standard_Real myPeriod[2]; // U and V period correspondingly.
+};
+
+//=======================================================================
+//function : computePeriodicity
+//purpose  : Compute period information on adaptor.
+//=======================================================================
+static void computePeriodicity(const Handle(Adaptor3d_HSurface)& theSurf,
+                               Standard_Real &theUPeriod,
+                               Standard_Real &theVPeriod)
+{
+  theUPeriod = 0.0;
+  theVPeriod = 0.0;
+
+  // Compute once information about periodicity.
+  // Param space may be reduced in case of rectangular trimmed surface,
+  // in this case really trimmed bounds should be set as unperiodic.
+  Standard_Real aTrimF, aTrimL, aBaseF, aBaseL, aDummyF, aDummyL;
+  Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theSurf->Surface(), Standard_False); // Not trim.
+  // U param space.
+  if (theSurf->IsUPeriodic())
+  {
+    theUPeriod = theSurf->UPeriod();
+  }
+  else if(theSurf->IsUClosed())
+  {
+    theUPeriod = theSurf->LastUParameter() - theSurf->FirstUParameter();
+  }
+  if (theUPeriod != 0.0)
+  {
+    aTrimF = theSurf->FirstUParameter(); // Trimmed first
+    aTrimL = theSurf->LastUParameter(); // Trimmed last
+    aS->Bounds(aBaseF, aBaseL, aDummyF, aDummyL); // Non-trimmed values.
+    if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
+    {
+      // Param space reduced.
+      theUPeriod = 0.0;
+    }
+  }
+
+  // V param space.
+  if (theSurf->IsVPeriodic())
+  {
+    theVPeriod = theSurf->VPeriod();
+  }
+  else if(theSurf->IsVClosed())
+  {
+    theVPeriod = theSurf->LastVParameter() - theSurf->FirstVParameter();
+  }
+  if (theVPeriod != 0.0)
+  {
+    aTrimF = theSurf->FirstVParameter(); // Trimmed first
+    aTrimL = theSurf->LastVParameter(); // Trimmed last
+    aS->Bounds(aDummyF, aDummyL, aBaseF, aBaseL); // Non-trimmed values.
+    if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
+    {
+      // Param space reduced.
+      theVPeriod = 0.0;
+    }
+  }
+}
+
+//=======================================================================
+//function : aFuncValue
+//purpose  : compute functional value in (theU,theV) point
+//=======================================================================
+static Standard_Real anOrthogSqValue(const gp_Pnt& aBasePnt,
+                                     const Handle(Adaptor3d_HSurface)& Surf,
+                                     const Standard_Real theU,
+                                     const Standard_Real theV)
+{
+  // Since find projection, formula is:
+  // F1 = Dot(S_U, Vec(aBasePnt, aProjPnt))
+  // F2 = Dot(S_V, Vec(aBasePnt, aProjPnt))
+
+  gp_Pnt aProjPnt;
+  gp_Vec aSu, aSv;
+
+  Surf->D1(theU, theV, aProjPnt, aSu, aSv);
+  gp_Vec aBaseVec(aBasePnt, aProjPnt);
+
+  if (aSu.SquareMagnitude() > Precision::SquareConfusion())
+    aSu.Normalize();
+
+  if (aSv.SquareMagnitude() > Precision::SquareConfusion())
+    aSv.Normalize();
+
+  Standard_Real aFirstPart = aSu.Dot(aBaseVec);
+  Standard_Real aSecondPart = aSv.Dot(aBaseVec);
+  return (aFirstPart * aFirstPart + aSecondPart * aSecondPart);
+}
+
 //=======================================================================
 //function : Value
 //purpose  : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
 //            <InitCurve2d> use for calculate start 2D point.
 //=======================================================================
-
-static gp_Pnt2d Function_Value(const Standard_Real U,
-                              const Handle(Adaptor3d_HSurface)& Surf,    
-                              const Handle(Adaptor3d_HCurve)& Curve,
-                              const Handle(Adaptor2d_HCurve2d)& InitCurve2d,
-                               //OCC217
-                              const Standard_Real DistTol3d, const Standard_Real tolU, const Standard_Real tolV)
-                              //const Standard_Real Tolerance)
+static gp_Pnt2d Function_Value(const Standard_Real theU,
+                               const aFuncStruct& theData)
 {
-  //OCC217
-  //Standard_Real Tol3d = 100*Tolerance;
+  gp_Pnt2d p2d = theData.myInitCurve2d->Value(theU) ;
+  gp_Pnt p = theData.myCurve->Value(theU);
+  gp_Pnt aSurfPnt = theData.mySurf->Value(p2d.X(), p2d.Y());
+  Standard_Real aSurfPntDist = aSurfPnt.SquareDistance(p);
 
-  gp_Pnt2d p2d = InitCurve2d->Value(U) ;
-  gp_Pnt p = Curve->Value(U);
-//  Curve->D0(U,p) ;
   Standard_Real Uinf, Usup, Vinf, Vsup;
-  Uinf = Surf->Surface().FirstUParameter();
-  Usup = Surf->Surface().LastUParameter();
-  Vinf = Surf->Surface().FirstVParameter();
-  Vsup = Surf->Surface().LastVParameter();
+  Uinf = theData.mySurf->Surface().FirstUParameter();
+  Usup = theData.mySurf->Surface().LastUParameter();
+  Vinf = theData.mySurf->Surface().FirstVParameter();
+  Vsup = theData.mySurf->Surface().LastVParameter();
+
+  // Check case when curve is close to co-parametrized isoline on surf.
+  if (Abs (p2d.X() - Uinf) < Precision::PConfusion() ||
+      Abs (p2d.X() - Usup) < Precision::PConfusion() )
+  {
+    // V isoline.
+    gp_Pnt aPnt;
+    theData.mySurf->D0(p2d.X(), theU, aPnt);
+    if (aPnt.SquareDistance(p) < aSurfPntDist)
+      p2d.SetY(theU);
+  }
+
+  if (Abs (p2d.Y() - Vinf) < Precision::PConfusion() ||
+      Abs (p2d.Y() - Vsup) < Precision::PConfusion() )
+  {
+    // U isoline.
+    gp_Pnt aPnt;
+    theData.mySurf->D0(theU, p2d.Y(), aPnt);
+    if (aPnt.SquareDistance(p) < aSurfPntDist)
+      p2d.SetX(theU);
+  }
+
   Standard_Integer decalU = 0, decalV = 0;
   Standard_Real U0 = p2d.X(), V0 = p2d.Y();
-  
-  GeomAbs_SurfaceType Type = Surf->GetType();
+
+  GeomAbs_SurfaceType Type = theData.mySurf->GetType();
   if((Type != GeomAbs_BSplineSurface) && 
      (Type != GeomAbs_BezierSurface)  &&
-     (Type != GeomAbs_OffsetSurface)    ) {
+     (Type != GeomAbs_OffsetSurface)    )
+  {
+    // Analytical cases.
     Standard_Real S = 0., T = 0.;
-    switch (Type) {
-//    case GeomAbs_Plane:
-//      {
-//     gp_Pln Plane = Surf->Plane();
-//     ElSLib::Parameters( Plane, p, S, T);
-//     break;
-//      }
+    switch (Type)
+    {
     case GeomAbs_Cylinder:
       {
-       gp_Cylinder Cylinder = Surf->Cylinder();
-       ElSLib::Parameters( Cylinder, p, S, T);
-       if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
-       if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
-       S += decalU*2*M_PI;
-       break;
+        gp_Cylinder Cylinder = theData.mySurf->Cylinder();
+        ElSLib::Parameters( Cylinder, p, S, T);
+        if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
+        if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
+        S += decalU*2*M_PI;
+        break;
       }
     case GeomAbs_Cone:
       {
-       gp_Cone Cone = Surf->Cone();
-       ElSLib::Parameters( Cone, p, S, T);
-       if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
-       if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
-       S += decalU*2*M_PI;
-       break;
+        gp_Cone Cone = theData.mySurf->Cone();
+        ElSLib::Parameters( Cone, p, S, T);
+        if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
+        if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
+        S += decalU*2*M_PI;
+        break;
       }
     case GeomAbs_Sphere:
       {
-       gp_Sphere Sphere = Surf->Sphere();
-       ElSLib::Parameters( Sphere, p, S, T);
-       if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
-       if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
-       S += decalU*2*M_PI;
-       if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
-       if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
-       T += decalV*2*M_PI;
-       if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI) {
-         T = M_PI - T;
-         if(U0 < S)
-           S -= M_PI;
-         else
-           S += M_PI;
-       }
-       break;
+        gp_Sphere Sphere = theData.mySurf->Sphere();
+        ElSLib::Parameters( Sphere, p, S, T);
+        if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
+        if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
+        S += decalU*2*M_PI;
+        if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
+        if(V0 > (Vsup+(Vsup-Vinf))) decalV =  int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
+        T += decalV*2*M_PI;
+        if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI)
+        {
+          T = M_PI - T;
+          if(U0 < S)
+            S -= M_PI;
+          else
+            S += M_PI;
+        }
+        break;
       }
     case GeomAbs_Torus:
       {
-       gp_Torus Torus = Surf->Torus();
-       ElSLib::Parameters( Torus, p, S, T);
-       if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
-       if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
-       if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
-       if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
-       S += decalU*2*M_PI; T += decalV*2*M_PI;
-       break;
+        gp_Torus Torus = theData.mySurf->Torus();
+        ElSLib::Parameters( Torus, p, S, T);
+        if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
+        if(U0 > Usup) decalU =  int((U0 - Usup)/(2*M_PI))+1;
+        if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
+        if(V0 > Vsup) decalV =  int((V0 - Vsup)/(2*M_PI))+1;
+        S += decalU*2*M_PI; T += decalV*2*M_PI;
+        break;
       }
     default:
       Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::Value");
     }
     return gp_Pnt2d(S, T);
   }
-  
-  //////////////////
+
+  // Non-analytical case.
   Standard_Real Dist2Min = RealLast();
-  //OCC217
-  //Standard_Real tolU,tolV ;
-  //tolU = Tolerance;
-  //tolV = Tolerance;
-  
-  Standard_Real uperiod =0, vperiod = 0, u, v;
-  // U0 and V0 are the points within the initialized period 
-  // (periode with u and v),
-  // U1 and V1 are the points for construction of tops
-  
-  if(Surf->IsUPeriodic() || Surf->IsUClosed()) {
-    uperiod =  Surf->LastUParameter() - Surf->FirstUParameter();
-  }
-  if(Surf->IsVPeriodic() || Surf->IsVClosed()) {
-    vperiod = Surf->LastVParameter() - Surf->FirstVParameter();
-  } 
-  if(U0 < Uinf) {
+  Standard_Real uperiod = theData.myPeriod[0],
+                vperiod = theData.myPeriod[1],
+                u, v;
+
+  // U0 and V0 are the points within the initialized period.
+  if(U0 < Uinf)
+  {
     if(!uperiod)
       U0 = Uinf;
-    else {
+    else
+    {
       decalU = int((Uinf - U0)/uperiod)+1;
       U0 += decalU*uperiod;
     }
   }
-  if(U0 > Usup) {
+  if(U0 > Usup)
+  {
     if(!uperiod)
       U0 = Usup;
-    else {
+    else
+    {
       decalU = -(int((U0 - Usup)/uperiod)+1);
       U0 += decalU*uperiod;
     }
   }
-  if(V0 < Vinf) {
+  if(V0 < Vinf)
+  {
     if(!vperiod)
       V0 = Vinf;
-    else {
+    else
+    {
       decalV = int((Vinf - V0)/vperiod)+1;
       V0 += decalV*vperiod;
     }
   }
-  if(V0 > Vsup) {
+  if(V0 > Vsup)
+  {
     if(!vperiod)
       V0 = Vsup;
-    else {
+    else
+    {
       decalV = -int((V0 - Vsup)/vperiod)-1;
       V0 += decalV*vperiod;
     }
   }
-  
-  // The surface around U0 is reduced
+
+  // The surface around (U0,V0) is reduced.
   Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
   Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
   if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
   if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
   if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
   if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
-  
-  //  const Adaptor3d_Surface GAS = Surf->Surface();
 
   GeomAdaptor_Surface SurfLittle;
-  if (Type == GeomAbs_BSplineSurface) {
-    Handle(Geom_Surface) GBSS(Surf->Surface().BSpline());
+  if (Type == GeomAbs_BSplineSurface)
+  {
+    Handle(Geom_Surface) GBSS(theData.mySurf->Surface().BSpline());
     SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
   }
-  else if (Type == GeomAbs_BezierSurface) {
-    Handle(Geom_Surface) GS(Surf->Surface().Bezier());
+  else if (Type == GeomAbs_BezierSurface)
+  {
+    Handle(Geom_Surface) GS(theData.mySurf->Surface().Bezier());
     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
   }
-  else if (Type == GeomAbs_OffsetSurface) {
-    Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(Surf->Surface());
+  else if (Type == GeomAbs_OffsetSurface)
+  {
+    Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(theData.mySurf->Surface());
     SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
   }
-  else {
+  else
+  {
     Standard_NoSuchObject::Raise("");
   }
 
-  Extrema_GenLocateExtPS  locext(p, SurfLittle, U0, V0, tolU, tolV);
-  if (locext.IsDone()) {
-    Dist2Min = locext.SquareDistance();
-    if (Dist2Min < DistTol3d * DistTol3d) {
-      (locext.Point()).Parameter(u, v);
+  // Try to run simple search with initial point (U0, V0).
+  Extrema_GenLocateExtPS  locext(p, SurfLittle, U0, V0, theData.myTolU, theData.myTolV);
+  if (locext.IsDone()) 
+  {
+    locext.Point().Parameter(u, v);
+    Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
+    if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
+        locext.SquareDistance() < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
+    {
       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
       return pnt;
     }
-  } 
-  
-  Extrema_ExtPS  ext(p, SurfLittle, tolU, tolV) ;
-  if (ext.IsDone() && ext.NbExt()>=1 ) {
+  }
+
+  // Perform whole param space search.
+  Extrema_ExtPS  ext(p, SurfLittle, theData.myTolU, theData.myTolV);
+  if (ext.IsDone() && ext.NbExt() >= 1)
+  {
     Dist2Min = ext.SquareDistance(1);
     Standard_Integer GoodValue = 1;
-    for ( Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ ) 
-      if( Dist2Min > ext.SquareDistance(i)) {
-       Dist2Min = ext.SquareDistance(i);
-       GoodValue = i;
+    for (Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
+    {
+      if( Dist2Min > ext.SquareDistance(i))
+      {
+        Dist2Min = ext.SquareDistance(i);
+        GoodValue = i;
       }
-    if (Dist2Min < DistTol3d * DistTol3d) {
-      (ext.Point(GoodValue)).Parameter(u,v);
+    }
+    ext.Point(GoodValue).Parameter(u, v);
+    Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
+    if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
+        ext.SquareDistance(GoodValue) < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
+    {
       gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
       return pnt;
     }
   }
-  
+
+  // Both searches return bad values, use point from initial 2dcurve.
   return p2d;
 }
 
@@ -288,49 +414,58 @@ static gp_Pnt2d Function_Value(const Standard_Real U,
 //purpose  : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
 //=======================================================================
 
-class ProjLib_PolarFunction : public AppCont_Function2d
+class ProjLib_PolarFunction : public AppCont_Function
 {
-  Handle(Adaptor3d_HCurve)           myCurve;
-  Handle(Adaptor2d_HCurve2d)         myInitialCurve2d ;
-  Handle(Adaptor3d_HSurface)         mySurface ;
-  //OCC217
-  Standard_Real                      myTolU, myTolV; 
-  Standard_Real                      myDistTol3d; 
-  //Standard_Real                    myTolerance ; 
-  
+  aFuncStruct myStruct;
+
   public :
-  
+
   ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C, 
-                       const Handle(Adaptor3d_HSurface)& Surf,
-                       const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
-                       //OCC217
-                       const Standard_Real Tol3d) :
-                       //const Standard_Real Tolerance) :
-  myCurve(C),
-  myInitialCurve2d(InitialCurve2d),
-  mySurface(Surf),
-  //OCC217             
-  myTolU(Surf->UResolution(Tol3d)),
-  myTolV(Surf->VResolution(Tol3d)),
-  myDistTol3d(100.0*Tol3d){} ;
-  //myTolerance(Tolerance){} ;
-  
+                        const Handle(Adaptor3d_HSurface)& Surf,
+                        const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
+                        const Standard_Real Tol3d)
+  {
+    myNbPnt = 0;
+    myNbPnt2d = 1;
+
+    computePeriodicity(Surf, myStruct.myPeriod[0], myStruct.myPeriod[1]);
+
+    myStruct.myCurve = C;
+    myStruct.myInitCurve2d = InitialCurve2d;
+    myStruct.mySurf = Surf;
+    myStruct.mySqProjOrtTol = 10000.0 * Tol3d * Tol3d;
+    myStruct.myTolU = Surf->UResolution(Tol3d);
+    myStruct.myTolV = Surf->VResolution(Tol3d);
+  }
+
   ~ProjLib_PolarFunction() {}
-  
+
   Standard_Real FirstParameter() const
-  {return (myCurve->FirstParameter()/*+1.e-9*/);}
-  
+  {
+    return myStruct.myCurve->FirstParameter();
+  }
+
   Standard_Real LastParameter() const
-  {return (myCurve->LastParameter()/*-1.e-9*/);}
-  
-  gp_Pnt2d Value(const Standard_Real t) const {
-    return Function_Value
-      (t,mySurface,myCurve,myInitialCurve2d,myDistTol3d,myTolU,myTolV) ;  //OCC217
-      //(t,mySurface,myCurve,myInitialCurve2d,myTolerance) ;
+  {
+    return myStruct.myCurve->LastParameter();
   }
-  
-//  Standard_Boolean D1(const Standard_Real t, gp_Pnt2d& P, gp_Vec2d& V) const
-  Standard_Boolean D1(const Standard_Real , gp_Pnt2d& , gp_Vec2d& ) const
+
+  gp_Pnt2d Value(const Standard_Real t) const
+  {
+    return Function_Value(t, myStruct);
+  }
+
+  Standard_Boolean Value(const Standard_Real   theT,
+                         NCollection_Array1<gp_Pnt2d>& thePnt2d,
+                         NCollection_Array1<gp_Pnt>&   /*thePnt*/) const
+  {
+    thePnt2d(1) = Function_Value(theT, myStruct);
+    return Standard_True;
+  }
+
+  Standard_Boolean D1(const Standard_Real   /*theT*/,
+                      NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
+                      NCollection_Array1<gp_Vec>&   /*theVec*/) const
     {return Standard_False;}
 };
 
@@ -339,7 +474,11 @@ class ProjLib_PolarFunction : public AppCont_Function2d
 //purpose  : 
 //=======================================================================
 
-ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface(){}
+ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
+: myProjIsDone(Standard_False),
+  myTolerance (-1.0)
+{
+}
 
 
 //=======================================================================
@@ -348,15 +487,14 @@ ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface(){}
 //=======================================================================
 
 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
-(const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
- const Handle(Adaptor3d_HCurve)& Curve,
- const Handle(Adaptor3d_HSurface)& S,
- const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
- //const Standard_Real tolerance):myProjIsDone(Standard_False)
+                    (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
+                     const Handle(Adaptor3d_HCurve)&   theCurve,
+                     const Handle(Adaptor3d_HSurface)& theSurface,
+                     const Standard_Real               theTolerance3D)
+: myProjIsDone(Standard_False),
+  myTolerance (theTolerance3D)
 {
-  myTolerance = tol3d; //OCC217
-  //myTolerance = Max(Tolerance,Precision::PApproximation());
-  myBSpline = Perform(InitialCurve2d,Curve,S);
+  myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
 }
 //=======================================================================
 //function : ProjLib_ComputeApproxOnPolarSurface
@@ -364,21 +502,23 @@ ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 //=======================================================================
 
 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
-(const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
- const Handle(Adaptor2d_HCurve2d)& InitialCurve2dBis,
- const Handle(Adaptor3d_HCurve)& Curve,
- const Handle(Adaptor3d_HSurface)& S,
- const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
- //const Standard_Real tolerance):myProjIsDone(Standard_False)
-{// InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing 
-  myTolerance = tol3d; //OCC217
-  //myTolerance = Max(tolerance,Precision::PApproximation());
-  Handle(Geom2d_BSplineCurve) bsc = Perform(InitialCurve2d,Curve,S);
+                (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
+                 const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
+                 const Handle(Adaptor3d_HCurve)&   theCurve,
+                 const Handle(Adaptor3d_HSurface)& theSurface,
+                 const Standard_Real               theTolerance3D)
+: myProjIsDone(Standard_False),
+  myTolerance (theTolerance3D)
+{
+  // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing 
+  Handle(Geom2d_BSplineCurve) bsc =
+    Perform(theInitialCurve2d, theCurve, theSurface);
+
   if(myProjIsDone) {
     gp_Pnt2d P2dproj, P2d, P2dBis;
     P2dproj = bsc->StartPoint();
-    P2d    = InitialCurve2d->Value(InitialCurve2d->FirstParameter());
-    P2dBis = InitialCurve2dBis->Value(InitialCurve2dBis->FirstParameter());
+    P2d    = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
+    P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
 
     Standard_Real Dist, DistBis;
     Dist    = P2dproj.Distance(P2d);
@@ -403,19 +543,25 @@ ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
 //=======================================================================
 
 ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
-(const Handle(Adaptor3d_HCurve)& Curve,
- const Handle(Adaptor3d_HSurface)& S,
- const Standard_Real tol3d):myProjIsDone(Standard_False)  //OCC217
- //const Standard_Real tolerance):myProjIsDone(Standard_False)
+                      (const Handle(Adaptor3d_HCurve)&   theCurve,
+                       const Handle(Adaptor3d_HSurface)& theSurface,
+                       const Standard_Real               theTolerance3D)
+: myProjIsDone(Standard_False),
+  myTolerance (theTolerance3D)
 {
-  myTolerance = tol3d; //OCC217
-  //myTolerance = Max(tolerance,Precision::PApproximation());
-  const Handle_Adaptor2d_HCurve2d InitCurve2d ;
-  myBSpline = Perform(InitCurve2d,Curve,S);  
+  const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
+  myBSpline = Perform(anInitCurve2d, theCurve, theSurface);  
 } 
 
+//=======================================================================
+//function : Concat
+//purpose  : 
+//=======================================================================
+
 static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
-                                         Handle(Geom2d_BSplineCurve) C2)
+                                          Handle(Geom2d_BSplineCurve) C2,
+                                          Standard_Real theUJump,
+                                          Standard_Real theVJump)
 {
   Standard_Integer deg, deg1, deg2;
   deg1 = C1->Degree();
@@ -472,7 +618,8 @@ static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
   }
   for (i = 2; i <= np2; i++) {
     count++;
-    P(count) = P2(i);
+    P(count).SetX(P2(i).X() + theUJump);
+    P(count).SetY(P2(i).Y() + theVJump);
   }
 
   Handle(Geom2d_BSplineCurve) BS = 
@@ -664,15 +811,38 @@ Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
        }
       }
 
+      Standard_Real anUPeriod, anVPeriod;
+      computePeriodicity(S, anUPeriod, anVPeriod);
       Standard_Integer NbC = LOfBSpline2d.Extent();
       Handle(Geom2d_BSplineCurve) CurBS;
       CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
       LOfBSpline2d.RemoveFirst();
-      for (Standard_Integer ii = 2; ii <= NbC; ii++) {
-       Handle(Geom2d_BSplineCurve) BS = 
-         Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
-       CurBS = Concat(CurBS,BS);
-       LOfBSpline2d.RemoveFirst();
+      for (Standard_Integer ii = 2; ii <= NbC; ii++)
+      {
+        Handle(Geom2d_BSplineCurve) BS = 
+          Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
+
+        //Check for period jump in point of contact.
+        gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
+        gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
+        Standard_Real anUJump = 0.0, anVJump = 0.0;
+
+        if (anUPeriod > 0.0 &&
+            Abs (aC1End.X() - aC2Beg.X()) > (anUPeriod ) / 2.01)
+        {
+          Standard_Real aMultCoeff =  aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
+          anUJump = (anUPeriod) * aMultCoeff;
+        }
+
+        if (anVPeriod &&
+            Abs (aC1End.Y() - aC2Beg.Y()) > (anVPeriod) / 2.01)
+        {
+          Standard_Real aMultCoeff =  aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
+          anVJump = (anVPeriod) * aMultCoeff;
+        }
+
+        CurBS = Concat(CurBS,BS, anUJump, anVJump);
+        LOfBSpline2d.RemoveFirst();
       }
       return CurBS;
     }
@@ -705,14 +875,9 @@ Handle(Adaptor2d_HCurve2d)
   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
   Standard_Real DistTol3d = 100.0*Tol3d;
 
-  Standard_Real uperiod = 0., vperiod = 0.;
-  if(Surf->IsUPeriodic() || Surf->IsUClosed())
-    uperiod = Surf->LastUParameter() - Surf->FirstUParameter(); 
-  
-  if(Surf->IsVPeriodic() || Surf->IsVClosed())
-    vperiod = Surf->LastVParameter() - Surf->FirstVParameter(); 
+  Standard_Real uperiod = 0.0, vperiod = 0.0;
+  computePeriodicity(Surf, uperiod, vperiod);
 
-  
   // NO myTol is Tol2d !!!!
   //Standard_Real TolU = myTolerance, TolV = myTolerance;
   //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
@@ -1106,18 +1271,50 @@ Handle(Adaptor2d_HCurve2d)
                   V1 = V0 + vsens*vperiod;
                   Pts2d(i).SetCoord(U1,V1);
                   myProjIsDone = Standard_True;
+
+                  if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
+                  {//Make 1st point more precise for periodic surfaces
+                    const Standard_Integer aSize = 3;
+                    const gp_Pnt2d aP(Pts2d(2)); 
+                    Standard_Real aUpar[aSize], aVpar[aSize];
+                    Pts2d(1).Coord(aUpar[1], aVpar[1]);
+                    aUpar[0] = aUpar[1] - uperiod;
+                    aUpar[2] = aUpar[1] + uperiod;
+                    aVpar[0] = aVpar[1] - vperiod;
+                    aVpar[2] = aVpar[1] + vperiod;
+
+                    Standard_Real aSQdistMin = RealLast();
+                    Standard_Integer aBestUInd = 1, aBestVInd = 1;
+                    const Standard_Integer  aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
+                                            aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
+                    for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
+                    {
+                      for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
+                      {
+                        Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
+                        if(aSQdist < aSQdistMin)
+                        {
+                          aSQdistMin = aSQdist;
+                          aBestUInd = uInd;
+                          aBestVInd = vInd;
+                        }
+                      }
+                    }
+
+                    Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
+                  }//if(i == 2) condition
                 }
               }
             }
           }
          if(!myProjIsDone && uperiod) {
-           Standard_Real Uinf, Usup, Uaux;
-           Uinf = Surf->Surface().FirstUParameter();
-           Usup = Surf->Surface().LastUParameter();
-           if((Usup - U0) > (U0 - Uinf)) 
-             Uaux = 2*Uinf - U0 + uperiod;
+           Standard_Real aUinf, aUsup, Uaux;
+           aUinf = Surf->Surface().FirstUParameter();
+           aUsup = Surf->Surface().LastUParameter();
+           if((aUsup - U0) > (U0 - aUinf)) 
+             Uaux = 2*aUinf - U0 + uperiod;
            else 
-             Uaux = 2*Usup - U0 - uperiod;
+             Uaux = 2*aUsup - U0 - uperiod;
            Extrema_GenLocateExtPS  locext(pntproj, 
                                           Surf->Surface(), 
                                           Uaux, V0, TolU, TolV);
@@ -1125,7 +1322,7 @@ Handle(Adaptor2d_HCurve2d)
              if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
              //if (locext.SquareDistance() < Tol3d * Tol3d) {
                (locext.Point()).Parameter(u,v);
-               if((Usup - U0) > (U0 - Uinf)) 
+               if((aUsup - U0) > (U0 - aUinf)) 
                  usens--;
                else 
                  usens++;
@@ -1137,13 +1334,13 @@ Handle(Adaptor2d_HCurve2d)
              }
          }
          if(!myProjIsDone && vperiod) {
-           Standard_Real Vinf, Vsup, Vaux;
-           Vinf = Surf->Surface().FirstVParameter();
-           Vsup = Surf->Surface().LastVParameter();
-           if((Vsup - V0) > (V0 - Vinf)) 
-             Vaux = 2*Vinf - V0 + vperiod;
+           Standard_Real aVinf, aVsup, Vaux;
+           aVinf = Surf->Surface().FirstVParameter();
+           aVsup = Surf->Surface().LastVParameter();
+           if((aVsup - V0) > (V0 - aVinf)) 
+             Vaux = 2*aVinf - V0 + vperiod;
            else 
-             Vaux = 2*Vsup - V0 - vperiod;
+             Vaux = 2*aVsup - V0 - vperiod;
            Extrema_GenLocateExtPS  locext(pntproj, 
                                           Surf->Surface(), 
                                           U0, Vaux, TolU, TolV) ;
@@ -1151,7 +1348,7 @@ Handle(Adaptor2d_HCurve2d)
              if (locext.SquareDistance() < DistTol3d * DistTol3d) {  //OCC217
              //if (locext.SquareDistance() < Tol3d * Tol3d) {
                (locext.Point()).Parameter(u,v);
-               if((Vsup - V0) > (V0 - Vinf)) 
+               if((aVsup - V0) > (V0 - aVinf)) 
                  vsens--;
                else 
                  vsens++;
@@ -1198,15 +1395,15 @@ Handle(Adaptor2d_HCurve2d)
            Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
            if (ext.IsDone()) {
              Dist2Min = ext.SquareDistance(1);
-             Standard_Integer GoodValue = 1;
+             Standard_Integer aGoodValue = 1;
              for ( j = 2 ; j <= ext.NbExt() ; j++ )
                if( Dist2Min > ext.SquareDistance(j)) {
                  Dist2Min = ext.SquareDistance(j);
-                 GoodValue = j;
+                 aGoodValue = j;
                }
              if (Dist2Min < DistTol3d * DistTol3d) {
              //if (Dist2Min < Tol3d * Tol3d) {
-               (ext.Point(GoodValue)).Parameter(u,v);
+               (ext.Point(aGoodValue)).Parameter(u,v);
                if(uperiod) {
                  if((U0 - u) > (2*uperiod/3)) {
                    usens++;
@@ -1270,7 +1467,7 @@ Handle(Adaptor2d_HCurve2d)
     //////////////////////////////////////////
     Geom2dAdaptor_Curve GAC(myBSpline);
     Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
-#ifdef DEB
+#ifdef OCCT_DEBUG
 //    char name [100];
 //    sprintf(name,"%s_%d","build",compteur++);
 //    DrawTrSurf::Set(name,myBSpline);
@@ -1306,7 +1503,7 @@ Handle(Geom2d_BSplineCurve)
   Standard_Real Tol3d = myTolerance;
   Standard_Real DistTol3d = 1.0*Tol3d;
   Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
-  Standard_Real Tol2d = Sqrt(TolU*TolU + TolV*TolV);
+  Standard_Real Tol2d = Max(Sqrt(TolU*TolU + TolV*TolV), Precision::PConfusion());
 
   Standard_Integer i;
   GeomAbs_SurfaceType TheTypeS = Surf->GetType();
@@ -1566,7 +1763,7 @@ Handle(Geom2d_BSplineCurve)
   ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ;  //OCC217
   //ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, myTolerance) ;
 
-#ifdef DEB
+#ifdef OCCT_DEBUG
   Standard_Integer Nb = 50;
   
   Standard_Real U, U1, U2;
@@ -1613,14 +1810,14 @@ Handle(Geom2d_BSplineCurve)
                            Standard_True);
 
   if(Fit.IsAllApproximated()) {
-    Standard_Integer i;
+    Standard_Integer j;
     Standard_Integer NbCurves = Fit.NbMultiCurves();
     Standard_Integer MaxDeg = 0;
     // To transform the MultiCurve into BSpline, it is required that all  
     // Bezier constituing it have the same degree -> Calculation of MaxDeg
     Standard_Integer NbPoles  = 1;
-    for (i = 1; i <= NbCurves; i++) {
-      Standard_Integer Deg = Fit.Value(i).Degree();
+    for (j = 1; j <= NbCurves; j++) {
+      Standard_Integer Deg = Fit.Value(j).Degree();
       MaxDeg = Max ( MaxDeg, Deg);
     }
 
@@ -1645,15 +1842,15 @@ Handle(Geom2d_BSplineCurve)
        BSplCLib::IncreaseDegree( MaxDeg, Poles2d, PLib::NoWeights(), 
                         TempPoles, PLib::NoWeights());
        //update of tops of the PCurve
-       for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
-         Poles.SetValue( Compt, TempPoles( j));
+       for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
+         Poles.SetValue( Compt, TempPoles( k));
          Compt++;
        }
       }
       else {
        //update of tops of the PCurve
-       for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
-         Poles.SetValue( Compt, Poles2d( j));
+       for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
+         Poles.SetValue( Compt, Poles2d( k));
          Compt++;
        }
       } 
@@ -1663,12 +1860,7 @@ Handle(Geom2d_BSplineCurve)
     
     //update of fields of  ProjLib_Approx
     Standard_Integer NbKnots = NbCurves + 1;
-    
-    // The start and end nodes are not correct : Cf: opening of the interval
-    //Knots( 1) -= 1.e-9;
-    //Knots(NbKnots) += 1.e-9; 
-    
-    
+
     TColStd_Array1OfInteger   Mults( 1, NbKnots);
     Mults.Init(MaxDeg);
     Mults.SetValue( 1, MaxDeg + 1);
@@ -1685,7 +1877,7 @@ Handle(Geom2d_BSplineCurve)
       OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,Tol3d);  //OCC217
       //OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,myTolerance);
     }
-#ifdef DEB
+#ifdef OCCT_DEBUG
     if (!OK) {
       cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<endl;
     }