0024023: Revamp the OCCT Handle -- downcast (automatic)
[occt.git] / src / TopOpeBRepTool / TopOpeBRepTool_CurveTool.cxx
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 #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
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
83 TopOpeBRepTool_CurveTool::TopOpeBRepTool_CurveTool()
84 {}
85
86 //=======================================================================
87 //function : CurveTool
88 //purpose  : 
89 //=======================================================================
90
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
103 TopOpeBRepTool_CurveTool::TopOpeBRepTool_CurveTool
104 (const TopOpeBRepTool_GeomTool& GT)
105 {
106   SetGeomTool(GT);
107 }
108
109 //=======================================================================
110 //function : ChangeGeomTool
111 //purpose  : 
112 //=======================================================================
113
114 TopOpeBRepTool_GeomTool& TopOpeBRepTool_CurveTool::ChangeGeomTool()
115 {
116   return myGeomTool;
117 }
118
119 //=======================================================================
120 //function : GetGeomTool
121 //purpose  : 
122 //=======================================================================
123
124 const TopOpeBRepTool_GeomTool& TopOpeBRepTool_CurveTool::GetGeomTool()const
125 {
126   return myGeomTool;
127 }
128
129 //=======================================================================
130 //function : SetGeomTool
131 //purpose  : 
132 //=======================================================================
133
134 void TopOpeBRepTool_CurveTool::SetGeomTool
135 (const TopOpeBRepTool_GeomTool& GT)
136 {
137   myGeomTool.Define(GT);
138 }
139
140 //-----------------------------------------------------------------------
141 //function : MakePCurve
142 //purpose  : 
143 //-----------------------------------------------------------------------
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   case GeomAbs_BezierCurve : case GeomAbs_OtherCurve : default :
155     Standard_NotImplemented::Raise("CurveTool::MakePCurve");
156     break;
157   }
158   return C2D;
159 }
160
161 //------------------------------------------------------------------
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 //------------------------------------------------------------------
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 //------------------------------------------------------------------
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 //------------------------------------------------------------------
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 //------------------------------------------------------------------
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
301 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   Standard_Boolean CompC3D = myGeomTool.CompC3D();
311
312   //cout << "MakeCurves begin" << endl;
313   if (!CompC3D) return Standard_False;
314
315   Standard_Boolean CompPC1 = myGeomTool.CompPC1();
316   Standard_Boolean CompPC2 = myGeomTool.CompPC2();
317   Standard_Real tol3d,tol2d;
318   Standard_Boolean RelativeTol;
319   myGeomTool.GetTolerances(tol3d,tol2d,RelativeTol);
320   Standard_Integer NbPntMax = myGeomTool.NbPntMax();
321
322 #ifdef OCCT_DEBUG
323   if (TopOpeBRepTool_GettraceKRO()) KRO_CURVETOOL_APPRO.Start();
324 #endif
325
326 //----------------------------------
327 ///*
328 #ifdef IFV
329   char name[16];
330   char *nm = &name[0];
331   sprintf(name,"C3D_%d",++NbCalls);
332   DrawTrSurf::Set(nm, C3D);
333   sprintf(name,"PC1_%d",NbCalls);
334   DrawTrSurf::Set(nm, PC1);
335   sprintf(name,"PC2_%d",NbCalls);
336   DrawTrSurf::Set(nm, PC2);
337 #endif
338 //*/
339 //---------------------------------------------
340
341   Standard_Integer iparmin = (Standard_Integer)parmin;
342   Standard_Integer iparmax = (Standard_Integer)parmax;
343
344   Handle(Geom_BSplineCurve)&   HC3D = Handle(Geom_BSplineCurve)::DownCast (C3D);
345   Handle(Geom2d_BSplineCurve)& HPC1 = Handle(Geom2d_BSplineCurve)::DownCast (PC1);
346   Handle(Geom2d_BSplineCurve)& HPC2 = Handle(Geom2d_BSplineCurve)::DownCast (PC2);
347
348 //--------------------- IFV - "improving" initial curves
349 #ifdef CurveImprovement
350   Standard_Integer nbpol = HC3D->NbPoles();
351   //cout <<"nbpol = " << nbpol << endl;
352   if(nbpol > 100) {
353     TColgp_Array1OfPnt PolC3D(1, nbpol);
354     TColgp_Array1OfPnt2d PolPC1(1, nbpol);
355     TColgp_Array1OfPnt2d PolPC2(1, nbpol);
356     TColStd_Array1OfBoolean IsValid(1, nbpol);
357     IsValid.Init(Standard_True);
358     Standard_Real tol = Max(1.e-10, 100.*tol3d*tol3d); //tol *= tol; - square distance
359     Standard_Real tl2d = tol*(tol2d*tol2d)/(tol3d*tol3d);
360     Standard_Real def = tol;
361     Standard_Real def2d = tol2d;
362     HC3D->Poles(PolC3D);
363     if(CompPC1) HPC1->Poles(PolPC1);
364     if(CompPC2) HPC2->Poles(PolPC2);
365     
366     Standard_Integer ip = 1, NbPol = 1;
367     Standard_Real d, d1, d2;
368     gp_Pnt P = PolC3D(ip);
369     gp_Pnt2d P1, P2;
370     if(CompPC1) P1 = PolPC1(ip);
371     if(CompPC2) P2 = PolPC2(ip);
372
373     for(ip = 2; ip <= nbpol; ip++) {
374       if( IsValid(ip) ) {
375         d = P.SquareDistance(PolC3D(ip));
376         if(CompPC1) {d1 = P1.SquareDistance(PolPC1(ip));} else {d1 = 0.;}
377         if(CompPC2) {d2 = P2.SquareDistance(PolPC2(ip));} else {d2 = 0.;}
378         if(d > tol || d1 > tl2d || d2 > tl2d ) {
379           Standard_Real dd=0.;
380           if(ip < nbpol ) dd = P.Distance(PolC3D(ip+1));
381           if(ip < nbpol && dd < 10.*tol) {
382             gce_MakeLin mkL(P, PolC3D(ip+1));
383             if(mkL.IsDone()) {
384               gp_Lin L = mkL.Value();
385               d = L.SquareDistance(PolC3D(ip));
386               if(CompPC1) {
387                 gp_Lin2d L1 = gce_MakeLin2d(P1, PolPC1(ip+1));
388                 d1 = L1.SquareDistance(PolPC1(ip));
389               }
390               else d1 = 0.;
391               if(CompPC2) {
392                 gp_Lin2d L2 = gce_MakeLin2d(P2, PolPC2(ip+1));
393                 d2 = L2.SquareDistance(PolPC2(ip));
394               }
395               else d2 = 0.;
396             
397               if(d > def || d1 > def2d || d2 > def2d ) {
398                 NbPol++;
399                 P = PolC3D(ip);
400                 if(CompPC1) P1 = PolPC1(ip);
401                 if(CompPC2) P2 = PolPC2(ip);
402               }
403               else {
404                 IsValid(ip) = Standard_False;
405               }
406             }
407             else {
408               IsValid(ip+1) = Standard_False;
409               NbPol++;
410               P = PolC3D(ip);
411               if(CompPC1) P1 = PolPC1(ip);
412               if(CompPC2) P2 = PolPC2(ip);
413             }
414           }
415           else {
416             NbPol++;
417             P = PolC3D(ip);
418             if(CompPC1) P1 = PolPC1(ip);
419             if(CompPC2) P2 = PolPC2(ip);
420           }
421         }
422         else {
423           IsValid(ip) = Standard_False;
424         }
425       }
426     }
427  
428     if(NbPol < 2) {IsValid(nbpol) = Standard_True; NbPol++;}
429     
430     if(NbPol < nbpol) {
431 #ifdef IFV
432       cout << "NbPol < nbpol : " << NbPol << " " <<  nbpol << endl;
433 #endif
434       TColgp_Array1OfPnt Polc3d(1, NbPol);
435       TColgp_Array1OfPnt2d Polpc1(1, NbPol);
436       TColgp_Array1OfPnt2d Polpc2(1, NbPol);
437       TColStd_Array1OfReal knots(1,NbPol);
438       TColStd_Array1OfInteger mults(1,NbPol);
439       mults.Init(1); mults(1) = 2; mults(NbPol) = 2;
440       Standard_Integer count = 0;
441       for(ip = 1; ip <= nbpol; ip++) {
442         if(IsValid(ip)) {
443           count++;
444           Polc3d(count) = PolC3D(ip);
445           if(CompPC1) Polpc1(count) = PolPC1(ip);
446           if(CompPC2) Polpc2(count) = PolPC2(ip);
447           knots(count) = count;
448         }
449       }
450
451       Polc3d(NbPol) = PolC3D(nbpol);
452       if(CompPC1) Polpc1(NbPol) = PolPC1(nbpol);
453       if(CompPC2) Polpc2(NbPol) = PolPC2(nbpol);
454       
455       HC3D = new Geom_BSplineCurve(Polc3d, knots, mults, 1);
456       if(CompPC1) HPC1 = new Geom2d_BSplineCurve(Polpc1, knots, mults, 1);
457       if(CompPC2) HPC2 = new Geom2d_BSplineCurve(Polpc2, knots, mults, 1);
458       iparmax = NbPol;
459
460 #ifdef IFV
461       sprintf(name,"C3Dmod_%d",NbCalls);
462       nm = &name[0];
463       DrawTrSurf::Set(nm, HC3D);
464       sprintf(name,"PC1mod_%d",NbCalls);
465       nm = &name[0];
466       DrawTrSurf::Set(nm, HPC1);
467       sprintf(name,"PC2mod_%d",NbCalls);
468       nm = &name[0];
469       DrawTrSurf::Set(nm, HPC2);
470 #endif
471
472     }
473   }
474 //--------------- IFV - end "improving"
475 #endif
476
477
478   BRepApprox_Approx Approx;
479
480   Standard_Integer degmin = 4;
481   Standard_Integer degmax = 8;
482   Approx_ParametrizationType parametrization = Approx_ChordLength;
483
484   Standard_Integer npol = HC3D->NbPoles();
485   TColgp_Array1OfPnt Polc3d(1, npol);
486   TColStd_Array1OfReal par(1, npol);
487   HC3D->Poles(Polc3d);
488   gp_Pnt P = Polc3d(1);
489
490   Standard_Boolean IsBad = Standard_False;
491   Standard_Integer ip;
492   Standard_Real ksi = 0., kc, kcprev = 0.;
493   Standard_Real dist;
494   par(1) = 0.;
495   for(ip = 2; ip <= npol; ip++) {
496     dist = P.Distance(Polc3d(ip));
497
498     if(dist < Precision::Confusion()) {
499       IsBad = Standard_True;
500       break;
501     }
502
503     par(ip) = par(ip-1) + dist; 
504     P = Polc3d(ip);
505   }
506
507   if(!IsBad) {
508
509     TColStd_Array1OfReal Curvature(1, npol);
510
511     if(npol > 3) {
512       Standard_Integer ic = 1;
513       for(ip = 2; ip <= npol-1; ip += npol-3) {
514         gp_Vec v1(Polc3d(ip-1),Polc3d(ip));
515         gp_Vec v2(Polc3d(ip),Polc3d(ip+1));
516         if(!v1.IsParallel(v2, 1.e-4)) {
517           gce_MakeCirc mc(Polc3d(ip-1),Polc3d(ip),Polc3d(ip+1));
518           gp_Circ cir = mc.Value();
519           kc = 1./cir.Radius();
520         }
521         else kc = 0.;
522     
523         Curvature(ic) = kc;
524         ic = npol;
525       }
526     }
527     else if(npol == 3) {
528       ip = 2;
529       gp_Vec v1(Polc3d(ip-1),Polc3d(ip));
530       gp_Vec v2(Polc3d(ip),Polc3d(ip+1));
531       if(!v1.IsParallel(v2, 1.e-4)) {
532         gce_MakeCirc mc(Polc3d(ip-1),Polc3d(ip),Polc3d(ip+1));
533         gp_Circ cir = mc.Value();
534         kc = 1./cir.Radius();
535       }
536       else kc = 0.;
537       Curvature(1) = Curvature(npol) = kc;
538     }
539     else {
540       Curvature(1) = Curvature(npol) = 0.;
541     }
542     
543     ip = 1;
544     Standard_Real dt = par(ip+1) - par(ip);
545     Standard_Real dx = (Polc3d(ip+1).X() - Polc3d(ip).X())/dt,
546                   dy = (Polc3d(ip+1).Y() - Polc3d(ip).Y())/dt,
547                   dz = (Polc3d(ip+1).Z() - Polc3d(ip).Z())/dt;
548     Standard_Real dx1, dy1, dz1, d2x, d2y, d2z, d2t;
549
550     for(ip = 2; ip <= npol-1; ip++) {
551       dt = par(ip+1) - par(ip);
552       dx1 = (Polc3d(ip+1).X() - Polc3d(ip).X())/dt;
553       dy1 = (Polc3d(ip+1).Y() - Polc3d(ip).Y())/dt,
554       dz1 = (Polc3d(ip+1).Z() - Polc3d(ip).Z())/dt;
555       d2t = 2./(par(ip+1) - par(ip-1));
556       d2x = d2t*(dx1 - dx);
557       d2y = d2t*(dy1 - dy);
558       d2z = d2t*(dz1 - dz);
559       Curvature(ip) = Sqrt(d2x*d2x + d2y*d2y + d2z*d2z);
560       dx = dx1; dy = dy1; dz = dz1;
561     }
562
563     Standard_Real crit = 1000.; // empirical criterion !!!
564
565     dt = par(2) - par(1);
566     kcprev = (Curvature(2) - Curvature(1))/dt; 
567     for(ip = 2; ip <= npol-1; ip++) {
568       dt = par(ip+1) - par(ip);
569       kc = (Curvature(ip+1) - Curvature(ip))/dt;
570       ksi = ksi + Abs(kc - kcprev);
571       if(ksi > crit) {IsBad = Standard_True;break;}
572       kc = kcprev; 
573     }
574
575   }
576   //cout << NbCalls << " ksi = " << ksi << endl;
577   //cout << "IsBad = " << IsBad << endl;
578
579   if(IsBad){
580     Standard_Real tt = Min(10.*tol3d,Precision::Approximation());
581     tol2d = tt * tol2d / tol3d;
582     tol3d = tt;
583     NbPntMax = 40;
584     degmax = 4;
585     parametrization = Approx_Centripetal;
586   }
587  
588   Standard_Integer nitmax = 0; // use projection only
589   Standard_Boolean withtangency = Standard_True;
590   
591   Standard_Boolean compminmaxUV = Standard_True;
592   BRepAdaptor_Surface BAS1(TopoDS::Face(S1),compminmaxUV);
593   BRepAdaptor_Surface BAS2(TopoDS::Face(S2),compminmaxUV);
594
595
596   Handle(BRepApprox_ApproxLine) AL;
597   AL = new BRepApprox_ApproxLine(HC3D,HPC1,HPC2);
598
599     Approx.SetParameters(tol3d,tol2d,RelativeTol,degmin,degmax,nitmax,NbPntMax,withtangency,
600                          parametrization);
601
602     if     (CompC3D && CompPC1 && BAS1.GetType() == GeomAbs_Plane) { 
603       //-- The curve X,Y,Z and U2,V2 is approximated
604       Approx.Perform(BAS1,BAS2,AL,CompC3D,Standard_False,CompPC2,iparmin,iparmax);
605     }
606     
607     else if(CompC3D && CompPC2 && BAS2.GetType() == GeomAbs_Plane) {
608       //-- The curve X,Y,Z and U1,V1 is approximated
609       Approx.Perform(BAS1,BAS2,AL,CompC3D,CompPC1,Standard_False,iparmin,iparmax);
610     }
611
612     else { 
613       Approx.Perform(BAS1,BAS2,AL,CompC3D,CompPC1,CompPC2,iparmin,iparmax);
614     }
615
616   // MSV Nov 9, 2001: do not raise exception in the case of failure,
617   //                  but return status
618
619   Standard_Boolean done = Approx.IsDone();
620   done = done && CheckApproxResults(Approx);
621
622   if (done) {
623     if     (CompC3D && CompPC1 && BAS1.GetType() == GeomAbs_Plane) { 
624       C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
625       PC1new = ::MakeCurve2DfromWLineApproxAndPlane(Approx,BAS1.Plane());
626       if (CompPC2) PC2new = ::MakeCurve2DfromWLineApprox(Approx,2);
627     }
628     else if(CompC3D && CompPC2 && BAS2.GetType() == GeomAbs_Plane) {
629       C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
630       if (CompPC1) PC1new = ::MakeCurve2DfromWLineApprox(Approx,2);
631       PC2new = ::MakeCurve2DfromWLineApproxAndPlane(Approx,BAS2.Plane());
632     }
633     else { 
634       if (CompC3D) C3Dnew = ::MakeCurve3DfromWLineApprox(Approx,1);
635       if (CompPC1) PC1new = ::MakeCurve2DfromWLineApprox(Approx,2);
636       if (CompPC2) {
637         Standard_Integer i32 = (CompPC1) ? 3 : 2;
638         PC2new = ::MakeCurve2DfromWLineApprox(Approx,i32);
639       }
640     }
641
642     // check the pcurves relatively the faces bounds
643     if (CompPC1)
644       done = done && CheckPCurve(PC1new,TopoDS::Face(S1));
645     if (CompPC2)
646       done = done && CheckPCurve(PC2new,TopoDS::Face(S2));
647   }
648
649   if (!done) {
650     if (CompC3D) C3Dnew.Nullify();
651     if (CompPC1) PC1new.Nullify();
652     if (CompPC2) PC2new.Nullify();
653     return Standard_False;
654   }
655
656 #ifdef IFV
657     sprintf(name,"C3Dnew_%d", NbCalls);
658     nm = &name[0];
659     DrawTrSurf::Set(nm, C3Dnew);
660     if (CompPC1) {
661       sprintf(name,"PC1new_%d", NbCalls);
662       nm = &name[0];
663       DrawTrSurf::Set(nm, PC1new);
664     }
665     if (CompPC2) {
666       sprintf(name,"PC2new_%d", NbCalls);
667       nm = &name[0];
668       DrawTrSurf::Set(nm, PC2new);
669     }
670
671 #endif
672
673   TolReached3d = Approx.TolReached3d();
674   TolReached2d = Approx.TolReached2d();
675 #ifdef IFV
676   cout << NbCalls << " : Tol = " << TolReached3d << " " << TolReached2d << endl;
677 #endif
678
679   Standard_Boolean bf, bl;
680
681   Handle(Geom_BSplineCurve) Curve (Handle(Geom_BSplineCurve)::DownCast(C3Dnew));
682   if(!Curve.IsNull()) {
683     GeomLib_CheckBSplineCurve cbsc(Curve, 1.e-7, 0.1);
684     cbsc.NeedTangentFix(bf, bl);
685
686 #ifdef OCCT_DEBUG
687     if (TopOpeBRepTool_GettraceCHKBSPL()) {
688       if(bf || bl) {
689         cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
690         cout<<" Last = "<<bl<<endl;
691       }
692     }
693 #endif
694     cbsc.FixTangent(bf, bl);
695   }
696   
697   Handle(Geom2d_BSplineCurve) Curve2df (Handle(Geom2d_BSplineCurve)::DownCast(PC1new));
698   if(!Curve2df.IsNull()) {
699     GeomLib_Check2dBSplineCurve cbsc2df(Curve2df, 1.e-7, 0.1);
700     cbsc2df.NeedTangentFix(bf, bl);
701 #ifdef OCCT_DEBUG
702     if (TopOpeBRepTool_GettraceCHKBSPL()) {
703       if(bf || bl) {
704         cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
705         cout<<" Last = "<<bl<<endl;
706       }
707     }
708 #endif
709     cbsc2df.FixTangent(bf, bl);
710   }
711   
712   Handle(Geom2d_BSplineCurve) Curve2ds (Handle(Geom2d_BSplineCurve)::DownCast(PC2new));
713   if(!Curve2ds.IsNull()) {
714     GeomLib_Check2dBSplineCurve cbsc2ds(Curve2ds, 1.e-7, 0.1);
715     cbsc2ds.NeedTangentFix(bf, bl);
716 #ifdef OCCT_DEBUG
717     if (TopOpeBRepTool_GettraceCHKBSPL()) {
718       if(bf || bl) {
719         cout<<"Problem orientation GeomLib_CheckBSplineCurve : First = "<<bf;
720         cout<<" Last = "<<bl<<endl;
721       }
722     }
723 #endif
724     cbsc2ds.FixTangent(bf, bl);
725   }
726   
727 #ifdef OCCT_DEBUG
728   if (TopOpeBRepTool_GettraceKRO()) KRO_CURVETOOL_APPRO.Stop();
729 #endif
730 //  cout << "MakeCurves end" << endl;
731
732   return Standard_True;
733 }
734
735
736 //=======================================================================
737 //function : MakeBSpline1fromPnt
738 //purpose  : 
739 //=======================================================================
740
741 Handle(Geom_Curve) TopOpeBRepTool_CurveTool::MakeBSpline1fromPnt
742 (const TColgp_Array1OfPnt& Points)
743 {
744   Standard_Integer Degree = 1;
745   
746   Standard_Integer i,nbpoints = Points.Length();
747   Standard_Integer nbknots = nbpoints - Degree +1;
748   
749   //  First compute the parameters
750   //  Standard_Real length = 0.;
751   //  TColStd_Array1OfReal parameters(1,nbpoints);
752   //  for (i = 1; i < nbpoints; i++) {
753   //    parameters(i) = length;
754   //    Standard_Real dist = Points(i).Distance(Points(i+1));
755   //    length += dist;
756   //  }
757   //  parameters(nbpoints) = length;
758   
759   // knots and multiplicities
760   TColStd_Array1OfReal knots(1,nbknots);
761   TColStd_Array1OfInteger mults(1,nbknots);
762   mults.Init(1);
763   mults(1) = mults(nbknots) = Degree + 1;
764   
765   //  knots(1) = 0;
766   //  for (i=2;i<nbknots;i++) knots(i) = (parameters(i) + parameters(i+1)) /2.;
767   //  knots(nbknots) = length;
768   
769   // take point index as parameter : JYL 01/AUG/94
770   for (i = 1; i <= nbknots; i++) knots(i) = (Standard_Real)i;
771   Handle(Geom_Curve) C = new Geom_BSplineCurve(Points,knots,mults,Degree);
772   return C;
773 }
774
775
776 //=======================================================================
777 //function : MakeBSpline1fromPnt2d
778 //purpose  : 
779 //=======================================================================
780
781 Handle(Geom2d_Curve) TopOpeBRepTool_CurveTool::MakeBSpline1fromPnt2d
782 (const TColgp_Array1OfPnt2d& Points)
783 {
784   Standard_Integer Degree = 1;
785   
786   Standard_Integer i,nbpoints = Points.Length();
787   Standard_Integer nbknots = nbpoints - Degree +1;
788   
789   //  First compute the parameters
790   //  Standard_Real length = 0;
791   //  TColStd_Array1OfReal parameters(1,nbpoints);
792   //  for (i = 1; i < nbpoints; i++) {
793   //    parameters(i) = length;
794   //    Standard_Real dist = Points(i).Distance(Points(i+1));
795   //    length += dist;
796   //  }
797   //  parameters(nbpoints) = length;
798   
799   // knots and multiplicities
800   TColStd_Array1OfReal knots(1,nbknots);
801   TColStd_Array1OfInteger mults(1,nbknots);
802   mults.Init(1);
803   mults(1) = mults(nbknots) = Degree + 1;
804   
805   //  knots(1) = 0;
806   //  for (i=2;i<nbknots;i++) knots(i) = (parameters(i) + parameters(i+1)) /2.;
807   //  knots(nbknots) = length;
808   
809   // take point index as parameter : JYL 01/AUG/94
810   for (i = 1; i <= nbknots; i++) knots(i) = (Standard_Real)i;
811   Handle(Geom2d_Curve) C = new Geom2d_BSplineCurve(Points,knots,mults,Degree);
812   return C;
813 }
814
815 //=======================================================================
816 //function : IsProjectable
817 //purpose  : 
818 //=======================================================================
819
820 Standard_Boolean TopOpeBRepTool_CurveTool::IsProjectable
821 (const TopoDS_Shape& S, const Handle(Geom_Curve)& C3D)
822 {
823   const TopoDS_Face& F = TopoDS::Face(S);
824   Standard_Boolean compminmaxUV = Standard_False;
825   BRepAdaptor_Surface BAS(F,compminmaxUV);
826   GeomAbs_SurfaceType suty = BAS.GetType();
827   GeomAdaptor_Curve GAC(C3D);
828   GeomAbs_CurveType cuty = GAC.GetType();
829
830   // --------
831   // avoid projection of 3d curve on surface in case
832   // of a quadric (ellipse,hyperbola,parabola) on a cone.
833   // Projection fails when the curve in not fully inside the UV domain 
834   // of the cone : only part of 2d curve is built.
835   // NYI : projection of quadric on cone (crossing cone domain)
836   // --------
837   
838   Standard_Boolean projectable = Standard_True;
839   if ( suty == GeomAbs_Cone ) {
840     if( (cuty == GeomAbs_Ellipse) || 
841        ( cuty == GeomAbs_Hyperbola) || 
842        ( cuty == GeomAbs_Parabola) ) {
843       projectable = Standard_False;
844     }
845   }
846   else if ( suty == GeomAbs_Cylinder ) {
847     if (cuty == GeomAbs_Ellipse) {
848       projectable = Standard_False;
849     }
850   }
851   else if ( suty == GeomAbs_Sphere ) {
852     if (cuty == GeomAbs_Circle) {
853       projectable = Standard_False;
854     }
855   }
856   else if ( suty == GeomAbs_Torus ) {
857     if (cuty == GeomAbs_Circle) {
858       projectable = Standard_False;
859     }
860   }
861   
862 #ifdef OCCT_DEBUG
863   if (TopOpeBRepTool_GettracePCURV()) {
864     cout<<"--- IsProjectable : "; 
865     if (projectable) cout<<"projectable"<<endl;
866     else cout<<"NOT projectable"<<endl;
867   }
868 #endif
869   
870   return projectable;
871 }
872
873
874 //=======================================================================
875 //function : MakePCurveOnFace
876 //purpose  : 
877 //=======================================================================
878
879 Handle(Geom2d_Curve) TopOpeBRepTool_CurveTool::MakePCurveOnFace
880 (const TopoDS_Shape& S,
881  const Handle(Geom_Curve)& C3D,
882  Standard_Real& TolReached2d,
883  const Standard_Real first, const Standard_Real last)
884  
885 {
886   Standard_Boolean trim = Standard_False;
887   if (first < last)
888     trim = Standard_True;
889
890   const TopoDS_Face& F = TopoDS::Face(S);
891   Standard_Boolean compminmaxUV = Standard_False;
892   BRepAdaptor_Surface BAS(F,compminmaxUV);
893   GeomAdaptor_Curve GAC;
894   if (trim) GAC.Load(C3D,first,last);
895   else      GAC.Load(C3D);
896   Handle(BRepAdaptor_HSurface) BAHS = new BRepAdaptor_HSurface(BAS);
897   Handle(GeomAdaptor_HCurve) BAHC = new GeomAdaptor_HCurve(GAC);
898   ProjLib_ProjectedCurve projcurv(BAHS,BAHC);
899   Handle(Geom2d_Curve) C2D = ::MakePCurve(projcurv);
900   TolReached2d = projcurv.GetTolerance();
901
902   Standard_Real UMin, UMax, VMin, VMax;
903   BRepTools::UVBounds(F,UMin, UMax, VMin, VMax);
904
905   Standard_Real f = GAC.FirstParameter();
906   Standard_Real l = GAC.LastParameter();
907   Standard_Real t =(f+l)*.5;
908   gp_Pnt2d pC2D; C2D->D0(t,pC2D);
909   Standard_Real u2 = pC2D.X();
910   Standard_Real v2 = pC2D.Y();
911
912   if (BAS.GetType() == GeomAbs_Sphere) {
913     // MSV: consider quasiperiodic shift of pcurve
914     Standard_Real VFirst = BAS.FirstVParameter();
915     Standard_Real VLast = BAS.LastVParameter();
916     Standard_Boolean mincond = v2 < VFirst;
917     Standard_Boolean maxcond = v2 > VLast;
918     if (mincond || maxcond) {
919       Handle(Geom2d_Curve) PCT = Handle(Geom2d_Curve)::DownCast(C2D->Copy());
920       // make mirror relative to the isoline of apex -PI/2 or PI/2
921       gp_Trsf2d aTrsf;
922       gp_Pnt2d po(0,-M_PI/2);
923       if (maxcond) po.SetY(M_PI/2);
924       aTrsf.SetMirror(gp_Ax2d(po, gp_Dir2d(1,0)));
925       PCT->Transform(aTrsf);
926       // add translation along U direction on PI
927       gp_Vec2d vec(M_PI,0);
928       Standard_Real UFirst = BAS.FirstUParameter();
929       if (u2-UFirst-M_PI > -1e-7) vec.Reverse();
930       PCT->Translate(vec);
931       C2D = PCT;
932       // recompute the test point
933       C2D->D0(t,pC2D);
934       u2 = pC2D.X();
935       v2 = pC2D.Y();
936     }
937   }
938
939   Standard_Real du = 0.;
940   if (BAHS->IsUPeriodic()) {
941     //modified by NIZHNY-MZV  Thu Mar 30 10:03:15 2000
942     Standard_Boolean mincond = (UMin - u2 > 1e-7) ? Standard_True : Standard_False;
943     Standard_Boolean maxcond = (u2 - UMax > 1e-7) ? Standard_True : Standard_False;
944     Standard_Boolean decalu = mincond || maxcond;
945     if (decalu) du = ( mincond ) ? BAHS->UPeriod() : -BAHS->UPeriod();
946     //Standard_Boolean decalu = ( u2 < UMin || u2 > UMax);
947     //if (decalu) du = ( u2 < UMin ) ? BAHS->UPeriod() : -BAHS->UPeriod();
948   }
949   Standard_Real dv = 0.;
950   if (BAHS->IsVPeriodic()) {
951     //modified by NIZHNY-MZV  Thu Mar 30 10:06:24 2000
952     Standard_Boolean mincond = (VMin - v2 > 1e-7) ? Standard_True : Standard_False;
953     Standard_Boolean maxcond = (v2 - VMax > 1e-7) ? Standard_True : Standard_False;
954     Standard_Boolean decalv = mincond || maxcond;
955     if (decalv) dv = ( mincond ) ? BAHS->VPeriod() : -BAHS->VPeriod();
956     //Standard_Boolean decalv = ( v2 < VMin || v2 > VMax);
957     //if (decalv) dv = ( v2 < VMin ) ? BAHS->VPeriod() : -BAHS->VPeriod();
958   }
959   
960   if ( du != 0. || dv != 0.) {
961     Handle(Geom2d_Curve) PCT = Handle(Geom2d_Curve)::DownCast(C2D->Copy());
962     PCT->Translate(gp_Vec2d(du,dv));
963     C2D = PCT;
964   }
965
966   return C2D;
967 }