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