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