projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031687: Draw Harness, ViewerTest - extend command vrenderparams with option updating...
[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
7fd59977 31
42cf5bc1 32#include <Approx_CurveOnSurface.hxx>
33#include <Geom2d_BSplineCurve.hxx>
9a6ea9c4 34#include <Geom2d_Circle.hxx>
42cf5bc1 35#include <Geom2d_Curve.hxx>
9a6ea9c4 36#include <Geom2d_Ellipse.hxx>
37#include <Geom2d_Hyperbola.hxx>
42cf5bc1 38#include <Geom2d_Line.hxx>
9a6ea9c4 39#include <Geom2d_Parabola.hxx>
7fd59977 40#include <Geom2d_TrimmedCurve.hxx>
42cf5bc1 41#include <Geom2dAdaptor.hxx>
42#include <Geom2dAPI_Interpolate.hxx>
43#include <Geom_BezierSurface.hxx>
44#include <Geom_BoundedCurve.hxx>
7fd59977 45#include <Geom_BSplineCurve.hxx>
42cf5bc1 46#include <Geom_BSplineSurface.hxx>
47#include <Geom_Curve.hxx>
7fd59977 48#include <Geom_OffsetSurface.hxx>
49#include <Geom_Plane.hxx>
42cf5bc1 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>
7fd59977 55#include <GeomAdaptor_HCurve.hxx>
42cf5bc1 56#include <GeomAdaptor_HSurface.hxx>
57#include <GeomAPI_Interpolate.hxx>
58#include <GeomAPI_PointsToBSpline.hxx>
59#include <GeomProjLib.hxx>
60#include <gp_Pnt2d.hxx>
2c26a53d 61#include <ElCLib.hxx>
42cf5bc1 62#include <NCollection_Sequence.hxx>
63#include <Precision.hxx>
64#include <ProjLib_CompProjectedCurve.hxx>
65#include <ProjLib_HCompProjectedCurve.hxx>
66#include <ProjLib_ProjectedCurve.hxx>
02fd709b 67#include <ShapeAnalysis.hxx>
7fd59977 68#include <ShapeAnalysis_Curve.hxx>
69#include <ShapeAnalysis_Surface.hxx>
42cf5bc1 70#include <ShapeConstruct_ProjectCurveOnSurface.hxx>
7fd59977 71#include <ShapeExtend.hxx>
42cf5bc1 72#include <Standard_ErrorHandler.hxx>
73#include <Standard_Failure.hxx>
74#include <Standard_Type.hxx>
75#include <TColgp_Array1OfPnt.hxx>
76#include <TColStd_Array1OfInteger.hxx>
2c26a53d 77#include <IntRes2d_Domain.hxx>
78#include <IntCurve_IntConicConic.hxx>
7fd59977 79
42cf5bc1 80#include <algorithm>
25e59720 81IMPLEMENT_STANDARD_RTTIEXT(ShapeConstruct_ProjectCurveOnSurface,Standard_Transient)
92efcf78 82
7fd59977 83#define NCONTROL 23
b311480e 84
2c26a53d 85
86static void AdjustSecondPointToFirstPoint(const gp_Pnt2d& theFirstPoint,
87 gp_Pnt2d& theSecondPoint,
88 const Handle(Geom_Surface)& theSurf)
89{
90 if (theSurf->IsUPeriodic())
91 {
92 Standard_Real UPeriod = theSurf->UPeriod();
93 Standard_Real NewU = ElCLib::InPeriod(theSecondPoint.X(),
94 theFirstPoint.X() - UPeriod/2,
95 theFirstPoint.X() + UPeriod/2);
96 theSecondPoint.SetX(NewU);
97 }
98 if (theSurf->IsVPeriodic())
99 {
100 Standard_Real VPeriod = theSurf->VPeriod();
101 Standard_Real NewV = ElCLib::InPeriod(theSecondPoint.Y(),
102 theFirstPoint.Y() - VPeriod/2,
103 theFirstPoint.Y() + VPeriod/2);
104 theSecondPoint.SetY(NewV);
105 }
106}
107
108
7fd59977 109//=======================================================================
110//function : ShapeConstruct_ProjectCurveOnSurface
111//purpose :
112//=======================================================================
113
114ShapeConstruct_ProjectCurveOnSurface::ShapeConstruct_ProjectCurveOnSurface()
115{
116 myPreci = Precision::Confusion();
117 myBuild = Standard_False;
118 myAdjustOverDegen = 1; //:c0 //szv#4:S4163:12Mar99 was boolean
119 myNbCashe = 0; //:q9
120}
121
122//=======================================================================
123//function : Init
124//purpose :
125//=======================================================================
126
127 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(Geom_Surface)& surf,const Standard_Real preci)
128{
129 Init (new ShapeAnalysis_Surface (surf), preci);
130}
131
132//=======================================================================
133//function : Init
134//purpose :
135//=======================================================================
136
137 void ShapeConstruct_ProjectCurveOnSurface::Init(const Handle(ShapeAnalysis_Surface)& surf,const Standard_Real preci)
138{
139 SetSurface (surf);
140 SetPrecision (preci);
141}
142
143//=======================================================================
144//function : SetSurface
145//purpose :
146//=======================================================================
147
148 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(Geom_Surface)& surf)
149{
150 SetSurface (new ShapeAnalysis_Surface (surf));
151}
152
153//=======================================================================
154//function : SetSurface
155//purpose :
156//=======================================================================
157
158 void ShapeConstruct_ProjectCurveOnSurface::SetSurface(const Handle(ShapeAnalysis_Surface)& surf)
159{
160 if ( mySurf == surf ) return;
161 mySurf = surf;
162 myNbCashe = 0; //:q9
163}
164
165//=======================================================================
166//function : SetPrecision
167//purpose :
168//=======================================================================
169
170 void ShapeConstruct_ProjectCurveOnSurface::SetPrecision(const Standard_Real preci)
171{
172 myPreci = preci;
173}
174
175//=======================================================================
176//function : BuildCurveMode
177//purpose :
178//=======================================================================
179
180 Standard_Boolean& ShapeConstruct_ProjectCurveOnSurface::BuildCurveMode()
181{
182 return myBuild;
183}
184
185//=======================================================================
186//function : AdjustOverDegenMode
187//purpose :
188//=======================================================================
189//:c0
190
191//szv#4:S4163:12Mar99 was Boolean
192 Standard_Integer& ShapeConstruct_ProjectCurveOnSurface::AdjustOverDegenMode()
193 {
194 return myAdjustOverDegen;
195 }
196
7fd59977 197//=======================================================================
198//function : Status
199//purpose :
200//=======================================================================
201
9fd2d2c3 202 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Status (const ShapeExtend_Status theStatus) const
7fd59977 203{
9fd2d2c3 204 return ShapeExtend::DecodeStatus (myStatus, theStatus);
7fd59977 205}
206
207//=======================================================================
208//function : Perform
209//purpose :
210//=======================================================================
7fd59977 211Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::Perform (Handle(Geom_Curve)& c3d,
212 const Standard_Real First,
213 const Standard_Real Last,
214 Handle(Geom2d_Curve)& c2d,
2c26a53d 215 const Standard_Real TolFirst,
216 const Standard_Real TolLast)
7fd59977 217{
218 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_OK);
219 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 not needed
220
221 if (mySurf.IsNull()) {
222 c2d.Nullify();
223 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
224 return Standard_False;
225 }
7fd59977 226// Projection Analytique
227 Handle(Geom_Curve) crv3dtrim = c3d;
228 if ( ! c3d->IsKind(STANDARD_TYPE(Geom_BoundedCurve)) )
229 crv3dtrim = new Geom_TrimmedCurve ( c3d, First, Last );
230 c2d = ProjectAnalytic ( crv3dtrim );
231 if (!c2d.IsNull()) {
232 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
233 return Standard_True;
234 }
235
236// Projection par approximation
237
238 // discretize the 3d curve
239
240 Standard_Integer nbrPnt;
241
242// $$$$ :92 abv 28 Jan 98 see PRO10107, big BSplineCurve C0
243 Standard_Integer nbPini = NCONTROL; // as in BRepCheck_Edge (RLN/Nijni)
244 // 20; // number of points for interpolation, should be "parametric dependent"
245
246 //:92 abv 28 Jan 98: if curve is BSpline with many intervals,
247 // increase number of points to provide at least Degree()+1 points per interval
248 Handle(Geom_BSplineCurve) bspl;
249 if ( c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve)) ) {
250 Handle(Geom_TrimmedCurve) ctrim = Handle(Geom_TrimmedCurve)::DownCast(c3d);
251 bspl = Handle(Geom_BSplineCurve)::DownCast ( ctrim->BasisCurve() );
252 }
253 else bspl = Handle(Geom_BSplineCurve)::DownCast ( c3d );
254 if ( ! bspl.IsNull() ) {
255 Standard_Integer nint = 0;
256 for ( Standard_Integer i=1; i < bspl->NbKnots(); i++ )
2c26a53d 257 {
258 if ( bspl->Knot(i+1) > First && bspl->Knot(i) < Last )
259 nint++;
260 }
7fd59977 261 Standard_Integer minPnt = nint * ( bspl->Degree() + 1 );
262 while ( nbPini < minPnt ) nbPini += NCONTROL - 1;
0797d9d3 263#ifdef OCCT_DEBUG
7fd59977 264 if ( nbPini > NCONTROL )
04232180 265 std::cout << "Warning: number of points for projecting is " << nbPini << std::endl;
7fd59977 266#endif
267 }
268
269// $$$$ end :92 (big BSplineCurve C0)
270
271 // this number should be "parametric dependent"
2c26a53d 272 TColgp_SequenceOfPnt points;
273 TColStd_SequenceOfReal params;
9a6ea9c4 274 NCollection_Sequence<Standard_Real> aKnotCoeffs;
7fd59977 275 gp_Pnt p3d;
9a6ea9c4 276 Standard_Integer iPnt;
277
278 // In case of bspline compute parametrization speed on each
279 // knot interval inside [aFirstParam, aLastParam].
280 // If quotient = (MaxSpeed / MinSpeed) >= aMaxQuotientCoeff then
281 // use PerformByProjLib algorithm.
282 if(!bspl.IsNull())
283 {
9a6ea9c4 284 Standard_Real aFirstParam = First; // First parameter of current interval.
285 Standard_Real aLastParam = Last; // Last parameter of current interval.
286
287 // First index computation.
288 Standard_Integer anIdx = 1;
289 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
290 {
291 if(bspl->Knot(anIdx) > First)
292 {
293 break;
294 }
295 }
296
6a2ee094 297 GeomAdaptor_Curve aC3DAdaptor(c3d);
67285996 298 Standard_Real aMinParSpeed = Precision::Infinite(); // Minimal parameterization speed.
9a6ea9c4 299 for(; anIdx <= bspl->NbKnots() && aFirstParam < Last; anIdx++)
300 {
301 // Fill current knot interval.
302 aLastParam = Min(Last, bspl->Knot(anIdx));
6a2ee094 303 Standard_Integer aNbIntPnts = NCONTROL;
304 // Number of inner points is adapted according to the length of the interval
305 // to avoid a lot of calculations on small range of parameters.
306 if (anIdx > 1)
307 {
308 const Standard_Real aLenThres = 1.e-2;
309 const Standard_Real aLenRatio =
310 (aLastParam - aFirstParam) / (bspl->Knot(anIdx) - bspl->Knot(anIdx - 1));
311 if (aLenRatio < aLenThres)
312 {
313 aNbIntPnts = Standard_Integer(aLenRatio / aLenThres * aNbIntPnts);
314 if (aNbIntPnts < 2)
315 aNbIntPnts = 2;
316 }
317 }
9a6ea9c4 318 Standard_Real aStep = (aLastParam - aFirstParam) / (aNbIntPnts - 1);
319 Standard_Integer anIntIdx;
320 gp_Pnt p3d1, p3d2;
6a2ee094 321 // Start filling from first point.
322 aC3DAdaptor.D0(aFirstParam, p3d1);
323
9a6ea9c4 324 Standard_Real aLength3d = 0.0;
6a2ee094 325 for(anIntIdx = 1; anIntIdx < aNbIntPnts; anIntIdx++)
9a6ea9c4 326 {
6a2ee094 327 Standard_Real aParam = aFirstParam + aStep * anIntIdx;
328 aC3DAdaptor.D0 (aParam, p3d2);
67285996 329 const Standard_Real aDist = p3d2.Distance(p3d1);
330
331 aLength3d += aDist;
6a2ee094 332 p3d1 = p3d2;
67285996 333
334 aMinParSpeed = Min(aMinParSpeed, aDist / aStep);
9a6ea9c4 335 }
6a2ee094 336 const Standard_Real aCoeff = aLength3d / (aLastParam - aFirstParam);
337 if (Abs(aCoeff) > gp::Resolution())
338 aKnotCoeffs.Append(aCoeff);
9a6ea9c4 339 aFirstParam = aLastParam;
340 }
341
342 Standard_Real anEvenlyCoeff = 0;
343 if (aKnotCoeffs.Size() > 0)
344 {
345 anEvenlyCoeff = *std::max_element(aKnotCoeffs.begin(), aKnotCoeffs.end()) /
346 *std::min_element(aKnotCoeffs.begin(), aKnotCoeffs.end());
347 }
348
349 const Standard_Real aMaxQuotientCoeff = 1500.0;
67285996 350 if (anEvenlyCoeff > aMaxQuotientCoeff &&
351 aMinParSpeed > Precision::Confusion() )
9a6ea9c4 352 {
353 PerformByProjLib(c3d, First, Last, c2d);
354 // PerformByProjLib fail detection:
355 if (!c2d.IsNull())
356 {
357 return Status (ShapeExtend_DONE);
358 }
359 }
360 }
361
7fd59977 362 Standard_Real deltaT, t;
363 deltaT = (Last - First) / (nbPini-1);
364 nbrPnt = nbPini;
9a6ea9c4 365 for (iPnt = 1; iPnt <= nbPini; iPnt ++)
366 {
7fd59977 367 if (iPnt == 1) t = First;
368 else if (iPnt == nbPini) t = Last;
369 else t = First + (iPnt - 1) * deltaT;
370
371 c3d->D0 (t, p3d);
2c26a53d 372 points.Append(p3d);
373 params.Append(t);
7fd59977 374 }
375
9a6ea9c4 376 // CALCUL par approximation
2c26a53d 377 TColgp_SequenceOfPnt2d pnt2d;
378 ApproxPCurve (nbrPnt,c3d,TolFirst,TolLast,
379 points,params,pnt2d,c2d); //szv#4:S4163:12Mar99 OK not needed
380 nbPini = points.Length();
7fd59977 381 if (!c2d.IsNull()) {
382 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE2);
383 return Standard_True;
384 }// cas particulier d iso
385
386// INTERPOLATION du resultat
387
388 if ( myBuild ) {
389 Handle(TColgp_HArray1OfPnt) thePnts = new TColgp_HArray1OfPnt (1, nbPini);
390 Handle(TColStd_HArray1OfReal) theParams = new TColStd_HArray1OfReal(1, nbPini);
391 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
392 thePnts->SetValue(iPnt, points(iPnt));
393 theParams->SetValue(iPnt, params(iPnt));
394 }
395
396 Handle(Geom_Curve) newc3d = InterpolateCurve3d (nbPini,thePnts,theParams, c3d);
397 if ( newc3d.IsNull() ) myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
398 else {
399 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE3);
400 c3d = newc3d;
401 }
402 }
403
404 Handle(TColgp_HArray1OfPnt2d) thePnts2d = new TColgp_HArray1OfPnt2d(1, nbPini);
405 Handle(TColStd_HArray1OfReal) theParams2d = new TColStd_HArray1OfReal(1, nbPini);
406 for (iPnt = 1; iPnt <= nbPini ; iPnt ++) {
407 theParams2d->SetValue(iPnt, params(iPnt));
408 thePnts2d->SetValue(iPnt, pnt2d(iPnt));
409 }
410
411 c2d = InterpolatePCurve (nbPini, thePnts2d, theParams2d, c3d);
412// c2d = ApproximatePCurve (nbPini, thePnts2d, theParams2d, c3d);
413// Faut-il aussi reprendre la C3D ?
7fd59977 414 myStatus |= ShapeExtend::EncodeStatus (c2d.IsNull() ? ShapeExtend_FAIL1 : ShapeExtend_DONE2);
415 return Status (ShapeExtend_DONE);
416}
417
418//=======================================================================
419//function : PerformByProjLib
420//purpose :
421//=======================================================================
422
423Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(Handle(Geom_Curve)& c3d,
424 const Standard_Real First,
425 const Standard_Real Last,
426 Handle(Geom2d_Curve)& c2d,
9a6ea9c4 427 const GeomAbs_Shape /*continuity*/,
428 const Standard_Integer /*maxdeg */,
429 const Standard_Integer /*nbinterval */)
7fd59977 430{
431 //Standard_Boolean OK = Standard_True; //szv#4:S4163:12Mar99 unused
432 c2d.Nullify();
433 if (mySurf.IsNull()) {
434 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL1);
435 return Standard_False;
436 }
437
9a6ea9c4 438 try
439 {
7fd59977 440 OCC_CATCH_SIGNALS
441 Handle(GeomAdaptor_HSurface) GAS = mySurf->Adaptor3d();
7fd59977 442 Handle(GeomAdaptor_HCurve) GAC = new GeomAdaptor_HCurve (c3d,First,Last);
9a6ea9c4 443 ProjLib_ProjectedCurve Projector(GAS, GAC);
444
445 switch (Projector.GetType())
446 {
447 case GeomAbs_Line :
448 c2d = new Geom2d_Line(Projector.Line());
449 break;
450 case GeomAbs_Circle :
451 c2d = new Geom2d_Circle(Projector.Circle());
452 break;
453 case GeomAbs_Ellipse :
454 c2d = new Geom2d_Ellipse(Projector.Ellipse());
455 break;
456 case GeomAbs_Parabola :
457 c2d = new Geom2d_Parabola(Projector.Parabola());
458 break;
459 case GeomAbs_Hyperbola :
460 c2d = new Geom2d_Hyperbola(Projector.Hyperbola());
461 break;
462 case GeomAbs_BSplineCurve :
463 c2d = Projector.BSpline();
464 break;
1aec3320 465 default:
9a6ea9c4 466 // Not possible, handling added to avoid gcc warning.
467 break;
7fd59977 468 }
9a6ea9c4 469
470 if(c2d.IsNull())
471 {
7fd59977 472 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL2);
473 return Standard_False;
474 }
9a6ea9c4 475 else
476 {
7fd59977 477 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_DONE1);
478 return Standard_True;
479 }
480
481 }
9775fa61 482 catch(Standard_Failure const& anException) {
0797d9d3 483#ifdef OCCT_DEBUG
04232180 484 std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::PerformByProjLib(): Exception: ";
485 anException.Print(std::cout); std::cout << std::endl;
7fd59977 486#endif
9775fa61 487 (void)anException;
7fd59977 488 myStatus = ShapeExtend::EncodeStatus (ShapeExtend_FAIL3);
489 c2d.Nullify();
490 }
491 return Standard_False;
492}
493
7fd59977 494//=======================================================================
495//function : ProjectAnalytic
496//purpose :
497//=======================================================================
498
499 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ProjectAnalytic(const Handle(Geom_Curve)& c3d) const
500{
501 Handle(Geom2d_Curve) result;
502
503 //:k1 abv 16 Dec 98: limit analytic cases by Plane surfaces only
504 // This is necessary for K4L since it fails on other surfaces
505 // when general method GeomProjLib::Curve2d() is used
506 // Projection is done as in BRep_Tool and BRepCheck_Edge
507 Handle(Geom_Surface) surf = mySurf->Surface();
508 Handle(Geom_Plane) Plane = Handle(Geom_Plane)::DownCast ( surf );
509 if ( Plane.IsNull() ) {
510 Handle(Geom_RectangularTrimmedSurface) RTS =
511 Handle(Geom_RectangularTrimmedSurface)::DownCast ( surf );
512 if ( ! RTS.IsNull() ) Plane = Handle(Geom_Plane)::DownCast ( RTS->BasisSurface() );
513 else {
514 Handle(Geom_OffsetSurface) OS =
515 Handle(Geom_OffsetSurface)::DownCast ( surf );
516 if ( ! OS.IsNull() )
517 Plane = Handle(Geom_Plane)::DownCast ( OS->BasisSurface() );
518 }
519 }
520 if ( ! Plane.IsNull() ) {
521 Handle(Geom_Curve) ProjOnPlane =
522 GeomProjLib::ProjectOnPlane (c3d, Plane,
523 Plane->Position().Direction(), Standard_True);
524 Handle(GeomAdaptor_HCurve) HC = new GeomAdaptor_HCurve ( ProjOnPlane );
525 ProjLib_ProjectedCurve Proj ( mySurf->Adaptor3d(), HC );
526
527 result = Geom2dAdaptor::MakeCurve(Proj);
528 if ( result.IsNull() ) return result;
529 if ( result->IsKind(STANDARD_TYPE(Geom2d_TrimmedCurve)) ) {
530 Handle(Geom2d_TrimmedCurve) TC = Handle(Geom2d_TrimmedCurve)::DownCast ( result );
531 result = TC->BasisCurve();
532 }
63c629aa 533
7fd59977 534 return result;
535 }
536
537 return result;
538}
539
15b54261 540 //! Fix possible period jump and handle walking period parameter.
541 static Standard_Boolean fixPeriodictyTroubles(gp_Pnt2d *thePnt, // pointer to gp_Pnt2d[4] beginning
542 Standard_Integer theIdx, // Index of objective coord: 1 ~ X, 2 ~ Y
02fd709b 543 Standard_Real thePeriod, // Period on objective coord
544 Standard_Integer theSavedPoint, // Point number to choose period
545 Standard_Real theSavedParam) // Param from cashe to choose period
546{
547 Standard_Real aSavedParam;
548 Standard_Integer aSavedPoint;
549 Standard_Real aMinParam = 0.0, aMaxParam = thePeriod;
550 if (theSavedPoint < 0) {
551 // normalize to first period by default
552 aSavedParam = 0.5 * thePeriod;
553 aSavedPoint = 0;
554 }
555 else {
556 aSavedParam = theSavedParam;
557 aSavedPoint = theSavedPoint;
558 while (aMinParam > aSavedParam) {
559 aMinParam -= thePeriod;
560 aMaxParam -= thePeriod;
561 }
562 while (aMaxParam < aSavedParam) {
563 aMinParam += thePeriod;
564 aMaxParam += thePeriod;
565 }
566 }
15b54261 567
02fd709b 568 Standard_Real aFixIsoParam = aMinParam;
569 Standard_Boolean isIsoLine = Standard_False;
570 if (aMaxParam - aSavedParam < Precision::PConfusion() ||
571 aSavedParam - aMinParam < Precision::PConfusion()) {
572 aFixIsoParam = aSavedParam;
573 isIsoLine = Standard_True;
574 }
575 // normalize all coordinates to [aMinParam, aMaxParam)
576 for (Standard_Integer i = 0; i < 4; i++) {
577 Standard_Real aParam = thePnt[i].Coord(theIdx);
578 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(aParam, aMinParam, aMaxParam);
579 aParam += aShift;
580 // Walk over period coord -> not walking on another isoline in parameter space.
581 if (isIsoLine) {
582 if (aMaxParam - aParam < Precision::PConfusion() || aParam - aMinParam < Precision::PConfusion())
583 aParam = aFixIsoParam;
584 }
585 else {
586 if (aMaxParam - aParam < Precision::PConfusion())
587 aParam = aMaxParam;
588 if (aParam - aMinParam < Precision::PConfusion())
589 aParam = aMinParam;
590 }
15b54261 591
02fd709b 592 thePnt[i].SetCoord(theIdx, aParam);
593 }
15b54261 594
02fd709b 595 // find possible period jump and increasing or decreasing coordinates vector
596 Standard_Boolean isJump = Standard_False;
597 Standard_Real aPrevDiff = 0.0;
598 Standard_Real aSumDiff = 1.0;
599 for (Standard_Integer i = 0; i < 3; i++) {
600 Standard_Real aDiff = thePnt[i + 1].Coord(theIdx) - thePnt[i].Coord(theIdx);
601 if (aDiff < -Precision::PConfusion()) {
602 aSumDiff *= -1.0;
603 }
604 //if first derivative changes its sign then period jump may exists in this place
605 if (aDiff * aPrevDiff < -Precision::PConfusion()) {
606 isJump = Standard_True;
607 }
608 aPrevDiff = aDiff;
609 }
15b54261 610
02fd709b 611 if (!isJump)
612 return Standard_False;
613
614 if (aSumDiff > 0) { // decreasing sequence (parameters decrease twice(--) and one period jump(+))
615 for (Standard_Integer i = aSavedPoint; i > 0; i--)
616 if (thePnt[i].Coord(theIdx) > thePnt[i - 1].Coord(theIdx)) {
617 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) + thePeriod);
618 }
619 for (Standard_Integer i = aSavedPoint; i < 3; i++)
620 if (thePnt[i].Coord(theIdx) < thePnt[i + 1].Coord(theIdx)) {
621 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) - thePeriod);
622 }
623 }
624 else {// increasing sequence (parameters increase twice(++) and one period jump(-))
625 for (Standard_Integer i = aSavedPoint; i > 0; i--)
626 if (thePnt[i].Coord(theIdx) < thePnt[i - 1].Coord(theIdx)) {
627 thePnt[i - 1].SetCoord(theIdx, thePnt[i - 1].Coord(theIdx) - thePeriod);
628 }
629 for (Standard_Integer i = aSavedPoint; i < 3; i++)
630 if (thePnt[i].Coord(theIdx) > thePnt[i + 1].Coord(theIdx)) {
631 thePnt[i + 1].SetCoord(theIdx, thePnt[i + 1].Coord(theIdx) + thePeriod);
632 }
633 }
15b54261 634
02fd709b 635 // Do not return false, because for nonlinear 2d curves vector of parameters
636 // may change its first derivative and shifted parameters will be broken for this case.
637 return Standard_True;
15b54261 638 }
639
640//=======================================================================
641//function : getLine
642//purpose :
643//=======================================================================
644
645 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::getLine(
2c26a53d 646 const TColgp_SequenceOfPnt& thepoints,
647 const TColStd_SequenceOfReal& theparams,
648 TColgp_SequenceOfPnt2d& thePnt2ds,
15b54261 649 Standard_Real theTol,
02fd709b 650 Standard_Boolean &isRecompute,
651 Standard_Boolean &isFromCashe) const
15b54261 652 {
653 Standard_Integer nb = thepoints.Length();
654 gp_Pnt aP[4];
655 aP[0] = thepoints(1);
656 aP[1] = thepoints(2);
657 aP[2] = thepoints(nb - 1);
658 aP[3] = thepoints(nb);
659 gp_Pnt2d aP2d[4];
660 Standard_Integer i = 0;
661
662 Standard_Real aTol2 = theTol * theTol;
663 Standard_Boolean isPeriodicU = mySurf->Surface()->IsUPeriodic();
664 Standard_Boolean isPeriodicV = mySurf->Surface()->IsVPeriodic();
665
666 // Workaround:
667 // Protection against bad "tolerance" shapes.
668 if (aTol2 > 1.0)
669 {
670 theTol = Precision::Confusion();
671 aTol2 = theTol * theTol;
672 }
02fd709b 673 if (aTol2 < Precision::SquareConfusion())
674 aTol2 = Precision::SquareConfusion();
15b54261 675 Standard_Real anOldTol2 = aTol2;
02fd709b 676 // auxiliary variables to choose period for connection with previous 2dcurve (if exist)
677 Standard_Integer aSavedPointNum = -1;
678 gp_Pnt2d aSavedPoint;
15b54261 679
680 // project first and last points
681 for( ; i < 4; i +=3)
682 {
683 Standard_Integer j;
2c26a53d 684 for (j = 0; j < myNbCashe; j++)
c48e2889 685 {
15b54261 686 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
687 {
c48e2889 688 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
02fd709b 689 aSavedPointNum = i;
690 aSavedPoint = myCashe2d[j];
691 if (i == 0)
692 isFromCashe = Standard_True;
15b54261 693 break;
694 }
c48e2889 695 }
2c26a53d 696
c48e2889 697 if (j >= myNbCashe)
c48e2889 698 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
2c26a53d 699
c48e2889 700 Standard_Real aDist = mySurf->Gap();
701 Standard_Real aCurDist = aDist * aDist;
2c26a53d 702 if( aTol2 < aDist * aDist)
c48e2889 703 aTol2 = aCurDist;
15b54261 704 }
2c26a53d 705
15b54261 706 if ( isPeriodicU || isPeriodicV )
707 {
708 // Compute second and last but one c2d points.
709 for(i = 1; i < 3; i++)
710 {
711 Standard_Integer j;
2c26a53d 712 for (j = 0; j < myNbCashe; j++)
c48e2889 713 {
15b54261 714 if ( myCashe3d[j].SquareDistance (aP[i] ) < aTol2)
715 {
716 aP2d[i] = mySurf->NextValueOfUV (myCashe2d[j], aP[i], theTol, theTol);
02fd709b 717 aSavedPointNum = i;
718 aSavedPoint = myCashe2d[j];
15b54261 719 break;
720 }
c48e2889 721 }
2c26a53d 722
c48e2889 723 if (j >= myNbCashe)
c48e2889 724 aP2d[i] = mySurf->ValueOfUV(aP[i], theTol);
2c26a53d 725
c48e2889 726 Standard_Real aDist = mySurf->Gap();
727 Standard_Real aCurDist = aDist * aDist;
2c26a53d 728 if( aTol2 < aDist * aDist)
c48e2889 729 aTol2 = aCurDist;
15b54261 730 }
731
732 if (isPeriodicU)
733 {
02fd709b 734 isRecompute = fixPeriodictyTroubles(&aP2d[0], 1 /* X Coord */, mySurf->Surface()->UPeriod(), aSavedPointNum, aSavedPoint.X());
15b54261 735 }
736
737 if (isPeriodicV)
738 {
02fd709b 739 isRecompute = fixPeriodictyTroubles(&aP2d[0], 2 /* Y Coord */, mySurf->Surface()->VPeriod(), aSavedPointNum, aSavedPoint.Y());
15b54261 740 }
741 }
742
02fd709b 743 if (isRecompute && mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
744 // Do not try to make line, because in this case may be very special case when 3d curve
745 // go over the pole of, e.g., sphere, and partly lies along seam
746 // (see ApproxPCurve() for more information).
747 return 0;
748 }
749
15b54261 750 thePnt2ds.SetValue(1, aP2d[0]);
751 thePnt2ds.SetValue(nb, aP2d[3]);
752
753 // Restore old tolerance in 2d space to avoid big gap cases.
754 aTol2 = anOldTol2;
755 // Check that straight line in 2d with parameterisation as in 3d will fit
756 // fit 3d curve at all points.
757 Standard_Real dPar = theparams(nb) - theparams(1);
758 if ( Abs(dPar) < Precision::PConfusion() )
759 return 0;
760 gp_Vec2d aVec0 (aP2d[0], aP2d[3]);
761 gp_Vec2d aVec = aVec0 / dPar;
02fd709b 762 Handle(Geom_Surface) aSurf = mySurf->Surface();
763 Standard_Boolean isNormalCheck = aSurf->IsCNu(1) && aSurf->IsCNv(1);
764 if (isNormalCheck) {
765 for(i = 1; i <= nb; i++)
766 {
767 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
768 gp_Pnt aCurP;
769 gp_Vec aNormalVec, aDu, aDv;
770 aSurf->D1(aCurPoint.X(), aCurPoint.Y(), aCurP, aDu, aDv);
771 aNormalVec = aDu ^ aDv;
772 if (aNormalVec.SquareMagnitude() < Precision::SquareConfusion()) {
773 isNormalCheck = Standard_False;
774 break;
775 }
776 gp_Lin aNormalLine(aCurP, gp_Dir(aNormalVec));
777 Standard_Real aDist = aNormalLine.Distance(thepoints(i));
778 if (aDist > theTol)
779 return 0;
780 }
781 }
782 if (!isNormalCheck) {
783 Standard_Real aFirstPointDist = mySurf->Surface()->Value(aP2d[0].X(), aP2d[0].Y()).
784 SquareDistance(thepoints(1));
785 aTol2 = Max(aTol2, aTol2 * 2 * aFirstPointDist);
786 for(i = 2; i < nb; i++)
787 {
788 gp_XY aCurPoint = aP2d[0].XY() + aVec.XY() * (theparams(i) - theparams(1));
789 gp_Pnt aCurP;
790 aSurf->D0(aCurPoint.X(), aCurPoint.Y(), aCurP);
791 Standard_Real aDist1 = aCurP.SquareDistance(thepoints(i));
792
793 if(Abs (aFirstPointDist - aDist1) > aTol2)
794 return 0;
795 }
15b54261 796 }
797
798 // check if pcurve can be represented by Geom2d_Line (parameterised by length)
799 Standard_Real aLLength = aVec0.Magnitude();
800 if ( Abs (aLLength - dPar) <= Precision::PConfusion() )
801 {
802 gp_XY aDirL = aVec0.XY() / aLLength;
803 gp_Pnt2d aPL (aP2d[0].XY() - theparams(1) * aDirL);
804 return new Geom2d_Line (aPL, gp_Dir2d(aDirL));
805 }
806
807 // create straight bspline
808 TColgp_Array1OfPnt2d aPoles(1, 2);
809 aPoles(1) = aP2d[0];
810 aPoles(2) = aP2d[3];
811
812 TColStd_Array1OfReal aKnots(1,2);
813 aKnots(1) = theparams(1);
814 aKnots(2) = theparams(theparams.Length());
815
816 TColStd_Array1OfInteger aMults(1,2);
817 aMults(1) = 2;
818 aMults(2) = 2;
819 Standard_Integer aDegree = 1;
820 Handle(Geom2d_BSplineCurve) abspl2d =
821 new Geom2d_BSplineCurve (aPoles, aKnots, aMults, aDegree);
822 return abspl2d;
823 }
824
7fd59977 825//=======================================================================
826//function : ApproxPCurve
827//purpose :
828//=======================================================================
829
15b54261 830 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::ApproxPCurve(const Standard_Integer nbrPnt,
2c26a53d 831 const Handle(Geom_Curve)& c3d,
832 const Standard_Real TolFirst,
833 const Standard_Real TolLast,
834 TColgp_SequenceOfPnt& points,
835 TColStd_SequenceOfReal& params,
836 TColgp_SequenceOfPnt2d& pnt2d,
837 Handle(Geom2d_Curve)& c2d)
7fd59977 838{
15b54261 839 // for performance, first try to handle typical case when pcurve is straight
840 Standard_Boolean isRecompute = Standard_False;
02fd709b 841 Standard_Boolean isFromCasheLine = Standard_False;
2c26a53d 842 for (Standard_Integer iseq = 1; iseq <= nbrPnt; iseq++)
843 {
844 gp_Pnt2d aP2d(0.,0.);
845 pnt2d.Append(aP2d);
846 }
02fd709b 847 c2d = getLine(points, params, pnt2d, myPreci, isRecompute, isFromCasheLine);
15b54261 848 if(!c2d.IsNull())
849 {
02fd709b 850 // fill cashe
851 Standard_Boolean ChangeCycle = Standard_False;
852 if(myNbCashe>0 && myCashe3d[0].Distance(points(1)) > myCashe3d[0].Distance(points(nbrPnt)) &&
853 myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
854 ChangeCycle = Standard_True;
855 myNbCashe = 2;
856 if(ChangeCycle) {
857 myCashe3d[0] = points(1);
858 myCashe3d[1] = points(nbrPnt);
859 myCashe2d[0] = pnt2d(1);
860 myCashe2d[1] = pnt2d(nbrPnt);
861 }
862 else {
863 myCashe3d[1] = points(1);
864 myCashe3d[0] = points(nbrPnt);
865 myCashe2d[1] = pnt2d(1);
866 myCashe2d[0] = pnt2d(nbrPnt);
867 }
15b54261 868 return Standard_True;
869 }
870 Standard_Boolean isDone = Standard_True;
7fd59977 871 // test if the curve 3d is a boundary of the surface
872 // (only for Bezier or BSpline surface)
873
874 Standard_Boolean isoParam, isoPar2d3d, isoTypeU, p1OnIso, p2OnIso, isoclosed;
875 gp_Pnt2d valueP1, valueP2;
876 Handle(Geom_Curve) cIso;
877 Standard_Real t1, t2;
878
879 Handle(Standard_Type) sType = mySurf->Surface()->DynamicType();
880 Standard_Boolean isAnalytic = Standard_True;
881 if (sType == STANDARD_TYPE(Geom_BezierSurface) || sType == STANDARD_TYPE(Geom_BSplineSurface)) isAnalytic = Standard_False;
882 Standard_Real uf, ul, vf, vl;
883 mySurf->Surface()->Bounds(uf, ul, vf, vl);
884 isoclosed = Standard_False;
885 TColStd_Array1OfReal pout(1, nbrPnt);
886
887 isoParam = IsAnIsoparametric(nbrPnt, points, params,
888 isoTypeU, p1OnIso, valueP1, p2OnIso, valueP2,
889 isoPar2d3d, cIso, t1, t2, pout);
890
891 // projection of the points on surfaces
892
893 gp_Pnt p3d;
894 gp_Pnt2d p2d;
7fd59977 895 Standard_Real isoValue=0., isoPar1=0., isoPar2=0., tPar=0., tdeb,tfin;
896 Standard_Real Cf, Cl, parf, parl; //szv#4:S4163:12Mar99 dist not needed
897
898 // Le calcul part-il dans le bon sens, c-a-d deb et fin dans le bon ordre ?
899 // Si uclosed et iso en V, attention isoPar1 ET/OU 2 peut toucher la fermeture
900 if(isoParam){
901 if(isoTypeU){
902 isoValue = valueP1.X();
903 isoPar1 = valueP1.Y();
904 isoPar2 = valueP2.Y();
905 isoclosed = mySurf->IsVClosed(myPreci);//#78 rln 12.03.99 S4135
906 parf = vf; parl = vl;
907 }
908 else {
909 isoValue = valueP1.Y();
910 isoPar1 = valueP1.X();
911 isoPar2 = valueP2.X();
912 isoclosed = mySurf->IsUClosed(myPreci);//#78 rln 12.03.99 S4135
913 parf = uf; parl = ul;
914 }
915 if (!isoPar2d3d && !isAnalytic) {
916 Cf = cIso->FirstParameter();
917 Cl = cIso->LastParameter();
918 if (Precision::IsInfinite(Cf)) Cf = -1000;
919 if (Precision::IsInfinite(Cl)) Cl = +1000;
920 //pdn S4030 optimizing and fix isopar case on PRO41323
921 tdeb = pout(2);
922 // dist = ShapeAnalysis_Curve().Project (cIso,points(2),myPreci,pt,tdeb,Cf,Cl);
923 // Chacun des par1 ou par2 est-il sur un bord. Attention first/last : recaler
924 if (isoclosed && (isoPar1 == parf || isoPar1 == parl)) {
925 if (Abs(tdeb-parf) < Abs(tdeb-parl)) isoPar1 = parf;
926 else isoPar1 = parl;
927 if (isoTypeU) valueP1.SetY (isoPar1);
928 else valueP1.SetX (isoPar1);
929 }
930 if (isoclosed && (isoPar2 == parf || isoPar2 == parl)) {
931 //pdn S4030 optimizing and fix isopar case on PRO41323
932 tfin = pout(nbrPnt-1);
933 //dist = ShapeAnalysis_Curve().Project (cIso,points(nbrPnt-1),myPreci,pt,tfin,Cf,Cl);
934 if (Abs(tfin-parf) < Abs(tfin-parl)) isoPar2 = parf;
935 else isoPar2 = parl;
936 if (isoTypeU) valueP2.SetY (isoPar2);
937 else valueP2.SetX (isoPar2);
938 }
939
940 // Interversion Par1/Par2 (ne veut que si les 2 sont sur les bords ...)
941 // Est-ce encore necessaire apres ce qui vient d etre fait ?
942
943 // PTV 05.02.02 fix for translation face from 12_hp_mouse (PARASOLID) face 24008
944 // if curve is periodic do not change the points
945 // skl change "if" for pout(nbrPnt-1) 19.11.2003
946 if (!isoclosed) {
947 if( (Abs(tdeb-isoPar1)>Abs(tdeb-isoPar2)) &&
948 (Abs(pout(nbrPnt-1)-isoPar2)>Abs(pout(nbrPnt-1)-isoPar1)) ) {
949 gp_Pnt2d valueTmp = valueP1;
950 valueP1 = valueP2; valueP2 = valueTmp;
951 if (isoTypeU) {
952 isoValue = valueP1.X();
953 isoPar1 = valueP1.Y();
954 isoPar2 = valueP2.Y();
955 }
956 else {
957 isoValue = valueP1.Y();
958 isoPar1 = valueP1.X();
959 isoPar2 = valueP2.X();
960 }
961 // Fin calcul sens de courbe iso
962 }
963 } // end of fix check 05.02.02
964 }
965 }
966
967 // Si pas isoParam, on a quand meme du p1OnIso/p2OnIso possible ... !!!
968 // (utile pour detromper bug de projection). Mais detromper aussi circularite
969 //else {
970 //if (p1OnIso) valueP1 =
971 //BestExtremum (valueP1,points(1),points(2));
972 //if (p2OnIso) valueP2 =
973 //BestExtremum (valueP2,points(nbrPnt),points(nbrPnt-1));
974 //}
975
976 Standard_Real gap = myPreci; //:q1
977 Standard_Boolean ChangeCycle = Standard_False; //skl for OCC3430
02fd709b 978 // auxiliaruy variables to shift 2dcurve, according to previous
979 Standard_Boolean isFromCashe = Standard_False;
980 gp_Pnt2d aSavedPoint;
7fd59977 981 if( myNbCashe>0 && myCashe3d[0].Distance(points(1))>myCashe3d[0].Distance(points(nbrPnt)) )
982 //if(myCashe3d[0].Distance(points(nbrPnt))<myPreci)
983 if(myCashe3d[0].Distance(points(nbrPnt))<Precision::Confusion())
984 ChangeCycle = Standard_True;
985 //for( i = 1; i <= nbrPnt; i ++) {
986 for(Standard_Integer ii=1; ii<=nbrPnt; ii++) {
2c26a53d 987 const Standard_Integer aPntIndex = ChangeCycle ? (nbrPnt - ii + 1) : ii;
988 p3d = points (aPntIndex);
7fd59977 989 if (isoParam) {
990
991 if (isoPar2d3d) {
2c26a53d 992 if (isoPar2 > isoPar1) tPar = params (aPntIndex);
993 else tPar = t1 + t2 - params(aPntIndex);
7fd59977 994 } else if (!isAnalytic) {
995 // projection to iso
2c26a53d 996 if (aPntIndex == 1) tPar = isoPar1;
997 else if (aPntIndex == nbrPnt) tPar = isoPar2;
7fd59977 998 else {
2c26a53d 999 tPar = pout(aPntIndex);
7fd59977 1000 //:S4030 ShapeAnalysis_Curve().Project (cIso,p3d,myPreci,pt,tPar,Cf,Cl); //szv#4:S4163:12Mar99 `dist=` not needed
1001 }
1002 }
1003
1004 if (!isoPar2d3d && isAnalytic) {
2c26a53d 1005 if (aPntIndex == 1) p2d = valueP1;
1006 else if (aPntIndex == nbrPnt) p2d = valueP2;
7fd59977 1007 else {
1008 p2d = mySurf->NextValueOfUV(p2d,p3d, myPreci, //%12 pdn 15.02.99 optimizing
1009 Precision::Confusion()+1000*gap); //:q1
1010 gap = mySurf->Gap();
1011 }
1012 } else {
1013 if(isoTypeU) { p2d.SetX(isoValue); p2d.SetY(tPar); }
1014 else { p2d.SetX(tPar); p2d.SetY(isoValue); }
1015 }
1016 }
1017
1018 else {
2c26a53d 1019 if ( (aPntIndex == 1) && p1OnIso) p2d = valueP1;
1020 else if( (aPntIndex == nbrPnt) && p2OnIso) p2d = valueP2;
7fd59977 1021 else {// general case (not an iso) mais attention aux singularites !
15b54261 1022 // first and last points are already computed by getLine()
2c26a53d 1023 if (aPntIndex == 1 || aPntIndex == nbrPnt)
15b54261 1024 {
1025 if (!isRecompute)
1026 {
2c26a53d 1027 p2d = pnt2d (aPntIndex);
15b54261 1028 gap = mySurf->Gap();
2c26a53d 1029 if (aPntIndex == 1) {
02fd709b 1030 isFromCashe = isFromCasheLine;
1031 aSavedPoint = p2d;
1032 }
15b54261 1033 continue;
1034 }
1035 else
1036 {
1037 //:q9 abv 23 Mar 99: use cashe as 1st approach
1038 Standard_Integer j; // svv #1
c48e2889 1039 for (j = 0; j < myNbCashe; ++j)
1040 {
15b54261 1041 if ( myCashe3d[j].SquareDistance ( p3d ) < myPreci*myPreci )
1042 {
c48e2889 1043 p2d = mySurf->NextValueOfUV (myCashe2d[j], p3d, myPreci, Precision::Confusion()+gap);
2c26a53d 1044 if (aPntIndex == 1)
c48e2889 1045 {
02fd709b 1046 isFromCashe = Standard_True;
1047 aSavedPoint = myCashe2d[j];
1048 }
15b54261 1049 break;
1050 }
c48e2889 1051 }
1052 if (j >= myNbCashe)
1053 {
1054 p2d = mySurf->ValueOfUV(p3d, myPreci);
1055 }
15b54261 1056 }
1057 }
7fd59977 1058 else {
1059 p2d = mySurf->NextValueOfUV (p2d, p3d, myPreci, //:S4030: optimizing
1060 Precision::Confusion()+1000*gap); //:q1
1061 }
1062 gap = mySurf->Gap();
1063 }
1064 }
2c26a53d 1065 pnt2d (aPntIndex) = p2d;
7fd59977 1066 if ( ii > 1 ) {
1067 if(ChangeCycle)
2c26a53d 1068 p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex + 1).XY() );
7fd59977 1069 else
2c26a53d 1070 p2d.SetXY ( 2. * p2d.XY() - pnt2d(aPntIndex - 1).XY() );
7fd59977 1071 }
1072 }
1073
1074 //pdn %12 11.02.99 PRO9234 entity 15402
1075 if (!isoPar2d3d) {
1076 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_True);
1077 mySurf->ProjectDegenerated(nbrPnt,points,pnt2d,myPreci,Standard_False);
1078 }
2c26a53d 1079
1080 //Check the extremities of 3d curve for coinciding with singularities of surf
1081 //Standard_Integer NbSing = mySurf->NbSingularities(Precision::Confusion());
1082 gp_Pnt PointFirst = points.First(), PointLast = points.Last();
1083 Standard_Real aTolFirst = (TolFirst == -1)? Precision::Confusion() : TolFirst;
1084 Standard_Real aTolLast = (TolLast == -1)? Precision::Confusion() : TolLast;
1085 for (Standard_Integer i = 1; ; i++)
1086 {
1087 Standard_Real aPreci, aFirstPar, aLastPar;
1088 gp_Pnt aP3d;
1089 gp_Pnt2d aFirstP2d, aLastP2d;
1090 Standard_Boolean IsUiso;
1091 if (!mySurf->Singularity(i, aPreci, aP3d, aFirstP2d, aLastP2d, aFirstPar, aLastPar, IsUiso))
1092 break;
1093 if (aPreci <= Precision::Confusion() &&
1094 PointFirst.Distance(aP3d) <= aTolFirst)
1095 {
1096 CorrectExtremity(c3d, params, pnt2d,
1097 Standard_True, //first point
1098 aFirstP2d,
1099 IsUiso);
1100 }
1101 if (aPreci <= Precision::Confusion() &&
1102 PointLast.Distance(aP3d) <= aTolLast)
1103 {
1104 CorrectExtremity(c3d, params, pnt2d,
1105 Standard_False, //last point
1106 aFirstP2d,
1107 IsUiso);
1108 }
1109 }
7fd59977 1110
1111 // attention aux singularites ... (hors cas iso qui les traite deja)
1112 // if (!isoParam) {
1113 // p2d = pnt2d (1);
1114 // if (mySurf->ProjectDegenerated (points(1),myPreci,pnt2d (2),p2d))
1115 // pnt2d (1) = p2d;
1116 // p2d = pnt2d (nbrPnt);
1117 // if (mySurf->ProjectDegenerated (points(nbrPnt),myPreci,pnt2d (nbrPnt-1),p2d))
1118 // pnt2d (nbrPnt) = p2d;
1119 // }
1120
1121 // Si la surface est UCLosed et VClosed, on recadre les points
1122 // algo un peu complique, on retarde l implementation
1123 Standard_Real Up = ul - uf;
1124 Standard_Real Vp = vl - vf;
1125 Standard_Real dist2d;
2c26a53d 1126 const Standard_Real TolOnUPeriod = Precision::Confusion() * Up;
1127 const Standard_Real TolOnVPeriod = Precision::Confusion() * Vp;
0797d9d3 1128#ifdef OCCT_DEBUG
7fd59977 1129 if (mySurf->IsUClosed(myPreci) && mySurf->IsVClosed(myPreci)) {//#78 rln 12.03.99 S4135
04232180 1130 std::cout << "WARNING : Recadrage incertain sur U & VClosed" << std::endl;
7fd59977 1131 }
1132#endif
1133 // Si la surface est UCLosed, on recadre les points
1134 if (mySurf->IsUClosed(myPreci)) {//#78 rln 12.03.99 S4135
1135 // Premier point dans le domain [uf, ul]
1136 Standard_Real prevX, firstX = pnt2d (1).X();
02fd709b 1137 if (!isFromCashe) {
1138 // do not shift 2dcurve, if it connects to previous
1139 while (firstX < uf) { firstX += Up; pnt2d (1).SetX(firstX); }
1140 while (firstX > ul) { firstX -= Up; pnt2d (1).SetX(firstX); }
1141 }
1142 // shift first point, according to cashe
1143 if (mySurf->Surface()->IsUPeriodic() && isFromCashe) {
1144 Standard_Real aMinParam = uf, aMaxParam = ul;
1145 while (aMinParam > aSavedPoint.X()) {
1146 aMinParam -= Up;
1147 aMaxParam -= Up;
1148 }
1149 while (aMaxParam < aSavedPoint.X()) {
1150 aMinParam += Up;
1151 aMaxParam += Up;
1152 }
1153 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstX, aMinParam, aMaxParam);
1154 firstX += aShift;
1155 pnt2d(1).SetX(firstX);
1156 }
7fd59977 1157 prevX = firstX;
1158
1159 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1160 Standard_Real minX = firstX, maxX = firstX;
2c26a53d 1161 Standard_Boolean ToAdjust = Standard_False;
7fd59977 1162
1163 // On decalle toujours le suivant
2c26a53d 1164 for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
1165 {
1166 // dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
1167 Standard_Real CurX = pnt2d (aPntIter).X();
7fd59977 1168 dist2d = Abs (CurX - prevX);
2c26a53d 1169 if (dist2d > ( Up / 2) )
1170 {
1171 InsertAdditionalPointOrAdjust(ToAdjust, 1, Up, TolOnUPeriod,
1172 CurX, prevX,
1173 c3d,
1174 aPntIter,
1175 points, params, pnt2d);
7fd59977 1176 }
1177 prevX = CurX;
1178 if ( minX > CurX ) minX = CurX; //:97
1179 else if ( maxX < CurX ) maxX = CurX; //:97
1180 }
1181
1182 //:97
02fd709b 1183 if (!isFromCashe) {
1184 // do not shift 2dcurve, if it connects to previous
1185 Standard_Real midX = 0.5 * ( minX + maxX );
1186 Standard_Real shiftX=0.;
1187 if ( midX > ul ) shiftX = -Up;
1188 else if ( midX < uf ) shiftX = Up;
1189 if ( shiftX != 0. )
2c26a53d 1190 for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
1191 pnt2d (aPntIter).SetX ( pnt2d (aPntIter).X() + shiftX );
02fd709b 1192 }
7fd59977 1193 }
1194 // Si la surface est VCLosed, on recadre les points
1195 // Same code as UClosed : optimisation souhaitable !!
1196 // CKY : d abord un code IDENTIQUE A UClosed; PUIS le special Seam ...
1197 // Si la surface est UCLosed, on recadre les points
1198 //
1199 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1200 //#78 rln 12.03.99 S4135
1201 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
1202 // Premier point dans le domain [vf, vl]
1203 Standard_Real prevY, firstY = pnt2d (1).Y();
02fd709b 1204 if (!isFromCashe) {
1205 // do not shift 2dcurve, if it connects to previous
1206 while (firstY < vf) { firstY += Vp; pnt2d (1).SetY(firstY); }
1207 while (firstY > vl) { firstY -= Vp; pnt2d (1).SetY(firstY); }
1208 }
1209 // shift first point, according to cashe
1210 if (mySurf->Surface()->IsVPeriodic() && isFromCashe) {
1211 Standard_Real aMinParam = vf, aMaxParam = vl;
1212 while (aMinParam > aSavedPoint.Y()) {
1213 aMinParam -= Vp;
1214 aMaxParam -= Vp;
1215 }
1216 while (aMaxParam < aSavedPoint.Y()) {
1217 aMinParam += Vp;
1218 aMaxParam += Vp;
1219 }
1220 Standard_Real aShift = ShapeAnalysis::AdjustToPeriod(firstY, aMinParam, aMaxParam);
1221 firstY += aShift;
1222 pnt2d(1).SetY(firstY);
1223 }
7fd59977 1224 prevY = firstY;
1225
1226 //:97 abv 1 Feb 98: treat case when curve is whole out of surface bounds
1227 Standard_Real minY = firstY, maxY = firstY;
2c26a53d 1228 Standard_Boolean ToAdjust = Standard_False;
7fd59977 1229
1230 // On decalle toujours le suivant
2c26a53d 1231 for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter)
1232 {
7fd59977 1233 // dist2d = pnt2d (i-1).Distance(pnt2d (i));
2c26a53d 1234 Standard_Real CurY = pnt2d (aPntIter).Y();
7fd59977 1235 dist2d = Abs (CurY - prevY);
2c26a53d 1236 if (dist2d > ( Vp / 2) )
1237 {
1238 InsertAdditionalPointOrAdjust(ToAdjust, 2, Vp, TolOnVPeriod,
1239 CurY, prevY,
1240 c3d,
1241 aPntIter,
1242 points, params, pnt2d);
7fd59977 1243 }
1244 prevY = CurY;
1245 if ( minY > CurY ) minY = CurY; //:97
1246 else if ( maxY < CurY ) maxY = CurY; //:97
1247 }
1248
1249 //:97
02fd709b 1250 if (!isFromCashe) {
1251 // do not shift 2dcurve, if it connects to previous
1252 Standard_Real midY = 0.5 * ( minY + maxY );
1253 Standard_Real shiftY=0.;
1254 if ( midY > vl ) shiftY = -Vp;
1255 else if ( midY < vf ) shiftY = Vp;
1256 if ( shiftY != 0. )
2c26a53d 1257 for (Standard_Integer aPntIter = 1; aPntIter <= pnt2d.Length(); ++aPntIter)
1258 pnt2d(aPntIter).SetY ( pnt2d(aPntIter).Y() + shiftY );
02fd709b 1259 }
7fd59977 1260 }
1261
1262 //#69 rln 01.03.99 S4135 bm2_sd_t4-A.stp entity 30
1263 //#78 rln 12.03.99 S4135
1264 if (mySurf->IsVClosed(myPreci) || mySurf->Surface()->IsKind (STANDARD_TYPE (Geom_SphericalSurface))) {
2c26a53d 1265 for (Standard_Integer aPntIter = 2; aPntIter <= pnt2d.Length(); ++aPntIter) {
7fd59977 1266 //#1 rln 11/02/98 ca_exhaust.stp entity #9869 dist2d = pnt2d (i-1).Distance(pnt2d (i));
2c26a53d 1267 dist2d = Abs (pnt2d (aPntIter).Y() - pnt2d (aPntIter - 1).Y());
7fd59977 1268 if (dist2d > ( Vp / 2) ) {
1269 // ATTENTION : il faut regarder ou le decalage se fait.
1270 // si plusieurs points sont decalles, il faut plusieurs passes
1271 // pour obtenir un resultat correct.
1272 // NOT YET IMPLEMENTED
1273
1274 // one of those point is incorrectly placed
1275 // i.e on the wrong side of the "seam"
1276 // on prend le point le plus pres des bords vf ou vl
1277 Standard_Boolean prevOnFirst = Standard_False;
1278 Standard_Boolean prevOnLast = Standard_False;
1279 Standard_Boolean currOnFirst = Standard_False;
1280 Standard_Boolean currOnLast = Standard_False;
1281
1282 // .X ? plutot .Y , non ?
2c26a53d 1283 Standard_Real distPrevVF = Abs(pnt2d (aPntIter - 1).Y() - vf);
1284 Standard_Real distPrevVL = Abs(pnt2d (aPntIter - 1).Y() - vl);
1285 Standard_Real distCurrVF = Abs(pnt2d (aPntIter).Y() - vf);
1286 Standard_Real distCurrVL = Abs(pnt2d (aPntIter).Y() - vl);
7fd59977 1287
1288 Standard_Real theMin = distPrevVF;
1289 prevOnFirst = Standard_True;
1290 if (distPrevVL < theMin) {
1291 theMin = distPrevVL;
1292 prevOnFirst = Standard_False;
1293 prevOnLast = Standard_True;
1294 }
1295 if (distCurrVF < theMin) {
1296 theMin = distCurrVF;
1297 prevOnFirst = Standard_False;
1298 prevOnLast = Standard_False;
1299 currOnFirst = Standard_True;
1300 }
1301 if (distCurrVL < theMin) {
1302 theMin = distCurrVL;
1303 prevOnFirst = Standard_False;
1304 prevOnLast = Standard_False;
1305 currOnFirst = Standard_False;
1306 currOnLast = Standard_True;
1307 }
1308 // Modifs RLN/Nijni 3-DEC-1997
1309 if (prevOnFirst) {
2c26a53d 1310 // on decalle le point (aPntIter - 1) en V Last
1311 gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vf); // instead of vl RLN/Nijni
1312 pnt2d (aPntIter - 1) = newPrev;
7fd59977 1313 }
1314 else if (prevOnLast) {
2c26a53d 1315 // on decalle le point (aPntIter - 1) en V first
1316 gp_Pnt2d newPrev(pnt2d (aPntIter - 1).X(), vl); // instead of vf RLN/Nijni
1317 pnt2d (aPntIter - 1) = newPrev;
7fd59977 1318 }
1319 else if (currOnFirst) {
2c26a53d 1320 // on decalle le point (aPntIter) en V Last
1321 gp_Pnt2d newCurr(pnt2d (aPntIter).X(),vf); // instead of vl RLN/Nijni
1322 pnt2d (aPntIter) = newCurr;
7fd59977 1323 }
1324 else if (currOnLast) {
2c26a53d 1325 // on decalle le point (aPntIter) en V First
1326 gp_Pnt2d newCurr(pnt2d (aPntIter).X(), vl); // instead of vf RLN/Nijni
1327 pnt2d (aPntIter) = newCurr;
7fd59977 1328 }
1329 // on verifie
0797d9d3 1330#ifdef OCCT_DEBUG
2c26a53d 1331 dist2d = pnt2d (aPntIter - 1).Distance(pnt2d (aPntIter));
7fd59977 1332 if (dist2d > ( Vp / 2) ) {
04232180 1333 std::cout << "Echec dans le recadrage" << std::endl;
7fd59977 1334 }
1335#endif
1336 }
1337 }
1338 }
1339
1340 //:c0 abv 20 Feb 98: treat very special case when 3d curve
1341 // go over the pole of, e.g., sphere, and partly lies along seam.
1342 // 2d representation of such a curve should consist of 3 parts - one on
1343 // regular part of surface (interior), one part along degenerated boundary
1344 // and one along seam.
1345 // Since it cannot be adjusted later by arranging pcurves (curve is single),
1346 // to fix it it is nesessary to have a possibility of adjusting seam
1347 // part of such curve either to left or right boundary of surface.
1348 // Test is performed only if flag AdjustOverDegen is not -1.
1349 // If AdjustOverDegen is True, seam part of curve is adjusted to
1350 // the left, and if False - to the right parametric boundary
1351 // If treated case is detected, flag DONE4 is set to status
1352 // NOTE: currently, precision is Precision::PConfusion() since it
1353 // is enough on encountered example
1354 // (ug_turbine-A.stp from ProSTEP Benchmark #3, entities ##2470 & 5680)
1355 // (r1001_ac.stp from Test Rally #10, face #35027 and others)
1356 if ( myAdjustOverDegen != -1 ) {
1357 if ( mySurf->IsUClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1358 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1359 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1360 // 1st, find gap point (degenerated pole)
1361 Standard_Real PrevX=0.;
1362 Standard_Integer OnBound=0, PrevOnBound=0;
1363 Standard_Integer ind; // svv #1
1364 Standard_Boolean start = Standard_True;
2c26a53d 1365 for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
7fd59977 1366 Standard_Real CurX = pnt2d(ind).X();
1367 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1368 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1369 continue;
1370 OnBound = ( Abs ( Abs ( CurX - 0.5 * ( ul + uf ) ) - Up/2 ) <=
1371 Precision::PConfusion() );
1372 if ( ! start && Abs ( Abs ( CurX - PrevX ) - Up/2 ) <= 0.01*Up )
1373 break;
1374 start = Standard_False;
1375 PrevX = CurX;
1376 PrevOnBound = OnBound;
1377 }
1378 // if found, adjust seam part
2c26a53d 1379 if ( ind <= pnt2d.Length() ) {
7fd59977 1380 PrevX = ( myAdjustOverDegen ? uf : ul );
1381 Standard_Real dU = Up/2 + Precision::PConfusion();
1382 if ( PrevOnBound ) {
1383 pnt2d(ind-1).SetX ( PrevX );
1384 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1385 Standard_Real CurX = pnt2d(j).X();
1386 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1387 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1388 }
1389 }
1390 else if ( OnBound ) {
1391 pnt2d(ind).SetX ( PrevX );
2c26a53d 1392 for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
7fd59977 1393 Standard_Real CurX = pnt2d(j).X();
1394 while ( CurX < PrevX - dU ) pnt2d(j).SetX ( CurX += Up );
1395 while ( CurX > PrevX + dU ) pnt2d(j).SetX ( CurX -= Up );
1396 }
1397 }
1398 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1399 }
1400 }
1401 }
1402 else if ( mySurf->IsVClosed(myPreci) ) {//#78 rln 12.03.99 S4135
1403 mySurf->IsDegenerated ( gp_Pnt(0,0,0), myPreci ); // pour calculer les dgnr
1404 if ( mySurf->NbSingularities(myPreci) > 0 ) { //rln S4135
1405 // 1st, find gap point (degenerated pole)
1406 Standard_Real PrevY=0.;
1407 Standard_Integer OnBound=0, PrevOnBound=0;
1408 Standard_Integer ind; // svv #1
1409 Standard_Boolean start = Standard_True;
2c26a53d 1410 for ( ind=1; ind <= pnt2d.Length(); ind++ ) {
7fd59977 1411 Standard_Real CurY = pnt2d(ind).Y();
1412 // abv 16 Mar 00: trj3_s1-ug.stp #697: ignore points in singularity
1413 if ( mySurf->IsDegenerated ( points(ind), Precision::Confusion() ) )
1414 continue;
1415 OnBound = ( Abs ( Abs ( CurY - 0.5 * ( vl + vf ) ) - Vp/2 ) <=
1416 Precision::PConfusion() );
1417 if ( ! start && Abs ( Abs ( CurY - PrevY ) - Vp/2 ) <= 0.01*Vp )
1418 break;
1419 start = Standard_False;
1420 PrevY = CurY;
1421 PrevOnBound = OnBound;
1422 }
1423 // if found, adjust seam part
2c26a53d 1424 if ( ind <= pnt2d.Length() ) {
7fd59977 1425 PrevY = ( myAdjustOverDegen ? vf : vl );
1426 Standard_Real dV = Vp/2 + Precision::PConfusion();
1427 if ( PrevOnBound ) {
1428 pnt2d(ind-1).SetY ( PrevY );
1429 for ( Standard_Integer j=ind-2; j >0; j-- ) {
1430 Standard_Real CurY = pnt2d(j).Y();
1431 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1432 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1433 }
1434 }
1435 else if ( OnBound ) {
1436 pnt2d(ind).SetY ( PrevY );
2c26a53d 1437 for ( Standard_Integer j=ind+1; j <= pnt2d.Length(); j++ ) {
7fd59977 1438 Standard_Real CurY = pnt2d(j).Y();
1439 while ( CurY < PrevY - dV ) pnt2d(j).SetY ( CurY += Vp );
1440 while ( CurY > PrevY + dV ) pnt2d(j).SetY ( CurY -= Vp );
1441 }
1442 }
1443 myStatus |= ShapeExtend::EncodeStatus (ShapeExtend_DONE4);
1444 }
1445 }
1446 }
1447 }
1448
1449 //:q9: fill cashe
1450 myNbCashe = 2;
1451 if(ChangeCycle) { // msv 10.08.04: avoid using of uninitialised field
1452 //if(myCashe3d[0].Distance(points(1))>Precision::Confusion() &&
1453 // myCashe3d[1].Distance(points(1))>Precision::Confusion()) {
1454 myCashe3d[0] = points(1);
2c26a53d 1455 myCashe3d[1] = points.Last();
7fd59977 1456 myCashe2d[0] = pnt2d(1);
2c26a53d 1457 myCashe2d[1] = pnt2d.Last();
7fd59977 1458 }
1459 else {
1460 myCashe3d[1] = points(1);
2c26a53d 1461 myCashe3d[0] = points.Last();
7fd59977 1462 myCashe2d[1] = pnt2d(1);
2c26a53d 1463 myCashe2d[0] = pnt2d.Last();
7fd59977 1464 }
7fd59977 1465 return isDone;
1466}
1467
1468//=======================================================================
1469//function : ApproximatePCurve
1470//purpose :
1471//=======================================================================
1472
1473Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(const Standard_Integer /*nbrPnt*/,
1474 Handle(TColgp_HArray1OfPnt2d)& points2d,
1475 Handle(TColStd_HArray1OfReal)& params,
1476 const Handle(Geom_Curve)& /*orig*/) const
1477{
1478// Standard_Real resol = Min(mySurf->Adaptor3d()->VResolution(myPreci), mySurf->Adaptor3d()->UResolution(myPreci));
1479 Standard_Real theTolerance2d = myPreci; // (100*nbrPnt);//resol;
1480 Handle(Geom2d_Curve) C2d;
1481 try {
1482 OCC_CATCH_SIGNALS
1483 CheckPoints2d (points2d, params, theTolerance2d);
1484 Standard_Integer numberPnt = points2d->Length();
1485
1486 TColgp_Array1OfPnt points3d(1,numberPnt);
1487 gp_Pnt2d pnt2d;
1488 gp_Pnt pnt;
1489 Standard_Integer i; // svv #1
1490 for( i = 1; i <= numberPnt; i++) {
1491 pnt2d = points2d->Value(i);
1492 pnt.SetCoord(pnt2d.X(),pnt2d.Y(),0);
1493 points3d(i) = pnt;
1494 }
1495
1496 GeomAPI_PointsToBSpline appr(points3d, params->Array1(), 1, 10, GeomAbs_C1, theTolerance2d);
1497 Handle(Geom_BSplineCurve) crv3d = appr.Curve();
1498 Standard_Integer NbPoles = crv3d->NbPoles();
1499 TColgp_Array1OfPnt poles3d (1, NbPoles);
1500 TColgp_Array1OfPnt2d poles2d (1, NbPoles);
1501 crv3d->Poles(poles3d);
1502 for( i = 1; i <= NbPoles; i++) {
1503 pnt2d.SetCoord(poles3d(i).X(),poles3d(i).Y());
1504 poles2d(i) = pnt2d;
1505 }
1506 TColStd_Array1OfReal weights (1,NbPoles);
1507 TColStd_Array1OfInteger multiplicities (1,crv3d->NbKnots());
1508 TColStd_Array1OfReal knots(1,crv3d->NbKnots());
1509 crv3d->Knots(knots);
1510 crv3d->Weights(weights);
1511 crv3d->Multiplicities(multiplicities);
1512 C2d = new Geom2d_BSplineCurve ( poles2d, weights, knots, multiplicities, crv3d->Degree(), crv3d->IsPeriodic());
1513 return C2d;
1514 }
9775fa61 1515 catch(Standard_Failure const& anException) {
1516#ifdef OCCT_DEBUG
1517 //:s5
7fd59977 1518 // debug ...
1519 Standard_Integer nbp = params->Length();
1520 Standard_Integer nb2 = points2d->Length();
04232180 1521 std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::ApproximatePCurve(): Exception: ";
1522 anException.Print(std::cout);
1523 std::cout<<"Pb Geom2dAPI_Approximate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
7fd59977 1524// if (nb2 > nbp) nb2 = nbp;
1525// Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1526// // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1527// dbl.AddReals (rb2,rbp);
1528// for (Standard_Integer i = 1; i <= nb2; i ++) {
1529// gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1530// dbl.AddXYZ (quoi);
1531// }
1532#endif
9775fa61 1533 (void)anException;
7fd59977 1534 C2d.Nullify();
1535 }
1536 return C2d;
1537}
1538
1539//=======================================================================
1540//function : InterpolatePCurve
1541//purpose :
1542//=======================================================================
1543
1544 Handle(Geom2d_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(const Standard_Integer nbrPnt,
1545 Handle(TColgp_HArray1OfPnt2d)& points2d,
1546 Handle(TColStd_HArray1OfReal)& params,
1547 const Handle(Geom_Curve)& /*orig*/) const
1548{
1549 Handle(Geom2d_Curve) C2d; // NULL si echec
1550 Standard_Real theTolerance2d = myPreci / (100 * nbrPnt);
1551 try {
1552 OCC_CATCH_SIGNALS
1553 // on verifie d abord s il n y a pas de points confondus
1554 // si besoin on retouche les valeurs ...
1555 CheckPoints2d (points2d, params, theTolerance2d);
1556 Geom2dAPI_Interpolate myInterPol2d (points2d, params,
1557 Standard_False, theTolerance2d);
1558 myInterPol2d.Perform();
1559 if (myInterPol2d.IsDone()) C2d = myInterPol2d.Curve();
1560 }
9775fa61 1561 catch(Standard_Failure const& anException) {
1562#ifdef OCCT_DEBUG
1563//:s5
7fd59977 1564// // debug ...
1565 Standard_Integer nbp = params->Length();
1566 Standard_Integer nb2 = points2d->Length();
04232180 1567 std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolatePCurve(): Exception: ";
1568 anException.Print(std::cout);
1569 std::cout<<"Pb Geom2dAPI_Interpolate, tol2d="<<theTolerance2d<<" NbParams="<<nbp<<" NbPnts="<<nb2<<std::endl;
7fd59977 1570// if (nb2 > nbp) nb2 = nbp;
1571// Standard_Real rbp,rb2; rbp = nbp; rb2 = nb2;
1572// // dbl.AddString ("NbP2d/NbParams puis X Y Param -> mini");
1573// dbl.AddReals (rb2,rbp);
1574// for (Standard_Integer i = 1; i <= nb2; i ++) {
1575// gp_XYZ quoi (points2d->Value(i).X(),points2d->Value(i).Y(),params->Value(i) );
1576// dbl.AddXYZ (quoi);
1577// }
1578#endif
9775fa61 1579 (void)anException;
7fd59977 1580 C2d.Nullify();
1581 }
1582 return C2d;
1583}
1584
1585//=======================================================================
1586//function : InterpolateCurve3d
1587//purpose :
1588//=======================================================================
1589
1590Handle(Geom_Curve) ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(const Standard_Integer,
1591 Handle(TColgp_HArray1OfPnt)& points,
1592 Handle(TColStd_HArray1OfReal)& params,
1593 const Handle(Geom_Curve)& /*orig*/) const
1594{
1595 Handle(Geom_Curve) C3d; // NULL si echec
1596 try {
1597 OCC_CATCH_SIGNALS
1598 Standard_Real Tol = myPreci;
1599 CheckPoints(points, params, Tol);
1600 GeomAPI_Interpolate myInterPol(points, params, Standard_False, Tol);
1601 myInterPol.Perform();
1602 if (myInterPol.IsDone()) C3d = myInterPol.Curve();
1603 }
9775fa61 1604 catch(Standard_Failure const& anException) {
1605#ifdef OCCT_DEBUG
1606 //:s5
04232180 1607 std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::InterpolateCurve3d(): Exception: ";
1608 anException.Print(std::cout); std::cout << std::endl;
7fd59977 1609#endif
9775fa61 1610 (void)anException;
1611 C3d.Nullify();
7fd59977 1612 }
1613 return C3d;
1614}
1615
2c26a53d 1616//============================================================================================
1617//function : CorrectExtremity
1618//purpose : corrects first or last 2d point of future curve
1619// in the case when it coincids with a singularity of surface
1620//============================================================================================
1621
1622void ShapeConstruct_ProjectCurveOnSurface::CorrectExtremity(const Handle(Geom_Curve)& theC3d,
1623 const TColStd_SequenceOfReal& theParams,
1624 TColgp_SequenceOfPnt2d& thePnt2d,
1625 const Standard_Boolean theIsFirstPoint,
1626 const gp_Pnt2d& thePointOnIsoLine,
1627 const Standard_Boolean theIsUiso)
1628{
1629 Standard_Integer NbPnt = thePnt2d.Length();
1630 Standard_Integer IndCoord = (theIsUiso)? 2 : 1;
1631 Standard_Real SingularityCoord = thePointOnIsoLine.Coord(3-IndCoord);
1632 gp_Pnt2d EndPoint = (theIsFirstPoint)? thePnt2d(1) : thePnt2d(NbPnt);
1633 Standard_Real FinishCoord = EndPoint.Coord(3-IndCoord); //the constant coord of isoline
1634
1635 gp_Dir2d aDir = (theIsUiso)? gp::DY2d() : gp::DX2d();
1636 gp_Lin2d anIsoLine(EndPoint, aDir);
1637 IntRes2d_Domain Dom1, Dom2;
1638
1639 Standard_Boolean IsPeriodic = (theIsUiso)?
1640 mySurf->Surface()->IsVPeriodic() : mySurf->Surface()->IsUPeriodic();
1641
1642 gp_Pnt2d FirstPointOfLine, SecondPointOfLine;
1643 Standard_Real FinishParam, FirstParam, SecondParam;
1644
1645 if (theIsFirstPoint)
1646 {
1647 FirstPointOfLine = thePnt2d(3);
1648 SecondPointOfLine = thePnt2d(2);
1649 FinishParam = theParams(1);
1650 FirstParam = theParams(3);
1651 SecondParam = theParams(2);
1652 }
1653 else //last point
1654 {
1655 FirstPointOfLine = thePnt2d(NbPnt-2);
1656 SecondPointOfLine = thePnt2d(NbPnt-1);
1657 FinishParam = theParams(NbPnt);
1658 FirstParam = theParams(NbPnt-2);
1659 SecondParam = theParams(NbPnt-1);
1660 }
1661
1662 if (SingularityCoord > FinishCoord &&
1663 SecondPointOfLine.Coord(3-IndCoord) > FinishCoord)
1664 return; //the curve passes through the singularity, do nothing
1665 if (SingularityCoord < FinishCoord &&
1666 SecondPointOfLine.Coord(3-IndCoord) < FinishCoord)
1667 return; //the curve passes through the singularity, do nothing
1668 //Check correctness of <EndPoint>
1669 {
1670 const Standard_Real aPrevDist = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
1671 const Standard_Real aCurDist = Abs(EndPoint.Coord(IndCoord) - SecondPointOfLine.Coord(IndCoord));
1672 if (aCurDist <= 2 * aPrevDist)
1673 return;
1674 }
1675
1676 gp_Pnt2d FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
1677 gp_Pnt2d(SecondPointOfLine.X(), FinishCoord); //first approximation of <FinishPoint>
1678
1679 for (;;)
1680 {
1681 if (Abs(SecondPointOfLine.Coord(3-IndCoord) - FinishCoord) <= 2*Precision::PConfusion())
1682 break;
1683
1684 gp_Vec2d aVec(FirstPointOfLine, SecondPointOfLine);
1685 Standard_Real aSqMagnitude = aVec.SquareMagnitude();
1686 if (aSqMagnitude <= 1.e-32)
1687 break;
1688 aDir.SetCoord(aVec.X(), aVec.Y());
1689
1690 gp_Lin2d aLine(FirstPointOfLine, aDir);
1691 IntCurve_IntConicConic Intersector(anIsoLine, Dom1,
1692 aLine, Dom2,
1693 1.e-10, 1.e-10);
1694 if (Intersector.IsDone() && !Intersector.IsEmpty())
1695 {
1696 IntRes2d_IntersectionPoint IntPoint = Intersector.Point(1);
1697 FinishPoint = IntPoint.Value();
1698 }
1699 else
1700 FinishPoint = (theIsUiso)? gp_Pnt2d(FinishCoord, SecondPointOfLine.Y()) :
1701 gp_Pnt2d(SecondPointOfLine.X(), FinishCoord);
1702
1703 gp_Pnt2d PrevPoint = FirstPointOfLine;
1704 FirstPointOfLine = SecondPointOfLine;
1705 FirstParam = SecondParam;
1706 SecondParam = (FirstParam + FinishParam)/2;
1707 if (Abs(SecondParam - FirstParam) <= 2*Precision::PConfusion())
1708 break;
1709 gp_Pnt aP3d;
1710 theC3d->D0(SecondParam, aP3d);
1711 SecondPointOfLine = mySurf->NextValueOfUV(FirstPointOfLine, aP3d,
1712 myPreci, Precision::Confusion());
1713 if (IsPeriodic)
1714 AdjustSecondPointToFirstPoint(FirstPointOfLine, SecondPointOfLine, mySurf->Surface());
1715
1716 //Check <SecondPointOfLine> to be enough close to <FirstPointOfLine>
1717 //because when a projected point is too close to singularity,
1718 //the non-constant coordinate becomes random.
1719 const Standard_Real aPrevDist = Abs(FirstPointOfLine.Coord(IndCoord) - PrevPoint.Coord(IndCoord));
1720 const Standard_Real aCurDist = Abs(SecondPointOfLine.Coord(IndCoord) - FirstPointOfLine.Coord(IndCoord));
1721 if (aCurDist > 2 * aPrevDist)
1722 break;
1723 }
1724
1725 if (theIsFirstPoint)
1726 thePnt2d(1) = FinishPoint;
1727 else
1728 thePnt2d(NbPnt) = FinishPoint;
1729}
1730
1731//============================================================================================
1732//function : InsertAdditionalPointOrAdjust
1733//purpose : If the current point is too far from the previous point
1734// (more than half-period of surface), it can happen in two cases:
1735// 1. Real current step on corresponding coordinate is small, all we need is adjust;
1736// 2. Current step on corresponding coordinate is really bigger than half-period of
1737// surface in this parametric direction, so we must add additional point to exclude
1738// such big intervals between points in 2d space.
1739//============================================================================================
1740
1741void ShapeConstruct_ProjectCurveOnSurface::
1742InsertAdditionalPointOrAdjust(Standard_Boolean& ToAdjust,
1743 const Standard_Integer theIndCoord,
1744 const Standard_Real Period,
1745 const Standard_Real TolOnPeriod,
1746 Standard_Real& CurCoord,
1747 const Standard_Real prevCoord,
1748 const Handle(Geom_Curve)& c3d,
1749 Standard_Integer& theIndex,
1750 TColgp_SequenceOfPnt& points,
1751 TColStd_SequenceOfReal& params,
1752 TColgp_SequenceOfPnt2d& pnt2d)
1753{
1754 Standard_Real CorrectedCurCoord = ElCLib::InPeriod(CurCoord,
1755 prevCoord - Period/2,
1756 prevCoord + Period/2);
1757 if (!ToAdjust)
1758 {
1759 Standard_Real CurPar = params(theIndex);
1760 Standard_Real PrevPar = params(theIndex-1);
1761 Standard_Real MidPar = (PrevPar + CurPar)/2;
1762 gp_Pnt MidP3d;
1763 c3d->D0(MidPar, MidP3d);
1764 gp_Pnt2d MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
1765 Standard_Real MidCoord = MidP2d.Coord(theIndCoord);
1766 MidCoord = ElCLib::InPeriod(MidCoord, prevCoord - Period/2, prevCoord + Period/2);
1767 Standard_Real FirstCoord = prevCoord, LastCoord = CorrectedCurCoord;
1768 if (LastCoord < FirstCoord)
1769 {Standard_Real tmp = FirstCoord; FirstCoord = LastCoord; LastCoord = tmp;}
1770 if (LastCoord - FirstCoord <= TolOnPeriod)
1771 ToAdjust = Standard_True;
1772 else if (FirstCoord <= MidCoord && MidCoord <= LastCoord)
1773 ToAdjust = Standard_True;
1774 else //add mid point
1775 {
1776 //Standard_Real RefU = prevX;
1777 Standard_Boolean Success = Standard_True;
1778 Standard_Real FirstT = PrevPar; //params(i-1)
1779 Standard_Real LastT = CurPar; //params(i)
1780 MidCoord = MidP2d.Coord(theIndCoord);
1781 while (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod ||
1782 Abs(CurCoord - MidCoord) >= Period/2 - TolOnPeriod)
1783 {
1784 if (MidPar - FirstT <= Precision::PConfusion() ||
1785 LastT - MidPar <= Precision::PConfusion())
1786 {
1787 Success = Standard_False;
1788 break; //wrong choice
1789 }
1790 if (Abs(MidCoord - prevCoord) >= Period/2 - TolOnPeriod)
1791 LastT = (FirstT + LastT)/2;
1792 else
1793 FirstT = (FirstT + LastT)/2;
1794 MidPar = (FirstT + LastT)/2;
1795 c3d->D0(MidPar, MidP3d);
1796 MidP2d = mySurf->ValueOfUV(MidP3d, myPreci);
1797 MidCoord = MidP2d.Coord(theIndCoord);
1798 }
1799 if (Success)
1800 {
1801 points.InsertBefore(theIndex, MidP3d);
1802 params.InsertBefore(theIndex, MidPar);
1803 pnt2d.InsertBefore(theIndex, MidP2d);
1804 theIndex++;
1805 }
1806 else
1807 ToAdjust = Standard_True;
1808 } //add mid point
1809 } //if (!ToAdjust)
1810 if (ToAdjust)
1811 {
1812 CurCoord = CorrectedCurCoord;
1813 pnt2d(theIndex).SetCoord (theIndCoord, CurCoord);
1814 }
1815}
1816
7fd59977 1817//=======================================================================
1818//function : CheckPoints
1819//purpose :
1820//=======================================================================
1821
1822 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints(Handle(TColgp_HArray1OfPnt)& points,Handle(TColStd_HArray1OfReal)& params,Standard_Real& preci) const
1823{
1824 Standard_Integer firstElem = points->Lower();
1825 Standard_Integer lastElem = points->Upper();
1826 Standard_Integer i;
1827 Standard_Integer nbPntDropped = 0;
1828 Standard_Integer lastValid = firstElem; // indice of last undropped point
1829
1830 // will store 0 when the point is to be removed, 1 otherwise
1831 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
2c26a53d 1832 for (i = firstElem; i<=lastElem ; i++)
1833 tmpParam.SetValue(i,1);
7ae65f0d 1834 Standard_Real DistMin2 = RealLast();
7fd59977 1835 gp_Pnt Prev = points->Value (lastValid);
1836 gp_Pnt Curr;
1837 for (i = firstElem + 1; i <= lastElem ; i ++) {
1838 Curr = points->Value(i);
7ae65f0d 1839 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1840 if (CurDist2 < gp::Resolution()) { // test 0
7fd59977 1841 nbPntDropped ++;
1842 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1843 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1844 } else {
7ae65f0d 1845 if (CurDist2 < DistMin2)
1846 DistMin2 = CurDist2;
7fd59977 1847 // lastValid becomes the current (i.e. i)
1848 lastValid = i;
1849 Prev = Curr;
1850 }
1851 }
7ae65f0d 1852 if (DistMin2 < RealLast())
1853 preci = 0.9 * Sqrt (DistMin2); // preci est la distance min entre les points on la reduit un peu
1854 if (nbPntDropped == 0)
7fd59977 1855 return;
7ae65f0d 1856
0797d9d3 1857#ifdef OCCT_DEBUG
04232180 1858 std::cout << "Warning : removing 3d points for interpolation" << std::endl;
7fd59977 1859#endif
1860 // Build new HArrays
1861 Standard_Integer newLast = lastElem - nbPntDropped;
1862 if ((newLast - firstElem + 1) < 2) {
0797d9d3 1863#ifdef OCCT_DEBUG
04232180 1864 std::cout << "Too many degenerated points for 3D interpolation" << std::endl;
7fd59977 1865#endif
1866 return;
1867 }
1868 Handle(TColgp_HArray1OfPnt) newPnts =
1869 new TColgp_HArray1OfPnt(firstElem, newLast);
1870 Handle(TColStd_HArray1OfReal) newParams =
1871 new TColStd_HArray1OfReal(firstElem, newLast);
1872 Standard_Integer newCurr = 1;
1873 for (i = firstElem; i<= lastElem ; i++) {
1874 if (tmpParam.Value(i) == 1) {
1875 newPnts->SetValue(newCurr, points->Value(i));
1876 newParams->SetValue(newCurr, params->Value(i));
1877 newCurr ++;
1878 }
1879 }
1880 points = newPnts;
1881 params = newParams;
7fd59977 1882 // on la reduit un peu
1883}
1884
1885//=======================================================================
1886//function : CheckPoints2d
1887//purpose :
1888//=======================================================================
1889
1890 void ShapeConstruct_ProjectCurveOnSurface::CheckPoints2d(Handle(TColgp_HArray1OfPnt2d)& points,
1891 Handle(TColStd_HArray1OfReal)& params,
1892 Standard_Real& preci) const
1893{
1894 Standard_Integer firstElem = points->Lower();
1895 Standard_Integer lastElem = points->Upper();
1896 Standard_Integer i;
1897 Standard_Integer nbPntDropped = 0;
1898 Standard_Integer lastValid = firstElem; // indice of last undropped point
1899
1900 // will store 0 when the point is to be removed, 1 otherwise
1901 TColStd_Array1OfInteger tmpParam(firstElem, lastElem);
1902 for (i = firstElem; i<=lastElem ; i++) {
1903 tmpParam.SetValue(i,1);
1904 }
7ae65f0d 1905 Standard_Real DistMin2 = RealLast();
7fd59977 1906 gp_Pnt2d Prev = points->Value(lastValid);
1907 gp_Pnt2d Curr;
1908 for (i = firstElem + 1; i<=lastElem ; i++) {
1909 Curr = points->Value(i);
7ae65f0d 1910 Standard_Real CurDist2 = Prev.SquareDistance(Curr);
1911 if (CurDist2 < gp::Resolution()) { // test 0
7fd59977 1912 nbPntDropped ++;
1913 if ( i == lastElem ) tmpParam.SetValue(lastValid, 0); // last point kept
1914 else tmpParam.SetValue(i, 0); // current dropped, lastValid unchanged
1915 } else {
7ae65f0d 1916 if (CurDist2 < DistMin2)
1917 DistMin2 = CurDist2;
7fd59977 1918 // lastValid becomes the current (i.e. i)
1919 lastValid = i;
1920 Prev = Curr;
1921 }
1922 }
7ae65f0d 1923 if (DistMin2 < RealLast())
1924 preci = 0.9 * Sqrt (DistMin2);
1925 if (nbPntDropped == 0)
7fd59977 1926 return;
7ae65f0d 1927
0797d9d3 1928#ifdef OCCT_DEBUG
04232180 1929 std::cout << "Warning : removing 2d points for interpolation" << std::endl;
7fd59977 1930#endif
1931 // Build new HArrays
1932 Standard_Integer newLast = lastElem - nbPntDropped;
1933 if ((newLast - firstElem + 1) < 2) {
0797d9d3 1934#ifdef OCCT_DEBUG
04232180 1935 std::cout << "Too many degenerated points for 2D interpolation" << std::endl;
7fd59977 1936#endif
1937 //pdn 12.02.99 S4135 Creating pcurve with minimal length.
1938 tmpParam.SetValue(firstElem,1);
1939 tmpParam.SetValue(lastElem,1);
1940 gp_XY lastPnt = points->Value(lastElem).XY();
1941 lastPnt.Add(gp_XY(preci,preci));
1942 points->SetValue(lastElem,lastPnt);
1943 newLast = firstElem+1;
1944 //return;
1945 }
1946 Handle(TColgp_HArray1OfPnt2d) newPnts =
1947 new TColgp_HArray1OfPnt2d(firstElem, newLast);
1948 Handle(TColStd_HArray1OfReal) newParams =
1949 new TColStd_HArray1OfReal(firstElem, newLast);
1950 Standard_Integer newCurr = 1;
1951 for (i = firstElem; i <= lastElem ; i++) {
1952 if (tmpParam.Value(i) == 1) {
0797d9d3 1953#ifdef OCCT_DEBUG
04232180 1954 std::cout << "Point " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << std::endl;
7fd59977 1955#endif
1956 newPnts->SetValue(newCurr, points->Value(i));
1957 newParams->SetValue(newCurr, params->Value(i));
1958 newCurr ++;
1959 }
1960 else {
0797d9d3 1961#ifdef OCCT_DEBUG
04232180 1962 std::cout << "Removed " << i << " : " << points->Value(i).X() << " " << points->Value(i).Y() << " at param " << params->Value(i) << std::endl;
7fd59977 1963#endif
1964 }
1965 }
1966 points = newPnts;
1967 params = newParams;
7fd59977 1968}
1969
1970//=======================================================================
1971//function : IsAnIsoparametric
1972//purpose :
1973//=======================================================================
1974//:S4030: modified for optimization
1975//:p9 abv 11 Mar 99: PRO7226 #489490: find nearest boundary instead of first one
1976
1977 Standard_Boolean ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(const Standard_Integer nbrPnt,
2c26a53d 1978 const TColgp_SequenceOfPnt& points,
1979 const TColStd_SequenceOfReal& params,
7fd59977 1980 Standard_Boolean& isoTypeU,
1981 Standard_Boolean& p1OnIso,
1982 gp_Pnt2d& valueP1,
1983 Standard_Boolean& p2OnIso,
1984 gp_Pnt2d& valueP2,
1985 Standard_Boolean& isoPar2d3d,
1986 Handle(Geom_Curve)& cIso,
1987 Standard_Real& t1,
1988 Standard_Real& t2,
1989 TColStd_Array1OfReal& pout) const
1990{
1991 try { // RAJOUT
1992 OCC_CATCH_SIGNALS
1993
1994 Standard_Real prec = Precision::Confusion();//myPreci;
1995
1996 Standard_Boolean isoParam = Standard_False;
1997 isoPar2d3d = Standard_False;
1998
1999 Standard_Real U1, U2, V1, V2;
2000 mySurf->Bounds(U1, U2, V1, V2);
2001
2002 if ( mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
2003 Handle(Geom_RectangularTrimmedSurface) sTrim =
2004 Handle(Geom_RectangularTrimmedSurface)::DownCast(mySurf->Surface());
2005 sTrim->Bounds(U1, U2, V1, V2);
2006 }
2007
2008 gp_Pnt pt;
2009 Standard_Integer mpt[2]; mpt[0] = mpt[1] = 0;
d20d815b 2010 Standard_Real t, tpar[2] = { 0.0, 0.0 }, isoValue=0.;
7fd59977 2011 Standard_Real mindist2;
2012 Standard_Real mind2[2];
2013 mindist2 = mind2[0] = mind2[1] = 4*prec*prec;
2014
2015 p1OnIso = Standard_False;
2016 p2OnIso = Standard_False;
2017 const Bnd_Box* aBox = 0;
2018
2019 for (Standard_Integer j=1; (j<=4) /*&& !isoParam*/; j++) {
2020 Standard_Real isoVal=0.;
2021 Standard_Boolean isoU=Standard_False; //szv#4:S4163:12Mar99 `isoU` must be Standard_Boolean
2022 Handle(Geom_Curve) cI;
2023 Standard_Real tt1, tt2;
2024
2025 if (j == 1 ) {
2026 if (Precision::IsInfinite(U1)) continue;
2027 cI = mySurf->UIso(U1);
2028 isoU = Standard_True;
2029 isoVal = U1;
2030 aBox = & mySurf->GetBoxUF();
2031 }
2032 else if (j == 2) {
2033 if (Precision::IsInfinite(U2)) continue;
2034 cI = mySurf->UIso(U2);
2035 isoU = Standard_True;
2036 isoVal = U2;
2037 aBox = & mySurf->GetBoxUL();
2038 }
2039 else if (j == 3) {
2040 if (Precision::IsInfinite(V1)) continue;
2041 cI = mySurf->VIso(V1);
2042 isoU = Standard_False;
2043 isoVal = V1;
2044 aBox = & mySurf->GetBoxVF();
2045 }
2046 else if (j == 4) {
2047 if (Precision::IsInfinite(V2)) continue;
2048 cI = mySurf->VIso(V2);
2049 isoU = Standard_False;
2050 isoVal = V2;
2051 aBox = & mySurf->GetBoxVL();
2052 }
2053 if(cI.IsNull())
2054 continue;
2055
2056 if (isoU) { tt1 = V1; tt2 = V2; }
2057 else { tt1 = U1; tt2 = U2; }
2058
2059 gp_Pnt ext1, ext2;
2060 cI->D0(tt1, ext1);
2061 cI->D0(tt2, ext2);
2062
2063// PATCH CKY 9-JUL-1998 : protection contre singularite
2064 gp_Pnt extmi;
2065 cI->D0( (tt1+tt2)/2,extmi);
2066 if (ext1.IsEqual(ext2,prec) && ext1.IsEqual(extmi,prec)) continue;
2067
2068 Standard_Boolean PtEQext1 = Standard_False;
2069 Standard_Boolean PtEQext2 = Standard_False;
2070
eb1ebea4 2071 Standard_Real currd2[2], tp[2] = {0, 0};
7fd59977 2072 Standard_Integer mp[2];
2073
2074 for (Standard_Integer i=0; i<2; i++) {
2075 mp[i] = 0;
2076 Standard_Integer k = (i == 0 ? 1 : nbrPnt);
2077
2078 // si ext1 == ext2 => valueP1 == valueP2 => vect null plus tard
2079 currd2[i] = points(k).SquareDistance ( ext1 );
2080 if ( currd2[i] <= prec*prec && !PtEQext1) {
2081 mp[i] = 1;
2082 tp[i] = tt1;
2083 PtEQext1 = Standard_True;
2084 continue;
2085 }
2086
2087 currd2[i] = points(k).SquareDistance ( ext2 );
2088 if ( currd2[i] <= prec*prec && !PtEQext2) {
2089 mp[i] = 2;
2090 tp[i] = tt2;
2091 PtEQext2 = Standard_True;
2092 continue;
2093 }
2094
2095 // On evite de projecter sur un iso degenere
2096 // on doit egalement le faire pour l apex du cone
2097 if (mySurf->Surface()->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) && !isoU) {
2098 continue;
2099 }
2100
2101 if(aBox->IsOut(points(k))) continue;
2102
2103 Standard_Real Cf = cI->FirstParameter();
2104 Standard_Real Cl = cI->LastParameter();
2105 if (Precision::IsInfinite(Cf)) Cf = -1000;
2106 if (Precision::IsInfinite(Cl)) Cl = +1000;
2107
2108 ShapeAnalysis_Curve sac;
2109 Standard_Real dist = sac.Project (cI,points(k),prec,pt,t,Cf,Cl);
2110 currd2[i] = dist * dist;
2111 if ((dist <= prec) && (t>= Cf) && (t<=Cl)) {
2112 mp[i] = 3;
2113 tp[i] = t;
2114 }
2115 }
2116
2117 //:e7 abv 21 Apr 98: ProSTEP TR8, r0501_pe #56679:
2118 // avoid possible null-length curves
2119 if ( mp[0] >0 && mp[1] >0 &&
2120 Abs ( tp[0] - tp[1] ) < Precision::PConfusion() ) continue;
2121
2122
2123 if (mp[0] > 0 &&
2124 ( ! p1OnIso || currd2[0] < mind2[0] ) ) {
2125 p1OnIso = Standard_True;
15b54261 2126 mind2[0] = currd2[0]; // LP2.stp #105899: FLT_INVALID_OPERATION on Windows 7 VC 9 Release mode on the whole file
7fd59977 2127 if (isoU) valueP1.SetCoord(isoVal, tp[0]);
2128 else valueP1.SetCoord(tp[0], isoVal);
2129 }
2130
2131 if (mp[1] > 0 &&
2132 ( ! p2OnIso || currd2[1] < mind2[1] ) ) {
2133 p2OnIso = Standard_True;
2134 mind2[1] = currd2[1];
2135 if (isoU) valueP2.SetCoord(isoVal, tp[1]);
2136 else valueP2.SetCoord(tp[1], isoVal);
2137 }
2138
2139 if ( mp[0] <=0 || mp[1] <=0 ) continue;
2140
2141 Standard_Real md2 = currd2[0] + currd2[1];
2142 if ( mindist2 <= md2 ) continue;
2143
2144 mindist2 = md2;
2145 mpt[0] = mp[0];
2146 mpt[1] = mp[1];
2147 tpar[0] = tp[0];
2148 tpar[1] = tp[1];
2149 isoTypeU = isoU;
2150 isoValue = isoVal;
2151 cIso = cI;
2152 t1 = tt1;
2153 t2 = tt2;
2154 }
2155
2156 // probablely it concerns an isoparametrics
2157 if ( mpt[0] >0 && mpt[1] >0 ) {
2158
2159 p1OnIso = p2OnIso = Standard_True;
2160 if (isoTypeU) {
2161 valueP1.SetCoord(isoValue, tpar[0]);
2162 valueP2.SetCoord(isoValue, tpar[1]);
2163 }
2164 else {
2165 valueP1.SetCoord(tpar[0], isoValue);
2166 valueP2.SetCoord(tpar[1], isoValue);
2167 }
2168
2169 if ( mpt[0] != 3 && mpt[1] != 3 ) {
2170 isoPar2d3d = Standard_True;
2171 for (Standard_Integer i=2; i < nbrPnt && isoPar2d3d; i++){
2172 if (tpar[1] > tpar[0]) t = params(i);
2173 else t = t1+t2-params(i);
2174 cIso->D0(t, pt);
2175 if (!points(i).IsEqual(pt, prec)) isoPar2d3d = Standard_False;
2176 }
2177 }
2178
2179 if (isoPar2d3d) isoParam = Standard_True;
2180 else {
2181 Standard_Real prevParam = tpar[0];
2182 Standard_Real Cf, Cl;
2183 Standard_Boolean isoByDistance = Standard_True;
2184 Cf = cIso->FirstParameter();
2185 Cl = cIso->LastParameter();
2186 if (Precision::IsInfinite(Cf)) Cf = -1000;
2187 if (Precision::IsInfinite(Cl)) Cl = +1000;
2188
2189 ShapeAnalysis_Curve sac;
2190 for (Standard_Integer i=2; i < nbrPnt && isoByDistance; i++) {
2191 Standard_Real dist = sac.NextProject (prevParam,cIso,points(i),
2192 prec,pt,t,Cf,Cl,
2193 Standard_False); //:j8 abv 10.12.98: TR10 r0501_db.stp #9423: avoid adjusting to ends
2194 prevParam = t;
2195 pout(i)=t;
2196 if( (dist > prec) || (t < Cf) || (t > Cl) )
2197 isoByDistance = Standard_False;
2198 }
2199 if (isoByDistance) isoParam = Standard_True;
2200 }
2201 }
2202/* if (!isoParam) { CKY 29-mai-1997 : garder tout ce qu on peut ?
2203 p1OnIso = Standard_False;
2204 p2OnIso = Standard_False;
2205 } */
2206 return isoParam;
2207 } // RAJOUT
9775fa61 2208 catch(Standard_Failure const& anException) {
0797d9d3 2209#ifdef OCCT_DEBUG
9775fa61 2210// pb : on affiche ce qu on peut
7fd59977 2211 for (Standard_Integer numpnt = 1; numpnt <= nbrPnt; numpnt ++) {
04232180 2212 std::cout<<"["<<numpnt<<"]param="<<params(numpnt)<<" point=("<<
2213 points(numpnt).X()<<" "<<points(numpnt).Y()<<" "<<points(numpnt).Z()<<")"<<std::endl;
7fd59977 2214 }
04232180 2215 std::cout << "Warning: ShapeConstruct_ProjectCurveOnSurface::IsAnIsoparametric(): Exception: ";
2216 anException.Print(std::cout); std::cout << std::endl;
7fd59977 2217#endif
9775fa61 2218 (void)anException;
7fd59977 2219 return Standard_False;
2220 }
7fd59977 2221}
2222
2223/* S4135 : BestExtremum is commented after IsAnIsoparametric works with Precision::Confusion()
2224//=======================================================================
2225//function : BestExtremum
2226//purpose : auxiliaire prenant le meilleur extremum si ISO car doute possible
2227//=======================================================================
2228
2229 gp_Pnt2d ShapeConstruct_ProjectCurveOnSurface::BestExtremum(const gp_Pnt2d& P2iso,const gp_Pnt& P3ext,const gp_Pnt& P3next) const
2230{
2231// P2iso a ete calcule depuis P3ext sur une iso externe de la surface
2232// En principe bon mais circularite possible ... et IsU/VClosed faillible
2233// (si baillement 1e-4 ou 1e-5, on est dedans !). DONC
2234// 1/ on privilegie l iso mais a tout hasard on verifie si Surf meilleur
2235// 2/ si iso, attention a la circularite (cas limite)
2236
2237// NB : si isoParam, on suppose que P2iso est bon (car il y en a 2). A voir...
2238
2239// D abord, calcul p2ext depuis la surface. choix surface/iso
2240 return P2iso;
2241 Standard_Real prec = Precision::Confusion();//myPreci;
2242 gp_Pnt2d P2cal = mySurf->ValueOfUV(P3ext, prec);
2243 gp_Pnt P3cal = mySurf->Value (P2cal);
2244 Standard_Real dcal = P3ext.Distance (P3cal);
2245 Standard_Real dnxt = P3ext.Distance (P3next);
2246 if (dcal > dnxt) return P2iso; // en fait protection sur BUG (PRO8468)
2247
2248// On choisit entre P2iso et P2cal, le plus proche de P2next ... !!!
2249 gp_Pnt2d P2next = mySurf->ValueOfUV(P3next, prec);
2250 if (P2next.Distance(P2cal) < P2next.Distance(P2iso)) return P2cal;
2251 return P2iso;
2252}
2253*/
Open CASCADE Technology repository