0024624: Lost word in license statement in source files
[occt.git] / src / Approx / Approx_SameParameter.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 6c9e58a..3fc41f3
@@ -1,23 +1,18 @@
 // Created on: 1995-06-06
 // Created by: Xavier BENVENISTE
 // Copyright (c) 1995-1999 Matra Datavision
-// Copyright (c) 1999-2012 OPEN CASCADE SAS
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
 //
-// The content of this file is subject to the Open CASCADE Technology Public
-// License Version 6.5 (the "License"). You may not use the content of this file
-// except in compliance with the License. Please obtain a copy of the License
-// at http://www.opencascade.org and read it completely before using this file.
+// This file is part of Open CASCADE Technology software library.
 //
-// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
-// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
 //
-// The Original Code and all software distributed under the License is
-// distributed on an "AS IS" basis, without warranty of any kind, and the
-// Initial Developer hereby disclaims all such warranties, including without
-// limitation, any warranties of merchantability, fitness for a particular
-// purpose or non-infringement. Please see the License for the specific terms
-// and conditions governing the rights and limitations under the License.
-
+// 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 <GeomAdaptor_HCurve.hxx>
 #include <GeomAdaptor_Surface.hxx>
 #include <GeomAdaptor_HSurface.hxx>
-//#include <GCPnts_UniformDeflection.hxx>
 #include <GCPnts_QuasiUniformDeflection.hxx>
 #include <Extrema_LocateExtPC.hxx>
 #include <AdvApprox_ApproxAFunction.hxx>
 #include <GeomLib_MakeCurvefromApprox.hxx>
 #include <Precision.hxx>
-
-#define MAX_ARRAY_SIZE 1000 // IFV, Jan 2000
+#include <Extrema_ExtPC.hxx>
 
 #ifdef DEB
 #ifdef DRAW
@@ -55,32 +48,32 @@ static Standard_Integer NbCurve = 0;
 
 
 static void ProjectPointOnCurve(const Standard_Real      InitValue,
-                               const gp_Pnt             APoint,
-                               const Standard_Real      Tolerance,
-                               const Standard_Integer   NumIteration,
-                               const Adaptor3d_Curve&     Curve,
-                               Standard_Boolean&        Status,
-                               Standard_Real&           Result)
+  const gp_Pnt             APoint,
+  const Standard_Real      Tolerance,
+  const Standard_Integer   NumIteration,
+  const Adaptor3d_Curve&     Curve,
+  Standard_Boolean&        Status,
+  Standard_Real&           Result)
 {
   Standard_Integer num_iter = 0,
-  not_done = 1,
-  ii ;
-  
+    not_done = 1,
+    ii ;
+
   gp_Pnt a_point ;
   gp_Vec   vector,
-  d1,
-  d2 ;
+    d1,
+    d2 ;
   Standard_Real func,
-  func_derivative,
-  param = InitValue ;
+    func_derivative,
+    param = InitValue ;
   Status = Standard_False ;
   Standard_Real Toler = 1.0e-12;
   do {
     num_iter += 1 ;
     Curve.D2(param,
-             a_point,
-            d1,
-            d2) ;
+      a_point,
+      d1,
+      d2) ;
     for (ii = 1 ; ii <= 3 ; ii++) {
       vector.SetCoord(ii, APoint.Coord(ii) - a_point.Coord(ii)) ;
     }
@@ -92,19 +85,19 @@ static void ProjectPointOnCurve(const Standard_Real      InitValue,
       Status = Standard_True ;
     }
     else
-      { // fixing a bug PRO18577 : avoid divizion by zero
-       if( Abs(func_derivative) > Toler )  {
-         param -= func / func_derivative ;
-       }
-       param = Max(param,Curve.FirstParameter()) ;
-       param = Min(param,Curve.LastParameter())  ;
-       Status = Standard_True ;
+    { // fixing a bug PRO18577 : avoid divizion by zero
+      if( Abs(func_derivative) > Toler )  {
+        param -= func / func_derivative ;
       }
+      param = Max(param,Curve.FirstParameter()) ;
+      param = Min(param,Curve.LastParameter())  ;
+      //Status = Standard_True ;
+    }
   } 
   while (not_done && num_iter <= NumIteration) ;
   Result = param ;
 }  
-     
+
 
 
 //=======================================================================
@@ -114,31 +107,31 @@ static void ProjectPointOnCurve(const Standard_Real      InitValue,
 
 class Approx_SameParameter_Evaluator : public AdvApprox_EvaluatorFunction
 {
- public:
+public:
   Approx_SameParameter_Evaluator (const TColStd_Array1OfReal& theFlatKnots, 
-                                  const TColStd_Array1OfReal& thePoles,
-                                  const Handle(Adaptor2d_HCurve2d)& theHCurve2d) 
+    const TColStd_Array1OfReal& thePoles,
+    const Handle(Adaptor2d_HCurve2d)& theHCurve2d) 
     : FlatKnots(theFlatKnots), Poles(thePoles), HCurve2d(theHCurve2d) {}
 
   virtual void Evaluate (Standard_Integer *Dimension,
-                        Standard_Real     StartEnd[2],
-                         Standard_Real    *Parameter,
-                         Standard_Integer *DerivativeRequest,
-                         Standard_Real    *Result, // [Dimension]
-                         Standard_Integer *ErrorCode);
-  
- private:
+    Standard_Real     StartEnd[2],
+    Standard_Real    *Parameter,
+    Standard_Integer *DerivativeRequest,
+    Standard_Real    *Result, // [Dimension]
+    Standard_Integer *ErrorCode);
+
+private:
   const TColStd_Array1OfReal& FlatKnots;
   const TColStd_Array1OfReal& Poles;
   Handle(Adaptor2d_HCurve2d) HCurve2d;
 };
 
 void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
-                                               Standard_Real    /*StartEnd*/[2],
-                                               Standard_Real    *Parameter,
-                                               Standard_Integer *DerivativeRequest,
-                                               Standard_Real    *Result,
-                                               Standard_Integer *ReturnCode) 
+  Standard_Real    /*StartEnd*/[2],
+  Standard_Real    *Parameter,
+  Standard_Integer *DerivativeRequest,
+  Standard_Real    *Result,
+  Standard_Integer *ReturnCode) 
 { 
   gp_Pnt2d Point ;
   gp_Vec2d Vector ;
@@ -151,16 +144,16 @@ void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
   // evaluate the 1D bspline that represents the change in parameterization
   //
   BSplCLib::Eval(*Parameter,
-                Standard_False,
-                *DerivativeRequest,
-                extrap_mode[0],
-                3,
-                FlatKnots,
-                1,
-                PolesArray[0],
-                 eval_result[0]) ;
-  
-  
+    Standard_False,
+    *DerivativeRequest,
+    extrap_mode[0],
+    3,
+    FlatKnots,
+    1,
+    PolesArray[0],
+    eval_result[0]) ;
+
+
   if (*DerivativeRequest == 0){
     HCurve2d->D0(eval_result[0],Point);
     Point.Coord(Result[0],Result[1]);
@@ -174,24 +167,22 @@ void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
 }
 
 static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
-                                      const Adaptor3d_CurveOnSurface& cons,
-                                      const Standard_Integer        nbp)
+  const Adaptor3d_CurveOnSurface& cons,
+  const Standard_Integer        nbp)
 {
   Standard_Real d2 = 0.;
-  Standard_Integer nn = nbp;
-  Standard_Real unsurnn = 1./nn;
-  Standard_Real first = c3d->FirstParameter();
-  Standard_Real last  = c3d->LastParameter();
-  for(Standard_Integer i = 0; i <= nn; i++){
-    Standard_Real t = unsurnn*i;
+  const Standard_Real first = c3d->FirstParameter();
+  const Standard_Real last  = c3d->LastParameter();
+  for(Standard_Integer i = 0; i <= nbp; i++){
+    Standard_Real t = IntToReal(i)/IntToReal(nbp);
     Standard_Real u = first*(1.-t) + last*t;
     gp_Pnt Pc3d = c3d->Value(u);
     gp_Pnt Pcons = cons.Value(u);
     if (Precision::IsInfinite(Pcons.X()) ||
-       Precision::IsInfinite(Pcons.Y()) ||
-       Precision::IsInfinite(Pcons.Z())) {
-      d2=Precision::Infinite();
-      break;
+      Precision::IsInfinite(Pcons.Y()) ||
+      Precision::IsInfinite(Pcons.Z())) {
+        d2=Precision::Infinite();
+        break;
     }
     Standard_Real temp = Pc3d.SquareDistance(Pcons);
     if(temp > d2) d2 = temp;
@@ -202,15 +193,15 @@ static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
 }
 
 static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots, 
-                              const TColStd_Array1OfReal& Poles,
-                              const Standard_Integer nbp,
-                             const TColStd_Array1OfReal& pc3d,
-//                           const TColStd_Array1OfReal& pcons,
-                             const TColStd_Array1OfReal& ,
-                             const Handle(Adaptor3d_HCurve)& c3d,
-                             const Adaptor3d_CurveOnSurface& cons,
-                             Standard_Real& tol,
-                             const Standard_Real oldtol)
+  const TColStd_Array1OfReal& Poles,
+  const Standard_Integer nbp,
+  const TColStd_Array1OfReal& pc3d,
+  //                         const TColStd_Array1OfReal& pcons,
+  const TColStd_Array1OfReal& ,
+  const Handle(Adaptor3d_HCurve)& c3d,
+  const Adaptor3d_CurveOnSurface& cons,
+  Standard_Real& tol,
+  const Standard_Real oldtol)
 {
   Standard_Real d = tol;
   Standard_Integer extrap_mode[2] ;
@@ -231,7 +222,7 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
     gp_Pnt Pc3d = c3d->Value(tc3d);
     Standard_Real tcons;
     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
-                  3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
+      3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
     gp_Pnt Pcons = cons.Value(tcons);
     Standard_Real temp = Pc3d.SquareDistance(Pcons);
     if(temp >= dglis) dglis = temp;
@@ -251,7 +242,7 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
     gp_Pnt Pc3d = c3d->Value(tc3d);
     Standard_Real tcons;
     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
-                  3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
+      3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
     gp_Pnt Pcons = cons.Value(tcons);
     Standard_Real temp = Pc3d.SquareDistance(Pcons);
     if(temp >= dglis) dglis = temp;
@@ -266,25 +257,25 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
   Standard_Real d2 = 0.;
   Standard_Integer nn = 2*nbp;
   Standard_Real unsurnn = 1./nn;
-//  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898 Begin
-// Correction of the interval of valid values. This condition has no sensible
-// grounds. But it is better then the old one (which is commented out) because
-// 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 firstborne= 2*pc3d(1)-pc3d(nbp);
-//   Standard_Real lastborne= 2*pc3d(nbp)-pc3d(1);
+  //  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898 Begin
+  // Correction of the interval of valid values. This condition has no sensible
+  // grounds. But it is better then the old one (which is commented out) because
+  // 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 firstborne= 2*pc3d(1)-pc3d(nbp);
+  //   Standard_Real lastborne= 2*pc3d(nbp)-pc3d(1);
   Standard_Real firstborne= 3.*pc3d(1)   - 2.*pc3d(nbp);
   Standard_Real lastborne = 3.*pc3d(nbp) - 2.*pc3d(1);
-//  Modified by skv - Wed Jun  2 11:50:03 2004 OCC5898 End
+  //  Modified by skv - Wed Jun  2 11:50:03 2004 OCC5898 End
   for(i = 0; i <= nn; i++){
     Standard_Real t = unsurnn*i;
     Standard_Real tc3d = pc3d(1)*(1.-t) + pc3d(nbp)*t;
     gp_Pnt Pc3d = c3d->Value(tc3d);
     Standard_Real tcons;
     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
-                  3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
+      3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
     if (tcons < firstborne || tcons > lastborne) {
       tol=Precision::Infinite();
       return Standard_False;
@@ -308,10 +299,10 @@ static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
 //=======================================================================
 
 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), myDone(Standard_False)
+  const Handle(Geom2d_Curve)& C2D,
+  const Handle(Geom_Surface)& S,
+  const Standard_Real         Tol):
+mySameParameter(Standard_True), myDone(Standard_False)
 {
   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
   myC3d      = new GeomAdaptor_HCurve(C3D);
@@ -326,10 +317,10 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)&   C3D,
 //=======================================================================
 
 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), myDone(Standard_False)
+  const Handle(Geom2d_Curve)&     C2D,
+  const Handle(Adaptor3d_HSurface)& S,
+  const Standard_Real            Tol):
+mySameParameter(Standard_True), myDone(Standard_False)
 {
   myC3d = C3D;
   mySurf = S;
@@ -344,10 +335,10 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D
 //=======================================================================
 
 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), myDone(Standard_False)
+  const Handle(Adaptor2d_HCurve2d)& C2D,
+  const Handle(Adaptor3d_HSurface)& S,
+  const Standard_Real            Tol):
+mySameParameter(Standard_True), myDone(Standard_False)
 {
   myC3d = C3D;
   mySurf = S;
@@ -360,9 +351,12 @@ Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D
 //function : Build
 //purpose  : 
 //=======================================================================
-
 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();
@@ -377,7 +371,7 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
   //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];
+  Standard_Real tangent[2] = { 0.0, 0.0 };
   gp_Pnt Pcons,Pc3d;
   gp_Vec Vcons,Vc3d;
 
@@ -417,11 +411,10 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
   //Take a multiple of the sample pof CheckShape,
   //at least the control points will be correct. No comment!!!
 
-  Standard_Integer NCONTROL = 22;
 #ifdef DEB
   Standard_Integer nbcoups = 0;
 #endif
-  
+
   Standard_Boolean interpolok = 0;
   Standard_Real tolsov = 1.e200;
   //Take parameters with  constant step on the curve on surface
@@ -433,9 +426,9 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 
   Standard_Real wcons = fcons;
   Standard_Real wc3d  = fc3d;
-  
-  Standard_Real qpcons[MAX_ARRAY_SIZE], qnewpcons[MAX_ARRAY_SIZE], 
-                qpc3d[MAX_ARRAY_SIZE], qnewpc3d[MAX_ARRAY_SIZE];
+
+  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;
 
@@ -450,51 +443,51 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 
   Standard_Integer New_NCONTROL = NCONTROL;
   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(Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
-     while(Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
-
-     while(inter <= NbInt || ii < NCONTROL) {
-       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;}
-        }
-        inter++;
-       }
-       else {
-        if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
-          new_par.Append(pcons[ii]);
-        }
-        ii++;
-       }
-     }
-     
-     new_par.Append(lcons);
-     New_NCONTROL = new_par.Length() - 1;
-     //simple protection if New_NCONTROL > allocated elements in array
-     if (New_NCONTROL > MAX_ARRAY_SIZE) {
-       mySameParameter = Standard_False;
-       return;
-     }
-     for(ii = 1; ii <= New_NCONTROL; ii++){
-       pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
-     }
-     pc3d[New_NCONTROL]  = lc3d;
-   }
-
-  
+    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(Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
+    while(Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
+
+    while(inter <= NbInt || ii < NCONTROL) {
+      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;}
+        }
+        inter++;
+      }
+      else {
+        if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
+          new_par.Append(pcons[ii]);
+        }
+        ii++;
+      }
+    }
+
+    new_par.Append(lcons);
+    New_NCONTROL = new_par.Length() - 1;
+    //simple protection if New_NCONTROL > allocated elements in array
+    if (New_NCONTROL > aMaxArraySize) {
+      mySameParameter = Standard_False;
+      return;
+    }
+    for(ii = 1; ii <= New_NCONTROL; ii++){
+      pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
+    }
+    pc3d[New_NCONTROL]  = lc3d;
+  }
+
+
   Extrema_LocateExtPC Projector;
   Projector.Initialize(myC3d->Curve(),fc3d,lc3d,Tol);
-  
+
   Standard_Integer count = 1;
   Standard_Real previousp = fc3d, initp=0, curp;//, deltamin = 50*Tolp;
   Standard_Real bornesup = lc3d - deltamin;
@@ -506,7 +499,7 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
     dist2 = Pcons.SquareDistance(Pc3d);
     use_parameter = (dist2 <= Tol2  && (pc3d[ii] > pc3d[count-1] + deltamin)) ;
     if(use_parameter) {
-      
+
       if(dist2 > dmax2) dmax2 = dist2;
       initp = previousp = pc3d[count] = pc3d[ii];
       pcons[count] = pcons[ii];
@@ -517,34 +510,72 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
       projok = mySameParameter = Standard_False;
       Projector.Perform(Pcons, initp);
       if (Projector.IsDone()) {
-       curp = Projector.Point().Parameter();
-       Standard_Real dist_2 = Projector.SquareDistance();
-       if(dist_2 > besttol2) besttol2 = dist_2;
-       projok = 1;
+        curp = Projector.Point().Parameter();
+        Standard_Real dist_2 = Projector.SquareDistance();
+        if(dist_2 > besttol2) besttol2 = dist_2;
+        projok = 1;
       }
-      else {
-       ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
+      else
+      {
+        ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
       }
-      if(projok){
-       if(curp > previousp + deltamin && curp < bornesup){
-         initp = previousp = pc3d[count] = curp;
-         pcons[count] = pcons[ii];
-         count++;
-       }
+      
+      if(projok)
+      {
+        if(curp > previousp + deltamin && curp < bornesup){
+          initp = previousp = pc3d[count] = curp;
+          pcons[count] = pcons[ii];
+          count++;
+        }
       }
-      else {
-#ifdef DEB 
-       // JAG
-       cout << "Projection not done" << endl;
+      else
+      {
+        Extrema_ExtPC PR(Pcons,myC3d->Curve(),fc3d,lc3d,Tol);
+        if(PR.IsDone())
+        {
+          const Standard_Integer aNbExt = PR.NbExt();
+          if(aNbExt > 0)
+          {
+            Standard_Integer anIndMin = 0;
+            Standard_Real aDistMin = 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 < aDistMin)
+              {
+                aDistMin = aDist2;
+                anIndMin = i;
+              }
+            }
+            curp = PR.Point(anIndMin).Parameter();
+            if(curp > previousp + deltamin && curp < bornesup)
+            {
+              initp = previousp = pc3d[count] = curp;
+              pcons[count] = pcons[ii];
+              count++;
+              projok = Standard_True;
+            }
+          }
+        }
+      }
+
+      if(!projok)
+      {
+        //Projector
+#ifdef DEB
+        // JAG
+        cout << "Projection not done" << endl;
 #endif
       }
     }
   }
+
   if(mySameParameter){
     myTolReached = 1.5*sqrt(dmax2);
     return;
   }
+
   if(!extrok) { // If not already SameP and tangent to mill, abandon.
     mySameParameter = Standard_False;
 #ifdef DEB
@@ -576,8 +607,11 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 #endif
   }
 #endif
-    
-  while(!interpolok){
+
+  Standard_Boolean hasCountChanged = Standard_False;
+
+  while(!interpolok)
+  {
     // The tables and their limits for the interpolation.
     Standard_Integer num_knots = count + 7;
     Standard_Integer num_poles = count + 3;
@@ -587,79 +621,79 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
     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);
+      1,Poles(1),inversion_problem);
     if(inversion_problem) {
       Standard_ConstructionError::Raise();
     }
 
     //-------------------------------------------
 #ifdef DEB
-  if (AffichFw) {
-    nbcoups ++;
-    char Name[17];
-    Name[0] = '\0';
-    Standard_Integer nnn = 100;
-    TColgp_Array1OfPnt2d    DEBP2d  (0,nnn);
-    TColStd_Array1OfInteger DEBMults(0,nnn); 
-    DEBMults.Init(1); DEBMults(0) = 2; DEBMults(nnn) = 2;
-    TColStd_Array1OfReal    DEBKnots(0,nnn);
-    Standard_Real du = (lc3d - fc3d) / nnn;
-    Standard_Real u3d = fc3d;
-    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()) ;
-
-    for (Standard_Integer DEBi = 0; DEBi <= nnn; DEBi++) {
-      DEBKnots(DEBi) = DEBi;
-      BSplCLib::Eval(u3d,
-                    Standard_False,
-                    DerivativeRequest,
-                    extrap_mode[0],
-                    3,
-                    FlatKnots,
-                    1,
-                    PolesArray[0],
-                    eval_result[0]) ;
-
-      DEBP2d  (DEBi) = gp_Pnt2d(u3d,eval_result[0]);
-      u3d += du;
-    }
+    if (AffichFw) {
+      nbcoups ++;
+      char Name[17];
+      Name[0] = '\0';
+      Standard_Integer nnn = 100;
+      TColgp_Array1OfPnt2d    DEBP2d  (0,nnn);
+      TColStd_Array1OfInteger DEBMults(0,nnn); 
+      DEBMults.Init(1); DEBMults(0) = 2; DEBMults(nnn) = 2;
+      TColStd_Array1OfReal    DEBKnots(0,nnn);
+      Standard_Real du = (lc3d - fc3d) / nnn;
+      Standard_Real u3d = fc3d;
+      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()) ;
+
+      for (Standard_Integer DEBi = 0; DEBi <= nnn; DEBi++) {
+        DEBKnots(DEBi) = DEBi;
+        BSplCLib::Eval(u3d,
+          Standard_False,
+          DerivativeRequest,
+          extrap_mode[0],
+          3,
+          FlatKnots,
+          1,
+          PolesArray[0],
+          eval_result[0]) ;
+
+        DEBP2d  (DEBi) = gp_Pnt2d(u3d,eval_result[0]);
+        u3d += du;
+      }
 
-    Handle(Geom2d_BSplineCurve) DEBBS = 
-      new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
-    Sprintf(Name,"DEBC2d_%d_%d",NbCurve,nbcoups );
+      Handle(Geom2d_BSplineCurve) DEBBS = 
+        new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
+      Sprintf(Name,"DEBC2d_%d_%d",NbCurve,nbcoups );
 #ifdef DRAW
-    DrawTrSurf::Set(Name,DEBBS);
+      DrawTrSurf::Set(Name,DEBBS);
 #endif
-  }
+    }
 #endif
-//-------------------------------------------    
+    //-------------------------------------------    
 
-//-------------------------------------------
+    //-------------------------------------------
     // Test if par2d(par3d) is monotonous function or not           ----- IFV, Jan 2000
     // and try to insert new point to improve BSpline interpolation
 
@@ -672,46 +706,46 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 
     Standard_Integer newcount = 0;
     for (ii = 0; ii < count; ii++) {
-      
+
       newpcons[newcount] = pcons[ii];
       newpc3d[newcount] = pc3d[ii];
       newcount++;
 
-      if(count - ii + newcount == MAX_ARRAY_SIZE) continue;
+      if(count - ii + newcount == aMaxArraySize) continue;
 
       BSplCLib::Eval(.5*(pc3d[ii]+pc3d[ii+1]), Standard_False, DerivativeRequest,
-                    extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
-                    
+        extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
+
       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]);
-       
-       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 {
+        Standard_Real ucons = 0.5*(pcons[ii]+pcons[ii+1]);
+        Standard_Real uc3d  = 0.5*(pc3d[ii]+pc3d[ii+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[ii] + deltamin && curp < pc3d[ii+1] - deltamin){
+            newpc3d[newcount] = curp;
+            newpcons[newcount] = ucons;
+            newcount ++;
+          }
+        }
+        else {
 #ifdef DEB 
-         // JAG
-         cout << "Projection not done" << endl;
+          // JAG
+          cout << "Projection not done" << endl;
 #endif
-       }
+        }
       }
-     
+
     }
 
     newpc3d[newcount] = pc3d[count];
@@ -724,14 +758,14 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
     pcons = newpcons;
     newpcons = temp;
 
-    if((count != newcount) && newcount < MAX_ARRAY_SIZE) { count = newcount; continue;}
+    if((count != newcount) && newcount < aMaxArraySize) { count = newcount; continue;}
 
     count = newcount;
 
     Standard_Real algtol = sqrt(besttol2);
 
     interpolok = Check (FlatKnots, Poles, count+1, Paramc3d, Paramcons, 
-                        myC3d, CurveOnSurface, algtol, tolsov);
+      myC3d, CurveOnSurface, algtol, tolsov);
 
     if (Precision::IsInfinite(algtol)) {
       mySameParameter = Standard_False;
@@ -743,15 +777,15 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 
     tolsov = algtol;
 
-    interpolok = (interpolok || count >= MAX_ARRAY_SIZE);
+    interpolok = (interpolok || count >= aMaxArraySize);
 
     if(interpolok) {
-        Standard_Real besttol = sqrt(besttol2);
+      Standard_Real besttol = sqrt(besttol2);
 #ifdef DEB
       if (Voir) {
-       if(algtol > besttol){
-         cout<<"SameParameter : Tol can't be reached before approx"<<endl;
-       }
+        if(algtol > besttol){
+          cout<<"SameParameter : Tol can't be reached before approx"<<endl;
+        }
       }
 #endif
       Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
@@ -761,58 +795,85 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
 
       Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d); 
       AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
-                                               Continuity,11,40,ev);
+        Continuity,11,40,ev);
 
       if (anApproximator.IsDone() || anApproximator.HasResult()) {
-       GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator) ;
-       myCurve2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
-       myHCurve2d = new Geom2dAdaptor_HCurve(myCurve2d);
-       CurveOnSurface.Load(myHCurve2d);
-       myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
-       myDone = Standard_True;
+        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);
+
+        myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
+
+        if(myTolReached > anErrorMAX)
+        {
+          //This tolerance is too big. Probably, we will not can get 
+          //edge with sameparameter in this case.
+
+          myDone = Standard_False;
+          return;
+        }
+
+        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  = new Geom2dAdaptor_HCurve(myCurve2d);
+          myDone = Standard_True;
+        }
+        else
+        {
+          interpolok = Standard_False;
+          CurveOnSurface = ACS;
+        }
       } 
     }
-    else {
+    
+    if(!interpolok)
+    {
 #ifdef DEB
       if (Voir)
-       cout<<"SameParameter : Not enough points, enrich"<<endl;
+        cout<<"SameParameter : Not enough points, enrich"<<endl;
 #endif
 
-      Standard_Integer newcount = 0;
+      newcount = 0;
       for(Standard_Integer n = 0; n < count; n++){
-       newpc3d[newcount] = pc3d[n];
-       newpcons[newcount] = pcons[n];
-       newcount ++;
-
-       if(count - n + newcount == MAX_ARRAY_SIZE) 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 {
+        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 DEB 
-         // JAG
-         cout << "Projection not done" << endl;
+          // JAG
+          cout << "Projection not done" << endl;
 #endif
-       }
+        }
       }
       newpc3d[newcount] = pc3d[count];
       newpcons[newcount] = pcons[count];
@@ -823,7 +884,16 @@ void Approx_SameParameter::Build(const Standard_Real Tolerance)
       tempx = pcons;
       pcons = newpcons;
       newpcons = tempx;
-      count = newcount;
+      
+      if(count != newcount)
+      {
+        count = newcount;
+        hasCountChanged = Standard_True;
+      }
+      else
+      {
+        hasCountChanged = Standard_False;
+      }
     }
   }
 }