1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
3 // This file is part of Open CASCADE Technology software library.
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
14 //:k1 abv 16.12.98 K4L PRO10107, PRO10108, PRO10109
15 //:j8 abv 10.12.98 TR10 r0501_db.stp #9423
16 //:S4030 abv, pdn: new methods - interface to standard ProjLib_CompProjectedCurve
17 //%12 pdn 15.02.99 PRO9234 optimizing
18 //%12 pdn 15.02.99 PRO9234 using improved ProjectDegenerated method
19 // rln 03.03.99 S4135: bm2_sd_t4-A.stp treatment of Geom_SphericalSurface together with V-closed surfaces
20 //:p9 abv 11.03.99 PRO7226 #489490: make IsAnIsoparametric to find nearest case
21 //:q1 abv 15.03.99 (pdn) PRO7226 #525030: limit NextValueOfUV() by tolerance
22 //:q5 abv 19.03.99 code improvement
23 //:q9 abv 23.03.99 PRO7226.stp #489490: cashe for projecting end points
24 //#78 rln 12.03.99 S4135: checking spatial closure with myPreci
25 // pdn 12.03.99 S4135: creating pcurve with minimal length in the case of densed points
26 // abv 29.03.99 IsAnIsoparametric with Precision::Confusion
27 // pdn 09.04.99 IsAnisoparametric uses already computed parameters (S4030, fix PRO14323)
29 //:s5 abv 22.04.99 Adding debug printouts in catch {} blocks
30 //#1 svv 11.01.00 Porting on DEC
32 #include <Approx_CurveOnSurface.hxx>
33 #include <Geom2d_BSplineCurve.hxx>
34 #include <Geom2d_Circle.hxx>
35 #include <Geom2d_Curve.hxx>
36 #include <Geom2d_Ellipse.hxx>
37 #include <Geom2d_Hyperbola.hxx>
38 #include <Geom2d_Line.hxx>
39 #include <Geom2d_Parabola.hxx>
40 #include <Geom2d_TrimmedCurve.hxx>
41 #include <Geom2dAdaptor.hxx>
42 #include <Geom2dAPI_Interpolate.hxx>
43 #include <Geom_BezierSurface.hxx>
44 #include <Geom_BoundedCurve.hxx>
45 #include <Geom_BSplineCurve.hxx>
46 #include <Geom_BSplineSurface.hxx>
47 #include <Geom_Curve.hxx>
48 #include <Geom_OffsetSurface.hxx>
49 #include <Geom_Plane.hxx>
50 #include <Geom_RectangularTrimmedSurface.hxx>
51 #include <Geom_SphericalSurface.hxx>
52 #include <Geom_Surface.hxx>
53 #include <Geom_SurfaceOfLinearExtrusion.hxx>
54 #include <Geom_TrimmedCurve.hxx>
55 #include <GeomAdaptor_HCurve.hxx>
56 #include <GeomAdaptor_HSurface.hxx>
57 #include <GeomAPI_Interpolate.hxx>
58 #include <GeomAPI_PointsToBSpline.hxx>
59 #include <GeomProjLib.hxx>
60 #include <gp_Pnt2d.hxx>
61 #include <NCollection_Sequence.hxx>
62 #include <Precision.hxx>
63 #include <ProjLib_CompProjectedCurve.hxx>
64 #include <ProjLib_HCompProjectedCurve.hxx>
65 #include <ProjLib_ProjectedCurve.hxx>
66 #include <ShapeAnalysis.hxx>
67 #include <ShapeAnalysis_Curve.hxx>
68 #include <ShapeAnalysis_Surface.hxx>
69 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
70 #include <ShapeExtend.hxx>
71 #include <Standard_ErrorHandler.hxx>
72 #include <Standard_Failure.hxx>
73 #include <Standard_Type.hxx>
74 #include <TColgp_Array1OfPnt.hxx>
75 #include <TColStd_Array1OfInteger.hxx>
78 IMPLEMENT_STANDARD_RTTIEXT(ShapeConstruct_ProjectCurveOnSurface,MMgt_TShared)
82 //=======================================================================
83 //function : ShapeConstruct_ProjectCurveOnSurface
85 //=======================================================================
87 ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
89 myPreci = Precision::Confusion();
90 myBuild = Standard_False;
91 myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
95 //=======================================================================
98 //=======================================================================
100 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
102 Init (new ShapeAnalysis_Surface (surf), preci);
105 //=======================================================================
108 //=======================================================================
110 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
113 SetPrecision (preci);
116 //=======================================================================
117 //function : SetSurface
119 //=======================================================================
121 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
123 SetSurface (new ShapeAnalysis_Surface (surf));
126 //=======================================================================
127 //function : SetSurface
129 //=======================================================================
131 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
133 if ( mySurf == surf ) return;
138 //=======================================================================
139 //function : SetPrecision
141 //=======================================================================
143 void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
148 //=======================================================================
149 //function : BuildCurveMode
151 //=======================================================================
153 Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
158 //=======================================================================
159 //function : AdjustOverDegenMode
161 //=======================================================================
164 //szv#4:S4163:12Mar99 was Boolean
165 Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
167 return myAdjustOverDegen;
171 //=======================================================================
172 //function : NbSurfIntervals
173 //purpose : work-around of bug in standard method
174 // GeomAdaptor_Surface->NbIntervals() (PRO16346)
175 //=======================================================================
177 static Standard_Integer NbSurfIntervals(const Handle(GeomAdaptor_HSurface)& GAS, const GeomAbs_Shape cont)
179 Standard_Integer NbU = 0;
180 if (GAS->GetType() == GeomAbs_SurfaceOfExtrusion) {
181 // extract the surface
182 Handle(Geom_SurfaceOfLinearExtrusion) surf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(GAS->ChangeSurface().Surface());
183 // build a 3d adaptor curve
184 GeomAdaptor_Curve Adaptor3dCurve(surf->BasisCurve(), GAS->FirstUParameter(), GAS->LastUParameter());
185 if (Adaptor3dCurve.GetType() == GeomAbs_BSplineCurve)
186 NbU = Adaptor3dCurve.NbIntervals(cont);
189 NbU = GAS->NbUIntervals(cont);
190 return NbU * (GAS->NbVIntervals(cont));
193 //=======================================================================
196 //=======================================================================
198 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status Status) const
200 return ShapeExtend::DecodeStatus (myStatus, Status);
203 //=======================================================================
206 //=======================================================================
207 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
208 const Standard_Real First,
209 const Standard_Real Last,
210 Handle(Geom2d_Curve)& c2d,
212 const Standard_Integer,
213 const Standard_Integer)
215 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
216 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
218 if (mySurf.IsNull()) {
220 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
221 return Standard_False;
223 // Projection Analytique
224 Handle(Geom_Curve) crv3dtrim = c3d;
225 if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
226 crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
227 c2d = ProjectAnalytic ( crv3dtrim );
229 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
230 return Standard_True;
233 // Projection par approximation
235 // discretize the 3d curve
237 Standard_Integer nbrPnt;
239 // $$$$ :92 abv 28 Jan 98 see PRO10107, big BSplineCurve C0
240 Standard_Integer nbPini = NCONTROL; // as in BRepCheck_Edge (RLN/Nijni)
241 // 20; // number of points for interpolation, should be "parametric dependent"
243 //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
244 // increase number of points to provide at least Degree()+1 points per interval
245 Handle(Geom_BSplineCurve) bspl;
246 if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
247 Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
248 bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
250 else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
251 if ( ! bspl.IsNull() ) {
252 Standard_Integer nint = 0;
253 for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
254 if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last ) nint++;
255 Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
256 while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
258 if ( nbPini > NCONTROL )
259 cout << "Warning: number of points for projecting is " << nbPini << endl;
263 // $$$$ end :92 (big BSplineCurve C0)
265 // this number should be "parametric dependent"
266 TColgp_Array1OfPnt points(1, nbPini);
267 TColStd_Array1OfReal params(1, nbPini);
268 NCollection_Sequence<Standard_Real> aKnotCoeffs;
270 Standard_Integer iPnt;
272 // In case of bspline compute parametrization speed on each
273 // knot interval inside [aFirstParam, aLastParam].
274 // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
275 // use PerformByProjLib algorithm.
278 Standard_Real aFirstParam = First; // First parameter of current interval.
279 Standard_Real aLastParam = Last; // Last parameter of current interval.
281 // First index computation.
282 Standard_Integer anIdx = 1;
283 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
285 if(bspl->Knot(anIdx) > First)
291 GeomAdaptor_Curve aC3DAdaptor(c3d);
292 Standard_Real aMinParSpeed = Precision::Infinite(); // Minimal parameterization speed.
293 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
295 // Fill current knot interval.
296 aLastParam = Min(Last, bspl->Knot(anIdx));
297 Standard_Integer aNbIntPnts = NCONTROL;
298 // Number of inner points is adapted according to the length of the interval
299 // to avoid a lot of calculations on small range of parameters.
302 const Standard_Real aLenThres = 1.e-2;
303 const Standard_Real aLenRatio =
304 (aLastParam - aFirstParam) / (bspl->Knot(anIdx) - bspl->Knot(anIdx - 1));
305 if (aLenRatio < aLenThres)
307 aNbIntPnts = Standard_Integer(aLenRatio / aLenThres * aNbIntPnts);
312 Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
313 Standard_Integer anIntIdx;
315 // Start filling from first point.
316 aC3DAdaptor.D0(aFirstParam, p3d1);
318 Standard_Real aLength3d = 0.0;
319 for(anIntIdx = 1; anIntIdx < aNbIntPnts; anIntIdx++)
321 Standard_Real aParam = aFirstParam + aStep * anIntIdx;
322 aC3DAdaptor.D0 (aParam, p3d2);
323 const Standard_Real aDist = p3d2.Distance(p3d1);
328 aMinParSpeed = Min(aMinParSpeed, aDist / aStep);
330 const Standard_Real aCoeff = aLength3d / (aLastParam - aFirstParam);
331 if (Abs(aCoeff) > gp::Resolution())
332 aKnotCoeffs.Append(aCoeff);
333 aFirstParam = aLastParam;
336 Standard_Real anEvenlyCoeff = 0;
337 if (aKnotCoeffs.Size() > 0)
339 anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
340 *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
343 const Standard_Real aMaxQuotientCoeff = 1500.0;
344 if (anEvenlyCoeff > aMaxQuotientCoeff &&
345 aMinParSpeed > Precision::Confusion() )
347 PerformByProjLib(c3d, First, Last, c2d);
348 // PerformByProjLib fail detection:
351 return Status (ShapeExtend_DONE);
356 Standard_Real deltaT, t;
357 deltaT = (Last - First) / (nbPini-1);
359 for (iPnt = 1; iPnt <= nbPini; iPnt ++)
361 if (iPnt == 1) t = First;
362 else if (iPnt == nbPini) t = Last;
363 else t = First + (iPnt - 1) * deltaT;
370 // CALCUL par approximation
371 TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
372 ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
374 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
375 return Standard_True;
376 }// cas particulier d iso
378 // INTERPOLATION du resultat
381 Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt (1, nbPini);
382 Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
383 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
384 thePnts->SetValue(iPnt, points(iPnt));
385 theParams->SetValue(iPnt, params(iPnt));
388 Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
389 if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
391 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
396 Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
397 Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
398 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
399 theParams2d->SetValue(iPnt, params(iPnt));
400 thePnts2d->SetValue(iPnt, pnt2d(iPnt));
403 c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
404 // c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
405 // Faut-il aussi reprendre la C3D ?
406 myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
407 return Status (ShapeExtend_DONE);
410 //=======================================================================
411 //function : PerformByProjLib
413 //=======================================================================
415 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
416 const Standard_Real First,
417 const Standard_Real Last,
418 Handle(Geom2d_Curve)& c2d,
419 const GeomAbs_Shape /*continuity*/,
420 const Standard_Integer /*maxdeg */,
421 const Standard_Integer /*nbinterval */)
423 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
425 if (mySurf.IsNull()) {
426 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
427 return Standard_False;
433 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
434 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
435 ProjLib_ProjectedCurve Projector(GAS, GAC);
437 switch (Projector.GetType())
440 c2d = new Geom2d_Line(Projector.Line());
442 case GeomAbs_Circle :
443 c2d = new Geom2d_Circle(Projector.Circle());
445 case GeomAbs_Ellipse :
446 c2d = new Geom2d_Ellipse(Projector.Ellipse());
448 case GeomAbs_Parabola :
449 c2d = new Geom2d_Parabola(Projector.Parabola());
451 case GeomAbs_Hyperbola :
452 c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
454 case GeomAbs_BSplineCurve :
455 c2d = Projector.BSpline();
458 // Not possible, handling added to avoid gcc warning.
464 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
465 return Standard_False;
469 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
470 return Standard_True;
474 catch(Standard_Failure const& anException) {
476 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
477 anException.Print(cout); cout << endl;
480 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
483 return Standard_False;
486 //=======================================================================
487 //function : PerformAdvanced
489 //=======================================================================
491 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
492 const Standard_Real First,
493 const Standard_Real Last,
494 Handle(Geom2d_Curve)& c2d)
496 Standard_Boolean hasResult = Standard_False;
497 Standard_Integer nbintervals;
499 Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
500 // && (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
502 if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
504 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
505 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
506 nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
507 isStandard = (nbintervals < 2);
510 hasResult = PerformByProjLib(c3d, First, Last, c2d);
512 if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
516 //=======================================================================
517 //function : ProjectAnalytic
519 //=======================================================================
521 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
523 Handle(Geom2d_Curve) result;
525 //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
526 // This is necessary for K4L since it fails on other surfaces
527 // when general method GeomProjLib::Curve2d() is used
528 // Projection is done as in BRep_Tool and BRepCheck_Edge
529 Handle(Geom_Surface) surf = mySurf->Surface();
530 Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
531 if ( Plane.IsNull() ) {
532 Handle(Geom_RectangularTrimmedSurface) RTS =
533 Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
534 if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
536 Handle(Geom_OffsetSurface) OS =
537 Handle(Geom_OffsetSurface)::DownCast ( surf );
539 Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
542 if ( ! Plane.IsNull() ) {
543 Handle(Geom_Curve) ProjOnPlane =
544 GeomProjLib::ProjectOnPlane (c3d, Plane,
545 Plane->Position().Direction(), Standard_True);
546 Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
547 ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
549 result = Geom2dAdaptor::MakeCurve(Proj);
550 if ( result.IsNull() ) return result;
551 if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
552 Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
553 result = TC->BasisCurve();
562 //! Fix possible period jump and handle walking period parameter.
563 static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
564 Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
565 Standard_Real thePeriod, // Period on objective coord
566 Standard_Integer theSavedPoint, // Point number to choose period
567 Standard_Real theSavedParam) // Param from cashe to choose period
569 Standard_Real aSavedParam;
570 Standard_Integer aSavedPoint;
571 Standard_Real aMinParam = 0.0, aMaxParam = thePeriod;
572 if (theSavedPoint < 0) {
573 // normalize to first period by default
574 aSavedParam = 0.5 * thePeriod;
578 aSavedParam = theSavedParam;
579 aSavedPoint = theSavedPoint;
580 while (aMinParam > aSavedParam) {
581 aMinParam -= thePeriod;
582 aMaxParam -= thePeriod;
584 while (aMaxParam < aSavedParam) {
585 aMinParam += thePeriod;
586 aMaxParam += thePeriod;
590 Standard_Real aFixIsoParam = aMinParam;
591 Standard_Boolean isIsoLine = Standard_False;
592 if (aMaxParam - aSavedParam < Precision::PConfusion() ||
593 aSavedParam - aMinParam < Precision::PConfusion()) {
594 aFixIsoParam = aSavedParam;
595 isIsoLine = Standard_True;
597 // normalize all coordinates to [aMinParam, aMaxParam)
598 for (Standard_Integer i = 0; i < 4; i++) {
599 Standard_Real aParam = thePnt[i].Coord(theIdx);
600 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(aParam, aMinParam, aMaxParam);
602 // Walk over period coord -> not walking on another isoline in parameter space.
604 if (aMaxParam - aParam < Precision::PConfusion() || aParam - aMinParam < Precision::PConfusion())
605 aParam = aFixIsoParam;
608 if (aMaxParam - aParam < Precision::PConfusion())
610 if (aParam - aMinParam < Precision::PConfusion())
614 thePnt[i].SetCoord(theIdx, aParam);
617 // find possible period jump and increasing or decreasing coordinates vector
618 Standard_Boolean isJump = Standard_False;
619 Standard_Real aPrevDiff = 0.0;
620 Standard_Real aSumDiff = 1.0;
621 for (Standard_Integer i = 0; i < 3; i++) {
622 Standard_Real aDiff = thePnt[i + 1].Coord(theIdx) - thePnt[i].Coord(theIdx);
623 if (aDiff < -Precision::PConfusion()) {
626 //if first derivative changes its sign then period jump may exists in this place
627 if (aDiff * aPrevDiff < -Precision::PConfusion()) {
628 isJump = Standard_True;
634 return Standard_False;
636 if (aSumDiff > 0) { // decreasing sequence (parameters decrease twice(--) and one period jump(+))
637 for (Standard_Integer i = aSavedPoint; i > 0; i--)
638 if (thePnt[i].Coord(theIdx) > thePnt[i - 1].Coord(theIdx)) {
639 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) + thePeriod);
641 for (Standard_Integer i = aSavedPoint; i < 3; i++)
642 if (thePnt[i].Coord(theIdx) < thePnt[i + 1].Coord(theIdx)) {
643 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) - thePeriod);
646 else {// increasing sequence (parameters increase twice(++) and one period jump(-))
647 for (Standard_Integer i = aSavedPoint; i > 0; i--)
648 if (thePnt[i].Coord(theIdx) < thePnt[i - 1].Coord(theIdx)) {
649 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) - thePeriod);
651 for (Standard_Integer i = aSavedPoint; i < 3; i++)
652 if (thePnt[i].Coord(theIdx) > thePnt[i + 1].Coord(theIdx)) {
653 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) + thePeriod);
657 // Do not return false, because for nonlinear 2d curves vector of parameters
658 // may change its first derivative and shifted parameters will be broken for this case.
659 return Standard_True;
662 //=======================================================================
665 //=======================================================================
667 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
668 const TColgp_Array1OfPnt& thepoints,
669 const TColStd_Array1OfReal& theparams,
670 TColgp_Array1OfPnt2d& thePnt2ds,
671 Standard_Real theTol,
672 Standard_Boolean &isRecompute,
673 Standard_Boolean &isFromCashe) const
675 Standard_Integer nb = thepoints.Length();
677 aP[0] = thepoints(1);
678 aP[1] = thepoints(2);
679 aP[2] = thepoints(nb - 1);
680 aP[3] = thepoints(nb);
682 Standard_Integer i = 0;
684 Standard_Real aTol2 = theTol * theTol;
685 Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
686 Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
689 // Protection against bad "tolerance" shapes.
692 theTol = Precision::Confusion();
693 aTol2 = theTol * theTol;
695 if (aTol2 < Precision::SquareConfusion())
696 aTol2 = Precision::SquareConfusion();
697 Standard_Real anOldTol2 = aTol2;
698 // auxiliary variables to choose period for connection with previous 2dcurve (if exist)
699 Standard_Integer aSavedPointNum = -1;
700 gp_Pnt2d aSavedPoint;
702 // project first and last points
706 for ( j=0; j < myNbCashe; j++ )
707 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
709 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol,
712 aSavedPoint = myCashe2d[j];
714 isFromCashe = Standard_True;
717 if ( j >= myNbCashe )
718 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
720 Standard_Real aDist = mySurf->Gap();
721 Standard_Real aCurDist = aDist * aDist;
722 if( aTol2 < aDist * aDist)
726 if ( isPeriodicU || isPeriodicV )
728 // Compute second and last but one c2d points.
729 for(i = 1; i < 3; i++)
732 for ( j=0; j < myNbCashe; j++ )
733 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
735 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
737 aSavedPoint = myCashe2d[j];
740 if ( j >= myNbCashe )
741 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
743 Standard_Real aDist = mySurf->Gap();
744 Standard_Real aCurDist = aDist * aDist;
745 if( aTol2 < aDist * aDist)
751 isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
756 isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
760 if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
761 // Do not try to make line, because in this case may be very special case when 3d curve
762 // go over the pole of, e.g., sphere, and partly lies along seam
763 // (see ApproxPCurve() for more information).
767 thePnt2ds.SetValue(1, aP2d[0]);
768 thePnt2ds.SetValue(nb, aP2d[3]);
770 // Restore old tolerance in 2d space to avoid big gap cases.
772 // Check that straight line in 2d with parameterisation as in 3d will fit
773 // fit 3d curve at all points.
774 Standard_Real dPar = theparams(nb) - theparams(1);
775 if ( Abs(dPar) < Precision::PConfusion() )
777 gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
778 gp_Vec2d aVec = aVec0 / dPar;
779 Handle(Geom_Surface) aSurf = mySurf->Surface();
780 Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
782 for(i = 1; i <= nb; i++)
784 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
786 gp_Vec aNormalVec, aDu, aDv;
787 aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
788 aNormalVec = aDu ^ aDv;
789 if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
790 isNormalCheck = Standard_False;
793 gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
794 Standard_Real aDist = aNormalLine.Distance(thepoints(i));
799 if (!isNormalCheck) {
800 Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
801 SquareDistance(thepoints(1));
802 aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
803 for(i = 2; i < nb; i++)
805 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
807 aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
808 Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
810 if(Abs (aFirstPointDist - aDist1) > aTol2)
815 // check if pcurve can be represented by Geom2d_Line (parameterised by length)
816 Standard_Real aLLength = aVec0.Magnitude();
817 if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
819 gp_XY aDirL = aVec0.XY() / aLLength;
820 gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
821 return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
824 // create straight bspline
825 TColgp_Array1OfPnt2d aPoles(1, 2);
829 TColStd_Array1OfReal aKnots(1,2);
830 aKnots(1) = theparams(1);
831 aKnots(2) = theparams(theparams.Length());
833 TColStd_Array1OfInteger aMults(1,2);
836 Standard_Integer aDegree = 1;
837 Handle(Geom2d_BSplineCurve) abspl2d =
838 new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
842 //=======================================================================
843 //function : ApproxPCurve
845 //=======================================================================
847 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
848 const TColgp_Array1OfPnt& points,
849 const TColStd_Array1OfReal& params,
850 TColgp_Array1OfPnt2d& pnt2d,
851 Handle(Geom2d_Curve)& c2d)
853 // for performance, first try to handle typical case when pcurve is straight
854 Standard_Boolean isRecompute = Standard_False;
855 Standard_Boolean isFromCasheLine = Standard_False;
856 c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
860 Standard_Boolean ChangeCycle = Standard_False;
861 if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
862 myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
863 ChangeCycle = Standard_True;
866 myCashe3d[0] = points(1);
867 myCashe3d[1] = points(nbrPnt);
868 myCashe2d[0] = pnt2d(1);
869 myCashe2d[1] = pnt2d(nbrPnt);
872 myCashe3d[1] = points(1);
873 myCashe3d[0] = points(nbrPnt);
874 myCashe2d[1] = pnt2d(1);
875 myCashe2d[0] = pnt2d(nbrPnt);
877 return Standard_True;
879 Standard_Boolean isDone = Standard_True;
880 // test if the curve 3d is a boundary of the surface
881 // (only for Bezier or BSpline surface)
883 Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
884 gp_Pnt2d valueP1, valueP2;
885 Handle(Geom_Curve) cIso;
886 Standard_Real t1, t2;
888 Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
889 Standard_Boolean isAnalytic = Standard_True;
890 if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
891 Standard_Real uf, ul, vf, vl;
892 mySurf->Surface()->Bounds(uf, ul, vf, vl);
893 isoclosed = Standard_False;
894 TColStd_Array1OfReal pout(1, nbrPnt);
896 isoParam = IsAnIsoparametric(nbrPnt, points, params,
897 isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
898 isoPar2d3d, cIso, t1, t2, pout);
900 // projection of the points on surfaces
905 Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
906 Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
908 // Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
909 // Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
912 isoValue = valueP1.X();
913 isoPar1 = valueP1.Y();
914 isoPar2 = valueP2.Y();
915 isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
916 parf = vf; parl = vl;
919 isoValue = valueP1.Y();
920 isoPar1 = valueP1.X();
921 isoPar2 = valueP2.X();
922 isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
923 parf = uf; parl = ul;
925 if (!isoPar2d3d && !isAnalytic) {
926 Cf = cIso->FirstParameter();
927 Cl = cIso->LastParameter();
928 if (Precision::IsInfinite(Cf)) Cf = -1000;
929 if (Precision::IsInfinite(Cl)) Cl = +1000;
930 //pdn S4030 optimizing and fix isopar case on PRO41323
932 // dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
933 // Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
934 if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
935 if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
937 if (isoTypeU) valueP1.SetY (isoPar1);
938 else valueP1.SetX (isoPar1);
940 if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
941 //pdn S4030 optimizing and fix isopar case on PRO41323
942 tfin = pout(nbrPnt-1);
943 //dist = ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
944 if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
946 if (isoTypeU) valueP2.SetY (isoPar2);
947 else valueP2.SetX (isoPar2);
950 // Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
951 // Est-ce encore necessaire apres ce qui vient d etre fait ?
953 // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
954 // if curve is periodic do not change the points
955 // skl change "if" for pout(nbrPnt-1) 19.11.2003
957 if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
958 (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
959 gp_Pnt2d valueTmp = valueP1;
960 valueP1 = valueP2; valueP2 = valueTmp;
962 isoValue = valueP1.X();
963 isoPar1 = valueP1.Y();
964 isoPar2 = valueP2.Y();
967 isoValue = valueP1.Y();
968 isoPar1 = valueP1.X();
969 isoPar2 = valueP2.X();
971 // Fin calcul sens de courbe iso
973 } // end of fix check 05.02.02
977 // Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
978 // (utile pour detromper bug de projection). Mais detromper aussi circularite
980 //if (p1OnIso) valueP1 =
981 //BestExtremum (valueP1,points(1),points(2));
982 //if (p2OnIso) valueP2 =
983 //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
986 Standard_Real gap = myPreci; //:q1
987 Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
988 // auxiliaruy variables to shift 2dcurve, according to previous
989 Standard_Boolean isFromCashe = Standard_False;
990 gp_Pnt2d aSavedPoint;
991 if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
992 //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
993 if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
994 ChangeCycle = Standard_True;
995 //for( i = 1; i <= nbrPnt; i ++) {
996 for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
997 if(ChangeCycle) //skl for OCC3430
1005 if (isoPar2 > isoPar1) tPar = params(i);
1006 else tPar = t1 + t2 - params(i);
1007 } else if (!isAnalytic) {
1008 // projection to iso
1009 if (i==1) tPar = isoPar1;
1010 else if (i==nbrPnt) tPar = isoPar2;
1013 //:S4030 ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
1017 if (!isoPar2d3d && isAnalytic) {
1018 if (i == 1) p2d = valueP1;
1019 else if (i == nbrPnt) p2d = valueP2;
1021 p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
1022 Precision::Confusion()+1000*gap); //:q1
1023 gap = mySurf->Gap();
1026 if(isoTypeU) { p2d.SetX(isoValue); p2d.SetY(tPar); }
1027 else { p2d.SetX(tPar); p2d.SetY(isoValue); }
1032 if ( (i == 1) && p1OnIso) p2d = valueP1;
1033 else if( (i == nbrPnt) && p2OnIso) p2d = valueP2;
1034 else {// general case (not an iso) mais attention aux singularites !
1035 // first and last points are already computed by getLine()
1036 if ( (i == 1 || i == nbrPnt))
1041 gap = mySurf->Gap();
1043 isFromCashe = isFromCasheLine;
1050 //:q9 abv 23 Mar 99: use cashe as 1st approach
1051 Standard_Integer j; // svv #1
1052 for ( j=0; j < myNbCashe; j++ )
1053 if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
1055 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci,
1056 Precision::Confusion()+gap);
1058 isFromCashe = Standard_True;
1059 aSavedPoint = myCashe2d[j];
1063 if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci);
1067 p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
1068 Precision::Confusion()+1000*gap); //:q1
1070 gap = mySurf->Gap();
1076 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
1078 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
1082 //pdn %12 11.02.99 PRO9234 entity 15402
1084 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
1085 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
1088 // attention aux singularites ... (hors cas iso qui les traite deja)
1091 // if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
1093 // p2d = pnt2d (nbrPnt);
1094 // if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
1095 // pnt2d (nbrPnt) = p2d;
1098 // Si la surface est UCLosed et VClosed, on recadre les points
1099 // algo un peu complique, on retarde l implementation
1100 Standard_Real Up = ul - uf;
1101 Standard_Real Vp = vl - vf;
1102 Standard_Real dist2d;
1104 if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
1105 cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
1108 // Si la surface est UCLosed, on recadre les points
1109 if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
1110 // Premier point dans le domain [uf, ul]
1111 Standard_Real prevX, firstX = pnt2d (1).X();
1113 // do not shift 2dcurve, if it connects to previous
1114 while (firstX < uf) { firstX += Up; pnt2d (1).SetX(firstX); }
1115 while (firstX > ul) { firstX -= Up; pnt2d (1).SetX(firstX); }
1117 // shift first point, according to cashe
1118 if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
1119 Standard_Real aMinParam = uf, aMaxParam = ul;
1120 while (aMinParam > aSavedPoint.X()) {
1124 while (aMaxParam < aSavedPoint.X()) {
1128 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
1130 pnt2d(1).SetX(firstX);
1134 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1135 Standard_Real minX = firstX, maxX = firstX;
1137 // On decalle toujours le suivant
1138 for (i = 2; i <= nbrPnt; i++) {
1139 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
1140 Standard_Real CurX = pnt2d (i).X();
1141 dist2d = Abs (CurX - prevX);
1142 if (dist2d > ( Up / 2) ) {
1143 if (CurX > prevX + Up/2) {
1144 while (CurX > prevX + Up/2) { CurX -= Up; pnt2d (i).SetX (CurX); }
1145 } else if (CurX < prevX - Up/2) {
1146 while (CurX < prevX - Up/2) { CurX += Up; pnt2d (i).SetX (CurX); }
1151 if ( minX > CurX ) minX = CurX; //:97
1152 else if ( maxX < CurX ) maxX = CurX; //:97
1157 // do not shift 2dcurve, if it connects to previous
1158 Standard_Real midX = 0.5 * ( minX + maxX );
1159 Standard_Real shiftX=0.;
1160 if ( midX > ul ) shiftX = -Up;
1161 else if ( midX < uf ) shiftX = Up;
1163 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
1166 // Si la surface est VCLosed, on recadre les points
1167 // Same code as UClosed : optimisation souhaitable !!
1168 // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
1169 // Si la surface est UCLosed, on recadre les points
1171 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1172 //#78 rln 12.03.99 S4135
1173 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1174 // Premier point dans le domain [vf, vl]
1175 Standard_Real prevY, firstY = pnt2d (1).Y();
1177 // do not shift 2dcurve, if it connects to previous
1178 while (firstY < vf) { firstY += Vp; pnt2d (1).SetY(firstY); }
1179 while (firstY > vl) { firstY -= Vp; pnt2d (1).SetY(firstY); }
1181 // shift first point, according to cashe
1182 if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
1183 Standard_Real aMinParam = vf, aMaxParam = vl;
1184 while (aMinParam > aSavedPoint.Y()) {
1188 while (aMaxParam < aSavedPoint.Y()) {
1192 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
1194 pnt2d(1).SetY(firstY);
1198 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1199 Standard_Real minY = firstY, maxY = firstY;
1201 // On decalle toujours le suivant
1202 for (i = 2; i <= nbrPnt; i ++) {
1203 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
1204 Standard_Real CurY = pnt2d (i).Y();
1205 dist2d = Abs (CurY - prevY);
1206 if (dist2d > ( Vp / 2) ) {
1207 if (CurY > prevY + Vp/2) {
1208 while (CurY > prevY + Vp/2) { CurY -= Vp; pnt2d (i).SetY (CurY); }
1209 } else if (CurY < prevY - Vp/2) {
1210 while (CurY < prevY - Vp/2) { CurY += Vp; pnt2d (i).SetY (CurY); }
1214 if ( minY > CurY ) minY = CurY; //:97
1215 else if ( maxY < CurY ) maxY = CurY; //:97
1220 // do not shift 2dcurve, if it connects to previous
1221 Standard_Real midY = 0.5 * ( minY + maxY );
1222 Standard_Real shiftY=0.;
1223 if ( midY > vl ) shiftY = -Vp;
1224 else if ( midY < vf ) shiftY = Vp;
1226 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
1230 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1231 //#78 rln 12.03.99 S4135
1232 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1233 for (i = 2; i <= nbrPnt; i++) {
1234 //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1235 dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
1236 if (dist2d > ( Vp / 2) ) {
1237 // ATTENTION : il faut regarder ou le decalage se fait.
1238 // si plusieurs points sont decalles, il faut plusieurs passes
1239 // pour obtenir un resultat correct.
1240 // NOT YET IMPLEMENTED
1242 // one of those point is incorrectly placed
1243 // i.e on the wrong side of the "seam"
1244 // on prend le point le plus pres des bords vf ou vl
1245 Standard_Boolean prevOnFirst = Standard_False;
1246 Standard_Boolean prevOnLast = Standard_False;
1247 Standard_Boolean currOnFirst = Standard_False;
1248 Standard_Boolean currOnLast = Standard_False;
1250 // .X ? plutot .Y , non ?
1251 Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
1252 Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
1253 Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
1254 Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
1256 Standard_Real theMin = distPrevVF;
1257 prevOnFirst = Standard_True;
1258 if (distPrevVL < theMin) {
1259 theMin = distPrevVL;
1260 prevOnFirst = Standard_False;
1261 prevOnLast = Standard_True;
1263 if (distCurrVF < theMin) {
1264 theMin = distCurrVF;
1265 prevOnFirst = Standard_False;
1266 prevOnLast = Standard_False;
1267 currOnFirst = Standard_True;
1269 if (distCurrVL < theMin) {
1270 theMin = distCurrVL;
1271 prevOnFirst = Standard_False;
1272 prevOnLast = Standard_False;
1273 currOnFirst = Standard_False;
1274 currOnLast = Standard_True;
1276 // Modifs RLN/Nijni 3-DEC-1997
1278 // on decalle le point (i-1) en V Last
1279 gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of vl RLN/Nijni
1280 pnt2d (i-1) = newPrev;
1282 else if (prevOnLast) {
1283 // on decalle le point (i-1) en V first
1284 gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of vf RLN/Nijni
1285 pnt2d (i-1) = newPrev;
1287 else if (currOnFirst) {
1288 // on decalle le point (i) en V Last
1289 gp_Pnt2d newCurr(pnt2d (i).X(),vf); // instead of vl RLN/Nijni
1290 pnt2d (i) = newCurr;
1292 else if (currOnLast) {
1293 // on decalle le point (i) en V First
1294 gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf RLN/Nijni
1295 pnt2d (i) = newCurr;
1299 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1300 if (dist2d > ( Vp / 2) ) {
1301 cout << "Echec dans le recadrage" << endl;
1308 //:c0 abv 20 Feb 98: treat very special case when 3d curve
1309 // go over the pole of, e.g., sphere, and partly lies along seam.
1310 // 2d representation of such a curve should consist of 3 parts - one on
1311 // regular part of surface (interior), one part along degenerated boundary
1312 // and one along seam.
1313 // Since it cannot be adjusted later by arranging pcurves (curve is single),
1314 // to fix it it is nesessary to have a possibility of adjusting seam
1315 // part of such curve either to left or right boundary of surface.
1316 // Test is performed only if flag AdjustOverDegen is not -1.
1317 // If AdjustOverDegen is True, seam part of curve is adjusted to
1318 // the left, and if False - to the right parametric boundary
1319 // If treated case is detected, flag DONE4 is set to status
1320 // NOTE: currently, precision is Precision::PConfusion() since it
1321 // is enough on encountered example
1322 // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1323 // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1324 if ( myAdjustOverDegen != -1 ) {
1325 if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1326 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1327 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1328 // 1st, find gap point (degenerated pole)
1329 Standard_Real PrevX=0.;
1330 Standard_Integer OnBound=0, PrevOnBound=0;
1331 Standard_Integer ind; // svv #1
1332 Standard_Boolean start = Standard_True;
1333 for ( ind=1; ind <= nbrPnt; ind++ ) {
1334 Standard_Real CurX = pnt2d(ind).X();
1335 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1336 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1338 OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1339 Precision::PConfusion() );
1340 if ( ! start && Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1342 start = Standard_False;
1344 PrevOnBound = OnBound;
1346 // if found, adjust seam part
1347 if ( ind <= nbrPnt ) {
1348 PrevX = ( myAdjustOverDegen ? uf : ul );
1349 Standard_Real dU = Up/2 + Precision::PConfusion();
1350 if ( PrevOnBound ) {
1351 pnt2d(ind-1).SetX ( PrevX );
1352 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1353 Standard_Real CurX = pnt2d(j).X();
1354 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1355 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1358 else if ( OnBound ) {
1359 pnt2d(ind).SetX ( PrevX );
1360 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1361 Standard_Real CurX = pnt2d(j).X();
1362 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1363 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1366 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1370 else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1371 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1372 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1373 // 1st, find gap point (degenerated pole)
1374 Standard_Real PrevY=0.;
1375 Standard_Integer OnBound=0, PrevOnBound=0;
1376 Standard_Integer ind; // svv #1
1377 Standard_Boolean start = Standard_True;
1378 for ( ind=1; ind <= nbrPnt; ind++ ) {
1379 Standard_Real CurY = pnt2d(ind).Y();
1380 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1381 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1383 OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1384 Precision::PConfusion() );
1385 if ( ! start && Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1387 start = Standard_False;
1389 PrevOnBound = OnBound;
1391 // if found, adjust seam part
1392 if ( ind <= nbrPnt ) {
1393 PrevY = ( myAdjustOverDegen ? vf : vl );
1394 Standard_Real dV = Vp/2 + Precision::PConfusion();
1395 if ( PrevOnBound ) {
1396 pnt2d(ind-1).SetY ( PrevY );
1397 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1398 Standard_Real CurY = pnt2d(j).Y();
1399 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1400 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1403 else if ( OnBound ) {
1404 pnt2d(ind).SetY ( PrevY );
1405 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1406 Standard_Real CurY = pnt2d(j).Y();
1407 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1408 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1411 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1419 if(ChangeCycle) { // msv 10.08.04: avoid using of uninitialised field
1420 //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1421 // myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1422 myCashe3d[0] = points(1);
1423 myCashe3d[1] = points(nbrPnt);
1424 myCashe2d[0] = pnt2d(1);
1425 myCashe2d[1] = pnt2d(nbrPnt);
1428 myCashe3d[1] = points(1);
1429 myCashe3d[0] = points(nbrPnt);
1430 myCashe2d[1] = pnt2d(1);
1431 myCashe2d[0] = pnt2d(nbrPnt);
1436 //=======================================================================
1437 //function : ApproximatePCurve
1439 //=======================================================================
1441 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1442 Handle(TColgp_HArray1OfPnt2d)& points2d,
1443 Handle(TColStd_HArray1OfReal)& params,
1444 const Handle(Geom_Curve)& /*orig*/) const
1446 // Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1447 Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1448 Handle(Geom2d_Curve) C2d;
1451 CheckPoints2d (points2d, params, theTolerance2d);
1452 Standard_Integer numberPnt = points2d->Length();
1454 TColgp_Array1OfPnt points3d(1,numberPnt);
1457 Standard_Integer i; // svv #1
1458 for( i = 1; i <= numberPnt; i++) {
1459 pnt2d = points2d->Value(i);
1460 pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1464 GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1465 Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1466 Standard_Integer NbPoles = crv3d->NbPoles();
1467 TColgp_Array1OfPnt poles3d (1, NbPoles);
1468 TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1469 crv3d->Poles(poles3d);
1470 for( i = 1; i <= NbPoles; i++) {
1471 pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1474 TColStd_Array1OfReal weights (1,NbPoles);
1475 TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1476 TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1477 crv3d->Knots(knots);
1478 crv3d->Weights(weights);
1479 crv3d->Multiplicities(multiplicities);
1480 C2d = new Geom2d_BSplineCurve ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1483 catch(Standard_Failure const& anException) {
1487 Standard_Integer nbp = params->Length();
1488 Standard_Integer nb2 = points2d->Length();
1489 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1490 anException.Print(cout);
1491 cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1492 // if (nb2 > nbp) nb2 = nbp;
1493 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1494 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1495 // dbl.AddReals (rb2,rbp);
1496 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1497 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1498 // dbl.AddXYZ (quoi);
1507 //=======================================================================
1508 //function : InterpolatePCurve
1510 //=======================================================================
1512 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1513 Handle(TColgp_HArray1OfPnt2d)& points2d,
1514 Handle(TColStd_HArray1OfReal)& params,
1515 const Handle(Geom_Curve)& /*orig*/) const
1517 Handle(Geom2d_Curve) C2d; // NULL si echec
1518 Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1521 // on verifie d abord s il n y a pas de points confondus
1522 // si besoin on retouche les valeurs ...
1523 CheckPoints2d (points2d, params, theTolerance2d);
1524 Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1525 Standard_False, theTolerance2d);
1526 myInterPol2d.Perform();
1527 if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1529 catch(Standard_Failure const& anException) {
1533 Standard_Integer nbp = params->Length();
1534 Standard_Integer nb2 = points2d->Length();
1535 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1536 anException.Print(cout);
1537 cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1538 // if (nb2 > nbp) nb2 = nbp;
1539 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1540 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1541 // dbl.AddReals (rb2,rbp);
1542 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1543 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1544 // dbl.AddXYZ (quoi);
1553 //=======================================================================
1554 //function : InterpolateCurve3d
1556 //=======================================================================
1558 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1559 Handle(TColgp_HArray1OfPnt)& points,
1560 Handle(TColStd_HArray1OfReal)& params,
1561 const Handle(Geom_Curve)& /*orig*/) const
1563 Handle(Geom_Curve) C3d; // NULL si echec
1566 Standard_Real Tol = myPreci;
1567 CheckPoints(points, params, Tol);
1568 GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1569 myInterPol.Perform();
1570 if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1572 catch(Standard_Failure const& anException) {
1575 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1576 anException.Print(cout); cout << endl;
1584 //=======================================================================
1585 //function : CheckPoints
1587 //=======================================================================
1589 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1591 Standard_Integer firstElem = points->Lower();
1592 Standard_Integer lastElem = points->Upper();
1594 Standard_Integer nbPntDropped = 0;
1595 Standard_Integer lastValid = firstElem; // indice of last undropped point
1597 // will store 0 when the point is to be removed, 1 otherwise
1598 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1599 for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
1600 Standard_Real DistMin2 = RealLast();
1601 gp_Pnt Prev = points->Value (lastValid);
1603 for (i = firstElem + 1; i <= lastElem ; i ++) {
1604 Curr = points->Value(i);
1605 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1606 if (CurDist2 < gp::Resolution()) { // test 0
1608 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1609 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1611 if (CurDist2 < DistMin2)
1612 DistMin2 = CurDist2;
1613 // lastValid becomes the current (i.e. i)
1618 if (DistMin2 < RealLast())
1619 preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1620 if (nbPntDropped == 0)
1624 cout << "Warning : removing 3d points for interpolation" << endl;
1626 // Build new HArrays
1627 Standard_Integer newLast = lastElem - nbPntDropped;
1628 if ((newLast - firstElem + 1) < 2) {
1630 cout << "Too many degenerated points for 3D interpolation" << endl;
1634 Handle(TColgp_HArray1OfPnt) newPnts =
1635 new TColgp_HArray1OfPnt(firstElem, newLast);
1636 Handle(TColStd_HArray1OfReal) newParams =
1637 new TColStd_HArray1OfReal(firstElem, newLast);
1638 Standard_Integer newCurr = 1;
1639 for (i = firstElem; i<= lastElem ; i++) {
1640 if (tmpParam.Value(i) == 1) {
1641 newPnts->SetValue(newCurr, points->Value(i));
1642 newParams->SetValue(newCurr, params->Value(i));
1648 // on la reduit un peu
1651 //=======================================================================
1652 //function : CheckPoints2d
1654 //=======================================================================
1656 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1657 Handle(TColStd_HArray1OfReal)& params,
1658 Standard_Real& preci) const
1660 Standard_Integer firstElem = points->Lower();
1661 Standard_Integer lastElem = points->Upper();
1663 Standard_Integer nbPntDropped = 0;
1664 Standard_Integer lastValid = firstElem; // indice of last undropped point
1666 // will store 0 when the point is to be removed, 1 otherwise
1667 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1668 for (i = firstElem; i<=lastElem ; i++) {
1669 tmpParam.SetValue(i,1);
1671 Standard_Real DistMin2 = RealLast();
1672 gp_Pnt2d Prev = points->Value(lastValid);
1674 for (i = firstElem + 1; i<=lastElem ; i++) {
1675 Curr = points->Value(i);
1676 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1677 if (CurDist2 < gp::Resolution()) { // test 0
1679 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1680 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1682 if (CurDist2 < DistMin2)
1683 DistMin2 = CurDist2;
1684 // lastValid becomes the current (i.e. i)
1689 if (DistMin2 < RealLast())
1690 preci = 0.9 * Sqrt (DistMin2);
1691 if (nbPntDropped == 0)
1695 cout << "Warning : removing 2d points for interpolation" << endl;
1697 // Build new HArrays
1698 Standard_Integer newLast = lastElem - nbPntDropped;
1699 if ((newLast - firstElem + 1) < 2) {
1701 cout << "Too many degenerated points for 2D interpolation" << endl;
1703 //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1704 tmpParam.SetValue(firstElem,1);
1705 tmpParam.SetValue(lastElem,1);
1706 gp_XY lastPnt = points->Value(lastElem).XY();
1707 lastPnt.Add(gp_XY(preci,preci));
1708 points->SetValue(lastElem,lastPnt);
1709 newLast = firstElem+1;
1712 Handle(TColgp_HArray1OfPnt2d) newPnts =
1713 new TColgp_HArray1OfPnt2d(firstElem, newLast);
1714 Handle(TColStd_HArray1OfReal) newParams =
1715 new TColStd_HArray1OfReal(firstElem, newLast);
1716 Standard_Integer newCurr = 1;
1717 for (i = firstElem; i <= lastElem ; i++) {
1718 if (tmpParam.Value(i) == 1) {
1720 cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1722 newPnts->SetValue(newCurr, points->Value(i));
1723 newParams->SetValue(newCurr, params->Value(i));
1728 cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1736 //=======================================================================
1737 //function : IsAnIsoparametric
1739 //=======================================================================
1740 //:S4030: modified for optimization
1741 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1743 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1744 const TColgp_Array1OfPnt& points,
1745 const TColStd_Array1OfReal& params,
1746 Standard_Boolean& isoTypeU,
1747 Standard_Boolean& p1OnIso,
1749 Standard_Boolean& p2OnIso,
1751 Standard_Boolean& isoPar2d3d,
1752 Handle(Geom_Curve)& cIso,
1755 TColStd_Array1OfReal& pout) const
1760 Standard_Real prec = Precision::Confusion();//myPreci;
1762 Standard_Boolean isoParam = Standard_False;
1763 isoPar2d3d = Standard_False;
1765 Standard_Real U1, U2, V1, V2;
1766 mySurf->Bounds(U1, U2, V1, V2);
1768 if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1769 Handle(Geom_RectangularTrimmedSurface) sTrim =
1770 Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1771 sTrim->Bounds(U1, U2, V1, V2);
1775 Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
1776 Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
1777 Standard_Real mindist2;
1778 Standard_Real mind2[2];
1779 mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1781 p1OnIso = Standard_False;
1782 p2OnIso = Standard_False;
1783 const Bnd_Box* aBox = 0;
1785 for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1786 Standard_Real isoVal=0.;
1787 Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1788 Handle(Geom_Curve) cI;
1789 Standard_Real tt1, tt2;
1792 if (Precision::IsInfinite(U1)) continue;
1793 cI = mySurf->UIso(U1);
1794 isoU = Standard_True;
1796 aBox = & mySurf->GetBoxUF();
1799 if (Precision::IsInfinite(U2)) continue;
1800 cI = mySurf->UIso(U2);
1801 isoU = Standard_True;
1803 aBox = & mySurf->GetBoxUL();
1806 if (Precision::IsInfinite(V1)) continue;
1807 cI = mySurf->VIso(V1);
1808 isoU = Standard_False;
1810 aBox = & mySurf->GetBoxVF();
1813 if (Precision::IsInfinite(V2)) continue;
1814 cI = mySurf->VIso(V2);
1815 isoU = Standard_False;
1817 aBox = & mySurf->GetBoxVL();
1822 if (isoU) { tt1 = V1; tt2 = V2; }
1823 else { tt1 = U1; tt2 = U2; }
1829 // PATCH CKY 9-JUL-1998 : protection contre singularite
1831 cI->D0( (tt1+tt2)/2,extmi);
1832 if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1834 Standard_Boolean PtEQext1 = Standard_False;
1835 Standard_Boolean PtEQext2 = Standard_False;
1837 Standard_Real currd2[2], tp[2] = {0, 0};
1838 Standard_Integer mp[2];
1840 for (Standard_Integer i=0; i<2; i++) {
1842 Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1844 // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1845 currd2[i] = points(k).SquareDistance ( ext1 );
1846 if ( currd2[i] <= prec*prec && !PtEQext1) {
1849 PtEQext1 = Standard_True;
1853 currd2[i] = points(k).SquareDistance ( ext2 );
1854 if ( currd2[i] <= prec*prec && !PtEQext2) {
1857 PtEQext2 = Standard_True;
1861 // On evite de projecter sur un iso degenere
1862 // on doit egalement le faire pour l apex du cone
1863 if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1867 if(aBox->IsOut(points(k))) continue;
1869 Standard_Real Cf = cI->FirstParameter();
1870 Standard_Real Cl = cI->LastParameter();
1871 if (Precision::IsInfinite(Cf)) Cf = -1000;
1872 if (Precision::IsInfinite(Cl)) Cl = +1000;
1874 ShapeAnalysis_Curve sac;
1875 Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1876 currd2[i] = dist * dist;
1877 if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1883 //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1884 // avoid possible null-length curves
1885 if ( mp[0] >0 && mp[1] >0 &&
1886 Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1890 ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1891 p1OnIso = Standard_True;
1892 mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
1893 if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1894 else valueP1.SetCoord(tp[0], isoVal);
1898 ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1899 p2OnIso = Standard_True;
1900 mind2[1] = currd2[1];
1901 if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1902 else valueP2.SetCoord(tp[1], isoVal);
1905 if ( mp[0] <=0 || mp[1] <=0 ) continue;
1907 Standard_Real md2 = currd2[0] + currd2[1];
1908 if ( mindist2 <= md2 ) continue;
1922 // probablely it concerns an isoparametrics
1923 if ( mpt[0] >0 && mpt[1] >0 ) {
1925 p1OnIso = p2OnIso = Standard_True;
1927 valueP1.SetCoord(isoValue, tpar[0]);
1928 valueP2.SetCoord(isoValue, tpar[1]);
1931 valueP1.SetCoord(tpar[0], isoValue);
1932 valueP2.SetCoord(tpar[1], isoValue);
1935 if ( mpt[0] != 3 && mpt[1] != 3 ) {
1936 isoPar2d3d = Standard_True;
1937 for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1938 if (tpar[1] > tpar[0]) t = params(i);
1939 else t = t1+t2-params(i);
1941 if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1945 if (isoPar2d3d) isoParam = Standard_True;
1947 Standard_Real prevParam = tpar[0];
1948 Standard_Real Cf, Cl;
1949 Standard_Boolean isoByDistance = Standard_True;
1950 Cf = cIso->FirstParameter();
1951 Cl = cIso->LastParameter();
1952 if (Precision::IsInfinite(Cf)) Cf = -1000;
1953 if (Precision::IsInfinite(Cl)) Cl = +1000;
1955 ShapeAnalysis_Curve sac;
1956 for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1957 Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1959 Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1962 if( (dist > prec) || (t < Cf) || (t > Cl) )
1963 isoByDistance = Standard_False;
1965 if (isoByDistance) isoParam = Standard_True;
1968 /* if (!isoParam) { CKY 29-mai-1997 : garder tout ce qu on peut ?
1969 p1OnIso = Standard_False;
1970 p2OnIso = Standard_False;
1974 catch(Standard_Failure const& anException) {
1976 // pb : on affiche ce qu on peut
1977 for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1978 cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1979 points(numpnt).X()<<" "<<points(numpnt).Y()<<" "<<points(numpnt).Z()<<")"<<endl;
1981 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
1982 anException.Print(cout); cout << endl;
1985 return Standard_False;
1989 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
1990 //=======================================================================
1991 //function : BestExtremum
1992 //purpose : auxiliaire prenant le meilleur extremum si ISO car doute possible
1993 //=======================================================================
1995 gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
1997 // P2iso a ete calcule depuis P3ext sur une iso externe de la surface
1998 // En principe bon mais circularite possible ... et IsU/VClosed faillible
1999 // (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
2000 // 1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
2001 // 2/ si iso, attention a la circularite (cas limite)
2003 // NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
2005 // D abord, calcul p2ext depuis la surface. choix surface/iso
2007 Standard_Real prec = Precision::Confusion();//myPreci;
2008 gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
2009 gp_Pnt P3cal = mySurf->Value (P2cal);
2010 Standard_Real dcal = P3ext.Distance (P3cal);
2011 Standard_Real dnxt = P3ext.Distance (P3next);
2012 if (dcal > dnxt) return P2iso; // en fait protection sur BUG (PRO8468)
2014 // On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
2015 gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
2016 if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;