src/IntStart/IntStart_SearchOnBoundaries.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #include <algorithm>
  16 #include <TopoDS_Edge.hxx>
  17 #include <Geom_Curve.hxx>
  18 #include <BRepAdaptor_Curve.hxx>
  19 #include <Adaptor3d_HSurface.hxx>
  20 #include <GeomAbs_SurfaceType.hxx>
  21 #include <BRep_Tool.hxx>
  22 #include <Geom_Line.hxx>
  23 #include <gp_Lin.hxx>
  24 #include <gp_Vec.hxx>
  25 #include <gp_Dir.hxx>
  26 #include <gp_Cylinder.hxx>
  27 #include <gp_Ax1.hxx>
  28 #include <gp_Lin.hxx>
  29
  30 #include <GeomAdaptor_Curve.hxx>
  31 #include <Precision.hxx>
  32 #include <Extrema_ExtCC.hxx>
  33 #include <Extrema_POnCurv.hxx>
  34
  35 #include <math_FunctionSample.hxx>
  36 #include <math_FunctionAllRoots.hxx>
  37 #include <TColgp_SequenceOfPnt.hxx>
  38
  39 //  Modified by skv - Tue Aug 31 12:13:51 2004 OCC569
  40
  41 #include <Precision.hxx>
  42 #include <IntSurf_Quadric.hxx>
  43 #include <math_Function.hxx>
  44 #include <math_BrentMinimum.hxx>
  45 #include <math_Matrix.hxx>
  46 #include <math_Vector.hxx>
  47 #include <NCollection_Array1.hxx>
  48
  49 static void FindVertex (const TheArc&,
  50                         const Handle(TheTopolTool)&,
  51                         TheFunction&,
  52                         IntStart_SequenceOfPathPoint&,
  53                         const Standard_Real);
  54
  55
  56 static void BoundedArc (const TheArc& A,
  57                         const Handle(TheTopolTool)& Domain,
  58                         const Standard_Real Pdeb,
  59                         const Standard_Real Pfin,
  60                         TheFunction& Func,
  61                         IntStart_SequenceOfPathPoint& pnt,
  62                         IntStart_SequenceOfSegment& seg,
  63                         const Standard_Real TolBoundary,
  64                         const Standard_Real TolTangency,
  65                         Standard_Boolean& Arcsol,
  66                         const Standard_Boolean RecheckOnRegularity);
  67
  68 static void PointProcess (const gp_Pnt&,
  69                           const Standard_Real,
  70                           const TheArc&,
  71                           const Handle(TheTopolTool)&,
  72                           IntStart_SequenceOfPathPoint&,
  73                           const Standard_Real,
  74                           Standard_Integer&);
  75
  76 static Standard_Integer TreatLC (const TheArc& A,
  77                                  const Handle(TheTopolTool)& aDomain,
  78                                  const IntSurf_Quadric& aQuadric,
  79                                  const Standard_Real TolBoundary,
  80                                  IntStart_SequenceOfPathPoint& pnt);
  81
  82 static Standard_Boolean IsRegularity(const TheArc& A,
  83                                      const Handle(TheTopolTool)& aDomain);
  84
  85 class MinFunction : public math_Function
  86 {
  87 public:
  88   MinFunction(TheFunction &theFunc) : myFunc(&theFunc) {};
  89
  90   //returns value of the one-dimension-function when parameter
  91   //is equal to theX
  92   virtual Standard_Boolean Value(const Standard_Real theX,
  93                                  Standard_Real& theFVal)
  94   {
  95     if(!myFunc->Value(theX, theFVal))
  96       return Standard_False;
  97
  98     theFVal *= theFVal;
  99     return Standard_True;
 100   }
 101
 102   //see analogical method for abstract owner class math_Function
 103   virtual Standard_Integer GetStateNumber()
 104   {
 105     return 0;
 106   }
 107
 108 private:
 109   TheFunction *myFunc;
 110 };
 111
 112
 113 //=======================================================================
 114 //function : FindVertex
 115 //purpose  :
 116 //=======================================================================
 117 void FindVertex (const TheArc& A,
 118                  const Handle(TheTopolTool)& Domain,
 119                  TheFunction& Func,
 120                  IntStart_SequenceOfPathPoint& pnt,
 121                  const Standard_Real Toler)
 122 {
 123
 124 // Find the vertex of the arc A restriction solutions. It stores
 125 // Vertex in the list solutions pnt.
 126
 127
 128   TheVertex vtx;
 129   Standard_Real param,valf;
 130   Standard_Integer itemp;
 131
 132   Domain->Initialize(A);
 133   Domain->InitVertexIterator();
 134   while (Domain->MoreVertex()) {
 135     vtx = Domain->Vertex();
 136     param = TheSOBTool::Parameter(vtx,A);
 137
 138     // Evaluate the function and look compared to tolerance of the
 139     // Vertex. If distance <= tolerance then add a vertex to the list of solutions.
 140     // The arc is already assumed in the load function.
 141
 142     Func.Value(param,valf);
 143     if (Abs(valf) <= Toler) {
 144       itemp = Func.GetStateNumber();
 145       pnt.Append(IntStart_ThePathPoint(Func.Valpoint(itemp),Toler, vtx,A,param));
 146       // Solution is added
 147     }
 148     Domain->NextVertex();
 149   }
 150 }
 151
 152 class SolInfo
 153 {
 154 public:
 155   SolInfo() : myMathIndex(-1), myValue(RealLast())
 156   {
 157   }
 158
 159   void Init(const math_FunctionAllRoots& theSolution, const Standard_Integer theIndex)
 160   {
 161     myMathIndex = theIndex;
 162     myValue = theSolution.GetPoint(theIndex);
 163   }
 164
 165   Standard_Real Value() const
 166   {
 167     return myValue;
 168   }
 169
 170   Standard_Integer Index() const
 171   {
 172     return myMathIndex;
 173   }
 174
 175   bool operator>(const SolInfo& theOther) const
 176   {
 177     return myValue > theOther.myValue;
 178   }
 179
 180   bool operator<(const SolInfo& theOther) const
 181   {
 182     return myValue < theOther.myValue;
 183   }
 184
 185   bool operator==(const SolInfo& theOther) const
 186   {
 187     return myValue == theOther.myValue;
 188   }
 189
 190   Standard_Real& ChangeValue()
 191   {
 192     return myValue;
 193   }
 194
 195 private:
 196   Standard_Integer myMathIndex;
 197   Standard_Real myValue;
 198 };
 199
 200 static
 201 void BoundedArc (const TheArc& A,
 202                  const Handle(TheTopolTool)& Domain,
 203                  const Standard_Real Pdeb,
 204                  const Standard_Real Pfin,
 205                  TheFunction& Func,
 206                  IntStart_SequenceOfPathPoint& pnt,
 207                  IntStart_SequenceOfSegment& seg,
 208                  const Standard_Real TolBoundary,
 209                  const Standard_Real TolTangency,
 210                  Standard_Boolean& Arcsol,
 211                  const Standard_Boolean RecheckOnRegularity)
 212 {
 213   // Recherche des points solutions et des bouts d arc solution sur un arc donne.
 214   // On utilise la fonction math_FunctionAllRoots. Ne convient donc que pour
 215   // des arcs ayant un point debut et un point de fin (intervalle ferme de
 216   // parametrage).
 217
 218   Standard_Integer i,Nbi,Nbp;
 219
 220   gp_Pnt ptdeb,ptfin;
 221   Standard_Real pardeb = 0., parfin = 0.;
 222   Standard_Integer ideb,ifin,range,ranged,rangef;
 223
 224
 225   // Creer l echantillonage (math_FunctionSample ou classe heritant)
 226   // Appel a math_FunctionAllRoots
 227
 228   Standard_Real EpsX = TheArcTool::Resolution(A,Precision::Confusion());
 229   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 230   //@@@ La Tolerance est asociee a l arc  ( Incoherence avec le cheminement )
 231   //@@@   ( EpsX ~ 1e-5   et ResolutionU et V ~ 1e-9 )
 232   //@@@   le vertex trouve ici n'est pas retrouve comme point d arret d une
 233   //@@@   ligne de cheminement
 234   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 235   EpsX = 0.0000000001;
 236   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 237   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 238   //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 239
 240   //  Standard_Integer NbEchant = TheSOBTool::NbSamplesOnArc(A);
 241   Standard_Integer NbEchant = Func.NbSamples();
 242
 243   //-- Modif 24  Aout 93 -----------------------------
 244   Standard_Real nTolTangency = TolTangency;
 245   if((Pfin - Pdeb) < (TolTangency*10.0)) {
 246     nTolTangency=(Pfin-Pdeb)*0.1;
 247   }
 248   if(EpsX>(nTolTangency+nTolTangency)) {
 249     EpsX = nTolTangency * 0.1;
 250   }
 251   //--------------------------------------------------
 252   //-- Plante avec un edge avec 2 Samples
 253   //-- dont les extremites son solutions (f=0)
 254   //-- et ou la derivee est nulle
 255   //-- Exemple : un segment diametre d une sphere
 256   //-- if(NbEchant<3) NbEchant = 3; //-- lbr le 19 Avril 95
 257   //--------------------------------------------------
 258   Standard_Real para=0,dist,maxdist;
 259   /*  if(NbEchant<20) NbEchant = 20; //-- lbr le 22 Avril 96
 260   //-- Toujours des pbs
 261   */
 262   if(NbEchant<100) NbEchant = 100; //-- lbr le 22 Avril 96
 263   //-- Toujours des pbs
 264
 265   //-------------------------------------------------------------- REJECTIONS le 15 oct 98
 266   Standard_Boolean Rejection=Standard_True;
 267   Standard_Real maxdr,maxr,minr,ur,dur;
 268   minr=RealLast();
 269   maxr=-minr;
 270   maxdr=-minr;
 271   dur=(Pfin-Pdeb)*0.2;
 272   for(i=1,ur=Pdeb;i<=6;i++) {
 273     Standard_Real F,D;
 274     if(Func.Values(ur,F,D)) {
 275       Standard_Real lminr,lmaxr;
 276       if(D<0.0) D=-D;
 277       D*=dur+dur;
 278       if(D>maxdr) maxdr=D;
 279       lminr=F-D;
 280       lmaxr=F+D;
 281       if(lminr<minr) minr=lminr;
 282       if(lmaxr>maxr) maxr=lmaxr;
 283       if(minr<0.0 && maxr>0.0)  {
 284         Rejection=Standard_False;
 285         break;
 286       }
 287     }
 288     ur+=dur;
 289   }
 290   if(Rejection)
 291   {
 292     dur=0.001+maxdr+(maxr-minr)*0.1;
 293     minr-=dur;
 294     maxr+=dur;
 295     if(minr<0.0 && maxr>0.0)  {
 296       Rejection=Standard_False;
 297     }
 298   }
 299
 300   Arcsol=Standard_False;
 301
 302   if(Rejection==Standard_False) {
 303     math_FunctionSample Echant(Pdeb,Pfin,NbEchant);
 304
 305     Standard_Boolean aelargir=Standard_True;
 306     //modified by NIZNHY-PKV Thu Apr 12 09:25:19 2001 f
 307     //
 308     //maxdist = 100.0*TolBoundary;
 309     maxdist = TolBoundary+TolTangency;
 310     //
 311     //modified by NIZNHY-PKV Thu Apr 12 09:25:23 2001 t
 312     for(i=1; i<=NbEchant && aelargir;i++) {
 313       Standard_Real u = Echant.GetParameter(i);
 314       if(Func.Value(u,dist)) {
 315         if(dist>maxdist || -dist>maxdist) {
 316           aelargir=Standard_False;
 317         }
 318       }
 319     }
 320     if(!(aelargir && maxdist<0.01)) {
 321       maxdist = TolBoundary;
 322     }
 323
 324     math_FunctionAllRoots Sol(Func,Echant,EpsX,maxdist,maxdist); //-- TolBoundary,nTolTangency);
 325
 326     if (!Sol.IsDone()) {throw Standard_Failure();}
 327
 328     Nbp=Sol.NbPoints();
 329     //
 330     //jgv: build solution on the whole boundary
 331     if (RecheckOnRegularity && Nbp > 0 && IsRegularity(A, Domain))
 332     {
 333       //Standard_Real theTol = Domain->MaxTolerance(A);
 334       //theTol += theTol;
 335       Standard_Real theTol = 5.e-4;
 336       math_FunctionAllRoots SolAgain(Func,Echant,EpsX,theTol,theTol); //-- TolBoundary,nTolTangency);
 337
 338       if (!SolAgain.IsDone()) {throw Standard_Failure();}
 339
 340       Standard_Integer Nbi_again = SolAgain.NbIntervals();
 341
 342       if (Nbi_again > 0)
 343       {
 344         Standard_Integer NbSamples = 10;
 345         Standard_Real delta = (Pfin - Pdeb)/NbSamples;
 346         Standard_Real GlobalTol = theTol*10;
 347         Standard_Boolean SolOnBoundary = Standard_True;
 348         for (i = 0; i <= NbSamples; i++)
 349         {
 350           Standard_Real aParam = Pdeb + i*delta;
 351           Standard_Real aValue;
 352           Func.Value(aParam, aValue);
 353           if (Abs(aValue) > GlobalTol)
 354           {
 355             SolOnBoundary = Standard_False;
 356             break;
 357           }
 358         }
 359
 360         if (SolOnBoundary)
 361         {
 362           for (i = 1; i <= Nbi_again; i++)
 363           {
 364             IntStart_TheSegment newseg;
 365             newseg.SetValue(A);
 366             // Recuperer point debut et fin, et leur parametre.
 367             SolAgain.GetInterval(i,pardeb,parfin);
 368
 369             if (Abs(pardeb - Pdeb) <= Precision::PConfusion())
 370               pardeb = Pdeb;
 371             if (Abs(parfin - Pfin) <= Precision::PConfusion())
 372               parfin = Pfin;
 373
 374             SolAgain.GetIntervalState(i,ideb,ifin);
 375
 376             //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<"  ParFin:"<<parfin<<endl;
 377
 378             ptdeb=Func.Valpoint(ideb);
 379             ptfin=Func.Valpoint(ifin);
 380
 381             PointProcess(ptdeb,pardeb,A,Domain,pnt,theTol,ranged);
 382             newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
 383             PointProcess(ptfin,parfin,A,Domain,pnt,theTol,rangef);
 384             newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
 385             seg.Append(newseg);
 386           }
 387           Arcsol=Standard_True;
 388           return;
 389         }
 390       }
 391     }
 392     ////////////////////////////////////////////
 393
 394     //-- detection du cas ou la fonction est quasi tangente et que les
 395     //-- zeros sont quasi confondus.
 396     //-- Dans ce cas on prend le point "milieu"
 397     //-- On suppose que les solutions sont triees.
 398
 399     if(Nbp) {
 400       NCollection_Array1<SolInfo> aSI(1, Nbp);
 401
 402       for(i=1;i<=Nbp;i++)
 403       {
 404         aSI(i).Init(Sol, i);
 405       }
 406
 407       std::sort(aSI.begin(), aSI.end());
 408
 409       //modified by NIZNHY-PKV Wed Mar 21 18:34:18 2001 f
 410       //////////////////////////////////////////////////////////
 411       // The treatment of the situation when line(arc) that is
 412       // tangent to cylinder(domain).
 413       // We should have only one solution i.e Nbp=1. Ok?
 414       // But we have 2,3,.. solutions.     That is wrong ersult.
 415       // The TreatLC(...) function is dedicated to solve the pb.
 416       //                           PKV Fri Mar 23 12:17:29 2001
 417       Standard_Integer ip;
 418       const IntSurf_Quadric& aQuadric=Func.Quadric();
 419
 420       ip=TreatLC (A, Domain, aQuadric, TolBoundary, pnt);
 421       if (ip) {
 422         //////////////////////////////////////////////////////////
 423         //modified by NIZNHY-PKV Wed Mar 21 18:34:23 2001 t
 424         //
 425         // Using of old usual way proposed by Laurent
 426         //
 427         for(i=1;i<Nbp;i++) {
 428           Standard_Real parap1 = aSI(i + 1).Value();
 429           para = aSI(i).Value();
 430
 431           Standard_Real param=(para+parap1)*0.5;
 432           Standard_Real ym;
 433           if(Func.Value(param,ym)) {
 434             if(Abs(ym)<maxdist) {
 435               //  Modified by skv - Tue Aug 31 12:13:51 2004 OCC569 Begin
 436               // Consider this interval as tangent one. Treat it to find
 437               // parameter with the lowest function value.
 438
 439               // Compute the number of nodes.
 440               Standard_Real    aTol = TolBoundary*1000.0;
 441               if(aTol > 0.001)
 442                 aTol = 0.001;
 443
 444               // fix floating point exception 569, chl-922-e9
 445               parap1 = (Abs(parap1) < 1.e9) ? parap1 : ((parap1 >= 0.) ? 1.e9 : -1.e9);
 446               para   = (Abs(para) < 1.e9) ? para : ((para >= 0.) ? 1.e9 : -1.e9);
 447
 448               Standard_Integer aNbNodes = RealToInt(Ceiling((parap1 - para)/aTol));
 449
 450               Standard_Real    aVal     = RealLast();
 451               //Standard_Integer aNbNodes = 23;
 452               Standard_Real    aDelta   = (parap1 - para)/(aNbNodes + 1.);
 453               Standard_Integer ii;
 454               Standard_Real    aCurPar;
 455               Standard_Real    aCurVal;
 456
 457               for (ii = 0; ii <= aNbNodes + 1; ii++) {
 458                 aCurPar = (ii < aNbNodes + 1) ? para + ii*aDelta : parap1;
 459
 460                 if (Func.Value(aCurPar, aCurVal)) {
 461                   //if (aCurVal < aVal) {
 462                   if (Abs(aCurVal) < aVal) {
 463                     //aVal  = aCurVal;
 464                     aVal  = Abs(aCurVal);
 465                     param = aCurPar;
 466                   }
 467                 }
 468               }
 469               //  Modified by skv - Tue Aug 31 12:13:51 2004 OCC569 End
 470               aSI(i).ChangeValue() = Pdeb - 1;
 471               aSI(i + 1).ChangeValue() = param;
 472             }
 473           }
 474         }
 475
 476         for (i=1; i<=Nbp; i++) {
 477           para = aSI(i).Value();
 478           if((para-Pdeb)<EpsX || (Pfin-para)<EpsX)
 479             continue;
 480
 481           if(!Func.Value(para,dist))
 482             continue;
 483
 484           dist = Abs(dist);
 485
 486           Standard_Integer anIndx = -1;
 487           const Standard_Real aParam = Sol.GetPoint(aSI(i).Index());
 488           if (dist < maxdist)
 489           {
 490             if (Abs(aParam - Pdeb) <= Precision::PConfusion() || Abs(aParam - Pfin) <= Precision::PConfusion())
 491             {
 492               Standard_Real aDistTemp = RealLast();
 493               if (Func.Value(aParam, aDistTemp))
 494               {
 495                 if (Abs(aDistTemp) < maxdist)
 496                 {
 497                   anIndx = Sol.GetPointState(aSI(i).Index());
 498                 }
 499               }
 500             }
 501           }
 502
 503           gp_Pnt aPnt(anIndx < 0 ? Func.LastComputedPoint() : Func.Valpoint(anIndx));
 504
 505           if (dist > 0.1*Precision::Confusion())
 506           {
 507             //Precise found points. It results in following:
 508             //  1. Make the vertex nearer to the intersection line
 509             //    (see description to issue #27252 in order to
 510             //    understand necessity).
 511             //  2. Merge two near vertices to single point.
 512
 513             //All members in TabSol array has already been sorted in increase order.
 514             //Now, we limit precise boundaries in order to avoid changing this order.
 515             const Standard_Real aFPar = (i == 1) ? Pdeb : (para + aSI(i - 1).Value()) / 2.0;
 516             const Standard_Real aLPar = (i == Nbp) ? Pfin : (para + aSI(i + 1).Value()) / 2.0;
 517
 518             MinFunction aNewFunc(Func);
 519             math_BrentMinimum aMin(Precision::Confusion());
 520
 521             aMin.Perform(aNewFunc, aFPar, para, aLPar);
 522             if(aMin.IsDone())
 523             {
 524               para = aMin.Location();
 525               const gp_Pnt2d aP2d(A->Value(para));
 526               aPnt = Func.Surface()->Value(aP2d.X(), aP2d.Y());
 527             }
 528           }
 529
 530           PointProcess(aPnt, para, A, Domain, pnt, TolBoundary, range);
 531         }
 532       }// end of if(ip)
 533     } // end of if(Nbp)
 534
 535     // Pour chaque intervalle trouve faire
 536     //   Traiter les extremites comme des points
 537     //   Ajouter intervalle dans la liste des segments
 538
 539     Nbi=Sol.NbIntervals();
 540
 541
 542     if (!RecheckOnRegularity && Nbp) {
 543       //--cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx :Nbp>0  0 <- Nbi "<<Nbi<<endl;
 544       Nbi=0;
 545     }
 546
 547     //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : Nbi : "<<Nbi<<endl;
 548
 549     for (i=1; i<=Nbi; i++) {
 550       IntStart_TheSegment newseg;
 551       newseg.SetValue(A);
 552       // Recuperer point debut et fin, et leur parametre.
 553       Sol.GetInterval(i,pardeb,parfin);
 554       Sol.GetIntervalState(i,ideb,ifin);
 555
 556
 557       //-- cout<<" Debug : IntStart_SearchOnBoundaries_1.gxx : i= "<<i<<" ParDeb:"<<pardeb<<"  ParFin:"<<parfin<<endl;
 558
 559       ptdeb=Func.Valpoint(ideb);
 560       ptfin=Func.Valpoint(ifin);
 561
 562       PointProcess(ptdeb,pardeb,A,Domain,pnt,TolBoundary,ranged);
 563       newseg.SetLimitPoint(pnt.Value(ranged),Standard_True);
 564       PointProcess(ptfin,parfin,A,Domain,pnt,TolBoundary,rangef);
 565       newseg.SetLimitPoint(pnt.Value(rangef),Standard_False);
 566       seg.Append(newseg);
 567     }
 568
 569     if (Nbi==1) {
 570       if((Abs(pardeb - Pdeb) < Precision::PConfusion()) &&
 571          (Abs(parfin - Pfin) < Precision::PConfusion()))
 572       {
 573         Arcsol=Standard_True;
 574       }
 575     }
 576   }
 577 }
 578
 579 //=======================================================================
 580 //function : ComputeBoundsfromInfinite
 581 //purpose  :
 582 //=======================================================================
 583 // - PROVISIONAL - TEMPORARY - NOT GOOD - NYI - TO DO
 584 // - Temporary - temporary - not good - nyi - to do
 585 void ComputeBoundsfromInfinite(TheFunction& Func,
 586                                Standard_Real& PDeb,
 587                                Standard_Real& PFin,
 588                                Standard_Integer& NbEchant)
 589 {
 590
 591   // - We are looking for parameters for start and end of the arc (2d curve)
 592   // - Infinity, a way to intersect the quadric with a portion of arc
 593   // - Finished.
 594   //
 595   // - The quadric is a plane, a cylinder, a cone and a sphere.
 596   // - Idea: We take any point on the arc and the fact grow
 597   // - Terminals to the signed distance function values or is likely
 598   // - S cancel.
 599   //
 600   // - WARNING: The following calculations provide a very estimated coarse parameters.
 601   // - This avoids the raises and allows a case of Boxes
 602   // - Inifinies walk. It will take this code
 603   // - With curve surface intersections.
 604
 605   NbEchant = 100;
 606
 607   Standard_Real U0    = 0.0;
 608   Standard_Real dU    = 0.001;
 609   Standard_Real Dist0,Dist1;
 610
 611   Func.Value(U0   , Dist0);
 612   Func.Value(U0+dU, Dist1);
 613   Standard_Real dDist = Dist1 - Dist0;
 614   if(dDist) {
 615     U0  -=  dU*Dist0 / dDist;
 616     PDeb = PFin = U0;
 617     Standard_Real Umin = U0 - 1e5;
 618     Func.Value(Umin   , Dist0);
 619     Func.Value(Umin+dU, Dist1);
 620     dDist = Dist1-Dist0;
 621     if(dDist) {
 622       Umin  -=  dU*Dist0 / dDist;
 623     }
 624     else {
 625       Umin-=10.0;
 626     }
 627     Standard_Real Umax = U0 + 1e8;
 628     Func.Value(Umax   , Dist0);
 629     Func.Value(Umax+dU, Dist1);
 630     dDist = Dist1-Dist0;
 631     if(dDist) {
 632       Umax  -=  dU*Dist0 / dDist;
 633     }
 634     else {
 635       Umax+=10.0;
 636     }
 637     if(Umin>U0) { Umin=U0-10.0; }
 638     if(Umax<U0) { Umax=U0+10.0; }
 639
 640     PFin = Umax + 10. * (Umax - Umin);
 641     PDeb = Umin - 10. * (Umax - Umin);
 642   }
 643   else {
 644     //-- Possibilite de Arc totalement inclu ds Quad
 645     PDeb = 1e10;
 646     PFin = -1e10;
 647   }
 648 }
 649
 650 //=======================================================================
 651 //function : PointProcess
 652 //purpose  :
 653 //=======================================================================
 654 void PointProcess (const gp_Pnt& Pt,
 655                    const Standard_Real Para,
 656                    const TheArc& A,
 657                    const Handle(TheTopolTool)& Domain,
 658                    IntStart_SequenceOfPathPoint& pnt,
 659                    const Standard_Real Tol,
 660                    Standard_Integer& Range)
 661 {
 662
 663 // Check to see if a solution point is coincident with a vertex.
 664 // If confused, you should find this vertex in the list of
 665 // Start. It then returns the position of this point in the list pnt.
 666 // Otherwise, add the point in the list.
 667
 668   Standard_Integer k;
 669   Standard_Boolean found,goon;
 670   Standard_Real dist,toler;
 671
 672   Standard_Integer Nbsol = pnt.Length();
 673   TheVertex vtx;
 674   IntStart_ThePathPoint ptsol;
 675
 676   Domain->Initialize(A);
 677   Domain->InitVertexIterator();
 678   found = Standard_False;
 679   goon = Domain->MoreVertex();
 680   while (goon) {
 681     vtx = Domain->Vertex();
 682     dist= Abs(Para-TheSOBTool::Parameter(vtx,A));
 683     toler = TheSOBTool::Tolerance(vtx,A);
 684 #ifdef OCCT_DEBUG
 685     if(toler>0.1) {
 686       cout<<"IntStart_SearchOnBoundaries_1.gxx  : ** WARNING ** Tol Vertex="<<toler<<endl;
 687       cout<<"                                     Ou Edge degenere Ou Kro pointu"<<endl;
 688       if(toler>10000) toler=1e-7;
 689     }
 690 #endif
 691
 692     if (dist <= toler) {
 693       // Locate the vertex in the list of solutions
 694       k=1;
 695       found = (k>Nbsol);
 696       while (!found) {
 697         ptsol = pnt.Value(k);
 698         if (!ptsol.IsNew()) {
 699         //jag 940608  if (ptsol.Vertex() == vtx && ptsol.Arc()    == A) {
 700           if (Domain->Identical(ptsol.Vertex(),vtx) &&
 701                     ptsol.Arc()    == A &&
 702                     Abs(ptsol.Parameter()-Para) <= toler) {
 703             found=Standard_True;
 704           }
 705           else {
 706             k=k+1;
 707             found=(k>Nbsol);
 708           }
 709         }
 710         else {
 711           k=k+1;
 712           found=(k>Nbsol);
 713         }
 714       }
 715       if (k<=Nbsol) {     // We find the vertex
 716         Range = k;
 717       }
 718       else {              // Otherwise
 719         ptsol.SetValue(Pt,Tol,vtx,A,Para);
 720         pnt.Append(ptsol);
 721         Range = pnt.Length();
 722       }
 723       found = Standard_True;
 724       goon = Standard_False;
 725     }
 726     else {
 727       Domain->NextVertex();
 728       goon = Domain->MoreVertex();
 729     }
 730   }
 731
 732   if (!found) {   // No one is falling on a vertex
 733     //jgv: do not add segment's extremities if they already exist
 734     Standard_Boolean found_internal = Standard_False;
 735     for (k = 1; k <= pnt.Length(); k++)
 736     {
 737       ptsol = pnt.Value(k);
 738       if (ptsol.Arc() != A ||
 739           !ptsol.IsNew()) //vertex
 740         continue;
 741       if (Abs(ptsol.Parameter()-Para) <= Precision::PConfusion())
 742       {
 743         found_internal = Standard_True;
 744         Range = k;
 745       }
 746     }
 747     /////////////////////////////////////////////////////////////
 748
 749     if (!found_internal)
 750     {
 751       Standard_Real TOL=Tol;
 752       TOL*=1000.0;
 753       //if(TOL>0.001) TOL=0.001;
 754       if(TOL>0.005) TOL=0.005; //#24643
 755
 756       ptsol.SetValue(Pt,TOL,A,Para);
 757       pnt.Append(ptsol);
 758       Range = pnt.Length();
 759     }
 760   }
 761 }
 762
 763 //=======================================================================
 764 //function : IsRegularity
 765 //purpose  :
 766 //=======================================================================
 767 Standard_Boolean IsRegularity(const TheArc& /*A*/,
 768                               const Handle(TheTopolTool)& aDomain)
 769 {
 770   Standard_Address anEAddress=aDomain->Edge();
 771   if (anEAddress==NULL) {
 772     return Standard_False;
 773   }
 774
 775   TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;
 776
 777   return (BRep_Tool::HasContinuity(*anE));
 778 }
 779
 780 //=======================================================================
 781 //function : TreatLC
 782 //purpose  :
 783 //=======================================================================
 784 Standard_Integer TreatLC (const TheArc& A,
 785                           const Handle(TheTopolTool)& aDomain,
 786                           const IntSurf_Quadric& aQuadric,
 787                           const Standard_Real TolBoundary,
 788                           IntStart_SequenceOfPathPoint& pnt)
 789 {
 790   Standard_Integer anExitCode=1, aNbExt;
 791
 792   Standard_Address anEAddress=aDomain->Edge();
 793   if (anEAddress==NULL) {
 794     return anExitCode;
 795   }
 796
 797   TopoDS_Edge* anE=(TopoDS_Edge*)anEAddress;
 798
 799   if (BRep_Tool::Degenerated(*anE)) {
 800     return anExitCode;
 801   }
 802
 803   GeomAbs_CurveType   aTypeE;
 804   BRepAdaptor_Curve aBAC(*anE);
 805   aTypeE=aBAC.GetType();
 806
 807   if (aTypeE!=GeomAbs_Line) {
 808     return anExitCode;
 809   }
 810
 811   GeomAbs_SurfaceType aTypeS;
 812   aTypeS=aQuadric.TypeQuadric();
 813
 814   if (aTypeS!=GeomAbs_Cylinder) {
 815     return anExitCode;
 816   }
 817
 818   Standard_Real f, l, U1f, U1l, U2f, U2l, UEgde, TOL, aDist, aR, aRRel, Tol;
 819   Handle(Geom_Curve) aCEdge=BRep_Tool::Curve(*anE, f, l);
 820
 821   gp_Cylinder aCyl=aQuadric.Cylinder();
 822   const gp_Ax1& anAx1=aCyl.Axis();
 823   gp_Lin aLin(anAx1);
 824   Handle(Geom_Line) aCAxis=new Geom_Line (aLin);
 825   aR=aCyl.Radius();
 826
 827   U1f = aCAxis->FirstParameter();
 828   U1l = aCAxis->LastParameter();
 829
 830   U2f = aCEdge->FirstParameter();
 831   U2l = aCEdge->LastParameter();
 832
 833
 834   GeomAdaptor_Curve C1, C2;
 835
 836   C1.Load(aCAxis);
 837   C2.Load(aCEdge);
 838
 839   Tol = Precision::PConfusion();
 840
 841   Extrema_ExtCC anExtCC(C1, C2, U1f, U1l, U2f, U2l, Tol, Tol);
 842
 843   aNbExt=anExtCC.NbExt();
 844   if (aNbExt!=1) {
 845     return anExitCode;
 846   }
 847
 848   gp_Pnt P1,PEdge;
 849   Extrema_POnCurv PC1, PC2;
 850
 851   anExtCC.Points(1, PC1, PC2);
 852
 853   P1   =PC1.Value();
 854   PEdge=PC2.Value();
 855
 856   UEgde=PC2.Parameter();
 857
 858   aDist=PEdge.Distance(P1);
 859   aRRel=fabs(aDist-aR)/aR;
 860   if (aRRel > TolBoundary) {
 861     return anExitCode;
 862   }
 863
 864   if (UEgde < (f+TolBoundary) || UEgde > (l-TolBoundary)) {
 865     return anExitCode;
 866   }
 867   //
 868   // Do not wonder !
 869   // It was done as into PointProcess(...) function
 870   //printf("TreatLC()=> tangent line is found\n");
 871   TOL=1000.*TolBoundary;
 872   if(TOL>0.001) TOL=0.001;
 873
 874   IntStart_ThePathPoint ptsol;
 875   ptsol.SetValue(PEdge, TOL, A, UEgde);
 876   pnt.Append(ptsol);
 877
 878   anExitCode=0;
 879   return anExitCode;
 880
 881 }
 882
 883
 884 //=======================================================================
 885 //function : IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries
 886 //purpose  :
 887 //=======================================================================
 888 IntStart_SearchOnBoundaries::IntStart_SearchOnBoundaries ()
 889 :  done(Standard_False)
 890 {
 891 }
 892
 893 //=======================================================================
 894 //function : Perform
 895 //purpose  :
 896 //=======================================================================
 897   void IntStart_SearchOnBoundaries::Perform (TheFunction& Func,
 898                                              const Handle(TheTopolTool)& Domain,
 899                                              const Standard_Real TolBoundary,
 900                                              const Standard_Real TolTangency,
 901                                              const Standard_Boolean RecheckOnRegularity)
 902 {
 903
 904   done = Standard_False;
 905   spnt.Clear();
 906   sseg.Clear();
 907
 908   Standard_Boolean Arcsol;
 909   Standard_Real PDeb,PFin, prm, tol;
 910   Standard_Integer i, nbknown, nbfound,index;
 911   gp_Pnt pt;
 912
 913   Domain->Init();
 914
 915   if (Domain->More()) {
 916     all  = Standard_True;
 917   }
 918   else {
 919     all = Standard_False;
 920   }
 921
 922   while (Domain->More()) {
 923     TheArc A = Domain->Value();
 924     if (!TheSOBTool::HasBeenSeen(A)) {
 925       Func.Set(A);
 926       FindVertex(A,Domain,Func,spnt,TolBoundary);
 927       TheSOBTool::Bounds(A,PDeb,PFin);
 928       if(Precision::IsNegativeInfinite(PDeb) ||
 929          Precision::IsPositiveInfinite(PFin)) {
 930         Standard_Integer NbEchant;
 931         ComputeBoundsfromInfinite(Func,PDeb,PFin,NbEchant);
 932       }
 933       BoundedArc(A,Domain,PDeb,PFin,Func,spnt,sseg,TolBoundary,TolTangency,Arcsol,RecheckOnRegularity);
 934       all = (all && Arcsol);
 935     }
 936
 937     else {
 938       // as it seems we'll never be here, because
 939       // TheSOBTool::HasBeenSeen(A) always returns FALSE
 940       nbfound = spnt.Length();
 941
 942       // On recupere les points connus
 943       nbknown = TheSOBTool::NbPoints(A);
 944       for (i=1; i<=nbknown; i++) {
 945         TheSOBTool::Value(A,i,pt,tol,prm);
 946         if (TheSOBTool::IsVertex(A,i)) {
 947           TheVertex vtx;
 948           TheSOBTool::Vertex(A,i,vtx);
 949           spnt.Append(IntStart_ThePathPoint(pt,tol,vtx,A,prm));
 950         }
 951         else {
 952           spnt.Append(IntStart_ThePathPoint(pt,tol,A,prm));
 953         }
 954       }
 955       // On recupere les arcs solutions
 956       nbknown = TheSOBTool::NbSegments(A);
 957       for (i=1; i<=nbknown; i++) {
 958         IntStart_TheSegment newseg;
 959         newseg.SetValue(A);
 960         if (TheSOBTool::HasFirstPoint(A,i,index)) {
 961           newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_True);
 962         }
 963         if (TheSOBTool::HasLastPoint(A,i,index)) {
 964           newseg.SetLimitPoint(spnt.Value(nbfound+index),Standard_False);
 965         }
 966         sseg.Append(newseg);
 967       }
 968       all = (all& TheSOBTool::IsAllSolution(A));
 969     }
 970     Domain->Next();
 971   }
 972   done = Standard_True;
 973 }