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