src/Approx/Approx_SameParameter.cxx

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