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