src/ShapeAnalysis/ShapeAnalysis_Curve.cxx

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