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