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