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