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