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