0023464: Projection algorithm produces wrong results.
[occt.git] / src / ProjLib / ProjLib_ProjectedCurve.cxx
CommitLineData
b311480e 1// Created on: 1993-08-25
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
21
7fd59977 22
23// Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272
24
25#include <GeomAbs_SurfaceType.hxx>
26#include <Standard_NoSuchObject.hxx>
27#include <Standard_NotImplemented.hxx>
28#include <ProjLib_ProjectedCurve.hxx>
29#include <ProjLib_CompProjectedCurve.hxx>
30#include <ProjLib_HCompProjectedCurve.hxx>
31#include <ProjLib_ComputeApproxOnPolarSurface.hxx>
32#include <ProjLib_ComputeApprox.hxx>
33#include <ProjLib_Projector.hxx>
34#include <Handle_Adaptor3d_HCurve.hxx>
35#include <Handle_Adaptor3d_HSurface.hxx>
36#include <Adaptor3d_HCurve.hxx>
37#include <Adaptor3d_HSurface.hxx>
38#include <Approx_CurveOnSurface.hxx>
39#include <ProjLib_Plane.hxx>
40#include <ProjLib_Cylinder.hxx>
41#include <ProjLib_Cone.hxx>
42#include <ProjLib_Sphere.hxx>
43#include <ProjLib_Torus.hxx>
44#include <Precision.hxx>
45#include <Handle_Geom_BSplineCurve.hxx>
46#include <Geom2d_BSplineCurve.hxx>
47#include <Handle_Geom2d_BSplineCurve.hxx>
48#include <Geom2d_BezierCurve.hxx>
49#include <Handle_Geom2d_BezierCurve.hxx>
50#include <Handle_Adaptor2d_HCurve2d.hxx>
51#include <gp_Vec2d.hxx>
52#include <StdFail_NotDone.hxx>
53#include <gp_XY.hxx>
54#include <TColgp_HArray1OfPnt2d.hxx>
55#include <TColStd_HArray1OfReal.hxx>
56#include <Geom2dConvert_CompCurveToBSplineCurve.hxx>
57#include <TColStd_Array1OfReal.hxx>
58#include <TColStd_Array1OfInteger.hxx>
59#include <TColgp_Array1OfPnt2d.hxx>
60#include <TColgp_HArray1OfVec2d.hxx>
61#include <TColStd_HArray1OfBoolean.hxx>
62#include <BSplCLib.hxx>
63#include <GeomAbs_IsoType.hxx>
bd82d4b2 64#include <Geom2d_Line.hxx>
65#include <Geom2d_TrimmedCurve.hxx>
66#include <ElCLib.hxx>
67#include <GeomLib.hxx>
7fd59977 68
69//=======================================================================
70//function : IsoIsDeg
71//purpose :
72//=======================================================================
73
74static Standard_Boolean IsoIsDeg (const Adaptor3d_Surface& S,
75 const Standard_Real Param,
76 const GeomAbs_IsoType IT,
77 const Standard_Real TolMin,
78 const Standard_Real TolMax)
79{
80 Standard_Real U1=0.,U2=0.,V1=0.,V2=0.,T;
81 Standard_Boolean Along = Standard_True;
82 U1 = S.FirstUParameter();
83 U2 = S.LastUParameter();
84 V1 = S.FirstVParameter();
85 V2 = S.LastVParameter();
86 gp_Vec D1U,D1V;
87 gp_Pnt P;
88 Standard_Real Step,D1NormMax;
89 if (IT == GeomAbs_IsoV)
90 {
91 Step = (U2 - U1)/10;
92 D1NormMax=0.;
93 for (T=U1;T<=U2;T=T+Step)
94 {
95 S.D1(T,Param,P,D1U,D1V);
96 D1NormMax=Max(D1NormMax,D1U.Magnitude());
97 }
98
99 if (D1NormMax >TolMax || D1NormMax < TolMin )
100 Along = Standard_False;
101 }
102 else
103 {
104 Step = (V2 - V1)/10;
105 D1NormMax=0.;
106 for (T=V1;T<=V2;T=T+Step)
107 {
108 S.D1(Param,T,P,D1U,D1V);
109 D1NormMax=Max(D1NormMax,D1V.Magnitude());
110 }
111
112 if (D1NormMax >TolMax || D1NormMax < TolMin )
113 Along = Standard_False;
114
115
116 }
117 return Along;
118}
119
120//=======================================================================
121//function : Interpolate
122//purpose :
123//=======================================================================
124
125static Handle(Geom2d_BSplineCurve) Interpolate(const Handle(TColgp_HArray1OfPnt2d)& myPoints,
126 const Handle(TColStd_HArray1OfReal)& myParameters,
127 const gp_Vec2d& InitialTangent,
128 const gp_Vec2d& FinalTangent)
129{
130 Handle(Geom2d_BSplineCurve) myCurve = NULL;
131
132// This code is extraction from Geom2dAPI_Interpolate with small correction
133// This is done to avoid of cyclic dependency if Geom2dAPI is used in ProjLib
134
135 Standard_Integer degree,
136 ii,
137 jj,
138 index,
139 index1,
140 index2,
141 index3,
142 mult_index,
143 inversion_problem,
144 num_points,
145 num_distinct_knots,
146 num_poles ;
147
148 gp_Pnt2d a_point ;
149
150 Standard_Boolean myTangentRequest = Standard_True;
151
152 Handle(TColgp_HArray1OfVec2d) myTangents =
153 new TColgp_HArray1OfVec2d(myPoints->Lower(),
154 myPoints->Upper()) ;
155 Handle(TColStd_HArray1OfBoolean) myTangentFlags =
156 new TColStd_HArray1OfBoolean(myPoints->Lower(),
157 myPoints->Upper()) ;
158 myTangentFlags->Init(Standard_False);
159
160 myTangentFlags->SetValue(1,Standard_True) ;
161 myTangentFlags->SetValue(myPoints->Length(),Standard_True) ;
162 myTangents->SetValue(1,InitialTangent) ;
163 myTangents->SetValue(myPoints->Length(),FinalTangent);
164
165 num_points =
166 num_distinct_knots =
167 num_poles = myPoints->Length() ;
168 if (num_poles == 2 && !myTangentRequest) {
169 degree = 1 ;
170 }
171 else if (num_poles == 3 && !myTangentRequest) {
172 degree = 2 ;
173 num_distinct_knots = 2 ;
174 }
175 else {
176 degree = 3 ;
177 num_poles += 2 ;
178 //if (myTangentRequest)
179 //for (ii = myTangentFlags->Lower() + 1 ;
180 // ii < myTangentFlags->Upper() ; ii++) {
181 //if (myTangentFlags->Value(ii)) {
182 //num_poles += 1 ;
183 //}
184 //}
185 }
186
187
188 TColStd_Array1OfReal parameters(1,num_poles) ;
189 TColStd_Array1OfReal flatknots(1,num_poles + degree + 1) ;
190 TColStd_Array1OfInteger mults(1,num_distinct_knots) ;
191 TColStd_Array1OfReal knots(1,num_distinct_knots) ;
192 TColStd_Array1OfInteger contact_order_array(1, num_poles) ;
193 TColgp_Array1OfPnt2d poles(1,num_poles) ;
194
195 for (ii = 1 ; ii <= degree + 1 ; ii++) {
196 flatknots.SetValue(ii,myParameters->Value(1)) ;
197 flatknots.SetValue(ii + num_poles,
198 myParameters->Value(num_points)) ;
199 }
200 for (ii = 1 ; ii <= num_poles ; ii++) {
201 contact_order_array.SetValue(ii,0) ;
202 }
203 for (ii = 2 ; ii < num_distinct_knots ; ii++) {
204 mults.SetValue(ii,1) ;
205 }
206 mults.SetValue(1,degree + 1) ;
207 mults.SetValue(num_distinct_knots ,degree + 1) ;
208
209 switch (degree) {
210 case 1:
211 for (ii = 1 ; ii <= num_poles ; ii++) {
212 poles.SetValue(ii ,myPoints->Value(ii)) ;
213 }
214 myCurve =
215 new Geom2d_BSplineCurve(poles,
216 myParameters->Array1(),
217 mults,
218 degree) ;
219 //myIsDone = Standard_True ;
220 break ;
221 case 2:
222 knots.SetValue(1,myParameters->Value(1)) ;
223 knots.SetValue(2,myParameters->Value(3)) ;
224 for (ii = 1 ; ii <= num_poles ; ii++) {
225 poles.SetValue(ii,myPoints->Value(ii)) ;
226
227 }
228 BSplCLib::Interpolate(degree,
229 flatknots,
230 myParameters->Array1(),
231 contact_order_array,
232 poles,
233 inversion_problem) ;
234 if (!inversion_problem) {
235 myCurve =
236 new Geom2d_BSplineCurve(poles,
237 knots,
238 mults,
239 degree) ;
240 //myIsDone = Standard_True ;
241 }
242 break ;
243 case 3:
244//
245// check if the boundary conditions are set
246//
247 //if (num_points >= 3) {
248//
249// cannot build the tangents with degree 3 with only 2 points
250// if those where not given in advance
251//
252 //BuildTangents(myPoints->Array1(),
253 //myTangents->ChangeArray1(),
254 //myTangentFlags->ChangeArray1(),
255 //myParameters->Array1()) ;
256 //}
257 contact_order_array.SetValue(2,1) ;
258 parameters.SetValue(1,myParameters->Value(1)) ;
259 parameters.SetValue(2,myParameters->Value(1)) ;
260 poles.SetValue(1,myPoints->Value(1)) ;
261 for (jj = 1 ; jj <= 2 ; jj++) {
262 a_point.SetCoord(jj,myTangents->Value(1).Coord(jj)) ;
263
264 }
265 poles.SetValue(2,a_point) ;
266 mult_index = 2 ;
267 index = 3 ;
268 index1 = 2 ;
269 index2 = myPoints->Lower() + 1 ;
270 index3 = degree + 2 ;
271 if (myTangentRequest) {
272 for (ii = myParameters->Lower() + 1 ;
273 ii < myParameters->Upper() ; ii++) {
274 parameters.SetValue(index,myParameters->Value(ii)) ;
275 poles.SetValue(index,myPoints->Value(index2)) ;
276 flatknots.SetValue(index3,myParameters->Value(ii)) ;
277 index += 1 ;
278 index3 += 1 ;
279 if (myTangentFlags->Value(index1)) {
280//
281// set the multiplicities, the order of the contact, the
282// the flatknots,
283//
284 mults.SetValue(mult_index,mults.Value(mult_index) + 1) ;
285 contact_order_array(index) = 1 ;
286 flatknots.SetValue(index3, myParameters->Value(ii)) ;
287 parameters.SetValue(index,
288 myParameters->Value(ii)) ;
289 for (jj = 1 ; jj <= 2 ; jj++) {
290 a_point.SetCoord(jj,myTangents->Value(ii).Coord(jj)) ;
291 }
292 poles.SetValue(index,a_point) ;
293 index += 1 ;
294 index3 += 1 ;
295 }
296 mult_index += 1 ;
297 index1 += 1 ;
298 index2 += 1 ;
299
300 }
301 }
302 else {
303 index1 = 2 ;
304 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
305 parameters.SetValue(index1,
306 myParameters->Value(ii)) ;
307 index1 += 1 ;
308 }
309 index = 3 ;
310 for (ii = myPoints->Lower() + 1 ; ii <= myPoints->Upper() - 1 ; ii++) {
311 poles.SetValue(index,
312 myPoints->Value(ii)) ;
313 index += 1 ;
314 }
315
316
317 index = degree + 1 ;
318 for(ii = myParameters->Lower() ; ii <= myParameters->Upper() ; ii++) {
319 flatknots.SetValue(index,
320 myParameters->Value(ii)) ;
321 index += 1 ;
322 }
323 }
324 for (jj = 1 ; jj <= 2 ; jj++) {
325 a_point.SetCoord(jj,
326 myTangents->Value(num_points).Coord(jj)) ;
327 }
328 poles.SetValue(num_poles-1 ,a_point) ;
329
330 contact_order_array.SetValue(num_poles - 1,1) ;
331 parameters.SetValue(num_poles,
332 myParameters->Value(myParameters->Upper())) ;
333 parameters.SetValue(num_poles -1,
334 myParameters->Value(myParameters->Upper())) ;
335
336 poles.SetValue(num_poles,
337 myPoints->Value(num_points)) ;
338
339 BSplCLib::Interpolate(degree,
340 flatknots,
341 parameters,
342 contact_order_array,
343 poles,
344 inversion_problem) ;
345 if (!inversion_problem) {
346 myCurve =
347 new Geom2d_BSplineCurve(poles,
348 myParameters->Array1(),
349 mults,
350 degree) ;
351 //myIsDone = Standard_True ;
352 }
353 break ;
354
355 }
356
357
358 return myCurve;
359}
360
361//=======================================================================
362//function : TrimC3d
363//purpose :
364//=======================================================================
365
366static void TrimC3d(Handle(Adaptor3d_HCurve)& myCurve,
367 Standard_Boolean* IsTrimmed,
368 const Standard_Real dt,
bd82d4b2 369 const gp_Pnt& Pole,
370 Standard_Integer* SingularCase,
371 const Standard_Integer NumberOfSingularCase)
7fd59977 372{
373 Standard_Real f = myCurve->FirstParameter();
374 Standard_Real l = myCurve->LastParameter();
375
376 gp_Pnt P = myCurve->Value(f);
377
378 if(P.Distance(Pole) < Precision::Confusion()) {
379 IsTrimmed[0] = Standard_True;
380 f = f+dt;
381 myCurve = myCurve->Trim(f, l, Precision::Confusion());
bd82d4b2 382 SingularCase[0] = NumberOfSingularCase;
7fd59977 383 }
384
385 P = myCurve->Value(l);
386 if(P.Distance(Pole) < Precision::Confusion()) {
387 IsTrimmed[1] = Standard_True;
388 l = l-dt;
389 myCurve = myCurve->Trim(f, l, Precision::Confusion());
bd82d4b2 390 SingularCase[1] = NumberOfSingularCase;
7fd59977 391 }
392}
393
394//=======================================================================
395//function : ExtendC2d
396//purpose :
397//=======================================================================
398
399static void ExtendC2d(Handle(Geom2d_BSplineCurve)& aRes,
400 const Standard_Real t,
401 const Standard_Real dt,
402 const Standard_Real u1,
403 const Standard_Real u2,
404 const Standard_Real v1,
bd82d4b2 405 const Standard_Real v2,
406 const Standard_Integer FirstOrLast,
407 const Standard_Integer NumberOfSingularCase)
7fd59977 408{
bd82d4b2 409 Standard_Real theParam = (FirstOrLast == 0)? aRes->FirstParameter()
410 : aRes->LastParameter();
411
412 gp_Pnt2d aPBnd;
413 gp_Vec2d aVBnd;
414 gp_Dir2d aDBnd;
415 Handle(Geom2d_TrimmedCurve) aSegment;
416 Geom2dConvert_CompCurveToBSplineCurve aCompCurve(aRes, Convert_RationalC1);
417 Standard_Real aTol = Precision::Confusion();
418
419 aRes->D1(theParam, aPBnd, aVBnd);
420 aDBnd.SetXY(aVBnd.XY());
421 gp_Lin2d aLin(aPBnd, aDBnd); //line in direction of derivative
422
423 gp_Pnt2d thePole;
424 gp_Dir2d theBoundDir;
425 switch (NumberOfSingularCase)
426 {
427 case 1:
428 {
429 thePole.SetCoord(u1, v1);
430 theBoundDir.SetCoord(0., 1.);
431 break;
432 }
433 case 2:
434 {
435 thePole.SetCoord(u2, v1);
436 theBoundDir.SetCoord(0., 1.);
437 break;
438 }
439 case 3:
440 {
441 thePole.SetCoord(u1, v1);
442 theBoundDir.SetCoord(1., 0.);
443 break;
444 }
445 case 4:
446 {
447 thePole.SetCoord(u1, v2);
448 theBoundDir.SetCoord(1., 0.);
449 break;
450 }
451 }
452 gp_Lin2d BoundLin(thePole, theBoundDir); //one of the bounds of rectangle
453
454 Standard_Real U1x = BoundLin.Direction().X();
455 Standard_Real U1y = BoundLin.Direction().Y();
456 Standard_Real U2x = aLin.Direction().X();
457 Standard_Real U2y = aLin.Direction().Y();
458 Standard_Real Uo21x = aLin.Location().X() - BoundLin.Location().X();
459 Standard_Real Uo21y = aLin.Location().Y() - BoundLin.Location().Y();
460
461 Standard_Real D = U1y*U2x-U1x*U2y;
462
463 Standard_Real ParOnLin = (Uo21y * U1x - Uo21x * U1y)/D; //parameter of intersection point
464
465 Handle(Geom2d_Line) aSegLine = new Geom2d_Line(aLin);
466 aSegment = (FirstOrLast == 0)?
467 new Geom2d_TrimmedCurve(aSegLine, ParOnLin, 0.) :
468 new Geom2d_TrimmedCurve(aSegLine, 0., ParOnLin);
469
470 aCompCurve.Add(aSegment, aTol);
471 aRes = aCompCurve.BSplineCurve();
472
473 /*
7fd59977 474 gp_Pnt2d P0;
475 gp_Vec2d V01, V02;
476 aRes->D2(t, P0, V01, V02);
477
478 gp_XY XYP1 = P0.XY() + V01.XY()*dt + .5*V02.XY()*dt*dt;
479
480 gp_Vec2d V11 = V01 + V02*dt;
481
482 if(XYP1.X() < u1) XYP1.SetX(u1);
483 if(XYP1.X() > u2) XYP1.SetX(u2);
484 if(XYP1.Y() < v1) XYP1.SetY(v1);
485 if(XYP1.Y() > v2) XYP1.SetY(v2);
486
487 Handle(TColgp_HArray1OfPnt2d) aPnts = new TColgp_HArray1OfPnt2d(1, 2);
488 Handle(TColStd_HArray1OfReal) aPars = new TColStd_HArray1OfReal(1, 2);
489
490 if(dt < 0.) {
491 aPnts->SetValue(1, gp_Pnt2d(XYP1));
492 aPnts->SetValue(2, P0);
493 aPars->SetValue(1, t + dt);
494 aPars->SetValue(2, t);
495 }
496 else {
497 aPnts->SetValue(2, gp_Pnt2d(XYP1));
498 aPnts->SetValue(1, P0);
499 aPars->SetValue(2, t + dt);
500 aPars->SetValue(1, t);
501 }
502
503 Handle(Geom2d_BSplineCurve) aC;
504
505 if(dt < 0.) {
506 aC = Interpolate(aPnts, aPars, V11, V01);
507 }
508 else {
509 aC = Interpolate(aPnts, aPars, V01, V11);
510 }
511
512
513 Geom2dConvert_CompCurveToBSplineCurve aConcat(aRes);
514 aConcat.Add(aC, Precision::PConfusion());
515
516 aRes = aConcat.BSplineCurve();
bd82d4b2 517 */
7fd59977 518}
519
520//=======================================================================
521//function : Project
522//purpose :
523//=======================================================================
524
525static void Project(ProjLib_Projector& P, Handle(Adaptor3d_HCurve)& C)
526{
527 GeomAbs_CurveType CType = C->GetType();
528 switch (CType) {
529 case GeomAbs_Line:
530 P.Project(C->Line());
531 break;
532 case GeomAbs_Circle:
533 P.Project(C->Circle());
534 break;
535 case GeomAbs_Ellipse:
536 P.Project(C->Ellipse());
537 break;
538 case GeomAbs_Hyperbola:
539 P.Project(C->Hyperbola());
540 break;
541 case GeomAbs_Parabola:
542 P.Project(C->Parabola());
543 break;
544 case GeomAbs_BSplineCurve:
545 case GeomAbs_BezierCurve:
546 case GeomAbs_OtherCurve: // try the approximation
547 break;
548 default:
549 Standard_NoSuchObject::Raise(" ");
550 }
551}
552
553//=======================================================================
554//function : ProjLib_ProjectedCurve
555//purpose :
556//=======================================================================
557
558ProjLib_ProjectedCurve::ProjLib_ProjectedCurve()
559
560{
561 myTolerance = Precision::PApproximation();
562}
563
564
565//=======================================================================
566//function : ProjLib_ProjectedCurve
567//purpose :
568//=======================================================================
569
570ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
571(const Handle(Adaptor3d_HSurface)& S)
572{
573 myTolerance = Precision::PApproximation();
574 Load(S);
575}
576
577
578//=======================================================================
579//function : ProjLib_ProjectedCurve
580//purpose :
581//=======================================================================
582
583ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
584(const Handle(Adaptor3d_HSurface)& S,
585 const Handle(Adaptor3d_HCurve)& C)
586{
587 myTolerance = Precision::PApproximation();
588 Load(S);
589 Load(C);
590}
591
592
593//=======================================================================
594//function : ProjLib_ProjectedCurve
595//purpose :
596//=======================================================================
597
598ProjLib_ProjectedCurve::ProjLib_ProjectedCurve
599(const Handle(Adaptor3d_HSurface)& S,
600 const Handle(Adaptor3d_HCurve)& C,
601 const Standard_Real Tol)
602{
603 myTolerance = Max(Tol, Precision::PApproximation());
604 Load(S);
605 Load(C);
606}
607
608
609//=======================================================================
610//function : Load
611//purpose :
612//=======================================================================
613
614void ProjLib_ProjectedCurve::Load(const Handle(Adaptor3d_HSurface)& S)
615{
616 mySurface = S ;
617}
618
619
620//=======================================================================
621//function : Load
622//purpose :
623//=======================================================================
624
625void ProjLib_ProjectedCurve::Load(const Handle(Adaptor3d_HCurve)& C)
626{
627 myTolerance = Max(myTolerance,Precision::PApproximation());
628 myCurve = C;
bd82d4b2 629 Standard_Real FirstPar = C->FirstParameter();
630 Standard_Real LastPar = C->LastParameter();
7fd59977 631 GeomAbs_SurfaceType SType = mySurface->GetType();
632 GeomAbs_CurveType CType = myCurve->GetType();
633
634 switch (SType) {
635
636 case GeomAbs_Plane:
637 {
638 ProjLib_Plane P(mySurface->Plane());
639 Project(P,myCurve);
640 myResult = P;
641 }
642 break;
643
644 case GeomAbs_Cylinder:
645 {
646 ProjLib_Cylinder P(mySurface->Cylinder());
647 Project(P,myCurve);
648 myResult = P;
649 }
650 break;
651
652 case GeomAbs_Cone:
653 {
654 ProjLib_Cone P(mySurface->Cone());
655 Project(P,myCurve);
656 myResult = P;
657 }
658 break;
659
660 case GeomAbs_Sphere:
661 {
662 ProjLib_Sphere P(mySurface->Sphere());
663 Project(P,myCurve);
664 if ( P.IsDone()) {
665 // on met dans la pseudo-periode ( car Sphere n'est pas
666 // periodique en V !)
667 P.SetInBounds(myCurve->FirstParameter());
668 }
669 myResult = P;
670 }
671 break;
672
673 case GeomAbs_Torus:
674 {
675 ProjLib_Torus P(mySurface->Torus());
676 Project(P,myCurve);
677 myResult = P;
678 }
679 break;
680
681 case GeomAbs_BezierSurface:
682 case GeomAbs_BSplineSurface:
683 {
684
685 Standard_Boolean IsTrimmed[2] = {Standard_False, Standard_False};
bd82d4b2 686 Standard_Integer SingularCase[2];
7fd59977 687 Standard_Real f, l, dt;
688 const Standard_Real eps = 0.01;
689 f = myCurve->FirstParameter();
690 l = myCurve->LastParameter();
691 dt = (l-f)*eps;
692
693 Standard_Real U1=0.,U2=0.,V1=0.,V2=0;
694 const Adaptor3d_Surface& S = mySurface->Surface();
695 U1 = S.FirstUParameter();
696 U2 = S.LastUParameter();
697 V1 = S.FirstVParameter();
698 V2 = S.LastVParameter();
699
074028c6 700 if(IsoIsDeg(S, U1, GeomAbs_IsoU, 0., myTolerance) ) {
7fd59977 701 //Surface has pole at U = Umin
702 gp_Pnt Pole = mySurface->Value(U1, V1);
bd82d4b2 703 TrimC3d(myCurve, IsTrimmed, dt, Pole, SingularCase, 1);
7fd59977 704 }
705
074028c6 706 if(IsoIsDeg(S, U2, GeomAbs_IsoU, 0., myTolerance) ) {
7fd59977 707 //Surface has pole at U = Umax
708 gp_Pnt Pole = mySurface->Value(U2, V1);
bd82d4b2 709 TrimC3d(myCurve, IsTrimmed, dt, Pole, SingularCase, 2);
7fd59977 710 }
711
074028c6 712 if(IsoIsDeg(S, V1, GeomAbs_IsoV, 0., myTolerance) ) {
7fd59977 713 //Surface has pole at V = Vmin
714 gp_Pnt Pole = mySurface->Value(U1, V1);
bd82d4b2 715 TrimC3d(myCurve, IsTrimmed, dt, Pole, SingularCase, 3);
7fd59977 716 }
717
074028c6 718 if(IsoIsDeg(S, V2, GeomAbs_IsoV, 0., myTolerance) ) {
7fd59977 719 //Surface has pole at V = Vmax
720 gp_Pnt Pole = mySurface->Value(U1, V2);
bd82d4b2 721 TrimC3d(myCurve, IsTrimmed, dt, Pole, SingularCase, 4);
7fd59977 722 }
723
724 ProjLib_ComputeApproxOnPolarSurface polar(myCurve,
725 mySurface,
726 myTolerance);
727
728 Handle(Geom2d_BSplineCurve) aRes = polar.BSpline();
729
730 if(IsTrimmed[0] || IsTrimmed[1]) {
731 if(IsTrimmed[0]) {
732 //Add segment before start of curve
733 f = myCurve->FirstParameter();
bd82d4b2 734 ExtendC2d(aRes, f, -dt, U1, U2, V1, V2, 0, SingularCase[0]);
7fd59977 735 }
736 if(IsTrimmed[1]) {
737 //Add segment after end of curve
738 l = myCurve->LastParameter();
bd82d4b2 739 ExtendC2d(aRes, l, dt, U1, U2, V1, V2, 1, SingularCase[1]);
7fd59977 740 }
bd82d4b2 741 Handle(Geom2d_Curve) NewCurve2d;
742 GeomLib::SameRange(Precision::PConfusion(), aRes,
743 aRes->FirstParameter(), aRes->LastParameter(),
744 FirstPar, LastPar,
745 NewCurve2d);
746 aRes = Handle(Geom2d_BSplineCurve)::DownCast(NewCurve2d);
747 }
7fd59977 748 myResult.SetBSpline(aRes);
749 myResult.Done();
750 myResult.SetType(GeomAbs_BSplineCurve);
751 }
752 break;
753
754 default:
755 {
756 Standard_Boolean IsTrimmed[2] = {Standard_False, Standard_False};
bd82d4b2 757 Standard_Real Vsingular[2]; //for surfaces of revolution
7fd59977 758 Standard_Real f, l, dt;
759 const Standard_Real eps = 0.01;
760
761 if(mySurface->GetType() == GeomAbs_SurfaceOfRevolution) {
762 //Check possible singularity
763
764 gp_Pnt P = mySurface->AxeOfRevolution().Location();
765 gp_Dir N = mySurface->AxeOfRevolution().Direction();
766
767 gp_Lin L(P, N);
768
769 f = myCurve->FirstParameter();
770 l = myCurve->LastParameter();
771 dt = (l-f)*eps;
772
773 P = myCurve->Value(f);
774 if(L.Distance(P) < Precision::Confusion()) {
775 IsTrimmed[0] = Standard_True;
776 f = f+dt;
777 myCurve = myCurve->Trim(f, l, Precision::Confusion());
bd82d4b2 778 Vsingular[0] = ElCLib::Parameter(L, P);
779 //SingularCase[0] = 3;
7fd59977 780 }
781
782 P = myCurve->Value(l);
783 if(L.Distance(P) < Precision::Confusion()) {
784 IsTrimmed[1] = Standard_True;
785 l = l-dt;
786 myCurve = myCurve->Trim(f, l, Precision::Confusion());
bd82d4b2 787 Vsingular[1] = ElCLib::Parameter(L, P);
788 //SingularCase[1] = 3;
7fd59977 789 }
790 }
791
792 ProjLib_CompProjectedCurve Projector(mySurface,myCurve,
793 myTolerance,myTolerance);
794 Handle(ProjLib_HCompProjectedCurve) HProjector =
795 new ProjLib_HCompProjectedCurve();
796 HProjector->Set(Projector);
797
798 // Normalement, dans le cadre de ProjLib, le resultat
799 // doit etre une et une seule courbe !!!
800 // De plus, cette courbe ne doit pas etre Single point
801 Standard_Integer NbCurves = Projector.NbCurves();
802 Standard_Real Udeb,Ufin;
803 if (NbCurves > 0) {
804 Projector.Bounds(1,Udeb,Ufin);
805 }
806 else {
807 StdFail_NotDone::Raise("ProjLib CompProjectedCurve Not Done");
808 }
809 // Approximons cette courbe algorithmique.
810 Standard_Boolean Only3d = Standard_False;
811 Standard_Boolean Only2d = Standard_True;
812 GeomAbs_Shape Continuity = GeomAbs_C1;
813 Standard_Integer MaxDegree = 14;
814 Standard_Integer MaxSeg = 16;
815
816 Approx_CurveOnSurface appr(HProjector, mySurface, Udeb, Ufin,
817 myTolerance,
818 Continuity, MaxDegree, MaxSeg,
819 Only3d, Only2d);
820
821 Handle(Geom2d_BSplineCurve) aRes = appr.Curve2d();
822
823 if(IsTrimmed[0] || IsTrimmed[1]) {
824 // Treatment only for surface of revolution
825 Standard_Real u1, u2, v1, v2;
826 u1 = mySurface->FirstUParameter();
827 u2 = mySurface->LastUParameter();
828 v1 = mySurface->FirstVParameter();
829 v2 = mySurface->LastVParameter();
830
831 if(IsTrimmed[0]) {
832 //Add segment before start of curve
bd82d4b2 833 ExtendC2d(aRes, f, -dt, u1, u2, Vsingular[0], v2, 0, 3);
7fd59977 834 }
835 if(IsTrimmed[1]) {
836 //Add segment after end of curve
bd82d4b2 837 ExtendC2d(aRes, l, dt, u1, u2, Vsingular[1], v2, 1, 3);
7fd59977 838 }
bd82d4b2 839 Handle(Geom2d_Curve) NewCurve2d;
840 GeomLib::SameRange(Precision::PConfusion(), aRes,
841 aRes->FirstParameter(), aRes->LastParameter(),
842 FirstPar, LastPar,
843 NewCurve2d);
844 aRes = Handle(Geom2d_BSplineCurve)::DownCast(NewCurve2d);
7fd59977 845 }
846
847 myResult.SetBSpline(aRes);
848 myResult.Done();
849 myResult.SetType(GeomAbs_BSplineCurve);
850 }
851 }
852 if ( !myResult.IsDone()) {
853 ProjLib_ComputeApprox Comp( myCurve, mySurface, myTolerance);
854 myResult.Done();
855
856 // set the type
857 if ( SType == GeomAbs_Plane && CType == GeomAbs_BezierCurve) {
858 myResult.SetType(GeomAbs_BezierCurve);
859 myResult.SetBezier(Comp.Bezier()) ;
860 }
861 else {
862 myResult.SetType(GeomAbs_BSplineCurve);
863 myResult.SetBSpline(Comp.BSpline()) ;
864 }
865 // set the periodicity flag
866 if ( SType == GeomAbs_Plane &&
867 CType == GeomAbs_BSplineCurve &&
868 myCurve->IsPeriodic() ) {
869 myResult.SetPeriodic();
870 }
871 myTolerance = Comp.Tolerance();
872 }
873
874 else {
875 // On remet arbitrairement la tol atteinte a une valeur
876 // petite en attendant mieux. dub lbo 11/03/97
877 myTolerance = Min(myTolerance,Precision::Confusion());
878
879 // Translate the projected curve to keep the first point
880 // In the canonical boundaries of periodic surfaces.
881 if (mySurface->IsUPeriodic()) {
882 // xf
883 Standard_Real aT1, aT2, aU1, aU2, aUPeriod, aUr, aUm, aUmid, dUm, dUr;
884 GeomAbs_CurveType aTypeR;
885 ProjLib_Projector aResult;
886 //
887 aT1=myCurve->FirstParameter();
888 aT2=myCurve->LastParameter();
889 aU1=mySurface->FirstUParameter();
890 aU2=mySurface->LastUParameter();
891 aUPeriod=mySurface->UPeriod();
892 //
893 aTypeR=myResult.GetType();
894 if ((aU2-aU1)<(aUPeriod-myTolerance) && aTypeR == GeomAbs_Line) {
895 aResult=myResult;
896 aResult.UFrame(aT1, aT2, aU1, aUPeriod);
897 //
898 gp_Lin2d &aLr = (gp_Lin2d &) aResult.Line();
899 aUr=aLr.Location().X();
900 gp_Lin2d &aLm = (gp_Lin2d &) myResult.Line();
901 aUm=aLm.Location().X();
902 //
903 aUmid=0.5*(aU2+aU1);
904 dUm=fabs(aUm-aUmid);
905 dUr=fabs(aUr-aUmid);
906 if (dUr<dUm) {
907 myResult=aResult;
908 }
909 }
910 else {
911 myResult.UFrame(aT1, aT2, aU1, aUPeriod);
912 }
913 //
914 /*
915 myResult.UFrame(myCurve->FirstParameter(),
916 myCurve->LastParameter(),
917 mySurface->FirstUParameter(),
918 mySurface->UPeriod());
919 */
920 //xt
921 // Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 Begin
922 // Correct the U isoline in periodical surface
923 // to be inside restriction boundaries.
924 if (myResult.GetType() == GeomAbs_Line) {
925 gp_Lin2d &aLine = (gp_Lin2d &) myResult.Line();
926
927 Standard_Real aPeriod = mySurface->UPeriod();
928 Standard_Real aFUPar = mySurface->FirstUParameter();
929 Standard_Real aLUPar = mySurface->LastUParameter();
930
931 // Check if the parametric range is lower then the period.
932 if (aLUPar - aFUPar < aPeriod - myTolerance) {
933 Standard_Real aU = aLine.Location().X();
934
935 if (Abs(aU + aPeriod - aFUPar) < myTolerance ||
936 Abs(aU - aPeriod - aFUPar) < myTolerance) {
937 gp_Pnt2d aNewLoc(aFUPar, aLine.Location().Y());
938
939 aLine.SetLocation(aNewLoc);
940 } else if (Abs(aU + aPeriod - aLUPar) < myTolerance ||
941 Abs(aU - aPeriod - aLUPar) < myTolerance) {
942 gp_Pnt2d aNewLoc(aLUPar, aLine.Location().Y());
943
944 aLine.SetLocation(aNewLoc);
945 }
946 }
947 }
948 }
949// Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 End
950
951 if (mySurface->IsVPeriodic()) {
952 myResult.VFrame(myCurve->FirstParameter(),
953 myCurve->LastParameter(),
954 mySurface->FirstVParameter(),
955 mySurface->VPeriod());
956// Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 Begin
957// Correct the V isoline in a periodical surface
958// to be inside restriction boundaries.
959 if (myResult.GetType() == GeomAbs_Line) {
960 gp_Lin2d &aLine = (gp_Lin2d &) myResult.Line();
961
962 Standard_Real aPeriod = mySurface->VPeriod();
963 Standard_Real aFVPar = mySurface->FirstVParameter();
964 Standard_Real aLVPar = mySurface->LastVParameter();
965
966 // Check if the parametric range is lower then the period.
967 if (aLVPar - aFVPar < aPeriod - myTolerance) {
968 Standard_Real aV = aLine.Location().Y();
969
970 if (Abs(aV + aPeriod - aFVPar) < myTolerance ||
971 Abs(aV - aPeriod - aFVPar) < myTolerance) {
972 gp_Pnt2d aNewLoc(aLine.Location().X(), aFVPar);
973
974 aLine.SetLocation(aNewLoc);
975 } else if (Abs(aV + aPeriod - aLVPar) < myTolerance ||
976 Abs(aV - aPeriod - aLVPar) < myTolerance) {
977 gp_Pnt2d aNewLoc(aLine.Location().X(), aLVPar);
978
979 aLine.SetLocation(aNewLoc);
980 }
981 }
982 }
983 }
984// Modified by skv - Wed Aug 11 15:45:58 2004 OCC6272 End
985 }
986}
987
988
989//=======================================================================
990//function : GetSurface
991//purpose :
992//=======================================================================
993
994const Handle(Adaptor3d_HSurface)& ProjLib_ProjectedCurve::GetSurface() const
995{
996 return mySurface;
997}
998
999
1000//=======================================================================
1001//function : GetCurve
1002//purpose :
1003//=======================================================================
1004
1005const Handle(Adaptor3d_HCurve)& ProjLib_ProjectedCurve::GetCurve() const
1006{
1007 return myCurve;
1008}
1009
1010
1011//=======================================================================
1012//function : GetTolerance
1013//purpose :
1014//=======================================================================
1015
1016Standard_Real ProjLib_ProjectedCurve::GetTolerance() const
1017{
1018 return myTolerance;
1019}
1020
1021
1022//=======================================================================
1023//function : FirstParameter
1024//purpose :
1025//=======================================================================
1026
1027Standard_Real ProjLib_ProjectedCurve::FirstParameter() const
1028{
1029 return myCurve->FirstParameter();
1030}
1031
1032
1033//=======================================================================
1034//function : LastParameter
1035//purpose :
1036//=======================================================================
1037
1038Standard_Real ProjLib_ProjectedCurve::LastParameter() const
1039{
1040 return myCurve->LastParameter();
1041}
1042
1043
1044//=======================================================================
1045//function : Continuity
1046//purpose :
1047//=======================================================================
1048
1049GeomAbs_Shape ProjLib_ProjectedCurve::Continuity() const
1050{
1051 Standard_NotImplemented::Raise("");
1052 return GeomAbs_C0;
1053}
1054
1055
1056//=======================================================================
1057//function : NbIntervals
1058//purpose :
1059//=======================================================================
1060
1061Standard_Integer ProjLib_ProjectedCurve::NbIntervals(const GeomAbs_Shape ) const
1062{
1063 Standard_NotImplemented::Raise("");
1064 return 0;
1065}
1066
1067
1068//=======================================================================
1069//function : Intervals
1070//purpose :
1071//=======================================================================
1072
1073//void ProjLib_ProjectedCurve::Intervals(TColStd_Array1OfReal& T,
1074void ProjLib_ProjectedCurve::Intervals(TColStd_Array1OfReal& ,
1075 const GeomAbs_Shape ) const
1076{
1077 Standard_NotImplemented::Raise("");
1078}
1079
1080
1081//=======================================================================
1082//function : IsClosed
1083//purpose :
1084//=======================================================================
1085
1086Standard_Boolean ProjLib_ProjectedCurve::IsClosed() const
1087{
1088 Standard_NotImplemented::Raise("");
1089 return Standard_True;
1090}
1091
1092
1093//=======================================================================
1094//function : IsPeriodic
1095//purpose :
1096//=======================================================================
1097
1098Standard_Boolean ProjLib_ProjectedCurve::IsPeriodic() const
1099{
1100 return myResult.IsPeriodic();
1101}
1102
1103
1104//=======================================================================
1105//function : Period
1106//purpose :
1107//=======================================================================
1108
1109Standard_Real ProjLib_ProjectedCurve::Period() const
1110{
1111 Standard_NotImplemented::Raise("");
1112 return 0.;
1113}
1114
1115
1116//=======================================================================
1117//function : Value
1118//purpose :
1119//=======================================================================
1120
1121gp_Pnt2d ProjLib_ProjectedCurve::Value(const Standard_Real ) const
1122{
1123 Standard_NotImplemented::Raise("");
1124 return gp_Pnt2d(0.,0.);
1125}
1126
1127
1128//=======================================================================
1129//function : D0
1130//purpose :
1131//=======================================================================
1132
1133void ProjLib_ProjectedCurve::D0(const Standard_Real , gp_Pnt2d& ) const
1134{
1135 Standard_NotImplemented::Raise("");
1136}
1137
1138
1139//=======================================================================
1140//function : D1
1141//purpose :
1142//=======================================================================
1143
1144void ProjLib_ProjectedCurve::D1(const Standard_Real ,
1145 gp_Pnt2d& ,
1146 gp_Vec2d& ) const
1147{
1148 Standard_NotImplemented::Raise("");
1149}
1150
1151
1152//=======================================================================
1153//function : D2
1154//purpose :
1155//=======================================================================
1156
1157void ProjLib_ProjectedCurve::D2(const Standard_Real ,
1158 gp_Pnt2d& ,
1159 gp_Vec2d& ,
1160 gp_Vec2d& ) const
1161{
1162 Standard_NotImplemented::Raise("");
1163}
1164
1165
1166//=======================================================================
1167//function : D3
1168//purpose :
1169//=======================================================================
1170
1171void ProjLib_ProjectedCurve::D3(const Standard_Real,
1172 gp_Pnt2d&,
1173 gp_Vec2d&,
1174 gp_Vec2d&,
1175 gp_Vec2d&) const
1176{
1177 Standard_NotImplemented::Raise("");
1178}
1179
1180
1181//=======================================================================
1182//function : DN
1183//purpose :
1184//=======================================================================
1185
1186gp_Vec2d ProjLib_ProjectedCurve::DN(const Standard_Real,
1187 const Standard_Integer) const
1188{
1189 Standard_NotImplemented::Raise("");
1190 return gp_Vec2d(0.,0.);
1191}
1192
1193
1194//=======================================================================
1195//function : Resolution
1196//purpose :
1197//=======================================================================
1198
1199Standard_Real ProjLib_ProjectedCurve::Resolution(const Standard_Real) const
1200{
1201 Standard_NotImplemented::Raise("");
1202 return 0.;
1203}
1204
1205
1206//=======================================================================
1207//function : GetType
1208//purpose :
1209//=======================================================================
1210
1211GeomAbs_CurveType ProjLib_ProjectedCurve::GetType() const
1212{
1213 return myResult.GetType();
1214}
1215
1216
1217//=======================================================================
1218//function : Line
1219//purpose :
1220//=======================================================================
1221
1222gp_Lin2d ProjLib_ProjectedCurve::Line() const
1223{
1224 return myResult.Line();
1225}
1226
1227
1228//=======================================================================
1229//function : Circle
1230//purpose :
1231//=======================================================================
1232
1233gp_Circ2d ProjLib_ProjectedCurve::Circle() const
1234{
1235 return myResult.Circle();
1236}
1237
1238
1239//=======================================================================
1240//function : Ellipse
1241//purpose :
1242//=======================================================================
1243
1244gp_Elips2d ProjLib_ProjectedCurve::Ellipse() const
1245{
1246 return myResult.Ellipse();
1247}
1248
1249
1250//=======================================================================
1251//function : Hyperbola
1252//purpose :
1253//=======================================================================
1254
1255gp_Hypr2d ProjLib_ProjectedCurve::Hyperbola() const
1256{
1257 return myResult.Hyperbola();
1258}
1259
1260
1261//=======================================================================
1262//function : Parabola
1263//purpose :
1264//=======================================================================
1265
1266gp_Parab2d ProjLib_ProjectedCurve::Parabola() const
1267{
1268 return myResult.Parabola();
1269}
1270
1271
1272
1273//=======================================================================
1274//function : Degree
1275//purpose :
1276//=======================================================================
1277
1278Standard_Integer ProjLib_ProjectedCurve::Degree() const
1279{
1280 Standard_NoSuchObject_Raise_if
1281 ( (GetType() != GeomAbs_BSplineCurve) &&
1282 (GetType() != GeomAbs_BezierCurve),
1283 "ProjLib_ProjectedCurve:Degree");
1284 if (GetType() == GeomAbs_BSplineCurve) {
1285 return myResult.BSpline()->Degree();
1286 }
1287 else if (GetType() == GeomAbs_BezierCurve) {
1288 return myResult.Bezier()->Degree();
1289 }
1290
1291 // portage WNT
1292 return 0;
1293}
1294
1295//=======================================================================
1296//function : IsRational
1297//purpose :
1298//=======================================================================
1299
1300Standard_Boolean ProjLib_ProjectedCurve::IsRational() const
1301{
1302 Standard_NoSuchObject_Raise_if
1303 ( (GetType() != GeomAbs_BSplineCurve) &&
1304 (GetType() != GeomAbs_BezierCurve),
1305 "ProjLib_ProjectedCurve:IsRational");
1306 if (GetType() == GeomAbs_BSplineCurve) {
1307 return myResult.BSpline()->IsRational();
1308 }
1309 else if (GetType() == GeomAbs_BezierCurve) {
1310 return myResult.Bezier()->IsRational();
1311 }
1312 // portage WNT
1313 return Standard_False;
1314}
1315
1316//=======================================================================
1317//function : NbPoles
1318//purpose :
1319//=======================================================================
1320
1321Standard_Integer ProjLib_ProjectedCurve::NbPoles() const
1322{
1323 Standard_NoSuchObject_Raise_if
1324 ( (GetType() != GeomAbs_BSplineCurve) &&
1325 (GetType() != GeomAbs_BezierCurve)
1326 ,"ProjLib_ProjectedCurve:NbPoles" );
1327 if (GetType() == GeomAbs_BSplineCurve) {
1328 return myResult.BSpline()->NbPoles();
1329 }
1330 else if (GetType() == GeomAbs_BezierCurve) {
1331 return myResult.Bezier()->NbPoles();
1332 }
1333
1334 // portage WNT
1335 return 0;
1336}
1337
1338//=======================================================================
1339//function : NbKnots
1340//purpose :
1341//=======================================================================
1342
1343Standard_Integer ProjLib_ProjectedCurve::NbKnots() const
1344{
1345 Standard_NoSuchObject_Raise_if ( GetType() != GeomAbs_BSplineCurve,
1346 "ProjLib_ProjectedCurve:NbKnots");
1347 return myResult.BSpline()->NbKnots();
1348}
1349
1350//=======================================================================
1351//function : Bezier
1352//purpose :
1353//=======================================================================
1354
1355Handle(Geom2d_BezierCurve) ProjLib_ProjectedCurve::Bezier() const
1356{
1357 return myResult.Bezier() ;
1358}
1359
1360//=======================================================================
1361//function : BSpline
1362//purpose :
1363//=======================================================================
1364
1365Handle(Geom2d_BSplineCurve) ProjLib_ProjectedCurve::BSpline() const
1366{
1367 return myResult.BSpline() ;
1368}
1369//=======================================================================
1370//function : Trim
1371//purpose :
1372//=======================================================================
1373
1374Handle(Adaptor2d_HCurve2d) ProjLib_ProjectedCurve::Trim
1375//(const Standard_Real First,
1376// const Standard_Real Last,
1377// const Standard_Real Tolerance) const
1378(const Standard_Real ,
1379 const Standard_Real ,
1380 const Standard_Real ) const
1381{
1382 Standard_NotImplemented::Raise("");
1383 return NULL ;
1384}
1385