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