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