36cf841c77c9460388f2ace78ca8f8bdeca138a3
[occt.git] / src / ShapeAnalysis / ShapeAnalysis_Curve.cxx
1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
2 //
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
7 //
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
10 //
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
17
18 // pdn 04.12.98 Add method using Adaptor_Curve
19 //:j8 abv 10.12.98 TR10 r0501_db.stp #9423
20 //pdn 25.12.98 private method ProjectAct
21 //szv#4 S4163
22 //:s5 abv 22.04.99 Adding debug printouts in catch {} blocks
23 //    abv 14.05.99 S4174: Adding method for exact computing of the boundary box 
24 //    gka 21.06.99 S4208: adding method NextProject(Adaptor_Curve)
25 //    msv 30.05.00 correct IsPlanar for a conic curve
26 #include <ShapeAnalysis_Curve.ixx>
27
28 #include <ElCLib.hxx>
29
30 #include <Geom2d_BoundedCurve.hxx>
31 #include <Geom2d_Line.hxx>
32 #include <Geom_BSplineCurve.hxx>
33 #include <GeomAdaptor_Curve.hxx>
34
35 #include <Precision.hxx>
36
37 #include <Standard_ErrorHandler.hxx>
38 #include <Standard_Failure.hxx>
39 #include <Adaptor3d_Curve.hxx>
40 #include <Extrema_ExtPC.hxx>
41 #include <ShapeAnalysis.hxx>
42 #include <TColgp_SequenceOfPnt.hxx>
43 #include <Geom_Line.hxx>
44 #include <Geom_Conic.hxx>
45 #include <Geom_TrimmedCurve.hxx>
46 #include <Geom_OffsetCurve.hxx>
47 #include <Geom_BezierCurve.hxx>
48 #include <ShapeExtend_ComplexCurve.hxx>
49 #include <Geom2d_Conic.hxx>
50 #include <Geom2d_TrimmedCurve.hxx>
51 #include <Geom2d_BSplineCurve.hxx>
52 #include <Geom2d_BezierCurve.hxx>
53
54 #include <Geom2d_OffsetCurve.hxx>
55 #include <Geom2dInt_Geom2dCurveTool.hxx>
56 #include <Geom2dAdaptor_Curve.hxx>
57 #include <Geom_Circle.hxx>
58
59
60 //=======================================================================
61 //function : ProjectOnSegments
62 //purpose  : 
63 //=======================================================================
64
65 static void ProjectOnSegments (const Adaptor3d_Curve& AC, const gp_Pnt& P3D,
66                                const Standard_Integer nbseg,
67                                Standard_Real& uMin, Standard_Real& uMax,
68                                Standard_Real& distmin, gp_Pnt& proj, Standard_Real& param)
69 {
70   //  On considere <nbseg> points sur [uMin,uMax]
71   //  Quel est le plus proche. Et quel est le nouvel intervalle
72   //  (il ne peut pas deborder l ancien)
73   Standard_Real u, dist, delta = (nbseg == 0)? 0 : (uMax-uMin)/nbseg; //szv#4:S4163:12Mar99 anti-exception
74   for (Standard_Integer i = 0; i <= nbseg; i ++) {
75     u = uMin + (delta * i);
76     gp_Pnt PU = AC.Value (u);
77     dist = PU.Distance (P3D);
78     if (dist < distmin)  {  distmin = dist;  proj = PU;  param = u;  }
79   }
80
81 #ifdef DEBUG
82   cout<<"ShapeAnalysis_Geom:Project, param="<<param<<" -> distmin="<<distmin<<endl;
83 #endif
84
85   uMax = Min (uMax, param+delta);
86   uMin = Max (uMin, param-delta);
87 }
88
89
90 //=======================================================================
91 //function : CurveNewton
92 //purpose  : Newton algo on curve (study S4030)
93 //=======================================================================
94
95 static Standard_Boolean CurveNewton(const Standard_Real paramPrev,
96                                     const Adaptor3d_Curve& Ad,
97                                     const gp_Pnt& P3D,
98                                     const Standard_Real /*preci*/,
99                                     Standard_Real& param,
100                                     const Standard_Real cf,
101                                     const Standard_Real cl)
102 {
103   //szv#4:S4163:12Mar99 optimized
104   Standard_Real uMin = cf, uMax = cl;
105
106   Standard_Real Tol = Precision::Confusion();
107   Standard_Real Tol2 = Tol*Tol, rs2p = 1.e+10;
108   Standard_Real X = paramPrev;
109
110   for ( Standard_Integer i=0; i<20; i++) {
111     gp_Vec v1, v2;
112     gp_Pnt pnt;
113     Ad.D2 (X, pnt, v1, v2);
114     Standard_Real nv1 = v1.SquareMagnitude();
115     if (nv1 < 1e-10) break;
116     
117     gp_Vec rs (P3D, pnt);
118     Standard_Real D = nv1 + (rs * v2);
119     if ( fabs( D ) <1e-10) break;
120     
121     Standard_Real dX = -(rs *v1)/D;
122     X += dX;
123     
124     Standard_Real rs2 = rs.SquareMagnitude();
125     if (rs2 > 4.*rs2p) break;
126     rs2p = rs2;
127     
128     if ( fabs(dX) > 1e-12) continue;
129      
130     if (X < uMin || X > uMax) break;
131     
132     Standard_Real rsn = rs * v1;
133     if (rsn*rsn/nv1 > Tol2) break;
134     param = X;
135     return Standard_True;
136   }
137
138   // cout <<"NewtonC: failed after " << i+1 << " iterations (fail " << fail << " )" << endl; 
139   return Standard_False;
140 }
141
142 //=======================================================================
143 //function : Project
144 //purpose  : 
145 //=======================================================================
146
147 Standard_Real ShapeAnalysis_Curve::Project(const Handle(Geom_Curve)& C3D,
148                                            const gp_Pnt& P3D,
149                                            const Standard_Real preci,
150                                            gp_Pnt& proj,
151                                            Standard_Real& param,
152                                            const Standard_Boolean AdjustToEnds) const
153 {
154   Standard_Real uMin = C3D->FirstParameter();
155   Standard_Real uMax = C3D->LastParameter();
156   if (uMin < uMax)  return Project (C3D,P3D,preci,proj,param,uMin,uMax,AdjustToEnds);
157   else              return Project (C3D,P3D,preci,proj,param,uMax,uMin,AdjustToEnds);
158 }
159
160 //=======================================================================
161 //function : Project
162 //purpose  : 
163 //=======================================================================
164
165 Standard_Real ShapeAnalysis_Curve::Project(const Handle(Geom_Curve)& C3D,
166                                            const gp_Pnt& P3D,
167                                            const Standard_Real preci,
168                                            gp_Pnt& proj,
169                                            Standard_Real& param,
170                                            const Standard_Real cf,
171                                            const Standard_Real cl,
172                                            const Standard_Boolean AdjustToEnds) const
173 {
174   Standard_Real distmin;
175   Standard_Real uMin = (cf < cl ? cf : cl);
176   Standard_Real uMax = (cf < cl ? cl : cf);
177   
178   if (C3D->IsKind(STANDARD_TYPE(Geom_BoundedCurve))) {
179     Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
180     gp_Pnt LowBound = C3D->Value(uMin);
181     gp_Pnt HigBound = C3D->Value(uMax);
182     distmin = LowBound.Distance(P3D);
183     if (distmin <= prec) {
184       param = uMin;
185       proj  = LowBound;
186       return distmin;
187     } 
188     distmin = HigBound.Distance(P3D);
189     if (distmin <= prec) {
190       param = uMax;
191       proj  = HigBound;
192       return distmin;
193     }
194   }
195
196   GeomAdaptor_Curve GAC(C3D, uMin, uMax);
197   if (!C3D->IsClosed()) {
198     //modified by rln on 16/12/97 after CSR# PRO11641 entity 20767
199     //the VIso was not closed (according to C3D->IsClosed()) while it "almost"
200     //was (the distance between ends of the curve was a little bit more than
201     //Precision::Confusion())
202     //in that case value 0.1 was too much and this method returned not correct parameter
203     //uMin = uMin - 0.1;
204     //uMax = uMax + 0.1;
205     // modified by pdn on 01.07.98 after BUC60195 entity 1952 (Min() added)
206     Standard_Real delta = Min (GAC.Resolution (preci), (uMax - uMin) * 0.1);
207     uMin -= delta;
208     uMax += delta;
209     GAC.Load(C3D, uMin, uMax);
210   }
211
212   return ProjectAct(GAC, P3D, preci, proj, param);
213 }
214
215 //=======================================================================
216 //function : Project
217 //purpose  : 
218 //=======================================================================
219
220 Standard_Real ShapeAnalysis_Curve::Project(const Adaptor3d_Curve& C3D,
221                                            const gp_Pnt& P3D,
222                                            const Standard_Real preci,
223                                            gp_Pnt& proj,
224                                            Standard_Real& param,
225                                            const Standard_Boolean AdjustToEnds) const
226
227 {
228   Standard_Real uMin = C3D.FirstParameter();
229   Standard_Real uMax = C3D.LastParameter();
230   Standard_Real distmin;
231   if (!Precision::IsInfinite(uMin)||!Precision::IsInfinite(uMax)) {
232     Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
233     gp_Pnt LowBound = C3D.Value(uMin);
234     gp_Pnt HigBound = C3D.Value(uMax);
235     distmin = LowBound.Distance(P3D);
236     if (distmin <= prec) {
237       param = uMin;
238       proj  = LowBound;
239       return distmin;
240     }
241     distmin = HigBound.Distance(P3D);
242     if (distmin <= prec) {
243       param = uMax;
244       proj  = HigBound;
245       return distmin;
246     } 
247   }
248   return ProjectAct(C3D, P3D, preci, proj, param);
249 }
250
251 //=======================================================================
252 //function : ProjectAct
253 //purpose  : 
254 //=======================================================================
255
256 Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
257                                               const gp_Pnt& P3D,
258                                               const Standard_Real preci,
259                                               gp_Pnt& proj,
260                                               Standard_Real& param) const
261        
262 {
263   Standard_Boolean useExtrema = Standard_True; //:i5
264
265 #ifdef WNT
266   //:i5 abv 11 Sep 98: UKI60195.igs on NT: hanging on null-length bsplines
267   if (C3D.IsClosed() && C3D.GetType()==GeomAbs_BSplineCurve) {
268     Handle(Geom_BSplineCurve) bspl = C3D.BSpline();
269     Standard_Real prec2 = gp::Resolution();
270     prec2 *= prec2;
271     gp_Pnt p0 = bspl->Pole(1), p1;
272     for ( Standard_Integer i=2; i <= bspl->NbPoles(); i++, p0 = p1 ) {
273       p1 = bspl->Pole(i);
274       if ( p0.SquareDistance ( p1 ) <= prec2 ) { 
275         useExtrema = Standard_False;
276         break;
277       }
278     }
279   }
280 #endif
281
282   Standard_Boolean OK = Standard_False;
283   if ( useExtrema ) { //:i5 //szv#4:S4163:12Mar99 try moved into if
284     try {
285       OCC_CATCH_SIGNALS
286       Extrema_ExtPC myExtPC(P3D,C3D);
287       //Standard_Boolean done = myExtPC.IsDone() && ( myExtPC.NbExt() > 0); //szv#4:S4163:12Mar99 not needed
288       if ( myExtPC.IsDone() && ( myExtPC.NbExt() > 0) ) {
289         Standard_Real dist2, dist2Min = myExtPC.SquareDistance(1);
290         Standard_Integer index = 1;
291         for ( Standard_Integer i = 2; i <= myExtPC.NbExt(); i++) {
292           dist2 = myExtPC.SquareDistance(i);
293           if ( dist2 < dist2Min) { dist2Min = dist2; index = i; }
294         }
295         param = (myExtPC.Point(index)).Parameter();
296         proj  = (myExtPC.Point(index)).Value();
297         OK = Standard_True;
298       }
299     }
300     catch(Standard_Failure) {
301       OK = Standard_False;
302 #ifdef DEB //:s5
303       cout << "\nWarning: ShapeAnalysis_Curve::ProjectAct(): Exception in Extrema_ExtPC: "; 
304       Standard_Failure::Caught()->Print(cout); cout << endl;
305 #endif
306     }
307   }
308
309   //szv#4:S4163:12Mar99 moved
310   Standard_Real uMin = C3D.FirstParameter(), uMax = C3D.LastParameter();
311   Standard_Boolean closed = Standard_False;  // si on franchit les bornes ...
312   Standard_Real distmin = RealLast(), valclosed = 0.;
313   Standard_Real aModParam = param;
314   Standard_Real aModMin = distmin;
315   
316   // PTV 29.05.2002 remember the old solution, cause it could be better
317   Standard_Real anOldParam =0.;
318   Standard_Boolean IsHaveOldSol = Standard_False;
319   gp_Pnt anOldProj;
320   if (OK) {
321     IsHaveOldSol = Standard_True;
322     anOldProj = proj;
323     anOldParam = param;
324     distmin = proj.Distance (P3D);
325     aModMin = distmin;
326     if (distmin > preci) OK = Standard_False;
327     // Cas TrimmedCurve a cheval. Voir la courbe de base.
328     // Si fermee, passer une periode
329     if (C3D.IsClosed()) {
330       closed = Standard_True;
331       valclosed = uMax - uMin; //szv#4:S4163:12Mar99 optimized
332     }
333   }
334
335   if (!OK) {
336     // BUG NICOLAS - Si le point est sur la courbe 0 Solutions
337     // Cela fonctionne avec ElCLib
338
339     // D une facon generale, on essaie de TOUJOURS retourner un resultat
340     //  MEME PAS BIEN BON. L appelant pourra decider alors quoi faire
341     param = 0.;
342
343     switch(C3D.GetType()) {
344     case GeomAbs_Circle:
345       {
346         const gp_Circ& aCirc = C3D.Circle();
347         proj = aCirc.Position().Location();
348         if(aCirc.Radius() <= gp::Resolution() ||
349            P3D.SquareDistance(proj) <= gp::Resolution() ) {
350           param = C3D.FirstParameter();
351           proj = proj.XYZ() + aCirc.XAxis().Direction().XYZ() * aCirc.Radius();
352         }
353         else {
354           param = ElCLib::Parameter(aCirc, P3D);
355           proj  = ElCLib::Value(param, aCirc);
356         }
357         closed = Standard_True;
358         valclosed = 2.*M_PI;
359       }
360       break;
361     case GeomAbs_Hyperbola:
362       {
363         param = ElCLib::Parameter(C3D.Hyperbola(), P3D);
364         proj  = ElCLib::Value(param, C3D.Hyperbola());
365       }
366       break;
367     case GeomAbs_Parabola:
368       {
369         param = ElCLib::Parameter(C3D.Parabola(), P3D);
370         proj  = ElCLib::Value(param, C3D.Parabola());
371       }
372       break;
373     case GeomAbs_Line:
374       {
375         param = ElCLib::Parameter(C3D.Line(), P3D);
376         proj  = ElCLib::Value(param, C3D.Line());
377       }
378       break;
379     case GeomAbs_Ellipse:
380       {
381         param = ElCLib::Parameter(C3D.Ellipse(), P3D);
382         proj  = ElCLib::Value(param, C3D.Ellipse());
383         closed = Standard_True;
384         valclosed = 2.*M_PI;
385
386       }
387       break;
388     default:
389       {
390         //  on ne va quand meme pas se laisser abattre ... ???
391         //  on tente ceci : 21 points sur la courbe, quel est le plus proche
392         distmin = Precision::Infinite();
393         ProjectOnSegments (C3D,P3D,25,uMin,uMax,distmin,proj,param);
394         if (distmin <= preci) return distmin;
395         if (CurveNewton(param, C3D, P3D, preci, param, uMin, uMax)) { //:S4030
396           C3D.D0(param, proj);
397           Standard_Real aDistNewton = P3D.Distance(proj);
398           if(aDistNewton < aModMin)
399             return aDistNewton;
400         }
401         // cout <<"newton failed"<<endl;    
402         ProjectOnSegments (C3D,P3D,40,uMin,uMax,distmin,proj,param);
403         if (distmin <= preci) return distmin;
404         ProjectOnSegments (C3D,P3D,20,uMin,uMax,distmin,proj,param);
405         if (distmin <= preci) return distmin;
406         ProjectOnSegments (C3D,P3D,25,uMin,uMax,distmin,proj,param);
407         if (distmin <= preci) return distmin;
408         ProjectOnSegments (C3D,P3D,40,uMin,uMax,distmin,proj,param);
409         if (distmin <= preci) return distmin;
410         //  soyons raisonnable ...
411         if(distmin > aModMin) {
412           distmin = aModMin;
413           param = aModParam;
414         }
415
416         return distmin;
417       }
418     }
419   }
420   
421   //  p = PPOC.LowerDistanceParameter();  cf try
422   if ( closed && ( param < uMin || param > uMax ) ) 
423     param += ShapeAnalysis::AdjustByPeriod ( param, 0.5 * ( uMin + uMax ), valclosed );
424   
425   if (IsHaveOldSol) {
426     // PTV 29.05.2002 Compare old solution and new;
427     Standard_Real adist1, adist2;
428     adist1 = anOldProj.SquareDistance(P3D);
429     adist2 = proj.SquareDistance (P3D);
430     if (adist1 < adist2) {
431       proj = anOldProj;
432       param = anOldParam;
433     }
434   }
435   return proj.Distance (P3D);
436 }
437
438
439 //=======================================================================
440 //function : NextProject
441 //purpose  : Newton algo for projecting point on curve (S4030)
442 //=======================================================================
443
444 Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
445                                                const Handle(Geom_Curve)& C3D,
446                                                const gp_Pnt& P3D,
447                                                const Standard_Real preci,
448                                                gp_Pnt& proj,
449                                                Standard_Real& param,
450                                                const Standard_Real cf,
451                                                const Standard_Real cl,
452                                                const Standard_Boolean AdjustToEnds) const
453 {
454   Standard_Real uMin = (cf < cl ? cf : cl);
455   Standard_Real uMax = (cf < cl ? cl : cf);
456   Standard_Real distmin;
457   if (C3D->IsKind(STANDARD_TYPE(Geom_BoundedCurve))) {
458     Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
459     gp_Pnt LowBound = C3D->Value(uMin);
460     gp_Pnt HigBound = C3D->Value(uMax);
461     distmin = LowBound.Distance(P3D);
462     if (distmin <= prec) {
463       param = uMin;
464       proj  = LowBound;
465       return distmin;
466     } 
467     distmin = HigBound.Distance(P3D);
468     if (distmin <= prec) {
469       param = uMax;
470       proj  = HigBound;
471       return distmin;
472     }
473   }
474
475   GeomAdaptor_Curve GAC(C3D, uMin, uMax);
476   if (!C3D->IsClosed()) {
477     //modified by rln on 16/12/97 after CSR# PRO11641 entity 20767
478     //the VIso was not closed (according to C3D->IsClosed()) while it "almost"
479     //was (the distance between ends of the curve was a little bit more than
480     //Precision::Confusion())
481     //in that case value 0.1 was too much and this method returned not correct parameter
482     //uMin = uMin - 0.1;
483     //uMax = uMax + 0.1;
484     // modified by pdn on 01.07.98 after BUC60195 entity 1952 (Min() added)
485     Standard_Real delta = Min (GAC.Resolution (preci), (uMax - uMin) * 0.1);
486     uMin -= delta;
487     uMax += delta;
488     GAC.Load(C3D, uMin, uMax);
489   }
490   return NextProject ( paramPrev, GAC, P3D, preci, proj, param );
491 }
492
493 //=======================================================================
494 //function : NextProject
495 //purpose  : 
496 //=======================================================================
497
498 Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
499                                                const Adaptor3d_Curve& C3D,
500                                                const gp_Pnt& P3D,
501                                                const Standard_Real preci,
502                                                gp_Pnt& proj,
503                                                Standard_Real& param) const
504 {
505   Standard_Real uMin = C3D.FirstParameter();
506   Standard_Real uMax = C3D.LastParameter();
507   
508   if (CurveNewton(paramPrev, C3D, P3D, preci, param, uMin, uMax)) { 
509     C3D.D0(param, proj);
510     return P3D.Distance(proj);
511   }
512
513   return Project(C3D, P3D, preci, proj, param);
514 }
515
516 //=======================================================================
517 //function : AdjustParameters
518 //purpose  : Copied from StepToTopoDS_GeometricTuul::UpdateParam3d (Aug 2001)
519 //=======================================================================
520
521 Standard_Boolean ShapeAnalysis_Curve::ValidateRange (const Handle(Geom_Curve)& theCurve, 
522                                                      Standard_Real& First, Standard_Real& Last,
523                                                      const Standard_Real preci) const
524 {
525   // First et/ou Last peuvent etre en dehors des bornes naturelles de la courbe.
526   // On donnera alors la valeur en bout a First et/ou Last
527   
528   Standard_Real cf = theCurve->FirstParameter();
529   Standard_Real cl = theCurve->LastParameter();
530 //  Standard_Real preci = BRepAPI::Precision();
531
532   if (theCurve->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) && !theCurve->IsClosed()) {
533     if (First < cf) {
534 #ifdef DEBUG
535       cout << "Update Edge First Parameter to Curve First Parameter" << endl;
536 #endif
537       First = cf;
538     }
539     else if (First > cl) {
540 #ifdef DEBUG
541       cout << "Update Edge First Parameter to Curve Last Parameter" << endl;
542 #endif
543       First = cl;
544     }
545     if (Last < cf) {
546 #ifdef DEBUG
547       cout << "Update Edge Last Parameter to Curve First Parameter" << endl;
548 #endif
549       Last = cf;
550     }
551     else if (Last > cl) {
552 #ifdef DEBUG
553       cout << "Update Edge Last Parameter to Curve Last Parameter" << endl;
554 #endif
555       Last = cl;
556     }
557   }
558
559   if (First < Last) return Standard_True;
560
561   // 15.11.2002 PTV OCC966
562   if (ShapeAnalysis_Curve::IsPeriodic(theCurve)) 
563     ElCLib::AdjustPeriodic(cf,cl,Precision::PConfusion(),First,Last); //:a7 abv 11 Feb 98: preci -> PConfusion()
564   else if (theCurve->IsClosed()) {
565     // l'un des points projecte se trouve sur l'origine du parametrage
566     // de la courbe 3D. L algo a donne cl +- preci au lieu de cf ou vice-versa
567     // DANGER precision 3d applique a une espace 1d
568     
569     // Last = cf au lieu de Last = cl
570     if      (Abs(Last - cf) < Precision::PConfusion() /*preci*/)  Last = cl ;
571     // First = cl au lieu de First = cf
572     else if (Abs(First - cl) < Precision::PConfusion() /*preci*/)  First = cf;
573
574     // on se trouve dans un cas ou l origine est traversee
575     // illegal sur une courbe fermee non periodique
576     // on inverse quand meme les parametres !!!!!!
577     else {
578       //:S4136 abv 20 Apr 99: r0701_ug.stp #6230: add check in 3d
579       if ( theCurve->Value(First).Distance(theCurve->Value(cf)) < preci ) First = cf;
580       if ( theCurve->Value(Last).Distance(theCurve->Value(cl)) < preci ) Last = cl;
581       if ( First > Last ) {
582 #ifdef DEBUG
583         cout << "Warning : parameter range of edge crossing non periodic curve origin" << endl;
584 #endif
585         Standard_Real tmp = First;
586         First = Last;
587         Last = tmp;
588       }
589     }
590   }
591   // The curve is closed within the 3D tolerance
592   else if (theCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
593     Handle(Geom_BSplineCurve) aBSpline = 
594       Handle(Geom_BSplineCurve)::DownCast(theCurve);
595     if (aBSpline->StartPoint().Distance(aBSpline->EndPoint()) <= preci ) {
596 //:S4136        <= BRepAPI::Precision()) {
597       // l'un des points projecte se trouve sur l'origine du parametrage
598       // de la courbe 3D. L algo a donne cl +- preci au lieu de cf ou vice-versa
599       // DANGER precision 3d applique a une espace 1d
600       
601       // Last = cf au lieu de Last = cl
602       if      (Abs(Last - cf) < Precision::PConfusion() /*preci*/) Last = cl ;
603       // First = cl au lieu de First = cf
604       else if (Abs(First - cl) < Precision::PConfusion() /*preci*/) First =  cf;
605
606       // on se trouve dans un cas ou l origine est traversee
607       // illegal sur une courbe fermee non periodique
608       // on inverse quand meme les parametres !!!!!!
609       else {
610 #ifdef DEBUG
611         cout << "Warning : parameter range of edge crossing non periodic curve origin" << endl;
612 #endif
613         Standard_Real tmp = First;
614         First = Last;
615         Last = tmp;
616       }
617     }
618     //abv 15.03.00 #72 bm1_pe_t4 protection of exceptions in draw
619     else if ( First > Last ) {
620 #ifdef DEBUG
621       cout << "Warning: parameter range is bad; curve reversed" << endl;
622 #endif
623       First = theCurve->ReversedParameter ( First );
624       Last = theCurve->ReversedParameter ( Last );
625       theCurve->Reverse();
626     }
627 //:j9 abv 11 Dec 98: PRO7747 #4875, after :j8:    else 
628     if (First == Last) {  //gka 10.07.1998 file PRO7656 entity 33334
629       First = cf; Last = cl;
630       return Standard_False;
631     }
632   }
633   else {
634 #ifdef DEBUG
635     cout << "UpdateParam3d Failed" << endl;
636     cout << "  - Curve Type : " << theCurve->DynamicType() << endl;
637     cout << "  - Param 1    : " << First << endl;
638     cout << "  - Param 2    : " << Last << endl;
639 #endif
640     //abv 15.03.00 #72 bm1_pe_t4 protection of exceptions in draw
641     if ( First > Last ) {
642 #ifdef DEBUG
643       cout << "Warning: parameter range is bad; curve reversed" << endl;
644 #endif
645       First = theCurve->ReversedParameter ( First );
646       Last = theCurve->ReversedParameter ( Last );
647       theCurve->Reverse();
648     }
649     //pdn 11.01.99 #144 bm1_pe_t4 protection of exceptions in draw
650     if (First == Last) {
651       First -= Precision::PConfusion();
652       Last += Precision::PConfusion();
653     }
654     return Standard_False;
655   }
656   return Standard_True;
657 }
658
659 //=======================================================================
660 //function : FillBndBox
661 //purpose  : WORK-AROUND for methods brom BndLib which give not exact bounds
662 //=======================================================================
663
664 // search for extremum using Newton
665 static Standard_Integer SearchForExtremum (const Handle(Geom2d_Curve)& C2d,
666                                            const Standard_Real First,
667                                            const Standard_Real Last,
668                                            const gp_Vec2d &dir,
669                                            Standard_Real &par,
670                                            gp_Pnt2d &res)
671 {
672   Standard_Real prevpar;
673   for ( Standard_Integer i=0; i <10; i++ ) {
674     prevpar = par;
675       
676     gp_Vec2d D1, D2;
677     C2d->D2 ( par, res, D1, D2 );
678     Standard_Real Det = ( D2 * dir );
679     if ( Abs ( Det ) < 1e-10 ) return Standard_True;
680     
681     par -= ( D1 * dir ) / Det;
682     if ( Abs ( par - prevpar ) < Precision::PConfusion() ) return Standard_True;
683     
684     if ( First - par >= Precision::PConfusion() || 
685          par - Last  >= Precision::PConfusion() ) return Standard_False;
686   }
687   return Standard_True;
688 }
689
690 void ShapeAnalysis_Curve::FillBndBox (const Handle(Geom2d_Curve)& C2d,
691                                       const Standard_Real First,
692                                       const Standard_Real Last,
693                                       const Standard_Integer NPoints,
694                                       const Standard_Boolean Exact,
695                                       Bnd_Box2d &Box) const
696 {
697   Standard_Integer nseg = ( NPoints <2 ? 1 : NPoints-1 );
698   Standard_Real step = ( Last - First ) / nseg;
699   for ( Standard_Integer i=0; i <= nseg; i++ ) {
700     Standard_Real par = First + i * step;
701     gp_Pnt2d pnt = C2d->Value ( par );
702     Box.Add ( pnt );
703     if ( ! Exact ) continue;
704     
705     gp_Pnt2d pextr;
706     Standard_Real parextr = par;
707     if ( SearchForExtremum ( C2d, Max(First,par-2.*step), Min(Last,par+2.*step),
708                              gp_Vec2d(1,0), parextr, pextr ) ) {
709       Box.Add ( pextr );
710     }
711     parextr = par;
712     if ( SearchForExtremum ( C2d, Max(First,par-2.*step), Min(Last,par+2.*step),
713                              gp_Vec2d(0,1), parextr, pextr ) ) {
714       Box.Add ( pextr );
715     }
716   }
717 }
718
719 //=======================================================================
720 //function : SelectForwardSeam
721 //purpose  : 
722 //=======================================================================
723
724 Standard_Integer ShapeAnalysis_Curve::SelectForwardSeam(const Handle(Geom2d_Curve)& C1,
725                                                         const Handle(Geom2d_Curve)& C2) const
726 {
727   //  SelectForward est destine a devenir un outil distinct
728   //  Il est sans doute optimisable !
729
730   Standard_Integer theCurveIndice = 0;
731
732   Handle(Geom2d_Line) L1 = Handle(Geom2d_Line)::DownCast(C1);
733   if (L1.IsNull()) {
734     // if we have BoundedCurve, create a line from C1
735     Handle(Geom2d_BoundedCurve) BC1 = Handle(Geom2d_BoundedCurve)::DownCast(C1);
736     if (BC1.IsNull()) return theCurveIndice;
737     gp_Pnt2d StartBC1 = BC1->StartPoint();
738     gp_Pnt2d EndBC1   = BC1->EndPoint();
739     gp_Vec2d VecBC1(StartBC1, EndBC1);
740     L1 = new Geom2d_Line(StartBC1, VecBC1);
741   }
742
743   Handle(Geom2d_Line) L2 = Handle(Geom2d_Line)::DownCast(C2);
744   if (L2.IsNull()) {
745     // if we have BoundedCurve, creates a line from C2
746     Handle(Geom2d_BoundedCurve) BC2 = Handle(Geom2d_BoundedCurve)::DownCast(C2);
747     if (BC2.IsNull()) return theCurveIndice;
748     gp_Pnt2d StartBC2 = BC2->StartPoint();
749     gp_Pnt2d EndBC2   = BC2->EndPoint();
750     gp_Vec2d VecBC2(StartBC2, EndBC2);
751     L2 = new Geom2d_Line(StartBC2, VecBC2);
752   }
753
754   Standard_Boolean UdirPos, UdirNeg, VdirPos, VdirNeg;
755   UdirPos = UdirNeg = VdirPos = VdirNeg = Standard_False;
756
757   gp_Dir2d theDir  = L1->Direction();
758   gp_Pnt2d theLoc1 = L1->Location();
759   gp_Pnt2d theLoc2 = L2->Location();
760
761   if        (theDir.X() > 0.) {
762     UdirPos = Standard_True; //szv#4:S4163:12Mar99 Udir unused
763   } else if (theDir.X() < 0.) {
764     UdirNeg = Standard_True; //szv#4:S4163:12Mar99 Udir unused
765   } else if (theDir.Y() > 0.) {
766     VdirPos = Standard_True; //szv#4:S4163:12Mar99 Vdir unused
767   } else if (theDir.Y() < 0.) {
768     VdirNeg = Standard_True; //szv#4:S4163:12Mar99 Vdir unused
769   }
770   
771   if        (VdirPos) {
772     // max of Loc1.X() Loc2.X()
773     if (theLoc1.X() > theLoc2.X())   theCurveIndice  = 1;
774     else                             theCurveIndice  = 2;
775   } else if (VdirNeg) {
776     if (theLoc1.X() > theLoc2.X())   theCurveIndice  = 2;
777     else                             theCurveIndice  = 1;
778   } else if (UdirPos) {
779     // min of Loc1.X() Loc2.X()
780     if (theLoc1.Y() < theLoc2.Y())   theCurveIndice  = 1;
781     else                             theCurveIndice  = 2;
782   } else if (UdirNeg) {
783     if (theLoc1.Y() < theLoc2.Y())   theCurveIndice  = 2;
784     else                             theCurveIndice  = 1;
785   }
786
787   return theCurveIndice;
788 }
789
790 //%11 pdn 12.01.98
791 //=============================================================================
792 // Static methods for IsPlanar
793 // IsPlanar
794 //=============================================================================
795
796 static gp_XYZ GetAnyNormal ( gp_XYZ orig )
797 {
798   gp_XYZ Norm;
799   if ( Abs ( orig.Z() ) < Precision::Confusion() )
800     Norm.SetCoord ( 0, 0, 1 );
801   else {
802     Norm.SetCoord ( orig.Z(), 0, -orig.X() );
803     Standard_Real nrm = Norm.Modulus();
804     if ( nrm < Precision::Confusion() ) Norm.SetCoord ( 0, 0, 1 );
805     else Norm = Norm / nrm;
806   }
807   return Norm;
808 }
809
810 //=======================================================================
811 //function : GetSamplePoints
812 //purpose  : 
813 //=======================================================================
814 static void AppendControlPoles (TColgp_SequenceOfPnt& seq,
815                                 const Handle(Geom_Curve) curve)
816 {
817   if ( curve->IsKind(STANDARD_TYPE(Geom_Line))) {
818     seq.Append(curve->Value(0));
819     seq.Append(curve->Value(1));
820   } else if ( curve->IsKind(STANDARD_TYPE(Geom_Conic))) {
821     seq.Append(curve->Value(0));
822     seq.Append(curve->Value(M_PI/2));
823     seq.Append(curve->Value(M_PI));
824   } else if ( curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
825     //DeclareAndCast(Geom_TrimmedCurve, Trimmed, curve);
826     Handle(Geom_TrimmedCurve) Trimmed = *((Handle(Geom_TrimmedCurve) *) &curve);
827 //     AppendControlPoles(seq,Trimmed->BasisCurve());
828     Handle(Geom_Curve) aBaseCrv = Trimmed->BasisCurve();
829     Standard_Boolean done = Standard_False;
830     if ( aBaseCrv->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
831       try {
832         OCC_CATCH_SIGNALS
833         Handle(Geom_Geometry) Ctmp = aBaseCrv->Copy();
834         Handle(Geom_BSplineCurve) bslp = Handle(Geom_BSplineCurve)::DownCast(Ctmp);
835         bslp->Segment(curve->FirstParameter(), curve->LastParameter());
836         AppendControlPoles(seq,bslp);
837         done = Standard_True;
838       }
839       catch (Standard_Failure) {
840       }
841     }
842     else if ( aBaseCrv->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
843       try {
844         OCC_CATCH_SIGNALS
845         Handle(Geom_Geometry) Ctmp = aBaseCrv->Copy();
846         Handle(Geom_BezierCurve) bz = Handle(Geom_BezierCurve)::DownCast(Ctmp);
847         bz->Segment(curve->FirstParameter(), curve->LastParameter());
848         AppendControlPoles(seq,bz);
849         done = Standard_True;
850       }
851       catch (Standard_Failure) {
852       }
853     }
854     if (!done) {
855       seq.Append(curve->Value(curve->FirstParameter()));
856       seq.Append(curve->Value((curve->FirstParameter() + curve->LastParameter())/2.));
857       seq.Append(curve->Value(curve->LastParameter()));
858     }
859   } else if ( curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
860     //DeclareAndCast(Geom_OffsetCurve, OffsetC, curve);
861     Handle(Geom_OffsetCurve) OffsetC = *((Handle(Geom_OffsetCurve) *) &curve);
862 //     AppendControlPoles(seq,OffsetC->BasisCurve());
863     seq.Append(curve->Value(curve->FirstParameter()));
864     seq.Append(curve->Value((curve->FirstParameter() + curve->LastParameter())/2.));
865     seq.Append(curve->Value(curve->LastParameter()));
866   } else if ( curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
867     //DeclareAndCast(Geom_BSplineCurve, BSpline, curve);
868     Handle(Geom_BSplineCurve) BSpline = *((Handle(Geom_BSplineCurve) *) &curve);
869     TColgp_Array1OfPnt Poles(1,BSpline->NbPoles());
870     BSpline->Poles(Poles);
871     for(Standard_Integer i = 1; i <= BSpline->NbPoles(); i++)
872       seq.Append(Poles(i));
873   } else if ( curve->IsKind(STANDARD_TYPE(Geom_BezierCurve)))  {
874     //DeclareAndCast(Geom_BezierCurve, Bezier, curve);
875     //Handle(Geom_BezierCurve) Bezier = Handle(Geom_BezierCurve)::DownCast(curve);
876     Handle(Geom_BezierCurve) Bezier = *((Handle(Geom_BezierCurve) *) &curve);
877     TColgp_Array1OfPnt Poles(1,Bezier->NbPoles());
878     Bezier->Poles(Poles);
879     for(Standard_Integer i = 1; i <= Bezier->NbPoles(); i++)
880       seq.Append(Poles(i));
881   }
882 }
883
884 //%11 pdn 12.01.98
885 //szv modified
886 //=======================================================================
887 //function : IsPlanar
888 //purpose  : Detects if points lie in some plane and returns normal
889 //=======================================================================
890
891 Standard_Boolean ShapeAnalysis_Curve::IsPlanar (const TColgp_Array1OfPnt& pnts,
892                                                 gp_XYZ& Normal,
893                                                 const Standard_Real preci)
894 {
895   Standard_Real precision = (preci > 0.0)? preci : Precision::Confusion();
896   Standard_Boolean noNorm = (Normal.SquareModulus() == 0);
897
898   if (pnts.Length() < 3) {
899     gp_XYZ N1 = pnts(1).XYZ()-pnts(2).XYZ();
900     if (noNorm) {
901       Normal = GetAnyNormal(N1);
902       return Standard_True;
903     }
904     return Abs(N1*Normal) < Precision::Confusion();
905   }
906   
907   gp_XYZ aMaxDir;
908   if (noNorm) {
909     //define a center point
910     gp_XYZ aCenter(0,0,0);
911     Standard_Integer i = 1;
912     for (; i <= pnts.Length(); i++) 
913       aCenter += pnts(i).XYZ();
914     aCenter/=pnts.Length();
915     
916     
917     aMaxDir = pnts(1).XYZ() - aCenter;
918     Normal = (pnts(pnts.Length()).XYZ() - aCenter) ^ aMaxDir;
919
920     for ( i = 1; i < pnts.Length(); i++) {
921       gp_XYZ aTmpDir = pnts(i+1).XYZ() - aCenter;
922       if(aTmpDir.SquareModulus() > aMaxDir.SquareModulus())
923         aMaxDir = aTmpDir;
924
925       gp_XYZ aDelta = (pnts(i).XYZ() - aCenter) ^ (pnts(i+1).XYZ() - aCenter);
926       if(Normal*aDelta < 0)
927         aDelta*=-1;
928       Normal += aDelta;
929     }
930   }
931
932   // check if points are linear
933   Standard_Real nrm = Normal.Modulus();
934   if ( nrm < Precision::Confusion() ) {
935     Normal = GetAnyNormal(aMaxDir);
936     return Standard_True;
937   }
938   Normal = Normal / nrm;
939
940   Standard_Real mind = RealLast(), maxd = -RealLast(), dev;
941   for (Standard_Integer i = 1; i <= pnts.Length(); i++) {
942     dev = pnts(i).XYZ() * Normal;
943     if (dev < mind) mind = dev;
944     if (dev > maxd) maxd = dev;
945   }
946
947   return ((maxd - mind) <= precision);
948 }
949
950
951 //=======================================================================
952 //function : IsPlanar
953 //purpose  : 
954 //=======================================================================
955
956  Standard_Boolean ShapeAnalysis_Curve::IsPlanar (const Handle(Geom_Curve)& curve,
957                                                  gp_XYZ& Normal,
958                                                  const Standard_Real preci)
959 {
960   Standard_Real precision = (preci > 0.0)? preci : Precision::Confusion();
961   Standard_Boolean noNorm = (Normal.SquareModulus() == 0);
962
963   if (curve->IsKind(STANDARD_TYPE(Geom_Line))) {
964     //DeclareAndCast(Geom_Line, Line, curve);
965     Handle(Geom_Line) Line = *((Handle(Geom_Line) *) &curve);
966     gp_XYZ N1 = Line->Position().Direction().XYZ();
967     if (noNorm) {
968       Normal = GetAnyNormal(N1);
969       return Standard_True;
970     }
971     return Abs(N1*Normal) < Precision::Confusion();
972   }
973
974   if (curve->IsKind(STANDARD_TYPE(Geom_Conic))) {
975     //DeclareAndCast(Geom_Conic, Conic, curve);
976     Handle(Geom_Conic) Conic = *((Handle(Geom_Conic) *) &curve);
977     gp_XYZ N1 = Conic->Axis().Direction().XYZ();
978     if (noNorm) {
979       Normal = N1;
980       return Standard_True;
981     }
982     gp_XYZ aVecMul = N1^Normal;
983     return aVecMul.SquareModulus() < Precision::Confusion()*Precision::Confusion();
984   }
985
986   if (curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
987     //DeclareAndCast(Geom_TrimmedCurve, Trimmed, curve);
988     Handle(Geom_TrimmedCurve) Trimmed = *((Handle(Geom_TrimmedCurve) *) &curve);
989     return IsPlanar(Trimmed->BasisCurve(),Normal,precision);
990   }
991
992   if (curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
993     //DeclareAndCast(Geom_OffsetCurve, OffsetC, curve);
994     Handle(Geom_OffsetCurve) OffsetC = *((Handle(Geom_OffsetCurve) *) &curve);
995     return IsPlanar(OffsetC->BasisCurve(),Normal,precision);
996   }
997
998   if (curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
999     //DeclareAndCast(Geom_BSplineCurve, BSpline, curve);
1000     Handle(Geom_BSplineCurve) BSpline = *((Handle(Geom_BSplineCurve) *) &curve);
1001     TColgp_Array1OfPnt Poles(1,BSpline->NbPoles());
1002     BSpline->Poles(Poles);
1003     return IsPlanar(Poles,Normal,precision);
1004   }
1005
1006   if (curve->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
1007     //DeclareAndCast(Geom_BezierCurve, Bezier, curve);
1008     Handle(Geom_BezierCurve) Bezier = *((Handle(Geom_BezierCurve) *) &curve);
1009     TColgp_Array1OfPnt Poles(1,Bezier->NbPoles());
1010     Bezier->Poles(Poles);
1011     return IsPlanar(Poles,Normal,precision);
1012   }
1013
1014   if (curve->IsKind(STANDARD_TYPE(ShapeExtend_ComplexCurve))) {
1015     //DeclareAndCast(ShapeExtend_ComplexCurve, Complex, curve);
1016     Handle(ShapeExtend_ComplexCurve) Complex = *((Handle(ShapeExtend_ComplexCurve) *) &curve);
1017     TColgp_SequenceOfPnt sequence;
1018     Standard_Integer i; // svv Jan11 2000 : porting on DEC
1019     for (i = 1; i <= Complex->NbCurves(); i++)
1020       AppendControlPoles(sequence,Complex->Curve(i));
1021     TColgp_Array1OfPnt Poles(1,sequence.Length());
1022     for (i=1; i <= sequence.Length(); i++) Poles(i) = sequence(i);
1023     return IsPlanar(Poles,Normal,precision);
1024   }
1025
1026   return Standard_False;
1027 }
1028
1029 //=======================================================================
1030 //function : GetSamplePoints
1031 //purpose  : 
1032 //=======================================================================
1033
1034 Standard_Boolean ShapeAnalysis_Curve::GetSamplePoints (const Handle(Geom_Curve)& curve,
1035                                                        const Standard_Real first,
1036                                                        const Standard_Real last,
1037                                                        TColgp_SequenceOfPnt& seq)
1038 {
1039   Standard_Real adelta = curve->LastParameter() - curve->FirstParameter();
1040   if(!adelta )
1041     return Standard_False;
1042   
1043   Standard_Integer aK = (Standard_Integer)ceil ((last - first) / adelta);
1044   Standard_Integer nbp =100*aK;
1045   if(curve->IsKind(STANDARD_TYPE(Geom_Line)))
1046     nbp =2;
1047   else if(curve->IsKind(STANDARD_TYPE(Geom_Circle)))
1048     nbp =360*aK;
1049   
1050   else if (curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
1051     Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(curve);
1052     
1053     nbp = aBspl->NbKnots() * aBspl->Degree()*aK;
1054     if(nbp < 2.0) nbp=2;
1055   }
1056   else if (curve->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
1057     Handle(Geom_BezierCurve) aB = Handle(Geom_BezierCurve)::DownCast(curve);
1058     nbp = 3 + aB->NbPoles();
1059   }
1060   else if(curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
1061     Handle(Geom_OffsetCurve) aC = Handle(Geom_OffsetCurve)::DownCast(curve);
1062     return GetSamplePoints(aC->BasisCurve(),first,last,seq);
1063   }
1064   else if(curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
1065     Handle(Geom_TrimmedCurve) aC = Handle(Geom_TrimmedCurve)::DownCast(curve);
1066     return GetSamplePoints(aC->BasisCurve(),first,last,seq);
1067   }
1068   Standard_Real step = ( last - first ) / (Standard_Real)( nbp - 1 );
1069   Standard_Real par = first, stop = last - 0.5 * step;
1070   for ( ; par < stop; par += step )
1071     seq.Append(curve->Value(par));
1072   seq.Append(curve->Value(last));
1073   return Standard_True;
1074 }
1075 //=======================================================================
1076 //function : GetSamplePoints
1077 //purpose  : 
1078 //=======================================================================
1079
1080 Standard_Boolean ShapeAnalysis_Curve::GetSamplePoints (const Handle(Geom2d_Curve)& curve,
1081                                                        const Standard_Real first,
1082                                                        const Standard_Real last,
1083                                                        TColgp_SequenceOfPnt2d& seq)
1084 {
1085   //:abv 05.06.02: TUBE.stp 
1086   // Use the same distribution of points as BRepTopAdaptor_FClass2d for consistency
1087   Geom2dAdaptor_Curve C ( curve, first, last );
1088   Standard_Integer nbs = Geom2dInt_Geom2dCurveTool::NbSamples(C);
1089   //-- Attention aux bsplines rationnelles de degree 3. (bouts de cercles entre autres)
1090   if (nbs > 2) nbs*=4;
1091   Standard_Real step = ( last - first ) / (Standard_Real)( nbs - 1 );
1092   Standard_Real par = first, stop = last - 0.5 * step;
1093   for ( ; par < stop; par += step )
1094     seq.Append(curve->Value(par));
1095   seq.Append(curve->Value(last));
1096   return Standard_True;
1097 /*
1098   Standard_Integer i;
1099   Standard_Real step;
1100   gp_Pnt2d Ptmp;
1101   if ( curve->IsKind(STANDARD_TYPE(Geom2d_Line))) {
1102     seq.Append(curve->Value(first));
1103     seq.Append(curve->Value(last));
1104     return Standard_True;
1105   }
1106   else if(curve->IsKind(STANDARD_TYPE(Geom2d_Conic))) {
1107     step = Min ( M_PI, last-first ) / 19; //:abv 05.06.02 TUBE.stp #19209...: M_PI/16
1108 //     if( step>(last-first) ) {
1109 //       seq.Append(curve->Value(first));
1110 //       seq.Append(curve->Value((last+first)/2));
1111 //       seq.Append(curve->Value(last));
1112 //       return Standard_True;
1113 //     }
1114 //     else {
1115       Standard_Real par=first;
1116       for(i=0; par<last; i++) {
1117         seq.Append(curve->Value(par));
1118         par += step;
1119       }
1120       seq.Append(curve->Value(last));
1121       return Standard_True;
1122 //     }
1123   }
1124   else if ( curve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) {
1125     DeclareAndCast(Geom2d_TrimmedCurve, Trimmed, curve);
1126     return GetControlPoints(Trimmed->BasisCurve(), seq, first, last);
1127   }
1128   else if ( curve->IsKind(STANDARD_TYPE(Geom2d_BSplineCurve))) {
1129     DeclareAndCast(Geom2d_BSplineCurve, aBs, curve);
1130     TColStd_SequenceOfReal aSeqParam;
1131     if(!aBs.IsNull()) {
1132       aSeqParam.Append(first);
1133       for(i=1; i<=aBs->NbKnots(); i++) {
1134         if( aBs->Knot(i)>first && aBs->Knot(i)<last ) 
1135           aSeqParam.Append(aBs->Knot(i));
1136       }
1137       aSeqParam.Append(last);
1138       Standard_Integer NbPoints=aBs->Degree();
1139       if( (aSeqParam.Length()-1)*NbPoints>10 ) {
1140         for(i=1; i<aSeqParam.Length(); i++) {
1141           Standard_Real FirstPar = aSeqParam.Value(i);
1142           Standard_Real LastPar = aSeqParam.Value(i+1);
1143           step = (LastPar-FirstPar)/NbPoints;
1144           for(Standard_Integer k=0; k<NbPoints; k++ ) {
1145             aBs->D0(FirstPar+step*k, Ptmp);
1146             seq.Append(Ptmp);
1147           }
1148         }
1149         aBs->D0(last, Ptmp);
1150         seq.Append(Ptmp);
1151         return Standard_True;
1152       }
1153       else {
1154         step = (last-first)/10;
1155         for(Standard_Integer k=0; k<=10; k++ ) {
1156           aBs->D0(first+step*k, Ptmp);
1157           seq.Append(Ptmp);
1158         }
1159         return Standard_True;
1160       }
1161     }
1162   }
1163   else if ( curve->IsKind(STANDARD_TYPE(Geom2d_BezierCurve)))  {
1164     DeclareAndCast(Geom2d_BezierCurve, aBz, curve);
1165     if(!aBz.IsNull()) {
1166       Standard_Integer NbPoints=aBz->Degree();
1167       step = (last-first)/NbPoints;
1168       for(Standard_Integer k=0; k<NbPoints; k++ ) {
1169         aBz->D0(first+step*k, Ptmp);
1170         seq.Append(Ptmp);
1171       }
1172       aBz->D0(last, Ptmp);
1173       seq.Append(Ptmp);
1174       return Standard_True;
1175     }
1176   }
1177   return Standard_False;
1178 */
1179 }
1180
1181 //=======================================================================
1182 //function : IsClosed
1183 //purpose  : 
1184 //=======================================================================
1185
1186 Standard_Boolean ShapeAnalysis_Curve::IsClosed(const Handle(Geom_Curve)& theCurve,
1187                                                const Standard_Real preci)
1188 {
1189   if (theCurve->IsClosed())
1190     return Standard_True;
1191   
1192   Standard_Real prec = Max (preci, Precision::Confusion());
1193
1194   Standard_Real f, l;
1195   f = theCurve->FirstParameter();
1196   l = theCurve->LastParameter();
1197   
1198   if (Precision::IsInfinite (f) || Precision::IsInfinite (l))
1199     return Standard_False;
1200   
1201   Standard_Real aClosedVal = theCurve->Value(f).SquareDistance(theCurve->Value(l));
1202   Standard_Real preci2 = prec*prec;
1203
1204   return (aClosedVal <= preci2);
1205 }
1206 //=======================================================================
1207 //function : IsPeriodic
1208 //purpose  : OCC996
1209 //=======================================================================
1210
1211 Standard_Boolean ShapeAnalysis_Curve::IsPeriodic(const Handle(Geom_Curve)& theCurve)
1212 {
1213   // 15.11.2002 PTV OCC966
1214   // remove regressions in DE tests (diva, divb, divc, toe3) in KAS:dev
1215   // ask IsPeriodic on BasisCurve
1216   Handle(Geom_Curve) aTmpCurve = theCurve;
1217   while ( (aTmpCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) ||
1218           (aTmpCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) ) {
1219     if (aTmpCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
1220       aTmpCurve = Handle(Geom_OffsetCurve)::DownCast(aTmpCurve)->BasisCurve();
1221     if (aTmpCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
1222       aTmpCurve = Handle(Geom_TrimmedCurve)::DownCast(aTmpCurve)->BasisCurve();
1223   }
1224   return aTmpCurve->IsPeriodic();
1225 }
1226
1227 Standard_Boolean ShapeAnalysis_Curve::IsPeriodic(const Handle(Geom2d_Curve)& theCurve)
1228 {
1229   // 15.11.2002 PTV OCC966
1230   // remove regressions in DE tests (diva, divb, divc, toe3) in KAS:dev
1231   // ask IsPeriodic on BasisCurve
1232   Handle(Geom2d_Curve) aTmpCurve = theCurve;
1233   while ( (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_OffsetCurve))) ||
1234           (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) ) {
1235     if (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_OffsetCurve)))
1236       aTmpCurve = Handle(Geom2d_OffsetCurve)::DownCast(aTmpCurve)->BasisCurve();
1237     if (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
1238       aTmpCurve = Handle(Geom2d_TrimmedCurve)::DownCast(aTmpCurve)->BasisCurve();
1239   }
1240   return aTmpCurve->IsPeriodic();
1241 }