0022598: Approximation of p-curve by 2D line
[occt.git] / src / ShapeConstruct / ShapeConstruct_ProjectCurveOnSurface.cxx
CommitLineData
973c2be1 1// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 2//
973c2be1 3// This file is part of Open CASCADE Technology software library.
b311480e 4//
d5f74e42 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
973c2be1 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.
b311480e 10//
973c2be1 11// Alternatively, this file may be used under the terms of Open CASCADE
12// commercial license or contractual agreement.
b311480e 13
7fd59977 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)
28//szv#4 S4163
29//:s5 abv 22.04.99 Adding debug printouts in catch {} blocks
30//#1 svv 11.01.00 Porting on DEC
31#include <ShapeConstruct_ProjectCurveOnSurface.ixx>
32
33#include <Standard_ErrorHandler.hxx>
34#include <Standard_Failure.hxx>
35
36#include <Precision.hxx>
37#include <TColStd_Array1OfInteger.hxx>
38#include <TColgp_Array1OfPnt.hxx>
39
40#include <GeomAPI_PointsToBSpline.hxx>
41#include <Geom2dAPI_Interpolate.hxx>
42#include <GeomAPI_Interpolate.hxx>
43#include <Geom2dAdaptor.hxx>
44#include <Geom2d_Line.hxx>
45#include <Geom2d_TrimmedCurve.hxx>
46#include <Geom2d_BSplineCurve.hxx>
47#include <Geom_BSplineCurve.hxx>
48#include <Geom_TrimmedCurve.hxx>
49#include <Geom_RectangularTrimmedSurface.hxx>
50#include <Geom_SurfaceOfLinearExtrusion.hxx>
51#include <Geom_SphericalSurface.hxx>
52#include <Geom_OffsetSurface.hxx>
53#include <Geom_Plane.hxx>
54#include <GeomProjLib.hxx>
55#include <GeomAdaptor_HSurface.hxx>
56#include <GeomAdaptor_HCurve.hxx>
57
58#include <ShapeAnalysis_Curve.hxx>
59#include <ShapeAnalysis_Surface.hxx>
60#include <ShapeExtend.hxx>
61
62#include <ProjLib_ProjectedCurve.hxx>
63#include <ProjLib_CompProjectedCurve.hxx>
64#include <ProjLib_HCompProjectedCurve.hxx>
65#include <Approx_CurveOnSurface.hxx>
66
67#define NCONTROL 23
b311480e 68
7fd59977 69//=======================================================================
70//function : ShapeConstruct_ProjectCurveOnSurface
71//purpose :
72//=======================================================================
73
74ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
75{
76 myPreci = Precision::Confusion();
77 myBuild = Standard_False;
78 myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
79 myNbCashe = 0; //:q9
80}
81
82//=======================================================================
83//function : Init
84//purpose :
85//=======================================================================
86
87 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
88{
89 Init (new ShapeAnalysis_Surface (surf), preci);
90}
91
92//=======================================================================
93//function : Init
94//purpose :
95//=======================================================================
96
97 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
98{
99 SetSurface (surf);
100 SetPrecision (preci);
101}
102
103//=======================================================================
104//function : SetSurface
105//purpose :
106//=======================================================================
107
108 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
109{
110 SetSurface (new ShapeAnalysis_Surface (surf));
111}
112
113//=======================================================================
114//function : SetSurface
115//purpose :
116//=======================================================================
117
118 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
119{
120 if ( mySurf == surf ) return;
121 mySurf = surf;
122 myNbCashe = 0; //:q9
123}
124
125//=======================================================================
126//function : SetPrecision
127//purpose :
128//=======================================================================
129
130 void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
131{
132 myPreci = preci;
133}
134
135//=======================================================================
136//function : BuildCurveMode
137//purpose :
138//=======================================================================
139
140 Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
141{
142 return myBuild;
143}
144
145//=======================================================================
146//function : AdjustOverDegenMode
147//purpose :
148//=======================================================================
149//:c0
150
151//szv#4:S4163:12Mar99 was Boolean
152 Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
153 {
154 return myAdjustOverDegen;
155 }
156
157
158//=======================================================================
159//function : NbSurfIntervals
160//purpose : work-around of bug in standard method
161// GeomAdaptor_Surface->NbIntervals() (PRO16346)
162//=======================================================================
163
164static Standard_Integer NbSurfIntervals(const Handle(GeomAdaptor_HSurface)& GAS, const GeomAbs_Shape cont)
165{
166 Standard_Integer NbU = 0;
167 if (GAS->GetType() == GeomAbs_SurfaceOfExtrusion) {
168 // extract the surface
169 Handle(Geom_SurfaceOfLinearExtrusion) surf = Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(GAS->ChangeSurface().Surface());
170 // build a 3d adaptor curve
171 GeomAdaptor_Curve Adaptor3dCurve(surf->BasisCurve(), GAS->FirstUParameter(), GAS->LastUParameter());
172 if (Adaptor3dCurve.GetType() == GeomAbs_BSplineCurve)
173 NbU = Adaptor3dCurve.NbIntervals(cont);
174 }
175 if (NbU == 0)
176 NbU = GAS->NbUIntervals(cont);
177 return NbU * (GAS->NbVIntervals(cont));
178}
179
180//=======================================================================
181//function : Status
182//purpose :
183//=======================================================================
184
185 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status Status) const
186{
187 return ShapeExtend::DecodeStatus (myStatus, Status);
188}
189
190//=======================================================================
191//function : Perform
192//purpose :
193//=======================================================================
194
195Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
196 const Standard_Real First,
197 const Standard_Real Last,
198 Handle(Geom2d_Curve)& c2d,
199 const GeomAbs_Shape,
200 const Standard_Integer,
201 const Standard_Integer)
202{
203 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
204 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
205
206 if (mySurf.IsNull()) {
207 c2d.Nullify();
208 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
209 return Standard_False;
210 }
211
212// Projection Analytique
213 Handle(Geom_Curve) crv3dtrim = c3d;
214 if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
215 crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
216 c2d = ProjectAnalytic ( crv3dtrim );
217 if (!c2d.IsNull()) {
218 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
219 return Standard_True;
220 }
221
222// Projection par approximation
223
224 // discretize the 3d curve
225
226 Standard_Integer nbrPnt;
227
228// $$$$ :92 abv 28 Jan 98 see PRO10107, big BSplineCurve C0
229 Standard_Integer nbPini = NCONTROL; // as in BRepCheck_Edge (RLN/Nijni)
230 // 20; // number of points for interpolation, should be "parametric dependent"
231
232 //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
233 // increase number of points to provide at least Degree()+1 points per interval
234 Handle(Geom_BSplineCurve) bspl;
235 if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
236 Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
237 bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
238 }
239 else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
240 if ( ! bspl.IsNull() ) {
241 Standard_Integer nint = 0;
242 for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
243 if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last ) nint++;
244 Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
245 while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
0797d9d3 246#ifdef OCCT_DEBUG
7fd59977 247 if ( nbPini > NCONTROL )
248 cout << "Warning: number of points for projecting is " << nbPini << endl;
249#endif
250 }
251
252// $$$$ end :92 (big BSplineCurve C0)
253
254 // this number should be "parametric dependent"
255 TColgp_Array1OfPnt points(1, nbPini);
256 TColStd_Array1OfReal params(1, nbPini);
257
258 Standard_Integer iPnt;
259 gp_Pnt p3d;
260 Standard_Real deltaT, t;
261 deltaT = (Last - First) / (nbPini-1);
262 nbrPnt = nbPini;
263 for (iPnt = 1; iPnt <= nbPini; iPnt ++) {
264 if (iPnt == 1) t = First;
265 else if (iPnt == nbPini) t = Last;
266 else t = First + (iPnt - 1) * deltaT;
267
268 c3d->D0 (t, p3d);
269 points(iPnt) = p3d;
270 params(iPnt) = t;
271 }
272
273// CALCUL par approximation
274
275 TColgp_Array1OfPnt2d pnt2d(1, nbrPnt);
276 ApproxPCurve (nbrPnt,points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
277 if (!c2d.IsNull()) {
278 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
279 return Standard_True;
280 }// cas particulier d iso
281
282// INTERPOLATION du resultat
283
284 if ( myBuild ) {
285 Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt (1, nbPini);
286 Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
287 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
288 thePnts->SetValue(iPnt, points(iPnt));
289 theParams->SetValue(iPnt, params(iPnt));
290 }
291
292 Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
293 if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
294 else {
295 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
296 c3d = newc3d;
297 }
298 }
299
300 Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
301 Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
302 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
303 theParams2d->SetValue(iPnt, params(iPnt));
304 thePnts2d->SetValue(iPnt, pnt2d(iPnt));
305 }
306
307 c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
308// c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
309// Faut-il aussi reprendre la C3D ?
7fd59977 310 myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
311 return Status (ShapeExtend_DONE);
312}
313
314//=======================================================================
315//function : PerformByProjLib
316//purpose :
317//=======================================================================
318
319Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
320 const Standard_Real First,
321 const Standard_Real Last,
322 Handle(Geom2d_Curve)& c2d,
323 const GeomAbs_Shape continuity,
324 const Standard_Integer maxdeg,
325 const Standard_Integer nbinterval)
326{
327 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
328 c2d.Nullify();
329 if (mySurf.IsNull()) {
330 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
331 return Standard_False;
332 }
333
334 try {
335 OCC_CATCH_SIGNALS
336 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
337 Standard_Real URes = GAS->ChangeSurface().UResolution ( myPreci );
338 Standard_Real VRes = GAS->ChangeSurface().VResolution ( myPreci );
339 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
340 ProjLib_CompProjectedCurve Projector ( GAS, GAC, URes, VRes );
341
342 Standard_Real ubeg, ufin;
343 Standard_Integer nbSol = Projector.NbCurves();
344 if (nbSol==1) {
345 Projector.Bounds(1, ubeg, ufin);
346 if((ubeg<=First)&&(ufin>=Last)) {
347 Standard_Integer nbintervals = ( nbinterval < 1 ?
348 NbSurfIntervals(GAS, GeomAbs_C3)+GAC->NbIntervals(GeomAbs_C3)+2:
349 nbinterval);
350 Handle(ProjLib_HCompProjectedCurve) HProjector = new ProjLib_HCompProjectedCurve();
351 HProjector->Set(Projector);
352 Handle(Adaptor2d_HCurve2d) HPCur = HProjector;
353 Approx_CurveOnSurface appr(HPCur, GAS, First, Last, myPreci,
354 continuity, maxdeg,
355 nbintervals,
356 Standard_False, Standard_True);
357 if ( appr.IsDone() )
358 c2d = appr.Curve2d();
359 }
0797d9d3 360#ifdef OCCT_DEBUG
7fd59977 361 else
362 cout<<"Warning: ProjLib cutting pcurve "<< First << " -> " << ubeg <<" ; "<< Last << " -> " << ufin << endl;
363#endif
364 }
0797d9d3 365#ifdef OCCT_DEBUG
7fd59977 366 else cout<<"Warning: ProjLib "<< nbSol << " curves in ProjLib"<<endl;
367#endif
368 if(c2d.IsNull()) {
369 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
370 return Standard_False;
371 }
372 else {
373 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
374 return Standard_True;
375 }
376
377 }
378 catch(Standard_Failure) {
0797d9d3 379#ifdef OCCT_DEBUG
7fd59977 380 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
381 Standard_Failure::Caught()->Print(cout); cout << endl;
382#endif
383 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
384 c2d.Nullify();
385 }
386 return Standard_False;
387}
388
389//=======================================================================
390//function : PerformAdvanced
391//purpose :
392//=======================================================================
393
394Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformAdvanced (Handle(Geom_Curve)& c3d,
395 const Standard_Real First,
396 const Standard_Real Last,
397 Handle(Geom2d_Curve)& c2d)
398{
399 Standard_Boolean hasResult = Standard_False;
400 Standard_Integer nbintervals;
401
402 Standard_Boolean isStandard = (mySurf->Adaptor3d()->GetType() != GeomAbs_Cylinder);
403// && (mySurf->Adaptor3d()->GetType() != GeomAbs_SurfaceOfRevolution);
404
405 if (isStandard) isStandard = !mySurf->HasSingularities(myPreci);
406 if (isStandard) {
407 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
408 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
409 nbintervals = NbSurfIntervals(GAS, GeomAbs_C1);//+GAC->NbIntervals(GeomAbs_C3);
410 isStandard = (nbintervals < 2);
411 }
412 if (isStandard) {
413 hasResult = PerformByProjLib(c3d, First, Last, c2d);
414 }
415 if (!hasResult) hasResult = Perform (c3d, First, Last, c2d);
416 return hasResult;
417}
418
419//=======================================================================
420//function : ProjectAnalytic
421//purpose :
422//=======================================================================
423
424 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
425{
426 Handle(Geom2d_Curve) result;
427
428 //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
429 // This is necessary for K4L since it fails on other surfaces
430 // when general method GeomProjLib::Curve2d() is used
431 // Projection is done as in BRep_Tool and BRepCheck_Edge
432 Handle(Geom_Surface) surf = mySurf->Surface();
433 Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
434 if ( Plane.IsNull() ) {
435 Handle(Geom_RectangularTrimmedSurface) RTS =
436 Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
437 if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
438 else {
439 Handle(Geom_OffsetSurface) OS =
440 Handle(Geom_OffsetSurface)::DownCast ( surf );
441 if ( ! OS.IsNull() )
442 Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
443 }
444 }
445 if ( ! Plane.IsNull() ) {
446 Handle(Geom_Curve) ProjOnPlane =
447 GeomProjLib::ProjectOnPlane (c3d, Plane,
448 Plane->Position().Direction(), Standard_True);
449 Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
450 ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
451
452 result = Geom2dAdaptor::MakeCurve(Proj);
453 if ( result.IsNull() ) return result;
454 if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
455 Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
456 result = TC->BasisCurve();
457 }
63c629aa 458
7fd59977 459 return result;
460 }
461
462 return result;
463}
464
15b54261 465 //! Fix possible period jump and handle walking period parameter.
466 static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
467 Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
468 Standard_Real thePeriod) // Period on objective coord
469 {
470 Standard_Integer i;
471
472 Standard_Boolean isNeedToFix = Standard_True;
473 for (i = 0; i < 3; i++)
474 {
475 Standard_Real aDiff = Abs (thePnt[i].Coord(theIdx) - thePnt[i + 1].Coord(theIdx));
476 if ( aDiff > Precision::PConfusion() &&
477 aDiff < thePeriod - Precision::PConfusion())
478 {
479 // Walk over period coord -> not walking on another isoline in parameter space.
480 isNeedToFix = Standard_False;
481 }
482 }
483
484 if (isNeedToFix)
485 {
486 // Walking on isoline on another parameter. Fix period paramter to obtained minimum.
487 Standard_Real aFixParam = Min (thePnt[0].Coord(theIdx), thePnt[3].Coord(theIdx));
488 for(i = 0; i < 4; i++)
489 thePnt[i].SetCoord(theIdx, aFixParam);
490 }
491
492 // Fix possible period jump on first point.
493 if ( Abs(thePnt[0].Coord(theIdx) - thePnt[1].Coord(theIdx) ) > thePeriod / 2.01)
494 {
495 Standard_Real aMult = thePnt[0].Coord(theIdx) < thePnt[1].Coord(theIdx) ? 1.0 : -1.0;
496 Standard_Real aNewParam = thePnt[0].Coord(theIdx) + aMult * thePeriod;
497 thePnt[0].SetCoord(theIdx, aNewParam);
498 return Standard_False;
499 }
500
501 // Fix possible period jump on last point.
502 if ( Abs(thePnt[2].Coord(theIdx) - thePnt[3].Coord(theIdx) ) > thePeriod / 2.01)
503 {
504 Standard_Real aMult = thePnt[3].Coord(theIdx) < thePnt[2].Coord(theIdx) ? 1.0 : -1.0;
505 Standard_Real aNewParam = thePnt[3].Coord(theIdx) + aMult * thePeriod;
506 thePnt[3].SetCoord(theIdx, aNewParam);
507 return Standard_False;
508 }
509
510 return Standard_True;
511 }
512
513//=======================================================================
514//function : getLine
515//purpose :
516//=======================================================================
517
518 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
519 const TColgp_Array1OfPnt& thepoints,
520 const TColStd_Array1OfReal& theparams,
521 TColgp_Array1OfPnt2d& thePnt2ds,
522 Standard_Real theTol,
523 Standard_Boolean &isRecompute) const
524 {
525 Standard_Integer nb = thepoints.Length();
526 gp_Pnt aP[4];
527 aP[0] = thepoints(1);
528 aP[1] = thepoints(2);
529 aP[2] = thepoints(nb - 1);
530 aP[3] = thepoints(nb);
531 gp_Pnt2d aP2d[4];
532 Standard_Integer i = 0;
533
534 Standard_Real aTol2 = theTol * theTol;
535 Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
536 Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
537
538 // Workaround:
539 // Protection against bad "tolerance" shapes.
540 if (aTol2 > 1.0)
541 {
542 theTol = Precision::Confusion();
543 aTol2 = theTol * theTol;
544 }
545 Standard_Real anOldTol2 = aTol2;
546
547 // project first and last points
548 for( ; i < 4; i +=3)
549 {
550 Standard_Integer j;
551 for ( j=0; j < myNbCashe; j++ )
552 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
553 {
554 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol,
555 theTol);
556 break;
557 }
558 if ( j >= myNbCashe )
559 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
560
561 Standard_Real aDist = mySurf->Gap();
562 Standard_Real aCurDist = aDist * aDist;
563 if( aTol2 < aDist * aDist)
564 aTol2 = aCurDist;
565 }
566
567 if ( isPeriodicU || isPeriodicV )
568 {
569 // Compute second and last but one c2d points.
570 for(i = 1; i < 3; i++)
571 {
572 Standard_Integer j;
573 for ( j=0; j < myNbCashe; j++ )
574 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
575 {
576 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
577 break;
578 }
579 if ( j >= myNbCashe )
580 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
581
582 Standard_Real aDist = mySurf->Gap();
583 Standard_Real aCurDist = aDist * aDist;
584 if( aTol2 < aDist * aDist)
585 aTol2 = aCurDist;
586 }
587
588 if (isPeriodicU)
589 {
590 isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod());
591 }
592
593 if (isPeriodicV)
594 {
595 isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod());
596 }
597 }
598
599 thePnt2ds.SetValue(1, aP2d[0]);
600 thePnt2ds.SetValue(nb, aP2d[3]);
601
602 // Restore old tolerance in 2d space to avoid big gap cases.
603 aTol2 = anOldTol2;
604 // Check that straight line in 2d with parameterisation as in 3d will fit
605 // fit 3d curve at all points.
606 Standard_Real dPar = theparams(nb) - theparams(1);
607 if ( Abs(dPar) < Precision::PConfusion() )
608 return 0;
609 gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
610 gp_Vec2d aVec = aVec0 / dPar;
611 Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
612 SquareDistance(thepoints(1));
613 for(i = 2; i < nb; i++)
614 {
615 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
616 gp_Pnt aCurP;
617 mySurf->Surface()->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
618 Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
619
620 if(Abs (aFirstPointDist - aDist1) > aTol2)
621 return 0;
622 }
623
624 // check if pcurve can be represented by Geom2d_Line (parameterised by length)
625 Standard_Real aLLength = aVec0.Magnitude();
626 if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
627 {
628 gp_XY aDirL = aVec0.XY() / aLLength;
629 gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
630 return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
631 }
632
633 // create straight bspline
634 TColgp_Array1OfPnt2d aPoles(1, 2);
635 aPoles(1) = aP2d[0];
636 aPoles(2) = aP2d[3];
637
638 TColStd_Array1OfReal aKnots(1,2);
639 aKnots(1) = theparams(1);
640 aKnots(2) = theparams(theparams.Length());
641
642 TColStd_Array1OfInteger aMults(1,2);
643 aMults(1) = 2;
644 aMults(2) = 2;
645 Standard_Integer aDegree = 1;
646 Handle(Geom2d_BSplineCurve) abspl2d =
647 new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
648 return abspl2d;
649 }
650
7fd59977 651//=======================================================================
652//function : ApproxPCurve
653//purpose :
654//=======================================================================
655
15b54261 656 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
7fd59977 657 const TColgp_Array1OfPnt& points,
658 const TColStd_Array1OfReal& params,
659 TColgp_Array1OfPnt2d& pnt2d,
660 Handle(Geom2d_Curve)& c2d)
661{
15b54261 662 // for performance, first try to handle typical case when pcurve is straight
663 Standard_Boolean isRecompute = Standard_False;
664 c2d = getLine(points, params, pnt2d, myPreci, isRecompute);
665 if(!c2d.IsNull())
666 {
667 return Standard_True;
668 }
669 Standard_Boolean isDone = Standard_True;
7fd59977 670 // test if the curve 3d is a boundary of the surface
671 // (only for Bezier or BSpline surface)
672
673 Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
674 gp_Pnt2d valueP1, valueP2;
675 Handle(Geom_Curve) cIso;
676 Standard_Real t1, t2;
677
678 Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
679 Standard_Boolean isAnalytic = Standard_True;
680 if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
681 Standard_Real uf, ul, vf, vl;
682 mySurf->Surface()->Bounds(uf, ul, vf, vl);
683 isoclosed = Standard_False;
684 TColStd_Array1OfReal pout(1, nbrPnt);
685
686 isoParam = IsAnIsoparametric(nbrPnt, points, params,
687 isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
688 isoPar2d3d, cIso, t1, t2, pout);
689
690 // projection of the points on surfaces
691
692 gp_Pnt p3d;
693 gp_Pnt2d p2d;
694 Standard_Integer i;
695 Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
696 Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
697
698 // Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
699 // Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
700 if(isoParam){
701 if(isoTypeU){
702 isoValue = valueP1.X();
703 isoPar1 = valueP1.Y();
704 isoPar2 = valueP2.Y();
705 isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
706 parf = vf; parl = vl;
707 }
708 else {
709 isoValue = valueP1.Y();
710 isoPar1 = valueP1.X();
711 isoPar2 = valueP2.X();
712 isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
713 parf = uf; parl = ul;
714 }
715 if (!isoPar2d3d && !isAnalytic) {
716 Cf = cIso->FirstParameter();
717 Cl = cIso->LastParameter();
718 if (Precision::IsInfinite(Cf)) Cf = -1000;
719 if (Precision::IsInfinite(Cl)) Cl = +1000;
720 //pdn S4030 optimizing and fix isopar case on PRO41323
721 tdeb = pout(2);
722 // dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
723 // Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
724 if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
725 if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
726 else isoPar1 = parl;
727 if (isoTypeU) valueP1.SetY (isoPar1);
728 else valueP1.SetX (isoPar1);
729 }
730 if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
731 //pdn S4030 optimizing and fix isopar case on PRO41323
732 tfin = pout(nbrPnt-1);
733 //dist = ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
734 if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
735 else isoPar2 = parl;
736 if (isoTypeU) valueP2.SetY (isoPar2);
737 else valueP2.SetX (isoPar2);
738 }
739
740 // Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
741 // Est-ce encore necessaire apres ce qui vient d etre fait ?
742
743 // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
744 // if curve is periodic do not change the points
745 // skl change "if" for pout(nbrPnt-1) 19.11.2003
746 if (!isoclosed) {
747 if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
748 (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
749 gp_Pnt2d valueTmp = valueP1;
750 valueP1 = valueP2; valueP2 = valueTmp;
751 if (isoTypeU) {
752 isoValue = valueP1.X();
753 isoPar1 = valueP1.Y();
754 isoPar2 = valueP2.Y();
755 }
756 else {
757 isoValue = valueP1.Y();
758 isoPar1 = valueP1.X();
759 isoPar2 = valueP2.X();
760 }
761 // Fin calcul sens de courbe iso
762 }
763 } // end of fix check 05.02.02
764 }
765 }
766
767 // Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
768 // (utile pour detromper bug de projection). Mais detromper aussi circularite
769 //else {
770 //if (p1OnIso) valueP1 =
771 //BestExtremum (valueP1,points(1),points(2));
772 //if (p2OnIso) valueP2 =
773 //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
774 //}
775
776 Standard_Real gap = myPreci; //:q1
777 Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
778 if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
779 //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
780 if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
781 ChangeCycle = Standard_True;
782 //for( i = 1; i <= nbrPnt; i ++) {
783 for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
784 if(ChangeCycle) //skl for OCC3430
785 i=nbrPnt-ii+1;
786 else
787 i=ii;
788 p3d = points(i);
789 if (isoParam) {
790
791 if (isoPar2d3d) {
792 if (isoPar2 > isoPar1) tPar = params(i);
793 else tPar = t1 + t2 - params(i);
794 } else if (!isAnalytic) {
795 // projection to iso
796 if (i==1) tPar = isoPar1;
797 else if (i==nbrPnt) tPar = isoPar2;
798 else {
799 tPar = pout(i);
800 //:S4030 ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
801 }
802 }
803
804 if (!isoPar2d3d && isAnalytic) {
805 if (i == 1) p2d = valueP1;
806 else if (i == nbrPnt) p2d = valueP2;
807 else {
808 p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
809 Precision::Confusion()+1000*gap); //:q1
810 gap = mySurf->Gap();
811 }
812 } else {
813 if(isoTypeU) { p2d.SetX(isoValue); p2d.SetY(tPar); }
814 else { p2d.SetX(tPar); p2d.SetY(isoValue); }
815 }
816 }
817
818 else {
819 if ( (i == 1) && p1OnIso) p2d = valueP1;
820 else if( (i == nbrPnt) && p2OnIso) p2d = valueP2;
821 else {// general case (not an iso) mais attention aux singularites !
15b54261 822 // first and last points are already computed by getLine()
823 if ( (i == 1 || i == nbrPnt))
824 {
825 if (!isRecompute)
826 {
827 p2d = pnt2d(i);
828 gap = mySurf->Gap();
829 continue;
830 }
831 else
832 {
833 //:q9 abv 23 Mar 99: use cashe as 1st approach
834 Standard_Integer j; // svv #1
835 for ( j=0; j < myNbCashe; j++ )
836 if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
837 {
838 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci,
839 Precision::Confusion()+gap);
840 break;
841 }
842 if ( j >= myNbCashe ) p2d = mySurf->ValueOfUV(p3d, myPreci);
843 }
844 }
7fd59977 845 else {
846 p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
847 Precision::Confusion()+1000*gap); //:q1
848 }
849 gap = mySurf->Gap();
850 }
851 }
852 pnt2d (i) = p2d;
853 if ( ii > 1 ) {
854 if(ChangeCycle)
855 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i+1).XY() );
856 else
857 p2d.SetXY ( 2. * p2d.XY() - pnt2d(i-1).XY() );
858 }
859 }
860
861 //pdn %12 11.02.99 PRO9234 entity 15402
862 if (!isoPar2d3d) {
863 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
864 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
865 }
866
867 // attention aux singularites ... (hors cas iso qui les traite deja)
868 // if (!isoParam) {
869 // p2d = pnt2d (1);
870 // if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
871 // pnt2d (1) = p2d;
872 // p2d = pnt2d (nbrPnt);
873 // if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
874 // pnt2d (nbrPnt) = p2d;
875 // }
876
877 // Si la surface est UCLosed et VClosed, on recadre les points
878 // algo un peu complique, on retarde l implementation
879 Standard_Real Up = ul - uf;
880 Standard_Real Vp = vl - vf;
881 Standard_Real dist2d;
0797d9d3 882#ifdef OCCT_DEBUG
7fd59977 883 if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
884 cout << "WARNING : Recadrage incertain sur U & VClosed" << endl;
885 }
886#endif
887 // Si la surface est UCLosed, on recadre les points
888 if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
889 // Premier point dans le domain [uf, ul]
890 Standard_Real prevX, firstX = pnt2d (1).X();
891 while (firstX < uf) { firstX += Up; pnt2d (1).SetX(firstX); }
892 while (firstX > ul) { firstX -= Up; pnt2d (1).SetX(firstX); }
893 prevX = firstX;
894
895 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
896 Standard_Real minX = firstX, maxX = firstX;
897
898 // On decalle toujours le suivant
899 for (i = 2; i <= nbrPnt; i++) {
900 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
901 Standard_Real CurX = pnt2d (i).X();
902 dist2d = Abs (CurX - prevX);
903 if (dist2d > ( Up / 2) ) {
904 if (CurX > prevX + Up/2) {
905 while (CurX > prevX + Up/2) { CurX -= Up; pnt2d (i).SetX (CurX); }
906 } else if (CurX < prevX - Up/2) {
907 while (CurX < prevX - Up/2) { CurX += Up; pnt2d (i).SetX (CurX); }
908 }
909
910 }
911 prevX = CurX;
912 if ( minX > CurX ) minX = CurX; //:97
913 else if ( maxX < CurX ) maxX = CurX; //:97
914 }
915
916 //:97
917 Standard_Real midX = 0.5 * ( minX + maxX );
918 Standard_Real shiftX=0.;
919 if ( midX > ul ) shiftX = -Up;
920 else if ( midX < uf ) shiftX = Up;
921 if ( shiftX != 0. )
922 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetX ( pnt2d(i).X() + shiftX );
923 }
924 // Si la surface est VCLosed, on recadre les points
925 // Same code as UClosed : optimisation souhaitable !!
926 // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
927 // Si la surface est UCLosed, on recadre les points
928 //
929 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
930 //#78 rln 12.03.99 S4135
931 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
932 // Premier point dans le domain [vf, vl]
933 Standard_Real prevY, firstY = pnt2d (1).Y();
934 while (firstY < vf) { firstY += Vp; pnt2d (1).SetY(firstY); }
935 while (firstY > vl) { firstY -= Vp; pnt2d (1).SetY(firstY); }
936 prevY = firstY;
937
938 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
939 Standard_Real minY = firstY, maxY = firstY;
940
941 // On decalle toujours le suivant
942 for (i = 2; i <= nbrPnt; i ++) {
943 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
944 Standard_Real CurY = pnt2d (i).Y();
945 dist2d = Abs (CurY - prevY);
946 if (dist2d > ( Vp / 2) ) {
947 if (CurY > prevY + Vp/2) {
948 while (CurY > prevY + Vp/2) { CurY -= Vp; pnt2d (i).SetY (CurY); }
949 } else if (CurY < prevY - Vp/2) {
950 while (CurY < prevY - Vp/2) { CurY += Vp; pnt2d (i).SetY (CurY); }
951 }
952 }
953 prevY = CurY;
954 if ( minY > CurY ) minY = CurY; //:97
955 else if ( maxY < CurY ) maxY = CurY; //:97
956 }
957
958 //:97
959 Standard_Real midY = 0.5 * ( minY + maxY );
960 Standard_Real shiftY=0.;
961 if ( midY > vl ) shiftY = -Vp;
962 else if ( midY < vf ) shiftY = Vp;
963 if ( shiftY != 0. )
964 for ( i=1; i <= nbrPnt; i++ ) pnt2d(i).SetY ( pnt2d(i).Y() + shiftY );
965 }
966
967 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
968 //#78 rln 12.03.99 S4135
969 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
970 for (i = 2; i <= nbrPnt; i++) {
971 //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
972 dist2d = Abs (pnt2d(i).Y() - pnt2d(i - 1).Y());
973 if (dist2d > ( Vp / 2) ) {
974 // ATTENTION : il faut regarder ou le decalage se fait.
975 // si plusieurs points sont decalles, il faut plusieurs passes
976 // pour obtenir un resultat correct.
977 // NOT YET IMPLEMENTED
978
979 // one of those point is incorrectly placed
980 // i.e on the wrong side of the "seam"
981 // on prend le point le plus pres des bords vf ou vl
982 Standard_Boolean prevOnFirst = Standard_False;
983 Standard_Boolean prevOnLast = Standard_False;
984 Standard_Boolean currOnFirst = Standard_False;
985 Standard_Boolean currOnLast = Standard_False;
986
987 // .X ? plutot .Y , non ?
988 Standard_Real distPrevVF = Abs(pnt2d (i-1).Y() - vf);
989 Standard_Real distPrevVL = Abs(pnt2d (i-1).Y() - vl);
990 Standard_Real distCurrVF = Abs(pnt2d (i).Y() - vf);
991 Standard_Real distCurrVL = Abs(pnt2d (i).Y() - vl);
992
993 Standard_Real theMin = distPrevVF;
994 prevOnFirst = Standard_True;
995 if (distPrevVL < theMin) {
996 theMin = distPrevVL;
997 prevOnFirst = Standard_False;
998 prevOnLast = Standard_True;
999 }
1000 if (distCurrVF < theMin) {
1001 theMin = distCurrVF;
1002 prevOnFirst = Standard_False;
1003 prevOnLast = Standard_False;
1004 currOnFirst = Standard_True;
1005 }
1006 if (distCurrVL < theMin) {
1007 theMin = distCurrVL;
1008 prevOnFirst = Standard_False;
1009 prevOnLast = Standard_False;
1010 currOnFirst = Standard_False;
1011 currOnLast = Standard_True;
1012 }
1013 // Modifs RLN/Nijni 3-DEC-1997
1014 if (prevOnFirst) {
1015 // on decalle le point (i-1) en V Last
1016 gp_Pnt2d newPrev(pnt2d (i-1).X(), vf); // instead of vl RLN/Nijni
1017 pnt2d (i-1) = newPrev;
1018 }
1019 else if (prevOnLast) {
1020 // on decalle le point (i-1) en V first
1021 gp_Pnt2d newPrev(pnt2d (i-1).X(), vl); // instead of vf RLN/Nijni
1022 pnt2d (i-1) = newPrev;
1023 }
1024 else if (currOnFirst) {
1025 // on decalle le point (i) en V Last
1026 gp_Pnt2d newCurr(pnt2d (i).X(),vf); // instead of vl RLN/Nijni
1027 pnt2d (i) = newCurr;
1028 }
1029 else if (currOnLast) {
1030 // on decalle le point (i) en V First
1031 gp_Pnt2d newCurr(pnt2d (i).X(), vl); // instead of vf RLN/Nijni
1032 pnt2d (i) = newCurr;
1033 }
1034 // on verifie
0797d9d3 1035#ifdef OCCT_DEBUG
7fd59977 1036 dist2d = pnt2d (i-1).Distance(pnt2d (i));
1037 if (dist2d > ( Vp / 2) ) {
1038 cout << "Echec dans le recadrage" << endl;
1039 }
1040#endif
1041 }
1042 }
1043 }
1044
1045 //:c0 abv 20 Feb 98: treat very special case when 3d curve
1046 // go over the pole of, e.g., sphere, and partly lies along seam.
1047 // 2d representation of such a curve should consist of 3 parts - one on
1048 // regular part of surface (interior), one part along degenerated boundary
1049 // and one along seam.
1050 // Since it cannot be adjusted later by arranging pcurves (curve is single),
1051 // to fix it it is nesessary to have a possibility of adjusting seam
1052 // part of such curve either to left or right boundary of surface.
1053 // Test is performed only if flag AdjustOverDegen is not -1.
1054 // If AdjustOverDegen is True, seam part of curve is adjusted to
1055 // the left, and if False - to the right parametric boundary
1056 // If treated case is detected, flag DONE4 is set to status
1057 // NOTE: currently, precision is Precision::PConfusion() since it
1058 // is enough on encountered example
1059 // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1060 // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1061 if ( myAdjustOverDegen != -1 ) {
1062 if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1063 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1064 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1065 // 1st, find gap point (degenerated pole)
1066 Standard_Real PrevX=0.;
1067 Standard_Integer OnBound=0, PrevOnBound=0;
1068 Standard_Integer ind; // svv #1
1069 Standard_Boolean start = Standard_True;
1070 for ( ind=1; ind <= nbrPnt; ind++ ) {
1071 Standard_Real CurX = pnt2d(ind).X();
1072 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1073 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1074 continue;
1075 OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1076 Precision::PConfusion() );
1077 if ( ! start && Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1078 break;
1079 start = Standard_False;
1080 PrevX = CurX;
1081 PrevOnBound = OnBound;
1082 }
1083 // if found, adjust seam part
1084 if ( ind <= nbrPnt ) {
1085 PrevX = ( myAdjustOverDegen ? uf : ul );
1086 Standard_Real dU = Up/2 + Precision::PConfusion();
1087 if ( PrevOnBound ) {
1088 pnt2d(ind-1).SetX ( PrevX );
1089 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1090 Standard_Real CurX = pnt2d(j).X();
1091 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1092 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1093 }
1094 }
1095 else if ( OnBound ) {
1096 pnt2d(ind).SetX ( PrevX );
1097 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1098 Standard_Real CurX = pnt2d(j).X();
1099 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1100 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1101 }
1102 }
1103 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1104 }
1105 }
1106 }
1107 else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1108 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1109 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1110 // 1st, find gap point (degenerated pole)
1111 Standard_Real PrevY=0.;
1112 Standard_Integer OnBound=0, PrevOnBound=0;
1113 Standard_Integer ind; // svv #1
1114 Standard_Boolean start = Standard_True;
1115 for ( ind=1; ind <= nbrPnt; ind++ ) {
1116 Standard_Real CurY = pnt2d(ind).Y();
1117 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1118 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1119 continue;
1120 OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1121 Precision::PConfusion() );
1122 if ( ! start && Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1123 break;
1124 start = Standard_False;
1125 PrevY = CurY;
1126 PrevOnBound = OnBound;
1127 }
1128 // if found, adjust seam part
1129 if ( ind <= nbrPnt ) {
1130 PrevY = ( myAdjustOverDegen ? vf : vl );
1131 Standard_Real dV = Vp/2 + Precision::PConfusion();
1132 if ( PrevOnBound ) {
1133 pnt2d(ind-1).SetY ( PrevY );
1134 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1135 Standard_Real CurY = pnt2d(j).Y();
1136 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1137 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1138 }
1139 }
1140 else if ( OnBound ) {
1141 pnt2d(ind).SetY ( PrevY );
1142 for ( Standard_Integer j=ind+1; j <= nbrPnt; j++ ) {
1143 Standard_Real CurY = pnt2d(j).Y();
1144 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1145 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1146 }
1147 }
1148 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1149 }
1150 }
1151 }
1152 }
1153
1154 //:q9: fill cashe
1155 myNbCashe = 2;
1156 if(ChangeCycle) { // msv 10.08.04: avoid using of uninitialised field
1157 //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1158 // myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1159 myCashe3d[0] = points(1);
1160 myCashe3d[1] = points(nbrPnt);
1161 myCashe2d[0] = pnt2d(1);
1162 myCashe2d[1] = pnt2d(nbrPnt);
1163 }
1164 else {
1165 myCashe3d[1] = points(1);
1166 myCashe3d[0] = points(nbrPnt);
1167 myCashe2d[1] = pnt2d(1);
1168 myCashe2d[0] = pnt2d(nbrPnt);
1169 }
7fd59977 1170 return isDone;
1171}
1172
1173//=======================================================================
1174//function : ApproximatePCurve
1175//purpose :
1176//=======================================================================
1177
1178Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1179 Handle(TColgp_HArray1OfPnt2d)& points2d,
1180 Handle(TColStd_HArray1OfReal)& params,
1181 const Handle(Geom_Curve)& /*orig*/) const
1182{
1183// Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1184 Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1185 Handle(Geom2d_Curve) C2d;
1186 try {
1187 OCC_CATCH_SIGNALS
1188 CheckPoints2d (points2d, params, theTolerance2d);
1189 Standard_Integer numberPnt = points2d->Length();
1190
1191 TColgp_Array1OfPnt points3d(1,numberPnt);
1192 gp_Pnt2d pnt2d;
1193 gp_Pnt pnt;
1194 Standard_Integer i; // svv #1
1195 for( i = 1; i <= numberPnt; i++) {
1196 pnt2d = points2d->Value(i);
1197 pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1198 points3d(i) = pnt;
1199 }
1200
1201 GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1202 Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1203 Standard_Integer NbPoles = crv3d->NbPoles();
1204 TColgp_Array1OfPnt poles3d (1, NbPoles);
1205 TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1206 crv3d->Poles(poles3d);
1207 for( i = 1; i <= NbPoles; i++) {
1208 pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1209 poles2d(i) = pnt2d;
1210 }
1211 TColStd_Array1OfReal weights (1,NbPoles);
1212 TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1213 TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1214 crv3d->Knots(knots);
1215 crv3d->Weights(weights);
1216 crv3d->Multiplicities(multiplicities);
1217 C2d = new Geom2d_BSplineCurve ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1218 return C2d;
1219 }
1220 catch(Standard_Failure) {
0797d9d3 1221#ifdef OCCT_DEBUG //:s5
7fd59977 1222 // debug ...
1223 Standard_Integer nbp = params->Length();
1224 Standard_Integer nb2 = points2d->Length();
1225 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1226 Standard_Failure::Caught()->Print(cout);
1227 cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1228// if (nb2 > nbp) nb2 = nbp;
1229// Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1230// // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1231// dbl.AddReals (rb2,rbp);
1232// for (Standard_Integer i = 1; i <= nb2; i ++) {
1233// gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1234// dbl.AddXYZ (quoi);
1235// }
1236#endif
1237 C2d.Nullify();
1238 }
1239 return C2d;
1240}
1241
1242//=======================================================================
1243//function : InterpolatePCurve
1244//purpose :
1245//=======================================================================
1246
1247 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1248 Handle(TColgp_HArray1OfPnt2d)& points2d,
1249 Handle(TColStd_HArray1OfReal)& params,
1250 const Handle(Geom_Curve)& /*orig*/) const
1251{
1252 Handle(Geom2d_Curve) C2d; // NULL si echec
1253 Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1254 try {
1255 OCC_CATCH_SIGNALS
1256 // on verifie d abord s il n y a pas de points confondus
1257 // si besoin on retouche les valeurs ...
1258 CheckPoints2d (points2d, params, theTolerance2d);
1259 Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1260 Standard_False, theTolerance2d);
1261 myInterPol2d.Perform();
1262 if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1263 }
1264 catch(Standard_Failure) {
0797d9d3 1265#ifdef OCCT_DEBUG //:s5
7fd59977 1266// // debug ...
1267 Standard_Integer nbp = params->Length();
1268 Standard_Integer nb2 = points2d->Length();
1269 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1270 Standard_Failure::Caught()->Print(cout);
1271 cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<endl;
1272// if (nb2 > nbp) nb2 = nbp;
1273// Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1274// // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1275// dbl.AddReals (rb2,rbp);
1276// for (Standard_Integer i = 1; i <= nb2; i ++) {
1277// gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1278// dbl.AddXYZ (quoi);
1279// }
1280#endif
1281 C2d.Nullify();
1282 }
1283 return C2d;
1284}
1285
1286//=======================================================================
1287//function : InterpolateCurve3d
1288//purpose :
1289//=======================================================================
1290
1291Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1292 Handle(TColgp_HArray1OfPnt)& points,
1293 Handle(TColStd_HArray1OfReal)& params,
1294 const Handle(Geom_Curve)& /*orig*/) const
1295{
1296 Handle(Geom_Curve) C3d; // NULL si echec
1297 try {
1298 OCC_CATCH_SIGNALS
1299 Standard_Real Tol = myPreci;
1300 CheckPoints(points, params, Tol);
1301 GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1302 myInterPol.Perform();
1303 if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1304 }
1305 catch(Standard_Failure) {
1306 C3d.Nullify();
0797d9d3 1307#ifdef OCCT_DEBUG //:s5
7fd59977 1308 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1309 Standard_Failure::Caught()->Print(cout); cout << endl;
1310#endif
1311 }
1312 return C3d;
1313}
1314
1315//=======================================================================
1316//function : CheckPoints
1317//purpose :
1318//=======================================================================
1319
1320 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1321{
1322 Standard_Integer firstElem = points->Lower();
1323 Standard_Integer lastElem = points->Upper();
1324 Standard_Integer i;
1325 Standard_Integer nbPntDropped = 0;
1326 Standard_Integer lastValid = firstElem; // indice of last undropped point
1327
1328 // will store 0 when the point is to be removed, 1 otherwise
1329 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1330 for (i = firstElem; i<=lastElem ; i++) tmpParam.SetValue(i,1);
7ae65f0d 1331 Standard_Real DistMin2 = RealLast();
7fd59977 1332 gp_Pnt Prev = points->Value (lastValid);
1333 gp_Pnt Curr;
1334 for (i = firstElem + 1; i <= lastElem ; i ++) {
1335 Curr = points->Value(i);
7ae65f0d 1336 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1337 if (CurDist2 < gp::Resolution()) { // test 0
7fd59977 1338 nbPntDropped ++;
1339 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1340 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1341 } else {
7ae65f0d 1342 if (CurDist2 < DistMin2)
1343 DistMin2 = CurDist2;
7fd59977 1344 // lastValid becomes the current (i.e. i)
1345 lastValid = i;
1346 Prev = Curr;
1347 }
1348 }
7ae65f0d 1349 if (DistMin2 < RealLast())
1350 preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1351 if (nbPntDropped == 0)
7fd59977 1352 return;
7ae65f0d 1353
0797d9d3 1354#ifdef OCCT_DEBUG
7fd59977 1355 cout << "Warning : removing 3d points for interpolation" << endl;
1356#endif
1357 // Build new HArrays
1358 Standard_Integer newLast = lastElem - nbPntDropped;
1359 if ((newLast - firstElem + 1) < 2) {
0797d9d3 1360#ifdef OCCT_DEBUG
7fd59977 1361 cout << "Too many degenerated points for 3D interpolation" << endl;
1362#endif
1363 return;
1364 }
1365 Handle(TColgp_HArray1OfPnt) newPnts =
1366 new TColgp_HArray1OfPnt(firstElem, newLast);
1367 Handle(TColStd_HArray1OfReal) newParams =
1368 new TColStd_HArray1OfReal(firstElem, newLast);
1369 Standard_Integer newCurr = 1;
1370 for (i = firstElem; i<= lastElem ; i++) {
1371 if (tmpParam.Value(i) == 1) {
1372 newPnts->SetValue(newCurr, points->Value(i));
1373 newParams->SetValue(newCurr, params->Value(i));
1374 newCurr ++;
1375 }
1376 }
1377 points = newPnts;
1378 params = newParams;
7fd59977 1379 // on la reduit un peu
1380}
1381
1382//=======================================================================
1383//function : CheckPoints2d
1384//purpose :
1385//=======================================================================
1386
1387 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1388 Handle(TColStd_HArray1OfReal)& params,
1389 Standard_Real& preci) const
1390{
1391 Standard_Integer firstElem = points->Lower();
1392 Standard_Integer lastElem = points->Upper();
1393 Standard_Integer i;
1394 Standard_Integer nbPntDropped = 0;
1395 Standard_Integer lastValid = firstElem; // indice of last undropped point
1396
1397 // will store 0 when the point is to be removed, 1 otherwise
1398 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1399 for (i = firstElem; i<=lastElem ; i++) {
1400 tmpParam.SetValue(i,1);
1401 }
7ae65f0d 1402 Standard_Real DistMin2 = RealLast();
7fd59977 1403 gp_Pnt2d Prev = points->Value(lastValid);
1404 gp_Pnt2d Curr;
1405 for (i = firstElem + 1; i<=lastElem ; i++) {
1406 Curr = points->Value(i);
7ae65f0d 1407 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1408 if (CurDist2 < gp::Resolution()) { // test 0
7fd59977 1409 nbPntDropped ++;
1410 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1411 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1412 } else {
7ae65f0d 1413 if (CurDist2 < DistMin2)
1414 DistMin2 = CurDist2;
7fd59977 1415 // lastValid becomes the current (i.e. i)
1416 lastValid = i;
1417 Prev = Curr;
1418 }
1419 }
7ae65f0d 1420 if (DistMin2 < RealLast())
1421 preci = 0.9 * Sqrt (DistMin2);
1422 if (nbPntDropped == 0)
7fd59977 1423 return;
7ae65f0d 1424
0797d9d3 1425#ifdef OCCT_DEBUG
7fd59977 1426 cout << "Warning : removing 2d points for interpolation" << endl;
1427#endif
1428 // Build new HArrays
1429 Standard_Integer newLast = lastElem - nbPntDropped;
1430 if ((newLast - firstElem + 1) < 2) {
0797d9d3 1431#ifdef OCCT_DEBUG
7fd59977 1432 cout << "Too many degenerated points for 2D interpolation" << endl;
1433#endif
1434 //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1435 tmpParam.SetValue(firstElem,1);
1436 tmpParam.SetValue(lastElem,1);
1437 gp_XY lastPnt = points->Value(lastElem).XY();
1438 lastPnt.Add(gp_XY(preci,preci));
1439 points->SetValue(lastElem,lastPnt);
1440 newLast = firstElem+1;
1441 //return;
1442 }
1443 Handle(TColgp_HArray1OfPnt2d) newPnts =
1444 new TColgp_HArray1OfPnt2d(firstElem, newLast);
1445 Handle(TColStd_HArray1OfReal) newParams =
1446 new TColStd_HArray1OfReal(firstElem, newLast);
1447 Standard_Integer newCurr = 1;
1448 for (i = firstElem; i <= lastElem ; i++) {
1449 if (tmpParam.Value(i) == 1) {
0797d9d3 1450#ifdef OCCT_DEBUG
7fd59977 1451 cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1452#endif
1453 newPnts->SetValue(newCurr, points->Value(i));
1454 newParams->SetValue(newCurr, params->Value(i));
1455 newCurr ++;
1456 }
1457 else {
0797d9d3 1458#ifdef OCCT_DEBUG
7fd59977 1459 cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << endl;
1460#endif
1461 }
1462 }
1463 points = newPnts;
1464 params = newParams;
7fd59977 1465}
1466
1467//=======================================================================
1468//function : IsAnIsoparametric
1469//purpose :
1470//=======================================================================
1471//:S4030: modified for optimization
1472//:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1473
1474 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
1475 const TColgp_Array1OfPnt& points,
1476 const TColStd_Array1OfReal& params,
1477 Standard_Boolean& isoTypeU,
1478 Standard_Boolean& p1OnIso,
1479 gp_Pnt2d& valueP1,
1480 Standard_Boolean& p2OnIso,
1481 gp_Pnt2d& valueP2,
1482 Standard_Boolean& isoPar2d3d,
1483 Handle(Geom_Curve)& cIso,
1484 Standard_Real& t1,
1485 Standard_Real& t2,
1486 TColStd_Array1OfReal& pout) const
1487{
1488 try { // RAJOUT
1489 OCC_CATCH_SIGNALS
1490
1491 Standard_Real prec = Precision::Confusion();//myPreci;
1492
1493 Standard_Boolean isoParam = Standard_False;
1494 isoPar2d3d = Standard_False;
1495
1496 Standard_Real U1, U2, V1, V2;
1497 mySurf->Bounds(U1, U2, V1, V2);
1498
1499 if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1500 Handle(Geom_RectangularTrimmedSurface) sTrim =
1501 Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
1502 sTrim->Bounds(U1, U2, V1, V2);
1503 }
1504
1505 gp_Pnt pt;
1506 Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
d20d815b 1507 Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
7fd59977 1508 Standard_Real mindist2;
1509 Standard_Real mind2[2];
1510 mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
1511
1512 p1OnIso = Standard_False;
1513 p2OnIso = Standard_False;
1514 const Bnd_Box* aBox = 0;
1515
1516 for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
1517 Standard_Real isoVal=0.;
1518 Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
1519 Handle(Geom_Curve) cI;
1520 Standard_Real tt1, tt2;
1521
1522 if (j == 1 ) {
1523 if (Precision::IsInfinite(U1)) continue;
1524 cI = mySurf->UIso(U1);
1525 isoU = Standard_True;
1526 isoVal = U1;
1527 aBox = & mySurf->GetBoxUF();
1528 }
1529 else if (j == 2) {
1530 if (Precision::IsInfinite(U2)) continue;
1531 cI = mySurf->UIso(U2);
1532 isoU = Standard_True;
1533 isoVal = U2;
1534 aBox = & mySurf->GetBoxUL();
1535 }
1536 else if (j == 3) {
1537 if (Precision::IsInfinite(V1)) continue;
1538 cI = mySurf->VIso(V1);
1539 isoU = Standard_False;
1540 isoVal = V1;
1541 aBox = & mySurf->GetBoxVF();
1542 }
1543 else if (j == 4) {
1544 if (Precision::IsInfinite(V2)) continue;
1545 cI = mySurf->VIso(V2);
1546 isoU = Standard_False;
1547 isoVal = V2;
1548 aBox = & mySurf->GetBoxVL();
1549 }
1550 if(cI.IsNull())
1551 continue;
1552
1553 if (isoU) { tt1 = V1; tt2 = V2; }
1554 else { tt1 = U1; tt2 = U2; }
1555
1556 gp_Pnt ext1, ext2;
1557 cI->D0(tt1, ext1);
1558 cI->D0(tt2, ext2);
1559
1560// PATCH CKY 9-JUL-1998 : protection contre singularite
1561 gp_Pnt extmi;
1562 cI->D0( (tt1+tt2)/2,extmi);
1563 if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
1564
1565 Standard_Boolean PtEQext1 = Standard_False;
1566 Standard_Boolean PtEQext2 = Standard_False;
1567
eb1ebea4 1568 Standard_Real currd2[2], tp[2] = {0, 0};
7fd59977 1569 Standard_Integer mp[2];
1570
1571 for (Standard_Integer i=0; i<2; i++) {
1572 mp[i] = 0;
1573 Standard_Integer k = (i == 0 ? 1 : nbrPnt);
1574
1575 // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
1576 currd2[i] = points(k).SquareDistance ( ext1 );
1577 if ( currd2[i] <= prec*prec && !PtEQext1) {
1578 mp[i] = 1;
1579 tp[i] = tt1;
1580 PtEQext1 = Standard_True;
1581 continue;
1582 }
1583
1584 currd2[i] = points(k).SquareDistance ( ext2 );
1585 if ( currd2[i] <= prec*prec && !PtEQext2) {
1586 mp[i] = 2;
1587 tp[i] = tt2;
1588 PtEQext2 = Standard_True;
1589 continue;
1590 }
1591
1592 // On evite de projecter sur un iso degenere
1593 // on doit egalement le faire pour l apex du cone
1594 if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
1595 continue;
1596 }
1597
1598 if(aBox->IsOut(points(k))) continue;
1599
1600 Standard_Real Cf = cI->FirstParameter();
1601 Standard_Real Cl = cI->LastParameter();
1602 if (Precision::IsInfinite(Cf)) Cf = -1000;
1603 if (Precision::IsInfinite(Cl)) Cl = +1000;
1604
1605 ShapeAnalysis_Curve sac;
1606 Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
1607 currd2[i] = dist * dist;
1608 if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
1609 mp[i] = 3;
1610 tp[i] = t;
1611 }
1612 }
1613
1614 //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
1615 // avoid possible null-length curves
1616 if ( mp[0] >0 && mp[1] >0 &&
1617 Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
1618
1619
1620 if (mp[0] > 0 &&
1621 ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
1622 p1OnIso = Standard_True;
15b54261 1623 mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
7fd59977 1624 if (isoU) valueP1.SetCoord(isoVal, tp[0]);
1625 else valueP1.SetCoord(tp[0], isoVal);
1626 }
1627
1628 if (mp[1] > 0 &&
1629 ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
1630 p2OnIso = Standard_True;
1631 mind2[1] = currd2[1];
1632 if (isoU) valueP2.SetCoord(isoVal, tp[1]);
1633 else valueP2.SetCoord(tp[1], isoVal);
1634 }
1635
1636 if ( mp[0] <=0 || mp[1] <=0 ) continue;
1637
1638 Standard_Real md2 = currd2[0] + currd2[1];
1639 if ( mindist2 <= md2 ) continue;
1640
1641 mindist2 = md2;
1642 mpt[0] = mp[0];
1643 mpt[1] = mp[1];
1644 tpar[0] = tp[0];
1645 tpar[1] = tp[1];
1646 isoTypeU = isoU;
1647 isoValue = isoVal;
1648 cIso = cI;
1649 t1 = tt1;
1650 t2 = tt2;
1651 }
1652
1653 // probablely it concerns an isoparametrics
1654 if ( mpt[0] >0 && mpt[1] >0 ) {
1655
1656 p1OnIso = p2OnIso = Standard_True;
1657 if (isoTypeU) {
1658 valueP1.SetCoord(isoValue, tpar[0]);
1659 valueP2.SetCoord(isoValue, tpar[1]);
1660 }
1661 else {
1662 valueP1.SetCoord(tpar[0], isoValue);
1663 valueP2.SetCoord(tpar[1], isoValue);
1664 }
1665
1666 if ( mpt[0] != 3 && mpt[1] != 3 ) {
1667 isoPar2d3d = Standard_True;
1668 for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
1669 if (tpar[1] > tpar[0]) t = params(i);
1670 else t = t1+t2-params(i);
1671 cIso->D0(t, pt);
1672 if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
1673 }
1674 }
1675
1676 if (isoPar2d3d) isoParam = Standard_True;
1677 else {
1678 Standard_Real prevParam = tpar[0];
1679 Standard_Real Cf, Cl;
1680 Standard_Boolean isoByDistance = Standard_True;
1681 Cf = cIso->FirstParameter();
1682 Cl = cIso->LastParameter();
1683 if (Precision::IsInfinite(Cf)) Cf = -1000;
1684 if (Precision::IsInfinite(Cl)) Cl = +1000;
1685
1686 ShapeAnalysis_Curve sac;
1687 for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
1688 Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
1689 prec,pt,t,Cf,Cl,
1690 Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
1691 prevParam = t;
1692 pout(i)=t;
1693 if( (dist > prec) || (t < Cf) || (t > Cl) )
1694 isoByDistance = Standard_False;
1695 }
1696 if (isoByDistance) isoParam = Standard_True;
1697 }
1698 }
1699/* if (!isoParam) { CKY 29-mai-1997 : garder tout ce qu on peut ?
1700 p1OnIso = Standard_False;
1701 p2OnIso = Standard_False;
1702 } */
1703 return isoParam;
1704 } // RAJOUT
1705 catch(Standard_Failure) {
1706// pb : on affiche ce qu on peut
0797d9d3 1707#ifdef OCCT_DEBUG
7fd59977 1708 for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
1709 cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
1710 points(numpnt).X()<<" "<<points(numpnt).Y()<<" "<<points(numpnt).Z()<<")"<<endl;
1711 }
1712#endif
0797d9d3 1713#ifdef OCCT_DEBUG //:s5
7fd59977 1714 cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
1715 Standard_Failure::Caught()->Print(cout); cout << endl;
1716#endif
1717 return Standard_False;
1718 }
7fd59977 1719}
1720
1721/* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
1722//=======================================================================
1723//function : BestExtremum
1724//purpose : auxiliaire prenant le meilleur extremum si ISO car doute possible
1725//=======================================================================
1726
1727 gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
1728{
1729// P2iso a ete calcule depuis P3ext sur une iso externe de la surface
1730// En principe bon mais circularite possible ... et IsU/VClosed faillible
1731// (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
1732// 1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
1733// 2/ si iso, attention a la circularite (cas limite)
1734
1735// NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
1736
1737// D abord, calcul p2ext depuis la surface. choix surface/iso
1738 return P2iso;
1739 Standard_Real prec = Precision::Confusion();//myPreci;
1740 gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
1741 gp_Pnt P3cal = mySurf->Value (P2cal);
1742 Standard_Real dcal = P3ext.Distance (P3cal);
1743 Standard_Real dnxt = P3ext.Distance (P3next);
1744 if (dcal > dnxt) return P2iso; // en fait protection sur BUG (PRO8468)
1745
1746// On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
1747 gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
1748 if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;
1749 return P2iso;
1750}
1751*/