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