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_Curve.hxx>
67 #include <ShapeAnalysis_Surface.hxx>
68 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
69 #include <ShapeExtend.hxx>
70 #include <Standard_ErrorHandler.hxx>
71 #include <Standard_Failure.hxx>
72 #include <Standard_Type.hxx>
73 #include <TColgp_Array1OfPnt.hxx>
74 #include <TColStd_Array1OfInteger.hxx>
79 //=======================================================================
80 //function : ShapeConstruct_ProjectCurveOnSurface
82 //=======================================================================
84 ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
86 myPreci = Precision::Confusion();
87 myBuild = Standard_False;
88 myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
92 //=======================================================================
95 //=======================================================================
97 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
99 Init (new ShapeAnalysis_Surface (surf), preci);
102 //=======================================================================
105 //=======================================================================
107 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
110 SetPrecision (preci);
113 //=======================================================================
114 //function : SetSurface
116 //=======================================================================
118 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
120 SetSurface (new ShapeAnalysis_Surface (surf));
123 //=======================================================================
124 //function : SetSurface
126 //=======================================================================
128 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
130 if ( mySurf == surf ) return;
135 //=======================================================================
136 //function : SetPrecision
138 //=======================================================================
140 void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
145 //=======================================================================
146 //function : BuildCurveMode
148 //=======================================================================
150 Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
155 //=======================================================================
156 //function : AdjustOverDegenMode
158 //=======================================================================
161 //szv#4:S4163:12Mar99 was Boolean
162 Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
164 return myAdjustOverDegen;
168 //=======================================================================
169 //function : NbSurfIntervals
170 //purpose : work-around of bug in standard method
171 // GeomAdaptor_Surface->NbIntervals() (PRO16346)
172 //=======================================================================
174 static Standard_Integer NbSurfIntervals(const Handle(GeomAdaptor_HSurface)& GAS, const GeomAbs_Shape cont)
176 Standard_Integer NbU = 0;
177 if (GAS->GetType() == GeomAbs_SurfaceOfExtrusion) {
178 // extract the surface
179 Handle(Geom_SurfaceOfLinearExtrusion) surf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(GAS->ChangeSurface().Surface());
180 // build a 3d adaptor curve
181 GeomAdaptor_Curve Adaptor3dCurve(surf->BasisCurve(), GAS->FirstUParameter(), GAS->LastUParameter());
182 if (Adaptor3dCurve.GetType() == GeomAbs_BSplineCurve)
183 NbU = Adaptor3dCurve.NbIntervals(cont);
186 NbU = GAS->NbUIntervals(cont);
187 return NbU * (GAS->NbVIntervals(cont));
190 //=======================================================================
193 //=======================================================================
195 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status Status) const
197 return ShapeExtend::DecodeStatus (myStatus, Status);
200 //=======================================================================
203 //=======================================================================
204 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
205 const Standard_Real First,
206 const Standard_Real Last,
207 Handle(Geom2d_Curve)& c2d,
209 const Standard_Integer,
210 const Standard_Integer)
212 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
213 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
215 if (mySurf.IsNull()) {
217 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
218 return Standard_False;
220 // Projection Analytique
221 Handle(Geom_Curve) crv3dtrim = c3d;
222 if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
223 crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
224 c2d = ProjectAnalytic ( crv3dtrim );
226 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
227 return Standard_True;
230 // Projection par approximation
232 // discretize the 3d curve
234 Standard_Integer nbrPnt;
236 // $$$$ :92 abv 28 Jan 98 see PRO10107, big BSplineCurve C0
237 Standard_Integer nbPini = NCONTROL; // as in BRepCheck_Edge (RLN/Nijni)
238 // 20; // number of points for interpolation, should be "parametric dependent"
240 //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
241 // increase number of points to provide at least Degree()+1 points per interval
242 Handle(Geom_BSplineCurve) bspl;
243 if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
244 Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
245 bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
247 else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
248 if ( ! bspl.IsNull() ) {
249 Standard_Integer nint = 0;
250 for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
251 if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last ) nint++;
252 Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
253 while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
255 if ( nbPini > NCONTROL )
256 cout << "Warning: number of points for projecting is " << nbPini << endl;
260 // $$$$ end :92 (big BSplineCurve C0)
262 // this number should be "parametric dependent"
263 TColgp_Array1OfPnt points(1, nbPini);
264 TColStd_Array1OfReal params(1, nbPini);
265 NCollection_Sequence<Standard_Real> aKnotCoeffs;
267 Standard_Integer iPnt;
269 // In case of bspline compute parametrization speed on each
270 // knot interval inside [aFirstParam, aLastParam].
271 // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
272 // use PerformByProjLib algorithm.
275 Standard_Integer aNbIntPnts = NCONTROL;
276 Standard_Real aFirstParam = First; // First parameter of current interval.
277 Standard_Real aLastParam = Last; // Last parameter of current interval.
279 // First index computation.
280 Standard_Integer anIdx = 1;
281 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
283 if(bspl->Knot(anIdx) > First)
289 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
291 // Fill current knot interval.
292 aLastParam = Min(Last, bspl->Knot(anIdx));
293 Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
294 Standard_Integer anIntIdx;
296 Standard_Real aLength3d = 0.0;
297 for(anIntIdx = 0; anIntIdx < aNbIntPnts - 1; anIntIdx++)
299 // Start filling from first point.
300 Standard_Real aParam1 = aFirstParam + aStep * anIntIdx;
301 Standard_Real aParam2 = aFirstParam + aStep * (anIntIdx + 1);
302 c3d->D0 (aParam1, p3d1);
303 c3d->D0 (aParam2, p3d2);
304 aLength3d += p3d2.Distance(p3d1);
306 aKnotCoeffs.Append(aLength3d / (aLastParam - aFirstParam));
307 aFirstParam = aLastParam;
310 Standard_Real anEvenlyCoeff = 0;
311 if (aKnotCoeffs.Size() > 0)
313 anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
314 *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
317 const Standard_Real aMaxQuotientCoeff = 1500.0;
318 if (anEvenlyCoeff > aMaxQuotientCoeff)
320 PerformByProjLib(c3d, First, Last, c2d);
321 // PerformByProjLib fail detection:
324 return Status (ShapeExtend_DONE);
329 Standard_Real deltaT, t;
330 deltaT = (Last - First) / (nbPini-1);
332 for (iPnt = 1; iPnt <= nbPini; iPnt ++)
334 if (iPnt == 1) t = First;
335 else if (iPnt == nbPini) t = Last;
336 else t = First + (iPnt - 1) * deltaT;
343 // CALCUL par approximation
344 TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
345 ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
347 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
348 return Standard_True;
349 }// cas particulier d iso
351 // INTERPOLATION du resultat
354 Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt (1, nbPini);
355 Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
356 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
357 thePnts->SetValue(iPnt, points(iPnt));
358 theParams->SetValue(iPnt, params(iPnt));
361 Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
362 if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
364 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
369 Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
370 Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
371 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
372 theParams2d->SetValue(iPnt, params(iPnt));
373 thePnts2d->SetValue(iPnt, pnt2d(iPnt));
376 c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
377 // c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
378 // Faut-il aussi reprendre la C3D ?
379 myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
380 return Status (ShapeExtend_DONE);
383 //=======================================================================
384 //function : PerformByProjLib
386 //=======================================================================
388 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
389 const Standard_Real First,
390 const Standard_Real Last,
391 Handle(Geom2d_Curve)& c2d,
392 const GeomAbs_Shape /*continuity*/,
393 const Standard_Integer /*maxdeg */,
394 const Standard_Integer /*nbinterval */)
396 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
398 if (mySurf.IsNull()) {
399 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
400 return Standard_False;
406 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
407 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
408 ProjLib_ProjectedCurve Projector(GAS, GAC);
410 switch (Projector.GetType())
413 c2d = new Geom2d_Line(Projector.Line());
415 case GeomAbs_Circle :
416 c2d = new Geom2d_Circle(Projector.Circle());
418 case GeomAbs_Ellipse :
419 c2d = new Geom2d_Ellipse(Projector.Ellipse());
421 case GeomAbs_Parabola :
422 c2d = new Geom2d_Parabola(Projector.Parabola());
424 case GeomAbs_Hyperbola :
425 c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
427 case GeomAbs_BSplineCurve :
428 c2d = Projector.BSpline();
431 // Not possible, handling added to avoid gcc warning.
437 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
438 return Standard_False;
442 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
443 return Standard_True;
447 catch(Standard_Failure)
450 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
451 Standard_Failure::Caught()->Print(cout); cout << endl;
453 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
456 return Standard_False;
459 //=======================================================================
460 //function : PerformAdvanced
462 //=======================================================================
464 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
465 const Standard_Real First,
466 const Standard_Real Last,
467 Handle(Geom2d_Curve)& c2d)
469 Standard_Boolean hasResult = Standard_False;
470 Standard_Integer nbintervals;
472 Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
473 // && (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
475 if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
477 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
478 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
479 nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
480 isStandard = (nbintervals < 2);
483 hasResult = PerformByProjLib(c3d, First, Last, c2d);
485 if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
489 //=======================================================================
490 //function : ProjectAnalytic
492 //=======================================================================
494 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
496 Handle(Geom2d_Curve) result;
498 //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
499 // This is necessary for K4L since it fails on other surfaces
500 // when general method GeomProjLib::Curve2d() is used
501 // Projection is done as in BRep_Tool and BRepCheck_Edge
502 Handle(Geom_Surface) surf = mySurf->Surface();
503 Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
504 if ( Plane.IsNull() ) {
505 Handle(Geom_RectangularTrimmedSurface) RTS =
506 Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
507 if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
509 Handle(Geom_OffsetSurface) OS =
510 Handle(Geom_OffsetSurface)::DownCast ( surf );
512 Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
515 if ( ! Plane.IsNull() ) {
516 Handle(Geom_Curve) ProjOnPlane =
517 GeomProjLib::ProjectOnPlane (c3d, Plane,
518 Plane->Position().Direction(), Standard_True);
519 Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
520 ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
522 result = Geom2dAdaptor::MakeCurve(Proj);
523 if ( result.IsNull() ) return result;
524 if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
525 Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
526 result = TC->BasisCurve();
535 //! Fix possible period jump and handle walking period parameter.
536 static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
537 Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
538 Standard_Real thePeriod) // Period on objective coord
542 Standard_Boolean isNeedToFix = Standard_True;
543 for (i = 0; i < 3; i++)
545 Standard_Real aDiff = Abs (thePnt[i].Coord(theIdx) - thePnt[i + 1].Coord(theIdx));
546 if ( aDiff > Precision::PConfusion() &&
547 aDiff < thePeriod - Precision::PConfusion())
549 // Walk over period coord -> not walking on another isoline in parameter space.
550 isNeedToFix = Standard_False;
556 // Walking on isoline on another parameter. Fix period paramter to obtained minimum.
557 Standard_Real aFixParam = Min (thePnt[0].Coord(theIdx), thePnt[3].Coord(theIdx));
558 for(i = 0; i < 4; i++)
559 thePnt[i].SetCoord(theIdx, aFixParam);
562 // Fix possible period jump on first point.
563 if ( Abs(thePnt[0].Coord(theIdx) - thePnt[1].Coord(theIdx) ) > thePeriod / 2.01)
565 Standard_Real aMult = thePnt[0].Coord(theIdx) < thePnt[1].Coord(theIdx) ? 1.0 : -1.0;
566 Standard_Real aNewParam = thePnt[0].Coord(theIdx) + aMult * thePeriod;
567 thePnt[0].SetCoord(theIdx, aNewParam);
568 return Standard_False;
571 // Fix possible period jump on last point.
572 if ( Abs(thePnt[2].Coord(theIdx) - thePnt[3].Coord(theIdx) ) > thePeriod / 2.01)
574 Standard_Real aMult = thePnt[3].Coord(theIdx) < thePnt[2].Coord(theIdx) ? 1.0 : -1.0;
575 Standard_Real aNewParam = thePnt[3].Coord(theIdx) + aMult * thePeriod;
576 thePnt[3].SetCoord(theIdx, aNewParam);
577 return Standard_False;
580 return Standard_True;
583 //=======================================================================
586 //=======================================================================
588 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
589 const TColgp_Array1OfPnt& thepoints,
590 const TColStd_Array1OfReal& theparams,
591 TColgp_Array1OfPnt2d& thePnt2ds,
592 Standard_Real theTol,
593 Standard_Boolean &isRecompute) const
595 Standard_Integer nb = thepoints.Length();
597 aP[0] = thepoints(1);
598 aP[1] = thepoints(2);
599 aP[2] = thepoints(nb - 1);
600 aP[3] = thepoints(nb);
602 Standard_Integer i = 0;
604 Standard_Real aTol2 = theTol * theTol;
605 Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
606 Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
609 // Protection against bad "tolerance" shapes.
612 theTol = Precision::Confusion();
613 aTol2 = theTol * theTol;
615 Standard_Real anOldTol2 = aTol2;
617 // project first and last points
621 for ( j=0; j < myNbCashe; j++ )
622 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
624 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol,
628 if ( j >= myNbCashe )
629 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
631 Standard_Real aDist = mySurf->Gap();
632 Standard_Real aCurDist = aDist * aDist;
633 if( aTol2 < aDist * aDist)
637 if ( isPeriodicU || isPeriodicV )
639 // Compute second and last but one c2d points.
640 for(i = 1; i < 3; i++)
643 for ( j=0; j < myNbCashe; j++ )
644 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
646 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
649 if ( j >= myNbCashe )
650 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
652 Standard_Real aDist = mySurf->Gap();
653 Standard_Real aCurDist = aDist * aDist;
654 if( aTol2 < aDist * aDist)
660 isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod());
665 isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod());
669 thePnt2ds.SetValue(1, aP2d[0]);
670 thePnt2ds.SetValue(nb, aP2d[3]);
672 // Restore old tolerance in 2d space to avoid big gap cases.
674 // Check that straight line in 2d with parameterisation as in 3d will fit
675 // fit 3d curve at all points.
676 Standard_Real dPar = theparams(nb) - theparams(1);
677 if ( Abs(dPar) < Precision::PConfusion() )
679 gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
680 gp_Vec2d aVec = aVec0 / dPar;
681 Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
682 SquareDistance(thepoints(1));
683 for(i = 2; i < nb; i++)
685 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
687 mySurf->Surface()->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
688 Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
690 if(Abs (aFirstPointDist - aDist1) > aTol2)
694 // check if pcurve can be represented by Geom2d_Line (parameterised by length)
695 Standard_Real aLLength = aVec0.Magnitude();
696 if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
698 gp_XY aDirL = aVec0.XY() / aLLength;
699 gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
700 return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
703 // create straight bspline
704 TColgp_Array1OfPnt2d aPoles(1, 2);
708 TColStd_Array1OfReal aKnots(1,2);
709 aKnots(1) = theparams(1);
710 aKnots(2) = theparams(theparams.Length());
712 TColStd_Array1OfInteger aMults(1,2);
715 Standard_Integer aDegree = 1;
716 Handle(Geom2d_BSplineCurve) abspl2d =
717 new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
721 //=======================================================================
722 //function : ApproxPCurve
724 //=======================================================================
726 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
727 const TColgp_Array1OfPnt& points,
728 const TColStd_Array1OfReal& params,
729 TColgp_Array1OfPnt2d& pnt2d,
730 Handle(Geom2d_Curve)& c2d)
732 // for performance, first try to handle typical case when pcurve is straight
733 Standard_Boolean isRecompute = Standard_False;
734 c2d = getLine(points, params, pnt2d, myPreci, isRecompute);
737 return Standard_True;
739 Standard_Boolean isDone = Standard_True;
740 // test if the curve 3d is a boundary of the surface
741 // (only for Bezier or BSpline surface)
743 Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
744 gp_Pnt2d valueP1, valueP2;
745 Handle(Geom_Curve) cIso;
746 Standard_Real t1, t2;
748 Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
749 Standard_Boolean isAnalytic = Standard_True;
750 if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
751 Standard_Real uf, ul, vf, vl;
752 mySurf->Surface()->Bounds(uf, ul, vf, vl);
753 isoclosed = Standard_False;
754 TColStd_Array1OfReal pout(1, nbrPnt);
756 isoParam = IsAnIsoparametric(nbrPnt, points, params,
757 isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
758 isoPar2d3d, cIso, t1, t2, pout);
760 // projection of the points on surfaces
765 Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
766 Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
768 // Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
769 // Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
772 isoValue = valueP1.X();
773 isoPar1 = valueP1.Y();
774 isoPar2 = valueP2.Y();
775 isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
776 parf = vf; parl = vl;
779 isoValue = valueP1.Y();
780 isoPar1 = valueP1.X();
781 isoPar2 = valueP2.X();
782 isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
783 parf = uf; parl = ul;
785 if (!isoPar2d3d && !isAnalytic) {
786 Cf = cIso->FirstParameter();
787 Cl = cIso->LastParameter();
788 if (Precision::IsInfinite(Cf)) Cf = -1000;
789 if (Precision::IsInfinite(Cl)) Cl = +1000;
790 //pdn S4030 optimizing and fix isopar case on PRO41323
792 // dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
793 // Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
794 if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
795 if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
797 if (isoTypeU) valueP1.SetY (isoPar1);
798 else valueP1.SetX (isoPar1);
800 if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
801 //pdn S4030 optimizing and fix isopar case on PRO41323
802 tfin = pout(nbrPnt-1);
803 //dist = ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
804 if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
806 if (isoTypeU) valueP2.SetY (isoPar2);
807 else valueP2.SetX (isoPar2);
810 // Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
811 // Est-ce encore necessaire apres ce qui vient d etre fait ?
813 // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
814 // if curve is periodic do not change the points
815 // skl change "if" for pout(nbrPnt-1) 19.11.2003
817 if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
818 (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
819 gp_Pnt2d valueTmp = valueP1;
820 valueP1 = valueP2; valueP2 = valueTmp;
822 isoValue = valueP1.X();
823 isoPar1 = valueP1.Y();
824 isoPar2 = valueP2.Y();
827 isoValue = valueP1.Y();
828 isoPar1 = valueP1.X();
829 isoPar2 = valueP2.X();
831 // Fin calcul sens de courbe iso
833 } // end of fix check 05.02.02
837 // Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
838 // (utile pour detromper bug de projection). Mais detromper aussi circularite
840 //if (p1OnIso) valueP1 =
841 //BestExtremum (valueP1,points(1),points(2));
842 //if (p2OnIso) valueP2 =
843 //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
846 Standard_Real gap = myPreci; //:q1
847 Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
848 if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
849 //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
850 if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
851 ChangeCycle = Standard_True;
852 //for( i = 1; i <= nbrPnt; i ++) {
853 for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
854 if(ChangeCycle) //skl for OCC3430
862 if (isoPar2 > isoPar1) tPar = params(i);
863 else tPar = t1 + t2 - params(i);
864 } else if (!isAnalytic) {
866 if (i==1) tPar = isoPar1;
867 else if (i==nbrPnt) tPar = isoPar2;
870 //:S4030 ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
874 if (!isoPar2d3d && isAnalytic) {
875 if (i == 1) p2d = valueP1;
876 else if (i == nbrPnt) p2d = valueP2;
878 p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
879 Precision::Confusion()+1000*gap); //:q1
883 if(isoTypeU) { p2d.SetX(isoValue); p2d.SetY(tPar); }
884 else { p2d.SetX(tPar); p2d.SetY(isoValue); }
889 if ( (i == 1) && p1OnIso) p2d = valueP1;
890 else if( (i == nbrPnt) && p2OnIso) p2d = valueP2;
891 else {// general case (not an iso) mais attention aux singularites !
892 // first and last points are already computed by getLine()
893 if ( (i == 1 || i == nbrPnt))
903 //:q9 abv 23 Mar 99: use cashe as 1st approach
904 Standard_Integer j; // svv #1
905 for ( j=0; j < myNbCashe; j++ )
906 if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
908 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci,
909 Precision::Confusion()+gap);
912 if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci);
916 p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
917 Precision::Confusion()+1000*gap); //:q1
925 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
927 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
931 //pdn %12 11.02.99 PRO9234 entity 15402
933 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
934 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
937 // attention aux singularites ... (hors cas iso qui les traite deja)
940 // if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
942 // p2d = pnt2d (nbrPnt);
943 // if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
944 // pnt2d (nbrPnt) = p2d;
947 // Si la surface est UCLosed et VClosed, on recadre les points
948 // algo un peu complique, on retarde l implementation
949 Standard_Real Up = ul - uf;
950 Standard_Real Vp = vl - vf;
951 Standard_Real dist2d;
953 if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
954 cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
957 // Si la surface est UCLosed, on recadre les points
958 if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
959 // Premier point dans le domain [uf, ul]
960 Standard_Real prevX, firstX = pnt2d (1).X();
961 while (firstX < uf) { firstX += Up; pnt2d (1).SetX(firstX); }
962 while (firstX > ul) { firstX -= Up; pnt2d (1).SetX(firstX); }
965 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
966 Standard_Real minX = firstX, maxX = firstX;
968 // On decalle toujours le suivant
969 for (i = 2; i <= nbrPnt; i++) {
970 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
971 Standard_Real CurX = pnt2d (i).X();
972 dist2d = Abs (CurX - prevX);
973 if (dist2d > ( Up / 2) ) {
974 if (CurX > prevX + Up/2) {
975 while (CurX > prevX + Up/2) { CurX -= Up; pnt2d (i).SetX (CurX); }
976 } else if (CurX < prevX - Up/2) {
977 while (CurX < prevX - Up/2) { CurX += Up; pnt2d (i).SetX (CurX); }
982 if ( minX > CurX ) minX = CurX; //:97
983 else if ( maxX < CurX ) maxX = CurX; //:97
987 Standard_Real midX = 0.5 * ( minX + maxX );
988 Standard_Real shiftX=0.;
989 if ( midX > ul ) shiftX = -Up;
990 else if ( midX < uf ) shiftX = Up;
992 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
994 // Si la surface est VCLosed, on recadre les points
995 // Same code as UClosed : optimisation souhaitable !!
996 // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
997 // Si la surface est UCLosed, on recadre les points
999 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1000 //#78 rln 12.03.99 S4135
1001 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1002 // Premier point dans le domain [vf, vl]
1003 Standard_Real prevY, firstY = pnt2d (1).Y();
1004 while (firstY < vf) { firstY += Vp; pnt2d (1).SetY(firstY); }
1005 while (firstY > vl) { firstY -= Vp; pnt2d (1).SetY(firstY); }
1008 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1009 Standard_Real minY = firstY, maxY = firstY;
1011 // On decalle toujours le suivant
1012 for (i = 2; i <= nbrPnt; i ++) {
1013 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
1014 Standard_Real CurY = pnt2d (i).Y();
1015 dist2d = Abs (CurY - prevY);
1016 if (dist2d > ( Vp / 2) ) {
1017 if (CurY > prevY + Vp/2) {
1018 while (CurY > prevY + Vp/2) { CurY -= Vp; pnt2d (i).SetY (CurY); }
1019 } else if (CurY < prevY - Vp/2) {
1020 while (CurY < prevY - Vp/2) { CurY += Vp; pnt2d (i).SetY (CurY); }
1024 if ( minY > CurY ) minY = CurY; //:97
1025 else if ( maxY < CurY ) maxY = CurY; //:97
1029 Standard_Real midY = 0.5 * ( minY + maxY );
1030 Standard_Real shiftY=0.;
1031 if ( midY > vl ) shiftY = -Vp;
1032 else if ( midY < vf ) shiftY = Vp;
1034 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
1037 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1038 //#78 rln 12.03.99 S4135
1039 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1040 for (i = 2; i <= nbrPnt; i++) {
1041 //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1042 dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
1043 if (dist2d > ( Vp / 2) ) {
1044 // ATTENTION : il faut regarder ou le decalage se fait.
1045 // si plusieurs points sont decalles, il faut plusieurs passes
1046 // pour obtenir un resultat correct.
1047 // NOT YET IMPLEMENTED
1049 // one of those point is incorrectly placed
1050 // i.e on the wrong side of the "seam"
1051 // on prend le point le plus pres des bords vf ou vl
1052 Standard_Boolean prevOnFirst = Standard_False;
1053 Standard_Boolean prevOnLast = Standard_False;
1054 Standard_Boolean currOnFirst = Standard_False;
1055 Standard_Boolean currOnLast = Standard_False;
1057 // .X ? plutot .Y , non ?
1058 Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
1059 Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
1060 Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
1061 Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
1063 Standard_Real theMin = distPrevVF;
1064 prevOnFirst = Standard_True;
1065 if (distPrevVL < theMin) {
1066 theMin = distPrevVL;
1067 prevOnFirst = Standard_False;
1068 prevOnLast = Standard_True;
1070 if (distCurrVF < theMin) {
1071 theMin = distCurrVF;
1072 prevOnFirst = Standard_False;
1073 prevOnLast = Standard_False;
1074 currOnFirst = Standard_True;
1076 if (distCurrVL < theMin) {
1077 theMin = distCurrVL;
1078 prevOnFirst = Standard_False;
1079 prevOnLast = Standard_False;
1080 currOnFirst = Standard_False;
1081 currOnLast = Standard_True;
1083 // Modifs RLN/Nijni 3-DEC-1997
1085 // on decalle le point (i-1) en V Last
1086 gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of vl RLN/Nijni
1087 pnt2d (i-1) = newPrev;
1089 else if (prevOnLast) {
1090 // on decalle le point (i-1) en V first
1091 gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of vf RLN/Nijni
1092 pnt2d (i-1) = newPrev;
1094 else if (currOnFirst) {
1095 // on decalle le point (i) en V Last
1096 gp_Pnt2d newCurr(pnt2d (i).X(),vf); // instead of vl RLN/Nijni
1097 pnt2d (i) = newCurr;
1099 else if (currOnLast) {
1100 // on decalle le point (i) en V First
1101 gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf RLN/Nijni
1102 pnt2d (i) = newCurr;
1106 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1107 if (dist2d > ( Vp / 2) ) {
1108 cout << "Echec dans le recadrage" << endl;
1115 //:c0 abv 20 Feb 98: treat very special case when 3d curve
1116 // go over the pole of, e.g., sphere, and partly lies along seam.
1117 // 2d representation of such a curve should consist of 3 parts - one on
1118 // regular part of surface (interior), one part along degenerated boundary
1119 // and one along seam.
1120 // Since it cannot be adjusted later by arranging pcurves (curve is single),
1121 // to fix it it is nesessary to have a possibility of adjusting seam
1122 // part of such curve either to left or right boundary of surface.
1123 // Test is performed only if flag AdjustOverDegen is not -1.
1124 // If AdjustOverDegen is True, seam part of curve is adjusted to
1125 // the left, and if False - to the right parametric boundary
1126 // If treated case is detected, flag DONE4 is set to status
1127 // NOTE: currently, precision is Precision::PConfusion() since it
1128 // is enough on encountered example
1129 // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1130 // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1131 if ( myAdjustOverDegen != -1 ) {
1132 if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1133 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1134 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1135 // 1st, find gap point (degenerated pole)
1136 Standard_Real PrevX=0.;
1137 Standard_Integer OnBound=0, PrevOnBound=0;
1138 Standard_Integer ind; // svv #1
1139 Standard_Boolean start = Standard_True;
1140 for ( ind=1; ind <= nbrPnt; ind++ ) {
1141 Standard_Real CurX = pnt2d(ind).X();
1142 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1143 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1145 OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1146 Precision::PConfusion() );
1147 if ( ! start && Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1149 start = Standard_False;
1151 PrevOnBound = OnBound;
1153 // if found, adjust seam part
1154 if ( ind <= nbrPnt ) {
1155 PrevX = ( myAdjustOverDegen ? uf : ul );
1156 Standard_Real dU = Up/2 + Precision::PConfusion();
1157 if ( PrevOnBound ) {
1158 pnt2d(ind-1).SetX ( PrevX );
1159 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1160 Standard_Real CurX = pnt2d(j).X();
1161 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1162 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1165 else if ( OnBound ) {
1166 pnt2d(ind).SetX ( PrevX );
1167 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1168 Standard_Real CurX = pnt2d(j).X();
1169 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1170 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1173 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1177 else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1178 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1179 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1180 // 1st, find gap point (degenerated pole)
1181 Standard_Real PrevY=0.;
1182 Standard_Integer OnBound=0, PrevOnBound=0;
1183 Standard_Integer ind; // svv #1
1184 Standard_Boolean start = Standard_True;
1185 for ( ind=1; ind <= nbrPnt; ind++ ) {
1186 Standard_Real CurY = pnt2d(ind).Y();
1187 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1188 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1190 OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1191 Precision::PConfusion() );
1192 if ( ! start && Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1194 start = Standard_False;
1196 PrevOnBound = OnBound;
1198 // if found, adjust seam part
1199 if ( ind <= nbrPnt ) {
1200 PrevY = ( myAdjustOverDegen ? vf : vl );
1201 Standard_Real dV = Vp/2 + Precision::PConfusion();
1202 if ( PrevOnBound ) {
1203 pnt2d(ind-1).SetY ( PrevY );
1204 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1205 Standard_Real CurY = pnt2d(j).Y();
1206 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1207 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1210 else if ( OnBound ) {
1211 pnt2d(ind).SetY ( PrevY );
1212 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1213 Standard_Real CurY = pnt2d(j).Y();
1214 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1215 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1218 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1226 if(ChangeCycle) { // msv 10.08.04: avoid using of uninitialised field
1227 //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1228 // myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1229 myCashe3d[0] = points(1);
1230 myCashe3d[1] = points(nbrPnt);
1231 myCashe2d[0] = pnt2d(1);
1232 myCashe2d[1] = pnt2d(nbrPnt);
1235 myCashe3d[1] = points(1);
1236 myCashe3d[0] = points(nbrPnt);
1237 myCashe2d[1] = pnt2d(1);
1238 myCashe2d[0] = pnt2d(nbrPnt);
1243 //=======================================================================
1244 //function : ApproximatePCurve
1246 //=======================================================================
1248 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1249 Handle(TColgp_HArray1OfPnt2d)& points2d,
1250 Handle(TColStd_HArray1OfReal)& params,
1251 const Handle(Geom_Curve)& /*orig*/) const
1253 // Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1254 Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1255 Handle(Geom2d_Curve) C2d;
1258 CheckPoints2d (points2d, params, theTolerance2d);
1259 Standard_Integer numberPnt = points2d->Length();
1261 TColgp_Array1OfPnt points3d(1,numberPnt);
1264 Standard_Integer i; // svv #1
1265 for( i = 1; i <= numberPnt; i++) {
1266 pnt2d = points2d->Value(i);
1267 pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1271 GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1272 Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1273 Standard_Integer NbPoles = crv3d->NbPoles();
1274 TColgp_Array1OfPnt poles3d (1, NbPoles);
1275 TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1276 crv3d->Poles(poles3d);
1277 for( i = 1; i <= NbPoles; i++) {
1278 pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1281 TColStd_Array1OfReal weights (1,NbPoles);
1282 TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1283 TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1284 crv3d->Knots(knots);
1285 crv3d->Weights(weights);
1286 crv3d->Multiplicities(multiplicities);
1287 C2d = new Geom2d_BSplineCurve ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1290 catch(Standard_Failure) {
1291 #ifdef OCCT_DEBUG //:s5
1293 Standard_Integer nbp = params->Length();
1294 Standard_Integer nb2 = points2d->Length();
1295 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1296 Standard_Failure::Caught()->Print(cout);
1297 cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1298 // if (nb2 > nbp) nb2 = nbp;
1299 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1300 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1301 // dbl.AddReals (rb2,rbp);
1302 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1303 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1304 // dbl.AddXYZ (quoi);
1312 //=======================================================================
1313 //function : InterpolatePCurve
1315 //=======================================================================
1317 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1318 Handle(TColgp_HArray1OfPnt2d)& points2d,
1319 Handle(TColStd_HArray1OfReal)& params,
1320 const Handle(Geom_Curve)& /*orig*/) const
1322 Handle(Geom2d_Curve) C2d; // NULL si echec
1323 Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1326 // on verifie d abord s il n y a pas de points confondus
1327 // si besoin on retouche les valeurs ...
1328 CheckPoints2d (points2d, params, theTolerance2d);
1329 Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1330 Standard_False, theTolerance2d);
1331 myInterPol2d.Perform();
1332 if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1334 catch(Standard_Failure) {
1335 #ifdef OCCT_DEBUG //:s5
1337 Standard_Integer nbp = params->Length();
1338 Standard_Integer nb2 = points2d->Length();
1339 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1340 Standard_Failure::Caught()->Print(cout);
1341 cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1342 // if (nb2 > nbp) nb2 = nbp;
1343 // Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1344 // // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1345 // dbl.AddReals (rb2,rbp);
1346 // for (Standard_Integer i = 1; i <= nb2; i ++) {
1347 // gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1348 // dbl.AddXYZ (quoi);
1356 //=======================================================================
1357 //function : InterpolateCurve3d
1359 //=======================================================================
1361 Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1362 Handle(TColgp_HArray1OfPnt)& points,
1363 Handle(TColStd_HArray1OfReal)& params,
1364 const Handle(Geom_Curve)& /*orig*/) const
1366 Handle(Geom_Curve) C3d; // NULL si echec
1369 Standard_Real Tol = myPreci;
1370 CheckPoints(points, params, Tol);
1371 GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1372 myInterPol.Perform();
1373 if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1375 catch(Standard_Failure) {
1377 #ifdef OCCT_DEBUG //:s5
1378 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1379 Standard_Failure::Caught()->Print(cout); cout << endl;
1385 //=======================================================================
1386 //function : CheckPoints
1388 //=======================================================================
1390 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1392 Standard_Integer firstElem = points->Lower();
1393 Standard_Integer lastElem = points->Upper();
1395 Standard_Integer nbPntDropped = 0;
1396 Standard_Integer lastValid = firstElem; // indice of last undropped point
1398 // will store 0 when the point is to be removed, 1 otherwise
1399 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1400 for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
1401 Standard_Real DistMin2 = RealLast();
1402 gp_Pnt Prev = points->Value (lastValid);
1404 for (i = firstElem + 1; i <= lastElem ; i ++) {
1405 Curr = points->Value(i);
1406 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1407 if (CurDist2 < gp::Resolution()) { // test 0
1409 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1410 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1412 if (CurDist2 < DistMin2)
1413 DistMin2 = CurDist2;
1414 // lastValid becomes the current (i.e. i)
1419 if (DistMin2 < RealLast())
1420 preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1421 if (nbPntDropped == 0)
1425 cout << "Warning : removing 3d points for interpolation" << endl;
1427 // Build new HArrays
1428 Standard_Integer newLast = lastElem - nbPntDropped;
1429 if ((newLast - firstElem + 1) < 2) {
1431 cout << "Too many degenerated points for 3D interpolation" << endl;
1435 Handle(TColgp_HArray1OfPnt) newPnts =
1436 new TColgp_HArray1OfPnt(firstElem, newLast);
1437 Handle(TColStd_HArray1OfReal) newParams =
1438 new TColStd_HArray1OfReal(firstElem, newLast);
1439 Standard_Integer newCurr = 1;
1440 for (i = firstElem; i<= lastElem ; i++) {
1441 if (tmpParam.Value(i) == 1) {
1442 newPnts->SetValue(newCurr, points->Value(i));
1443 newParams->SetValue(newCurr, params->Value(i));
1449 // on la reduit un peu
1452 //=======================================================================
1453 //function : CheckPoints2d
1455 //=======================================================================
1457 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1458 Handle(TColStd_HArray1OfReal)& params,
1459 Standard_Real& preci) const
1461 Standard_Integer firstElem = points->Lower();
1462 Standard_Integer lastElem = points->Upper();
1464 Standard_Integer nbPntDropped = 0;
1465 Standard_Integer lastValid = firstElem; // indice of last undropped point
1467 // will store 0 when the point is to be removed, 1 otherwise
1468 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1469 for (i = firstElem; i<=lastElem ; i++) {
1470 tmpParam.SetValue(i,1);
1472 Standard_Real DistMin2 = RealLast();
1473 gp_Pnt2d Prev = points->Value(lastValid);
1475 for (i = firstElem + 1; i<=lastElem ; i++) {
1476 Curr = points->Value(i);
1477 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1478 if (CurDist2 < gp::Resolution()) { // test 0
1480 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1481 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1483 if (CurDist2 < DistMin2)
1484 DistMin2 = CurDist2;
1485 // lastValid becomes the current (i.e. i)
1490 if (DistMin2 < RealLast())
1491 preci = 0.9 * Sqrt (DistMin2);
1492 if (nbPntDropped == 0)
1496 cout << "Warning : removing 2d points for interpolation" << endl;
1498 // Build new HArrays
1499 Standard_Integer newLast = lastElem - nbPntDropped;
1500 if ((newLast - firstElem + 1) < 2) {
1502 cout << "Too many degenerated points for 2D interpolation" << endl;
1504 //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1505 tmpParam.SetValue(firstElem,1);
1506 tmpParam.SetValue(lastElem,1);
1507 gp_XY lastPnt = points->Value(lastElem).XY();
1508 lastPnt.Add(gp_XY(preci,preci));
1509 points->SetValue(lastElem,lastPnt);
1510 newLast = firstElem+1;
1513 Handle(TColgp_HArray1OfPnt2d) newPnts =
1514 new TColgp_HArray1OfPnt2d(firstElem, newLast);
1515 Handle(TColStd_HArray1OfReal) newParams =
1516 new TColStd_HArray1OfReal(firstElem, newLast);
1517 Standard_Integer newCurr = 1;
1518 for (i = firstElem; i <= lastElem ; i++) {
1519 if (tmpParam.Value(i) == 1) {
1521 cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1523 newPnts->SetValue(newCurr, points->Value(i));
1524 newParams->SetValue(newCurr, params->Value(i));
1529 cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1537 //=======================================================================
1538 //function : IsAnIsoparametric
1540 //=======================================================================
1541 //:S4030: modified for optimization
1542 //:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1544 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1545 const TColgp_Array1OfPnt& points,
1546 const TColStd_Array1OfReal& params,
1547 Standard_Boolean& isoTypeU,
1548 Standard_Boolean& p1OnIso,
1550 Standard_Boolean& p2OnIso,
1552 Standard_Boolean& isoPar2d3d,
1553 Handle(Geom_Curve)& cIso,
1556 TColStd_Array1OfReal& pout) const
1561 Standard_Real prec = Precision::Confusion();//myPreci;
1563 Standard_Boolean isoParam = Standard_False;
1564 isoPar2d3d = Standard_False;
1566 Standard_Real U1, U2, V1, V2;
1567 mySurf->Bounds(U1, U2, V1, V2);
1569 if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1570 Handle(Geom_RectangularTrimmedSurface) sTrim =
1571 Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1572 sTrim->Bounds(U1, U2, V1, V2);
1576 Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
1577 Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
1578 Standard_Real mindist2;
1579 Standard_Real mind2[2];
1580 mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1582 p1OnIso = Standard_False;
1583 p2OnIso = Standard_False;
1584 const Bnd_Box* aBox = 0;
1586 for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1587 Standard_Real isoVal=0.;
1588 Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1589 Handle(Geom_Curve) cI;
1590 Standard_Real tt1, tt2;
1593 if (Precision::IsInfinite(U1)) continue;
1594 cI = mySurf->UIso(U1);
1595 isoU = Standard_True;
1597 aBox = & mySurf->GetBoxUF();
1600 if (Precision::IsInfinite(U2)) continue;
1601 cI = mySurf->UIso(U2);
1602 isoU = Standard_True;
1604 aBox = & mySurf->GetBoxUL();
1607 if (Precision::IsInfinite(V1)) continue;
1608 cI = mySurf->VIso(V1);
1609 isoU = Standard_False;
1611 aBox = & mySurf->GetBoxVF();
1614 if (Precision::IsInfinite(V2)) continue;
1615 cI = mySurf->VIso(V2);
1616 isoU = Standard_False;
1618 aBox = & mySurf->GetBoxVL();
1623 if (isoU) { tt1 = V1; tt2 = V2; }
1624 else { tt1 = U1; tt2 = U2; }
1630 // PATCH CKY 9-JUL-1998 : protection contre singularite
1632 cI->D0( (tt1+tt2)/2,extmi);
1633 if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1635 Standard_Boolean PtEQext1 = Standard_False;
1636 Standard_Boolean PtEQext2 = Standard_False;
1638 Standard_Real currd2[2], tp[2] = {0, 0};
1639 Standard_Integer mp[2];
1641 for (Standard_Integer i=0; i<2; i++) {
1643 Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1645 // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1646 currd2[i] = points(k).SquareDistance ( ext1 );
1647 if ( currd2[i] <= prec*prec && !PtEQext1) {
1650 PtEQext1 = Standard_True;
1654 currd2[i] = points(k).SquareDistance ( ext2 );
1655 if ( currd2[i] <= prec*prec && !PtEQext2) {
1658 PtEQext2 = Standard_True;
1662 // On evite de projecter sur un iso degenere
1663 // on doit egalement le faire pour l apex du cone
1664 if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1668 if(aBox->IsOut(points(k))) continue;
1670 Standard_Real Cf = cI->FirstParameter();
1671 Standard_Real Cl = cI->LastParameter();
1672 if (Precision::IsInfinite(Cf)) Cf = -1000;
1673 if (Precision::IsInfinite(Cl)) Cl = +1000;
1675 ShapeAnalysis_Curve sac;
1676 Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1677 currd2[i] = dist * dist;
1678 if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1684 //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1685 // avoid possible null-length curves
1686 if ( mp[0] >0 && mp[1] >0 &&
1687 Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1691 ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1692 p1OnIso = Standard_True;
1693 mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
1694 if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1695 else valueP1.SetCoord(tp[0], isoVal);
1699 ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1700 p2OnIso = Standard_True;
1701 mind2[1] = currd2[1];
1702 if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1703 else valueP2.SetCoord(tp[1], isoVal);
1706 if ( mp[0] <=0 || mp[1] <=0 ) continue;
1708 Standard_Real md2 = currd2[0] + currd2[1];
1709 if ( mindist2 <= md2 ) continue;
1723 // probablely it concerns an isoparametrics
1724 if ( mpt[0] >0 && mpt[1] >0 ) {
1726 p1OnIso = p2OnIso = Standard_True;
1728 valueP1.SetCoord(isoValue, tpar[0]);
1729 valueP2.SetCoord(isoValue, tpar[1]);
1732 valueP1.SetCoord(tpar[0], isoValue);
1733 valueP2.SetCoord(tpar[1], isoValue);
1736 if ( mpt[0] != 3 && mpt[1] != 3 ) {
1737 isoPar2d3d = Standard_True;
1738 for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1739 if (tpar[1] > tpar[0]) t = params(i);
1740 else t = t1+t2-params(i);
1742 if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1746 if (isoPar2d3d) isoParam = Standard_True;
1748 Standard_Real prevParam = tpar[0];
1749 Standard_Real Cf, Cl;
1750 Standard_Boolean isoByDistance = Standard_True;
1751 Cf = cIso->FirstParameter();
1752 Cl = cIso->LastParameter();
1753 if (Precision::IsInfinite(Cf)) Cf = -1000;
1754 if (Precision::IsInfinite(Cl)) Cl = +1000;
1756 ShapeAnalysis_Curve sac;
1757 for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1758 Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1760 Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1763 if( (dist > prec) || (t < Cf) || (t > Cl) )
1764 isoByDistance = Standard_False;
1766 if (isoByDistance) isoParam = Standard_True;
1769 /* if (!isoParam) { CKY 29-mai-1997 : garder tout ce qu on peut ?
1770 p1OnIso = Standard_False;
1771 p2OnIso = Standard_False;
1775 catch(Standard_Failure) {
1776 // pb : on affiche ce qu on peut
1778 for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1779 cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1780 points(numpnt).X()<<" "<<points(numpnt).Y()<<" "<<points(numpnt).Z()<<")"<<endl;
1783 #ifdef OCCT_DEBUG //:s5
1784 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
1785 Standard_Failure::Caught()->Print(cout); cout << endl;
1787 return Standard_False;
1791 /* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
1792 //=======================================================================
1793 //function : BestExtremum
1794 //purpose : auxiliaire prenant le meilleur extremum si ISO car doute possible
1795 //=======================================================================
1797 gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
1799 // P2iso a ete calcule depuis P3ext sur une iso externe de la surface
1800 // En principe bon mais circularite possible ... et IsU/VClosed faillible
1801 // (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
1802 // 1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
1803 // 2/ si iso, attention a la circularite (cas limite)
1805 // NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
1807 // D abord, calcul p2ext depuis la surface. choix surface/iso
1809 Standard_Real prec = Precision::Confusion();//myPreci;
1810 gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
1811 gp_Pnt P3cal = mySurf->Value (P2cal);
1812 Standard_Real dcal = P3ext.Distance (P3cal);
1813 Standard_Real dnxt = P3ext.Distance (P3next);
1814 if (dcal > dnxt) return P2iso; // en fait protection sur BUG (PRO8468)
1816 // On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
1817 gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
1818 if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;