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
  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
  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   default :
 155     throw Standard_NotImplemented("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   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
 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
 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
 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
 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 }