src/Approx/Approx_SameParameter.cxx

   1 // Created on: 1995-06-06
   2 // Created by: Xavier BENVENISTE
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   5 //
   6 // The content of this file is subject to the Open CASCADE Technology Public
   7 // License Version 6.5 (the "License"). You may not use the content of this file
   8 // except in compliance with the License. Please obtain a copy of the License
   9 // at http://www.opencascade.org and read it completely before using this file.
  10 //
  11 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  12 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  13 //
  14 // The Original Code and all software distributed under the License is
  15 // distributed on an "AS IS" basis, without warranty of any kind, and the
  16 // Initial Developer hereby disclaims all such warranties, including without
  17 // limitation, any warranties of merchantability, fitness for a particular
  18 // purpose or non-infringement. Please see the License for the specific terms
  19 // and conditions governing the rights and limitations under the License.
  20
  21
  22 //  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898
  23
  24 #include <Approx_SameParameter.ixx>
  25 #include <TColStd_Array1OfReal.hxx>
  26 #include <BSplCLib.hxx>
  27 #include <Adaptor3d_CurveOnSurface.hxx>
  28 #include <Geom2dAdaptor_Curve.hxx>
  29 #include <Geom2dAdaptor_HCurve.hxx>
  30 #include <GeomAdaptor_Curve.hxx>
  31 #include <GeomAdaptor_HCurve.hxx>
  32 #include <GeomAdaptor_Surface.hxx>
  33 #include <GeomAdaptor_HSurface.hxx>
  34 //#include <GCPnts_UniformDeflection.hxx>
  35 #include <GCPnts_QuasiUniformDeflection.hxx>
  36 #include <Extrema_LocateExtPC.hxx>
  37 #include <AdvApprox_ApproxAFunction.hxx>
  38 #include <GeomLib_MakeCurvefromApprox.hxx>
  39 #include <Precision.hxx>
  40
  41 #define MAX_ARRAY_SIZE 1000 // IFV, Jan 2000
  42
  43 #ifdef DEB
  44 #ifdef DRAW
  45 #include <DrawTrSurf.hxx>
  46 #endif
  47 #include <Geom2d_BSplineCurve.hxx>
  48 #include <stdio.h>
  49 static Standard_Boolean Voir     = Standard_False;
  50 static Standard_Boolean AffichFw = Standard_False;
  51 static Standard_Integer NbCurve = 0;
  52 #endif
  53 //
  54 //   allows testing if Extrema produces correct results/
  55
  56
  57 static void ProjectPointOnCurve(const Standard_Real      InitValue,
  58                                 const gp_Pnt             APoint,
  59                                 const Standard_Real      Tolerance,
  60                                 const Standard_Integer   NumIteration,
  61                                 const Adaptor3d_Curve&     Curve,
  62                                 Standard_Boolean&        Status,
  63                                 Standard_Real&           Result)
  64 {
  65   Standard_Integer num_iter = 0,
  66   not_done = 1,
  67   ii ;
  68
  69   gp_Pnt a_point ;
  70   gp_Vec   vector,
  71   d1,
  72   d2 ;
  73   Standard_Real func,
  74   func_derivative,
  75   param = InitValue ;
  76   Status = Standard_False ;
  77   Standard_Real Toler = 1.0e-12;
  78   do {
  79     num_iter += 1 ;
  80     Curve.D2(param,
  81               a_point,
  82              d1,
  83              d2) ;
  84     for (ii = 1 ; ii <= 3 ; ii++) {
  85       vector.SetCoord(ii, APoint.Coord(ii) - a_point.Coord(ii)) ;
  86     }
  87     func = vector.Dot(d1) ;
  88     func_derivative = vector.Dot(d2) ;
  89     func_derivative -= d1.Dot(d1) ;
  90     if ( Abs(func) < Tolerance * d1.Magnitude()) {
  91       not_done = 0 ;
  92       Status = Standard_True ;
  93     }
  94     else
  95       { // fixing a bug PRO18577 : avoid divizion by zero
  96         if( Abs(func_derivative) > Toler )  {
  97           param -= func / func_derivative ;
  98         }
  99         param = Max(param,Curve.FirstParameter()) ;
 100         param = Min(param,Curve.LastParameter())  ;
 101         Status = Standard_True ;
 102       }
 103   }
 104   while (not_done && num_iter <= NumIteration) ;
 105   Result = param ;
 106 }
 107
 108
 109
 110 //=======================================================================
 111 //class : Approx_SameParameter_Evaluator
 112 //purpose  :
 113 //=======================================================================
 114
 115 class Approx_SameParameter_Evaluator : public AdvApprox_EvaluatorFunction
 116 {
 117  public:
 118   Approx_SameParameter_Evaluator (const TColStd_Array1OfReal& theFlatKnots,
 119                                   const TColStd_Array1OfReal& thePoles,
 120                                   const Handle(Adaptor2d_HCurve2d)& theHCurve2d)
 121     : FlatKnots(theFlatKnots), Poles(thePoles), HCurve2d(theHCurve2d) {}
 122
 123   virtual void Evaluate (Standard_Integer *Dimension,
 124                          Standard_Real     StartEnd[2],
 125                          Standard_Real    *Parameter,
 126                          Standard_Integer *DerivativeRequest,
 127                          Standard_Real    *Result, // [Dimension]
 128                          Standard_Integer *ErrorCode);
 129
 130  private:
 131   const TColStd_Array1OfReal& FlatKnots;
 132   const TColStd_Array1OfReal& Poles;
 133   Handle(Adaptor2d_HCurve2d) HCurve2d;
 134 };
 135
 136 void Approx_SameParameter_Evaluator::Evaluate (Standard_Integer *,/*Dimension*/
 137                                                Standard_Real    /*StartEnd*/[2],
 138                                                Standard_Real    *Parameter,
 139                                                Standard_Integer *DerivativeRequest,
 140                                                Standard_Real    *Result,
 141                                                Standard_Integer *ReturnCode)
 142 {
 143   gp_Pnt2d Point ;
 144   gp_Vec2d Vector ;
 145   Standard_Integer extrap_mode[2] ;
 146   extrap_mode[0] = extrap_mode[1] = 3;
 147   Standard_Real eval_result[2] ;
 148   Standard_Real *PolesArray =
 149     (Standard_Real *) &Poles(Poles.Lower()) ;
 150   //
 151   // evaluate the 1D bspline that represents the change in parameterization
 152   //
 153   BSplCLib::Eval(*Parameter,
 154                  Standard_False,
 155                  *DerivativeRequest,
 156                  extrap_mode[0],
 157                  3,
 158                  FlatKnots,
 159                  1,
 160                  PolesArray[0],
 161                  eval_result[0]) ;
 162
 163
 164   if (*DerivativeRequest == 0){
 165     HCurve2d->D0(eval_result[0],Point);
 166     Point.Coord(Result[0],Result[1]);
 167   }
 168   else if (*DerivativeRequest == 1){
 169     HCurve2d->D1(eval_result[0], Point, Vector);
 170     Vector.Multiply(eval_result[1]);
 171     Vector.Coord(Result[0],Result[1]);
 172   }
 173   ReturnCode[0] = 0 ;
 174 }
 175
 176 static Standard_Real ComputeTolReached(const Handle(Adaptor3d_HCurve)& c3d,
 177                                        const Adaptor3d_CurveOnSurface& cons,
 178                                        const Standard_Integer        nbp)
 179 {
 180   Standard_Real d2 = 0.;
 181   Standard_Integer nn = nbp;
 182   Standard_Real unsurnn = 1./nn;
 183   Standard_Real first = c3d->FirstParameter();
 184   Standard_Real last  = c3d->LastParameter();
 185   for(Standard_Integer i = 0; i <= nn; i++){
 186     Standard_Real t = unsurnn*i;
 187     Standard_Real u = first*(1.-t) + last*t;
 188     gp_Pnt Pc3d = c3d->Value(u);
 189     gp_Pnt Pcons = cons.Value(u);
 190     if (Precision::IsInfinite(Pcons.X()) ||
 191         Precision::IsInfinite(Pcons.Y()) ||
 192         Precision::IsInfinite(Pcons.Z())) {
 193       d2=Precision::Infinite();
 194       break;
 195     }
 196     Standard_Real temp = Pc3d.SquareDistance(Pcons);
 197     if(temp > d2) d2 = temp;
 198   }
 199   d2 = 1.5*sqrt(d2);
 200   if(d2<1.e-7) d2 = 1.e-7;
 201   return d2;
 202 }
 203
 204 static Standard_Boolean Check(const TColStd_Array1OfReal& FlatKnots,
 205                               const TColStd_Array1OfReal& Poles,
 206                               const Standard_Integer nbp,
 207                               const TColStd_Array1OfReal& pc3d,
 208 //                            const TColStd_Array1OfReal& pcons,
 209                               const TColStd_Array1OfReal& ,
 210                               const Handle(Adaptor3d_HCurve)& c3d,
 211                               const Adaptor3d_CurveOnSurface& cons,
 212                               Standard_Real& tol,
 213                               const Standard_Real oldtol)
 214 {
 215   Standard_Real d = tol;
 216   Standard_Integer extrap_mode[2] ;
 217   extrap_mode[0] = extrap_mode[1] = 3;
 218   Standard_Integer i;
 219 #ifdef DEB
 220   if (Voir) {
 221     cout<<endl;
 222     cout<<"Control the change of variable : "<<endl;
 223     cout<<"yawn mesured by projection : "<<d<<endl;
 224     cout<<"Number of points : "<<nbp<<endl;
 225   }
 226 #endif
 227 #if 0
 228   Standard_Real glis = 0., dglis = 0.;
 229   for(i = 1; i <= nbp; i++){
 230     Standard_Real tc3d = pc3d(i);
 231     gp_Pnt Pc3d = c3d->Value(tc3d);
 232     Standard_Real tcons;
 233     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
 234                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
 235     gp_Pnt Pcons = cons.Value(tcons);
 236     Standard_Real temp = Pc3d.SquareDistance(Pcons);
 237     if(temp >= dglis) dglis = temp;
 238     temp = Abs(tcons-pcons(i));
 239     if(temp >= glis) glis = temp;
 240   }
 241   dglis = sqrt(dglis);
 242 #ifdef DEB
 243   if ( Voir) {
 244     cout<<"shift of parameter to the imposed points : "<<glis<<endl;
 245     cout<<"shift distance at the imposed points : "<<dglis<<endl;
 246   }
 247 #endif
 248   dglis = 0.;
 249   for(i = 1; i < nbp; i++){
 250     Standard_Real tc3d = 0.5*(pc3d(i)+pc3d(i+1));
 251     gp_Pnt Pc3d = c3d->Value(tc3d);
 252     Standard_Real tcons;
 253     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
 254                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
 255     gp_Pnt Pcons = cons.Value(tcons);
 256     Standard_Real temp = Pc3d.SquareDistance(Pcons);
 257     if(temp >= dglis) dglis = temp;
 258   }
 259   dglis = sqrt(dglis);
 260 #ifdef DEB
 261   if (Voir)
 262     cout<<"distance de glissement en milieu d intervals : "<<dglis<<endl;
 263 #endif
 264 #endif
 265
 266   Standard_Real d2 = 0.;
 267   Standard_Integer nn = 2*nbp;
 268   Standard_Real unsurnn = 1./nn;
 269 //  Modified by skv - Wed Jun  2 11:49:59 2004 OCC5898 Begin
 270 // Correction of the interval of valid values. This condition has no sensible
 271 // grounds. But it is better then the old one (which is commented out) because
 272 // it fixes the bug OCC5898. To develop more or less sensible criterion it is
 273 // necessary to deeply investigate this problem which is not possible in frames
 274 // of debugging.
 275
 276 //   Standard_Real firstborne= 2*pc3d(1)-pc3d(nbp);
 277 //   Standard_Real lastborne= 2*pc3d(nbp)-pc3d(1);
 278   Standard_Real firstborne= 3.*pc3d(1)   - 2.*pc3d(nbp);
 279   Standard_Real lastborne = 3.*pc3d(nbp) - 2.*pc3d(1);
 280 //  Modified by skv - Wed Jun  2 11:50:03 2004 OCC5898 End
 281   for(i = 0; i <= nn; i++){
 282     Standard_Real t = unsurnn*i;
 283     Standard_Real tc3d = pc3d(1)*(1.-t) + pc3d(nbp)*t;
 284     gp_Pnt Pc3d = c3d->Value(tc3d);
 285     Standard_Real tcons;
 286     BSplCLib::Eval(tc3d,Standard_False,0,extrap_mode[0],
 287                    3,FlatKnots,1, (Standard_Real&)Poles(1),tcons);
 288     if (tcons < firstborne || tcons > lastborne) {
 289       tol=Precision::Infinite();
 290       return Standard_False;
 291     }
 292     gp_Pnt Pcons = cons.Value(tcons);
 293     Standard_Real temp = Pc3d.SquareDistance(Pcons);
 294     if(temp > d2) d2 = temp;
 295   }
 296   tol = sqrt(d2);
 297 #ifdef DEB
 298   if (Voir)
 299     cout<<"distance max on "<<nn<<" points : "<<tol<<endl<<endl;
 300 #endif
 301   return ((tol <= d) || (tol > 0.8 * oldtol));
 302 }
 303
 304
 305 //=======================================================================
 306 //function : Approx_SameParameter
 307 //purpose  :
 308 //=======================================================================
 309
 310 Approx_SameParameter::Approx_SameParameter(const Handle(Geom_Curve)&   C3D,
 311                                            const Handle(Geom2d_Curve)& C2D,
 312                                            const Handle(Geom_Surface)& S,
 313                                            const Standard_Real         Tol):
 314  mySameParameter(Standard_True), myDone(Standard_False)
 315 {
 316   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
 317   myC3d      = new GeomAdaptor_HCurve(C3D);
 318   mySurf     = new GeomAdaptor_HSurface(S);
 319   Build(Tol);
 320 }
 321
 322
 323 //=======================================================================
 324 //function : Approx_SameParameter
 325 //purpose  :
 326 //=======================================================================
 327
 328 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
 329                                            const Handle(Geom2d_Curve)&     C2D,
 330                                            const Handle(Adaptor3d_HSurface)& S,
 331                                            const Standard_Real            Tol):
 332  mySameParameter(Standard_True), myDone(Standard_False)
 333 {
 334   myC3d = C3D;
 335   mySurf = S;
 336   myHCurve2d = new Geom2dAdaptor_HCurve(C2D);
 337   Build(Tol);
 338 }
 339
 340
 341 //=======================================================================
 342 //function : Approx_SameParameter
 343 //purpose  :
 344 //=======================================================================
 345
 346 Approx_SameParameter::Approx_SameParameter(const Handle(Adaptor3d_HCurve)&   C3D,
 347                                            const Handle(Adaptor2d_HCurve2d)& C2D,
 348                                            const Handle(Adaptor3d_HSurface)& S,
 349                                            const Standard_Real            Tol):
 350  mySameParameter(Standard_True), myDone(Standard_False)
 351 {
 352   myC3d = C3D;
 353   mySurf = S;
 354   myHCurve2d = C2D;
 355   Build(Tol);
 356 }
 357
 358
 359 //=======================================================================
 360 //function : Build
 361 //purpose  :
 362 //=======================================================================
 363
 364 void Approx_SameParameter::Build(const Standard_Real Tolerance)
 365 {
 366   Standard_Integer ii ;
 367   Adaptor3d_CurveOnSurface CurveOnSurface(myHCurve2d,mySurf);
 368   Standard_Real fcons = CurveOnSurface.FirstParameter();
 369   Standard_Real lcons = CurveOnSurface.LastParameter();
 370   Standard_Real fc3d = myC3d->FirstParameter();
 371   Standard_Real lc3d = myC3d->LastParameter();
 372
 373   GeomAbs_Shape Continuity = myHCurve2d->Continuity();
 374
 375   if(Continuity > GeomAbs_C1) Continuity = GeomAbs_C1;
 376
 377   //Control tangents at the extremities to know if the
 378   //reparametring is possible and calculate the tangents
 379   //at the extremities of the function of change of variable.
 380   Standard_Real tangent[2] = { 0.0, 0.0 };
 381   gp_Pnt Pcons,Pc3d;
 382   gp_Vec Vcons,Vc3d;
 383
 384   Standard_Real Tol = Tolerance;
 385   Standard_Real Tol2 = Tol * Tol;
 386   Standard_Real Tolp = myC3d->Resolution(Tol), deltamin = 50*Tolp;
 387
 388   Standard_Real besttol2 = Tol2;
 389   Standard_Boolean extrok = 0;
 390
 391   extrok = 1;
 392   CurveOnSurface.D1(fcons,Pcons,Vcons);
 393   myC3d->D1(fc3d,Pc3d,Vc3d);
 394   Standard_Real dist2 = Pcons.SquareDistance(Pc3d);
 395   Standard_Real dmax2 = dist2;
 396
 397   Standard_Real magVcons = Vcons.Magnitude();
 398   if (magVcons > 1.e-12){
 399     tangent[0] = Vc3d.Magnitude() / magVcons;
 400   }
 401   else extrok = 0;
 402
 403   CurveOnSurface.D1(lcons,Pcons,Vcons);
 404   myC3d->D1(lc3d,Pc3d,Vc3d);
 405   dist2 = Pcons.SquareDistance(Pc3d);
 406
 407   if(dist2 > dmax2) dmax2 = dist2;
 408   magVcons = Vcons.Magnitude();
 409   if (magVcons > 1.e-12){
 410     tangent[1] = Vc3d.Magnitude() / magVcons;
 411   }
 412   else extrok = 0;
 413
 414
 415   if(dmax2 > besttol2) besttol2 = dmax2;
 416
 417   //Take a multiple of the sample pof CheckShape,
 418   //at least the control points will be correct. No comment!!!
 419
 420   Standard_Integer NCONTROL = 22;
 421 #ifdef DEB
 422   Standard_Integer nbcoups = 0;
 423 #endif
 424
 425   Standard_Boolean interpolok = 0;
 426   Standard_Real tolsov = 1.e200;
 427   //Take parameters with  constant step on the curve on surface
 428   //and on curve 3d.
 429   Standard_Real deltacons = lcons - fcons;
 430   deltacons /= (NCONTROL);
 431   Standard_Real deltac3d = lc3d - fc3d;
 432   deltac3d /= (NCONTROL);
 433
 434   Standard_Real wcons = fcons;
 435   Standard_Real wc3d  = fc3d;
 436
 437   Standard_Real qpcons[MAX_ARRAY_SIZE], qnewpcons[MAX_ARRAY_SIZE],
 438                 qpc3d[MAX_ARRAY_SIZE], qnewpc3d[MAX_ARRAY_SIZE];
 439   Standard_Real * pcons = qpcons; Standard_Real * newpcons = qnewpcons;
 440   Standard_Real * pc3d = qpc3d; Standard_Real * newpc3d = qnewpc3d;
 441
 442   for ( ii = 0 ; ii < NCONTROL; ii++) {
 443     pcons[ii] = wcons;
 444     pc3d[ii]  = wc3d;
 445     wcons += deltacons;
 446     wc3d  += deltac3d;
 447   }
 448   pcons[NCONTROL] = lcons;
 449   pc3d[NCONTROL]  = lc3d;
 450
 451   Standard_Integer New_NCONTROL = NCONTROL;
 452   if(Continuity < GeomAbs_C1) {
 453      Standard_Integer NbInt = myHCurve2d->NbIntervals(GeomAbs_C1) + 1;
 454      TColStd_Array1OfReal Param_de_decoupeC1 (1, NbInt);
 455      myHCurve2d->Intervals(Param_de_decoupeC1, GeomAbs_C1);
 456      TColStd_SequenceOfReal new_par;
 457      Standard_Integer inter = 1;
 458      ii =1;
 459      new_par.Append(fcons);
 460
 461      while(Param_de_decoupeC1(inter) <= fcons + deltamin) inter++;
 462      while(Param_de_decoupeC1(NbInt) >= lcons - deltamin) NbInt--;
 463
 464      while(inter <= NbInt || ii < NCONTROL) {
 465        if(Param_de_decoupeC1(inter) < pcons[ii]) {
 466          new_par.Append(Param_de_decoupeC1(inter));
 467          if((pcons[ii] - Param_de_decoupeC1(inter)) <= deltamin) {
 468            ii++;
 469            if(ii > NCONTROL) {ii = NCONTROL;}
 470          }
 471          inter++;
 472        }
 473        else {
 474          if((Param_de_decoupeC1(inter) - pcons[ii]) > deltamin) {
 475            new_par.Append(pcons[ii]);
 476          }
 477          ii++;
 478        }
 479      }
 480
 481      new_par.Append(lcons);
 482      New_NCONTROL = new_par.Length() - 1;
 483      //simple protection if New_NCONTROL > allocated elements in array
 484      if (New_NCONTROL > MAX_ARRAY_SIZE) {
 485        mySameParameter = Standard_False;
 486        return;
 487      }
 488      for(ii = 1; ii <= New_NCONTROL; ii++){
 489        pcons[ii] = pc3d[ii] = new_par.Value(ii + 1);
 490      }
 491      pc3d[New_NCONTROL]  = lc3d;
 492    }
 493
 494
 495   Extrema_LocateExtPC Projector;
 496   Projector.Initialize(myC3d->Curve(),fc3d,lc3d,Tol);
 497
 498   Standard_Integer count = 1;
 499   Standard_Real previousp = fc3d, initp=0, curp;//, deltamin = 50*Tolp;
 500   Standard_Real bornesup = lc3d - deltamin;
 501   Standard_Boolean projok = 0,
 502     use_parameter ;
 503   for (ii = 1; ii < New_NCONTROL; ii++){
 504     CurveOnSurface.D0(pcons[ii],Pcons);
 505     myC3d->D0(pc3d[ii],Pc3d);
 506     dist2 = Pcons.SquareDistance(Pc3d);
 507     use_parameter = (dist2 <= Tol2  && (pc3d[ii] > pc3d[count-1] + deltamin)) ;
 508     if(use_parameter) {
 509
 510       if(dist2 > dmax2) dmax2 = dist2;
 511       initp = previousp = pc3d[count] = pc3d[ii];
 512       pcons[count] = pcons[ii];
 513       count++;
 514     }
 515     else {
 516       if(!projok) initp = pc3d[ii];
 517       projok = mySameParameter = Standard_False;
 518       Projector.Perform(Pcons, initp);
 519       if (Projector.IsDone()) {
 520         curp = Projector.Point().Parameter();
 521         Standard_Real dist_2 = Projector.SquareDistance();
 522         if(dist_2 > besttol2) besttol2 = dist_2;
 523         projok = 1;
 524       }
 525       else {
 526         ProjectPointOnCurve(initp,Pcons,Tol,30,myC3d->Curve(),projok,curp);
 527       }
 528       if(projok){
 529         if(curp > previousp + deltamin && curp < bornesup){
 530           initp = previousp = pc3d[count] = curp;
 531           pcons[count] = pcons[ii];
 532           count++;
 533         }
 534       }
 535       else {
 536 #ifdef DEB
 537         // JAG
 538         cout << "Projection not done" << endl;
 539 #endif
 540       }
 541     }
 542   }
 543   if(mySameParameter){
 544     myTolReached = 1.5*sqrt(dmax2);
 545     return;
 546   }
 547
 548   if(!extrok) { // If not already SameP and tangent to mill, abandon.
 549     mySameParameter = Standard_False;
 550 #ifdef DEB
 551     cout<<"SameParameter problem  : zero tangent to extremities"<<endl;
 552 #endif
 553     return;
 554   }
 555
 556   pcons[count] = lcons;
 557   pc3d[count]  = lc3d;
 558
 559 #ifdef DEB
 560   if (AffichFw) {
 561     char Name[17];
 562     Name[0]='\0';
 563     TColgp_Array1OfPnt2d    DEBP2d  (0,count);
 564     TColStd_Array1OfInteger DEBMults(0,count);
 565     DEBMults.Init(1); DEBMults(0) = 2; DEBMults(count) = 2;
 566     TColStd_Array1OfReal    DEBKnots(0,count);
 567     for (Standard_Integer DEBi = 0; DEBi <= count; DEBi++) {
 568       DEBKnots(DEBi) = DEBi;
 569       DEBP2d  (DEBi) = gp_Pnt2d(pc3d[DEBi],pcons[DEBi]);
 570     }
 571     Handle(Geom2d_BSplineCurve) DEBBS =
 572       new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
 573     Sprintf(Name,"DEBC2d_%d",++NbCurve);
 574 #ifdef DRAW
 575     DrawTrSurf::Set(Name,DEBBS);
 576 #endif
 577   }
 578 #endif
 579
 580   while(!interpolok){
 581     // The tables and their limits for the interpolation.
 582     Standard_Integer num_knots = count + 7;
 583     Standard_Integer num_poles = count + 3;
 584     TColStd_Array1OfReal    Paramc3d(*pc3d,1,count+1);
 585     TColStd_Array1OfReal    Paramcons(*pcons,1,count+1);
 586     TColStd_Array1OfInteger ContactOrder(1,num_poles) ;
 587     TColStd_Array1OfReal    Poles(1,num_poles) ;
 588     TColStd_Array1OfReal    InterpolationParameters(1,num_poles) ;
 589     TColStd_Array1OfReal    FlatKnots(1,num_knots) ;
 590
 591     // Fill tables taking attention to end values.
 592     ContactOrder.Init(0);
 593     ContactOrder(2) = ContactOrder(num_poles - 1) = 1;
 594
 595     FlatKnots(1) = FlatKnots(2) = FlatKnots(3) = FlatKnots(4) = fc3d;
 596     FlatKnots(num_poles + 1) = FlatKnots(num_poles + 2) =
 597       FlatKnots(num_poles + 3) = FlatKnots(num_poles + 4) = lc3d;
 598
 599     Poles(1) = fcons; Poles(num_poles) = lcons;
 600     Poles(2) = tangent[0]; Poles(num_poles - 1) = tangent[1];
 601
 602     InterpolationParameters(1) = InterpolationParameters(2) = fc3d;
 603     InterpolationParameters(num_poles - 1) = InterpolationParameters(num_poles) = lc3d;
 604
 605     for (ii = 3; ii <= num_poles - 2; ii++) {
 606       Poles(ii) = Paramcons(ii - 1);
 607       InterpolationParameters(ii) = FlatKnots(ii+2) = Paramc3d(ii - 1);
 608     }
 609     Standard_Integer inversion_problem;
 610     BSplCLib::Interpolate(3,FlatKnots,InterpolationParameters,ContactOrder,
 611                           1,Poles(1),inversion_problem);
 612     if(inversion_problem) {
 613       Standard_ConstructionError::Raise();
 614     }
 615
 616     //-------------------------------------------
 617 #ifdef DEB
 618   if (AffichFw) {
 619     nbcoups ++;
 620     char Name[17];
 621     Name[0] = '\0';
 622     Standard_Integer nnn = 100;
 623     TColgp_Array1OfPnt2d    DEBP2d  (0,nnn);
 624     TColStd_Array1OfInteger DEBMults(0,nnn);
 625     DEBMults.Init(1); DEBMults(0) = 2; DEBMults(nnn) = 2;
 626     TColStd_Array1OfReal    DEBKnots(0,nnn);
 627     Standard_Real du = (lc3d - fc3d) / nnn;
 628     Standard_Real u3d = fc3d;
 629     Standard_Integer extrap_mode[2] ;
 630     extrap_mode[0] = extrap_mode[1] = 3;
 631     Standard_Real eval_result[2] ;
 632     Standard_Integer DerivativeRequest = 0;
 633     Standard_Real *PolesArray =
 634       (Standard_Real *) &Poles(Poles.Lower()) ;
 635
 636     for (Standard_Integer DEBi = 0; DEBi <= nnn; DEBi++) {
 637       DEBKnots(DEBi) = DEBi;
 638       BSplCLib::Eval(u3d,
 639                      Standard_False,
 640                      DerivativeRequest,
 641                      extrap_mode[0],
 642                      3,
 643                      FlatKnots,
 644                      1,
 645                      PolesArray[0],
 646                      eval_result[0]) ;
 647
 648       DEBP2d  (DEBi) = gp_Pnt2d(u3d,eval_result[0]);
 649       u3d += du;
 650     }
 651
 652     Handle(Geom2d_BSplineCurve) DEBBS =
 653       new Geom2d_BSplineCurve(DEBP2d,DEBKnots,DEBMults,1);
 654     Sprintf(Name,"DEBC2d_%d_%d",NbCurve,nbcoups );
 655 #ifdef DRAW
 656     DrawTrSurf::Set(Name,DEBBS);
 657 #endif
 658   }
 659 #endif
 660 //-------------------------------------------
 661
 662 //-------------------------------------------
 663     // Test if par2d(par3d) is monotonous function or not           ----- IFV, Jan 2000
 664     // and try to insert new point to improve BSpline interpolation
 665
 666     Standard_Integer extrap_mode[2] ;
 667     extrap_mode[0] = extrap_mode[1] = 3;
 668     Standard_Real eval_result[2] ;
 669     Standard_Integer DerivativeRequest = 0;
 670     Standard_Real *PolesArray =
 671       (Standard_Real *) &Poles(Poles.Lower()) ;
 672
 673     Standard_Integer newcount = 0;
 674     for (ii = 0; ii < count; ii++) {
 675
 676       newpcons[newcount] = pcons[ii];
 677       newpc3d[newcount] = pc3d[ii];
 678       newcount++;
 679
 680       if(count - ii + newcount == MAX_ARRAY_SIZE) continue;
 681
 682       BSplCLib::Eval(.5*(pc3d[ii]+pc3d[ii+1]), Standard_False, DerivativeRequest,
 683                      extrap_mode[0], 3, FlatKnots, 1, PolesArray[0], eval_result[0]);
 684
 685       if(eval_result[0] < pcons[ii] || eval_result[0] > pcons[ii+1]) {
 686         Standard_Real ucons = 0.5*(pcons[ii]+pcons[ii+1]);
 687         Standard_Real uc3d  = 0.5*(pc3d[ii]+pc3d[ii+1]);
 688
 689         CurveOnSurface.D0(ucons,Pcons);
 690         Projector.Perform(Pcons, uc3d);
 691         if (Projector.IsDone()) {
 692           curp = Projector.Point().Parameter();
 693           Standard_Real dist_2 = Projector.SquareDistance();
 694           if(dist_2 > besttol2) besttol2 = dist_2;
 695           projok = 1;
 696         }
 697         else {
 698           ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
 699         }
 700         if(projok){
 701           if(curp > pc3d[ii] + deltamin && curp < pc3d[ii+1] - deltamin){
 702             newpc3d[newcount] = curp;
 703             newpcons[newcount] = ucons;
 704             newcount ++;
 705           }
 706         }
 707         else {
 708 #ifdef DEB
 709           // JAG
 710           cout << "Projection not done" << endl;
 711 #endif
 712         }
 713       }
 714
 715     }
 716
 717     newpc3d[newcount] = pc3d[count];
 718     newpcons[newcount] = pcons[count];
 719     Standard_Real * temp;
 720     temp = pc3d;
 721     pc3d = newpc3d;
 722     newpc3d = temp;
 723     temp = pcons;
 724     pcons = newpcons;
 725     newpcons = temp;
 726
 727     if((count != newcount) && newcount < MAX_ARRAY_SIZE) { count = newcount; continue;}
 728
 729     count = newcount;
 730
 731     Standard_Real algtol = sqrt(besttol2);
 732
 733     interpolok = Check (FlatKnots, Poles, count+1, Paramc3d, Paramcons,
 734                         myC3d, CurveOnSurface, algtol, tolsov);
 735
 736     if (Precision::IsInfinite(algtol)) {
 737       mySameParameter = Standard_False;
 738 #ifdef DEB
 739       cout<<"SameParameter problem  : function of interpolation of parametration at mills !!"<<endl;
 740 #endif
 741       return;
 742     }
 743
 744     tolsov = algtol;
 745
 746     interpolok = (interpolok || count >= MAX_ARRAY_SIZE);
 747
 748     if(interpolok) {
 749         Standard_Real besttol = sqrt(besttol2);
 750 #ifdef DEB
 751       if (Voir) {
 752         if(algtol > besttol){
 753           cout<<"SameParameter : Tol can't be reached before approx"<<endl;
 754         }
 755       }
 756 #endif
 757       Handle(TColStd_HArray1OfReal) tol1d,tol2d,tol3d;
 758       tol1d = new TColStd_HArray1OfReal(1,2) ;
 759       tol1d->SetValue(1, mySurf->UResolution(besttol));
 760       tol1d->SetValue(2, mySurf->VResolution(besttol));
 761
 762       Approx_SameParameter_Evaluator ev (FlatKnots, Poles, myHCurve2d);
 763       AdvApprox_ApproxAFunction  anApproximator(2,0,0,tol1d,tol2d,tol3d,fc3d,lc3d,
 764                                                 Continuity,11,40,ev);
 765
 766       if (anApproximator.IsDone() || anApproximator.HasResult()) {
 767         GeomLib_MakeCurvefromApprox  aCurveBuilder(anApproximator) ;
 768         myCurve2d = aCurveBuilder.Curve2dFromTwo1d(1,2) ;
 769         myHCurve2d = new Geom2dAdaptor_HCurve(myCurve2d);
 770         CurveOnSurface.Load(myHCurve2d);
 771         myTolReached = ComputeTolReached(myC3d,CurveOnSurface,NCONTROL);
 772         myDone = Standard_True;
 773       }
 774     }
 775     else {
 776 #ifdef DEB
 777       if (Voir)
 778         cout<<"SameParameter : Not enough points, enrich"<<endl;
 779 #endif
 780
 781       Standard_Integer newcount = 0;
 782       for(Standard_Integer n = 0; n < count; n++){
 783         newpc3d[newcount] = pc3d[n];
 784         newpcons[newcount] = pcons[n];
 785         newcount ++;
 786
 787         if(count - n + newcount == MAX_ARRAY_SIZE) continue;
 788
 789         Standard_Real ucons = 0.5*(pcons[n]+pcons[n+1]);
 790         Standard_Real uc3d  = 0.5*(pc3d[n]+pc3d[n+1]);
 791
 792         CurveOnSurface.D0(ucons,Pcons);
 793         Projector.Perform(Pcons, uc3d);
 794         if (Projector.IsDone()) {
 795           curp = Projector.Point().Parameter();
 796           Standard_Real dist_2 = Projector.SquareDistance();
 797           if(dist_2 > besttol2) besttol2 = dist_2;
 798           projok = 1;
 799         }
 800         else {
 801           ProjectPointOnCurve(uc3d,Pcons,Tol,30,myC3d->Curve(),projok,curp);
 802         }
 803         if(projok){
 804           if(curp > pc3d[n] + deltamin && curp < pc3d[n+1] - deltamin){
 805             newpc3d[newcount] = curp;
 806             newpcons[newcount] = ucons;
 807             newcount ++;
 808           }
 809         }
 810         else {
 811 #ifdef DEB
 812           // JAG
 813           cout << "Projection not done" << endl;
 814 #endif
 815         }
 816       }
 817       newpc3d[newcount] = pc3d[count];
 818       newpcons[newcount] = pcons[count];
 819       Standard_Real * tempx;
 820       tempx = pc3d;
 821       pc3d = newpc3d;
 822       newpc3d = tempx;
 823       tempx = pcons;
 824       pcons = newpcons;
 825       newpcons = tempx;
 826       count = newcount;
 827     }
 828   }
 829 }