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