projects / occt.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
0026196: Wrong result obtained by projection algorithm.
[occt.git] / src / ProjLib / ProjLib_ComputeApproxOnPolarSurface.cxx
... / ...
CommitLineData
1// Created by: Bruno DUMORTIER
2// Copyright (c) 1995-1999 Matra Datavision
3// Copyright (c) 1999-2014 OPEN CASCADE SAS
4//
5// This file is part of Open CASCADE Technology software library.
6//
7// This library is free software; you can redistribute it and/or modify it under
8// the terms of the GNU Lesser General Public License version 2.1 as published
9// by the Free Software Foundation, with special exception defined in the file
10// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11// distribution for complete text of the license and disclaimer of any warranty.
12//
13// Alternatively, this file may be used under the terms of Open CASCADE
14// commercial license or contractual agreement.
15
16#include <ProjLib_ComputeApproxOnPolarSurface.hxx>
17#include <AppCont_Function.hxx>
18#include <ElSLib.hxx>
19#include <ElCLib.hxx>
20#include <BSplCLib.hxx>
21#include <PLib.hxx>
22#include <Standard_NoSuchObject.hxx>
23#include <Geom_UndefinedDerivative.hxx>
24#include <gp_Trsf.hxx>
25#include <Precision.hxx>
26#include <Approx_FitAndDivide2d.hxx>
27#include <math.hxx>
28#include <AppParCurves_MultiCurve.hxx>
29#include <Geom_Surface.hxx>
30#include <Geom2d_BSplineCurve.hxx>
31#include <Geom2d_BezierCurve.hxx>
32#include <Geom2d_Line.hxx>
33#include <Geom2d_Circle.hxx>
34#include <Geom2d_Ellipse.hxx>
35#include <Geom2d_Hyperbola.hxx>
36#include <Geom2d_Parabola.hxx>
37#include <Geom2d_TrimmedCurve.hxx>
38#include <Geom_BSplineSurface.hxx>
39#include <Geom_BezierSurface.hxx>
40#include <Geom_BSplineCurve.hxx>
41#include <Geom_BezierCurve.hxx>
42#include <Geom_TrimmedCurve.hxx>
43
44#include <TColgp_Array1OfPnt2d.hxx>
45#include <TColgp_Array2OfPnt2d.hxx>
46#include <TColgp_Array1OfPnt.hxx>
47#include <TColgp_SequenceOfPnt2d.hxx>
48#include <TColStd_Array1OfReal.hxx>
49#include <TColStd_Array1OfInteger.hxx>
50#include <TColStd_SequenceOfReal.hxx>
51#include <TColStd_ListOfTransient.hxx>
52
53#include <GeomAbs_SurfaceType.hxx>
54#include <GeomAbs_CurveType.hxx>
55#include <Adaptor3d_Surface.hxx>
56#include <Adaptor3d_Curve.hxx>
57#include <Adaptor3d_HSurface.hxx>
58#include <Adaptor3d_HCurve.hxx>
59#include <Adaptor2d_HCurve2d.hxx>
60#include <Geom2dAdaptor_Curve.hxx>
61#include <Geom2dAdaptor_HCurve.hxx>
62#include <GeomAdaptor_HCurve.hxx>
63#include <GeomAdaptor.hxx>
64#include <GeomAdaptor_Surface.hxx>
65#include <TColgp_SequenceOfPnt.hxx>
66
67#include <gp_Pnt.hxx>
68#include <gp_Pnt2d.hxx>
69#include <gp_Vec2d.hxx>
70#include <Extrema_GenLocateExtPS.hxx>
71#include <Extrema_ExtPS.hxx>
72#include <GCPnts_QuasiUniformAbscissa.hxx>
73#include <Standard_DomainError.hxx>
74//#include <GeomLib_IsIso.hxx>
75//#include <GeomLib_CheckSameParameter.hxx>
76
77#ifdef OCCT_DEBUG
78#ifdef DRAW
79#include <DrawTrSurf.hxx>
80#include <Geom2d_Curve.hxx>
81#endif
82//static Standard_Integer compteur = 0;
83#endif
84
85struct aFuncStruct
86{
87 Handle(Adaptor3d_HSurface) mySurf; // Surface where to project.
88 Handle(Adaptor3d_HCurve) myCurve; // Curve to project.
89 Handle(Adaptor2d_HCurve2d) myInitCurve2d; // Initial 2dcurve projection.
90 Standard_Real mySqProjOrtTol; // Used to filter non-orthogonal projected point.
91 Standard_Real myTolU;
92 Standard_Real myTolV;
93 Standard_Real myPeriod[2]; // U and V period correspondingly.
94};
95
96//=======================================================================
97//function : computePeriodicity
98//purpose : Compute period information on adaptor.
99//=======================================================================
100static void computePeriodicity(const Handle(Adaptor3d_HSurface)& theSurf,
101 Standard_Real &theUPeriod,
102 Standard_Real &theVPeriod)
103{
104 theUPeriod = 0.0;
105 theVPeriod = 0.0;
106
107 // Compute once information about periodicity.
108 // Param space may be reduced in case of rectangular trimmed surface,
109 // in this case really trimmed bounds should be set as unperiodic.
110 Standard_Real aTrimF, aTrimL, aBaseF, aBaseL, aDummyF, aDummyL;
111 Handle(Geom_Surface) aS = GeomAdaptor::MakeSurface(theSurf->Surface(), Standard_False); // Not trim.
112 // U param space.
113 if (theSurf->IsUPeriodic())
114 {
115 theUPeriod = theSurf->UPeriod();
116 }
117 else if(theSurf->IsUClosed())
118 {
119 theUPeriod = theSurf->LastUParameter() - theSurf->FirstUParameter();
120 }
121 if (theUPeriod != 0.0)
122 {
123 aTrimF = theSurf->FirstUParameter(); // Trimmed first
124 aTrimL = theSurf->LastUParameter(); // Trimmed last
125 aS->Bounds(aBaseF, aBaseL, aDummyF, aDummyL); // Non-trimmed values.
126 if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
127 {
128 // Param space reduced.
129 theUPeriod = 0.0;
130 }
131 }
132
133 // V param space.
134 if (theSurf->IsVPeriodic())
135 {
136 theVPeriod = theSurf->VPeriod();
137 }
138 else if(theSurf->IsVClosed())
139 {
140 theVPeriod = theSurf->LastVParameter() - theSurf->FirstVParameter();
141 }
142 if (theVPeriod != 0.0)
143 {
144 aTrimF = theSurf->FirstVParameter(); // Trimmed first
145 aTrimL = theSurf->LastVParameter(); // Trimmed last
146 aS->Bounds(aDummyF, aDummyL, aBaseF, aBaseL); // Non-trimmed values.
147 if (Abs (aBaseF - aTrimF) + Abs (aBaseL - aTrimL) > Precision::PConfusion())
148 {
149 // Param space reduced.
150 theVPeriod = 0.0;
151 }
152 }
153}
154
155//=======================================================================
156//function : aFuncValue
157//purpose : compute functional value in (theU,theV) point
158//=======================================================================
159static Standard_Real anOrthogSqValue(const gp_Pnt& aBasePnt,
160 const Handle(Adaptor3d_HSurface)& Surf,
161 const Standard_Real theU,
162 const Standard_Real theV)
163{
164 // Since find projection, formula is:
165 // F1 = Dot(S_U, Vec(aBasePnt, aProjPnt))
166 // F2 = Dot(S_V, Vec(aBasePnt, aProjPnt))
167
168 gp_Pnt aProjPnt;
169 gp_Vec aSu, aSv;
170
171 Surf->D1(theU, theV, aProjPnt, aSu, aSv);
172 gp_Vec aBaseVec(aBasePnt, aProjPnt);
173
174 if (aSu.SquareMagnitude() > Precision::SquareConfusion())
175 aSu.Normalize();
176
177 if (aSv.SquareMagnitude() > Precision::SquareConfusion())
178 aSv.Normalize();
179
180 Standard_Real aFirstPart = aSu.Dot(aBaseVec);
181 Standard_Real aSecondPart = aSv.Dot(aBaseVec);
182 return (aFirstPart * aFirstPart + aSecondPart * aSecondPart);
183}
184
185//=======================================================================
186//function : Value
187//purpose : (OCC217 - apo)- Compute Point2d that project on polar surface(<Surf>) 3D<Curve>
188// <InitCurve2d> use for calculate start 2D point.
189//=======================================================================
190static gp_Pnt2d Function_Value(const Standard_Real theU,
191 const aFuncStruct& theData)
192{
193 gp_Pnt2d p2d = theData.myInitCurve2d->Value(theU) ;
194 gp_Pnt p = theData.myCurve->Value(theU);
195 gp_Pnt aSurfPnt = theData.mySurf->Value(p2d.X(), p2d.Y());
196 Standard_Real aSurfPntDist = aSurfPnt.SquareDistance(p);
197
198 Standard_Real Uinf, Usup, Vinf, Vsup;
199 Uinf = theData.mySurf->Surface().FirstUParameter();
200 Usup = theData.mySurf->Surface().LastUParameter();
201 Vinf = theData.mySurf->Surface().FirstVParameter();
202 Vsup = theData.mySurf->Surface().LastVParameter();
203
204 // Check case when curve is close to co-parametrized isoline on surf.
205 if (Abs (p2d.X() - Uinf) < Precision::PConfusion() ||
206 Abs (p2d.X() - Usup) < Precision::PConfusion() )
207 {
208 // V isoline.
209 gp_Pnt aPnt;
210 theData.mySurf->D0(p2d.X(), theU, aPnt);
211 if (aPnt.SquareDistance(p) < aSurfPntDist)
212 p2d.SetY(theU);
213 }
214
215 if (Abs (p2d.Y() - Vinf) < Precision::PConfusion() ||
216 Abs (p2d.Y() - Vsup) < Precision::PConfusion() )
217 {
218 // U isoline.
219 gp_Pnt aPnt;
220 theData.mySurf->D0(theU, p2d.Y(), aPnt);
221 if (aPnt.SquareDistance(p) < aSurfPntDist)
222 p2d.SetX(theU);
223 }
224
225 Standard_Integer decalU = 0, decalV = 0;
226 Standard_Real U0 = p2d.X(), V0 = p2d.Y();
227
228 GeomAbs_SurfaceType Type = theData.mySurf->GetType();
229 if((Type != GeomAbs_BSplineSurface) &&
230 (Type != GeomAbs_BezierSurface) &&
231 (Type != GeomAbs_OffsetSurface) )
232 {
233 // Analytical cases.
234 Standard_Real S = 0., T = 0.;
235 switch (Type)
236 {
237 case GeomAbs_Cylinder:
238 {
239 gp_Cylinder Cylinder = theData.mySurf->Cylinder();
240 ElSLib::Parameters( Cylinder, p, S, T);
241 if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
242 if(U0 > Usup) decalU = int((U0 - Usup)/(2*M_PI))+1;
243 S += decalU*2*M_PI;
244 break;
245 }
246 case GeomAbs_Cone:
247 {
248 gp_Cone Cone = theData.mySurf->Cone();
249 ElSLib::Parameters( Cone, p, S, T);
250 if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
251 if(U0 > Usup) decalU = int((U0 - Usup)/(2*M_PI))+1;
252 S += decalU*2*M_PI;
253 break;
254 }
255 case GeomAbs_Sphere:
256 {
257 gp_Sphere Sphere = theData.mySurf->Sphere();
258 ElSLib::Parameters( Sphere, p, S, T);
259 if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
260 if(U0 > Usup) decalU = int((U0 - Usup)/(2*M_PI))+1;
261 S += decalU*2*M_PI;
262 if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
263 if(V0 > (Vsup+(Vsup-Vinf))) decalV = int((V0 - Vsup+(Vsup-Vinf))/(2*M_PI))+1;
264 T += decalV*2*M_PI;
265 if(0.4*M_PI < Abs(U0 - S) && Abs(U0 - S) < 1.6*M_PI)
266 {
267 T = M_PI - T;
268 if(U0 < S)
269 S -= M_PI;
270 else
271 S += M_PI;
272 }
273 break;
274 }
275 case GeomAbs_Torus:
276 {
277 gp_Torus Torus = theData.mySurf->Torus();
278 ElSLib::Parameters( Torus, p, S, T);
279 if(U0 < Uinf) decalU = -int((Uinf - U0)/(2*M_PI))-1;
280 if(U0 > Usup) decalU = int((U0 - Usup)/(2*M_PI))+1;
281 if(V0 < Vinf) decalV = -int((Vinf - V0)/(2*M_PI))-1;
282 if(V0 > Vsup) decalV = int((V0 - Vsup)/(2*M_PI))+1;
283 S += decalU*2*M_PI; T += decalV*2*M_PI;
284 break;
285 }
286 default:
287 Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::Value");
288 }
289 return gp_Pnt2d(S, T);
290 }
291
292 // Non-analytical case.
293 Standard_Real Dist2Min = RealLast();
294 Standard_Real uperiod = theData.myPeriod[0],
295 vperiod = theData.myPeriod[1],
296 u, v;
297
298 // U0 and V0 are the points within the initialized period.
299 if(U0 < Uinf)
300 {
301 if(!uperiod)
302 U0 = Uinf;
303 else
304 {
305 decalU = int((Uinf - U0)/uperiod)+1;
306 U0 += decalU*uperiod;
307 }
308 }
309 if(U0 > Usup)
310 {
311 if(!uperiod)
312 U0 = Usup;
313 else
314 {
315 decalU = -(int((U0 - Usup)/uperiod)+1);
316 U0 += decalU*uperiod;
317 }
318 }
319 if(V0 < Vinf)
320 {
321 if(!vperiod)
322 V0 = Vinf;
323 else
324 {
325 decalV = int((Vinf - V0)/vperiod)+1;
326 V0 += decalV*vperiod;
327 }
328 }
329 if(V0 > Vsup)
330 {
331 if(!vperiod)
332 V0 = Vsup;
333 else
334 {
335 decalV = -int((V0 - Vsup)/vperiod)-1;
336 V0 += decalV*vperiod;
337 }
338 }
339
340 // The surface around (U0,V0) is reduced.
341 Standard_Real uLittle = (Usup - Uinf)/10, vLittle = (Vsup - Vinf)/10;
342 Standard_Real uInfLi = 0, vInfLi = 0,uSupLi = 0, vSupLi = 0;
343 if((U0 - Uinf) > uLittle) uInfLi = U0 - uLittle; else uInfLi = Uinf;
344 if((V0 - Vinf) > vLittle) vInfLi = V0 - vLittle; else vInfLi = Vinf;
345 if((Usup - U0) > uLittle) uSupLi = U0 + uLittle; else uSupLi = Usup;
346 if((Vsup - V0) > vLittle) vSupLi = V0 + vLittle; else vSupLi = Vsup;
347
348 GeomAdaptor_Surface SurfLittle;
349 if (Type == GeomAbs_BSplineSurface)
350 {
351 Handle(Geom_Surface) GBSS(theData.mySurf->Surface().BSpline());
352 SurfLittle.Load(GBSS, uInfLi, uSupLi, vInfLi, vSupLi);
353 }
354 else if (Type == GeomAbs_BezierSurface)
355 {
356 Handle(Geom_Surface) GS(theData.mySurf->Surface().Bezier());
357 SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
358 }
359 else if (Type == GeomAbs_OffsetSurface)
360 {
361 Handle(Geom_Surface) GS = GeomAdaptor::MakeSurface(theData.mySurf->Surface());
362 SurfLittle.Load(GS, uInfLi, uSupLi, vInfLi, vSupLi);
363 }
364 else
365 {
366 Standard_NoSuchObject::Raise("");
367 }
368
369 // Try to run simple search with initial point (U0, V0).
370 Extrema_GenLocateExtPS locext(p, SurfLittle, U0, V0, theData.myTolU, theData.myTolV);
371 if (locext.IsDone())
372 {
373 locext.Point().Parameter(u, v);
374 Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
375 if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
376 locext.SquareDistance() < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
377 {
378 gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
379 return pnt;
380 }
381 }
382
383 // Perform whole param space search.
384 Extrema_ExtPS ext(p, SurfLittle, theData.myTolU, theData.myTolV);
385 if (ext.IsDone() && ext.NbExt() >= 1)
386 {
387 Dist2Min = ext.SquareDistance(1);
388 Standard_Integer GoodValue = 1;
389 for (Standard_Integer i = 2 ; i <= ext.NbExt() ; i++ )
390 {
391 if( Dist2Min > ext.SquareDistance(i))
392 {
393 Dist2Min = ext.SquareDistance(i);
394 GoodValue = i;
395 }
396 }
397 ext.Point(GoodValue).Parameter(u, v);
398 Dist2Min = anOrthogSqValue(p, theData.mySurf, u, v);
399 if (Dist2Min < theData.mySqProjOrtTol && // Point is projection.
400 ext.SquareDistance(GoodValue) < aSurfPntDist + Precision::SquareConfusion()) // Point better than initial.
401 {
402 gp_Pnt2d pnt(u - decalU*uperiod,v - decalV*vperiod);
403 return pnt;
404 }
405 }
406
407 // Both searches return bad values, use point from initial 2dcurve.
408 return p2d;
409}
410
411
412//=======================================================================
413//function : ProjLib_PolarFunction
414//purpose : (OCC217 - apo)- This class produce interface to call "gp_Pnt2d Function_Value(...)"
415//=======================================================================
416
417class ProjLib_PolarFunction : public AppCont_Function
418{
419 aFuncStruct myStruct;
420
421 public :
422
423 ProjLib_PolarFunction(const Handle(Adaptor3d_HCurve) & C,
424 const Handle(Adaptor3d_HSurface)& Surf,
425 const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
426 const Standard_Real Tol3d)
427 {
428 myNbPnt = 0;
429 myNbPnt2d = 1;
430
431 computePeriodicity(Surf, myStruct.myPeriod[0], myStruct.myPeriod[1]);
432
433 myStruct.myCurve = C;
434 myStruct.myInitCurve2d = InitialCurve2d;
435 myStruct.mySurf = Surf;
436 myStruct.mySqProjOrtTol = 10000.0 * Tol3d * Tol3d;
437 myStruct.myTolU = Surf->UResolution(Tol3d);
438 myStruct.myTolV = Surf->VResolution(Tol3d);
439 }
440
441 ~ProjLib_PolarFunction() {}
442
443 Standard_Real FirstParameter() const
444 {
445 return myStruct.myCurve->FirstParameter();
446 }
447
448 Standard_Real LastParameter() const
449 {
450 return myStruct.myCurve->LastParameter();
451 }
452
453 gp_Pnt2d Value(const Standard_Real t) const
454 {
455 return Function_Value(t, myStruct);
456 }
457
458 Standard_Boolean Value(const Standard_Real theT,
459 NCollection_Array1<gp_Pnt2d>& thePnt2d,
460 NCollection_Array1<gp_Pnt>& /*thePnt*/) const
461 {
462 thePnt2d(1) = Function_Value(theT, myStruct);
463 return Standard_True;
464 }
465
466 Standard_Boolean D1(const Standard_Real /*theT*/,
467 NCollection_Array1<gp_Vec2d>& /*theVec2d*/,
468 NCollection_Array1<gp_Vec>& /*theVec*/) const
469 {return Standard_False;}
470};
471
472//=======================================================================
473//function : ProjLib_ComputeApproxOnPolarSurface
474//purpose :
475//=======================================================================
476
477ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface()
478: myProjIsDone(Standard_False),
479 myTolerance (-1.0)
480{
481}
482
483
484//=======================================================================
485//function : ProjLib_ComputeApproxOnPolarSurface
486//purpose :
487//=======================================================================
488
489ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
490 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
491 const Handle(Adaptor3d_HCurve)& theCurve,
492 const Handle(Adaptor3d_HSurface)& theSurface,
493 const Standard_Real theTolerance3D)
494: myProjIsDone(Standard_False),
495 myTolerance (theTolerance3D)
496{
497 myBSpline = Perform(theInitialCurve2d, theCurve, theSurface);
498}
499//=======================================================================
500//function : ProjLib_ComputeApproxOnPolarSurface
501//purpose : Process the case of sewing
502//=======================================================================
503
504ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
505 (const Handle(Adaptor2d_HCurve2d)& theInitialCurve2d,
506 const Handle(Adaptor2d_HCurve2d)& theInitialCurve2dBis,
507 const Handle(Adaptor3d_HCurve)& theCurve,
508 const Handle(Adaptor3d_HSurface)& theSurface,
509 const Standard_Real theTolerance3D)
510: myProjIsDone(Standard_False),
511 myTolerance (theTolerance3D)
512{
513 // InitialCurve2d and InitialCurve2dBis are two pcurves of the sewing
514 Handle(Geom2d_BSplineCurve) bsc =
515 Perform(theInitialCurve2d, theCurve, theSurface);
516
517 if(myProjIsDone) {
518 gp_Pnt2d P2dproj, P2d, P2dBis;
519 P2dproj = bsc->StartPoint();
520 P2d = theInitialCurve2d->Value(theInitialCurve2d->FirstParameter());
521 P2dBis = theInitialCurve2dBis->Value(theInitialCurve2dBis->FirstParameter());
522
523 Standard_Real Dist, DistBis;
524 Dist = P2dproj.Distance(P2d);
525 DistBis = P2dproj.Distance(P2dBis);
526 if( Dist < DistBis) {
527 // myBSpline2d is the pcurve that is found. It is translated to obtain myCurve2d
528 myBSpline = bsc;
529 Handle(Geom2d_Geometry) GG = myBSpline->Translated(P2d, P2dBis);
530 my2ndCurve = Handle(Geom2d_Curve)::DownCast(GG);
531 }
532 else {
533 my2ndCurve = bsc;
534 Handle(Geom2d_Geometry) GG = my2ndCurve->Translated(P2dBis, P2d);
535 myBSpline = Handle(Geom2d_BSplineCurve)::DownCast(GG);
536 }
537 }
538}
539
540//=======================================================================
541//function : ProjLib_ComputeApproxOnPolarSurface
542//purpose : case without curve of initialization
543//=======================================================================
544
545ProjLib_ComputeApproxOnPolarSurface::ProjLib_ComputeApproxOnPolarSurface
546 (const Handle(Adaptor3d_HCurve)& theCurve,
547 const Handle(Adaptor3d_HSurface)& theSurface,
548 const Standard_Real theTolerance3D)
549: myProjIsDone(Standard_False),
550 myTolerance (theTolerance3D)
551{
552 const Handle(Adaptor2d_HCurve2d) anInitCurve2d;
553 myBSpline = Perform(anInitCurve2d, theCurve, theSurface);
554}
555
556//=======================================================================
557//function : Concat
558//purpose :
559//=======================================================================
560
561static Handle(Geom2d_BSplineCurve) Concat(Handle(Geom2d_BSplineCurve) C1,
562 Handle(Geom2d_BSplineCurve) C2,
563 Standard_Real theUJump,
564 Standard_Real theVJump)
565{
566 Standard_Integer deg, deg1, deg2;
567 deg1 = C1->Degree();
568 deg2 = C2->Degree();
569
570 if ( deg1 < deg2) {
571 C1->IncreaseDegree(deg2);
572 deg = deg2;
573 }
574 else if ( deg2 < deg1) {
575 C2->IncreaseDegree(deg1);
576 deg = deg1;
577 }
578 else deg = deg1;
579
580 Standard_Integer np1,np2,nk1,nk2,np,nk;
581 np1 = C1->NbPoles();
582 nk1 = C1->NbKnots();
583 np2 = C2->NbPoles();
584 nk2 = C2->NbKnots();
585 nk = nk1 + nk2 -1;
586 np = np1 + np2 -1;
587
588 TColStd_Array1OfReal K1(1,nk1); C1->Knots(K1);
589 TColStd_Array1OfInteger M1(1,nk1); C1->Multiplicities(M1);
590 TColgp_Array1OfPnt2d P1(1,np1); C1->Poles(P1);
591 TColStd_Array1OfReal K2(1,nk2); C2->Knots(K2);
592 TColStd_Array1OfInteger M2(1,nk2); C2->Multiplicities(M2);
593 TColgp_Array1OfPnt2d P2(1,np2); C2->Poles(P2);
594
595 // Compute the new BSplineCurve
596 TColStd_Array1OfReal K(1,nk);
597 TColStd_Array1OfInteger M(1,nk);
598 TColgp_Array1OfPnt2d P(1,np);
599
600 Standard_Integer i, count = 0;
601 // Set Knots and Mults
602 for ( i = 1; i <= nk1; i++) {
603 count++;
604 K(count) = K1(i);
605 M(count) = M1(i);
606 }
607 M(count) = deg;
608 for ( i = 2; i <= nk2; i++) {
609 count++;
610 K(count) = K2(i);
611 M(count) = M2(i);
612 }
613 // Set the Poles
614 count = 0;
615 for (i = 1; i <= np1; i++) {
616 count++;
617 P(count) = P1(i);
618 }
619 for (i = 2; i <= np2; i++) {
620 count++;
621 P(count).SetX(P2(i).X() + theUJump);
622 P(count).SetY(P2(i).Y() + theVJump);
623 }
624
625 Handle(Geom2d_BSplineCurve) BS =
626 new Geom2d_BSplineCurve(P,K,M,deg);
627 return BS;
628}
629
630
631//=======================================================================
632//function : Perform
633//purpose :
634//=======================================================================
635Handle(Geom2d_BSplineCurve) ProjLib_ComputeApproxOnPolarSurface::Perform
636(const Handle(Adaptor2d_HCurve2d)& InitialCurve2d,
637 const Handle(Adaptor3d_HCurve)& Curve,
638 const Handle(Adaptor3d_HSurface)& S)
639{
640 //OCC217
641 Standard_Real Tol3d = myTolerance;
642 Standard_Real ParamTol = Precision::PApproximation();
643
644 Handle(Adaptor2d_HCurve2d) AHC2d = InitialCurve2d;
645 Handle(Adaptor3d_HCurve) AHC = Curve;
646
647// if the curve 3d is a BSpline with degree C0, it is cut into sections with degree C1
648// -> bug cts18237
649 GeomAbs_CurveType typeCurve = Curve->GetType();
650 if(typeCurve == GeomAbs_BSplineCurve) {
651 TColStd_ListOfTransient LOfBSpline2d;
652 Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
653 Standard_Integer nbInter = Curve->NbIntervals(GeomAbs_C1);
654 if(nbInter > 1) {
655 Standard_Integer i, j;
656 Handle(Geom_TrimmedCurve) GTC;
657 Handle(Geom2d_TrimmedCurve) G2dTC;
658 TColStd_Array1OfReal Inter(1,nbInter+1);
659 Curve->Intervals(Inter,GeomAbs_C1);
660 Standard_Real firstinter = Inter.Value(1), secondinter = Inter.Value(2);
661 // initialization 3d
662 GTC = new Geom_TrimmedCurve(BSC, firstinter, secondinter);
663 AHC = new GeomAdaptor_HCurve(GTC);
664
665 // if there is an initialization curve:
666 // - either this is a BSpline C0, with discontinuity at the same parameters of nodes
667 // and the sections C1 are taken
668 // - or this is a curve C1 and the sections of intrest are taken otherwise the curve is created.
669
670 // initialization 2d
671 Standard_Integer nbInter2d;
672 Standard_Boolean C2dIsToCompute;
673 C2dIsToCompute = InitialCurve2d.IsNull();
674 Handle(Geom2d_BSplineCurve) BSC2d;
675 Handle(Geom2d_Curve) G2dC;
676
677 if(!C2dIsToCompute) {
678 nbInter2d = InitialCurve2d->NbIntervals(GeomAbs_C1);
679 TColStd_Array1OfReal Inter2d(1,nbInter2d+1);
680 InitialCurve2d->Intervals(Inter2d,GeomAbs_C1);
681 j = 1;
682 for(i = 1,j = 1;i <= nbInter;i++)
683 if(Abs(Inter.Value(i) - Inter2d.Value(j)) < ParamTol) { //OCC217
684 //if(Abs(Inter.Value(i) - Inter2d.Value(j)) < myTolerance) {
685 if (j > nbInter2d) break;
686 j++;
687 }
688 if(j != (nbInter2d+1)) {
689 C2dIsToCompute = Standard_True;
690 }
691 }
692
693 if(C2dIsToCompute) {
694 AHC2d = BuildInitialCurve2d(AHC, S);
695 }
696 else {
697 typeCurve = InitialCurve2d->GetType();
698 switch (typeCurve) {
699 case GeomAbs_Line: {
700 G2dC = new Geom2d_Line(InitialCurve2d->Line());
701 break;
702 }
703 case GeomAbs_Circle: {
704 G2dC = new Geom2d_Circle(InitialCurve2d->Circle());
705 break;
706 }
707 case GeomAbs_Ellipse: {
708 G2dC = new Geom2d_Ellipse(InitialCurve2d->Ellipse());
709 break;
710 }
711 case GeomAbs_Hyperbola: {
712 G2dC = new Geom2d_Hyperbola(InitialCurve2d->Hyperbola());
713 break;
714 }
715 case GeomAbs_Parabola: {
716 G2dC = new Geom2d_Parabola(InitialCurve2d->Parabola());
717 break;
718 }
719 case GeomAbs_BezierCurve: {
720 G2dC = InitialCurve2d->Bezier();
721 break;
722 }
723 case GeomAbs_BSplineCurve: {
724 G2dC = InitialCurve2d->BSpline();
725 break;
726 }
727 case GeomAbs_OtherCurve:
728 default:
729 break;
730 }
731 gp_Pnt2d fp2d = G2dC->Value(firstinter), lp2d = G2dC->Value(secondinter);
732 gp_Pnt fps, lps, fpc, lpc;
733 S->D0(fp2d.X(), fp2d.Y(), fps);
734 S->D0(lp2d.X(), lp2d.Y(), lps);
735 Curve->D0(firstinter, fpc);
736 Curve->D0(secondinter, lpc);
737 //OCC217
738 if((fps.IsEqual(fpc, Tol3d)) &&
739 (lps.IsEqual(lpc, Tol3d))) {
740 //if((fps.IsEqual(fpc, myTolerance)) &&
741 // (lps.IsEqual(lpc, myTolerance))) {
742 G2dTC = new Geom2d_TrimmedCurve(G2dC, firstinter, secondinter);
743 Geom2dAdaptor_Curve G2dAC(G2dTC);
744 AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
745 myProjIsDone = Standard_True;
746 }
747 else {
748 AHC2d = BuildInitialCurve2d(AHC, S);
749 C2dIsToCompute = Standard_True;
750 }
751 }
752
753 if(myProjIsDone) {
754 BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
755 if(BSC2d.IsNull()) return Handle(Geom2d_BSplineCurve)(); //IFV
756 LOfBSpline2d.Append(BSC2d);
757 }
758 else {
759 return Handle(Geom2d_BSplineCurve)();
760 }
761
762
763
764 Standard_Real iinter, ip1inter;
765 Standard_Integer nbK2d, deg;
766 nbK2d = BSC2d->NbKnots(); deg = BSC2d->Degree();
767
768 for(i = 2;i <= nbInter;i++) {
769 iinter = Inter.Value(i);
770 ip1inter = Inter.Value(i+1);
771 // general case 3d
772 GTC->SetTrim(iinter, ip1inter);
773 AHC = new GeomAdaptor_HCurve(GTC);
774
775 // general case 2d
776 if(C2dIsToCompute) {
777 AHC2d = BuildInitialCurve2d(AHC, S);
778 }
779 else {
780 gp_Pnt2d fp2d = G2dC->Value(iinter), lp2d = G2dC->Value(ip1inter);
781 gp_Pnt fps, lps, fpc, lpc;
782 S->D0(fp2d.X(), fp2d.Y(), fps);
783 S->D0(lp2d.X(), lp2d.Y(), lps);
784 Curve->D0(iinter, fpc);
785 Curve->D0(ip1inter, lpc);
786 //OCC217
787 if((fps.IsEqual(fpc, Tol3d)) &&
788 (lps.IsEqual(lpc, Tol3d))) {
789 //if((fps.IsEqual(fpc, myTolerance)) &&
790 // (lps.IsEqual(lpc, myTolerance))) {
791 G2dTC->SetTrim(iinter, ip1inter);
792 Geom2dAdaptor_Curve G2dAC(G2dTC);
793 AHC2d = new Geom2dAdaptor_HCurve(G2dAC);
794 myProjIsDone = Standard_True;
795 }
796 else {
797 AHC2d = BuildInitialCurve2d(AHC, S);
798 }
799 }
800 if(myProjIsDone) {
801 BSC2d = ProjectUsingInitialCurve2d(AHC, S, AHC2d);
802 if(BSC2d.IsNull()) {
803 return Handle(Geom2d_BSplineCurve)();
804 }
805 LOfBSpline2d.Append(BSC2d);
806 nbK2d += BSC2d->NbKnots() - 1;
807 deg = Max(deg, BSC2d->Degree());
808 }
809 else {
810 return Handle(Geom2d_BSplineCurve)();
811 }
812 }
813
814 Standard_Real anUPeriod, anVPeriod;
815 computePeriodicity(S, anUPeriod, anVPeriod);
816 Standard_Integer NbC = LOfBSpline2d.Extent();
817 Handle(Geom2d_BSplineCurve) CurBS;
818 CurBS = Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
819 LOfBSpline2d.RemoveFirst();
820 for (Standard_Integer ii = 2; ii <= NbC; ii++)
821 {
822 Handle(Geom2d_BSplineCurve) BS =
823 Handle(Geom2d_BSplineCurve)::DownCast(LOfBSpline2d.First());
824
825 //Check for period jump in point of contact.
826 gp_Pnt2d aC1End = CurBS->Pole(CurBS->NbPoles()); // End of C1.
827 gp_Pnt2d aC2Beg = BS->Pole(1); // Beginning of C2.
828 Standard_Real anUJump = 0.0, anVJump = 0.0;
829
830 if (anUPeriod > 0.0 &&
831 Abs (aC1End.X() - aC2Beg.X()) > (anUPeriod ) / 2.01)
832 {
833 Standard_Real aMultCoeff = aC2Beg.X() < aC1End.X() ? 1.0 : -1.0;
834 anUJump = (anUPeriod) * aMultCoeff;
835 }
836
837 if (anVPeriod &&
838 Abs (aC1End.Y() - aC2Beg.Y()) > (anVPeriod) / 2.01)
839 {
840 Standard_Real aMultCoeff = aC2Beg.Y() < aC1End.Y() ? 1.0 : -1.0;
841 anVJump = (anVPeriod) * aMultCoeff;
842 }
843
844 CurBS = Concat(CurBS,BS, anUJump, anVJump);
845 LOfBSpline2d.RemoveFirst();
846 }
847 return CurBS;
848 }
849 }
850
851 if(InitialCurve2d.IsNull()) {
852 AHC2d = BuildInitialCurve2d(Curve, S);
853 if(!myProjIsDone)
854 return Handle(Geom2d_BSplineCurve)();
855 }
856 return ProjectUsingInitialCurve2d(AHC, S, AHC2d);
857}
858
859//=======================================================================
860//function : ProjLib_BuildInitialCurve2d
861//purpose :
862//=======================================================================
863
864Handle(Adaptor2d_HCurve2d)
865 ProjLib_ComputeApproxOnPolarSurface::
866 BuildInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
867 const Handle(Adaptor3d_HSurface)& Surf)
868{
869 // discretize the Curve with quasiuniform deflection
870 // density at least NbOfPnts points
871 myProjIsDone = Standard_False;
872
873 //OCC217
874 Standard_Real Tol3d = myTolerance;
875 Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
876 Standard_Real DistTol3d = 100.0*Tol3d;
877
878 Standard_Real uperiod = 0.0, vperiod = 0.0;
879 computePeriodicity(Surf, uperiod, vperiod);
880
881 // NO myTol is Tol2d !!!!
882 //Standard_Real TolU = myTolerance, TolV = myTolerance;
883 //Standard_Real Tol3d = 100*myTolerance; // At random Balthazar.
884
885 Standard_Integer NbOfPnts = 61;
886 GCPnts_QuasiUniformAbscissa QUA(Curve->GetCurve(),NbOfPnts);
887 TColgp_Array1OfPnt Pts(1,NbOfPnts);
888 TColStd_Array1OfReal Param(1,NbOfPnts);
889 Standard_Integer i, j;
890 for( i = 1; i <= NbOfPnts ; i++ ) {
891 Param(i) = QUA.Parameter(i);
892 Pts(i) = Curve->Value(Param(i));
893 }
894
895 TColgp_Array1OfPnt2d Pts2d(1,NbOfPnts);
896 TColStd_Array1OfInteger Mult(1,NbOfPnts);
897 Mult.Init(1);
898 Mult(1) = Mult(NbOfPnts) = 2;
899
900 Standard_Real Uinf, Usup, Vinf, Vsup;
901 Uinf = Surf->Surface().FirstUParameter();
902 Usup = Surf->Surface().LastUParameter();
903 Vinf = Surf->Surface().FirstVParameter();
904 Vsup = Surf->Surface().LastVParameter();
905 GeomAbs_SurfaceType Type = Surf->GetType();
906 if((Type != GeomAbs_BSplineSurface) && (Type != GeomAbs_BezierSurface) &&
907 (Type != GeomAbs_OffsetSurface)) {
908 Standard_Real S, T;
909// Standard_Integer usens = 0, vsens = 0;
910 // to know the position relatively to the period
911 switch (Type) {
912// case GeomAbs_Plane:
913// {
914// gp_Pln Plane = Surf->Plane();
915// for ( i = 1 ; i <= NbOfPnts ; i++) {
916// ElSLib::Parameters( Plane, Pts(i), S, T);
917// Pts2d(i).SetCoord(S,T);
918// }
919// myProjIsDone = Standard_True;
920// break;
921// }
922 case GeomAbs_Cylinder:
923 {
924// Standard_Real Sloc, Tloc;
925 Standard_Real Sloc;
926 Standard_Integer usens = 0;
927 gp_Cylinder Cylinder = Surf->Cylinder();
928 ElSLib::Parameters( Cylinder, Pts(1), S, T);
929 Pts2d(1).SetCoord(S,T);
930 for ( i = 2 ; i <= NbOfPnts ; i++) {
931 Sloc = S;
932 ElSLib::Parameters( Cylinder, Pts(i), S, T);
933 if(Abs(Sloc - S) > M_PI) {
934 if(Sloc > S)
935 usens++;
936 else
937 usens--;
938 }
939 Pts2d(i).SetCoord(S+usens*2*M_PI,T);
940 }
941 myProjIsDone = Standard_True;
942 break;
943 }
944 case GeomAbs_Cone:
945 {
946// Standard_Real Sloc, Tloc;
947 Standard_Real Sloc;
948 Standard_Integer usens = 0;
949 gp_Cone Cone = Surf->Cone();
950 ElSLib::Parameters( Cone, Pts(1), S, T);
951 Pts2d(1).SetCoord(S,T);
952 for ( i = 2 ; i <= NbOfPnts ; i++) {
953 Sloc = S;
954 ElSLib::Parameters( Cone, Pts(i), S, T);
955 if(Abs(Sloc - S) > M_PI) {
956 if(Sloc > S)
957 usens++;
958 else
959 usens--;
960 }
961 Pts2d(i).SetCoord(S+usens*2*M_PI,T);
962 }
963 myProjIsDone = Standard_True;
964 break;
965 }
966 case GeomAbs_Sphere:
967 {
968 Standard_Real Sloc, Tloc;
969 Standard_Integer usens = 0, vsens = 0; //usens steps by half-period
970 Standard_Boolean vparit = Standard_False;
971 gp_Sphere Sphere = Surf->Sphere();
972 ElSLib::Parameters( Sphere, Pts(1), S, T);
973 Pts2d(1).SetCoord(S,T);
974 for ( i = 2 ; i <= NbOfPnts ; i++) {
975 Sloc = S;Tloc = T;
976 ElSLib::Parameters( Sphere, Pts(i), S, T);
977 if(1.6*M_PI < Abs(Sloc - S)) {
978 if(Sloc > S)
979 usens += 2;
980 else
981 usens -= 2;
982 }
983 if(1.6*M_PI > Abs(Sloc - S) && Abs(Sloc - S) > 0.4*M_PI) {
984 vparit = !vparit;
985 if(Sloc > S)
986 usens++;
987 else
988 usens--;
989 if(Abs(Tloc - Vsup) < (Vsup - Vinf)/5)
990 vsens++;
991 else
992 vsens--;
993 }
994 if(vparit) {
995 Pts2d(i).SetCoord(S+usens*M_PI,(M_PI - T)*(vsens-1));
996 }
997 else {
998 Pts2d(i).SetCoord(S+usens*M_PI,T+vsens*M_PI);
999
1000 }
1001 }
1002 myProjIsDone = Standard_True;
1003 break;
1004 }
1005 case GeomAbs_Torus:
1006 {
1007 Standard_Real Sloc, Tloc;
1008 Standard_Integer usens = 0, vsens = 0;
1009 gp_Torus Torus = Surf->Torus();
1010 ElSLib::Parameters( Torus, Pts(1), S, T);
1011 Pts2d(1).SetCoord(S,T);
1012 for ( i = 2 ; i <= NbOfPnts ; i++) {
1013 Sloc = S; Tloc = T;
1014 ElSLib::Parameters( Torus, Pts(i), S, T);
1015 if(Abs(Sloc - S) > M_PI) {
1016 if(Sloc > S)
1017 usens++;
1018 else
1019 usens--;
1020 }
1021 if(Abs(Tloc - T) > M_PI) {
1022 if(Tloc > T)
1023 vsens++;
1024 else
1025 vsens--;
1026 }
1027 Pts2d(i).SetCoord(S+usens*2*M_PI,T+vsens*2*M_PI);
1028 }
1029 myProjIsDone = Standard_True;
1030 break;
1031 }
1032 default:
1033 Standard_NoSuchObject::Raise("ProjLib_ComputeApproxOnPolarSurface::BuildInitialCurve2d");
1034 }
1035 }
1036 else {
1037 myProjIsDone = Standard_False;
1038 Standard_Real Dist2Min = 1.e+200, u = 0., v = 0.;
1039 gp_Pnt pntproj;
1040
1041 TColgp_SequenceOfPnt2d Sols;
1042 Standard_Boolean areManyZeros = Standard_False;
1043
1044 Curve->D0(Param.Value(1), pntproj) ;
1045 Extrema_ExtPS aExtPS(pntproj, Surf->Surface(), TolU, TolV) ;
1046 Standard_Real aMinSqDist = RealLast();
1047 if (aExtPS.IsDone())
1048 {
1049 for (i = 1; i <= aExtPS.NbExt(); i++)
1050 {
1051 Standard_Real aSqDist = aExtPS.SquareDistance(i);
1052 if (aSqDist < aMinSqDist)
1053 aMinSqDist = aSqDist;
1054 }
1055 }
1056 if (aMinSqDist > DistTol3d * DistTol3d) //try to project with less tolerance
1057 {
1058 TolU = Min(TolU, Precision::PConfusion());
1059 TolV = Min(TolV, Precision::PConfusion());
1060 aExtPS.Initialize(Surf->Surface(),
1061 Surf->Surface().FirstUParameter(), Surf->Surface().LastUParameter(),
1062 Surf->Surface().FirstVParameter(), Surf->Surface().LastVParameter(),
1063 TolU, TolV);
1064 aExtPS.Perform(pntproj);
1065 }
1066
1067 if( aExtPS.IsDone() && aExtPS.NbExt() >= 1 ) {
1068
1069 Standard_Integer GoodValue = 1;
1070
1071 for ( i = 1 ; i <= aExtPS.NbExt() ; i++ ) {
1072 if( aExtPS.SquareDistance(i) < DistTol3d * DistTol3d ) {
1073 if( aExtPS.SquareDistance(i) <= 1.e-18 ) {
1074 aExtPS.Point(i).Parameter(u,v);
1075 gp_Pnt2d p2d(u,v);
1076 Standard_Boolean isSame = Standard_False;
1077 for( j = 1; j <= Sols.Length(); j++ ) {
1078 if( p2d.SquareDistance( Sols.Value(j) ) <= 1.e-18 ) {
1079 isSame = Standard_True;
1080 break;
1081 }
1082 }
1083 if( !isSame ) Sols.Append( p2d );
1084 }
1085 if( Dist2Min > aExtPS.SquareDistance(i) ) {
1086 Dist2Min = aExtPS.SquareDistance(i);
1087 GoodValue = i;
1088 }
1089 }
1090 }
1091
1092 if( Sols.Length() > 1 ) areManyZeros = Standard_True;
1093
1094 if( Dist2Min <= DistTol3d * DistTol3d) {
1095 if( !areManyZeros ) {
1096 aExtPS.Point(GoodValue).Parameter(u,v);
1097 Pts2d(1).SetCoord(u,v);
1098 myProjIsDone = Standard_True;
1099 }
1100 else {
1101 Standard_Integer nbSols = Sols.Length();
1102 Standard_Real Dist2Max = -1.e+200;
1103 for( i = 1; i <= nbSols; i++ ) {
1104 const gp_Pnt2d& aP1 = Sols.Value(i);
1105 for( j = i+1; j <= nbSols; j++ ) {
1106 const gp_Pnt2d& aP2 = Sols.Value(j);
1107 Standard_Real aDist2 = aP1.SquareDistance(aP2);
1108 if( aDist2 > Dist2Max ) Dist2Max = aDist2;
1109 }
1110 }
1111 Standard_Real aMaxT2 = Max(TolU,TolV);
1112 aMaxT2 *= aMaxT2;
1113 if( Dist2Max > aMaxT2 ) {
1114 Standard_Integer tPp = 0;
1115 for( i = 1; i <= 5; i++ ) {
1116 Standard_Integer nbExtOk = 0;
1117 Standard_Integer indExt = 0;
1118 Standard_Integer iT = 1 + (NbOfPnts - 1)/5*i;
1119 Curve->D0( Param.Value(iT), pntproj );
1120 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1121 Dist2Min = 1.e+200;
1122 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1123 for( j = 1 ; j <= aTPS.NbExt() ; j++ ) {
1124 if( aTPS.SquareDistance(j) < DistTol3d * DistTol3d ) {
1125 nbExtOk++;
1126 if( aTPS.SquareDistance(j) < Dist2Min ) {
1127 Dist2Min = aTPS.SquareDistance(j);
1128 indExt = j;
1129 }
1130 }
1131 }
1132 }
1133 if( nbExtOk == 1 ) {
1134 tPp = iT;
1135 aTPS.Point(indExt).Parameter(u,v);
1136 break;
1137 }
1138 }
1139
1140 if( tPp != 0 ) {
1141 gp_Pnt2d aPp = gp_Pnt2d(u,v);
1142 gp_Pnt2d aPn;
1143 j = 1;
1144 Standard_Boolean isFound = Standard_False;
1145 while( !isFound ) {
1146 Curve->D0( Param.Value(tPp+j), pntproj );
1147 Extrema_ExtPS aTPS( pntproj, Surf->Surface(), TolU, TolV );
1148 Dist2Min = 1.e+200;
1149 Standard_Integer indExt = 0;
1150 if( aTPS.IsDone() && aTPS.NbExt() >= 1 ) {
1151 for( i = 1 ; i <= aTPS.NbExt() ; i++ ) {
1152 if( aTPS.SquareDistance(i) < DistTol3d * DistTol3d && aTPS.SquareDistance(i) < Dist2Min ) {
1153 Dist2Min = aTPS.SquareDistance(i);
1154 indExt = i;
1155 isFound = Standard_True;
1156 }
1157 }
1158 }
1159 if( isFound ) {
1160 aTPS.Point(indExt).Parameter(u,v);
1161 aPn = gp_Pnt2d(u,v);
1162 break;
1163 }
1164 j++;
1165 if( (tPp+j) > NbOfPnts ) break;
1166 }
1167
1168 if( isFound ) {
1169 gp_Vec2d atV(aPp,aPn);
1170 Standard_Boolean isChosen = Standard_False;
1171 for( i = 1; i <= nbSols; i++ ) {
1172 const gp_Pnt2d& aP1 = Sols.Value(i);
1173 gp_Vec2d asV(aP1,aPp);
1174 if( asV.Dot(atV) > 0. ) {
1175 isChosen = Standard_True;
1176 Pts2d(1).SetCoord(aP1.X(),aP1.Y());
1177 myProjIsDone = Standard_True;
1178 break;
1179 }
1180 }
1181 if( !isChosen ) {
1182 aExtPS.Point(GoodValue).Parameter(u,v);
1183 Pts2d(1).SetCoord(u,v);
1184 myProjIsDone = Standard_True;
1185 }
1186 }
1187 else {
1188 aExtPS.Point(GoodValue).Parameter(u,v);
1189 Pts2d(1).SetCoord(u,v);
1190 myProjIsDone = Standard_True;
1191 }
1192 }
1193 else {
1194 aExtPS.Point(GoodValue).Parameter(u,v);
1195 Pts2d(1).SetCoord(u,v);
1196 myProjIsDone = Standard_True;
1197 }
1198 }
1199 else {
1200 aExtPS.Point(GoodValue).Parameter(u,v);
1201 Pts2d(1).SetCoord(u,v);
1202 myProjIsDone = Standard_True;
1203 }
1204 }
1205 }
1206
1207 // calculate the following points with GenLocate_ExtPS
1208 // (and store the result and each parameter in a sequence)
1209 Standard_Integer usens = 0, vsens = 0;
1210 // to know the position relatively to the period
1211 Standard_Real U0 = u, V0 = v, U1 = u, V1 = v;
1212 // U0 and V0 are the points in the initialized period
1213 // (period with u and v),
1214 // U1 and V1 are the points for construction of poles
1215
1216 for ( i = 2 ; i <= NbOfPnts ; i++)
1217 if(myProjIsDone) {
1218 myProjIsDone = Standard_False;
1219 Dist2Min = RealLast();
1220 Curve->D0(Param.Value(i), pntproj);
1221 Extrema_GenLocateExtPS aLocateExtPS
1222 (pntproj, Surf->Surface(), U0, V0, TolU, TolV) ;
1223
1224 if (aLocateExtPS.IsDone())
1225 {
1226 if (aLocateExtPS.SquareDistance() < DistTol3d * DistTol3d)
1227 { //OCC217
1228 //if (aLocateExtPS.SquareDistance() < Tol3d * Tol3d) {
1229 (aLocateExtPS.Point()).Parameter(U0,V0);
1230 U1 = U0 + usens*uperiod;
1231 V1 = V0 + vsens*vperiod;
1232 Pts2d(i).SetCoord(U1,V1);
1233 myProjIsDone = Standard_True;
1234 }
1235 else
1236 {
1237 Extrema_ExtPS aGlobalExtr(pntproj, Surf->Surface(), TolU, TolV);
1238 if (aGlobalExtr.IsDone())
1239 {
1240 Standard_Real LocalMinSqDist = RealLast();
1241 Standard_Integer imin = 0;
1242 for (Standard_Integer isol = 1; isol <= aGlobalExtr.NbExt(); isol++)
1243 {
1244 Standard_Real aSqDist = aGlobalExtr.SquareDistance(isol);
1245 if (aSqDist < LocalMinSqDist)
1246 {
1247 LocalMinSqDist = aSqDist;
1248 imin = isol;
1249 }
1250 }
1251 if (LocalMinSqDist < DistTol3d * DistTol3d)
1252 {
1253 Standard_Real LocalU, LocalV;
1254 aGlobalExtr.Point(imin).Parameter(LocalU, LocalV);
1255 if (uperiod > 0. && Abs(U0 - LocalU) >= uperiod/2.)
1256 {
1257 if (LocalU > U0)
1258 usens = -1;
1259 else
1260 usens = 1;
1261 }
1262 if (vperiod > 0. && Abs(V0 - LocalV) >= vperiod/2.)
1263 {
1264 if (LocalV > V0)
1265 vsens = -1;
1266 else
1267 vsens = 1;
1268 }
1269 U0 = LocalU; V0 = LocalV;
1270 U1 = U0 + usens*uperiod;
1271 V1 = V0 + vsens*vperiod;
1272 Pts2d(i).SetCoord(U1,V1);
1273 myProjIsDone = Standard_True;
1274
1275 if((i == 2) && (!IsEqual(uperiod, 0.0) || !IsEqual(vperiod, 0.0)))
1276 {//Make 1st point more precise for periodic surfaces
1277 const Standard_Integer aSize = 3;
1278 const gp_Pnt2d aP(Pts2d(2));
1279 Standard_Real aUpar[aSize], aVpar[aSize];
1280 Pts2d(1).Coord(aUpar[1], aVpar[1]);
1281 aUpar[0] = aUpar[1] - uperiod;
1282 aUpar[2] = aUpar[1] + uperiod;
1283 aVpar[0] = aVpar[1] - vperiod;
1284 aVpar[2] = aVpar[1] + vperiod;
1285
1286 Standard_Real aSQdistMin = RealLast();
1287 Standard_Integer aBestUInd = 1, aBestVInd = 1;
1288 const Standard_Integer aSizeU = IsEqual(uperiod, 0.0) ? 1 : aSize,
1289 aSizeV = IsEqual(vperiod, 0.0) ? 1 : aSize;
1290 for(Standard_Integer uInd = 0; uInd < aSizeU; uInd++)
1291 {
1292 for(Standard_Integer vInd = 0; vInd < aSizeV; vInd++)
1293 {
1294 Standard_Real aSQdist = aP.SquareDistance(gp_Pnt2d(aUpar[uInd], aVpar[vInd]));
1295 if(aSQdist < aSQdistMin)
1296 {
1297 aSQdistMin = aSQdist;
1298 aBestUInd = uInd;
1299 aBestVInd = vInd;
1300 }
1301 }
1302 }
1303
1304 Pts2d(1).SetCoord(aUpar[aBestUInd], aVpar[aBestVInd]);
1305 }//if(i == 2) condition
1306 }
1307 }
1308 }
1309 }
1310 if(!myProjIsDone && uperiod) {
1311 Standard_Real aUinf, aUsup, Uaux;
1312 aUinf = Surf->Surface().FirstUParameter();
1313 aUsup = Surf->Surface().LastUParameter();
1314 if((aUsup - U0) > (U0 - aUinf))
1315 Uaux = 2*aUinf - U0 + uperiod;
1316 else
1317 Uaux = 2*aUsup - U0 - uperiod;
1318 Extrema_GenLocateExtPS locext(pntproj,
1319 Surf->Surface(),
1320 Uaux, V0, TolU, TolV);
1321 if (locext.IsDone())
1322 if (locext.SquareDistance() < DistTol3d * DistTol3d) { //OCC217
1323 //if (locext.SquareDistance() < Tol3d * Tol3d) {
1324 (locext.Point()).Parameter(u,v);
1325 if((aUsup - U0) > (U0 - aUinf))
1326 usens--;
1327 else
1328 usens++;
1329 U0 = u; V0 = v;
1330 U1 = U0 + usens*uperiod;
1331 V1 = V0 + vsens*vperiod;
1332 Pts2d(i).SetCoord(U1,V1);
1333 myProjIsDone = Standard_True;
1334 }
1335 }
1336 if(!myProjIsDone && vperiod) {
1337 Standard_Real aVinf, aVsup, Vaux;
1338 aVinf = Surf->Surface().FirstVParameter();
1339 aVsup = Surf->Surface().LastVParameter();
1340 if((aVsup - V0) > (V0 - aVinf))
1341 Vaux = 2*aVinf - V0 + vperiod;
1342 else
1343 Vaux = 2*aVsup - V0 - vperiod;
1344 Extrema_GenLocateExtPS locext(pntproj,
1345 Surf->Surface(),
1346 U0, Vaux, TolU, TolV) ;
1347 if (locext.IsDone())
1348 if (locext.SquareDistance() < DistTol3d * DistTol3d) { //OCC217
1349 //if (locext.SquareDistance() < Tol3d * Tol3d) {
1350 (locext.Point()).Parameter(u,v);
1351 if((aVsup - V0) > (V0 - aVinf))
1352 vsens--;
1353 else
1354 vsens++;
1355 U0 = u; V0 = v;
1356 U1 = U0 + usens*uperiod;
1357 V1 = V0 + vsens*vperiod;
1358 Pts2d(i).SetCoord(U1,V1);
1359 myProjIsDone = Standard_True;
1360 }
1361 }
1362 if(!myProjIsDone && uperiod && vperiod) {
1363 Standard_Real Uaux, Vaux;
1364 if((Usup - U0) > (U0 - Uinf))
1365 Uaux = 2*Uinf - U0 + uperiod;
1366 else
1367 Uaux = 2*Usup - U0 - uperiod;
1368 if((Vsup - V0) > (V0 - Vinf))
1369 Vaux = 2*Vinf - V0 + vperiod;
1370 else
1371 Vaux = 2*Vsup - V0 - vperiod;
1372 Extrema_GenLocateExtPS locext(pntproj,
1373 Surf->Surface(),
1374 Uaux, Vaux, TolU, TolV);
1375 if (locext.IsDone())
1376 if (locext.SquareDistance() < DistTol3d * DistTol3d) {
1377 //if (locext.SquareDistance() < Tol3d * Tol3d) {
1378 (locext.Point()).Parameter(u,v);
1379 if((Usup - U0) > (U0 - Uinf))
1380 usens--;
1381 else
1382 usens++;
1383 if((Vsup - V0) > (V0 - Vinf))
1384 vsens--;
1385 else
1386 vsens++;
1387 U0 = u; V0 = v;
1388 U1 = U0 + usens*uperiod;
1389 V1 = V0 + vsens*vperiod;
1390 Pts2d(i).SetCoord(U1,V1);
1391 myProjIsDone = Standard_True;
1392 }
1393 }
1394 if(!myProjIsDone) {
1395 Extrema_ExtPS ext(pntproj, Surf->Surface(), TolU, TolV) ;
1396 if (ext.IsDone()) {
1397 Dist2Min = ext.SquareDistance(1);
1398 Standard_Integer aGoodValue = 1;
1399 for ( j = 2 ; j <= ext.NbExt() ; j++ )
1400 if( Dist2Min > ext.SquareDistance(j)) {
1401 Dist2Min = ext.SquareDistance(j);
1402 aGoodValue = j;
1403 }
1404 if (Dist2Min < DistTol3d * DistTol3d) {
1405 //if (Dist2Min < Tol3d * Tol3d) {
1406 (ext.Point(aGoodValue)).Parameter(u,v);
1407 if(uperiod) {
1408 if((U0 - u) > (2*uperiod/3)) {
1409 usens++;
1410 }
1411 else
1412 if((u - U0) > (2*uperiod/3)) {
1413 usens--;
1414 }
1415 }
1416 if(vperiod) {
1417 if((V0 - v) > (vperiod/2)) {
1418 vsens++;
1419 }
1420 else
1421 if((v - V0) > (vperiod/2)) {
1422 vsens--;
1423 }
1424 }
1425 U0 = u; V0 = v;
1426 U1 = U0 + usens*uperiod;
1427 V1 = V0 + vsens*vperiod;
1428 Pts2d(i).SetCoord(U1,V1);
1429 myProjIsDone = Standard_True;
1430 }
1431 }
1432 }
1433 }
1434 else break;
1435 }
1436 }
1437 // -- Pnts2d is transformed into Geom2d_BSplineCurve, with the help of Param and Mult
1438 if(myProjIsDone) {
1439 myBSpline = new Geom2d_BSplineCurve(Pts2d,Param,Mult,1);
1440 //jgv: put the curve into parametric range
1441 gp_Pnt2d MidPoint = myBSpline->Value(0.5*(myBSpline->FirstParameter() + myBSpline->LastParameter()));
1442 Standard_Real TestU = MidPoint.X(), TestV = MidPoint.Y();
1443 Standard_Real sense = 0.;
1444 if (uperiod)
1445 {
1446 if (TestU < Uinf - TolU)
1447 sense = 1.;
1448 else if (TestU > Usup + TolU)
1449 sense = -1;
1450 while (TestU < Uinf - TolU || TestU > Usup + TolU)
1451 TestU += sense * uperiod;
1452 }
1453 if (vperiod)
1454 {
1455 sense = 0.;
1456 if (TestV < Vinf - TolV)
1457 sense = 1.;
1458 else if (TestV > Vsup + TolV)
1459 sense = -1.;
1460 while (TestV < Vinf - TolV || TestV > Vsup + TolV)
1461 TestV += sense * vperiod;
1462 }
1463 gp_Vec2d Offset(TestU - MidPoint.X(), TestV - MidPoint.Y());
1464 if (Abs(Offset.X()) > gp::Resolution() ||
1465 Abs(Offset.Y()) > gp::Resolution())
1466 myBSpline->Translate(Offset);
1467 //////////////////////////////////////////
1468 Geom2dAdaptor_Curve GAC(myBSpline);
1469 Handle(Adaptor2d_HCurve2d) IC2d = new Geom2dAdaptor_HCurve(GAC);
1470#ifdef OCCT_DEBUG
1471// char name [100];
1472// sprintf(name,"%s_%d","build",compteur++);
1473// DrawTrSurf::Set(name,myBSpline);
1474#endif
1475 return IC2d;
1476 }
1477 else {
1478// Modified by Sergey KHROMOV - Thu Apr 18 10:57:50 2002 Begin
1479// Standard_NoSuchObject_Raise_if(1,"ProjLib_Compu: build echec");
1480// Modified by Sergey KHROMOV - Thu Apr 18 10:57:51 2002 End
1481 return Handle(Adaptor2d_HCurve2d)();
1482 }
1483// myProjIsDone = Standard_False;
1484// Modified by Sergey KHROMOV - Thu Apr 18 10:58:01 2002 Begin
1485// Standard_NoSuchObject_Raise_if(1,"ProjLib_ComputeOnPS: build echec");
1486// Modified by Sergey KHROMOV - Thu Apr 18 10:58:02 2002 End
1487}
1488
1489
1490
1491
1492//=======================================================================
1493//function : ProjLib_ProjectUsingInitialCurve2d
1494//purpose :
1495//=======================================================================
1496Handle(Geom2d_BSplineCurve)
1497 ProjLib_ComputeApproxOnPolarSurface::
1498 ProjectUsingInitialCurve2d(const Handle(Adaptor3d_HCurve)& Curve,
1499 const Handle(Adaptor3d_HSurface)& Surf,
1500 const Handle(Adaptor2d_HCurve2d)& InitCurve2d)
1501{
1502 //OCC217
1503 Standard_Real Tol3d = myTolerance;
1504 Standard_Real DistTol3d = 1.0*Tol3d;
1505 Standard_Real TolU = Surf->UResolution(Tol3d), TolV = Surf->VResolution(Tol3d);
1506 Standard_Real Tol2d = Max(Sqrt(TolU*TolU + TolV*TolV), Precision::PConfusion());
1507
1508 Standard_Integer i;
1509 GeomAbs_SurfaceType TheTypeS = Surf->GetType();
1510 GeomAbs_CurveType TheTypeC = Curve->GetType();
1511// Handle(Standard_Type) TheTypeS = Surf->DynamicType();
1512// Handle(Standard_Type) TheTypeC = Curve->DynamicType(); // si on a :
1513// if(TheTypeS == STANDARD_TYPE(Geom_BSplineSurface)) {
1514 if(TheTypeS == GeomAbs_Plane) {
1515 Standard_Real S, T;
1516 gp_Pln Plane = Surf->Plane();
1517 if(TheTypeC == GeomAbs_BSplineCurve) {
1518 Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1519 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1520 for(i = 1;i <= Curve->NbPoles();i++) {
1521 ElSLib::Parameters( Plane, BSC->Pole(i), S, T);
1522 Poles2d(i).SetCoord(S,T);
1523 }
1524 TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1525 BSC->Knots(Knots);
1526 TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1527 BSC->Multiplicities(Mults);
1528 if(BSC->IsRational()) {
1529 TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1530 BSC->Weights(Weights);
1531 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1532 BSC->Degree(), BSC->IsPeriodic()) ;
1533 }
1534 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1535 BSC->Degree(), BSC->IsPeriodic()) ;
1536
1537 }
1538 if(TheTypeC == GeomAbs_BezierCurve) {
1539 Handle(Geom_BezierCurve) BC = Curve->Bezier();
1540 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1541 for(i = 1;i <= Curve->NbPoles();i++) {
1542 ElSLib::Parameters( Plane, BC->Pole(i), S, T);
1543 Poles2d(i).SetCoord(S,T);
1544 }
1545 TColStd_Array1OfReal Knots(1, 2);
1546 Knots.SetValue(1,0.0);
1547 Knots.SetValue(2,1.0);
1548 TColStd_Array1OfInteger Mults(1, 2);
1549 Mults.Init(BC->NbPoles());
1550 if(BC->IsRational()) {
1551 TColStd_Array1OfReal Weights(1, BC->NbPoles());
1552 BC->Weights(Weights);
1553 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1554 BC->Degree(), BC->IsPeriodic()) ;
1555 }
1556 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1557 BC->Degree(), BC->IsPeriodic()) ;
1558 }
1559 }
1560 if(TheTypeS == GeomAbs_BSplineSurface) {
1561 Handle(Geom_BSplineSurface) BSS = Surf->BSpline();
1562 if((BSS->MaxDegree() == 1) &&
1563 (BSS->NbUPoles() == 2) &&
1564 (BSS->NbVPoles() == 2)) {
1565 gp_Pnt p11 = BSS->Pole(1,1);
1566 gp_Pnt p12 = BSS->Pole(1,2);
1567 gp_Pnt p21 = BSS->Pole(2,1);
1568 gp_Pnt p22 = BSS->Pole(2,2);
1569 gp_Vec V1(p11,p12);
1570 gp_Vec V2(p21,p22);
1571 if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ //OCC217
1572 //if(V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1573 // so the polar surface is plane
1574 // and if it is enough to projet the poles of Curve
1575 Standard_Integer Dist2Min = IntegerLast();
1576 Standard_Real u,v;
1577 //OCC217
1578 //Standard_Real TolU = Surf->UResolution(myTolerance)
1579 // , TolV = Surf->VResolution(myTolerance);
1580// gp_Pnt pntproj;
1581 if(TheTypeC == GeomAbs_BSplineCurve) {
1582 Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1583 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1584 for(i = 1;i <= Curve->NbPoles();i++) {
1585 myProjIsDone = Standard_False;
1586 Dist2Min = IntegerLast();
1587 Extrema_GenLocateExtPS extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1588 (p11.Y()+p22.Y())/2,TolU,TolV) ;
1589 if (extrloc.IsDone()) {
1590 Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1591 if (Dist2Min < DistTol3d * DistTol3d) { //OCC217
1592 //if (Dist2Min < myTolerance * myTolerance) {
1593 (extrloc.Point()).Parameter(u,v);
1594 Poles2d(i).SetCoord(u,v);
1595 myProjIsDone = Standard_True;
1596 }
1597 else break;
1598 }
1599 else break;
1600 if(!myProjIsDone)
1601 break;
1602 }
1603 if(myProjIsDone) {
1604 TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1605 BSC->Knots(Knots);
1606 TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1607 BSC->Multiplicities(Mults);
1608 if(BSC->IsRational()) {
1609 TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1610 BSC->Weights(Weights);
1611 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1612 BSC->Degree(), BSC->IsPeriodic()) ;
1613 }
1614 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1615 BSC->Degree(), BSC->IsPeriodic()) ;
1616
1617
1618 }
1619 }
1620 if(TheTypeC == GeomAbs_BezierCurve) {
1621 Handle(Geom_BezierCurve) BC = Curve->Bezier();
1622 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1623 for(i = 1;i <= Curve->NbPoles();i++) {
1624 Dist2Min = IntegerLast();
1625 Extrema_GenLocateExtPS extrloc(BC->Pole(i),Surf->Surface(),0.5,
1626 0.5,TolU,TolV) ;
1627 if (extrloc.IsDone()) {
1628 Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1629 if (Dist2Min < DistTol3d * DistTol3d) { //OCC217
1630 //if (Dist2Min < myTolerance * myTolerance) {
1631 (extrloc.Point()).Parameter(u,v);
1632 Poles2d(i).SetCoord(u,v);
1633 myProjIsDone = Standard_True;
1634 }
1635 else break;
1636 }
1637 else break;
1638 if(myProjIsDone)
1639 myProjIsDone = Standard_False;
1640 else break;
1641 }
1642 if(myProjIsDone) {
1643 TColStd_Array1OfReal Knots(1, 2);
1644 Knots.SetValue(1,0.0);
1645 Knots.SetValue(2,1.0);
1646 TColStd_Array1OfInteger Mults(1, 2);
1647 Mults.Init(BC->NbPoles());
1648 if(BC->IsRational()) {
1649 TColStd_Array1OfReal Weights(1, BC->NbPoles());
1650 BC->Weights(Weights);
1651 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1652 BC->Degree(), BC->IsPeriodic()) ;
1653 }
1654 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1655 BC->Degree(), BC->IsPeriodic()) ;
1656 }
1657 }
1658 }
1659 }
1660 }
1661 else if(TheTypeS == GeomAbs_BezierSurface) {
1662 Handle(Geom_BezierSurface) BS = Surf->Bezier();
1663 if((BS->MaxDegree() == 1) &&
1664 (BS->NbUPoles() == 2) &&
1665 (BS->NbVPoles() == 2)) {
1666 gp_Pnt p11 = BS->Pole(1,1);
1667 gp_Pnt p12 = BS->Pole(1,2);
1668 gp_Pnt p21 = BS->Pole(2,1);
1669 gp_Pnt p22 = BS->Pole(2,2);
1670 gp_Vec V1(p11,p12);
1671 gp_Vec V2(p21,p22);
1672 if(V1.IsEqual(V2,Tol3d,Tol3d/(p11.Distance(p12)*180/M_PI))){ //OCC217
1673 //if (V1.IsEqual(V2,myTolerance,myTolerance/(p11.Distance(p12)*180/M_PI))){
1674 // and if it is enough to project the poles of Curve
1675 Standard_Integer Dist2Min = IntegerLast();
1676 Standard_Real u,v;
1677 //OCC217
1678 //Standard_Real TolU = Surf->UResolution(myTolerance)
1679 // , TolV = Surf->VResolution(myTolerance);
1680
1681// gp_Pnt pntproj;
1682 if(TheTypeC == GeomAbs_BSplineCurve) {
1683 Handle(Geom_BSplineCurve) BSC = Curve->BSpline();
1684 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1685 for(i = 1;i <= Curve->NbPoles();i++) {
1686 myProjIsDone = Standard_False;
1687 Dist2Min = IntegerLast();
1688 Extrema_GenLocateExtPS extrloc(BSC->Pole(i),Surf->Surface(),(p11.X()+p22.X())/2,
1689 (p11.Y()+p22.Y())/2,TolU,TolV) ;
1690 if (extrloc.IsDone()) {
1691 Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1692 if (Dist2Min < DistTol3d * DistTol3d) { //OCC217
1693 //if (Dist2Min < myTolerance * myTolerance) {
1694 (extrloc.Point()).Parameter(u,v);
1695 Poles2d(i).SetCoord(u,v);
1696 myProjIsDone = Standard_True;
1697 }
1698 else break;
1699 }
1700 else break;
1701 if(!myProjIsDone)
1702 break;
1703 }
1704 if(myProjIsDone) {
1705 TColStd_Array1OfReal Knots(1, BSC->NbKnots());
1706 BSC->Knots(Knots);
1707 TColStd_Array1OfInteger Mults(1, BSC->NbKnots());
1708 BSC->Multiplicities(Mults);
1709 if(BSC->IsRational()) {
1710 TColStd_Array1OfReal Weights(1, BSC->NbPoles());
1711 BSC->Weights(Weights);
1712 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1713 BSC->Degree(), BSC->IsPeriodic()) ;
1714 }
1715 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1716 BSC->Degree(), BSC->IsPeriodic()) ;
1717
1718
1719 }
1720 }
1721 if(TheTypeC == GeomAbs_BezierCurve) {
1722 Handle(Geom_BezierCurve) BC = Curve->Bezier();
1723 TColgp_Array1OfPnt2d Poles2d(1,Curve->NbPoles());
1724 for(i = 1;i <= Curve->NbPoles();i++) {
1725 Dist2Min = IntegerLast();
1726 Extrema_GenLocateExtPS extrloc(BC->Pole(i),Surf->Surface(),0.5,
1727 0.5,TolU,TolV) ;
1728 if (extrloc.IsDone()) {
1729 Dist2Min = (Standard_Integer ) extrloc.SquareDistance();
1730 if (Dist2Min < DistTol3d * DistTol3d) { //OCC217
1731 //if (Dist2Min < myTolerance * myTolerance) {
1732 (extrloc.Point()).Parameter(u,v);
1733 Poles2d(i).SetCoord(u,v);
1734 myProjIsDone = Standard_True;
1735 }
1736 else break;
1737 }
1738 else break;
1739 if(myProjIsDone)
1740 myProjIsDone = Standard_False;
1741 else break;
1742 }
1743 if(myProjIsDone) {
1744 TColStd_Array1OfReal Knots(1, 2);
1745 Knots.SetValue(1,0.0);
1746 Knots.SetValue(2,1.0);
1747 TColStd_Array1OfInteger Mults(1, 2);
1748 Mults.Init(BC->NbPoles());
1749 if(BC->IsRational()) {
1750 TColStd_Array1OfReal Weights(1, BC->NbPoles());
1751 BC->Weights(Weights);
1752 return new Geom2d_BSplineCurve(Poles2d, Weights, Knots, Mults,
1753 BC->Degree(), BC->IsPeriodic()) ;
1754 }
1755 return new Geom2d_BSplineCurve(Poles2d, Knots, Mults,
1756 BC->Degree(), BC->IsPeriodic()) ;
1757 }
1758 }
1759 }
1760 }
1761 }
1762
1763 ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, Tol3d) ; //OCC217
1764 //ProjLib_PolarFunction F(Curve, Surf, InitCurve2d, myTolerance) ;
1765
1766#ifdef OCCT_DEBUG
1767 Standard_Integer Nb = 50;
1768
1769 Standard_Real U, U1, U2;
1770 U1 = F.FirstParameter();
1771 U2 = F.LastParameter();
1772
1773 TColgp_Array1OfPnt2d DummyPoles(1,Nb+1);
1774 TColStd_Array1OfReal DummyKnots(1,Nb+1);
1775 TColStd_Array1OfInteger DummyMults(1,Nb+1);
1776 DummyMults.Init(1);
1777 DummyMults(1) = 2;
1778 DummyMults(Nb+1) = 2;
1779 for (Standard_Integer ij = 0; ij <= Nb; ij++) {
1780 U = (Nb-ij)*U1 + ij*U2;
1781 U /= Nb;
1782 DummyPoles(ij+1) = F.Value(U);
1783 DummyKnots(ij+1) = ij;
1784 }
1785 Handle(Geom2d_BSplineCurve) DummyC2d =
1786 new Geom2d_BSplineCurve(DummyPoles, DummyKnots, DummyMults, 1);
1787#ifdef DRAW
1788 Standard_CString Temp = "bs2d";
1789 DrawTrSurf::Set(Temp,DummyC2d);
1790#endif
1791// DrawTrSurf::Set((Standard_CString ) "bs2d",DummyC2d);
1792 Handle(Geom2dAdaptor_HCurve) DDD =
1793 Handle(Geom2dAdaptor_HCurve)::DownCast(InitCurve2d);
1794
1795#ifdef DRAW
1796 Temp = "initc2d";
1797 DrawTrSurf::Set(Temp,DDD->ChangeCurve2d().Curve());
1798#endif
1799// DrawTrSurf::Set((Standard_CString ) "initc2d",DDD->ChangeCurve2d().Curve());
1800#endif
1801
1802 Standard_Integer Deg1,Deg2;
1803// Deg1 = 8;
1804// Deg2 = 8;
1805 Deg1 = 2; //IFV
1806 Deg2 = 8; //IFV
1807
1808 Approx_FitAndDivide2d Fit(F,Deg1,Deg2,Tol3d,Tol2d, //OCC217
1809 //Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
1810 Standard_True);
1811
1812 if(Fit.IsAllApproximated()) {
1813 Standard_Integer j;
1814 Standard_Integer NbCurves = Fit.NbMultiCurves();
1815 Standard_Integer MaxDeg = 0;
1816 // To transform the MultiCurve into BSpline, it is required that all
1817 // Bezier constituing it have the same degree -> Calculation of MaxDeg
1818 Standard_Integer NbPoles = 1;
1819 for (j = 1; j <= NbCurves; j++) {
1820 Standard_Integer Deg = Fit.Value(j).Degree();
1821 MaxDeg = Max ( MaxDeg, Deg);
1822 }
1823
1824 NbPoles = MaxDeg * NbCurves + 1; //Tops on the BSpline
1825 TColgp_Array1OfPnt2d Poles( 1, NbPoles);
1826
1827 TColgp_Array1OfPnt2d TempPoles( 1, MaxDeg + 1);//to augment the degree
1828
1829 TColStd_Array1OfReal Knots( 1, NbCurves + 1); //Nodes of the BSpline
1830
1831 Standard_Integer Compt = 1;
1832 for (i = 1; i <= NbCurves; i++) {
1833 Fit.Parameters(i, Knots(i), Knots(i+1));
1834 AppParCurves_MultiCurve MC = Fit.Value( i); //Load the Ith Curve
1835 TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Retrieve the tops
1836 MC.Curve(1, Poles2d);
1837
1838 //Eventual augmentation of the degree
1839 Standard_Integer Inc = MaxDeg - MC.Degree();
1840 if ( Inc > 0) {
1841// BSplCLib::IncreaseDegree( Inc, Poles2d, PLib::NoWeights(),
1842 BSplCLib::IncreaseDegree( MaxDeg, Poles2d, PLib::NoWeights(),
1843 TempPoles, PLib::NoWeights());
1844 //update of tops of the PCurve
1845 for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1846 Poles.SetValue( Compt, TempPoles( k));
1847 Compt++;
1848 }
1849 }
1850 else {
1851 //update of tops of the PCurve
1852 for (Standard_Integer k = 1 ; k <= MaxDeg + 1; k++) {
1853 Poles.SetValue( Compt, Poles2d( k));
1854 Compt++;
1855 }
1856 }
1857
1858 Compt--;
1859 }
1860
1861 //update of fields of ProjLib_Approx
1862 Standard_Integer NbKnots = NbCurves + 1;
1863
1864 TColStd_Array1OfInteger Mults( 1, NbKnots);
1865 Mults.Init(MaxDeg);
1866 Mults.SetValue( 1, MaxDeg + 1);
1867 Mults.SetValue(NbKnots, MaxDeg + 1);
1868 myProjIsDone = Standard_True;
1869 Handle(Geom2d_BSplineCurve) Dummy =
1870 new Geom2d_BSplineCurve(Poles,Knots,Mults,MaxDeg);
1871
1872 // try to smoother the Curve GeomAbs_C1.
1873
1874 Standard_Boolean OK = Standard_True;
1875
1876 for (Standard_Integer ij = 2; ij < NbKnots; ij++) {
1877 OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,Tol3d); //OCC217
1878 //OK = OK && Dummy->RemoveKnot(ij,MaxDeg-1,myTolerance);
1879 }
1880#ifdef OCCT_DEBUG
1881 if (!OK) {
1882 cout << "ProjLib_ComputeApproxOnPolarSurface : Smoothing echoue"<<endl;
1883 }
1884#endif
1885 return Dummy;
1886 }
1887 return Handle(Geom2d_BSplineCurve)();
1888}
1889
1890//=======================================================================
1891//function : BSpline
1892//purpose :
1893//=======================================================================
1894
1895Handle(Geom2d_BSplineCurve)
1896 ProjLib_ComputeApproxOnPolarSurface::BSpline() const
1897
1898{
1899// Modified by Sergey KHROMOV - Thu Apr 18 11:16:46 2002 End
1900// Standard_NoSuchObject_Raise_if
1901// (!myProjIsDone,
1902// "ProjLib_ComputeApproxOnPolarSurface:BSpline");
1903// Modified by Sergey KHROMOV - Thu Apr 18 11:16:47 2002 End
1904 return myBSpline ;
1905}
1906
1907//=======================================================================
1908//function : Curve2d
1909//purpose :
1910//=======================================================================
1911
1912Handle(Geom2d_Curve)
1913 ProjLib_ComputeApproxOnPolarSurface::Curve2d() const
1914
1915{
1916 Standard_NoSuchObject_Raise_if
1917 (!myProjIsDone,
1918 "ProjLib_ComputeApproxOnPolarSurface:2ndCurve2d");
1919 return my2ndCurve ;
1920}
1921
1922
1923//=======================================================================
1924//function : IsDone
1925//purpose :
1926//=======================================================================
1927
1928Standard_Boolean ProjLib_ComputeApproxOnPolarSurface::IsDone() const
1929
1930{
1931 return myProjIsDone;
1932}
Open CASCADE Technology repository