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