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