1 // Copyright (c) 1999-2012 OPEN CASCADE SAS
3 // The content of this file is subject to the Open CASCADE Technology Public
4 // License Version 6.5 (the "License"). You may not use the content of this file
5 // except in compliance with the License. Please obtain a copy of the License
6 // at http://www.opencascade.org and read it completely before using this file.
8 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
9 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
11 // The Original Code and all software distributed under the License is
12 // distributed on an "AS IS" basis, without warranty of any kind, and the
13 // Initial Developer hereby disclaims all such warranties, including without
14 // limitation, any warranties of merchantability, fitness for a particular
15 // purpose or non-infringement. Please see the License for the specific terms
16 // and conditions governing the rights and limitations under the License.
18 // pdn 04.12.98 Add method using Adaptor_Curve
19 //:j8 abv 10.12.98 TR10 r0501_db.stp #9423
20 //pdn 25.12.98 private method ProjectAct
22 //:s5 abv 22.04.99 Adding debug printouts in catch {} blocks
23 // abv 14.05.99 S4174: Adding method for exact computing of the boundary box
24 // gka 21.06.99 S4208: adding method NextProject(Adaptor_Curve)
25 // msv 30.05.00 correct IsPlanar for a conic curve
26 #include <ShapeAnalysis_Curve.ixx>
30 #include <Geom2d_BoundedCurve.hxx>
31 #include <Geom2d_Line.hxx>
32 #include <Geom_BSplineCurve.hxx>
33 #include <GeomAdaptor_Curve.hxx>
35 #include <Precision.hxx>
37 #include <Standard_ErrorHandler.hxx>
38 #include <Standard_Failure.hxx>
39 #include <Adaptor3d_Curve.hxx>
40 #include <Extrema_ExtPC.hxx>
41 #include <ShapeAnalysis.hxx>
42 #include <TColgp_SequenceOfPnt.hxx>
43 #include <Geom_Line.hxx>
44 #include <Geom_Conic.hxx>
45 #include <Geom_TrimmedCurve.hxx>
46 #include <Geom_OffsetCurve.hxx>
47 #include <Geom_BezierCurve.hxx>
48 #include <ShapeExtend_ComplexCurve.hxx>
49 #include <Geom2d_Conic.hxx>
50 #include <Geom2d_TrimmedCurve.hxx>
51 #include <Geom2d_BSplineCurve.hxx>
52 #include <Geom2d_BezierCurve.hxx>
54 #include <Geom2d_OffsetCurve.hxx>
55 #include <Geom2dInt_Geom2dCurveTool.hxx>
56 #include <Geom2dAdaptor_Curve.hxx>
57 #include <Geom_Circle.hxx>
60 //=======================================================================
61 //function : ProjectOnSegments
63 //=======================================================================
65 static void ProjectOnSegments (const Adaptor3d_Curve& AC, const gp_Pnt& P3D,
66 const Standard_Integer nbseg,
67 Standard_Real& uMin, Standard_Real& uMax,
68 Standard_Real& distmin, gp_Pnt& proj, Standard_Real& param)
70 // On considere <nbseg> points sur [uMin,uMax]
71 // Quel est le plus proche. Et quel est le nouvel intervalle
72 // (il ne peut pas deborder l ancien)
73 Standard_Real u, dist, delta = (nbseg == 0)? 0 : (uMax-uMin)/nbseg; //szv#4:S4163:12Mar99 anti-exception
74 for (Standard_Integer i = 0; i <= nbseg; i ++) {
75 u = uMin + (delta * i);
76 gp_Pnt PU = AC.Value (u);
77 dist = PU.Distance (P3D);
78 if (dist < distmin) { distmin = dist; proj = PU; param = u; }
82 cout<<"ShapeAnalysis_Geom:Project, param="<<param<<" -> distmin="<<distmin<<endl;
85 uMax = Min (uMax, param+delta);
86 uMin = Max (uMin, param-delta);
90 //=======================================================================
91 //function : CurveNewton
92 //purpose : Newton algo on curve (study S4030)
93 //=======================================================================
95 static Standard_Boolean CurveNewton(const Standard_Real paramPrev,
96 const Adaptor3d_Curve& Ad,
98 const Standard_Real /*preci*/,
100 const Standard_Real cf,
101 const Standard_Real cl)
103 //szv#4:S4163:12Mar99 optimized
104 Standard_Real uMin = cf, uMax = cl;
106 Standard_Real Tol = Precision::Confusion();
107 Standard_Real Tol2 = Tol*Tol, rs2p = 1.e+10;
108 Standard_Real X = paramPrev;
110 for ( Standard_Integer i=0; i<20; i++) {
113 Ad.D2 (X, pnt, v1, v2);
114 Standard_Real nv1 = v1.SquareMagnitude();
115 if (nv1 < 1e-10) break;
117 gp_Vec rs (P3D, pnt);
118 Standard_Real D = nv1 + (rs * v2);
119 if ( fabs( D ) <1e-10) break;
121 Standard_Real dX = -(rs *v1)/D;
124 Standard_Real rs2 = rs.SquareMagnitude();
125 if (rs2 > 4.*rs2p) break;
128 if ( fabs(dX) > 1e-12) continue;
130 if (X < uMin || X > uMax) break;
132 Standard_Real rsn = rs * v1;
133 if (rsn*rsn/nv1 > Tol2) break;
135 return Standard_True;
138 // cout <<"NewtonC: failed after " << i+1 << " iterations (fail " << fail << " )" << endl;
139 return Standard_False;
142 //=======================================================================
145 //=======================================================================
147 Standard_Real ShapeAnalysis_Curve::Project(const Handle(Geom_Curve)& C3D,
149 const Standard_Real preci,
151 Standard_Real& param,
152 const Standard_Boolean AdjustToEnds) const
154 Standard_Real uMin = C3D->FirstParameter();
155 Standard_Real uMax = C3D->LastParameter();
156 if (uMin < uMax) return Project (C3D,P3D,preci,proj,param,uMin,uMax,AdjustToEnds);
157 else return Project (C3D,P3D,preci,proj,param,uMax,uMin,AdjustToEnds);
160 //=======================================================================
163 //=======================================================================
165 Standard_Real ShapeAnalysis_Curve::Project(const Handle(Geom_Curve)& C3D,
167 const Standard_Real preci,
169 Standard_Real& param,
170 const Standard_Real cf,
171 const Standard_Real cl,
172 const Standard_Boolean AdjustToEnds) const
174 Standard_Real distmin;
175 Standard_Real uMin = (cf < cl ? cf : cl);
176 Standard_Real uMax = (cf < cl ? cl : cf);
178 if (C3D->IsKind(STANDARD_TYPE(Geom_BoundedCurve))) {
179 Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
180 gp_Pnt LowBound = C3D->Value(uMin);
181 gp_Pnt HigBound = C3D->Value(uMax);
182 distmin = LowBound.Distance(P3D);
183 if (distmin <= prec) {
188 distmin = HigBound.Distance(P3D);
189 if (distmin <= prec) {
196 GeomAdaptor_Curve GAC(C3D, uMin, uMax);
197 if (!C3D->IsClosed()) {
198 //modified by rln on 16/12/97 after CSR# PRO11641 entity 20767
199 //the VIso was not closed (according to C3D->IsClosed()) while it "almost"
200 //was (the distance between ends of the curve was a little bit more than
201 //Precision::Confusion())
202 //in that case value 0.1 was too much and this method returned not correct parameter
205 // modified by pdn on 01.07.98 after BUC60195 entity 1952 (Min() added)
206 Standard_Real delta = Min (GAC.Resolution (preci), (uMax - uMin) * 0.1);
209 GAC.Load(C3D, uMin, uMax);
212 return ProjectAct(GAC, P3D, preci, proj, param);
215 //=======================================================================
218 //=======================================================================
220 Standard_Real ShapeAnalysis_Curve::Project(const Adaptor3d_Curve& C3D,
222 const Standard_Real preci,
224 Standard_Real& param,
225 const Standard_Boolean AdjustToEnds) const
228 Standard_Real uMin = C3D.FirstParameter();
229 Standard_Real uMax = C3D.LastParameter();
230 Standard_Real distmin;
231 if (!Precision::IsInfinite(uMin)||!Precision::IsInfinite(uMax)) {
232 Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
233 gp_Pnt LowBound = C3D.Value(uMin);
234 gp_Pnt HigBound = C3D.Value(uMax);
235 distmin = LowBound.Distance(P3D);
236 if (distmin <= prec) {
241 distmin = HigBound.Distance(P3D);
242 if (distmin <= prec) {
248 return ProjectAct(C3D, P3D, preci, proj, param);
251 //=======================================================================
252 //function : ProjectAct
254 //=======================================================================
256 Standard_Real ShapeAnalysis_Curve::ProjectAct(const Adaptor3d_Curve& C3D,
258 const Standard_Real preci,
260 Standard_Real& param) const
263 Standard_Boolean useExtrema = Standard_True; //:i5
266 //:i5 abv 11 Sep 98: UKI60195.igs on NT: hanging on null-length bsplines
267 if (C3D.IsClosed() && C3D.GetType()==GeomAbs_BSplineCurve) {
268 Handle(Geom_BSplineCurve) bspl = C3D.BSpline();
269 Standard_Real prec2 = gp::Resolution();
271 gp_Pnt p0 = bspl->Pole(1), p1;
272 for ( Standard_Integer i=2; i <= bspl->NbPoles(); i++, p0 = p1 ) {
274 if ( p0.SquareDistance ( p1 ) <= prec2 ) {
275 useExtrema = Standard_False;
282 Standard_Boolean OK = Standard_False;
283 if ( useExtrema ) { //:i5 //szv#4:S4163:12Mar99 try moved into if
286 Extrema_ExtPC myExtPC(P3D,C3D);
287 //Standard_Boolean done = myExtPC.IsDone() && ( myExtPC.NbExt() > 0); //szv#4:S4163:12Mar99 not needed
288 if ( myExtPC.IsDone() && ( myExtPC.NbExt() > 0) ) {
289 Standard_Real dist2, dist2Min = myExtPC.SquareDistance(1);
290 Standard_Integer index = 1;
291 for ( Standard_Integer i = 2; i <= myExtPC.NbExt(); i++) {
292 dist2 = myExtPC.SquareDistance(i);
293 if ( dist2 < dist2Min) { dist2Min = dist2; index = i; }
295 param = (myExtPC.Point(index)).Parameter();
296 proj = (myExtPC.Point(index)).Value();
300 catch(Standard_Failure) {
303 cout << "\nWarning: ShapeAnalysis_Curve::ProjectAct(): Exception in Extrema_ExtPC: ";
304 Standard_Failure::Caught()->Print(cout); cout << endl;
309 //szv#4:S4163:12Mar99 moved
310 Standard_Real uMin = C3D.FirstParameter(), uMax = C3D.LastParameter();
311 Standard_Boolean closed = Standard_False; // si on franchit les bornes ...
312 Standard_Real distmin = RealLast(), valclosed = 0.;
313 Standard_Real aModParam = param;
314 Standard_Real aModMin = distmin;
316 // PTV 29.05.2002 remember the old solution, cause it could be better
317 Standard_Real anOldParam =0.;
318 Standard_Boolean IsHaveOldSol = Standard_False;
321 IsHaveOldSol = Standard_True;
324 distmin = proj.Distance (P3D);
326 if (distmin > preci) OK = Standard_False;
327 // Cas TrimmedCurve a cheval. Voir la courbe de base.
328 // Si fermee, passer une periode
329 if (C3D.IsClosed()) {
330 closed = Standard_True;
331 valclosed = uMax - uMin; //szv#4:S4163:12Mar99 optimized
336 // BUG NICOLAS - Si le point est sur la courbe 0 Solutions
337 // Cela fonctionne avec ElCLib
339 // D une facon generale, on essaie de TOUJOURS retourner un resultat
340 // MEME PAS BIEN BON. L appelant pourra decider alors quoi faire
343 switch(C3D.GetType()) {
346 const gp_Circ& aCirc = C3D.Circle();
347 proj = aCirc.Position().Location();
348 if(aCirc.Radius() <= gp::Resolution() ||
349 P3D.SquareDistance(proj) <= gp::Resolution() ) {
350 param = C3D.FirstParameter();
351 proj = proj.XYZ() + aCirc.XAxis().Direction().XYZ() * aCirc.Radius();
354 param = ElCLib::Parameter(aCirc, P3D);
355 proj = ElCLib::Value(param, aCirc);
357 closed = Standard_True;
361 case GeomAbs_Hyperbola:
363 param = ElCLib::Parameter(C3D.Hyperbola(), P3D);
364 proj = ElCLib::Value(param, C3D.Hyperbola());
367 case GeomAbs_Parabola:
369 param = ElCLib::Parameter(C3D.Parabola(), P3D);
370 proj = ElCLib::Value(param, C3D.Parabola());
375 param = ElCLib::Parameter(C3D.Line(), P3D);
376 proj = ElCLib::Value(param, C3D.Line());
379 case GeomAbs_Ellipse:
381 param = ElCLib::Parameter(C3D.Ellipse(), P3D);
382 proj = ElCLib::Value(param, C3D.Ellipse());
383 closed = Standard_True;
390 // on ne va quand meme pas se laisser abattre ... ???
391 // on tente ceci : 21 points sur la courbe, quel est le plus proche
392 distmin = Precision::Infinite();
393 ProjectOnSegments (C3D,P3D,25,uMin,uMax,distmin,proj,param);
394 if (distmin <= preci) return distmin;
395 if (CurveNewton(param, C3D, P3D, preci, param, uMin, uMax)) { //:S4030
397 Standard_Real aDistNewton = P3D.Distance(proj);
398 if(aDistNewton < aModMin)
401 // cout <<"newton failed"<<endl;
402 ProjectOnSegments (C3D,P3D,40,uMin,uMax,distmin,proj,param);
403 if (distmin <= preci) return distmin;
404 ProjectOnSegments (C3D,P3D,20,uMin,uMax,distmin,proj,param);
405 if (distmin <= preci) return distmin;
406 ProjectOnSegments (C3D,P3D,25,uMin,uMax,distmin,proj,param);
407 if (distmin <= preci) return distmin;
408 ProjectOnSegments (C3D,P3D,40,uMin,uMax,distmin,proj,param);
409 if (distmin <= preci) return distmin;
410 // soyons raisonnable ...
411 if(distmin > aModMin) {
421 // p = PPOC.LowerDistanceParameter(); cf try
422 if ( closed && ( param < uMin || param > uMax ) )
423 param += ShapeAnalysis::AdjustByPeriod ( param, 0.5 * ( uMin + uMax ), valclosed );
426 // PTV 29.05.2002 Compare old solution and new;
427 Standard_Real adist1, adist2;
428 adist1 = anOldProj.SquareDistance(P3D);
429 adist2 = proj.SquareDistance (P3D);
430 if (adist1 < adist2) {
435 return proj.Distance (P3D);
439 //=======================================================================
440 //function : NextProject
441 //purpose : Newton algo for projecting point on curve (S4030)
442 //=======================================================================
444 Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
445 const Handle(Geom_Curve)& C3D,
447 const Standard_Real preci,
449 Standard_Real& param,
450 const Standard_Real cf,
451 const Standard_Real cl,
452 const Standard_Boolean AdjustToEnds) const
454 Standard_Real uMin = (cf < cl ? cf : cl);
455 Standard_Real uMax = (cf < cl ? cl : cf);
456 Standard_Real distmin;
457 if (C3D->IsKind(STANDARD_TYPE(Geom_BoundedCurve))) {
458 Standard_Real prec = ( AdjustToEnds ? preci : Precision::Confusion() ); //:j8 abv 10 Dec 98: tr10_r0501_db.stp #9423: protection against densing of points near one end
459 gp_Pnt LowBound = C3D->Value(uMin);
460 gp_Pnt HigBound = C3D->Value(uMax);
461 distmin = LowBound.Distance(P3D);
462 if (distmin <= prec) {
467 distmin = HigBound.Distance(P3D);
468 if (distmin <= prec) {
475 GeomAdaptor_Curve GAC(C3D, uMin, uMax);
476 if (!C3D->IsClosed()) {
477 //modified by rln on 16/12/97 after CSR# PRO11641 entity 20767
478 //the VIso was not closed (according to C3D->IsClosed()) while it "almost"
479 //was (the distance between ends of the curve was a little bit more than
480 //Precision::Confusion())
481 //in that case value 0.1 was too much and this method returned not correct parameter
484 // modified by pdn on 01.07.98 after BUC60195 entity 1952 (Min() added)
485 Standard_Real delta = Min (GAC.Resolution (preci), (uMax - uMin) * 0.1);
488 GAC.Load(C3D, uMin, uMax);
490 return NextProject ( paramPrev, GAC, P3D, preci, proj, param );
493 //=======================================================================
494 //function : NextProject
496 //=======================================================================
498 Standard_Real ShapeAnalysis_Curve::NextProject(const Standard_Real paramPrev,
499 const Adaptor3d_Curve& C3D,
501 const Standard_Real preci,
503 Standard_Real& param) const
505 Standard_Real uMin = C3D.FirstParameter();
506 Standard_Real uMax = C3D.LastParameter();
508 if (CurveNewton(paramPrev, C3D, P3D, preci, param, uMin, uMax)) {
510 return P3D.Distance(proj);
513 return Project(C3D, P3D, preci, proj, param);
516 //=======================================================================
517 //function : AdjustParameters
518 //purpose : Copied from StepToTopoDS_GeometricTuul::UpdateParam3d (Aug 2001)
519 //=======================================================================
521 Standard_Boolean ShapeAnalysis_Curve::ValidateRange (const Handle(Geom_Curve)& theCurve,
522 Standard_Real& First, Standard_Real& Last,
523 const Standard_Real preci) const
525 // First et/ou Last peuvent etre en dehors des bornes naturelles de la courbe.
526 // On donnera alors la valeur en bout a First et/ou Last
528 Standard_Real cf = theCurve->FirstParameter();
529 Standard_Real cl = theCurve->LastParameter();
530 // Standard_Real preci = BRepAPI::Precision();
532 if (theCurve->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) && !theCurve->IsClosed()) {
535 cout << "Update Edge First Parameter to Curve First Parameter" << endl;
539 else if (First > cl) {
541 cout << "Update Edge First Parameter to Curve Last Parameter" << endl;
547 cout << "Update Edge Last Parameter to Curve First Parameter" << endl;
551 else if (Last > cl) {
553 cout << "Update Edge Last Parameter to Curve Last Parameter" << endl;
559 if (First < Last) return Standard_True;
561 // 15.11.2002 PTV OCC966
562 if (ShapeAnalysis_Curve::IsPeriodic(theCurve))
563 ElCLib::AdjustPeriodic(cf,cl,Precision::PConfusion(),First,Last); //:a7 abv 11 Feb 98: preci -> PConfusion()
564 else if (theCurve->IsClosed()) {
565 // l'un des points projecte se trouve sur l'origine du parametrage
566 // de la courbe 3D. L algo a donne cl +- preci au lieu de cf ou vice-versa
567 // DANGER precision 3d applique a une espace 1d
569 // Last = cf au lieu de Last = cl
570 if (Abs(Last - cf) < Precision::PConfusion() /*preci*/) Last = cl ;
571 // First = cl au lieu de First = cf
572 else if (Abs(First - cl) < Precision::PConfusion() /*preci*/) First = cf;
574 // on se trouve dans un cas ou l origine est traversee
575 // illegal sur une courbe fermee non periodique
576 // on inverse quand meme les parametres !!!!!!
578 //:S4136 abv 20 Apr 99: r0701_ug.stp #6230: add check in 3d
579 if ( theCurve->Value(First).Distance(theCurve->Value(cf)) < preci ) First = cf;
580 if ( theCurve->Value(Last).Distance(theCurve->Value(cl)) < preci ) Last = cl;
581 if ( First > Last ) {
583 cout << "Warning : parameter range of edge crossing non periodic curve origin" << endl;
585 Standard_Real tmp = First;
591 // The curve is closed within the 3D tolerance
592 else if (theCurve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
593 Handle(Geom_BSplineCurve) aBSpline =
594 Handle(Geom_BSplineCurve)::DownCast(theCurve);
595 if (aBSpline->StartPoint().Distance(aBSpline->EndPoint()) <= preci ) {
596 //:S4136 <= BRepAPI::Precision()) {
597 // l'un des points projecte se trouve sur l'origine du parametrage
598 // de la courbe 3D. L algo a donne cl +- preci au lieu de cf ou vice-versa
599 // DANGER precision 3d applique a une espace 1d
601 // Last = cf au lieu de Last = cl
602 if (Abs(Last - cf) < Precision::PConfusion() /*preci*/) Last = cl ;
603 // First = cl au lieu de First = cf
604 else if (Abs(First - cl) < Precision::PConfusion() /*preci*/) First = cf;
606 // on se trouve dans un cas ou l origine est traversee
607 // illegal sur une courbe fermee non periodique
608 // on inverse quand meme les parametres !!!!!!
611 cout << "Warning : parameter range of edge crossing non periodic curve origin" << endl;
613 Standard_Real tmp = First;
618 //abv 15.03.00 #72 bm1_pe_t4 protection of exceptions in draw
619 else if ( First > Last ) {
621 cout << "Warning: parameter range is bad; curve reversed" << endl;
623 First = theCurve->ReversedParameter ( First );
624 Last = theCurve->ReversedParameter ( Last );
627 //:j9 abv 11 Dec 98: PRO7747 #4875, after :j8: else
628 if (First == Last) { //gka 10.07.1998 file PRO7656 entity 33334
629 First = cf; Last = cl;
630 return Standard_False;
635 cout << "UpdateParam3d Failed" << endl;
636 cout << " - Curve Type : " << theCurve->DynamicType() << endl;
637 cout << " - Param 1 : " << First << endl;
638 cout << " - Param 2 : " << Last << endl;
640 //abv 15.03.00 #72 bm1_pe_t4 protection of exceptions in draw
641 if ( First > Last ) {
643 cout << "Warning: parameter range is bad; curve reversed" << endl;
645 First = theCurve->ReversedParameter ( First );
646 Last = theCurve->ReversedParameter ( Last );
649 //pdn 11.01.99 #144 bm1_pe_t4 protection of exceptions in draw
651 First -= Precision::PConfusion();
652 Last += Precision::PConfusion();
654 return Standard_False;
656 return Standard_True;
659 //=======================================================================
660 //function : FillBndBox
661 //purpose : WORK-AROUND for methods brom BndLib which give not exact bounds
662 //=======================================================================
664 // search for extremum using Newton
665 static Standard_Integer SearchForExtremum (const Handle(Geom2d_Curve)& C2d,
666 const Standard_Real First,
667 const Standard_Real Last,
672 Standard_Real prevpar;
673 for ( Standard_Integer i=0; i <10; i++ ) {
677 C2d->D2 ( par, res, D1, D2 );
678 Standard_Real Det = ( D2 * dir );
679 if ( Abs ( Det ) < 1e-10 ) return Standard_True;
681 par -= ( D1 * dir ) / Det;
682 if ( Abs ( par - prevpar ) < Precision::PConfusion() ) return Standard_True;
684 if ( First - par >= Precision::PConfusion() ||
685 par - Last >= Precision::PConfusion() ) return Standard_False;
687 return Standard_True;
690 void ShapeAnalysis_Curve::FillBndBox (const Handle(Geom2d_Curve)& C2d,
691 const Standard_Real First,
692 const Standard_Real Last,
693 const Standard_Integer NPoints,
694 const Standard_Boolean Exact,
695 Bnd_Box2d &Box) const
697 Standard_Integer nseg = ( NPoints <2 ? 1 : NPoints-1 );
698 Standard_Real step = ( Last - First ) / nseg;
699 for ( Standard_Integer i=0; i <= nseg; i++ ) {
700 Standard_Real par = First + i * step;
701 gp_Pnt2d pnt = C2d->Value ( par );
703 if ( ! Exact ) continue;
706 Standard_Real parextr = par;
707 if ( SearchForExtremum ( C2d, Max(First,par-2.*step), Min(Last,par+2.*step),
708 gp_Vec2d(1,0), parextr, pextr ) ) {
712 if ( SearchForExtremum ( C2d, Max(First,par-2.*step), Min(Last,par+2.*step),
713 gp_Vec2d(0,1), parextr, pextr ) ) {
719 //=======================================================================
720 //function : SelectForwardSeam
722 //=======================================================================
724 Standard_Integer ShapeAnalysis_Curve::SelectForwardSeam(const Handle(Geom2d_Curve)& C1,
725 const Handle(Geom2d_Curve)& C2) const
727 // SelectForward est destine a devenir un outil distinct
728 // Il est sans doute optimisable !
730 Standard_Integer theCurveIndice = 0;
732 Handle(Geom2d_Line) L1 = Handle(Geom2d_Line)::DownCast(C1);
734 // if we have BoundedCurve, create a line from C1
735 Handle(Geom2d_BoundedCurve) BC1 = Handle(Geom2d_BoundedCurve)::DownCast(C1);
736 if (BC1.IsNull()) return theCurveIndice;
737 gp_Pnt2d StartBC1 = BC1->StartPoint();
738 gp_Pnt2d EndBC1 = BC1->EndPoint();
739 gp_Vec2d VecBC1(StartBC1, EndBC1);
740 L1 = new Geom2d_Line(StartBC1, VecBC1);
743 Handle(Geom2d_Line) L2 = Handle(Geom2d_Line)::DownCast(C2);
745 // if we have BoundedCurve, creates a line from C2
746 Handle(Geom2d_BoundedCurve) BC2 = Handle(Geom2d_BoundedCurve)::DownCast(C2);
747 if (BC2.IsNull()) return theCurveIndice;
748 gp_Pnt2d StartBC2 = BC2->StartPoint();
749 gp_Pnt2d EndBC2 = BC2->EndPoint();
750 gp_Vec2d VecBC2(StartBC2, EndBC2);
751 L2 = new Geom2d_Line(StartBC2, VecBC2);
754 Standard_Boolean UdirPos, UdirNeg, VdirPos, VdirNeg;
755 UdirPos = UdirNeg = VdirPos = VdirNeg = Standard_False;
757 gp_Dir2d theDir = L1->Direction();
758 gp_Pnt2d theLoc1 = L1->Location();
759 gp_Pnt2d theLoc2 = L2->Location();
761 if (theDir.X() > 0.) {
762 UdirPos = Standard_True; //szv#4:S4163:12Mar99 Udir unused
763 } else if (theDir.X() < 0.) {
764 UdirNeg = Standard_True; //szv#4:S4163:12Mar99 Udir unused
765 } else if (theDir.Y() > 0.) {
766 VdirPos = Standard_True; //szv#4:S4163:12Mar99 Vdir unused
767 } else if (theDir.Y() < 0.) {
768 VdirNeg = Standard_True; //szv#4:S4163:12Mar99 Vdir unused
772 // max of Loc1.X() Loc2.X()
773 if (theLoc1.X() > theLoc2.X()) theCurveIndice = 1;
774 else theCurveIndice = 2;
775 } else if (VdirNeg) {
776 if (theLoc1.X() > theLoc2.X()) theCurveIndice = 2;
777 else theCurveIndice = 1;
778 } else if (UdirPos) {
779 // min of Loc1.X() Loc2.X()
780 if (theLoc1.Y() < theLoc2.Y()) theCurveIndice = 1;
781 else theCurveIndice = 2;
782 } else if (UdirNeg) {
783 if (theLoc1.Y() < theLoc2.Y()) theCurveIndice = 2;
784 else theCurveIndice = 1;
787 return theCurveIndice;
791 //=============================================================================
792 // Static methods for IsPlanar
794 //=============================================================================
796 static gp_XYZ GetAnyNormal ( gp_XYZ orig )
799 if ( Abs ( orig.Z() ) < Precision::Confusion() )
800 Norm.SetCoord ( 0, 0, 1 );
802 Norm.SetCoord ( orig.Z(), 0, -orig.X() );
803 Standard_Real nrm = Norm.Modulus();
804 if ( nrm < Precision::Confusion() ) Norm.SetCoord ( 0, 0, 1 );
805 else Norm = Norm / nrm;
810 //=======================================================================
811 //function : GetSamplePoints
813 //=======================================================================
814 static void AppendControlPoles (TColgp_SequenceOfPnt& seq,
815 const Handle(Geom_Curve) curve)
817 if ( curve->IsKind(STANDARD_TYPE(Geom_Line))) {
818 seq.Append(curve->Value(0));
819 seq.Append(curve->Value(1));
820 } else if ( curve->IsKind(STANDARD_TYPE(Geom_Conic))) {
821 seq.Append(curve->Value(0));
822 seq.Append(curve->Value(M_PI/2));
823 seq.Append(curve->Value(M_PI));
824 } else if ( curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
825 //DeclareAndCast(Geom_TrimmedCurve, Trimmed, curve);
826 Handle(Geom_TrimmedCurve) Trimmed = *((Handle(Geom_TrimmedCurve) *) &curve);
827 // AppendControlPoles(seq,Trimmed->BasisCurve());
828 Handle(Geom_Curve) aBaseCrv = Trimmed->BasisCurve();
829 Standard_Boolean done = Standard_False;
830 if ( aBaseCrv->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
833 Handle(Geom_Geometry) Ctmp = aBaseCrv->Copy();
834 Handle(Geom_BSplineCurve) bslp = Handle(Geom_BSplineCurve)::DownCast(Ctmp);
835 bslp->Segment(curve->FirstParameter(), curve->LastParameter());
836 AppendControlPoles(seq,bslp);
837 done = Standard_True;
839 catch (Standard_Failure) {
842 else if ( aBaseCrv->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
845 Handle(Geom_Geometry) Ctmp = aBaseCrv->Copy();
846 Handle(Geom_BezierCurve) bz = Handle(Geom_BezierCurve)::DownCast(Ctmp);
847 bz->Segment(curve->FirstParameter(), curve->LastParameter());
848 AppendControlPoles(seq,bz);
849 done = Standard_True;
851 catch (Standard_Failure) {
855 seq.Append(curve->Value(curve->FirstParameter()));
856 seq.Append(curve->Value((curve->FirstParameter() + curve->LastParameter())/2.));
857 seq.Append(curve->Value(curve->LastParameter()));
859 } else if ( curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
860 //DeclareAndCast(Geom_OffsetCurve, OffsetC, curve);
861 Handle(Geom_OffsetCurve) OffsetC = *((Handle(Geom_OffsetCurve) *) &curve);
862 // AppendControlPoles(seq,OffsetC->BasisCurve());
863 seq.Append(curve->Value(curve->FirstParameter()));
864 seq.Append(curve->Value((curve->FirstParameter() + curve->LastParameter())/2.));
865 seq.Append(curve->Value(curve->LastParameter()));
866 } else if ( curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
867 //DeclareAndCast(Geom_BSplineCurve, BSpline, curve);
868 Handle(Geom_BSplineCurve) BSpline = *((Handle(Geom_BSplineCurve) *) &curve);
869 TColgp_Array1OfPnt Poles(1,BSpline->NbPoles());
870 BSpline->Poles(Poles);
871 for(Standard_Integer i = 1; i <= BSpline->NbPoles(); i++)
872 seq.Append(Poles(i));
873 } else if ( curve->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
874 //DeclareAndCast(Geom_BezierCurve, Bezier, curve);
875 //Handle(Geom_BezierCurve) Bezier = Handle(Geom_BezierCurve)::DownCast(curve);
876 Handle(Geom_BezierCurve) Bezier = *((Handle(Geom_BezierCurve) *) &curve);
877 TColgp_Array1OfPnt Poles(1,Bezier->NbPoles());
878 Bezier->Poles(Poles);
879 for(Standard_Integer i = 1; i <= Bezier->NbPoles(); i++)
880 seq.Append(Poles(i));
886 //=======================================================================
887 //function : IsPlanar
888 //purpose : Detects if points lie in some plane and returns normal
889 //=======================================================================
891 Standard_Boolean ShapeAnalysis_Curve::IsPlanar (const TColgp_Array1OfPnt& pnts,
893 const Standard_Real preci)
895 Standard_Real precision = (preci > 0.0)? preci : Precision::Confusion();
896 Standard_Boolean noNorm = (Normal.SquareModulus() == 0);
898 if (pnts.Length() < 3) {
899 gp_XYZ N1 = pnts(1).XYZ()-pnts(2).XYZ();
901 Normal = GetAnyNormal(N1);
902 return Standard_True;
904 return Abs(N1*Normal) < Precision::Confusion();
909 //define a center point
910 gp_XYZ aCenter(0,0,0);
911 Standard_Integer i = 1;
912 for (; i <= pnts.Length(); i++)
913 aCenter += pnts(i).XYZ();
914 aCenter/=pnts.Length();
917 aMaxDir = pnts(1).XYZ() - aCenter;
918 Normal = (pnts(pnts.Length()).XYZ() - aCenter) ^ aMaxDir;
920 for ( i = 1; i < pnts.Length(); i++) {
921 gp_XYZ aTmpDir = pnts(i+1).XYZ() - aCenter;
922 if(aTmpDir.SquareModulus() > aMaxDir.SquareModulus())
925 gp_XYZ aDelta = (pnts(i).XYZ() - aCenter) ^ (pnts(i+1).XYZ() - aCenter);
926 if(Normal*aDelta < 0)
932 // check if points are linear
933 Standard_Real nrm = Normal.Modulus();
934 if ( nrm < Precision::Confusion() ) {
935 Normal = GetAnyNormal(aMaxDir);
936 return Standard_True;
938 Normal = Normal / nrm;
940 Standard_Real mind = RealLast(), maxd = -RealLast(), dev;
941 for (Standard_Integer i = 1; i <= pnts.Length(); i++) {
942 dev = pnts(i).XYZ() * Normal;
943 if (dev < mind) mind = dev;
944 if (dev > maxd) maxd = dev;
947 return ((maxd - mind) <= precision);
951 //=======================================================================
952 //function : IsPlanar
954 //=======================================================================
956 Standard_Boolean ShapeAnalysis_Curve::IsPlanar (const Handle(Geom_Curve)& curve,
958 const Standard_Real preci)
960 Standard_Real precision = (preci > 0.0)? preci : Precision::Confusion();
961 Standard_Boolean noNorm = (Normal.SquareModulus() == 0);
963 if (curve->IsKind(STANDARD_TYPE(Geom_Line))) {
964 //DeclareAndCast(Geom_Line, Line, curve);
965 Handle(Geom_Line) Line = *((Handle(Geom_Line) *) &curve);
966 gp_XYZ N1 = Line->Position().Direction().XYZ();
968 Normal = GetAnyNormal(N1);
969 return Standard_True;
971 return Abs(N1*Normal) < Precision::Confusion();
974 if (curve->IsKind(STANDARD_TYPE(Geom_Conic))) {
975 //DeclareAndCast(Geom_Conic, Conic, curve);
976 Handle(Geom_Conic) Conic = *((Handle(Geom_Conic) *) &curve);
977 gp_XYZ N1 = Conic->Axis().Direction().XYZ();
980 return Standard_True;
982 gp_XYZ aVecMul = N1^Normal;
983 return aVecMul.SquareModulus() < Precision::SquareConfusion();
986 if (curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
987 //DeclareAndCast(Geom_TrimmedCurve, Trimmed, curve);
988 Handle(Geom_TrimmedCurve) Trimmed = *((Handle(Geom_TrimmedCurve) *) &curve);
989 return IsPlanar(Trimmed->BasisCurve(),Normal,precision);
992 if (curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
993 //DeclareAndCast(Geom_OffsetCurve, OffsetC, curve);
994 Handle(Geom_OffsetCurve) OffsetC = *((Handle(Geom_OffsetCurve) *) &curve);
995 return IsPlanar(OffsetC->BasisCurve(),Normal,precision);
998 if (curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
999 //DeclareAndCast(Geom_BSplineCurve, BSpline, curve);
1000 Handle(Geom_BSplineCurve) BSpline = *((Handle(Geom_BSplineCurve) *) &curve);
1001 TColgp_Array1OfPnt Poles(1,BSpline->NbPoles());
1002 BSpline->Poles(Poles);
1003 return IsPlanar(Poles,Normal,precision);
1006 if (curve->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
1007 //DeclareAndCast(Geom_BezierCurve, Bezier, curve);
1008 Handle(Geom_BezierCurve) Bezier = *((Handle(Geom_BezierCurve) *) &curve);
1009 TColgp_Array1OfPnt Poles(1,Bezier->NbPoles());
1010 Bezier->Poles(Poles);
1011 return IsPlanar(Poles,Normal,precision);
1014 if (curve->IsKind(STANDARD_TYPE(ShapeExtend_ComplexCurve))) {
1015 //DeclareAndCast(ShapeExtend_ComplexCurve, Complex, curve);
1016 Handle(ShapeExtend_ComplexCurve) Complex = *((Handle(ShapeExtend_ComplexCurve) *) &curve);
1017 TColgp_SequenceOfPnt sequence;
1018 Standard_Integer i; // svv Jan11 2000 : porting on DEC
1019 for (i = 1; i <= Complex->NbCurves(); i++)
1020 AppendControlPoles(sequence,Complex->Curve(i));
1021 TColgp_Array1OfPnt Poles(1,sequence.Length());
1022 for (i=1; i <= sequence.Length(); i++) Poles(i) = sequence(i);
1023 return IsPlanar(Poles,Normal,precision);
1026 return Standard_False;
1029 //=======================================================================
1030 //function : GetSamplePoints
1032 //=======================================================================
1034 Standard_Boolean ShapeAnalysis_Curve::GetSamplePoints (const Handle(Geom_Curve)& curve,
1035 const Standard_Real first,
1036 const Standard_Real last,
1037 TColgp_SequenceOfPnt& seq)
1039 Standard_Real adelta = curve->LastParameter() - curve->FirstParameter();
1041 return Standard_False;
1043 Standard_Integer aK = (Standard_Integer)ceil ((last - first) / adelta);
1044 Standard_Integer nbp =100*aK;
1045 if(curve->IsKind(STANDARD_TYPE(Geom_Line)))
1047 else if(curve->IsKind(STANDARD_TYPE(Geom_Circle)))
1050 else if (curve->IsKind(STANDARD_TYPE(Geom_BSplineCurve))) {
1051 Handle(Geom_BSplineCurve) aBspl = Handle(Geom_BSplineCurve)::DownCast(curve);
1053 nbp = aBspl->NbKnots() * aBspl->Degree()*aK;
1054 if(nbp < 2.0) nbp=2;
1056 else if (curve->IsKind(STANDARD_TYPE(Geom_BezierCurve))) {
1057 Handle(Geom_BezierCurve) aB = Handle(Geom_BezierCurve)::DownCast(curve);
1058 nbp = 3 + aB->NbPoles();
1060 else if(curve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) {
1061 Handle(Geom_OffsetCurve) aC = Handle(Geom_OffsetCurve)::DownCast(curve);
1062 return GetSamplePoints(aC->BasisCurve(),first,last,seq);
1064 else if(curve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
1065 Handle(Geom_TrimmedCurve) aC = Handle(Geom_TrimmedCurve)::DownCast(curve);
1066 return GetSamplePoints(aC->BasisCurve(),first,last,seq);
1068 Standard_Real step = ( last - first ) / (Standard_Real)( nbp - 1 );
1069 Standard_Real par = first, stop = last - 0.5 * step;
1070 for ( ; par < stop; par += step )
1071 seq.Append(curve->Value(par));
1072 seq.Append(curve->Value(last));
1073 return Standard_True;
1075 //=======================================================================
1076 //function : GetSamplePoints
1078 //=======================================================================
1080 Standard_Boolean ShapeAnalysis_Curve::GetSamplePoints (const Handle(Geom2d_Curve)& curve,
1081 const Standard_Real first,
1082 const Standard_Real last,
1083 TColgp_SequenceOfPnt2d& seq)
1085 //:abv 05.06.02: TUBE.stp
1086 // Use the same distribution of points as BRepTopAdaptor_FClass2d for consistency
1087 Geom2dAdaptor_Curve C ( curve, first, last );
1088 Standard_Integer nbs = Geom2dInt_Geom2dCurveTool::NbSamples(C);
1089 //-- Attention aux bsplines rationnelles de degree 3. (bouts de cercles entre autres)
1090 if (nbs > 2) nbs*=4;
1091 Standard_Real step = ( last - first ) / (Standard_Real)( nbs - 1 );
1092 Standard_Real par = first, stop = last - 0.5 * step;
1093 for ( ; par < stop; par += step )
1094 seq.Append(curve->Value(par));
1095 seq.Append(curve->Value(last));
1096 return Standard_True;
1101 if ( curve->IsKind(STANDARD_TYPE(Geom2d_Line))) {
1102 seq.Append(curve->Value(first));
1103 seq.Append(curve->Value(last));
1104 return Standard_True;
1106 else if(curve->IsKind(STANDARD_TYPE(Geom2d_Conic))) {
1107 step = Min ( M_PI, last-first ) / 19; //:abv 05.06.02 TUBE.stp #19209...: M_PI/16
1108 // if( step>(last-first) ) {
1109 // seq.Append(curve->Value(first));
1110 // seq.Append(curve->Value((last+first)/2));
1111 // seq.Append(curve->Value(last));
1112 // return Standard_True;
1115 Standard_Real par=first;
1116 for(i=0; par<last; i++) {
1117 seq.Append(curve->Value(par));
1120 seq.Append(curve->Value(last));
1121 return Standard_True;
1124 else if ( curve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) {
1125 DeclareAndCast(Geom2d_TrimmedCurve, Trimmed, curve);
1126 return GetControlPoints(Trimmed->BasisCurve(), seq, first, last);
1128 else if ( curve->IsKind(STANDARD_TYPE(Geom2d_BSplineCurve))) {
1129 DeclareAndCast(Geom2d_BSplineCurve, aBs, curve);
1130 TColStd_SequenceOfReal aSeqParam;
1132 aSeqParam.Append(first);
1133 for(i=1; i<=aBs->NbKnots(); i++) {
1134 if( aBs->Knot(i)>first && aBs->Knot(i)<last )
1135 aSeqParam.Append(aBs->Knot(i));
1137 aSeqParam.Append(last);
1138 Standard_Integer NbPoints=aBs->Degree();
1139 if( (aSeqParam.Length()-1)*NbPoints>10 ) {
1140 for(i=1; i<aSeqParam.Length(); i++) {
1141 Standard_Real FirstPar = aSeqParam.Value(i);
1142 Standard_Real LastPar = aSeqParam.Value(i+1);
1143 step = (LastPar-FirstPar)/NbPoints;
1144 for(Standard_Integer k=0; k<NbPoints; k++ ) {
1145 aBs->D0(FirstPar+step*k, Ptmp);
1149 aBs->D0(last, Ptmp);
1151 return Standard_True;
1154 step = (last-first)/10;
1155 for(Standard_Integer k=0; k<=10; k++ ) {
1156 aBs->D0(first+step*k, Ptmp);
1159 return Standard_True;
1163 else if ( curve->IsKind(STANDARD_TYPE(Geom2d_BezierCurve))) {
1164 DeclareAndCast(Geom2d_BezierCurve, aBz, curve);
1166 Standard_Integer NbPoints=aBz->Degree();
1167 step = (last-first)/NbPoints;
1168 for(Standard_Integer k=0; k<NbPoints; k++ ) {
1169 aBz->D0(first+step*k, Ptmp);
1172 aBz->D0(last, Ptmp);
1174 return Standard_True;
1177 return Standard_False;
1181 //=======================================================================
1182 //function : IsClosed
1184 //=======================================================================
1186 Standard_Boolean ShapeAnalysis_Curve::IsClosed(const Handle(Geom_Curve)& theCurve,
1187 const Standard_Real preci)
1189 if (theCurve->IsClosed())
1190 return Standard_True;
1192 Standard_Real prec = Max (preci, Precision::Confusion());
1195 f = theCurve->FirstParameter();
1196 l = theCurve->LastParameter();
1198 if (Precision::IsInfinite (f) || Precision::IsInfinite (l))
1199 return Standard_False;
1201 Standard_Real aClosedVal = theCurve->Value(f).SquareDistance(theCurve->Value(l));
1202 Standard_Real preci2 = prec*prec;
1204 return (aClosedVal <= preci2);
1206 //=======================================================================
1207 //function : IsPeriodic
1209 //=======================================================================
1211 Standard_Boolean ShapeAnalysis_Curve::IsPeriodic(const Handle(Geom_Curve)& theCurve)
1213 // 15.11.2002 PTV OCC966
1214 // remove regressions in DE tests (diva, divb, divc, toe3) in KAS:dev
1215 // ask IsPeriodic on BasisCurve
1216 Handle(Geom_Curve) aTmpCurve = theCurve;
1217 while ( (aTmpCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve))) ||
1218 (aTmpCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) ) {
1219 if (aTmpCurve->IsKind(STANDARD_TYPE(Geom_OffsetCurve)))
1220 aTmpCurve = Handle(Geom_OffsetCurve)::DownCast(aTmpCurve)->BasisCurve();
1221 if (aTmpCurve->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)))
1222 aTmpCurve = Handle(Geom_TrimmedCurve)::DownCast(aTmpCurve)->BasisCurve();
1224 return aTmpCurve->IsPeriodic();
1227 Standard_Boolean ShapeAnalysis_Curve::IsPeriodic(const Handle(Geom2d_Curve)& theCurve)
1229 // 15.11.2002 PTV OCC966
1230 // remove regressions in DE tests (diva, divb, divc, toe3) in KAS:dev
1231 // ask IsPeriodic on BasisCurve
1232 Handle(Geom2d_Curve) aTmpCurve = theCurve;
1233 while ( (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_OffsetCurve))) ||
1234 (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve))) ) {
1235 if (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_OffsetCurve)))
1236 aTmpCurve = Handle(Geom2d_OffsetCurve)::DownCast(aTmpCurve)->BasisCurve();
1237 if (aTmpCurve->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)))
1238 aTmpCurve = Handle(Geom2d_TrimmedCurve)::DownCast(aTmpCurve)->BasisCurve();
1240 return aTmpCurve->IsPeriodic();