0027531: Modeling Algorithms - Make the algorithm Approx_SameParameter more clear...
[occt.git] / src / Approx / Approx_SameParameter.cxx
index 54a8b56..7258312 100644 (file)
@@ -14,7 +14,6 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-//  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898
 
 #include <Adaptor2d_HCurve2d.hxx>
 #include <Adaptor3d_CurveOnSurface.hxx>
 #include <BSplCLib.hxx>
 #include <Extrema_ExtPC.hxx>
 #include <Extrema_LocateExtPC.hxx>
-#include <GCPnts_QuasiUniformDeflection.hxx>
 #include <Geom2d_BSplineCurve.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom2dAdaptor_Curve.hxx>
 #include <Geom2dAdaptor_HCurve.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Surface.hxx>
-#include <GeomAdaptor_Curve.hxx>
 #include <GeomAdaptor_HCurve.hxx>
 #include <GeomAdaptor_HSurface.hxx>
-#include <GeomAdaptor_Surface.hxx>
 #include <GeomLib_MakeCurvefromApprox.hxx>
 #include <Precision.hxx>
 #include <Standard_ConstructionError.hxx>
 #include <Standard_OutOfRange.hxx>
 #include <TColStd_Array1OfReal.hxx>
+#include <Geom2dAdaptor.hxx>
 
 //=======================================================================
 //class    : Approx_SameParameter_Evaluator
@@ -125,7 +122,7 @@ static void ProjectPointOnCurve(const Standard_Real      InitValue,
                                 Standard_Boolean&        Status,
                                 Standard_Real&           Result)
 {
-  Standard_Integer num_iter = 0, not_done = 1, ii;
+  Standard_Integer num_iter = 0, not_done = 1;
 
   gp_Pnt a_point;
   gp_Vec vector, d1, d2;
@@ -136,8 +133,7 @@ static void ProjectPointOnCurve(const Standard_Real      InitValue,
   {
     num_iter++;
     Curve.D2(param, a_point, d1, d2);
-    for (ii = 1 ; ii <= 3 ; ii++) 
-      vector.SetCoord(ii, APoint.Coord(ii) - a_point.Coord(ii));
+    vector = gp_Vec(a_point,APoint);
 
     func = vector.Dot(d1);
     if ( Abs(func) < Tolerance * d1.Magnitude())
@@ -189,7 +185,7 @@ static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
     d2 = Max(d2, Pc3d.SquareDistance(Pcons));
   }
 
-  const Standard_Real aMult = 1.5; // To be tolerant to discrete tolerance computing.
+  const Standard_Real aMult = 1. + 0.05; // 
   Standard_Real aDeviation = aMult * sqrt(d2);
   aDeviation = Max(aDeviation, Precision::Confusion()); // Tolerance in modeling space.
   return aDeviation;
@@ -199,11 +195,10 @@ static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
 //function : Check
 //purpose  : Check current interpolation for validity.
 //=======================================================================
-static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots, 
+static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
                               const TColStd_Array1OfReal& Poles,
                               const Standard_Integer nbp,
-                              const TColStd_Array1OfReal& pc3d,
-                              const TColStd_Array1OfReal& ,
+                              const Standard_Real *pc3d,
                               const Handle(Adaptor3d_HCurve)& c3d,
                               const Adaptor3d_CurveOnSurface& cons,
                               Standard_Real& tol,
@@ -217,8 +212,8 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
   // it fixes the bug OCC5898. To develop more or less sensible criterion it is
   // necessary to deeply investigate this problem which is not possible in frames
   // of debugging.
-  Standard_Real aParamFirst = 3.0 * pc3d(1)   - 2.0 * pc3d(nbp);
-  Standard_Real aParamLast = 3.0 * pc3d(nbp) - 2.0 * pc3d(1);
+  Standard_Real aParamFirst = 3.0 * pc3d[0]   - 2.0 * pc3d[nbp - 1];
+  Standard_Real aParamLast = 3.0 * pc3d[nbp - 1] - 2.0 * pc3d[0];
 
   Standard_Real FirstPar = cons.FirstParameter();
   Standard_Real LastPar  = cons.LastParameter();
@@ -227,27 +222,29 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
   if (aParamLast > LastPar)
     aParamLast = LastPar;
 
-
   Standard_Real d2 = 0.0; // Maximum square deviation on the samples.
   const Standard_Real d = tol;
   const Standard_Integer nn = 2 * nbp;
-  const Standard_Real unsurnn = 1.0/nn;
+  const Standard_Real unsurnn = 1.0 / nn;
+  Standard_Real tprev = aParamFirst;
   for(Standard_Integer i = 0; i <= nn; i++)
   {
     // Compute corresponding parameter on 2d curve.
     // It should be inside of 3d curve parameter space.
     Standard_Real t = unsurnn*i;
-    Standard_Real tc3d = pc3d(1)*(1.-t) + pc3d(nbp)*t;
+    Standard_Real tc3d = pc3d[0]*(1.0 - t) + pc3d[nbp - 1] * t; // weight function.
     gp_Pnt Pc3d = c3d->Value(tc3d);
     Standard_Real tcons;
-    BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
-                   aDegree,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
-    if (tcons < aParamFirst ||
+    BSplCLib::Eval(tc3d, Standard_False, 0, extrap_mode[0],
+                   aDegree, FlatKnots, 1, (Standard_Real&)Poles(1), tcons);
+
+    if (tcons < tprev ||
         tcons > aParamLast)
     {
       tol = Precision::Infinite();
       return Standard_False;
     }
+    tprev = tcons;
     gp_Pnt Pcons = cons.Value(tcons);
     Standard_Real temp = Pc3d.SquareDistance(Pcons);
     if(temp > d2) d2 = temp;
@@ -275,7 +272,8 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)&   C3D,
                                            const Handle(Geom2d_Curve)& C2D,
                                            const Handle(Geom_Surface)& S,
                                            const Standard_Real         Tol)
-: mySameParameter(Standard_True),
+: myDeltaMin(Precision::PConfusion()),
+  mySameParameter(Standard_True),
   myDone(Standard_False)
 {
   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
@@ -292,7 +290,8 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D
                                            const Handle(Geom2d_Curve)&       C2D,
                                            const Handle(Adaptor3d_HSurface)& S,
                                            const Standard_Real               Tol)
-: mySameParameter(Standard_True),
+: myDeltaMin(Precision::PConfusion()),
+  mySameParameter(Standard_True),
   myDone(Standard_False)
 {
   myC3d = C3D;
@@ -309,7 +308,8 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D
                                            const Handle(Adaptor2d_HCurve2d)& C2D,
                                            const Handle(Adaptor3d_HSurface)& S,
                                            const Standard_Real               Tol)
-: mySameParameter(Standard_True),
+: myDeltaMin(Precision::PConfusion()),
+  mySameParameter(Standard_True),
   myDone(Standard_False)
 {
   myC3d = C3D;
@@ -324,560 +324,620 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D
 //=======================================================================
 void Approx_SameParameter::Build(const Standard_Real Tolerance)
 {
-  const Standard_Real anErrorMAX = 1.0e15;
-  const Standard_Integer aMaxArraySize = 1000;
-  const Standard_Integer NCONTROL = 22;
-
-  Standard_Integer ii ;
-  Adaptor3d_CurveOnSurface CurveOnSurface(myHCurve2d,mySurf);
-  Standard_Real fcons = CurveOnSurface.FirstParameter();
-  Standard_Real lcons = CurveOnSurface.LastParameter();
-  Standard_Real fc3d = myC3d->FirstParameter();
-  Standard_Real lc3d = myC3d->LastParameter();
-
-  //Control tangents at the extremities to know if the
-  //reparametring is possible and calculate the tangents 
-  //at the extremities of the function of change of variable.
+  // Algorithm:
+  // 1) Build initial distribution. Increase number of samples in case of C0 continuity.
+  // 2.1) Check same parameter state on samples.
+  // 2.2) Compute parameters in 2d space if not same parameter.
+  // 3) Loop over poles number and try to interpolate 2d curve.
+  // 4) If loop is failed build 2d curve forcibly or use original pcurve.
+  Standard_Real qpcons[myMaxArraySize], qnewpcons[myMaxArraySize],
+                 qpc3d[myMaxArraySize], qnewpc3d[myMaxArraySize];
+
+  // Create and fill data structure.
+  Approx_SameParameter_Data aData;
+  aData.myCOnS = Adaptor3d_CurveOnSurface(myHCurve2d,mySurf);
+  aData.myC2dPF = aData.myCOnS.FirstParameter();
+  aData.myC2dPL = aData.myCOnS.LastParameter();
+  aData.myC3dPF = myC3d->FirstParameter();
+  aData.myC3dPL = myC3d->LastParameter();
+  aData.myNbPnt = 0; // No points initially.
+  aData.myPC2d = qpcons;
+  aData.myPC3d = qpc3d;
+  aData.myNewPC2d = qnewpcons;
+  aData.myNewPC3d = qnewpc3d;
+  aData.myTol = Tolerance;
+
+  // Build initial distribution.
+  if (!BuildInitialDistribution(aData))
+  {
+    mySameParameter = Standard_False;
+    myDone = Standard_False;
+    return;
+  }
+
+  // Check same parameter state on distribution.
+  Standard_Real aMaxSqDeviation = 0.0;
+  const Standard_Real aPercentOfBadProj = 0.3;
+  Standard_Integer aNbPnt = aData.myNbPnt - RealToInt(aPercentOfBadProj * aData.myNbPnt);
+  mySameParameter = CheckSameParameter(aData, aMaxSqDeviation);
+  if(mySameParameter)
+  {
+    myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
+    myDone = Standard_True;
+    return;
+  }
+  else
+  {
+    // Control number of sample points after checking sameparameter
+    // If number of points is less then initial one, it means that there are 
+    // problems with projection
+    if(aData.myNbPnt < aNbPnt )
+    {
+      myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
+      myCurve2d = Geom2dAdaptor::MakeCurve(myHCurve2d->Curve2d());
+      myDone = Standard_False;
+      return;
+    }
+  }
+
+  // Control tangents at the extremities to know if the
+  // reparametring is possible and calculate the tangents
+  // at the extremities of the function of change of variable.
   Standard_Real tangent[2] = { 0.0, 0.0 };
-  gp_Pnt Pcons,Pc3d;
-  gp_Vec Vcons,Vc3d;
+  if(!ComputeTangents(aData.myCOnS, tangent[0], tangent[1]))
+  {
+    // Cannot compute tangents.
+    mySameParameter = Standard_False;
+    myDone = Standard_False;
+    myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
+    return;
+  }
 
-  const Standard_Real Tol = Tolerance;
-  const Standard_Real Tol2 = Tol * Tol;
-  Standard_Real deltamin = Precision::PConfusion();
+  // There is at least one point where same parameter is broken.
+  // Try to build B-spline approximation curve using interpolation with degree 3.
+  // The loop is organized over number of poles.
+  GeomAbs_Shape aContinuity = myHCurve2d->Continuity();
+  if(aContinuity > GeomAbs_C1) aContinuity = GeomAbs_C1;
 
-  Standard_Real besttol2 = Tol2;
+  Standard_Real besttol2 = aData.myTol * aData.myTol,
+                tolsov = Precision::Infinite();
+  Standard_Boolean interpolok = Standard_False,
+                   hasCountChanged = Standard_False;
+  do
+  {
+    // Interpolation data.
+    Standard_Integer num_knots = aData.myNbPnt + 7;
+    Standard_Integer num_poles = aData.myNbPnt + 3;
+    TColStd_Array1OfReal    Poles(1, num_poles);
+    TColStd_Array1OfReal    FlatKnots(1 ,num_knots);
 
-  // Check tangency on curve border.
-  Standard_Boolean extrok = 1;
-  CurveOnSurface.D1(fcons,Pcons,Vcons);
-  myC3d->D1(fc3d,Pc3d,Vc3d);
-  Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
-  Standard_Real dmax2 = dist2;
+    if (!Interpolate(aData, tangent[0], tangent[1], Poles, FlatKnots))
+    {
+      // Interpolation fails.
+      mySameParameter = Standard_False;
+      myDone = Standard_False;
+      myTolReached = ComputeTolReached(myC3d, aData.myCOnS, 2 * myNbSamples);
+      return;
+    }
+
+    Standard_Real algtol = sqrt(besttol2);
+    interpolok = Check (FlatKnots, Poles, aData.myNbPnt+1, aData.myPC3d,
+                        myC3d, aData.myCOnS, algtol, tolsov);
+    tolsov = algtol;
 
-  Standard_Real magVcons = Vcons.Magnitude();
-  if (magVcons > 1.e-12)
-    tangent[0] = Vc3d.Magnitude() / magVcons;
-  else extrok = 0;
+    // Try to build 2d curve and check it for validity.
+    if(interpolok)
+    {
+      Standard_Real besttol = sqrt(besttol2);
 
-  CurveOnSurface.D1(lcons,Pcons,Vcons);
-  myC3d->D1(lc3d,Pc3d,Vc3d);
-  dist2 = Pcons.SquareDistance(Pc3d);
+      Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
+      tol1d = new TColStd_HArray1OfReal(1,2);
+      tol1d->SetValue(1, mySurf->UResolution(besttol));
+      tol1d->SetValue(2, mySurf->VResolution(besttol));
 
-  dmax2 = Max(dmax2, dist2);
-  magVcons = Vcons.Magnitude();
-  if (magVcons > 1.e-12)
-    tangent[1] = Vc3d.Magnitude() / magVcons;
-  else extrok = 0;
-
-
-  //Take a multiple of the sample pof CheckShape,
-  //at least the control points will be correct. No comment!!!
-
-  Standard_Boolean interpolok = 0;
-  Standard_Real tolsov = 1.e200;
-  //Take parameters with  constant step on the curve on surface
-  //and on curve 3d.
-  Standard_Real deltacons = lcons - fcons;
-  deltacons /= (NCONTROL);
-  Standard_Real deltac3d = lc3d - fc3d;
-  deltac3d /= (NCONTROL);
-
-  Standard_Real wcons = fcons;
-  Standard_Real wc3d  = fc3d;
-
-  Standard_Real qpcons[aMaxArraySize], qnewpcons[aMaxArraySize], 
-                 qpc3d[aMaxArraySize], qnewpc3d[aMaxArraySize];
-  Standard_Real * pcons = qpcons; Standard_Real * newpcons = qnewpcons;
-  Standard_Real * pc3d = qpc3d; Standard_Real * newpc3d = qnewpc3d;
-
-  for ( ii = 0 ; ii < NCONTROL; ii++) {
-    pcons[ii] = wcons;
-    pc3d[ii]  = wc3d;
-    wcons += deltacons;
-    wc3d  += deltac3d;
-  }
-  pcons[NCONTROL] = lcons;
-  pc3d[NCONTROL]  = lc3d;
+      Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d);
+      Standard_Integer aMaxDeg = 11, aMaxSeg = 1000;
+      AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,aData.myC3dPF,aData.myC3dPL,
+                                                aContinuity,aMaxDeg,aMaxSeg,ev);
 
-  // Change number of points in case of C0 continuity.
-  Standard_Integer New_NCONTROL = NCONTROL;
-  GeomAbs_Shape Continuity = myHCurve2d->Continuity();
-  if(Continuity > GeomAbs_C1) Continuity = GeomAbs_C1;
-  if(Continuity < GeomAbs_C1)
-  {
-    Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
-    TColStd_Array1OfReal Param_de_decoupeC1 (1, NbInt);
-    myHCurve2d->Intervals(Param_de_decoupeC1, GeomAbs_C1);
-    TColStd_SequenceOfReal new_par;
-    Standard_Integer inter = 1;
-    ii =1;
-    new_par.Append(fcons);
-
-    while(inter <= NbInt && Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
-    while(NbInt > 0 && Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
-
-    while(inter <= NbInt || (ii < NCONTROL && inter <= Param_de_decoupeC1.Length()) ) {
-      if(Param_de_decoupeC1(inter) < pcons[ii]) {
-        new_par.Append(Param_de_decoupeC1(inter));
-        if((pcons[ii] - Param_de_decoupeC1(inter)) <= deltamin) {
-          ii++;
-          if(ii > NCONTROL) {ii = NCONTROL;}
+      if (anApproximator.IsDone() || anApproximator.HasResult())
+      {
+        Adaptor3d_CurveOnSurface ACS = aData.myCOnS;
+        GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator);
+        Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2);
+        Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
+        aData.myCOnS.Load(aHCurve2d);
+        myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
+
+        const Standard_Real aMult = 250.0; // To be tolerant with discrete tolerance.
+        if (myTolReached < aMult * besttol ) 
+        {
+          myCurve2d = aC2d;
+          myHCurve2d = aHCurve2d;
+          myDone = Standard_True;
+          break;
         }
-        inter++;
-      }
-      else {
-        if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
-          new_par.Append(pcons[ii]);
+        else if(aData.myNbPnt < myMaxArraySize - 1)
+        {
+          interpolok = Standard_False;
+          aData.myCOnS = ACS;
+        }
+        else
+        {
+          break;
         }
-        ii++;
       }
     }
 
-    new_par.Append(lcons);
-    New_NCONTROL = new_par.Length() - 1;
-    // Simple protection if New_NCONTROL > allocated elements in array but one
-    // aMaxArraySize - 1 index may be filled after projection.
-    if (New_NCONTROL > aMaxArraySize - 1) {
-      mySameParameter = Standard_False;
+    if (!interpolok)
+      hasCountChanged = IncreaseNbPoles(Poles, FlatKnots, aData, besttol2);
+
+  } while(!interpolok && hasCountChanged);
+
+  if (!myDone)
+  {
+    // Loop is finished unsuccessfully. Fix tolerance by maximal deviation,
+    // using data from the last loop iteration or initial data. Use data set with minimal deflection.
+
+    // Original 2d curve.
+    aData.myCOnS.Load(myHCurve2d);
+    myTolReached = ComputeTolReached(myC3d,aData.myCOnS, 2 * myNbSamples);
+    myCurve2d = Geom2dAdaptor::MakeCurve(myHCurve2d->Curve2d());
+
+    // Approximation curve.
+    Standard_Integer num_knots = aData.myNbPnt + 7;
+    Standard_Integer num_poles = aData.myNbPnt + 3;
+    TColStd_Array1OfReal    Poles(1, num_poles);
+    TColStd_Array1OfReal    FlatKnots(1 ,num_knots);
+
+    Interpolate(aData, tangent[0], tangent[1],
+                Poles, FlatKnots);
+
+    Standard_Real besttol = sqrt(besttol2);
+    Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
+    tol1d = new TColStd_HArray1OfReal(1,2) ;
+    tol1d->SetValue(1, mySurf->UResolution(besttol));
+    tol1d->SetValue(2, mySurf->VResolution(besttol));
+
+    Approx_SameParameter_Evaluator ev(FlatKnots, Poles, myHCurve2d);
+    AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,aData.myC3dPF,aData.myC3dPL,
+                                              aContinuity,11,40,ev);
+
+    if (!anApproximator.IsDone() &&
+        !anApproximator.HasResult() )
+    {
+      myDone = Standard_False;
       return;
     }
-    for(ii = 1; ii <= New_NCONTROL; ii++){
-      pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
+
+    GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator);
+    Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2);
+    Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
+    aData.myCOnS.Load(aHCurve2d);
+
+    Standard_Real anApproxTol = ComputeTolReached(myC3d,aData.myCOnS,2 * myNbSamples);
+    if (anApproxTol < myTolReached)
+    {
+      myTolReached = anApproxTol;
+      myCurve2d = aC2d;
+      myHCurve2d  = aHCurve2d;
     }
-    pc3d[New_NCONTROL]  = lc3d;
+    myDone = Standard_True;
   }
+}
 
-  // Check existing same parameter state.
-  Extrema_LocateExtPC Projector;
-  Projector.Initialize(myC3d->Curve(),fc3d,lc3d,Tol);
-
-  Standard_Integer count = 1;
-  Standard_Real previousp = fc3d, initp=0, curp;
-  Standard_Real bornesup = lc3d - deltamin;
-  Standard_Boolean projok = 0, 
-    use_parameter ;
-  for (ii = 1; ii < New_NCONTROL; ii++){    
-    CurveOnSurface.D0(pcons[ii],Pcons);
-    myC3d->D0(pc3d[ii],Pc3d);
-    dist2 = Pcons.SquareDistance(Pc3d);
-    use_parameter = (dist2 <= Tol2  && (pc3d[ii] > pc3d[count-1] + deltamin)) ;
-    Standard_Real aDistMin = RealLast();
-    if(use_parameter) {
+//=======================================================================
+//function : BuildInitialDistribution
+//purpose  : Sub-method in Build.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::BuildInitialDistribution(Approx_SameParameter_Data &theData) const
+{
+  // Take a multiple of the sample pof CheckShape,
+  // at least the control points will be correct.
+  // Take parameters with constant step on the curve on surface
+  // and on curve 3d.
+  const Standard_Real deltacons = (theData.myC2dPL - theData.myC2dPF) / myNbSamples;
+  const Standard_Real deltac3d  = (theData.myC3dPL - theData.myC3dPF) / myNbSamples;
+  Standard_Real wcons = theData.myC2dPF;
+  Standard_Real wc3d  = theData.myC3dPF;
+  for (Standard_Integer ii = 0 ; ii < myNbSamples; ii++)
+  {
+    theData.myPC2d[ii] = wcons;
+    theData.myPC3d[ii]  = wc3d;
+    wcons += deltacons;
+    wc3d  += deltac3d;
+  }
+  theData.myNbPnt = myNbSamples;
+  theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
+  theData.myPC3d[theData.myNbPnt]  = theData.myC3dPL;
 
-      if(dist2 > dmax2) dmax2 = dist2;
-      initp = previousp = pc3d[count] = pc3d[ii];
-      pcons[count] = pcons[ii];
-      count++;
-      
+  // Change number of points in case of C0 continuity.
+  GeomAbs_Shape Continuity = myHCurve2d->Continuity();
+  if(Continuity < GeomAbs_C1)
+  {
+    if (!IncreaseInitialNbSamples(theData))
+    {
+      // Number of samples is too big.
+      return Standard_False;
     }
-    else {
-      if(!projok) initp = pc3d[ii];
-      projok = mySameParameter = Standard_False;
-      Projector.Perform(Pcons, initp);
-      if (Projector.IsDone()) {
-        curp = Projector.Point().Parameter();
-        Standard_Real dist_2 = Projector.SquareDistance();
-        projok = Standard_True;
-        aDistMin = dist_2; 
-      }
-      else
-      {
-        ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
-        if(projok)
-        {
-          const gp_Pnt& ap1 =myC3d->Value(curp);
-          aDistMin = Pcons.SquareDistance(ap1);
-        }
-      }
-      projok = (projok && (curp > previousp + deltamin && curp < bornesup));
-      if(projok)
-      {
-        initp = previousp = pc3d[count] = curp;
-        pcons[count] = pcons[ii];
-        count++;
-       
-      }
-      else
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : IncreaseInitialNbSamples
+//purpose  : Get number of C1 intervals and build new distribution on them.
+//           Sub-method in BuildInitialDistribution.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::IncreaseInitialNbSamples(Approx_SameParameter_Data &theData) const
+{
+  Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
+  TColStd_Array1OfReal aC1Intervals (1, NbInt);
+  myHCurve2d->Intervals(aC1Intervals, GeomAbs_C1);
+
+  Standard_Integer inter = 1;
+  while(inter <= NbInt && aC1Intervals(inter) <= theData.myC3dPF + myDeltaMin) inter++;
+  while(NbInt > 0 && aC1Intervals(NbInt) >= theData.myC3dPL - myDeltaMin) NbInt--;
+
+  // Compute new parameters.
+  TColStd_SequenceOfReal aNewPar;
+  aNewPar.Append(theData.myC3dPF);
+  Standard_Integer ii = 1;
+  while(inter <= NbInt || (ii < myNbSamples && inter <= aC1Intervals.Length()) )
+  {
+    if(aC1Intervals(inter) < theData.myPC2d[ii])
+    {
+      aNewPar.Append(aC1Intervals(inter));
+      if((theData.myPC2d[ii] - aC1Intervals(inter)) <= myDeltaMin)
       {
-        Extrema_ExtPC PR(Pcons,myC3d->Curve(),fc3d,lc3d,Tol);
-        if(PR.IsDone())
+        ii++;
+        if(ii > myNbSamples)
         {
-          const Standard_Integer aNbExt = PR.NbExt();
-          if(aNbExt > 0)
-          {
-            Standard_Integer anIndMin = 0;
-            Standard_Real aCurDistMin = RealLast();
-            for(Standard_Integer i = 1; i <= aNbExt; i++)
-            {
-              const gp_Pnt &aP = PR.Point(i).Value();
-              Standard_Real aDist2 = aP.SquareDistance(Pcons);
-              if(aDist2 < aCurDistMin)
-              {
-                aCurDistMin = aDist2;
-                anIndMin = i;
-              }
-            }
-            if(anIndMin)
-            {
-              curp = PR.Point(anIndMin).Parameter();
-              if( curp > previousp + deltamin && curp < bornesup)
-              {
-                aDistMin = aCurDistMin;
-                initp = previousp = pc3d[count] = curp;
-                pcons[count] = pcons[ii];
-                count++;
-                projok = Standard_True;
-                
-              }
-            }
-         
-          }
+          ii = myNbSamples;
         }
       }
-      if(projok && besttol2 < aDistMin)
-        besttol2 = aDistMin;
-        
-      else if(!projok)
+      inter++;
+    }
+    else
+    {
+      if((aC1Intervals(inter) - theData.myPC2d[ii]) > myDeltaMin)
       {
-        //Projector
-#ifdef OCCT_DEBUG
-        std::cout << "Projection not done" << std::endl;
-#endif
+        aNewPar.Append(theData.myPC2d[ii]);
       }
+      ii++;
     }
   }
-
-  if(mySameParameter){
-    myTolReached = 1.5*sqrt(dmax2);
-    return;
+  // Simple protection if theNewNbPoints > allocated elements in array but one
+  // myMaxArraySize - 1 index may be filled after projection.
+  theData.myNbPnt = aNewPar.Length();
+  if (theData.myNbPnt > myMaxArraySize - 1)
+  {
+    return Standard_False;
   }
 
-  if(!extrok)
+  for(ii = 1; ii < theData.myNbPnt; ii++)
   {
-    // If not already SameP and tangent to mill, abandon.
-    mySameParameter = Standard_False;
-#ifdef OCCT_DEBUG
-    std::cout<<"SameParameter problem  : zero tangent to extremities"<<std::endl;
-#endif
-    return;
+    // Copy only internal points.
+    theData.myPC2d[ii] = theData.myPC3d[ii] = aNewPar.Value(ii + 1);
   }
+  theData.myPC3d[theData.myNbPnt]  = theData.myC3dPL;
+  theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
 
-  pcons[count] = lcons;
-  pc3d[count]  = lc3d;
+  return Standard_True;
+}
 
-  // There is at least one point where same parameter is broken.
-  // Try to build B-spline interpolation curve with degree 3.
-  // The loop is organized over number of poles.
-  Standard_Boolean hasCountChanged = Standard_False;
-  do
-  {
-    // The tables and their limits for the interpolation.
-    Standard_Integer num_knots = count + 7;
-    Standard_Integer num_poles = count + 3;
-    TColStd_Array1OfReal    Paramc3d(*pc3d,1,count+1);
-    TColStd_Array1OfReal    Paramcons(*pcons,1,count+1);
-    TColStd_Array1OfInteger ContactOrder(1,num_poles) ;
-    TColStd_Array1OfReal    Poles(1,num_poles) ;
-    TColStd_Array1OfReal    InterpolationParameters(1,num_poles) ;
-    TColStd_Array1OfReal    FlatKnots(1,num_knots) ; 
-
-    // Fill tables taking attention to end values.
-    ContactOrder.Init(0);
-    ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
-
-    FlatKnots(1) = FlatKnots(2) = FlatKnots(3) = FlatKnots(4) = fc3d;
-    FlatKnots(num_poles + 1) = FlatKnots(num_poles + 2) = 
-      FlatKnots(num_poles + 3) = FlatKnots(num_poles + 4) = lc3d;
-
-    Poles(1) = fcons; Poles(num_poles) = lcons;
-    Poles(2) = tangent[0]; Poles(num_poles - 1) = tangent[1];
-
-    InterpolationParameters(1) = InterpolationParameters(2) = fc3d;
-    InterpolationParameters(num_poles - 1) = InterpolationParameters(num_poles) = lc3d;
-
-    for (ii = 3; ii <= num_poles - 2; ii++) {
-      Poles(ii) = Paramcons(ii - 1);
-      InterpolationParameters(ii) = FlatKnots(ii+2) = Paramc3d(ii - 1);
-    }
-    Standard_Integer inversion_problem;
-    BSplCLib::Interpolate(3,FlatKnots,InterpolationParameters,ContactOrder,
-      1,Poles(1),inversion_problem);
-    if(inversion_problem) {
-      throw Standard_ConstructionError();
-    }
+//=======================================================================
+//function : CheckSameParameter
+//purpose  : Sub-method in Build.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::CheckSameParameter(Approx_SameParameter_Data &theData,
+                                                          Standard_Real &theSqDist) const
+{
+  const Standard_Real Tol2 = theData.myTol * theData.myTol;
+  Standard_Boolean isSameParam = Standard_True;
 
-    // Test if par2d(par3d) is monotonous function or not           ----- IFV, Jan 2000
-    // and try to insert new point to improve BSpline interpolation
+  // Compute initial distance on boundary points.
+  gp_Pnt Pcons, Pc3d;
+  theData.myCOnS.D0(theData.myC2dPF, Pcons);
+  myC3d->D0(theData.myC3dPF, Pc3d);
+  Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
+  Standard_Real dmax2 = dist2;
 
-    Standard_Integer extrap_mode[2] ;
-    extrap_mode[0] = extrap_mode[1] = 3;
-    Standard_Real eval_result[2] ;
-    Standard_Integer DerivativeRequest = 0;
-    Standard_Real *PolesArray =
-      (Standard_Real *) &Poles(Poles.Lower()) ;
+  theData.myCOnS.D0(theData.myC2dPL, Pcons);
+  myC3d->D0(theData.myC3dPL, Pc3d);
+  dist2 = Pcons.SquareDistance(Pc3d);
+  dmax2 = Max(dmax2, dist2);
 
-    Standard_Integer newcount = 0;
-    for (ii = 0; ii < count; ii++) {
+  Extrema_LocateExtPC Projector;
+  Projector.Initialize(myC3d->Curve(), theData.myC3dPF, theData.myC3dPL, theData.myTol);
 
-      newpcons[newcount] = pcons[ii];
-      newpc3d[newcount] = pc3d[ii];
-      newcount++;
+  Standard_Integer count = 1;
+  Standard_Real previousp = theData.myC3dPF, initp=0, curp;
+  Standard_Real bornesup = theData.myC3dPL - myDeltaMin;
+  Standard_Boolean isProjOk = Standard_False;
+  for (Standard_Integer ii = 1; ii < theData.myNbPnt; ii++)
+  {
+    theData.myCOnS.D0(theData.myPC2d[ii],Pcons);
+    myC3d->D0(theData.myPC3d[ii],Pc3d);
+    dist2 = Pcons.SquareDistance(Pc3d);
 
-      if(count - ii + newcount == aMaxArraySize) continue;
+    // Same parameter point.
+    Standard_Boolean isUseParam = (dist2 <= Tol2  && // Good distance.
+      (theData.myPC3d[ii] > theData.myPC3d[count-1] + myDeltaMin)); // Point is separated from previous.
+    if(isUseParam)
+    {
+      if(dmax2 < dist2)
+        dmax2 = dist2;
+      initp = previousp = theData.myPC3d[count] = theData.myPC3d[ii];
+      theData.myPC2d[count] = theData.myPC2d[ii];
+      count++;
+      continue;
+    }
 
-      BSplCLib::Eval(.5*(pc3d[ii]+pc3d[ii+1]), Standard_False, DerivativeRequest,
-        extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
+    // Local search: local extrema and iterative projection algorithm.
+    if(!isProjOk)
+      initp = theData.myPC3d[ii];
+    isProjOk = isSameParam = Standard_False;
+    Projector.Perform(Pcons, initp);
+    if (Projector.IsDone())
+    {
+      // Local extrema is found.
+      curp = Projector.Point().Parameter();
+      isProjOk = Standard_True;
+    }
+    else
+    {
+      ProjectPointOnCurve(initp,Pcons,theData.myTol,30,myC3d->Curve(),isProjOk,curp);
+    }
+    isProjOk = isProjOk && // Good projection.
+               curp > previousp + myDeltaMin && // Point is separated from previous.
+               curp < bornesup; // Inside of parameter space.
+    if(isProjOk)
+    {
+      initp = previousp = theData.myPC3d[count] = curp;
+      theData.myPC2d[count] = theData.myPC2d[ii];
+      count++;
+      continue;
+    }
 
-      if(eval_result[0] < pcons[ii] || eval_result[0] > pcons[ii+1]) {
-        Standard_Real ucons = 0.5*(pcons[ii]+pcons[ii+1]);
-        Standard_Real uc3d  = 0.5*(pc3d[ii]+pc3d[ii+1]);
+    // Whole parameter space search using general extrema.
+    Extrema_ExtPC PR(Pcons,myC3d->Curve(),theData.myC3dPF, theData.myC3dPL,theData.myTol);
+    if (!PR.IsDone() || PR.NbExt() == 0) // Lazy evaluation is used.
+      continue;
 
-        CurveOnSurface.D0(ucons,Pcons);
-        Projector.Perform(Pcons, uc3d);
-        if (Projector.IsDone()) {
-          curp = Projector.Point().Parameter();
-          Standard_Real dist_2 = Projector.SquareDistance();
-          if(dist_2 > besttol2) besttol2 = dist_2;
-          projok = 1;
-        }
-        else {
-          ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
-        }
-        if(projok){
-          if(curp > pc3d[ii] + deltamin && curp < pc3d[ii+1] - deltamin){
-            newpc3d[newcount] = curp;
-            newpcons[newcount] = ucons;
-            newcount ++;
-          }
-        }
-        else {
-#ifdef OCCT_DEBUG
-          std::cout << "Projection not done" << std::endl;
-#endif
-        }
+    const Standard_Integer aNbExt = PR.NbExt();
+    Standard_Integer anIndMin = 0;
+    Standard_Real aCurDistMin = RealLast();
+    for(Standard_Integer i = 1; i <= aNbExt; i++)
+    {
+      const gp_Pnt &aP = PR.Point(i).Value();
+      Standard_Real aDist2 = aP.SquareDistance(Pcons);
+      if(aDist2 < aCurDistMin)
+      {
+        aCurDistMin = aDist2;
+        anIndMin = i;
       }
-
     }
-
-    newpc3d[newcount] = pc3d[count];
-    newpcons[newcount] = pcons[count];
-    Standard_Real * temp;
-    temp = pc3d;
-    pc3d = newpc3d;
-    newpc3d = temp;
-    temp = pcons;
-    pcons = newpcons;
-    newpcons = temp;
-
-    if((count != newcount) && newcount < aMaxArraySize)
+    if(anIndMin)
     {
-      hasCountChanged = Standard_True;
-      count = newcount;
-      continue;
+      curp = PR.Point(anIndMin).Parameter();
+      if( curp > previousp + myDeltaMin && curp < bornesup)
+      {
+        initp = previousp = theData.myPC3d[count] = curp;
+        theData.myPC2d[count] = theData.myPC2d[ii];
+        count++;
+        isProjOk = Standard_True;
+      }
     }
+  }
+  theData.myNbPnt = count;
+  theData.myPC2d[theData.myNbPnt] = theData.myC2dPL;
+  theData.myPC3d[theData.myNbPnt] = theData.myC3dPL;
 
-    count = newcount;
+  theSqDist = dmax2;
+  return isSameParam;
+}
 
-    Standard_Real algtol = sqrt(besttol2);
+//=======================================================================
+//function : ComputeTangents
+//purpose  : Sub-method in Build.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::ComputeTangents(const Adaptor3d_CurveOnSurface & theCOnS,
+                                                       Standard_Real &theFirstTangent,
+                                                       Standard_Real &theLastTangent) const
+{
+  const Standard_Real aSmallMagnitude = 1.0e-12;
+  // Check tangency on curve border.
+  gp_Pnt aPnt, aPntCOnS;
+  gp_Vec aVec, aVecConS;
+
+  // First point.
+  const Standard_Real aParamFirst = myC3d->FirstParameter();
+  theCOnS.D1(aParamFirst, aPntCOnS, aVecConS);
+  myC3d->D1(aParamFirst, aPnt, aVec);
+  Standard_Real aMagnitude = aVecConS.Magnitude();
+  if (aMagnitude > aSmallMagnitude)
+    theFirstTangent = aVec.Magnitude() / aMagnitude;
+  else
+    return Standard_False;
+
+  // Last point.
+  const Standard_Real aParamLast = myC3d->LastParameter();
+  theCOnS.D1(aParamLast,aPntCOnS,aVecConS);
+  myC3d->D1(aParamLast, aPnt, aVec);
+
+  aMagnitude = aVecConS.Magnitude();
+  if (aMagnitude > aSmallMagnitude)
+    theLastTangent = aVec.Magnitude() / aMagnitude;
+  else
+    return Standard_False;
+
+  return Standard_True;
+}
 
-    interpolok = Check (FlatKnots, Poles, count+1, Paramc3d, Paramcons, 
-                        myC3d, CurveOnSurface, algtol, tolsov);
+//=======================================================================
+//function : Interpolate
+//purpose  : Sub-method in Build.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::Interpolate(const Approx_SameParameter_Data & theData,
+                                       const Standard_Real aTangFirst,
+                                       const Standard_Real aTangLast,
+                                       TColStd_Array1OfReal & thePoles,
+                                       TColStd_Array1OfReal & theFlatKnots) const
+{
+  Standard_Integer num_poles = theData.myNbPnt + 3;
+  TColStd_Array1OfInteger ContactOrder(1,num_poles);
+  TColStd_Array1OfReal    aParameters(1, num_poles);
 
-    if (Precision::IsInfinite(algtol)) {
-      mySameParameter = Standard_False;
-#ifdef OCCT_DEBUG
-      std::cout<<"SameParameter problem  : function of interpolation of parametration at mills !!"<<std::endl;
-#endif
-      return;
-    }
+  // Fill tables taking attention to end values.
+  ContactOrder.Init(0);
+  ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
 
-    tolsov = algtol;
+  theFlatKnots(1) = theFlatKnots(2) = theFlatKnots(3) = theFlatKnots(4) = theData.myC3dPF;
+  theFlatKnots(num_poles + 1) = theFlatKnots(num_poles + 2) = 
+    theFlatKnots(num_poles + 3) = theFlatKnots(num_poles + 4) = theData.myC3dPL;
 
-    interpolok = (interpolok || // Good result.
-                  count >= aMaxArraySize - 1 ); // Number of points.
+  thePoles(1) = theData.myC2dPF; thePoles(num_poles) = theData.myC2dPL;
+  thePoles(2) = aTangFirst; thePoles(num_poles - 1) = aTangLast;
 
-    if(interpolok) {
-      Standard_Real besttol = sqrt(besttol2);
+  aParameters(1) = aParameters(2) = theData.myC3dPF;
+  aParameters(num_poles - 1) = aParameters(num_poles) = theData.myC3dPL;
 
-      Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
-      tol1d = new TColStd_HArray1OfReal(1,2) ;
-      tol1d->SetValue(1, mySurf->UResolution(besttol));
-      tol1d->SetValue(2, mySurf->VResolution(besttol));
+  for (Standard_Integer ii = 3; ii <= num_poles - 2; ii++)
+  {
+    thePoles(ii) = theData.myPC2d[ii - 2];
+    aParameters(ii) = theFlatKnots(ii+2) = theData.myPC3d[ii - 2];
+  }
+  Standard_Integer inversion_problem;
+  BSplCLib::Interpolate(3,theFlatKnots,aParameters,ContactOrder,
+                        1,thePoles(1),inversion_problem);
+  if(inversion_problem)
+  {
+    return Standard_False;
+  }
 
-      Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d); 
-      AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
-        Continuity,11,40,ev);
+  return Standard_True;
+}
 
-      if (anApproximator.IsDone() || anApproximator.HasResult()) {
-        Adaptor3d_CurveOnSurface ACS = CurveOnSurface;
-        GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator) ;
-        Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
-        Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
-        CurveOnSurface.Load(aHCurve2d);
+//=======================================================================
+//function : IncreaseNbPoles
+//purpose  : Sub-method in Build.
+//=======================================================================
+Standard_Boolean Approx_SameParameter::IncreaseNbPoles(const TColStd_Array1OfReal & thePoles,
+                                                       const TColStd_Array1OfReal & theFlatKnots,
+                                                       Approx_SameParameter_Data & theData,
+                                                       Standard_Real &theBestSqTol) const
+{
+  Extrema_LocateExtPC Projector;
+  Projector.Initialize(myC3d->Curve(), myC3d->FirstParameter(), myC3d->LastParameter(), theData.myTol);
+  Standard_Real curp = 0.0;
+  Standard_Boolean projok = Standard_False;
 
-        myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
+  // Project middle point to fix parameterization and check projection existence.
+  const Standard_Integer aDegree = 3;
+  const Standard_Integer DerivativeRequest = 0;
+  Standard_Integer extrap_mode[2] = {aDegree, aDegree};
+  Standard_Real eval_result;
+  Standard_Real *PolesArray = (Standard_Real *) &thePoles(thePoles.Lower());
+  Standard_Integer newcount = 0;
+  for (Standard_Integer ii = 0; ii < theData.myNbPnt; ii++)
+  {
+    theData.myNewPC2d[newcount] = theData.myPC2d[ii];
+    theData.myNewPC3d[newcount] = theData.myPC3d[ii];
+    newcount++;
 
-        if(myTolReached > anErrorMAX)
-        {
-          //This tolerance is too big. Probably, we will not be able to get 
-          //edge with sameparameter in this case.
+    if(theData.myNbPnt - ii + newcount == myMaxArraySize) continue;
 
-          myDone = Standard_False;
-          return;
-        }
+    BSplCLib::Eval(0.5*(theData.myPC3d[ii]+theData.myPC3d[ii+1]), Standard_False, DerivativeRequest,
+      extrap_mode[0], 3, theFlatKnots, 1, PolesArray[0], eval_result);
 
-        if( (myTolReached < 250.0*besttol) || 
-            (count >= aMaxArraySize-2) || 
-            !hasCountChanged) //if count does not change after adding new point
-                              //(else we can have circularity)
-        {
-          myCurve2d = aC2d;
-          myHCurve2d  = aHCurve2d;
-          myDone = Standard_True;
-        }
-        else
-        {
-          interpolok = Standard_False;
-          CurveOnSurface = ACS;
-        }
-      } 
-    }
-    
-    if(!interpolok)
+    if(eval_result < theData.myPC2d[ii] || eval_result > theData.myPC2d[ii+1])
     {
+      Standard_Real ucons = 0.5*(theData.myPC2d[ii]+theData.myPC2d[ii+1]);
+      Standard_Real uc3d  = 0.5*(theData.myPC3d[ii]+theData.myPC3d[ii+1]);
 
-      newcount = 0;
-      for(Standard_Integer n = 0; n < count; n++){
-        newpc3d[newcount] = pc3d[n];
-        newpcons[newcount] = pcons[n];
-        newcount ++;
-
-        if(count - n + newcount == aMaxArraySize) continue;
-
-        Standard_Real ucons = 0.5*(pcons[n]+pcons[n+1]);
-        Standard_Real uc3d  = 0.5*(pc3d[n]+pc3d[n+1]);
-
-        CurveOnSurface.D0(ucons,Pcons);
-        Projector.Perform(Pcons, uc3d);
-        if (Projector.IsDone()) {
-          curp = Projector.Point().Parameter();
-          Standard_Real dist_2 = Projector.SquareDistance();
-          if(dist_2 > besttol2) besttol2 = dist_2;
-          projok = 1;
-        }
-        else {
-          ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
-        }
-        if(projok){
-          if(curp > pc3d[n] + deltamin && curp < pc3d[n+1] - deltamin){
-            newpc3d[newcount] = curp;
-            newpcons[newcount] = ucons;
-            newcount ++;
-          }
-        }
-        else {
-#ifdef OCCT_DEBUG
-          std::cout << "Projection not done" << std::endl;
-#endif
-        }
-      }
-      newpc3d[newcount] = pc3d[count];
-      newpcons[newcount] = pcons[count];
-      Standard_Real * tempx;
-      tempx = pc3d;
-      pc3d = newpc3d;
-      newpc3d = tempx;
-      tempx = pcons;
-      pcons = newpcons;
-      newpcons = tempx;
-      
-      if(count != newcount)
+      gp_Pnt Pcons;
+      theData.myCOnS.D0(ucons,Pcons);
+      Projector.Perform(Pcons, uc3d);
+      if (Projector.IsDone())
       {
-        count = newcount;
-        hasCountChanged = Standard_True;
+        curp = Projector.Point().Parameter();
+        Standard_Real dist_2 = Projector.SquareDistance();
+        if(dist_2 > theBestSqTol) theBestSqTol = dist_2;
+        projok = 1;
       }
       else
       {
-        hasCountChanged = Standard_False;
+        ProjectPointOnCurve(uc3d,Pcons,theData.myTol,30,myC3d->Curve(),projok,curp);
+      }
+      if(projok)
+      {
+        if(curp > theData.myPC3d[ii] + myDeltaMin && curp < theData.myPC3d[ii+1] - myDeltaMin)
+        {
+          theData.myNewPC3d[newcount] = curp;
+          theData.myNewPC2d[newcount] = ucons;
+          newcount ++;
+        }
       }
     }
-  } while(!interpolok && hasCountChanged);
+  } //     for (ii = 0; ii < count; ii++)
+  theData.myNewPC3d[newcount] = theData.myPC3d[theData.myNbPnt];
+  theData.myNewPC2d[newcount] = theData.myPC2d[theData.myNbPnt];
 
-  if (!myDone)
+  if((theData.myNbPnt != newcount) && newcount < myMaxArraySize - 1)
   {
-    // Loop is finished unsuccessfully. Fix tolerance by maximal deviation,
-    // using data from the last loop iteration.
-    Standard_Integer num_knots = count + 7;
-    Standard_Integer num_poles = count + 3;
-    TColStd_Array1OfReal    Paramc3d(*pc3d,1,count + 1);
-    TColStd_Array1OfReal    Paramcons(*pcons,1,count + 1);
-    TColStd_Array1OfInteger ContactOrder(1,num_poles) ;
-    TColStd_Array1OfReal    Poles(1,num_poles) ;
-    TColStd_Array1OfReal    InterpolationParameters(1,num_poles) ;
-    TColStd_Array1OfReal    FlatKnots(1,num_knots) ; 
-
-    // Fill tables taking attention to end values.
-    ContactOrder.Init(0);
-    ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
-
-    FlatKnots(1) = FlatKnots(2) = FlatKnots(3) = FlatKnots(4) = fc3d;
-    FlatKnots(num_poles + 1) = FlatKnots(num_poles + 2) = 
-      FlatKnots(num_poles + 3) = FlatKnots(num_poles + 4) = lc3d;
-
-    Poles(1) = fcons; Poles(num_poles) = lcons;
-    Poles(2) = tangent[0]; Poles(num_poles - 1) = tangent[1];
-
-    InterpolationParameters(1) = InterpolationParameters(2) = fc3d;
-    InterpolationParameters(num_poles - 1) = InterpolationParameters(num_poles) = lc3d;
-
-    for (ii = 3; ii <= num_poles - 2; ii++)
-    {
-      Poles(ii) = Paramcons(ii - 1);
-      InterpolationParameters(ii) = FlatKnots(ii+2) = Paramc3d(ii - 1);
-    }
-    Standard_Integer inversion_problem;
-    BSplCLib::Interpolate(3,FlatKnots,InterpolationParameters,ContactOrder,
-                          1,Poles(1),inversion_problem);
-    if(inversion_problem)
-    {
-      throw Standard_ConstructionError();
-    }
+    // Distribution is changed.
+    theData.Swap(newcount);
+    return Standard_True;
+  }
 
-    Standard_Real besttol = sqrt(besttol2);
-    Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
-    tol1d = new TColStd_HArray1OfReal(1,2) ;
-    tol1d->SetValue(1, mySurf->UResolution(besttol));
-    tol1d->SetValue(2, mySurf->VResolution(besttol));
+  // Increase number of samples in two times.
+  newcount = 0;
+  for(Standard_Integer n = 0; n < theData.myNbPnt; n++)
+  {
+    theData.myNewPC3d[newcount] = theData.myPC3d[n];
+    theData.myNewPC2d[newcount] = theData.myPC2d[n];
+    newcount ++;
 
-    Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d); 
-    AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
-                                              Continuity,11,40,ev);
+    if(theData.myNbPnt - n + newcount == myMaxArraySize) continue;
 
-    if (!anApproximator.IsDone() &&
-        !anApproximator.HasResult() )
+    Standard_Real ucons = 0.5*(theData.myPC2d[n]+theData.myPC2d[n+1]);
+    Standard_Real uc3d  = 0.5*(theData.myPC3d[n]+theData.myPC3d[n+1]);
+
+    gp_Pnt Pcons;
+    theData.myCOnS.D0(ucons,Pcons);
+    Projector.Perform(Pcons, uc3d);
+    if (Projector.IsDone()) 
     {
-      myDone = Standard_False;
-      return;
+      curp = Projector.Point().Parameter();
+      Standard_Real dist_2 = Projector.SquareDistance();
+      if(dist_2 > theBestSqTol) theBestSqTol = dist_2;
+      projok = 1;
     }
-
-    GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator) ;
-    Handle(Geom2d_BSplineCurve) aC2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
-    Handle(Adaptor2d_HCurve2d) aHCurve2d = new Geom2dAdaptor_HCurve(aC2d);
-    CurveOnSurface.Load(aHCurve2d);
-
-    myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
-
-    if(myTolReached > anErrorMAX)
+    else 
     {
-      //This tolerance is too big. Probably, we will not be able to get
-      //edge with sameparameter in this case.
-      myDone = Standard_False;
-      return;
+      ProjectPointOnCurve(uc3d,Pcons,theData.myTol,30,myC3d->Curve(),projok,curp);
     }
+    if(projok)
+    {
+      if(curp > theData.myPC3d[n] + myDeltaMin && curp < theData.myPC3d[n+1] - myDeltaMin)
+      {
+        theData.myNewPC3d[newcount] = curp;
+        theData.myNewPC2d[newcount] = ucons;
+        newcount ++;
+      }
+    }
+  }
+  theData.myNewPC3d[newcount] = theData.myPC3d[theData.myNbPnt];
+  theData.myNewPC2d[newcount] = theData.myPC2d[theData.myNbPnt];
 
-    myCurve2d = aC2d;
-    myHCurve2d  = aHCurve2d;
-    myDone = Standard_True;
+  if(theData.myNbPnt != newcount)
+  {
+    // Distribution is changed.
+    theData.Swap(newcount);
+    return Standard_True;
   }
+
+  return Standard_False;
 }