1 // Created on: 1994-09-05
2 // Created by: Bruno DUMORTIER
3 // Copyright (c) 1994-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 // 09-Aug-95 : xab : changed the ProjLib_ProjectOnPlane in the case
18 // of the line and the parameteriation is kept
19 #include <ProjLib_ProjectOnPlane.hxx>
20 #include <AppCont_Function.hxx>
21 #include <Approx_FitAndDivide.hxx>
22 #include <AppParCurves_MultiCurve.hxx>
23 #include <Standard_ConstructionError.hxx>
24 #include <Standard_NoSuchObject.hxx>
25 #include <Standard_NotImplemented.hxx>
26 #include <Precision.hxx>
27 #include <BSplCLib.hxx>
29 #include <Geom_BSplineCurve.hxx>
30 #include <Geom_BezierCurve.hxx>
31 #include <TColgp_HArray1OfPnt.hxx>
32 #include <TColStd_HArray1OfReal.hxx>
33 #include <TColStd_HArray1OfInteger.hxx>
34 #include <Precision.hxx>
36 #include <Adaptor3d_Curve.hxx>
37 #include <Adaptor3d_HCurve.hxx>
38 #include <GeomAdaptor_HCurve.hxx>
39 #include <Geom_Line.hxx>
40 #include <GeomConvert.hxx>
41 #include <BSplCLib.hxx>
42 #include <Geom_TrimmedCurve.hxx>
43 #include <Geom_Circle.hxx>
44 #include <Geom_Parabola.hxx>
45 #include <Geom_Hyperbola.hxx>
46 #include <Geom_Ellipse.hxx>
47 #include <GeomLib_Tool.hxx>
48 #include <math_Jacobi.hxx>
49 #include <math_Matrix.hxx>
53 //=======================================================================
54 //function : OnPlane_Value
55 //purpose : Evaluate current point of the projected curve
56 //=======================================================================
58 static gp_Pnt OnPlane_Value(const Standard_Real U,
59 const Handle(Adaptor3d_HCurve)& aCurvePtr,
63 // PO . Z / Z = Pl.Direction()
64 // Proj(u) = P(u) + ------- * D avec \ O = Pl.Location()
67 gp_Pnt Point = aCurvePtr->Value(U);
69 gp_Vec PO(Point,Pl.Location());
71 Standard_Real Alpha = PO * gp_Vec(Pl.Direction());
72 Alpha /= D * Pl.Direction();
73 Point.SetXYZ(Point.XYZ() + Alpha * D.XYZ());
78 //=======================================================================
79 //function : OnPlane_DN
80 //purpose : Evaluate current point of the projected curve
81 //=======================================================================
83 static gp_Vec OnPlane_DN(const Standard_Real U,
84 const Standard_Integer DerivativeRequest,
85 const Handle(Adaptor3d_HCurve)& aCurvePtr,
89 // PO . Z / Z = Pl.Direction()
90 // Proj(u) = P(u) + ------- * D avec \ O = Pl.Location()
93 gp_Vec Vector = aCurvePtr->DN(U,DerivativeRequest);
95 gp_Dir Z = Pl.Direction();
98 Alpha = Vector * gp_Vec(Z);
101 Vector.SetXYZ( Vector.XYZ() - Alpha * D.XYZ());
105 //=======================================================================
106 //function : OnPlane_D1
108 //=======================================================================
110 static Standard_Boolean OnPlane_D1(const Standard_Real U,
113 const Handle(Adaptor3d_HCurve)& aCurvePtr,
121 gp_Dir Z = Pl.Direction();
123 aCurvePtr->D1(U, Point, Vector);
125 // evaluate the point as in `OnPlane_Value`
126 gp_Vec PO(Point,Pl.Location());
127 Alpha = PO * gp_Vec(Z);
129 P.SetXYZ(Point.XYZ() + Alpha * D.XYZ());
132 // evaluate the derivative.
134 // d(Proj) d(P) 1 d(P)
135 // ------ = --- - -------- * ( --- . Z ) * D
139 Alpha = Vector * gp_Vec(Z);
142 V.SetXYZ( Vector.XYZ() - Alpha * D.XYZ());
144 return Standard_True;
146 //=======================================================================
147 //function : OnPlane_D2
149 //=======================================================================
151 static Standard_Boolean OnPlane_D2(const Standard_Real U,
155 const Handle(Adaptor3d_HCurve) & aCurvePtr,
164 gp_Dir Z = Pl.Direction();
166 aCurvePtr->D2(U, Point, Vector1, Vector2);
168 // evaluate the point as in `OnPlane_Value`
169 gp_Vec PO(Point,Pl.Location());
170 Alpha = PO * gp_Vec(Z);
172 P.SetXYZ(Point.XYZ() + Alpha * D.XYZ());
174 // evaluate the derivative.
176 // d(Proj) d(P) 1 d(P)
177 // ------ = --- - -------- * ( --- . Z ) * D
181 Alpha = Vector1 * gp_Vec(Z);
184 V1.SetXYZ( Vector1.XYZ() - Alpha * D.XYZ());
186 Alpha = Vector2 * gp_Vec(Z);
189 V2.SetXYZ( Vector2.XYZ() - Alpha * D.XYZ());
190 return Standard_True;
193 //=======================================================================
194 //function : OnPlane_D3
196 //=======================================================================
198 static Standard_Boolean OnPlane_D3(const Standard_Real U,
203 const Handle(Adaptor3d_HCurve)& aCurvePtr,
213 gp_Dir Z = Pl.Direction();
215 aCurvePtr->D3(U, Point, Vector1, Vector2, Vector3);
217 // evaluate the point as in `OnPlane_Value`
218 gp_Vec PO(Point,Pl.Location());
219 Alpha = PO * gp_Vec(Z);
221 P.SetXYZ(Point.XYZ() + Alpha * D.XYZ());
223 // evaluate the derivative.
225 // d(Proj) d(P) 1 d(P)
226 // ------ = --- - -------- * ( --- . Z ) * D
230 Alpha = Vector1 * gp_Vec(Z);
233 V1.SetXYZ( Vector1.XYZ() - Alpha * D.XYZ());
235 Alpha = Vector2 * gp_Vec(Z);
238 V2.SetXYZ( Vector2.XYZ() - Alpha * D.XYZ());
239 Alpha = Vector3 * gp_Vec(Z);
242 V3.SetXYZ( Vector3.XYZ() - Alpha * D.XYZ());
243 return Standard_True;
246 //=======================================================================
247 // class : ProjLib_OnPlane
248 //purpose : Use to approximate the projection on a plane
249 //=======================================================================
251 class ProjLib_OnPlane : public AppCont_Function
254 Handle(Adaptor3d_HCurve) myCurve;
260 ProjLib_OnPlane(const Handle(Adaptor3d_HCurve)& C,
271 Standard_Real FirstParameter() const
272 {return myCurve->FirstParameter();}
274 Standard_Real LastParameter() const
275 {return myCurve->LastParameter();}
277 Standard_Boolean Value(const Standard_Real theT,
278 NCollection_Array1<gp_Pnt2d>& /*thePnt2d*/,
279 NCollection_Array1<gp_Pnt>& thePnt) const
281 thePnt(1) = OnPlane_Value(theT, myCurve, myPlane, myDirection);
282 return Standard_True;
285 Standard_Boolean D1(const Standard_Real theT,
286 NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
287 NCollection_Array1<gp_Vec>& theVec) const
290 return OnPlane_D1(theT, aDummyPnt, theVec(1),myCurve,myPlane,myDirection);
297 //=====================================================================//
299 // D E S C R I P T I O N O F T H E C L A S S : //
301 // P r o j L i b _ A p p r o x P r o j e c t O n P l a n e //
303 //=====================================================================//
307 //=======================================================================
308 //function : PerformApprox
310 //=======================================================================
312 static void PerformApprox (const Handle(Adaptor3d_HCurve)& C,
315 Handle(Geom_BSplineCurve) &BSplineCurvePtr)
318 ProjLib_OnPlane F(C,Pl,D);
320 Standard_Integer Deg1, Deg2;
323 Approx_FitAndDivide Fit(Deg1,Deg2,Precision::Approximation(),
324 Precision::PApproximation(),Standard_True);
325 Fit.SetMaxSegments(100);
327 if (!Fit.IsAllApproximated())
332 Standard_Integer NbCurves = Fit.NbMultiCurves();
333 Standard_Integer MaxDeg = 0;
335 // Pour transformer la MultiCurve en BSpline, il faut que toutes
336 // les Bezier la constituant aient le meme degre -> Calcul de MaxDeg
337 Standard_Integer NbPoles = 1;
338 for (i = 1; i <= NbCurves; i++) {
339 Standard_Integer Deg = Fit.Value(i).Degree();
340 MaxDeg = Max ( MaxDeg, Deg);
342 NbPoles = MaxDeg * NbCurves + 1; //Poles sur la BSpline
344 TColgp_Array1OfPnt Poles( 1, NbPoles);
346 TColgp_Array1OfPnt TempPoles( 1, MaxDeg + 1); //pour augmentation du degre
348 TColStd_Array1OfReal Knots( 1, NbCurves + 1); //Noeuds de la BSpline
350 Standard_Integer Compt = 1;
351 for (i = 1; i <= Fit.NbMultiCurves(); i++) {
352 Fit.Parameters(i, Knots(i), Knots(i+1));
354 AppParCurves_MultiCurve MC = Fit.Value( i); //Charge la Ieme Curve
355 TColgp_Array1OfPnt LocalPoles( 1, MC.Degree() + 1);//Recupere les poles
356 MC.Curve(1, LocalPoles);
358 //Augmentation eventuelle du degre
359 if (MaxDeg > MC.Degree() ) {
360 BSplCLib::IncreaseDegree(MaxDeg, LocalPoles, BSplCLib::NoWeights(),
361 TempPoles, BSplCLib::NoWeights());
362 //mise a jour des poles de la PCurve
363 for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
364 Poles.SetValue( Compt, TempPoles( j));
369 //mise a jour des poles de la PCurve
370 for (Standard_Integer j = 1 ; j <= MaxDeg + 1; j++) {
371 Poles.SetValue( Compt, LocalPoles( j));
379 //mise a jour des fields de ProjLib_Approx
382 NbKnots = NbCurves + 1;
384 TColStd_Array1OfInteger Mults(1, NbKnots);
385 Mults.SetValue( 1, MaxDeg + 1);
386 for ( i = 2; i <= NbCurves; i++) {
387 Mults.SetValue( i, MaxDeg);
389 Mults.SetValue( NbKnots, MaxDeg + 1);
391 new Geom_BSplineCurve(Poles,Knots,Mults,MaxDeg,Standard_False);
395 //=======================================================================
396 //function : ProjectOnPlane
398 //=======================================================================
400 ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane() :
401 myKeepParam(Standard_False),
405 myType (GeomAbs_OtherCurve),
406 myIsApprox (Standard_False)
410 //=======================================================================
411 //function : ProjLib_ProjectOnPlane
413 //=======================================================================
415 ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane(const gp_Ax3& Pl) :
417 myDirection (Pl.Direction()) ,
418 myKeepParam(Standard_False),
422 myType (GeomAbs_OtherCurve),
423 myIsApprox (Standard_False)
427 //=======================================================================
428 //function : ProjLib_ProjectOnPlane
430 //=======================================================================
432 ProjLib_ProjectOnPlane::ProjLib_ProjectOnPlane(const gp_Ax3& Pl,
436 myKeepParam(Standard_False),
440 myType (GeomAbs_OtherCurve),
441 myIsApprox (Standard_False)
443 // if ( Abs(D * Pl.Direction()) < Precision::Confusion()) {
444 // throw Standard_ConstructionError
445 // ("ProjLib_ProjectOnPlane: The Direction and the Plane are parallel");
449 //=======================================================================
451 //purpose : Returns the projection of a point <Point> on a plane
452 // <ThePlane> along a direction <TheDir>.
453 //=======================================================================
455 static gp_Pnt ProjectPnt(const gp_Ax3& ThePlane,
456 const gp_Dir& TheDir,
459 gp_Vec PO(Point,ThePlane.Location());
461 Standard_Real Alpha = PO * gp_Vec(ThePlane.Direction());
462 Alpha /= TheDir * ThePlane.Direction();
465 P.SetXYZ(Point.XYZ() + Alpha * TheDir.XYZ());
471 //=======================================================================
473 //purpose : Returns the projection of a Vector <Vec> on a plane
474 // <ThePlane> along a direction <TheDir>.
475 //=======================================================================
477 static gp_Vec ProjectVec(const gp_Ax3& ThePlane,
478 const gp_Dir& TheDir,
482 gp_Vec Z = ThePlane.Direction();
484 D -= ( (Vec * Z) / (TheDir * Z)) * TheDir;
489 //=======================================================================
492 //=======================================================================
494 void ProjLib_ProjectOnPlane::Load(const Handle(Adaptor3d_HCurve)& C,
495 const Standard_Real Tolerance,
496 const Standard_Boolean KeepParametrization)
500 myType = GeomAbs_OtherCurve;
501 myIsApprox = Standard_False;
502 myTolerance = Tolerance ;
504 Handle(Geom_BSplineCurve) ApproxCurve;
505 Handle(GeomAdaptor_HCurve) aGAHCurve;
507 Handle(Geom_Line) GeomLinePtr;
508 Handle(Geom_Circle) GeomCirclePtr ;
509 Handle(Geom_Ellipse) GeomEllipsePtr ;
510 Handle(Geom_Hyperbola) GeomHyperbolaPtr ;
516 Standard_Integer num_knots ;
517 GeomAbs_CurveType Type = C->GetType();
520 Standard_Real R1 =0., R2 =0.;
522 myKeepParam = KeepParametrization;
528 // ==> Q(u) = f(P(u))
529 // = f(O) + u * f(Xc)
531 gp_Lin L = myCurve->Line();
532 gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(L.Direction()));
534 if ( Xc.Magnitude() < Precision::Confusion()) { // line orthog au plan
535 myType = GeomAbs_BSplineCurve;
536 gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location());
537 TColStd_Array1OfInteger Mults(1,2); Mults.Init(2);
538 TColgp_Array1OfPnt Poles(1,2); Poles.Init(P);
539 TColStd_Array1OfReal Knots(1,2);
540 Knots(1) = myCurve->FirstParameter();
541 Knots(2) = myCurve->LastParameter();
542 Handle(Geom_BSplineCurve) BSP =
543 new Geom_BSplineCurve(Poles,Knots,Mults,1);
545 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
546 GeomAdaptor_Curve aGACurve(BSP);
547 myResult = new GeomAdaptor_HCurve(aGACurve);
548 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
550 else if ( Abs( Xc.Magnitude() - 1.) < Precision::Confusion()) {
551 myType = GeomAbs_Line;
552 gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location());
553 myFirstPar = myCurve->FirstParameter();
554 myLastPar = myCurve->LastParameter();
555 aLine = gp_Lin(P,gp_Dir(Xc));
556 GeomLinePtr = new Geom_Line(aLine);
558 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
559 GeomAdaptor_Curve aGACurve(GeomLinePtr,
560 myCurve->FirstParameter(),
561 myCurve->LastParameter() );
562 myResult = new GeomAdaptor_HCurve(aGACurve);
563 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
566 myType = GeomAbs_Line;
567 gp_Pnt P = ProjectPnt(myPlane,myDirection,L.Location());
568 aLine = gp_Lin(P,gp_Dir(Xc));
569 Standard_Real Udeb, Ufin;
571 // eval the first and last parameters of the projected curve
572 Udeb = myCurve->FirstParameter();
573 Ufin = myCurve->LastParameter();
574 gp_Pnt P1 = ProjectPnt(myPlane,myDirection,
575 myCurve->Value(Udeb));
576 gp_Pnt P2 = ProjectPnt(myPlane,myDirection,
577 myCurve->Value(Ufin));
578 myFirstPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P,P1));
579 myLastPar = gp_Vec(aLine.Direction()).Dot(gp_Vec(P,P2));
580 GeomLinePtr = new Geom_Line(aLine);
582 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
583 GeomAdaptor_Curve aGACurve(GeomLinePtr,
586 myResult = new GeomAdaptor_HCurve(aGACurve);
587 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
590 myType = GeomAbs_BSplineCurve;
592 // make a linear BSpline of degree 1 between the end points of
593 // the projected line
595 Handle(Geom_TrimmedCurve) NewTrimCurvePtr =
596 new Geom_TrimmedCurve(GeomLinePtr,
600 Handle(Geom_BSplineCurve) NewCurvePtr =
601 GeomConvert::CurveToBSplineCurve(NewTrimCurvePtr) ;
602 num_knots = NewCurvePtr->NbKnots() ;
603 TColStd_Array1OfReal BsplineKnots(1,num_knots) ;
604 NewCurvePtr->Knots(BsplineKnots) ;
606 BSplCLib::Reparametrize(myCurve->FirstParameter(),
607 myCurve->LastParameter(),
610 NewCurvePtr->SetKnots(BsplineKnots) ;
611 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
612 GeomAdaptor_Curve aGACurve(NewCurvePtr);
613 myResult = new GeomAdaptor_HCurve(aGACurve);
614 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
621 // Pour le cercle et l ellipse on a les relations suivantes:
622 // ( Rem : pour le cercle R1 = R2 = R)
623 // P(u) = O + R1 * Cos(u) * Xc + R2 * Sin(u) * Yc
624 // ==> Q(u) = f(P(u))
625 // = f(O) + R1 * Cos(u) * f(Xc) + R2 * Sin(u) * f(Yc)
627 gp_Circ Circ = myCurve->Circle();
628 Axis = Circ.Position();
629 R1 = R2 = Circ.Radius();
633 case GeomAbs_Ellipse:
635 if ( Type == GeomAbs_Ellipse) {
636 gp_Elips E = myCurve->Ellipse();
638 R1 = E.MajorRadius();
639 R2 = E.MinorRadius();
642 // Common Code for CIRCLE & ELLIPSE begin here
643 gp_Dir X = Axis.XDirection();
644 gp_Dir Y = Axis.YDirection();
645 gp_Vec VDx = ProjectVec(myPlane,myDirection,X);
646 gp_Vec VDy = ProjectVec(myPlane,myDirection,Y);
649 Standard_Real Tol2 = myTolerance*myTolerance;
650 if (VDx.SquareMagnitude() < Tol2 ||
651 VDy.SquareMagnitude() < Tol2 ||
652 VDx.CrossSquareMagnitude(VDy) < Tol2)
654 myIsApprox = Standard_True;
661 gp_Pnt O = Axis.Location();
662 gp_Pnt P = ProjectPnt(myPlane,myDirection,O);
663 gp_Pnt Px = ProjectPnt(myPlane,myDirection,O.Translated(R1*gp_Vec(X)));
664 gp_Pnt Py = ProjectPnt(myPlane,myDirection,O.Translated(R2*gp_Vec(Y)));
665 Standard_Real Major = P.Distance(Px);
666 Standard_Real Minor = P.Distance(Py);
670 myIsApprox = !gp_Dir(VDx).IsNormal(gp_Dir(VDy), Precision::Angular());
674 // Since it is not necessary to keep the same parameter for the point on the original and on the projected curves,
675 // we will use the following approach to find axes of the projected ellipse and provide the canonical curve:
676 // https://www.geometrictools.com/Documentation/ParallelProjectionEllipse.pdf
677 math_Matrix aMatrA(1, 2, 1, 2);
678 // A = Jp^T * Pr(Je), where
679 // Pr(Je) - projection of axes of original ellipse to the target plane
680 // Jp - X and Y axes of the target plane
681 aMatrA(1, 1) = myPlane.XDirection().XYZ().Dot(VDx.XYZ());
682 aMatrA(1, 2) = myPlane.XDirection().XYZ().Dot(VDy.XYZ());
683 aMatrA(2, 1) = myPlane.YDirection().XYZ().Dot(VDx.XYZ());
684 aMatrA(2, 2) = myPlane.YDirection().XYZ().Dot(VDy.XYZ());
686 math_Matrix aMatrDelta2(1, 2, 1, 2, 0.0);
687 // | 1/MajorRad^2 0 |
689 // | 0 1/MajorRad^2 |
690 aMatrDelta2(1, 1) = 1.0 / (R1 * R1);
691 aMatrDelta2(2, 2) = 1.0 / (R2 * R2);
693 math_Matrix aMatrAInv = aMatrA.Inverse();
694 math_Matrix aMatrM = aMatrAInv.Transposed() * aMatrDelta2 * aMatrAInv;
696 // perform eigenvalues calculation
697 math_Jacobi anEigenCalc(aMatrM);
698 if (anEigenCalc.IsDone())
700 // radii of the projected ellipse
701 Minor = 1.0 / Sqrt(anEigenCalc.Value(1));
702 Major = 1.0 / Sqrt(anEigenCalc.Value(2));
704 // calculate the rotation angle for the plane axes to meet the correct axes of the projected ellipse
705 // (swap eigenvectors in respect to major and minor axes)
706 const math_Matrix& anEigenVec = anEigenCalc.Vectors();
707 gp_Trsf2d aTrsfInPlane;
708 aTrsfInPlane.SetValues(anEigenVec(1, 2), anEigenVec(1, 1), 0.0,
709 anEigenVec(2, 2), anEigenVec(2, 1), 0.0);
711 aRot.SetRotation(gp_Ax1(P, myPlane.Direction()), aTrsfInPlane.RotationPart());
713 Dx = myPlane.XDirection().Transformed(aRot);
714 Dy = myPlane.YDirection().Transformed(aRot);
718 myIsApprox = Standard_True;
724 gp_Ax2 Axe(P, Dx^Dy, Dx);
726 if (Abs(Major - Minor) < Precision::Confusion()) {
727 myType = GeomAbs_Circle;
728 gp_Circ Circ(Axe, Major);
729 GeomCirclePtr = new Geom_Circle(Circ);
730 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
731 GeomAdaptor_Curve aGACurve(GeomCirclePtr);
732 myResult = new GeomAdaptor_HCurve(aGACurve);
733 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
735 else if ( Major > Minor) {
736 myType = GeomAbs_Ellipse;
737 Elips = gp_Elips( Axe, Major, Minor);
739 GeomEllipsePtr = new Geom_Ellipse(Elips);
740 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
741 GeomAdaptor_Curve aGACurve(GeomEllipsePtr);
742 myResult = new GeomAdaptor_HCurve(aGACurve);
743 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
746 myIsApprox = Standard_True;
751 // No way to build the canonical curve, approximate as B-spline
754 myType = GeomAbs_BSplineCurve;
755 PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
756 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
757 GeomAdaptor_Curve aGACurve(ApproxCurve);
758 myResult = new GeomAdaptor_HCurve(aGACurve);
759 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
761 else if (GeomCirclePtr || GeomEllipsePtr)
763 Handle(Geom_Curve) aResultCurve = GeomCirclePtr;
764 if (aResultCurve.IsNull())
765 aResultCurve = GeomEllipsePtr;
766 // start and end parameters of the projected curve
767 Standard_Real aParFirst = myCurve->FirstParameter();
768 Standard_Real aParLast = myCurve->LastParameter();
769 gp_Pnt aPntFirst = ProjectPnt(myPlane, myDirection, myCurve->Value(aParFirst));
770 gp_Pnt aPntLast = ProjectPnt(myPlane, myDirection, myCurve->Value(aParLast));
771 GeomLib_Tool::Parameter(aResultCurve, aPntFirst, Precision::Confusion(), myFirstPar);
772 GeomLib_Tool::Parameter(aResultCurve, aPntLast, Precision::Confusion(), myLastPar);
773 while (myLastPar <= myFirstPar)
774 myLastPar += myResult->Period();
778 case GeomAbs_Parabola:
780 myKeepParam = Standard_True;
781 // P(u) = O + (u*u)/(4*f) * Xc + u * Yc
782 // ==> Q(u) = f(P(u))
783 // = f(O) + (u*u)/(4*f) * f(Xc) + u * f(Yc)
785 gp_Parab Parab = myCurve->Parabola();
786 gp_Ax2 AxeRef = Parab.Position();
787 gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.XDirection()));
788 gp_Vec Yc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.YDirection()));
790 // fix for case when no one action is done. 28.03.2002
791 Standard_Boolean alocalIsDone = Standard_False;
793 if ( Abs( Yc.Magnitude() - 1.) < Precision::Confusion()) {
794 gp_Pnt P = ProjectPnt(myPlane,myDirection,AxeRef.Location());
795 if ( Xc.Magnitude() < Precision::Confusion()) {
796 myType = GeomAbs_Line;
797 aLine = gp_Lin( P, gp_Dir(Yc));
799 GeomLinePtr = new Geom_Line(aLine) ;
800 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
801 GeomAdaptor_Curve aGACurve(GeomLinePtr);
802 myResult = new GeomAdaptor_HCurve(aGACurve);
803 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
804 alocalIsDone = Standard_True;
806 else if ( Xc.IsNormal(Yc,Precision::Angular())) {
807 myType = GeomAbs_Parabola;
808 Standard_Real F = Parab.Focal() / Xc.Magnitude();
809 gp_Parab aParab = gp_Parab( gp_Ax2(P,Xc^Yc,Xc), F);
810 Handle(Geom_Parabola) GeomParabolaPtr =
811 new Geom_Parabola(aParab) ;
812 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
813 GeomAdaptor_Curve aGACurve(GeomParabolaPtr);
814 myResult = new GeomAdaptor_HCurve(aGACurve);
815 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
816 alocalIsDone = Standard_True;
819 if (!alocalIsDone)/*else*/ {
820 myIsApprox = Standard_True;
821 myType = GeomAbs_BSplineCurve;
822 PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
823 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
824 GeomAdaptor_Curve aGACurve(ApproxCurve);
825 myResult = new GeomAdaptor_HCurve(aGACurve);
826 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
830 case GeomAbs_Hyperbola:
832 myKeepParam = Standard_True;
833 // P(u) = O + R1 * Cosh(u) * Xc + R2 * Sinh(u) * Yc
834 // ==> Q(u) = f(P(u))
835 // = f(O) + R1 * Cosh(u) * f(Xc) + R2 * Sinh(u) * f(Yc)
837 gp_Hypr Hypr = myCurve->Hyperbola();
838 gp_Ax2 AxeRef = Hypr.Position();
839 gp_Vec Xc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.XDirection()));
840 gp_Vec Yc = ProjectVec(myPlane,myDirection,gp_Vec(AxeRef.YDirection()));
841 gp_Pnt P = ProjectPnt(myPlane,myDirection,AxeRef.Location());
842 Standard_Real aR1 = Hypr.MajorRadius();
843 Standard_Real aR2 = Hypr.MinorRadius();
844 gp_Dir Z = myPlane.Direction();
846 if ( Xc.Magnitude() < Precision::Confusion()) {
847 myType = GeomAbs_Hyperbola;
848 gp_Dir X = gp_Dir(Yc) ^ Z;
849 Hypr = gp_Hypr(gp_Ax2( P, Z, X), 0., aR2 * Yc.Magnitude());
851 new Geom_Hyperbola(Hypr) ;
852 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
853 GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
854 myResult = new GeomAdaptor_HCurve(aGACurve);
855 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
857 else if ( Yc.Magnitude() < Precision::Confusion()) {
858 myType = GeomAbs_Hyperbola;
860 gp_Hypr(gp_Ax2(P, Z, gp_Dir(Xc)), aR1 * Xc.Magnitude(), 0.);
862 new Geom_Hyperbola(Hypr) ;
863 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
864 GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
865 myResult = new GeomAdaptor_HCurve(aGACurve);
866 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
868 else if ( Xc.IsNormal(Yc,Precision::Angular())) {
869 myType = GeomAbs_Hyperbola;
870 Hypr = gp_Hypr( gp_Ax2( P, gp_Dir( Xc ^ Yc), gp_Dir( Xc)),
871 aR1 * Xc.Magnitude(), aR2 * Yc.Magnitude() );
873 new Geom_Hyperbola(Hypr) ;
874 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
875 GeomAdaptor_Curve aGACurve(GeomHyperbolaPtr);
876 myResult = new GeomAdaptor_HCurve(aGACurve);
877 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
880 myIsApprox = Standard_True;
881 myType = GeomAbs_BSplineCurve;
882 PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
883 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
884 GeomAdaptor_Curve aGACurve(ApproxCurve);
885 myResult = new GeomAdaptor_HCurve(aGACurve);
886 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
890 case GeomAbs_BezierCurve:
892 Handle(Geom_BezierCurve) BezierCurvePtr =
894 Standard_Integer NbPoles =
895 BezierCurvePtr->NbPoles() ;
897 Handle(Geom_BezierCurve) ProjCu =
898 Handle(Geom_BezierCurve)::DownCast(BezierCurvePtr->Copy());
900 myKeepParam = Standard_True;
901 myIsApprox = Standard_False;
903 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
905 (i,ProjectPnt(myPlane,myDirection,BezierCurvePtr->Pole(i)));
908 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
909 GeomAdaptor_Curve aGACurve(ProjCu);
910 myResult = new GeomAdaptor_HCurve(aGACurve);
911 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
914 case GeomAbs_BSplineCurve:
916 Handle(Geom_BSplineCurve) BSplineCurvePtr =
919 // make a copy of the curve and projects its poles
921 Handle(Geom_BSplineCurve) ProjectedBSplinePtr =
922 Handle(Geom_BSplineCurve)::DownCast(BSplineCurvePtr->Copy()) ;
924 myKeepParam = Standard_True;
925 myIsApprox = Standard_False;
927 for ( Standard_Integer i = 1; i <= BSplineCurvePtr->NbPoles(); i++) {
928 ProjectedBSplinePtr->SetPole
929 (i,ProjectPnt(myPlane,myDirection,BSplineCurvePtr->Pole(i)));
932 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
933 GeomAdaptor_Curve aGACurve(ProjectedBSplinePtr);
934 myResult = new GeomAdaptor_HCurve(aGACurve);
935 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
940 myKeepParam = Standard_True;
941 myIsApprox = Standard_True;
942 myType = GeomAbs_BSplineCurve;
943 PerformApprox(myCurve,myPlane,myDirection,ApproxCurve);
944 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:29 2002 Begin
945 GeomAdaptor_Curve aGACurve(ApproxCurve);
946 myResult = new GeomAdaptor_HCurve(aGACurve);
947 // Modified by Sergey KHROMOV - Tue Jan 29 16:57:30 2002 End
953 //=======================================================================
954 //function : GetPlane
956 //=======================================================================
958 const gp_Ax3& ProjLib_ProjectOnPlane::GetPlane() const
963 //=======================================================================
964 //function : GetDirection
966 //=======================================================================
968 const gp_Dir& ProjLib_ProjectOnPlane::GetDirection() const
973 //=======================================================================
974 //function : GetCurve
976 //=======================================================================
978 const Handle(Adaptor3d_HCurve)& ProjLib_ProjectOnPlane::GetCurve() const
983 //=======================================================================
984 //function : GetResult
986 //=======================================================================
988 const Handle(GeomAdaptor_HCurve)& ProjLib_ProjectOnPlane::GetResult() const
994 //=======================================================================
995 //function : FirstParameter
997 //=======================================================================
999 Standard_Real ProjLib_ProjectOnPlane::FirstParameter() const
1001 if ( myKeepParam || myIsApprox)
1002 return myCurve->FirstParameter();
1008 //=======================================================================
1009 //function : LastParameter
1011 //=======================================================================
1013 Standard_Real ProjLib_ProjectOnPlane::LastParameter() const
1015 if ( myKeepParam || myIsApprox)
1016 return myCurve->LastParameter();
1022 //=======================================================================
1023 //function : Continuity
1025 //=======================================================================
1027 GeomAbs_Shape ProjLib_ProjectOnPlane::Continuity() const
1029 return myCurve->Continuity() ;
1033 //=======================================================================
1034 //function : NbIntervals
1036 //=======================================================================
1038 Standard_Integer ProjLib_ProjectOnPlane::NbIntervals(const GeomAbs_Shape S) const
1040 return myCurve->NbIntervals(S) ;
1044 //=======================================================================
1045 //function : Intervals
1047 //=======================================================================
1049 void ProjLib_ProjectOnPlane::Intervals(TColStd_Array1OfReal& T,
1050 const GeomAbs_Shape S) const
1052 myCurve->Intervals(T,S) ;
1055 //=======================================================================
1058 //=======================================================================
1060 Handle(Adaptor3d_HCurve)
1061 ProjLib_ProjectOnPlane::Trim(const Standard_Real First,
1062 const Standard_Real Last,
1063 const Standard_Real Tolerance) const
1065 if (myType != GeomAbs_OtherCurve){
1066 return myResult->Trim(First,Last,Tolerance) ;
1069 throw Standard_NotImplemented ("ProjLib_ProjectOnPlane::Trim() - curve of unsupported type");
1074 //=======================================================================
1075 //function : IsClosed
1077 //=======================================================================
1079 Standard_Boolean ProjLib_ProjectOnPlane::IsClosed() const
1081 return myCurve->IsClosed() ;
1085 //=======================================================================
1086 //function : IsPeriodic
1088 //=======================================================================
1090 Standard_Boolean ProjLib_ProjectOnPlane::IsPeriodic() const
1093 return Standard_False;
1095 return myCurve->IsPeriodic();
1099 //=======================================================================
1102 //=======================================================================
1104 Standard_Real ProjLib_ProjectOnPlane::Period() const
1106 if ( !IsPeriodic()) {
1107 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane::Period");
1111 return Standard_False;
1113 return myCurve->Period();
1117 //=======================================================================
1120 //=======================================================================
1122 gp_Pnt ProjLib_ProjectOnPlane::Value(const Standard_Real U) const
1124 if (myType != GeomAbs_OtherCurve) {
1125 return myResult->Value(U);
1128 return OnPlane_Value(U,
1137 //=======================================================================
1140 //=======================================================================
1142 void ProjLib_ProjectOnPlane::D0(const Standard_Real U , gp_Pnt& P) const
1144 if (myType != GeomAbs_OtherCurve) {
1148 P = OnPlane_Value(U,
1156 //=======================================================================
1159 //=======================================================================
1161 void ProjLib_ProjectOnPlane::D1(const Standard_Real U,
1165 if (myType != GeomAbs_OtherCurve) {
1166 myResult->D1(U,P,V) ;
1179 //=======================================================================
1182 //=======================================================================
1184 void ProjLib_ProjectOnPlane::D2(const Standard_Real U,
1189 if (myType != GeomAbs_OtherCurve) {
1190 myResult->D2(U,P,V1,V2) ;
1204 //=======================================================================
1207 //=======================================================================
1209 void ProjLib_ProjectOnPlane::D3(const Standard_Real U,
1215 if (myType != GeomAbs_OtherCurve) {
1216 myResult->D3(U,P,V1,V2,V3) ;
1231 //=======================================================================
1234 //=======================================================================
1236 gp_Vec ProjLib_ProjectOnPlane::DN(const Standard_Real U,
1237 const Standard_Integer DerivativeRequest)
1240 if (myType != GeomAbs_OtherCurve) {
1241 return myResult->DN(U,DerivativeRequest) ;
1244 return OnPlane_DN(U,
1253 //=======================================================================
1254 //function : Resolution
1256 //=======================================================================
1258 Standard_Real ProjLib_ProjectOnPlane::Resolution
1259 (const Standard_Real Tolerance) const
1261 if (myType != GeomAbs_OtherCurve) {
1262 return myResult->Resolution(Tolerance) ;
1270 //=======================================================================
1271 //function : GetType
1273 //=======================================================================
1275 GeomAbs_CurveType ProjLib_ProjectOnPlane::GetType() const
1281 //=======================================================================
1284 //=======================================================================
1286 gp_Lin ProjLib_ProjectOnPlane::Line() const
1288 if (myType != GeomAbs_Line)
1289 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Line");
1291 return myResult->Line();
1295 //=======================================================================
1298 //=======================================================================
1300 gp_Circ ProjLib_ProjectOnPlane::Circle() const
1302 if (myType != GeomAbs_Circle)
1303 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Circle");
1305 return myResult->Circle();
1309 //=======================================================================
1310 //function : Ellipse
1312 //=======================================================================
1314 gp_Elips ProjLib_ProjectOnPlane::Ellipse() const
1316 if (myType != GeomAbs_Ellipse)
1317 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Ellipse");
1319 return myResult->Ellipse();
1323 //=======================================================================
1324 //function : Hyperbola
1326 //=======================================================================
1328 gp_Hypr ProjLib_ProjectOnPlane::Hyperbola() const
1330 if (myType != GeomAbs_Hyperbola)
1331 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Hyperbola");
1333 return myResult->Hyperbola() ;
1337 //=======================================================================
1338 //function : Parabola
1340 //=======================================================================
1342 gp_Parab ProjLib_ProjectOnPlane::Parabola() const
1344 if (myType != GeomAbs_Parabola)
1345 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Parabola");
1347 return myResult->Parabola() ;
1350 //=======================================================================
1353 //=======================================================================
1355 Standard_Integer ProjLib_ProjectOnPlane::Degree() const
1357 if ((GetType() != GeomAbs_BSplineCurve) &&
1358 (GetType() != GeomAbs_BezierCurve))
1359 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Degree");
1362 return myResult->Degree();
1364 return myCurve->Degree();
1367 //=======================================================================
1368 //function : IsRational
1370 //=======================================================================
1372 Standard_Boolean ProjLib_ProjectOnPlane::IsRational() const
1374 if ((GetType() != GeomAbs_BSplineCurve) &&
1375 (GetType() != GeomAbs_BezierCurve))
1376 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:IsRational");
1379 return myResult->IsRational();
1381 return myCurve->IsRational();
1384 //=======================================================================
1385 //function : NbPoles
1387 //=======================================================================
1389 Standard_Integer ProjLib_ProjectOnPlane::NbPoles() const
1391 if ((GetType() != GeomAbs_BSplineCurve) &&
1392 (GetType() != GeomAbs_BezierCurve))
1393 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:NbPoles");
1396 return myResult->NbPoles();
1398 return myCurve->NbPoles();
1401 //=======================================================================
1402 //function : NbKnots
1404 //=======================================================================
1406 Standard_Integer ProjLib_ProjectOnPlane::NbKnots() const
1408 if ( GetType() != GeomAbs_BSplineCurve)
1409 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:NbKnots");
1412 return myResult->NbKnots();
1414 return myCurve->NbKnots();
1418 //=======================================================================
1421 //=======================================================================
1423 Handle(Geom_BezierCurve) ProjLib_ProjectOnPlane::Bezier() const
1425 if (myType != GeomAbs_BezierCurve)
1426 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:Bezier");
1428 return myResult->Bezier() ;
1431 //=======================================================================
1434 //=======================================================================
1436 Handle(Geom_BSplineCurve) ProjLib_ProjectOnPlane::BSpline() const
1438 if (myType != GeomAbs_BSplineCurve)
1439 throw Standard_NoSuchObject("ProjLib_ProjectOnPlane:BSpline");
1441 return myResult->BSpline() ;