src/ShapeAnalysis/ShapeAnalysis_Curve.cxx

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