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