0024947: Redesign OCCT legacy type system
[occt.git] / src / TopOpeBRepTool / TopOpeBRepTool_CurveTool.cxx
CommitLineData
b311480e 1// Created on: 1993-06-24
2// Created by: Jean Yves LEBEY
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.
7fd59977 16
17#include <TopOpeBRepTool_CurveTool.ixx>
18
19#include <ProjLib_ProjectedCurve.hxx>
20#include <BRepApprox_Approx.hxx>
21#include <BRepApprox_ApproxLine.hxx>
22#include <TopoDS.hxx>
23#include <TopoDS_Face.hxx>
24#include <Geom_BSplineCurve.hxx>
25#include <Geom2d_Line.hxx>
26#include <Geom2d_Circle.hxx>
27#include <Geom2d_Ellipse.hxx>
28#include <Geom2d_Parabola.hxx>
29#include <Geom2d_Hyperbola.hxx>
30#include <Geom2d_BSplineCurve.hxx>
31#include <TColStd_Array1OfReal.hxx>
32#include <GeomLib_CheckBSplineCurve.hxx>
33#include <GeomLib_Check2dBSplineCurve.hxx>
34#include <TColStd_Array1OfInteger.hxx>
35#include <BRepAdaptor_Surface.hxx>
36#include <GeomAdaptor_Curve.hxx>
37#include <BRepAdaptor_HSurface.hxx>
38#include <GeomAdaptor_HCurve.hxx>
39#include <GeomAbs_SurfaceType.hxx>
40#include <gp_Pln.hxx>
41#include <gp_Pnt2d.hxx>
42#include <ElSLib.hxx>
43#include <Standard_NotImplemented.hxx>
44#include <Standard_ProgramError.hxx>
45#include <OSD_Chronometer.hxx>
46
47#include <GeomTools_Curve2dSet.hxx>
48#include <TColStd_Array1OfBoolean.hxx>
49#include <gp_Lin.hxx>
50#include <gce_MakeLin.hxx>
51#include <gp_Lin2d.hxx>
52#include <gce_MakeLin2d.hxx>
53#include <BRepTools.hxx>
54#include <Precision.hxx>
55#include <gce_MakeCirc.hxx>
56#include <gp_Vec.hxx>
57//#include <Approx.hxx>
58#include <Geom_Surface.hxx>
59#include <TopLoc_Location.hxx>
60#include <Geom_RectangularTrimmedSurface.hxx>
61#include <BRep_Tool.hxx>
62
0797d9d3 63#ifdef OCCT_DEBUG
7fd59977 64#include <TopOpeBRepTool_KRO.hxx>
65TOPKRO KRO_CURVETOOL_APPRO("approximation");
1d0a9d4d 66extern Standard_Boolean TopOpeBRepTool_GettraceKRO();
67extern Standard_Boolean TopOpeBRepTool_GettracePCURV();
68extern Standard_Boolean TopOpeBRepTool_GettraceCHKBSPL();
7fd59977 69#endif
70//#define DRAW
71//#define IFV
72#define CurveImprovement
73#ifdef DRAW
74#include <DrawTrSurf.hxx>
75static Standard_Integer NbCalls = 0;
76#endif
77//=======================================================================
78//function : CurveTool
79//purpose :
80//=======================================================================
81
82TopOpeBRepTool_CurveTool::TopOpeBRepTool_CurveTool()
83{}
84
85//=======================================================================
86//function : CurveTool
87//purpose :
88//=======================================================================
89
90TopOpeBRepTool_CurveTool::TopOpeBRepTool_CurveTool
91(const TopOpeBRepTool_OutCurveType O)
92{
93 TopOpeBRepTool_GeomTool GT(O);
94 SetGeomTool(GT);
95}
96
97//=======================================================================
98//function : CurveTool
99//purpose :
100//=======================================================================
101
102TopOpeBRepTool_CurveTool::TopOpeBRepTool_CurveTool
103(const TopOpeBRepTool_GeomTool& GT)
104{
105 SetGeomTool(GT);
106}
107
108//=======================================================================
109//function : ChangeGeomTool
110//purpose :
111//=======================================================================
112
113TopOpeBRepTool_GeomTool& TopOpeBRepTool_CurveTool::ChangeGeomTool()
114{
115 return myGeomTool;
116}
117
118//=======================================================================
119//function : GetGeomTool
120//purpose :
121//=======================================================================
122
123const TopOpeBRepTool_GeomTool& TopOpeBRepTool_CurveTool::GetGeomTool()const
124{
125 return myGeomTool;
126}
127
128//=======================================================================
129//function : SetGeomTool
130//purpose :
131//=======================================================================
132
133void TopOpeBRepTool_CurveTool::SetGeomTool
134(const TopOpeBRepTool_GeomTool& GT)
135{
136 myGeomTool.Define(GT);
137}
138
139//-----------------------------------------------------------------------
140//function : MakePCurve
141//purpose :
142//-----------------------------------------------------------------------
143Standard_EXPORT Handle(Geom2d_Curve) MakePCurve(const ProjLib_ProjectedCurve& PC)
144{
145 Handle(Geom2d_Curve) C2D;
146 switch (PC.GetType()) {
147 case GeomAbs_Line : C2D = new Geom2d_Line(PC.Line()); break;
148 case GeomAbs_Circle : C2D = new Geom2d_Circle(PC.Circle()); break;
149 case GeomAbs_Ellipse : C2D = new Geom2d_Ellipse(PC.Ellipse()); break;
150 case GeomAbs_Parabola : C2D = new Geom2d_Parabola(PC.Parabola()); break;
151 case GeomAbs_Hyperbola : C2D = new Geom2d_Hyperbola(PC.Hyperbola()); break;
152 case GeomAbs_BSplineCurve : C2D = PC.BSpline(); break;
153 case GeomAbs_BezierCurve : case GeomAbs_OtherCurve : default :
154 Standard_NotImplemented::Raise("CurveTool::MakePCurve");
155 break;
156 }
157 return C2D;
158}
159
160//------------------------------------------------------------------
161static Standard_Boolean CheckApproxResults
162 (const BRepApprox_Approx& Approx)
163//------------------------------------------------------------------
164{
165 const AppParCurves_MultiBSpCurve& amc = Approx.Value(1);
166 Standard_Integer np = amc.NbPoles();
167 Standard_Integer nc = amc.NbCurves();
168 if (np < 2 || nc < 1) return Standard_False;
169
170 // check the knots for coincidence
171 const TColStd_Array1OfReal& knots = amc.Knots();
172 for (Standard_Integer i = knots.Lower(); i < knots.Upper(); i++) {
173 if (knots(i+1) - knots(i) <= Epsilon(knots(i))) {
174 return Standard_False;
175 }
176 }
177 return Standard_True;
178}
179
180//------------------------------------------------------------------
181static Standard_Boolean CheckPCurve
182 (const Handle(Geom2d_Curve)& aPC, const TopoDS_Face& aFace)
183//------------------------------------------------------------------
184// check if points of the pcurve are out of the face bounds
185{
186 const Standard_Integer NPoints = 23;
187 Standard_Real umin,umax,vmin,vmax;
188
189 BRepTools::UVBounds(aFace, umin, umax, vmin, vmax);
190 Standard_Real tolU = Max ((umax-umin)*0.01, Precision::Confusion());
191 Standard_Real tolV = Max ((vmax-vmin)*0.01, Precision::Confusion());
192 Standard_Real fp = aPC->FirstParameter();
193 Standard_Real lp = aPC->LastParameter();
194 Standard_Real step = (lp-fp)/(NPoints+1);
195
196 // adjust domain for periodic surfaces
197 TopLoc_Location aLoc;
198 Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace, aLoc);
199 if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
200 aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
201
202 gp_Pnt2d pnt = aPC->Value((fp+lp)/2);
203 Standard_Real u,v;
204 pnt.Coord(u,v);
205
206 if (aSurf->IsUPeriodic()) {
207 Standard_Real aPer = aSurf->UPeriod();
208 Standard_Integer nshift = (Standard_Integer) ((u-umin)/aPer);
209 if (u < umin+aPer*nshift) nshift--;
210 umin += aPer*nshift;
211 umax += aPer*nshift;
212 }
213 if (aSurf->IsVPeriodic()) {
214 Standard_Real aPer = aSurf->VPeriod();
215 Standard_Integer nshift = (Standard_Integer) ((v-vmin)/aPer);
216 if (v < vmin+aPer*nshift) nshift--;
217 vmin += aPer*nshift;
218 vmax += aPer*nshift;
219 }
220
221 Standard_Integer i;
222 for (i=1; i <= NPoints; i++) {
223 Standard_Real p = fp + i * step;
224 pnt = aPC->Value(p);
225 pnt.Coord(u,v);
226 if (umin-u > tolU || u-umax > tolU ||
227 vmin-v > tolV || v-vmax > tolV)
228 return Standard_False;
229 }
230 return Standard_True;
231}
232
233//------------------------------------------------------------------
234static Handle(Geom_Curve) MakeCurve3DfromWLineApprox
235(const BRepApprox_Approx& Approx,const Standard_Integer )
236//------------------------------------------------------------------
237{
238 const AppParCurves_MultiBSpCurve& amc = Approx.Value(1);
239 Standard_Integer np = amc.NbPoles();
240 //Standard_Integer nc = amc.NbCurves();
241 TColgp_Array1OfPnt poles3d(1,np);
242 Standard_Integer ic = 1;
243 //for (ic=1; ic<=nc; ic++) {
244 //if (ic == CI) {
245 amc.Curve(ic,poles3d);
246 //}
247 //}
248
249 const TColStd_Array1OfReal& knots = amc.Knots();
250 const TColStd_Array1OfInteger& mults = amc.Multiplicities();
251 Standard_Integer degree = amc.Degree();
252 Handle(Geom_Curve) C3D = new Geom_BSplineCurve(poles3d,knots,mults,degree);
253 return C3D;
254}
255
256//------------------------------------------------------------------
257static Handle(Geom2d_Curve) MakeCurve2DfromWLineApproxAndPlane
258(const BRepApprox_Approx& Approx,const gp_Pln& Pl)
259//------------------------------------------------------------------
260{
261 const AppParCurves_MultiBSpCurve& amc = Approx.Value(1);
262 Standard_Integer np = amc.NbPoles();
263 TColgp_Array1OfPnt2d poles2d(1,np);
264 TColgp_Array1OfPnt poles3d(1,np);
265 amc.Curve(1,poles3d);
266 for(Standard_Integer i=1; i<=np; i++) {
267 Standard_Real U,V; ElSLib::Parameters(Pl,poles3d.Value(i),U,V);
268 poles2d.SetValue(i,gp_Pnt2d(U,V));
269 }
270 const TColStd_Array1OfReal& knots = amc.Knots();
271 const TColStd_Array1OfInteger& mults = amc.Multiplicities();
272 Standard_Integer degree = amc.Degree();
273 Handle(Geom2d_Curve) C2D = new Geom2d_BSplineCurve(poles2d,knots,mults,degree);
274 return C2D;
275}
276
277//------------------------------------------------------------------
278static Handle(Geom2d_Curve) MakeCurve2DfromWLineApprox
279(const BRepApprox_Approx& Approx,const Standard_Integer CI)
280//------------------------------------------------------------------
281{
282 const AppParCurves_MultiBSpCurve& amc = Approx.Value(1);
283 Standard_Integer np = amc.NbPoles();
284 TColgp_Array1OfPnt2d poles2d(1,np);
285 Standard_Integer nc = amc.NbCurves();
286 for (Standard_Integer ic=1; ic<=nc; ic++) if (ic == CI) amc.Curve(ic,poles2d);
287 const TColStd_Array1OfReal& knots = amc.Knots();
288 const TColStd_Array1OfInteger& mults = amc.Multiplicities();
289 Standard_Integer degree = amc.Degree();
290 Handle(Geom2d_Curve) C2D = new Geom2d_BSplineCurve(poles2d,knots,mults,degree);
291 return C2D;
292}
293
294
295//=======================================================================
296//function : MakeCurves
297//purpose :
298//=======================================================================
299
300Standard_Boolean TopOpeBRepTool_CurveTool::MakeCurves
301(const Standard_Real parmin, const Standard_Real parmax,
302 const Handle(Geom_Curve)& C3D,
303 const Handle(Geom2d_Curve)& PC1, const Handle(Geom2d_Curve)& PC2,
304 const TopoDS_Shape& S1, const TopoDS_Shape& S2,
305 Handle(Geom_Curve)& C3Dnew,
306 Handle(Geom2d_Curve)& PC1new, Handle(Geom2d_Curve)& PC2new,
307 Standard_Real& TolReached3d, Standard_Real& TolReached2d) const
308{
309 Standard_Boolean CompC3D = myGeomTool.CompC3D();
310
311 //cout << "MakeCurves begin" << endl;
312 if (!CompC3D) return Standard_False;
313
314 Standard_Boolean CompPC1 = myGeomTool.CompPC1();
315 Standard_Boolean CompPC2 = myGeomTool.CompPC2();
316 Standard_Real tol3d,tol2d;
317 Standard_Boolean RelativeTol;
318 myGeomTool.GetTolerances(tol3d,tol2d,RelativeTol);
319 Standard_Integer NbPntMax = myGeomTool.NbPntMax();
320
0797d9d3 321#ifdef OCCT_DEBUG
7fd59977 322 if (TopOpeBRepTool_GettraceKRO()) KRO_CURVETOOL_APPRO.Start();
323#endif
324
325//----------------------------------
326///*
327#ifdef IFV
328 char name[16];
329 char *nm = &name[0];
330 sprintf(name,"C3D_%d",++NbCalls);
331 DrawTrSurf::Set(nm, C3D);
332 sprintf(name,"PC1_%d",NbCalls);
333 DrawTrSurf::Set(nm, PC1);
334 sprintf(name,"PC2_%d",NbCalls);
335 DrawTrSurf::Set(nm, PC2);
336#endif
337//*/
338//---------------------------------------------
339
340 Standard_Integer iparmin = (Standard_Integer)parmin;
341 Standard_Integer iparmax = (Standard_Integer)parmax;
342
343 Handle(Geom_BSplineCurve)& HC3D = *((Handle(Geom_BSplineCurve)*)&C3D);
344 Handle(Geom2d_BSplineCurve)& HPC1 = *((Handle(Geom2d_BSplineCurve)*)&PC1);
345 Handle(Geom2d_BSplineCurve)& HPC2 = *((Handle(Geom2d_BSplineCurve)*)&PC2);
346
347//--------------------- IFV - "improving" initial curves
348#ifdef CurveImprovement
349 Standard_Integer nbpol = HC3D->NbPoles();
350 //cout <<"nbpol = " << nbpol << endl;
351 if(nbpol > 100) {
352 TColgp_Array1OfPnt PolC3D(1, nbpol);
353 TColgp_Array1OfPnt2d PolPC1(1, nbpol);
354 TColgp_Array1OfPnt2d PolPC2(1, nbpol);
355 TColStd_Array1OfBoolean IsValid(1, nbpol);
356 IsValid.Init(Standard_True);
357 Standard_Real tol = Max(1.e-10, 100.*tol3d*tol3d); //tol *= tol; - square distance
358 Standard_Real tl2d = tol*(tol2d*tol2d)/(tol3d*tol3d);
359 Standard_Real def = tol;
360 Standard_Real def2d = tol2d;
361 HC3D->Poles(PolC3D);
362 if(CompPC1) HPC1->Poles(PolPC1);
363 if(CompPC2) HPC2->Poles(PolPC2);
364
365 Standard_Integer ip = 1, NbPol = 1;
366 Standard_Real d, d1, d2;
367 gp_Pnt P = PolC3D(ip);
368 gp_Pnt2d P1, P2;
369 if(CompPC1) P1 = PolPC1(ip);
370 if(CompPC2) P2 = PolPC2(ip);
371
372 for(ip = 2; ip <= nbpol; ip++) {
373 if( IsValid(ip) ) {
374 d = P.SquareDistance(PolC3D(ip));
375 if(CompPC1) {d1 = P1.SquareDistance(PolPC1(ip));} else {d1 = 0.;}
376 if(CompPC2) {d2 = P2.SquareDistance(PolPC2(ip));} else {d2 = 0.;}
377 if(d > tol || d1 > tl2d || d2 > tl2d ) {
378 Standard_Real dd=0.;
379 if(ip < nbpol ) dd = P.Distance(PolC3D(ip+1));
380 if(ip < nbpol && dd < 10.*tol) {
381 gce_MakeLin mkL(P, PolC3D(ip+1));
382 if(mkL.IsDone()) {
383 gp_Lin L = mkL.Value();
384 d = L.SquareDistance(PolC3D(ip));
385 if(CompPC1) {
386 gp_Lin2d L1 = gce_MakeLin2d(P1, PolPC1(ip+1));
387 d1 = L1.SquareDistance(PolPC1(ip));
388 }
389 else d1 = 0.;
390 if(CompPC2) {
391 gp_Lin2d L2 = gce_MakeLin2d(P2, PolPC2(ip+1));
392 d2 = L2.SquareDistance(PolPC2(ip));
393 }
394 else d2 = 0.;
395
396 if(d > def || d1 > def2d || d2 > def2d ) {
397 NbPol++;
398 P = PolC3D(ip);
399 if(CompPC1) P1 = PolPC1(ip);
400 if(CompPC2) P2 = PolPC2(ip);
401 }
402 else {
403 IsValid(ip) = Standard_False;
404 }
405 }
406 else {
407 IsValid(ip+1) = Standard_False;
408 NbPol++;
409 P = PolC3D(ip);
410 if(CompPC1) P1 = PolPC1(ip);
411 if(CompPC2) P2 = PolPC2(ip);
412 }
413 }
414 else {
415 NbPol++;
416 P = PolC3D(ip);
417 if(CompPC1) P1 = PolPC1(ip);
418 if(CompPC2) P2 = PolPC2(ip);
419 }
420 }
421 else {
422 IsValid(ip) = Standard_False;
423 }
424 }
425 }
426
427 if(NbPol < 2) {IsValid(nbpol) = Standard_True; NbPol++;}
428
429 if(NbPol < nbpol) {
430#ifdef IFV
431 cout << "NbPol < nbpol : " << NbPol << " " << nbpol << endl;
432#endif
433 TColgp_Array1OfPnt Polc3d(1, NbPol);
434 TColgp_Array1OfPnt2d Polpc1(1, NbPol);
435 TColgp_Array1OfPnt2d Polpc2(1, NbPol);
436 TColStd_Array1OfReal knots(1,NbPol);
437 TColStd_Array1OfInteger mults(1,NbPol);
438 mults.Init(1); mults(1) = 2; mults(NbPol) = 2;
439 Standard_Integer count = 0;
440 for(ip = 1; ip <= nbpol; ip++) {
441 if(IsValid(ip)) {
442 count++;
443 Polc3d(count) = PolC3D(ip);
444 if(CompPC1) Polpc1(count) = PolPC1(ip);
445 if(CompPC2) Polpc2(count) = PolPC2(ip);
446 knots(count) = count;
447 }
448 }
449
450 Polc3d(NbPol) = PolC3D(nbpol);
451 if(CompPC1) Polpc1(NbPol) = PolPC1(nbpol);
452 if(CompPC2) Polpc2(NbPol) = PolPC2(nbpol);
453
454 HC3D = new Geom_BSplineCurve(Polc3d, knots, mults, 1);
455 if(CompPC1) HPC1 = new Geom2d_BSplineCurve(Polpc1, knots, mults, 1);
456 if(CompPC2) HPC2 = new Geom2d_BSplineCurve(Polpc2, knots, mults, 1);
457 iparmax = NbPol;
458
459#ifdef IFV
460 sprintf(name,"C3Dmod_%d",NbCalls);
461 nm = &name[0];
462 DrawTrSurf::Set(nm, HC3D);
463 sprintf(name,"PC1mod_%d",NbCalls);
464 nm = &name[0];
465 DrawTrSurf::Set(nm, HPC1);
466 sprintf(name,"PC2mod_%d",NbCalls);
467 nm = &name[0];
468 DrawTrSurf::Set(nm, HPC2);
469#endif
470
471 }
472 }
473//--------------- IFV - end "improving"
474#endif
475
476
477 BRepApprox_Approx Approx;
478
479 Standard_Integer degmin = 4;
480 Standard_Integer degmax = 8;
481 Approx_ParametrizationType parametrization = Approx_ChordLength;
482
483 Standard_Integer npol = HC3D->NbPoles();
484 TColgp_Array1OfPnt Polc3d(1, npol);
485 TColStd_Array1OfReal par(1, npol);
486 HC3D->Poles(Polc3d);
487 gp_Pnt P = Polc3d(1);
488
489 Standard_Boolean IsBad = Standard_False;
490 Standard_Integer ip;
491 Standard_Real ksi = 0., kc, kcprev = 0.;
492 Standard_Real dist;
493 par(1) = 0.;
494 for(ip = 2; ip <= npol; ip++) {
495 dist = P.Distance(Polc3d(ip));
496
497 if(dist < Precision::Confusion()) {
498 IsBad = Standard_True;
499 break;
500 }
501
502 par(ip) = par(ip-1) + dist;
503 P = Polc3d(ip);
504 }
505
506 if(!IsBad) {
507
508 TColStd_Array1OfReal Curvature(1, npol);
509
510 if(npol > 3) {
511 Standard_Integer ic = 1;
512 for(ip = 2; ip <= npol-1; ip += npol-3) {
513 gp_Vec v1(Polc3d(ip-1),Polc3d(ip));
514 gp_Vec v2(Polc3d(ip),Polc3d(ip+1));
515 if(!v1.IsParallel(v2, 1.e-4)) {
516 gce_MakeCirc mc(Polc3d(ip-1),Polc3d(ip),Polc3d(ip+1));
517 gp_Circ cir = mc.Value();
518 kc = 1./cir.Radius();
519 }
520 else kc = 0.;
521
522 Curvature(ic) = kc;
523 ic = npol;
524 }
525 }
526 else if(npol == 3) {
527 ip = 2;
528 gp_Vec v1(Polc3d(ip-1),Polc3d(ip));
529 gp_Vec v2(Polc3d(ip),Polc3d(ip+1));
530 if(!v1.IsParallel(v2, 1.e-4)) {
531 gce_MakeCirc mc(Polc3d(ip-1),Polc3d(ip),Polc3d(ip+1));
532 gp_Circ cir = mc.Value();
533 kc = 1./cir.Radius();
534 }
535 else kc = 0.;
536 Curvature(1) = Curvature(npol) = kc;
537 }
538 else {
539 Curvature(1) = Curvature(npol) = 0.;
540 }
541
542 ip = 1;
543 Standard_Real dt = par(ip+1) - par(ip);
544 Standard_Real dx = (Polc3d(ip+1).X() - Polc3d(ip).X())/dt,
545 dy = (Polc3d(ip+1).Y() - Polc3d(ip).Y())/dt,
546 dz = (Polc3d(ip+1).Z() - Polc3d(ip).Z())/dt;
547 Standard_Real dx1, dy1, dz1, d2x, d2y, d2z, d2t;
548
549 for(ip = 2; ip <= npol-1; ip++) {
550 dt = par(ip+1) - par(ip);
551 dx1 = (Polc3d(ip+1).X() - Polc3d(ip).X())/dt;
552 dy1 = (Polc3d(ip+1).Y() - Polc3d(ip).Y())/dt,
553 dz1 = (Polc3d(ip+1).Z() - Polc3d(ip).Z())/dt;
554 d2t = 2./(par(ip+1) - par(ip-1));
555 d2x = d2t*(dx1 - dx);
556 d2y = d2t*(dy1 - dy);
557 d2z = d2t*(dz1 - dz);
558 Curvature(ip) = Sqrt(d2x*d2x + d2y*d2y + d2z*d2z);
559 dx = dx1; dy = dy1; dz = dz1;
560 }
561
562 Standard_Real crit = 1000.; // empirical criterion !!!
563
564 dt = par(2) - par(1);
565 kcprev = (Curvature(2) - Curvature(1))/dt;
566 for(ip = 2; ip <= npol-1; ip++) {
567 dt = par(ip+1) - par(ip);
568 kc = (Curvature(ip+1) - Curvature(ip))/dt;
569 ksi = ksi + Abs(kc - kcprev);
570 if(ksi > crit) {IsBad = Standard_True;break;}
571 kc = kcprev;
572 }
573
574 }
575 //cout << NbCalls << " ksi = " << ksi << endl;
576 //cout << "IsBad = " << IsBad << endl;
577
578 if(IsBad){
579 Standard_Real tt = Min(10.*tol3d,Precision::Approximation());
580 tol2d = tt * tol2d / tol3d;
581 tol3d = tt;
582 NbPntMax = 40;
583 degmax = 4;
584 parametrization = Approx_Centripetal;
585 }
586
587 Standard_Integer nitmax = 0; // use projection only
588 Standard_Boolean withtangency = Standard_True;
589
590 Standard_Boolean compminmaxUV = Standard_True;
591 BRepAdaptor_Surface BAS1(TopoDS::Face(S1),compminmaxUV);
592 BRepAdaptor_Surface BAS2(TopoDS::Face(S2),compminmaxUV);
593
594
595 Handle(BRepApprox_ApproxLine) AL;
596 AL = new BRepApprox_ApproxLine(HC3D,HPC1,HPC2);
597
598 Approx.SetParameters(tol3d,tol2d,RelativeTol,degmin,degmax,nitmax,NbPntMax,withtangency,
599 parametrization);
600
601 if (CompC3D && CompPC1 && BAS1.GetType() == GeomAbs_Plane) {
602 //-- The curve X,Y,Z and U2,V2 is approximated
603 Approx.Perform(BAS1,BAS2,AL,CompC3D,Standard_False,CompPC2,iparmin,iparmax);
604 }
605
606 else if(CompC3D && CompPC2 && BAS2.GetType() == GeomAbs_Plane) {
607 //-- The curve X,Y,Z and U1,V1 is approximated
608 Approx.Perform(BAS1,BAS2,AL,CompC3D,CompPC1,Standard_False,iparmin,iparmax);
609 }
610
611 else {
612 Approx.Perform(BAS1,BAS2,AL,CompC3D,CompPC1,CompPC2,iparmin,iparmax);
613 }
614
615 // MSV Nov 9, 2001: do not raise exception in the case of failure,
616 // but return status
617
618 Standard_Boolean done = Approx.IsDone();
619 done = done && CheckApproxResults(Approx);
620
621 if (done) {
622 if (CompC3D && CompPC1 && BAS1.GetType() == GeomAbs_Plane) {
623 C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
624 PC1new = ::MakeCurve2DfromWLineApproxAndPlane(Approx,BAS1.Plane());
625 if (CompPC2) PC2new = ::MakeCurve2DfromWLineApprox(Approx,2);
626 }
627 else if(CompC3D && CompPC2 && BAS2.GetType() == GeomAbs_Plane) {
628 C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
629 if (CompPC1) PC1new = ::MakeCurve2DfromWLineApprox(Approx,2);
630 PC2new = ::MakeCurve2DfromWLineApproxAndPlane(Approx,BAS2.Plane());
631 }
632 else {
633 if (CompC3D) C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
634 if (CompPC1) PC1new = ::MakeCurve2DfromWLineApprox(Approx,2);
635 if (CompPC2) {
636 Standard_Integer i32 = (CompPC1) ? 3 : 2;
637 PC2new = ::MakeCurve2DfromWLineApprox(Approx,i32);
638 }
639 }
640
641 // check the pcurves relatively the faces bounds
642 if (CompPC1)
643 done = done && CheckPCurve(PC1new,TopoDS::Face(S1));
644 if (CompPC2)
645 done = done && CheckPCurve(PC2new,TopoDS::Face(S2));
646 }
647
648 if (!done) {
649 if (CompC3D) C3Dnew.Nullify();
650 if (CompPC1) PC1new.Nullify();
651 if (CompPC2) PC2new.Nullify();
652 return Standard_False;
653 }
654
655#ifdef IFV
656 sprintf(name,"C3Dnew_%d", NbCalls);
657 nm = &name[0];
658 DrawTrSurf::Set(nm, C3Dnew);
659 if (CompPC1) {
660 sprintf(name,"PC1new_%d", NbCalls);
661 nm = &name[0];
662 DrawTrSurf::Set(nm, PC1new);
663 }
664 if (CompPC2) {
665 sprintf(name,"PC2new_%d", NbCalls);
666 nm = &name[0];
667 DrawTrSurf::Set(nm, PC2new);
668 }
669
670#endif
671
672 TolReached3d = Approx.TolReached3d();
673 TolReached2d = Approx.TolReached2d();
674#ifdef IFV
675 cout << NbCalls << " : Tol = " << TolReached3d << " " << TolReached2d << endl;
676#endif
677
678 Standard_Boolean bf, bl;
679
680 const Handle(Geom_BSplineCurve)& Curve = Handle(Geom_BSplineCurve)::DownCast(C3Dnew);
681 if(!Curve.IsNull()) {
682 GeomLib_CheckBSplineCurve cbsc(Curve, 1.e-7, 0.1);
683 cbsc.NeedTangentFix(bf, bl);
684
0797d9d3 685#ifdef OCCT_DEBUG
7fd59977 686 if (TopOpeBRepTool_GettraceCHKBSPL()) {
687 if(bf || bl) {
688 cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
689 cout<<" Last = "<<bl<<endl;
690 }
691 }
692#endif
693 cbsc.FixTangent(bf, bl);
694 }
695
696 const Handle(Geom2d_BSplineCurve)& Curve2df = Handle(Geom2d_BSplineCurve)::DownCast(PC1new);
697 if(!Curve2df.IsNull()) {
698 GeomLib_Check2dBSplineCurve cbsc2df(Curve2df, 1.e-7, 0.1);
699 cbsc2df.NeedTangentFix(bf, bl);
0797d9d3 700#ifdef OCCT_DEBUG
7fd59977 701 if (TopOpeBRepTool_GettraceCHKBSPL()) {
702 if(bf || bl) {
703 cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
704 cout<<" Last = "<<bl<<endl;
705 }
706 }
707#endif
708 cbsc2df.FixTangent(bf, bl);
709 }
710
711 const Handle(Geom2d_BSplineCurve)& Curve2ds = Handle(Geom2d_BSplineCurve)::DownCast(PC2new);
712 if(!Curve2ds.IsNull()) {
713 GeomLib_Check2dBSplineCurve cbsc2ds(Curve2ds, 1.e-7, 0.1);
714 cbsc2ds.NeedTangentFix(bf, bl);
0797d9d3 715#ifdef OCCT_DEBUG
7fd59977 716 if (TopOpeBRepTool_GettraceCHKBSPL()) {
717 if(bf || bl) {
718 cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
719 cout<<" Last = "<<bl<<endl;
720 }
721 }
722#endif
723 cbsc2ds.FixTangent(bf, bl);
724 }
725
0797d9d3 726#ifdef OCCT_DEBUG
7fd59977 727 if (TopOpeBRepTool_GettraceKRO()) KRO_CURVETOOL_APPRO.Stop();
728#endif
729// cout << "MakeCurves end" << endl;
730
731 return Standard_True;
732}
733
734
735//=======================================================================
736//function : MakeBSpline1fromPnt
737//purpose :
738//=======================================================================
739
740Handle(Geom_Curve) TopOpeBRepTool_CurveTool::MakeBSpline1fromPnt
741(const TColgp_Array1OfPnt& Points)
742{
743 Standard_Integer Degree = 1;
744
745 Standard_Integer i,nbpoints = Points.Length();
746 Standard_Integer nbknots = nbpoints - Degree +1;
747
748 // First compute the parameters
749 // Standard_Real length = 0.;
750 // TColStd_Array1OfReal parameters(1,nbpoints);
751 // for (i = 1; i < nbpoints; i++) {
752 // parameters(i) = length;
753 // Standard_Real dist = Points(i).Distance(Points(i+1));
754 // length += dist;
755 // }
756 // parameters(nbpoints) = length;
757
758 // knots and multiplicities
759 TColStd_Array1OfReal knots(1,nbknots);
760 TColStd_Array1OfInteger mults(1,nbknots);
761 mults.Init(1);
762 mults(1) = mults(nbknots) = Degree + 1;
763
764 // knots(1) = 0;
765 // for (i=2;i<nbknots;i++) knots(i) = (parameters(i) + parameters(i+1)) /2.;
766 // knots(nbknots) = length;
767
768 // take point index as parameter : JYL 01/AUG/94
769 for (i = 1; i <= nbknots; i++) knots(i) = (Standard_Real)i;
770 Handle(Geom_Curve) C = new Geom_BSplineCurve(Points,knots,mults,Degree);
771 return C;
772}
773
774
775//=======================================================================
776//function : MakeBSpline1fromPnt2d
777//purpose :
778//=======================================================================
779
780Handle(Geom2d_Curve) TopOpeBRepTool_CurveTool::MakeBSpline1fromPnt2d
781(const TColgp_Array1OfPnt2d& Points)
782{
783 Standard_Integer Degree = 1;
784
785 Standard_Integer i,nbpoints = Points.Length();
786 Standard_Integer nbknots = nbpoints - Degree +1;
787
788 // First compute the parameters
789 // Standard_Real length = 0;
790 // TColStd_Array1OfReal parameters(1,nbpoints);
791 // for (i = 1; i < nbpoints; i++) {
792 // parameters(i) = length;
793 // Standard_Real dist = Points(i).Distance(Points(i+1));
794 // length += dist;
795 // }
796 // parameters(nbpoints) = length;
797
798 // knots and multiplicities
799 TColStd_Array1OfReal knots(1,nbknots);
800 TColStd_Array1OfInteger mults(1,nbknots);
801 mults.Init(1);
802 mults(1) = mults(nbknots) = Degree + 1;
803
804 // knots(1) = 0;
805 // for (i=2;i<nbknots;i++) knots(i) = (parameters(i) + parameters(i+1)) /2.;
806 // knots(nbknots) = length;
807
808 // take point index as parameter : JYL 01/AUG/94
809 for (i = 1; i <= nbknots; i++) knots(i) = (Standard_Real)i;
810 Handle(Geom2d_Curve) C = new Geom2d_BSplineCurve(Points,knots,mults,Degree);
811 return C;
812}
813
814//=======================================================================
815//function : IsProjectable
816//purpose :
817//=======================================================================
818
819Standard_Boolean TopOpeBRepTool_CurveTool::IsProjectable
820(const TopoDS_Shape& S, const Handle(Geom_Curve)& C3D)
821{
822 const TopoDS_Face& F = TopoDS::Face(S);
823 Standard_Boolean compminmaxUV = Standard_False;
824 BRepAdaptor_Surface BAS(F,compminmaxUV);
825 GeomAbs_SurfaceType suty = BAS.GetType();
826 GeomAdaptor_Curve GAC(C3D);
827 GeomAbs_CurveType cuty = GAC.GetType();
828
829 // --------
830 // avoid projection of 3d curve on surface in case
831 // of a quadric (ellipse,hyperbola,parabola) on a cone.
832 // Projection fails when the curve in not fully inside the UV domain
833 // of the cone : only part of 2d curve is built.
834 // NYI : projection of quadric on cone (crossing cone domain)
835 // --------
836
837 Standard_Boolean projectable = Standard_True;
838 if ( suty == GeomAbs_Cone ) {
839 if( (cuty == GeomAbs_Ellipse) ||
840 ( cuty == GeomAbs_Hyperbola) ||
841 ( cuty == GeomAbs_Parabola) ) {
842 projectable = Standard_False;
843 }
844 }
845 else if ( suty == GeomAbs_Cylinder ) {
846 if (cuty == GeomAbs_Ellipse) {
847 projectable = Standard_False;
848 }
849 }
850 else if ( suty == GeomAbs_Sphere ) {
851 if (cuty == GeomAbs_Circle) {
852 projectable = Standard_False;
853 }
854 }
855 else if ( suty == GeomAbs_Torus ) {
856 if (cuty == GeomAbs_Circle) {
857 projectable = Standard_False;
858 }
859 }
860
0797d9d3 861#ifdef OCCT_DEBUG
7fd59977 862 if (TopOpeBRepTool_GettracePCURV()) {
863 cout<<"--- IsProjectable : ";
864 if (projectable) cout<<"projectable"<<endl;
865 else cout<<"NOT projectable"<<endl;
866 }
867#endif
868
869 return projectable;
870}
871
872
873//=======================================================================
874//function : MakePCurveOnFace
875//purpose :
876//=======================================================================
877
878Handle(Geom2d_Curve) TopOpeBRepTool_CurveTool::MakePCurveOnFace
879(const TopoDS_Shape& S,
880 const Handle(Geom_Curve)& C3D,
881 Standard_Real& TolReached2d,
882 const Standard_Real first, const Standard_Real last)
883
884{
885 Standard_Boolean trim = Standard_False;
886 if (first < last)
887 trim = Standard_True;
888
889 const TopoDS_Face& F = TopoDS::Face(S);
890 Standard_Boolean compminmaxUV = Standard_False;
891 BRepAdaptor_Surface BAS(F,compminmaxUV);
892 GeomAdaptor_Curve GAC;
893 if (trim) GAC.Load(C3D,first,last);
894 else GAC.Load(C3D);
895 Handle(BRepAdaptor_HSurface) BAHS = new BRepAdaptor_HSurface(BAS);
896 Handle(GeomAdaptor_HCurve) BAHC = new GeomAdaptor_HCurve(GAC);
897 ProjLib_ProjectedCurve projcurv(BAHS,BAHC);
898 Handle(Geom2d_Curve) C2D = ::MakePCurve(projcurv);
899 TolReached2d = projcurv.GetTolerance();
900
901 Standard_Real UMin, UMax, VMin, VMax;
902 BRepTools::UVBounds(F,UMin, UMax, VMin, VMax);
903
904 Standard_Real f = GAC.FirstParameter();
905 Standard_Real l = GAC.LastParameter();
906 Standard_Real t =(f+l)*.5;
907 gp_Pnt2d pC2D; C2D->D0(t,pC2D);
908 Standard_Real u2 = pC2D.X();
909 Standard_Real v2 = pC2D.Y();
910
911 if (BAS.GetType() == GeomAbs_Sphere) {
912 // MSV: consider quasiperiodic shift of pcurve
913 Standard_Real VFirst = BAS.FirstVParameter();
914 Standard_Real VLast = BAS.LastVParameter();
915 Standard_Boolean mincond = v2 < VFirst;
916 Standard_Boolean maxcond = v2 > VLast;
917 if (mincond || maxcond) {
918 Handle(Geom2d_Curve) PCT = Handle(Geom2d_Curve)::DownCast(C2D->Copy());
919 // make mirror relative to the isoline of apex -PI/2 or PI/2
920 gp_Trsf2d aTrsf;
c6541a0c
D
921 gp_Pnt2d po(0,-M_PI/2);
922 if (maxcond) po.SetY(M_PI/2);
7fd59977 923 aTrsf.SetMirror(gp_Ax2d(po, gp_Dir2d(1,0)));
924 PCT->Transform(aTrsf);
925 // add translation along U direction on PI
c6541a0c 926 gp_Vec2d vec(M_PI,0);
7fd59977 927 Standard_Real UFirst = BAS.FirstUParameter();
c6541a0c 928 if (u2-UFirst-M_PI > -1e-7) vec.Reverse();
7fd59977 929 PCT->Translate(vec);
930 C2D = PCT;
931 // recompute the test point
932 C2D->D0(t,pC2D);
933 u2 = pC2D.X();
934 v2 = pC2D.Y();
935 }
936 }
937
938 Standard_Real du = 0.;
939 if (BAHS->IsUPeriodic()) {
940 //modified by NIZHNY-MZV Thu Mar 30 10:03:15 2000
941 Standard_Boolean mincond = (UMin - u2 > 1e-7) ? Standard_True : Standard_False;
942 Standard_Boolean maxcond = (u2 - UMax > 1e-7) ? Standard_True : Standard_False;
943 Standard_Boolean decalu = mincond || maxcond;
944 if (decalu) du = ( mincond ) ? BAHS->UPeriod() : -BAHS->UPeriod();
945 //Standard_Boolean decalu = ( u2 < UMin || u2 > UMax);
946 //if (decalu) du = ( u2 < UMin ) ? BAHS->UPeriod() : -BAHS->UPeriod();
947 }
948 Standard_Real dv = 0.;
949 if (BAHS->IsVPeriodic()) {
950 //modified by NIZHNY-MZV Thu Mar 30 10:06:24 2000
951 Standard_Boolean mincond = (VMin - v2 > 1e-7) ? Standard_True : Standard_False;
952 Standard_Boolean maxcond = (v2 - VMax > 1e-7) ? Standard_True : Standard_False;
953 Standard_Boolean decalv = mincond || maxcond;
954 if (decalv) dv = ( mincond ) ? BAHS->VPeriod() : -BAHS->VPeriod();
955 //Standard_Boolean decalv = ( v2 < VMin || v2 > VMax);
956 //if (decalv) dv = ( v2 < VMin ) ? BAHS->VPeriod() : -BAHS->VPeriod();
957 }
958
959 if ( du != 0. || dv != 0.) {
960 Handle(Geom2d_Curve) PCT = Handle(Geom2d_Curve)::DownCast(C2D->Copy());
961 PCT->Translate(gp_Vec2d(du,dv));
962 C2D = PCT;
963 }
964
965 return C2D;
966}