0022048: Visualization, AIS_InteractiveContext - single object selection should alway...
[occt.git] / 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 }