projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022627: Change OCCT memory management defaults
[occt.git] / src / ProjLib / ProjLib_ComputeApprox.cxx
CommitLineData
7fd59977 1// File: ProjLib_ComputeApprox.gxx
2// Created: Tue Sep 7 16:37:49 1993
3// Author: Bruno DUMORTIER
4// <dub@topsn3>
5// modified by NIZHNY-OFV Thu Jan 20 11:04:19 2005
6
7#include <ProjLib_ComputeApprox.hxx>
8
9#include <GeomAbs_SurfaceType.hxx>
10#include <GeomAbs_CurveType.hxx>
11#include <AppCont_Function2d.hxx>
12#include <Convert_CompBezierCurves2dToBSplineCurve2d.hxx>
13#include <ElSLib.hxx>
14#include <ElCLib.hxx>
15#include <BSplCLib.hxx>
16#include <Standard_NoSuchObject.hxx>
17#include <Geom_UndefinedDerivative.hxx>
18#include <gp.hxx>
19#include <gp_Trsf.hxx>
20#include <Precision.hxx>
21#include <Approx_FitAndDivide2d.hxx>
22#include <AppParCurves_MultiCurve.hxx>
23#include <Handle_Adaptor3d_HCurve.hxx>
24#include <Adaptor3d_HCurve.hxx>
25#include <Handle_Adaptor3d_HSurface.hxx>
26#include <Adaptor3d_HSurface.hxx>
27#include <TColgp_Array1OfPnt2d.hxx>
28#include <TColgp_Array1OfPnt.hxx>
29#include <TColStd_Array1OfReal.hxx>
30#include <TColStd_Array1OfInteger.hxx>
31#include <Geom_BSplineCurve.hxx>
32#include <Geom_BezierCurve.hxx>
33#include <Geom2d_BSplineCurve.hxx>
34#include <Geom2d_BezierCurve.hxx>
35
36#ifdef DRAW
37#include <DrawTrSurf.hxx>
38#endif
39#ifdef DEB
40static Standard_Boolean AffichValue = Standard_False;
41#endif
42
43static
44 void Parameters(const Handle(Adaptor3d_HCurve)& myCurve,
45 const Handle(Adaptor3d_HSurface)& mySurface,
46 const gp_Pnt& aP1,
47 const Standard_Integer iFirst,
48 const Standard_Real aTolU,
49 Standard_Real& aU,
50 Standard_Real& aV);
51
52//=======================================================================
53//function : IsEqual
54//purpose :
55//=======================================================================
56// OFV:
57static inline Standard_Boolean IsEqual(Standard_Real Check,Standard_Real With,Standard_Real Toler)
58{
59 return ((Abs(Check - With) < Toler) ? Standard_True : Standard_False);
60}
61
62
63//=======================================================================
64//function : Value
65//purpose :
66//=======================================================================
67
68static gp_Pnt2d Function_Value(const Standard_Real U,
69 const Handle(Adaptor3d_HCurve)& myCurve,
70 const Handle(Adaptor3d_HSurface)& mySurface,
71 const Standard_Real U1,
72 const Standard_Real U2,
73 const Standard_Real V1,
74 const Standard_Real V2,
75 const Standard_Boolean UCouture,
76 const Standard_Boolean VCouture )
77{
78 Standard_Real S, T;
79
80 gp_Pnt P3d = myCurve->Value(U);
81 GeomAbs_SurfaceType SType = mySurface->GetType();
82
83 switch ( SType ) {
84
85 case GeomAbs_Plane:
86 {
87 gp_Pln Plane = mySurface->Plane();
88 ElSLib::Parameters( Plane, P3d, S, T);
89 break;
90 }
91 case GeomAbs_Cylinder:
92 {
93 gp_Cylinder Cylinder = mySurface->Cylinder();
94 ElSLib::Parameters( Cylinder, P3d, S, T);
95 break;
96 }
97 case GeomAbs_Cone:
98 {
99 gp_Cone Cone = mySurface->Cone();
100 ElSLib::Parameters( Cone, P3d, S, T);
101 break;
102 }
103 case GeomAbs_Sphere:
104 {
105 gp_Sphere Sphere = mySurface->Sphere();
106 ElSLib::Parameters(Sphere, P3d, S, T);
107 break;
108 }
109 case GeomAbs_Torus:
110 {
111 gp_Torus Torus = mySurface->Torus();
112 ElSLib::Parameters( Torus, P3d, S, T);
113 break;
114 }
115 default:
116 Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
117 }
118
119 if ( UCouture) {
120 S = ElCLib::InPeriod(S, U1, U2);
121 }
122
123 if ( VCouture) {
124 if(SType == GeomAbs_Sphere) {
c6541a0c
D		 125 if ( Abs( S - U1 ) > M_PI ) {
126 T = M_PI - T;
127 S = M_PI + S;
7fd59977 128 }
129 S = ElCLib::InPeriod(S, U1, U2);
130 }
131 T = ElCLib::InPeriod(T, V1, V2);
132 }
133
134 return gp_Pnt2d(S, T);
135}
136//=======================================================================
137//function : D1
138//purpose :
139//=======================================================================
140static Standard_Boolean Function_D1( const Standard_Real U,
141 gp_Pnt2d& P,
142 gp_Vec2d& D,
143 const Handle(Adaptor3d_HCurve)& myCurve,
144 const Handle(Adaptor3d_HSurface)& mySurface,
145 const Standard_Real U1,
146 const Standard_Real U2,
147 const Standard_Real V1,
148 const Standard_Real V2,
149 const Standard_Boolean UCouture,
150 const Standard_Boolean VCouture )
151{
152 gp_Pnt P3d;
153 Standard_Real dU, dV;
154
155 P = Function_Value(U,myCurve,mySurface,U1,U2,V1,V2,UCouture,VCouture);
156
157 GeomAbs_SurfaceType Type = mySurface->GetType();
158
159 switch ( Type) {
160 case GeomAbs_Plane:
161 case GeomAbs_Cone:
162 case GeomAbs_Cylinder:
163 case GeomAbs_Sphere:
164 case GeomAbs_Torus:
165 {
166 gp_Vec D1U, D1V;
167 gp_Vec T;
168 myCurve->D1(U,P3d,T);
169 mySurface->D1(P.X(),P.Y(),P3d,D1U,D1V);
170
171 dU = T.Dot(D1U);
172 dV = T.Dot(D1V);
173 Standard_Real Nu = D1U.SquareMagnitude();
174 Standard_Real Nv = D1V.SquareMagnitude();
175
176 if ( Nu < Epsilon(1.) || Nv < Epsilon(1.))
177 return Standard_False;
178
179 dU /= Nu;
180 dV /= Nv;
181 D = gp_Vec2d( dU, dV);
182 }
183 break;
184
185 default:
186 return Standard_False;
187 }
188
189 return Standard_True;
190}
191
192//=======================================================================
193//function : Function_SetUVBounds
194//purpose :
195//=======================================================================
196static void Function_SetUVBounds(Standard_Real& myU1,
197 Standard_Real& myU2,
198 Standard_Real& myV1,
199 Standard_Real& myV2,
200 Standard_Boolean& UCouture,
201 Standard_Boolean& VCouture,
202 const Handle(Adaptor3d_HCurve)& myCurve,
203 const Handle(Adaptor3d_HSurface)& mySurface)
204{
205 Standard_Real W1, W2, W;
206 gp_Pnt P1, P2, P;
207 //
208 W1 = myCurve->FirstParameter();
209 W2 = myCurve->LastParameter ();
210 W = 0.5*(W1+W2);
211 // on ouvre l`intervalle
212 // W1 += 1.0e-9;
213 // W2 -= 1.0e-9;
214 P1 = myCurve->Value(W1);
215 P2 = myCurve->Value(W2);
216 P = myCurve->Value(W);
217
218 switch ( mySurface->GetType()) {
219
220 case GeomAbs_Cone: {
221 gp_Cone Cone = mySurface->Cone();
222 VCouture = Standard_False;
223
224 switch( myCurve->GetType() ){
225 case GeomAbs_Parabola:
226 case GeomAbs_Hyperbola:
227 case GeomAbs_Ellipse:{
228 Standard_Real U1, U2, V1, V2, U , V;
229 ElSLib::Parameters( Cone, P1, U1, V1);
230 ElSLib::Parameters( Cone, P2, U2, V2);
231 ElSLib::Parameters( Cone, P , U , V );
232 myU1 = Min(U1,U2);
233 myU2 = Max(U1,U2);
234 if ( ( U1 < U && U < U2 ) && !myCurve->IsClosed() ) {
235 UCouture = Standard_False;
236 }
237 else {
238 UCouture = Standard_True;
c6541a0c 239 myU2 = myU1 + 2*M_PI;
7fd59977 240 }
241
242 }
243 break;
244 default: {
245 Standard_Real U1, V1, U , V, Delta = 0., d = 0., pmin = W1, pmax = W1, dmax = 0., Uf, Ul;
246 ElSLib::Parameters( Cone, P1, U1, V1);
247 ElSLib::Parameters( Cone, P2, Ul, V1);
248 myU1 = U1; myU2 = U1; Uf = U1;
249 Standard_Real Step = .1;
250 Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
251 nbp = Max(nbp, 3);
252 Step = (W2 - W1) / (nbp - 1);
253 Standard_Boolean isclandper = (!(myCurve->IsClosed()) && !(myCurve->IsPeriodic()));
254 for(Standard_Real par = W1 + Step; par <= W2; par += Step) {
255 if(!isclandper) par += Step;
256 P = myCurve->Value(par);
257 ElSLib::Parameters( Cone, P, U, V);
258 U += Delta;
259 d = U - U1;
c6541a0c
D		 260 if(d > M_PI) {
261 if( ( (IsEqual(U,(2*M_PI),1.e-10) && (U1 >= 0. && U1 <= M_PI)) &&
7fd59977 262 (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,0.,1.e-10)) ) && isclandper ) U = 0.;
c6541a0c 263 else Delta -= 2*M_PI;
7fd59977 264 U += Delta;
265 d = U - U1;
266 }
c6541a0c
D		 267 else if(d < -M_PI) {
268 if( ( (IsEqual(U,0.,1.e-10) && (U1 >= M_PI && U1 <= (2*M_PI))) &&
269 (IsEqual(U,Ul,1.e-10) && !IsEqual(Uf,(2*M_PI),1.e-10)) ) && isclandper ) U = 2*M_PI;
270 else Delta += 2*M_PI;
7fd59977 271 U += Delta;
272 d = U - U1;
273 }
274 dmax = Max(dmax, Abs(d));
275 if(U < myU1) {myU1 = U; pmin = par;}
276 if(U > myU2) {myU2 = U; pmax = par;}
277 U1 = U;
278 }
279
280 if(!(Abs(pmin - W1) <= Precision::PConfusion() || Abs(pmin - W2) <= Precision::PConfusion()) ) myU1 -= dmax*.5;
281 if(!(Abs(pmax - W1) <= Precision::PConfusion() || Abs(pmax - W2) <= Precision::PConfusion()) ) myU2 += dmax*.5;
282
c6541a0c 283 if((myU1 >=0. && myU1 <= 2*M_PI) && (myU2 >=0. && myU2 <= 2*M_PI) ) UCouture = Standard_False;
7fd59977 284 else{
285 U = ( myU1 + myU2 ) /2.;
c6541a0c
D		 286 myU1 = U - M_PI;
287 myU2 = U + M_PI;
7fd59977 288 UCouture = Standard_True;
289 }
290 }
291 break;
292 }// switch curve type
293 }// case Cone
294 break;
295
296 case GeomAbs_Cylinder: {
297 gp_Cylinder Cylinder = mySurface->Cylinder();
298 VCouture = Standard_False;
299
300 if (myCurve->GetType() == GeomAbs_Ellipse) {
301
302 Standard_Real U1, U2, V1, V2, U , V;
303 ElSLib::Parameters( Cylinder, P1, U1, V1);
304 ElSLib::Parameters( Cylinder, P2, U2, V2);
305 ElSLib::Parameters( Cylinder, P , U , V );
306 myU1 = Min(U1,U2);
307 myU2 = Max(U1,U2);
308
309 if ( !myCurve->IsClosed()) {
310 if ( myU1 < U && U < myU2) {
311 U = ( myU1 + myU2 ) /2.;
c6541a0c
D		 312 myU1 = U - M_PI;
313 myU2 = U + M_PI;
7fd59977 314 }
315 else {
316 U = ( myU1 + myU2 ) /2.;
317 if ( myU1 < U) {
c6541a0c 318 myU1 = U - 2*M_PI;
7fd59977 319 myU2 = U;
320 }
321 else {
322 myU1 = U;
c6541a0c 323 myU2 = U + 2*M_PI;
7fd59977 324 }
325 }
326 UCouture = Standard_True;
327 }
328 else {
329 gp_Vec D1U, D1V;
330 gp_Vec T;
331 gp_Pnt P3d;
332 myCurve->D1(W1,P3d,T);
333 mySurface->D1(U1,U2,P3d,D1U,D1V);
334 Standard_Real dU = T.Dot(D1U);
335
336 UCouture = Standard_True;
337 if ( dU > 0.) {
c6541a0c 338 myU2 = myU1 + 2*M_PI;
7fd59977 339 }
340 else {
341 myU2 = myU1;
c6541a0c 342 myU1 -= 2*M_PI;
7fd59977 343 }
344 }
345 }
346 else {
347 Standard_Real U1, V1, U , V;
348 ElSLib::Parameters( Cylinder, P1, U1, V1);
349 Standard_Real Step = .1, Delta = 0.;
c6541a0c 350 Standard_Real eps = M_PI, dmax = 0., d = 0.;
7fd59977 351 Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
352 nbp = Max(nbp, 3);
353 Step = (W2 - W1) / (nbp - 1);
354 myU1 = U1; myU2 = U1;
355 Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
356 for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
357 P = myCurve->Value(par);
358 ElSLib::Parameters( Cylinder, P, U, V);
359 U += Delta;
360 d = U - U1;
361 if(d > eps) {
362 U -= Delta;
c6541a0c 363 Delta -= 2*M_PI;
7fd59977 364 U += Delta;
365 d = U - U1;
366 }
367 else if(d < -eps) {
368 U -= Delta;
c6541a0c 369 Delta += 2*M_PI;
7fd59977 370 U += Delta;
371 d = U - U1;
372 }
373 dmax = Max(dmax, Abs(d));
374 if(U < myU1) {myU1 = U; pmin = par;}
375 if(U > myU2) {myU2 = U; pmax = par;}
376 U1 = U;
377 }
378
379 if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
380 Abs(pmin - W2) <= Precision::PConfusion()) ) myU1 -= dmax*.5;
381 if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
382 Abs(pmax - W2) <= Precision::PConfusion()) ) myU2 += dmax*.5;
383
c6541a0c
D		 384 if((myU1 >=0. && myU1 <= 2*M_PI) &&
385 (myU2 >=0. && myU2 <= 2*M_PI) ) {
7fd59977 386 UCouture = Standard_False;
387 }
388 else {
389 U = ( myU1 + myU2 ) /2.;
c6541a0c
D		 390 myU1 = U - M_PI;
391 myU2 = U + M_PI;
7fd59977 392 UCouture = Standard_True;
393 }
394 }
395 }
396 break;
397 //
398 case GeomAbs_Sphere:{
399 VCouture = Standard_False;
400 gp_Sphere SP = mySurface->Sphere();
401 if ( myCurve->GetType() == GeomAbs_Circle) {
402 UCouture = Standard_True;
403
404 // on cherche a savoir le nombre de fois que la couture est
405 // traversee.
406 // si 0 ou 2 fois : la PCurve est fermee et dans l`intervalle
407 // [Uc-PI, Uc+PI] (Uc: U du centre du cercle)
408 // si 1 fois : la PCurve est ouverte et dans l`intervalle
409 // [U1, U1 +/- 2*PI]
410
411 // pour determiner le nombre de solution, on resoud le systeme
412 // x^2 + y^2 + z^2 = R^2 (1)
413 // A x + B y + C z + D = 0 (2)
414 // x > 0 (3)
415 // y = 0 (4)
416 // REM : (1) (2) : equation du cercle
417 // (1) (3) (4) : equation de la couture.
418 Standard_Integer NbSolutions = 0;
419 Standard_Real A, B, C, D, R, Tol = 1.e-10;
420 Standard_Real U1, U2, V1, V2, aTPC;
421 gp_Trsf Trsf;
422 //
423 aTPC=Precision::PConfusion();
424 //
425 gp_Circ Circle = myCurve->Circle();
426 Trsf.SetTransformation(SP.Position());
427 Circle.Transform(Trsf);
428 //
429 R = SP.Radius();
430 gp_Pln Plane( gp_Ax3(Circle.Position()));
431 Plane.Coefficients(A,B,C,D);
432 //
433 if ( Abs(C) < Tol) {
434 if ( Abs(A) > Tol) {
435 if ( (D/A) < 0.) {
436 if ( ( R - Abs(D/A)) > Tol) NbSolutions = 2;
437 else if ( Abs(R - Abs(D/A))< Tol) NbSolutions = 1;
438 else NbSolutions = 0;
439 }
440 }
441 }
442 else {
443 Standard_Real delta = R*R*(A*A+C*C) - D*D;
444 delta *= C*C;
445 if ( Abs(delta) < Tol*Tol) {
446 if ( A*D > 0.) NbSolutions = 1;
447 }
448 else if ( delta > 0) {
449 Standard_Real xx;
450 delta = Sqrt(delta);
451 xx = -A*D+delta;
452 //
453 if ( xx > Tol) NbSolutions++;
454 xx = -A*D-delta;
455 //
456 if ( xx > Tol) NbSolutions++;
457 }
458 }
459 //
460
461 // box+sphere >>
462 Standard_Real UU = 0.;
463 ElSLib::Parameters(SP, P1, U1, V1);
464 ElSLib::Parameters(SP, P2, U2, V1);
465 ElSLib::Parameters(SP, P, UU, V1);
466 Standard_Real UUmi = Min(Min(U1,UU),Min(UU,U2));
467 Standard_Real UUma = Max(Max(U1,UU),Max(UU,U2));
c6541a0c 468 Standard_Boolean reCalc = ((UUmi >= 0. && UUmi <= M_PI) && (UUma >= 0. && UUma <= M_PI));
7fd59977 469 // box+sphere <<
470
471 ElSLib::Parameters(SP, P1, U1, V1);//*
472 //
473 Parameters(myCurve, mySurface, P1, 1, aTPC, U1, V1);
474 //
475 //
c6541a0c 476 P2 = myCurve->Value(W1+M_PI/8);
7fd59977 477 ElSLib::Parameters(SP,P2,U2,V2);
478 //
479 if ( NbSolutions == 1) {
c6541a0c
D		 480 if ( Abs(U1-U2) > M_PI) { // on traverse la couture
481 if ( U1 > M_PI) {
7fd59977 482 myU1 = U1;
c6541a0c 483 myU2 = U1+2*M_PI;
7fd59977 484 }
485 else {
486 myU2 = U1;
c6541a0c 487 myU1 = U1-2*M_PI;
7fd59977 488 }
489 }
490 else { // on ne traverse pas la couture
491 if ( U1 > U2) {
492 myU2 = U1;
c6541a0c 493 myU1 = U1-2*M_PI;
7fd59977 494 }
495 else {
496 myU1 = U1;
c6541a0c 497 myU2 = U1+2*M_PI;
7fd59977 498 }
499 }
500 }
501 else { // 0 ou 2 solutions
502 gp_Pnt Center = Circle.Location();
503 Standard_Real U,V;
504 ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
c6541a0c
D		 505 myU1 = U-M_PI;
506 myU2 = U+M_PI;
7fd59977 507 }
508 //
509 // eval the VCouture.
510 if ( (C==0) || Abs(Abs(D/C)-R) > 1.e-10) {
511 VCouture = Standard_False;
512 }
513 else {
514 VCouture = Standard_True;
515 UCouture = Standard_True;
516
517 if ( D/C < 0.) {
c6541a0c
D		 518 myV1 = - M_PI / 2.;
519 myV2 = 3 * M_PI / 2.;
7fd59977 520 }
521 else {
c6541a0c
D		 522 myV1 = -3 * M_PI / 2.;
523 myV2 = M_PI / 2.;
7fd59977 524 }
525
526 // si P1.Z() vaut +/- R on est sur le sommet : pas significatif.
527 gp_Pnt pp = P1.Transformed(Trsf);
528
529 if ( Abs( Abs(pp.Z()) - R) < Tol) {
530 gp_Pnt Center = Circle.Location();
531 Standard_Real U,V;
532 ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
c6541a0c
D		 533 myU1 = U-M_PI;
534 myU2 = U+M_PI;
7fd59977 535 VCouture = Standard_False;
536 }
537 else {
538 ElSLib::Parameters(SP,P1,U1,V1);//*
539 //
540 Parameters(myCurve, mySurface, P1, 1, aTPC, U1, V1);
541 //
c6541a0c 542 P2 = myCurve->Value(W1+M_PI/8);
7fd59977 543 ElSLib::Parameters(SP,P2,U2,V2);
544
c6541a0c
D		 545 if ( Abs(U1-U2) > M_PI) { // on traverse la couture
546 if ( U1 > M_PI) {
7fd59977 547 myU1 = U1;
c6541a0c 548 myU2 = U1+2*M_PI;
7fd59977 549 }
550 else {
551 myU2 = U1;
c6541a0c 552 myU1 = U1-2*M_PI;
7fd59977 553 }
554 }
555 else { // on ne traverse pas la couture
556 if ( U1 > U2) {
557 myU2 = U1;
c6541a0c 558 myU1 = U1-2*M_PI;
7fd59977 559 }
560 else {
561 myU1 = U1;
c6541a0c 562 myU2 = U1+2*M_PI;
7fd59977 563 }
564 }
565 }
566 }
567
568 // box+sphere >>
569 myV1 = -1.e+100; myV2 = 1.e+100;
570 Standard_Real UU1 = myU1, UU2 = myU2;
c6541a0c 571 if((Abs(UU1) <= (2.*M_PI) && Abs(UU2) <= (2.*M_PI)) && NbSolutions == 1 && reCalc) {
7fd59977 572 gp_Pnt Center = Circle.Location();
573 Standard_Real U,V;
574 ElSLib::SphereParameters(gp_Ax3(gp::XOY()),1,Center, U, V);
c6541a0c
D		 575 myU1 = U-M_PI;
576 myU2 = U+M_PI;
7fd59977 577 myU1 = Min(UU1,myU1);
578 myU2 = Max(UU2,myU2);
579 }
580 // box+sphere <<
581
582 }//if ( myCurve->GetType() == GeomAbs_Circle)
583
584 else {
585 Standard_Real U1, V1, U , V;
586 ElSLib::Parameters( SP, P1, U1, V1);
587 Standard_Real Step = .1, Delta = 0.;
c6541a0c 588 Standard_Real eps = M_PI, dmax = 0., d = 0.;
7fd59977 589 Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
590 nbp = Max(nbp, 3);
591 Step = (W2 - W1) / (nbp - 1);
592 myU1 = U1; myU2 = U1;
593 Standard_Real pmin = W1, pmax = W1, plim = W2+.1*Step;
594 for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
595 P = myCurve->Value(par);
596 ElSLib::Parameters( SP, P, U, V);
597 U += Delta;
598 d = U - U1;
599 if(d > eps) {
600 U -= Delta;
c6541a0c 601 Delta -= 2*M_PI;
7fd59977 602 U += Delta;
603 d = U - U1;
604 }
605 else if(d < -eps) {
606 U -= Delta;
c6541a0c 607 Delta += 2*M_PI;
7fd59977 608 U += Delta;
609 d = U - U1;
610 }
611 dmax = Max(dmax, Abs(d));
612 if(U < myU1) {myU1 = U; pmin = par;}
613 if(U > myU2) {myU2 = U; pmax = par;}
614 U1 = U;
615 }
616
617 if(!(Abs(pmin - W1) <= Precision::PConfusion() ||
618 Abs(pmin - W2) <= Precision::PConfusion()) ) myU1 -= dmax*.5;
619 if(!(Abs(pmax - W1) <= Precision::PConfusion() ||
620 Abs(pmax - W2) <= Precision::PConfusion()) ) myU2 += dmax*.5;
621
c6541a0c
D		 622 if((myU1 >=0. && myU1 <= 2*M_PI) &&
623 (myU2 >=0. && myU2 <= 2*M_PI) ) {
7fd59977 624 myU1 = 0.;
c6541a0c 625 myU2 = 2.*M_PI;
7fd59977 626 UCouture = Standard_False;
627 }
628 else {
629 U = ( myU1 + myU2 ) /2.;
c6541a0c
D		 630 myU1 = U - M_PI;
631 myU2 = U + M_PI;
7fd59977 632 UCouture = Standard_True;
633 }
634
635 VCouture = Standard_False;
636 }
637 }
638 break;
639 //
640 case GeomAbs_Torus:{
641 gp_Torus TR = mySurface->Torus();
642 Standard_Real U1, V1, U , V;
643 ElSLib::Parameters( TR, P1, U1, V1);
644 Standard_Real Step = .1, DeltaU = 0., DeltaV = 0.;
c6541a0c 645 Standard_Real eps = M_PI, dmaxU = 0., dU = 0., dmaxV = 0., dV = 0.;
7fd59977 646 Standard_Integer nbp = (Standard_Integer)((W2 - W1) / Step + 1);
647 nbp = Max(nbp, 3);
648 Step = (W2 - W1) / (nbp - 1);
649 myU1 = U1; myU2 = U1;
650 myV1 = V1; myV2 = V1;
651 Standard_Real pminU = W1, pmaxU = W1, pminV = W1, pmaxV = W1,
652 plim = W2+.1*Step;
653 for(Standard_Real par = W1 + Step; par <= plim; par += Step) {
654 P = myCurve->Value(par);
655 ElSLib::Parameters( TR, P, U, V);
656 U += DeltaU;
657 V += DeltaV;
658 dU = U - U1;
659 dV = V - V1;
660 if(dU > eps) {
661 U -= DeltaU;
c6541a0c 662 DeltaU -= 2*M_PI;
7fd59977 663 U += DeltaU;
664 dU = U - U1;
665 }
666 else if(dU < -eps) {
667 U -= DeltaU;
c6541a0c 668 DeltaU += 2*M_PI;
7fd59977 669 U += DeltaU;
670 dU = U - U1;
671 }
672 if(dV > eps) {
673 V -= DeltaV;
c6541a0c 674 DeltaV -= 2*M_PI;
7fd59977 675 V += DeltaV;
676 dV = V - V1;
677 }
678 else if(dV < -eps) {
679 V -= DeltaV;
c6541a0c 680 DeltaV += 2*M_PI;
7fd59977 681 V += DeltaV;
682 dV = V - V1;
683 }
684 dmaxU = Max(dmaxU, Abs(dU));
685 dmaxV = Max(dmaxV, Abs(dV));
686 if(U < myU1) {myU1 = U; pminU = par;}
687 if(U > myU2) {myU2 = U; pmaxU = par;}
688 if(V < myV1) {myV1 = V; pminV = par;}
689 if(V > myV2) {myV2 = V; pmaxV = par;}
690 U1 = U;
691 V1 = V;
692 }
693
694 if(!(Abs(pminU - W1) <= Precision::PConfusion() ||
695 Abs(pminU - W2) <= Precision::PConfusion()) ) myU1 -= dmaxU*.5;
696 if(!(Abs(pmaxU - W1) <= Precision::PConfusion() ||
697 Abs(pmaxU - W2) <= Precision::PConfusion()) ) myU2 += dmaxU*.5;
698 if(!(Abs(pminV - W1) <= Precision::PConfusion() ||
699 Abs(pminV - W2) <= Precision::PConfusion()) ) myV1 -= dmaxV*.5;
700 if(!(Abs(pmaxV - W1) <= Precision::PConfusion() ||
701 Abs(pmaxV - W2) <= Precision::PConfusion()) ) myV2 += dmaxV*.5;
702
c6541a0c
D		 703 if((myU1 >=0. && myU1 <= 2*M_PI) &&
704 (myU2 >=0. && myU2 <= 2*M_PI) ) {
7fd59977 705 myU1 = 0.;
c6541a0c 706 myU2 = 2.*M_PI;
7fd59977 707 UCouture = Standard_False;
708 }
709 else {
710 U = ( myU1 + myU2 ) /2.;
c6541a0c
D		 711 myU1 = U - M_PI;
712 myU2 = U + M_PI;
7fd59977 713 UCouture = Standard_True;
714 }
c6541a0c
D		 715 if((myV1 >=0. && myV1 <= 2*M_PI) &&
716 (myV2 >=0. && myV2 <= 2*M_PI) ) {
7fd59977 717 VCouture = Standard_False;
718 }
719 else {
720 V = ( myV1 + myV2 ) /2.;
c6541a0c
D		 721 myV1 = V - M_PI;
722 myV2 = V + M_PI;
7fd59977 723 VCouture = Standard_True;
724 }
725
726 }
727 break;
728
729 default:
730 {
731 UCouture = Standard_False;
732 VCouture = Standard_False;
733 }
734 break;
735 }
736}
737//
738//=======================================================================
739//function : Parameters
740//purpose :
741//=======================================================================
742void Parameters(const Handle(Adaptor3d_HCurve)& myCurve,
743 const Handle(Adaptor3d_HSurface)& mySurface,
744 const gp_Pnt& aP1,
745 const Standard_Integer iFirst,
746 const Standard_Real aTolU,
747 Standard_Real& aU,
748 Standard_Real& aV)
749{
750 Standard_Real aTwoPI, aU1, aV1, aU2, aV2, aRSp, aTol3D;
751 Standard_Real aTF, aTL, aT2, dT;
752 GeomAbs_SurfaceType aSType;
753 GeomAbs_CurveType aCType;
754 gp_Pnt aP2;
755 //
c6541a0c 756 aTwoPI=2.*M_PI;
7fd59977 757 //
758 aSType=mySurface->GetType();
759 aCType=myCurve->GetType();
760 //
761 if (aSType==GeomAbs_Sphere && aCType==GeomAbs_Circle) {
762 gp_Sphere aSp=mySurface->Sphere();
763 //
764 aRSp=aSp.Radius();
765 aTol3D=aRSp*aTolU;
766 //
767 aTF = myCurve->FirstParameter();
768 aTL = myCurve->LastParameter ();
769 dT=myCurve->Resolution(aTol3D);
770 //
771 ElSLib::Parameters(aSp, aP1, aU1, aV1);
772 if (fabs(aU)<aTolU || fabs(aU-aTwoPI)<aTolU){
773 aT2=aTF+dT;
774 if (!iFirst) {
775 aT2=aTL-dT;
776 }
777 //
778 aP2=myCurve->Value(aT2);
779 ElSLib::Parameters(aSp, aP2, aU2, aV2);
780 //
781 aU1=0.;
c6541a0c 782 if (aU2>M_PI) {
7fd59977 783 aU1=aTwoPI;
784 }
785 }
786 aU=aU1;
787 aV=aV1;
788 }
789}
790//
791//=======================================================================
792//classn : ProjLib_Function
793//purpose :
794//=======================================================================
795class ProjLib_Function : public AppCont_Function2d
796{
797 Handle(Adaptor3d_HCurve) myCurve;
798 Handle(Adaptor3d_HSurface) mySurface;
799
800 public :
801
802 Standard_Real myU1,myU2,myV1,myV2;
803 Standard_Boolean UCouture,VCouture;
804
805 ProjLib_Function(const Handle(Adaptor3d_HCurve)& C,
806 const Handle(Adaptor3d_HSurface)& S) :
807 myCurve(C), mySurface(S),
808 myU1(0.0),
809 myU2(0.0),
810 myV1(0.0),
811 myV2(0.0),
812 UCouture(Standard_False),
813 VCouture(Standard_False)
814 {Function_SetUVBounds(myU1,myU2,myV1,myV2,UCouture,VCouture,myCurve,mySurface);}
815
816 Standard_Real FirstParameter() const
817 {return (myCurve->FirstParameter() + 1.e-9);}
818
819 Standard_Real LastParameter() const
820 {return (myCurve->LastParameter() -1.e-9);}
821
822
823 gp_Pnt2d Value(const Standard_Real t) const
824 {return Function_Value(t,myCurve,mySurface,myU1,myU2,myV1,myV2,UCouture,VCouture);}
825
826 Standard_Boolean D1(const Standard_Real t, gp_Pnt2d& P, gp_Vec2d& V) const
827 {return Function_D1(t,P,V,myCurve,mySurface,myU1,myU2,myV1,myV2,UCouture,VCouture);}
828};
829
830//=======================================================================
831//function : ProjLib_ComputeApprox
832//purpose :
833//=======================================================================
834
835ProjLib_ComputeApprox::ProjLib_ComputeApprox
836 (const Handle(Adaptor3d_HCurve) & C,
837 const Handle(Adaptor3d_HSurface) & S,
838 const Standard_Real Tol )
839{
840 // if the surface is a plane and the curve a BSpline or a BezierCurve,
841 // don`t make an Approx but only the projection of the poles.
842
843 myTolerance = Max(Precision::PApproximation(),Tol);
844 Standard_Integer NbKnots, NbPoles ;
845 GeomAbs_CurveType CType = C->GetType();
846 GeomAbs_SurfaceType SType = S->GetType();
847
848 Standard_Boolean SurfIsAnal = (SType != GeomAbs_BSplineSurface) &&
849 (SType != GeomAbs_BezierSurface) &&
850 (SType != GeomAbs_OtherSurface) ;
851
852 Standard_Boolean CurvIsAnal = (CType != GeomAbs_BSplineCurve) &&
853 (CType != GeomAbs_BezierCurve) &&
854 (CType != GeomAbs_OtherCurve) ;
855
856 Standard_Boolean simplecase = SurfIsAnal && CurvIsAnal;
857
858 if (CType == GeomAbs_BSplineCurve &&
859 SType == GeomAbs_Plane ) {
860
861 // get the poles and eventually the weights
862 Handle(Geom_BSplineCurve) BS = C->BSpline();
863 NbPoles = BS->NbPoles();
864 TColgp_Array1OfPnt P3d( 1, NbPoles);
865 TColgp_Array1OfPnt2d Poles( 1, NbPoles);
866 TColStd_Array1OfReal Weights( 1, NbPoles);
867 if ( BS->IsRational()) BS->Weights(Weights);
868 BS->Poles( P3d);
869 gp_Pln Plane = S->Plane();
870 Standard_Real U,V;
871 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
872 ElSLib::Parameters( Plane, P3d(i), U, V);
873 Poles.SetValue(i,gp_Pnt2d(U,V));
874 }
875 NbKnots = BS->NbKnots();
876 TColStd_Array1OfReal Knots(1,NbKnots);
877 TColStd_Array1OfInteger Mults(1,NbKnots);
878 BS->Knots(Knots) ;
879 BS->Multiplicities(Mults) ;
880 // get the knots and mults if BSplineCurve
881 if ( BS->IsRational()) {
882 myBSpline = new Geom2d_BSplineCurve(Poles,
883 Weights,
884 Knots,
885 Mults,
886 BS->Degree(),
887 BS->IsPeriodic());
888 }
889 else {
890 myBSpline = new Geom2d_BSplineCurve(Poles,
891 Knots,
892 Mults,
893 BS->Degree(),
894 BS->IsPeriodic());
895 }
896 }
897 else if (CType == GeomAbs_BezierCurve &&
898 SType == GeomAbs_Plane ) {
899
900 // get the poles and eventually the weights
901 Handle(Geom_BezierCurve) BezierCurvePtr = C->Bezier() ;
902 NbPoles = BezierCurvePtr->NbPoles();
903 TColgp_Array1OfPnt P3d( 1, NbPoles);
904 TColgp_Array1OfPnt2d Poles( 1, NbPoles);
905 TColStd_Array1OfReal Weights( 1, NbPoles);
906 if ( BezierCurvePtr->IsRational()) {
907 BezierCurvePtr->Weights(Weights);
908 }
909 BezierCurvePtr->Poles( P3d);
910
911 // project the 3D-Poles on the plane
912
913 gp_Pln Plane = S->Plane();
914 Standard_Real U,V;
915 for ( Standard_Integer i = 1; i <= NbPoles; i++) {
916 ElSLib::Parameters( Plane, P3d(i), U, V);
917 Poles.SetValue(i,gp_Pnt2d(U,V));
918 }
919 if ( BezierCurvePtr->IsRational()) {
920 myBezier = new Geom2d_BezierCurve(Poles, Weights);
921 }
922 else {
923 myBezier = new Geom2d_BezierCurve(Poles);
924 }
925 }
926 else {
927 ProjLib_Function F( C, S);
928
929#ifdef DEB
930 if ( AffichValue) {
931 Standard_Integer Nb = 20;
932 Standard_Real U1, U2, dU, U;
933 U1 = F.FirstParameter();
934 U2 = F.LastParameter();
935 dU = ( U2 - U1) / Nb;
936 TColStd_Array1OfInteger Mults(1,Nb+1);
937 TColStd_Array1OfReal Knots(1,Nb+1);
938 TColgp_Array1OfPnt2d Poles(1,Nb+1);
939 for ( Standard_Integer i = 1; i <= Nb+1; i++) {
940 U = U1 + (i-1)*dU;
941 Poles(i) = F.Value(U);
942 Knots(i) = i;
943 Mults(i) = 1;
944 }
945 Mults(1) = 2;
946 Mults(Nb+1) = 2;
947#ifdef DRAW
948// POP pour NT
949 char* ResultName = "Result";
950 DrawTrSurf::Set(ResultName,new Geom2d_BSplineCurve(Poles,Knots,Mults,1));
951// DrawTrSurf::Set("Result",new Geom2d_BSplineCurve(Poles,Knots,Mults,1));
952#endif
953 }
954#endif
955
956//-----------
957 Standard_Integer Deg1, Deg2;
958 if(simplecase) {
959 Deg1 = 8;
960 Deg2 = 10;
961 }
962 else {
963 Deg1 = 8;
964 Deg2 = 12;
965 }
966//-------------
967 Approx_FitAndDivide2d Fit(F,Deg1,Deg2,myTolerance,myTolerance,
968 Standard_True);
969 if(Fit.IsAllApproximated()) {
970 Standard_Integer i;
971 Standard_Integer NbCurves = Fit.NbMultiCurves();
972
973 // on essaie de rendre la courbe au moins C1
974 Convert_CompBezierCurves2dToBSplineCurve2d Conv;
975
976 myTolerance = 0;
977 Standard_Real Tol3d,Tol2d;
978 for (i = 1; i <= NbCurves; i++) {
979 Fit.Error(i,Tol3d, Tol2d);
980 myTolerance = Max(myTolerance, Tol2d);
981 AppParCurves_MultiCurve MC = Fit.Value( i); //Charge la Ieme Curve
982 TColgp_Array1OfPnt2d Poles2d( 1, MC.Degree() + 1);//Recupere les poles
983 MC.Curve(1, Poles2d);
984
985 Conv.AddCurve(Poles2d);
986 }
987
988 //mise a jour des fields de ProjLib_Approx
989 Conv.Perform();
990
991 NbPoles = Conv.NbPoles();
992 NbKnots = Conv.NbKnots();
993
994 //7626
995 if(NbPoles <= 0 || NbPoles > 100000)
996 return;
997 if(NbKnots <= 0 || NbKnots > 100000)
998 return;
999
1000 TColgp_Array1OfPnt2d NewPoles(1,NbPoles);
1001 TColStd_Array1OfReal NewKnots(1,NbKnots);
1002 TColStd_Array1OfInteger NewMults(1,NbKnots);
1003
1004 Conv.KnotsAndMults(NewKnots,NewMults);
1005 Conv.Poles(NewPoles);
1006
1007 BSplCLib::Reparametrize(C->FirstParameter(),
1008 C->LastParameter(),
1009 NewKnots);
1010
1011 // il faut recadrer les poles de debut et de fin:
1012 // ( Car pour les problemes de couture, on a du ouvrir l`intervalle
1013 // de definition de la courbe.)
1014 // On choisit de calculer ces poles par prolongement de la courbe
1015 // approximee.
1016
1017 gp_Pnt2d P;
1018 Standard_Real U;
1019
1020 U = C->FirstParameter() - 1.e-9;
1021 BSplCLib::D0(U,
1022 0,
1023 Conv.Degree(),
1024 Standard_False,
1025 NewPoles,
1026 BSplCLib::NoWeights(),
1027 NewKnots,
1028 NewMults,
1029 P);
1030 NewPoles.SetValue(1,P);
1031 U = C->LastParameter() + 1.e-9;
1032 BSplCLib::D0(U,
1033 0,
1034 Conv.Degree(),
1035 Standard_False,
1036 NewPoles,
1037 BSplCLib::NoWeights(),
1038 NewKnots,
1039 NewMults,
1040 P);
1041 NewPoles.SetValue(NbPoles,P);
1042 myBSpline = new Geom2d_BSplineCurve (NewPoles,
1043 NewKnots,
1044 NewMults,
1045 Conv.Degree());
1046 }
1047 else {
1048 Standard_Integer NbCurves = Fit.NbMultiCurves();
1049 if(NbCurves != 0) {
1050 Standard_Real Tol3d,Tol2d;
1051 Fit.Error(NbCurves,Tol3d, Tol2d);
1052 myTolerance = Tol2d;
1053 }
1054 }
1055
1056 //Return curve home
1057 Standard_Real UFirst = F.FirstParameter();
1058 gp_Pnt P3d = C->Value( UFirst );
1059 Standard_Real u, v;
1060 switch (SType)
1061 {
1062 case GeomAbs_Plane:
1063 {
1064 gp_Pln Plane = S->Plane();
1065 ElSLib::Parameters( Plane, P3d, u, v );
1066 break;
1067 }
1068 case GeomAbs_Cylinder:
1069 {
1070 gp_Cylinder Cylinder = S->Cylinder();
1071 ElSLib::Parameters( Cylinder, P3d, u, v );
1072 break;
1073 }
1074 case GeomAbs_Cone:
1075 {
1076 gp_Cone Cone = S->Cone();
1077 ElSLib::Parameters( Cone, P3d, u, v );
1078 break;
1079 }
1080 case GeomAbs_Sphere:
1081 {
1082 gp_Sphere Sphere = S->Sphere();
1083 ElSLib::Parameters( Sphere, P3d, u, v );
1084 break;
1085 }
1086 case GeomAbs_Torus:
1087 {
1088 gp_Torus Torus = S->Torus();
1089 ElSLib::Parameters( Torus, P3d, u, v );
1090 break;
1091 }
1092 default:
1093 Standard_NoSuchObject::Raise("ProjLib_ComputeApprox::Value");
1094 }
1095 Standard_Boolean ToMirror = Standard_False;
1096 Standard_Real du = 0., dv = 0.;
1097 Standard_Integer number;
1098 if (F.VCouture)
1099 {
c6541a0c 1100 if (SType == GeomAbs_Sphere && Abs(u-F.myU1) > M_PI)
7fd59977 1101 {
1102 ToMirror = Standard_True;
c6541a0c
D		 1103 dv = -M_PI;
1104 v = M_PI - v;
7fd59977 1105 }
1106 Standard_Real newV = ElCLib::InPeriod( v, F.myV1, F.myV2 );
1107 number = (Standard_Integer) (Floor((newV-v)/(F.myV2-F.myV1)));
1108 dv -= number*(F.myV2-F.myV1);
1109 }
1110 if (F.UCouture || F.VCouture && SType == GeomAbs_Sphere)
1111 {
1112 gp_Pnt2d P2d = F.Value( UFirst );
c6541a0c
D		 1113 number = (Standard_Integer) (Floor((P2d.X()-u)/M_PI + Epsilon(M_PI)));
1114 du = -number*M_PI;
7fd59977 1115 }
1116
1117 if (!myBSpline.IsNull())
1118 {
1119 if (du != 0. || dv != 0.)
1120 myBSpline->Translate( gp_Vec2d(du,dv) );
1121 if (ToMirror)
1122 {
1123 gp_Ax2d Axe( gp_Pnt2d(0.,0.), gp_Dir2d(1.,0.) );
1124 myBSpline->Mirror( Axe );
1125 }
1126 }
1127 }
1128}
1129
1130//=======================================================================
1131//function : BSpline
1132//purpose :
1133//=======================================================================
1134
1135Handle(Geom2d_BSplineCurve) ProjLib_ComputeApprox::BSpline() const
1136
1137{
1138 return myBSpline ;
1139}
1140
1141//=======================================================================
1142//function : Bezier
1143//purpose :
1144//=======================================================================
1145
1146Handle(Geom2d_BezierCurve) ProjLib_ComputeApprox::Bezier() const
1147
1148{
1149 return myBezier ;
1150}
1151
1152
1153//=======================================================================
1154//function : Tolerance
1155//purpose :
1156//=======================================================================
1157
1158Standard_Real ProjLib_ComputeApprox::Tolerance() const
1159{
1160 return myTolerance;
1161}
1162
1163
Open CASCADE Technology repository