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