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