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);
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 aLength3d += p3d2.Distance(p3d1);
326 const Standard_Real aCoeff = aLength3d / (aLastParam - aFirstParam);
327 if (Abs(aCoeff) > gp::Resolution())
328 aKnotCoeffs.Append(aCoeff);
329 aFirstParam = aLastParam;
332 Standard_Real anEvenlyCoeff = 0;
333 if (aKnotCoeffs.Size() > 0)
335 anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
336 *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
339 const Standard_Real aMaxQuotientCoeff = 1500.0;
340 if (anEvenlyCoeff > aMaxQuotientCoeff)
342 PerformByProjLib(c3d, First, Last, c2d);
343 // PerformByProjLib fail detection:
346 return Status (ShapeExtend_DONE);
351 Standard_Real deltaT, t;
352 deltaT = (Last - First) / (nbPini-1);
354 for (iPnt = 1; iPnt <= nbPini; iPnt ++)
356 if (iPnt == 1) t = First;
357 else if (iPnt == nbPini) t = Last;
358 else t = First + (iPnt - 1) * deltaT;
365 // CALCUL par approximation
366 TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
367 ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
369 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
370 return Standard_True;
371 }// cas particulier d iso
373 // INTERPOLATION du resultat
376 Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt (1, nbPini);
377 Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
378 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
379 thePnts->SetValue(iPnt, points(iPnt));
380 theParams->SetValue(iPnt, params(iPnt));
383 Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
384 if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
386 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
391 Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
392 Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
393 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
394 theParams2d->SetValue(iPnt, params(iPnt));
395 thePnts2d->SetValue(iPnt, pnt2d(iPnt));
398 c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
399 // c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
400 // Faut-il aussi reprendre la C3D ?
401 myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
402 return Status (ShapeExtend_DONE);
405 //=======================================================================
406 //function : PerformByProjLib
408 //=======================================================================
410 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
411 const Standard_Real First,
412 const Standard_Real Last,
413 Handle(Geom2d_Curve)& c2d,
414 const GeomAbs_Shape /*continuity*/,
415 const Standard_Integer /*maxdeg */,
416 const Standard_Integer /*nbinterval */)
418 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
420 if (mySurf.IsNull()) {
421 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
422 return Standard_False;
428 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
429 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
430 ProjLib_ProjectedCurve Projector(GAS, GAC);
432 switch (Projector.GetType())
435 c2d = new Geom2d_Line(Projector.Line());
437 case GeomAbs_Circle :
438 c2d = new Geom2d_Circle(Projector.Circle());
440 case GeomAbs_Ellipse :
441 c2d = new Geom2d_Ellipse(Projector.Ellipse());
443 case GeomAbs_Parabola :
444 c2d = new Geom2d_Parabola(Projector.Parabola());
446 case GeomAbs_Hyperbola :
447 c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
449 case GeomAbs_BSplineCurve :
450 c2d = Projector.BSpline();
453 // Not possible, handling added to avoid gcc warning.
459 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
460 return Standard_False;
464 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
465 return Standard_True;
469 catch(Standard_Failure)
472 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
473 Standard_Failure::Caught()->Print(cout); cout << endl;
475 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
478 return Standard_False;
481 //=======================================================================
482 //function : PerformAdvanced
484 //=======================================================================
486 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
487 const Standard_Real First,
488 const Standard_Real Last,
489 Handle(Geom2d_Curve)& c2d)
491 Standard_Boolean hasResult = Standard_False;
492 Standard_Integer nbintervals;
494 Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
495 // && (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
497 if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
499 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
500 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
501 nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
502 isStandard = (nbintervals < 2);
505 hasResult = PerformByProjLib(c3d, First, Last, c2d);
507 if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
511 //=======================================================================
512 //function : ProjectAnalytic
514 //=======================================================================
516 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
518 Handle(Geom2d_Curve) result;
520 //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
521 // This is necessary for K4L since it fails on other surfaces
522 // when general method GeomProjLib::Curve2d() is used
523 // Projection is done as in BRep_Tool and BRepCheck_Edge
524 Handle(Geom_Surface) surf = mySurf->Surface();
525 Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
526 if ( Plane.IsNull() ) {
527 Handle(Geom_RectangularTrimmedSurface) RTS =
528 Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
529 if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
531 Handle(Geom_OffsetSurface) OS =
532 Handle(Geom_OffsetSurface)::DownCast ( surf );
534 Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
537 if ( ! Plane.IsNull() ) {
538 Handle(Geom_Curve) ProjOnPlane =
539 GeomProjLib::ProjectOnPlane (c3d, Plane,
540 Plane->Position().Direction(), Standard_True);
541 Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
542 ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
544 result = Geom2dAdaptor::MakeCurve(Proj);
545 if ( result.IsNull() ) return result;
546 if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
547 Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
548 result = TC->BasisCurve();
557 //! Fix possible period jump and handle walking period parameter.
558 static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
559 Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
560 Standard_Real thePeriod, // Period on objective coord
561 Standard_Integer theSavedPoint, // Point number to choose period
562 Standard_Real theSavedParam) // Param from cashe to choose period
564 Standard_Real aSavedParam;
565 Standard_Integer aSavedPoint;
566 Standard_Real aMinParam = 0.0, aMaxParam = thePeriod;
567 if (theSavedPoint < 0) {
568 // normalize to first period by default
569 aSavedParam = 0.5 * thePeriod;
573 aSavedParam = theSavedParam;
574 aSavedPoint = theSavedPoint;
575 while (aMinParam > aSavedParam) {
576 aMinParam -= thePeriod;
577 aMaxParam -= thePeriod;
579 while (aMaxParam < aSavedParam) {
580 aMinParam += thePeriod;
581 aMaxParam += thePeriod;
585 Standard_Real aFixIsoParam = aMinParam;
586 Standard_Boolean isIsoLine = Standard_False;
587 if (aMaxParam - aSavedParam < Precision::PConfusion() ||
588 aSavedParam - aMinParam < Precision::PConfusion()) {
589 aFixIsoParam = aSavedParam;
590 isIsoLine = Standard_True;
592 // normalize all coordinates to [aMinParam, aMaxParam)
593 for (Standard_Integer i = 0; i < 4; i++) {
594 Standard_Real aParam = thePnt[i].Coord(theIdx);
595 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(aParam, aMinParam, aMaxParam);
597 // Walk over period coord -> not walking on another isoline in parameter space.
599 if (aMaxParam - aParam < Precision::PConfusion() || aParam - aMinParam < Precision::PConfusion())
600 aParam = aFixIsoParam;
603 if (aMaxParam - aParam < Precision::PConfusion())
605 if (aParam - aMinParam < Precision::PConfusion())
609 thePnt[i].SetCoord(theIdx, aParam);
612 // find possible period jump and increasing or decreasing coordinates vector
613 Standard_Boolean isJump = Standard_False;
614 Standard_Real aPrevDiff = 0.0;
615 Standard_Real aSumDiff = 1.0;
616 for (Standard_Integer i = 0; i < 3; i++) {
617 Standard_Real aDiff = thePnt[i + 1].Coord(theIdx) - thePnt[i].Coord(theIdx);
618 if (aDiff < -Precision::PConfusion()) {
621 //if first derivative changes its sign then period jump may exists in this place
622 if (aDiff * aPrevDiff < -Precision::PConfusion()) {
623 isJump = Standard_True;
629 return Standard_False;
631 if (aSumDiff > 0) { // decreasing sequence (parameters decrease twice(--) and one period jump(+))
632 for (Standard_Integer i = aSavedPoint; i > 0; i--)
633 if (thePnt[i].Coord(theIdx) > thePnt[i - 1].Coord(theIdx)) {
634 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) + thePeriod);
636 for (Standard_Integer i = aSavedPoint; i < 3; i++)
637 if (thePnt[i].Coord(theIdx) < thePnt[i + 1].Coord(theIdx)) {
638 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) - thePeriod);
641 else {// increasing sequence (parameters increase twice(++) and one period jump(-))
642 for (Standard_Integer i = aSavedPoint; i > 0; i--)
643 if (thePnt[i].Coord(theIdx) < thePnt[i - 1].Coord(theIdx)) {
644 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) - thePeriod);
646 for (Standard_Integer i = aSavedPoint; i < 3; i++)
647 if (thePnt[i].Coord(theIdx) > thePnt[i + 1].Coord(theIdx)) {
648 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) + thePeriod);
652 // Do not return false, because for nonlinear 2d curves vector of parameters
653 // may change its first derivative and shifted parameters will be broken for this case.
654 return Standard_True;
657 //=======================================================================
660 //=======================================================================
662 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
663 const TColgp_Array1OfPnt& thepoints,
664 const TColStd_Array1OfReal& theparams,
665 TColgp_Array1OfPnt2d& thePnt2ds,
666 Standard_Real theTol,
667 Standard_Boolean &isRecompute,
668 Standard_Boolean &isFromCashe) const
670 Standard_Integer nb = thepoints.Length();
672 aP[0] = thepoints(1);
673 aP[1] = thepoints(2);
674 aP[2] = thepoints(nb - 1);
675 aP[3] = thepoints(nb);
677 Standard_Integer i = 0;
679 Standard_Real aTol2 = theTol * theTol;
680 Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
681 Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
684 // Protection against bad "tolerance" shapes.
687 theTol = Precision::Confusion();
688 aTol2 = theTol * theTol;
690 if (aTol2 < Precision::SquareConfusion())
691 aTol2 = Precision::SquareConfusion();
692 Standard_Real anOldTol2 = aTol2;
693 // auxiliary variables to choose period for connection with previous 2dcurve (if exist)
694 Standard_Integer aSavedPointNum = -1;
695 gp_Pnt2d aSavedPoint;
697 // project first and last points
701 for ( j=0; j < myNbCashe; j++ )
702 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
704 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol,
707 aSavedPoint = myCashe2d[j];
709 isFromCashe = Standard_True;
712 if ( j >= myNbCashe )
713 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
715 Standard_Real aDist = mySurf->Gap();
716 Standard_Real aCurDist = aDist * aDist;
717 if( aTol2 < aDist * aDist)
721 if ( isPeriodicU || isPeriodicV )
723 // Compute second and last but one c2d points.
724 for(i = 1; i < 3; i++)
727 for ( j=0; j < myNbCashe; j++ )
728 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
730 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
732 aSavedPoint = myCashe2d[j];
735 if ( j >= myNbCashe )
736 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
738 Standard_Real aDist = mySurf->Gap();
739 Standard_Real aCurDist = aDist * aDist;
740 if( aTol2 < aDist * aDist)
746 isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
751 isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
755 if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
756 // Do not try to make line, because in this case may be very special case when 3d curve
757 // go over the pole of, e.g., sphere, and partly lies along seam
758 // (see ApproxPCurve() for more information).
762 thePnt2ds.SetValue(1, aP2d[0]);
763 thePnt2ds.SetValue(nb, aP2d[3]);
765 // Restore old tolerance in 2d space to avoid big gap cases.
767 // Check that straight line in 2d with parameterisation as in 3d will fit
768 // fit 3d curve at all points.
769 Standard_Real dPar = theparams(nb) - theparams(1);
770 if ( Abs(dPar) < Precision::PConfusion() )
772 gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
773 gp_Vec2d aVec = aVec0 / dPar;
774 Handle(Geom_Surface) aSurf = mySurf->Surface();
775 Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
777 for(i = 1; i <= nb; i++)
779 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
781 gp_Vec aNormalVec, aDu, aDv;
782 aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
783 aNormalVec = aDu ^ aDv;
784 if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
785 isNormalCheck = Standard_False;
788 gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
789 Standard_Real aDist = aNormalLine.Distance(thepoints(i));
794 if (!isNormalCheck) {
795 Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
796 SquareDistance(thepoints(1));
797 aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
798 for(i = 2; i < nb; i++)
800 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
802 aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
803 Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
805 if(Abs (aFirstPointDist - aDist1) > aTol2)
810 // check if pcurve can be represented by Geom2d_Line (parameterised by length)
811 Standard_Real aLLength = aVec0.Magnitude();
812 if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
814 gp_XY aDirL = aVec0.XY() / aLLength;
815 gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
816 return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
819 // create straight bspline
820 TColgp_Array1OfPnt2d aPoles(1, 2);
824 TColStd_Array1OfReal aKnots(1,2);
825 aKnots(1) = theparams(1);
826 aKnots(2) = theparams(theparams.Length());
828 TColStd_Array1OfInteger aMults(1,2);
831 Standard_Integer aDegree = 1;
832 Handle(Geom2d_BSplineCurve) abspl2d =
833 new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
837 //=======================================================================
838 //function : ApproxPCurve
840 //=======================================================================
842 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
843 const TColgp_Array1OfPnt& points,
844 const TColStd_Array1OfReal& params,
845 TColgp_Array1OfPnt2d& pnt2d,
846 Handle(Geom2d_Curve)& c2d)
848 // for performance, first try to handle typical case when pcurve is straight
849 Standard_Boolean isRecompute = Standard_False;
850 Standard_Boolean isFromCasheLine = Standard_False;
851 c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
855 Standard_Boolean ChangeCycle = Standard_False;
856 if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
857 myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
858 ChangeCycle = Standard_True;
861 myCashe3d[0] = points(1);
862 myCashe3d[1] = points(nbrPnt);
863 myCashe2d[0] = pnt2d(1);
864 myCashe2d[1] = pnt2d(nbrPnt);
867 myCashe3d[1] = points(1);
868 myCashe3d[0] = points(nbrPnt);
869 myCashe2d[1] = pnt2d(1);
870 myCashe2d[0] = pnt2d(nbrPnt);
872 return Standard_True;
874 Standard_Boolean isDone = Standard_True;
875 // test if the curve 3d is a boundary of the surface
876 // (only for Bezier or BSpline surface)
878 Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
879 gp_Pnt2d valueP1, valueP2;
880 Handle(Geom_Curve) cIso;
881 Standard_Real t1, t2;
883 Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
884 Standard_Boolean isAnalytic = Standard_True;
885 if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
886 Standard_Real uf, ul, vf, vl;
887 mySurf->Surface()->Bounds(uf, ul, vf, vl);
888 isoclosed = Standard_False;
889 TColStd_Array1OfReal pout(1, nbrPnt);
891 isoParam = IsAnIsoparametric(nbrPnt, points, params,
892 isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
893 isoPar2d3d, cIso, t1, t2, pout);
895 // projection of the points on surfaces
900 Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
901 Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
903 // Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
904 // Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
907 isoValue = valueP1.X();
908 isoPar1 = valueP1.Y();
909 isoPar2 = valueP2.Y();
910 isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
911 parf = vf; parl = vl;
914 isoValue = valueP1.Y();
915 isoPar1 = valueP1.X();
916 isoPar2 = valueP2.X();
917 isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
918 parf = uf; parl = ul;
920 if (!isoPar2d3d && !isAnalytic) {
921 Cf = cIso->FirstParameter();
922 Cl = cIso->LastParameter();
923 if (Precision::IsInfinite(Cf)) Cf = -1000;
924 if (Precision::IsInfinite(Cl)) Cl = +1000;
925 //pdn S4030 optimizing and fix isopar case on PRO41323
927 // dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
928 // Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
929 if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
930 if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
932 if (isoTypeU) valueP1.SetY (isoPar1);
933 else valueP1.SetX (isoPar1);
935 if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
936 //pdn S4030 optimizing and fix isopar case on PRO41323
937 tfin = pout(nbrPnt-1);
938 //dist = ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
939 if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
941 if (isoTypeU) valueP2.SetY (isoPar2);
942 else valueP2.SetX (isoPar2);
945 // Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
946 // Est-ce encore necessaire apres ce qui vient d etre fait ?
948 // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
949 // if curve is periodic do not change the points
950 // skl change "if" for pout(nbrPnt-1) 19.11.2003
952 if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
953 (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
954 gp_Pnt2d valueTmp = valueP1;
955 valueP1 = valueP2; valueP2 = valueTmp;
957 isoValue = valueP1.X();
958 isoPar1 = valueP1.Y();
959 isoPar2 = valueP2.Y();
962 isoValue = valueP1.Y();
963 isoPar1 = valueP1.X();
964 isoPar2 = valueP2.X();
966 // Fin calcul sens de courbe iso
968 } // end of fix check 05.02.02
972 // Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
973 // (utile pour detromper bug de projection). Mais detromper aussi circularite
975 //if (p1OnIso) valueP1 =
976 //BestExtremum (valueP1,points(1),points(2));
977 //if (p2OnIso) valueP2 =
978 //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
981 Standard_Real gap = myPreci; //:q1
982 Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
983 // auxiliaruy variables to shift 2dcurve, according to previous
984 Standard_Boolean isFromCashe = Standard_False;
985 gp_Pnt2d aSavedPoint;
986 if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
987 //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
988 if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
989 ChangeCycle = Standard_True;
990 //for( i = 1; i <= nbrPnt; i ++) {
991 for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
992 if(ChangeCycle) //skl for OCC3430
1000 if (isoPar2 > isoPar1) tPar = params(i);
1001 else tPar = t1 + t2 - params(i);
1002 } else if (!isAnalytic) {
1003 // projection to iso
1004 if (i==1) tPar = isoPar1;
1005 else if (i==nbrPnt) tPar = isoPar2;
1008 //:S4030 ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
1012 if (!isoPar2d3d && isAnalytic) {
1013 if (i == 1) p2d = valueP1;
1014 else if (i == nbrPnt) p2d = valueP2;
1016 p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
1017 Precision::Confusion()+1000*gap); //:q1
1018 gap = mySurf->Gap();
1021 if(isoTypeU) { p2d.SetX(isoValue); p2d.SetY(tPar); }
1022 else { p2d.SetX(tPar); p2d.SetY(isoValue); }
1027 if ( (i == 1) && p1OnIso) p2d = valueP1;
1028 else if( (i == nbrPnt) && p2OnIso) p2d = valueP2;
1029 else {// general case (not an iso) mais attention aux singularites !
1030 // first and last points are already computed by getLine()
1031 if ( (i == 1 || i == nbrPnt))
1036 gap = mySurf->Gap();
1038 isFromCashe = isFromCasheLine;
1045 //:q9 abv 23 Mar 99: use cashe as 1st approach
1046 Standard_Integer j; // svv #1
1047 for ( j=0; j < myNbCashe; j++ )
1048 if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
1050 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci,
1051 Precision::Confusion()+gap);
1053 isFromCashe = Standard_True;
1054 aSavedPoint = myCashe2d[j];
1058 if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci);
1062 p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
1063 Precision::Confusion()+1000*gap); //:q1
1065 gap = mySurf->Gap();
1071 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
1073 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
1077 //pdn %12 11.02.99 PRO9234 entity 15402
1079 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
1080 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
1083 // attention aux singularites ... (hors cas iso qui les traite deja)
1086 // if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
1088 // p2d = pnt2d (nbrPnt);
1089 // if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
1090 // pnt2d (nbrPnt) = p2d;
1093 // Si la surface est UCLosed et VClosed, on recadre les points
1094 // algo un peu complique, on retarde l implementation
1095 Standard_Real Up = ul - uf;
1096 Standard_Real Vp = vl - vf;
1097 Standard_Real dist2d;
1099 if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
1100 cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
1103 // Si la surface est UCLosed, on recadre les points
1104 if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
1105 // Premier point dans le domain [uf, ul]
1106 Standard_Real prevX, firstX = pnt2d (1).X();
1108 // do not shift 2dcurve, if it connects to previous
1109 while (firstX < uf) { firstX += Up; pnt2d (1).SetX(firstX); }
1110 while (firstX > ul) { firstX -= Up; pnt2d (1).SetX(firstX); }
1112 // shift first point, according to cashe
1113 if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
1114 Standard_Real aMinParam = uf, aMaxParam = ul;
1115 while (aMinParam > aSavedPoint.X()) {
1119 while (aMaxParam < aSavedPoint.X()) {
1123 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
1125 pnt2d(1).SetX(firstX);
1129 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1130 Standard_Real minX = firstX, maxX = firstX;
1132 // On decalle toujours le suivant
1133 for (i = 2; i <= nbrPnt; i++) {
1134 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
1135 Standard_Real CurX = pnt2d (i).X();
1136 dist2d = Abs (CurX - prevX);
1137 if (dist2d > ( Up / 2) ) {
1138 if (CurX > prevX + Up/2) {
1139 while (CurX > prevX + Up/2) { CurX -= Up; pnt2d (i).SetX (CurX); }
1140 } else if (CurX < prevX - Up/2) {
1141 while (CurX < prevX - Up/2) { CurX += Up; pnt2d (i).SetX (CurX); }
1146 if ( minX > CurX ) minX = CurX; //:97
1147 else if ( maxX < CurX ) maxX = CurX; //:97
1152 // do not shift 2dcurve, if it connects to previous
1153 Standard_Real midX = 0.5 * ( minX + maxX );
1154 Standard_Real shiftX=0.;
1155 if ( midX > ul ) shiftX = -Up;
1156 else if ( midX < uf ) shiftX = Up;
1158 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
1161 // Si la surface est VCLosed, on recadre les points
1162 // Same code as UClosed : optimisation souhaitable !!
1163 // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
1164 // Si la surface est UCLosed, on recadre les points
1166 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1167 //#78 rln 12.03.99 S4135
1168 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1169 // Premier point dans le domain [vf, vl]
1170 Standard_Real prevY, firstY = pnt2d (1).Y();
1172 // do not shift 2dcurve, if it connects to previous
1173 while (firstY < vf) { firstY += Vp; pnt2d (1).SetY(firstY); }
1174 while (firstY > vl) { firstY -= Vp; pnt2d (1).SetY(firstY); }
1176 // shift first point, according to cashe
1177 if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
1178 Standard_Real aMinParam = vf, aMaxParam = vl;
1179 while (aMinParam > aSavedPoint.Y()) {
1183 while (aMaxParam < aSavedPoint.Y()) {
1187 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
1189 pnt2d(1).SetY(firstY);
1193 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1194 Standard_Real minY = firstY, maxY = firstY;
1196 // On decalle toujours le suivant
1197 for (i = 2; i <= nbrPnt; i ++) {
1198 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
1199 Standard_Real CurY = pnt2d (i).Y();
1200 dist2d = Abs (CurY - prevY);
1201 if (dist2d > ( Vp / 2) ) {
1202 if (CurY > prevY + Vp/2) {
1203 while (CurY > prevY + Vp/2) { CurY -= Vp; pnt2d (i).SetY (CurY); }
1204 } else if (CurY < prevY - Vp/2) {
1205 while (CurY < prevY - Vp/2) { CurY += Vp; pnt2d (i).SetY (CurY); }
1209 if ( minY > CurY ) minY = CurY; //:97
1210 else if ( maxY < CurY ) maxY = CurY; //:97
1215 // do not shift 2dcurve, if it connects to previous
1216 Standard_Real midY = 0.5 * ( minY + maxY );
1217 Standard_Real shiftY=0.;
1218 if ( midY > vl ) shiftY = -Vp;
1219 else if ( midY < vf ) shiftY = Vp;
1221 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
1225 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1226 //#78 rln 12.03.99 S4135
1227 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1228 for (i = 2; i <= nbrPnt; i++) {
1229 //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1230 dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
1231 if (dist2d > ( Vp / 2) ) {
1232 // ATTENTION : il faut regarder ou le decalage se fait.
1233 // si plusieurs points sont decalles, il faut plusieurs passes
1234 // pour obtenir un resultat correct.
1235 // NOT YET IMPLEMENTED
1237 // one of those point is incorrectly placed
1238 // i.e on the wrong side of the "seam"
1239 // on prend le point le plus pres des bords vf ou vl
1240 Standard_Boolean prevOnFirst = Standard_False;
1241 Standard_Boolean prevOnLast = Standard_False;
1242 Standard_Boolean currOnFirst = Standard_False;
1243 Standard_Boolean currOnLast = Standard_False;
1245 // .X ? plutot .Y , non ?
1246 Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
1247 Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
1248 Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
1249 Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
1251 Standard_Real theMin = distPrevVF;
1252 prevOnFirst = Standard_True;
1253 if (distPrevVL < theMin) {
1254 theMin = distPrevVL;
1255 prevOnFirst = Standard_False;
1256 prevOnLast = Standard_True;
1258 if (distCurrVF < theMin) {
1259 theMin = distCurrVF;
1260 prevOnFirst = Standard_False;
1261 prevOnLast = Standard_False;
1262 currOnFirst = Standard_True;
1264 if (distCurrVL < theMin) {
1265 theMin = distCurrVL;
1266 prevOnFirst = Standard_False;
1267 prevOnLast = Standard_False;
1268 currOnFirst = Standard_False;
1269 currOnLast = Standard_True;
1271 // Modifs RLN/Nijni 3-DEC-1997
1273 // on decalle le point (i-1) en V Last
1274 gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of vl RLN/Nijni
1275 pnt2d (i-1) = newPrev;
1277 else if (prevOnLast) {
1278 // on decalle le point (i-1) en V first
1279 gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of vf RLN/Nijni
1280 pnt2d (i-1) = newPrev;
1282 else if (currOnFirst) {
1283 // on decalle le point (i) en V Last
1284 gp_Pnt2d newCurr(pnt2d (i).X(),vf); // instead of vl RLN/Nijni
1285 pnt2d (i) = newCurr;
1287 else if (currOnLast) {
1288 // on decalle le point (i) en V First
1289 gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf RLN/Nijni
1290 pnt2d (i) = newCurr;
1294 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1295 if (dist2d > ( Vp / 2) ) {
1296 cout << "Echec dans le recadrage" << endl;
1303 //:c0 abv 20 Feb 98: treat very special case when 3d curve
1304 // go over the pole of, e.g., sphere, and partly lies along seam.
1305 // 2d representation of such a curve should consist of 3 parts - one on
1306 // regular part of surface (interior), one part along degenerated boundary
1307 // and one along seam.
1308 // Since it cannot be adjusted later by arranging pcurves (curve is single),
1309 // to fix it it is nesessary to have a possibility of adjusting seam
1310 // part of such curve either to left or right boundary of surface.
1311 // Test is performed only if flag AdjustOverDegen is not -1.
1312 // If AdjustOverDegen is True, seam part of curve is adjusted to
1313 // the left, and if False - to the right parametric boundary
1314 // If treated case is detected, flag DONE4 is set to status
1315 // NOTE: currently, precision is Precision::PConfusion() since it
1316 // is enough on encountered example
1317 // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1318 // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1319 if ( myAdjustOverDegen != -1 ) {
1320 if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1321 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1322 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1323 // 1st, find gap point (degenerated pole)
1324 Standard_Real PrevX=0.;
1325 Standard_Integer OnBound=0, PrevOnBound=0;
1326 Standard_Integer ind; // svv #1
1327 Standard_Boolean start = Standard_True;
1328 for ( ind=1; ind <= nbrPnt; ind++ ) {
1329 Standard_Real CurX = pnt2d(ind).X();
1330 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1331 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1333 OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1334 Precision::PConfusion() );
1335 if ( ! start && Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1337 start = Standard_False;
1339 PrevOnBound = OnBound;
1341 // if found, adjust seam part
1342 if ( ind <= nbrPnt ) {
1343 PrevX = ( myAdjustOverDegen ? uf : ul );
1344 Standard_Real dU = Up/2 + Precision::PConfusion();
1345 if ( PrevOnBound ) {
1346 pnt2d(ind-1).SetX ( PrevX );
1347 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1348 Standard_Real CurX = pnt2d(j).X();
1349 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1350 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1353 else if ( OnBound ) {
1354 pnt2d(ind).SetX ( PrevX );
1355 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1356 Standard_Real CurX = pnt2d(j).X();
1357 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1358 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1361 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1365 else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1366 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1367 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1368 // 1st, find gap point (degenerated pole)
1369 Standard_Real PrevY=0.;
1370 Standard_Integer OnBound=0, PrevOnBound=0;
1371 Standard_Integer ind; // svv #1
1372 Standard_Boolean start = Standard_True;
1373 for ( ind=1; ind <= nbrPnt; ind++ ) {
1374 Standard_Real CurY = pnt2d(ind).Y();
1375 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1376 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1378 OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1379 Precision::PConfusion() );
1380 if ( ! start && Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1382 start = Standard_False;
1384 PrevOnBound = OnBound;
1386 // if found, adjust seam part
1387 if ( ind <= nbrPnt ) {
1388 PrevY = ( myAdjustOverDegen ? vf : vl );
1389 Standard_Real dV = Vp/2 + Precision::PConfusion();
1390 if ( PrevOnBound ) {
1391 pnt2d(ind-1).SetY ( PrevY );
1392 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1393 Standard_Real CurY = pnt2d(j).Y();
1394 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1395 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1398 else if ( OnBound ) {
1399 pnt2d(ind).SetY ( PrevY );
1400 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1401 Standard_Real CurY = pnt2d(j).Y();
1402 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1403 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1406 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1414 if(ChangeCycle) { // msv 10.08.04: avoid using of uninitialised field
1415 //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1416 // myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1417 myCashe3d[0] = points(1);
1418 myCashe3d[1] = points(nbrPnt);
1419 myCashe2d[0] = pnt2d(1);
1420 myCashe2d[1] = pnt2d(nbrPnt);
1423 myCashe3d[1] = points(1);
1424 myCashe3d[0] = points(nbrPnt);
1425 myCashe2d[1] = pnt2d(1);
1426 myCashe2d[0] = pnt2d(nbrPnt);
1431 //=======================================================================
1432 //function : ApproximatePCurve
1434 //=======================================================================
1436 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1437 Handle(TColgp_HArray1OfPnt2d)& points2d,
1438 Handle(TColStd_HArray1OfReal)& params,
1439 const Handle(Geom_Curve)& /*orig*/) const
1441 // Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1442 Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1443 Handle(Geom2d_Curve) C2d;
1446 CheckPoints2d (points2d, params, theTolerance2d);
1447 Standard_Integer numberPnt = points2d->Length();
1449 TColgp_Array1OfPnt points3d(1,numberPnt);
1452 Standard_Integer i; // svv #1
1453 for( i = 1; i <= numberPnt; i++) {
1454 pnt2d = points2d->Value(i);
1455 pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1459 GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1460 Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1461 Standard_Integer NbPoles = crv3d->NbPoles();
1462 TColgp_Array1OfPnt poles3d (1, NbPoles);
1463 TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1464 crv3d->Poles(poles3d);
1465 for( i = 1; i <= NbPoles; i++) {
1466 pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1469 TColStd_Array1OfReal weights (1,NbPoles);
1470 TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1471 TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1472 crv3d->Knots(knots);
1473 crv3d->Weights(weights);
1474 crv3d->Multiplicities(multiplicities);
1475 C2d = new Geom2d_BSplineCurve ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1478 catch(Standard_Failure) {
1479 #ifdef OCCT_DEBUG //:s5
1481 Standard_Integer nbp = params->Length();
1482 Standard_Integer nb2 = points2d->Length();
1483 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1484 Standard_Failure::Caught()->Print(cout);
1485 cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1486 // if (nb2 > nbp) nb2 = nbp;
1487 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1488 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1489 // dbl.AddReals (rb2,rbp);
1490 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1491 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1492 // dbl.AddXYZ (quoi);
1500 //=======================================================================
1501 //function : InterpolatePCurve
1503 //=======================================================================
1505 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1506 Handle(TColgp_HArray1OfPnt2d)& points2d,
1507 Handle(TColStd_HArray1OfReal)& params,
1508 const Handle(Geom_Curve)& /*orig*/) const
1510 Handle(Geom2d_Curve) C2d; // NULL si echec
1511 Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1514 // on verifie d abord s il n y a pas de points confondus
1515 // si besoin on retouche les valeurs ...
1516 CheckPoints2d (points2d, params, theTolerance2d);
1517 Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1518 Standard_False, theTolerance2d);
1519 myInterPol2d.Perform();
1520 if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1522 catch(Standard_Failure) {
1523 #ifdef OCCT_DEBUG //:s5
1525 Standard_Integer nbp = params->Length();
1526 Standard_Integer nb2 = points2d->Length();
1527 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1528 Standard_Failure::Caught()->Print(cout);
1529 cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1530 // if (nb2 > nbp) nb2 = nbp;
1531 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1532 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1533 // dbl.AddReals (rb2,rbp);
1534 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1535 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1536 // dbl.AddXYZ (quoi);
1544 //=======================================================================
1545 //function : InterpolateCurve3d
1547 //=======================================================================
1549 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1550 Handle(TColgp_HArray1OfPnt)& points,
1551 Handle(TColStd_HArray1OfReal)& params,
1552 const Handle(Geom_Curve)& /*orig*/) const
1554 Handle(Geom_Curve) C3d; // NULL si echec
1557 Standard_Real Tol = myPreci;
1558 CheckPoints(points, params, Tol);
1559 GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1560 myInterPol.Perform();
1561 if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1563 catch(Standard_Failure) {
1565 #ifdef OCCT_DEBUG //:s5
1566 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1567 Standard_Failure::Caught()->Print(cout); cout << endl;
1573 //=======================================================================
1574 //function : CheckPoints
1576 //=======================================================================
1578 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1580 Standard_Integer firstElem = points->Lower();
1581 Standard_Integer lastElem = points->Upper();
1583 Standard_Integer nbPntDropped = 0;
1584 Standard_Integer lastValid = firstElem; // indice of last undropped point
1586 // will store 0 when the point is to be removed, 1 otherwise
1587 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1588 for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
1589 Standard_Real DistMin2 = RealLast();
1590 gp_Pnt Prev = points->Value (lastValid);
1592 for (i = firstElem + 1; i <= lastElem ; i ++) {
1593 Curr = points->Value(i);
1594 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1595 if (CurDist2 < gp::Resolution()) { // test 0
1597 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1598 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1600 if (CurDist2 < DistMin2)
1601 DistMin2 = CurDist2;
1602 // lastValid becomes the current (i.e. i)
1607 if (DistMin2 < RealLast())
1608 preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1609 if (nbPntDropped == 0)
1613 cout << "Warning : removing 3d points for interpolation" << endl;
1615 // Build new HArrays
1616 Standard_Integer newLast = lastElem - nbPntDropped;
1617 if ((newLast - firstElem + 1) < 2) {
1619 cout << "Too many degenerated points for 3D interpolation" << endl;
1623 Handle(TColgp_HArray1OfPnt) newPnts =
1624 new TColgp_HArray1OfPnt(firstElem, newLast);
1625 Handle(TColStd_HArray1OfReal) newParams =
1626 new TColStd_HArray1OfReal(firstElem, newLast);
1627 Standard_Integer newCurr = 1;
1628 for (i = firstElem; i<= lastElem ; i++) {
1629 if (tmpParam.Value(i) == 1) {
1630 newPnts->SetValue(newCurr, points->Value(i));
1631 newParams->SetValue(newCurr, params->Value(i));
1637 // on la reduit un peu
1640 //=======================================================================
1641 //function : CheckPoints2d
1643 //=======================================================================
1645 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1646 Handle(TColStd_HArray1OfReal)& params,
1647 Standard_Real& preci) const
1649 Standard_Integer firstElem = points->Lower();
1650 Standard_Integer lastElem = points->Upper();
1652 Standard_Integer nbPntDropped = 0;
1653 Standard_Integer lastValid = firstElem; // indice of last undropped point
1655 // will store 0 when the point is to be removed, 1 otherwise
1656 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1657 for (i = firstElem; i<=lastElem ; i++) {
1658 tmpParam.SetValue(i,1);
1660 Standard_Real DistMin2 = RealLast();
1661 gp_Pnt2d Prev = points->Value(lastValid);
1663 for (i = firstElem + 1; i<=lastElem ; i++) {
1664 Curr = points->Value(i);
1665 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1666 if (CurDist2 < gp::Resolution()) { // test 0
1668 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1669 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1671 if (CurDist2 < DistMin2)
1672 DistMin2 = CurDist2;
1673 // lastValid becomes the current (i.e. i)
1678 if (DistMin2 < RealLast())
1679 preci = 0.9 * Sqrt (DistMin2);
1680 if (nbPntDropped == 0)
1684 cout << "Warning : removing 2d points for interpolation" << endl;
1686 // Build new HArrays
1687 Standard_Integer newLast = lastElem - nbPntDropped;
1688 if ((newLast - firstElem + 1) < 2) {
1690 cout << "Too many degenerated points for 2D interpolation" << endl;
1692 //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1693 tmpParam.SetValue(firstElem,1);
1694 tmpParam.SetValue(lastElem,1);
1695 gp_XY lastPnt = points->Value(lastElem).XY();
1696 lastPnt.Add(gp_XY(preci,preci));
1697 points->SetValue(lastElem,lastPnt);
1698 newLast = firstElem+1;
1701 Handle(TColgp_HArray1OfPnt2d) newPnts =
1702 new TColgp_HArray1OfPnt2d(firstElem, newLast);
1703 Handle(TColStd_HArray1OfReal) newParams =
1704 new TColStd_HArray1OfReal(firstElem, newLast);
1705 Standard_Integer newCurr = 1;
1706 for (i = firstElem; i <= lastElem ; i++) {
1707 if (tmpParam.Value(i) == 1) {
1709 cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1711 newPnts->SetValue(newCurr, points->Value(i));
1712 newParams->SetValue(newCurr, params->Value(i));
1717 cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1725 //=======================================================================
1726 //function : IsAnIsoparametric
1728 //=======================================================================
1729 //:S4030: modified for optimization
1730 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1732 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1733 const TColgp_Array1OfPnt& points,
1734 const TColStd_Array1OfReal& params,
1735 Standard_Boolean& isoTypeU,
1736 Standard_Boolean& p1OnIso,
1738 Standard_Boolean& p2OnIso,
1740 Standard_Boolean& isoPar2d3d,
1741 Handle(Geom_Curve)& cIso,
1744 TColStd_Array1OfReal& pout) const
1749 Standard_Real prec = Precision::Confusion();//myPreci;
1751 Standard_Boolean isoParam = Standard_False;
1752 isoPar2d3d = Standard_False;
1754 Standard_Real U1, U2, V1, V2;
1755 mySurf->Bounds(U1, U2, V1, V2);
1757 if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1758 Handle(Geom_RectangularTrimmedSurface) sTrim =
1759 Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1760 sTrim->Bounds(U1, U2, V1, V2);
1764 Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
1765 Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
1766 Standard_Real mindist2;
1767 Standard_Real mind2[2];
1768 mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1770 p1OnIso = Standard_False;
1771 p2OnIso = Standard_False;
1772 const Bnd_Box* aBox = 0;
1774 for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1775 Standard_Real isoVal=0.;
1776 Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1777 Handle(Geom_Curve) cI;
1778 Standard_Real tt1, tt2;
1781 if (Precision::IsInfinite(U1)) continue;
1782 cI = mySurf->UIso(U1);
1783 isoU = Standard_True;
1785 aBox = & mySurf->GetBoxUF();
1788 if (Precision::IsInfinite(U2)) continue;
1789 cI = mySurf->UIso(U2);
1790 isoU = Standard_True;
1792 aBox = & mySurf->GetBoxUL();
1795 if (Precision::IsInfinite(V1)) continue;
1796 cI = mySurf->VIso(V1);
1797 isoU = Standard_False;
1799 aBox = & mySurf->GetBoxVF();
1802 if (Precision::IsInfinite(V2)) continue;
1803 cI = mySurf->VIso(V2);
1804 isoU = Standard_False;
1806 aBox = & mySurf->GetBoxVL();
1811 if (isoU) { tt1 = V1; tt2 = V2; }
1812 else { tt1 = U1; tt2 = U2; }
1818 // PATCH CKY 9-JUL-1998 : protection contre singularite
1820 cI->D0( (tt1+tt2)/2,extmi);
1821 if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1823 Standard_Boolean PtEQext1 = Standard_False;
1824 Standard_Boolean PtEQext2 = Standard_False;
1826 Standard_Real currd2[2], tp[2] = {0, 0};
1827 Standard_Integer mp[2];
1829 for (Standard_Integer i=0; i<2; i++) {
1831 Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1833 // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1834 currd2[i] = points(k).SquareDistance ( ext1 );
1835 if ( currd2[i] <= prec*prec && !PtEQext1) {
1838 PtEQext1 = Standard_True;
1842 currd2[i] = points(k).SquareDistance ( ext2 );
1843 if ( currd2[i] <= prec*prec && !PtEQext2) {
1846 PtEQext2 = Standard_True;
1850 // On evite de projecter sur un iso degenere
1851 // on doit egalement le faire pour l apex du cone
1852 if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1856 if(aBox->IsOut(points(k))) continue;
1858 Standard_Real Cf = cI->FirstParameter();
1859 Standard_Real Cl = cI->LastParameter();
1860 if (Precision::IsInfinite(Cf)) Cf = -1000;
1861 if (Precision::IsInfinite(Cl)) Cl = +1000;
1863 ShapeAnalysis_Curve sac;
1864 Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1865 currd2[i] = dist * dist;
1866 if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1872 //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1873 // avoid possible null-length curves
1874 if ( mp[0] >0 && mp[1] >0 &&
1875 Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1879 ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1880 p1OnIso = Standard_True;
1881 mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
1882 if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1883 else valueP1.SetCoord(tp[0], isoVal);
1887 ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1888 p2OnIso = Standard_True;
1889 mind2[1] = currd2[1];
1890 if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1891 else valueP2.SetCoord(tp[1], isoVal);
1894 if ( mp[0] <=0 || mp[1] <=0 ) continue;
1896 Standard_Real md2 = currd2[0] + currd2[1];
1897 if ( mindist2 <= md2 ) continue;
1911 // probablely it concerns an isoparametrics
1912 if ( mpt[0] >0 && mpt[1] >0 ) {
1914 p1OnIso = p2OnIso = Standard_True;
1916 valueP1.SetCoord(isoValue, tpar[0]);
1917 valueP2.SetCoord(isoValue, tpar[1]);
1920 valueP1.SetCoord(tpar[0], isoValue);
1921 valueP2.SetCoord(tpar[1], isoValue);
1924 if ( mpt[0] != 3 && mpt[1] != 3 ) {
1925 isoPar2d3d = Standard_True;
1926 for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1927 if (tpar[1] > tpar[0]) t = params(i);
1928 else t = t1+t2-params(i);
1930 if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1934 if (isoPar2d3d) isoParam = Standard_True;
1936 Standard_Real prevParam = tpar[0];
1937 Standard_Real Cf, Cl;
1938 Standard_Boolean isoByDistance = Standard_True;
1939 Cf = cIso->FirstParameter();
1940 Cl = cIso->LastParameter();
1941 if (Precision::IsInfinite(Cf)) Cf = -1000;
1942 if (Precision::IsInfinite(Cl)) Cl = +1000;
1944 ShapeAnalysis_Curve sac;
1945 for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1946 Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1948 Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1951 if( (dist > prec) || (t < Cf) || (t > Cl) )
1952 isoByDistance = Standard_False;
1954 if (isoByDistance) isoParam = Standard_True;
1957 /* if (!isoParam) { CKY 29-mai-1997 : garder tout ce qu on peut ?
1958 p1OnIso = Standard_False;
1959 p2OnIso = Standard_False;
1963 catch(Standard_Failure) {
1964 // pb : on affiche ce qu on peut
1966 for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1967 cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1968 points(numpnt).X()<<" "<<points(numpnt).Y()<<" "<<points(numpnt).Z()<<")"<<endl;
1971 #ifdef OCCT_DEBUG //:s5
1972 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
1973 Standard_Failure::Caught()->Print(cout); cout << endl;
1975 return Standard_False;
1979 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
1980 //=======================================================================
1981 //function : BestExtremum
1982 //purpose : auxiliaire prenant le meilleur extremum si ISO car doute possible
1983 //=======================================================================
1985 gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
1987 // P2iso a ete calcule depuis P3ext sur une iso externe de la surface
1988 // En principe bon mais circularite possible ... et IsU/VClosed faillible
1989 // (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
1990 // 1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
1991 // 2/ si iso, attention a la circularite (cas limite)
1993 // NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
1995 // D abord, calcul p2ext depuis la surface. choix surface/iso
1997 Standard_Real prec = Precision::Confusion();//myPreci;
1998 gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
1999 gp_Pnt P3cal = mySurf->Value (P2cal);
2000 Standard_Real dcal = P3ext.Distance (P3cal);
2001 Standard_Real dnxt = P3ext.Distance (P3next);
2002 if (dcal > dnxt) return P2iso; // en fait protection sur BUG (PRO8468)
2004 // On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
2005 gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
2006 if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;