src/Extrema/Extrema_GenExtCS.cxx

   1 // Created on: 1995-07-18
   2 // Created by: Modelistation
   3 // Copyright (c) 1995-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <Extrema_GenExtCS.ixx>
  18
  19 #include <math_Vector.hxx>
  20 #include <math_FunctionSetRoot.hxx>
  21 #include <Precision.hxx>
  22
  23 #include <Geom_Line.hxx>
  24 #include <Adaptor3d_HCurve.hxx>
  25 #include <GeomAdaptor_Curve.hxx>
  26
  27 #include <Extrema_ExtCC.hxx>
  28 #include <Extrema_POnCurv.hxx>
  29 #include <Extrema_ExtPS.hxx>
  30
  31 //=======================================================================
  32 //function : Extrema_GenExtCS
  33 //purpose  :
  34 //=======================================================================
  35
  36 Extrema_GenExtCS::Extrema_GenExtCS()
  37 {
  38   myDone = Standard_False;
  39   myInit = Standard_False;
  40 }
  41
  42 //=======================================================================
  43 //function : Extrema_GenExtCS
  44 //purpose  :
  45 //=======================================================================
  46
  47 Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C,
  48   const Adaptor3d_Surface& S,
  49   const Standard_Integer NbT,
  50   const Standard_Integer NbU,
  51   const Standard_Integer NbV,
  52   const Standard_Real Tol1,
  53   const Standard_Real Tol2)
  54 {
  55   Initialize(S, NbU, NbV, Tol2);
  56   Perform(C, NbT, Tol1);
  57 }
  58
  59 //=======================================================================
  60 //function : Extrema_GenExtCS
  61 //purpose  :
  62 //=======================================================================
  63
  64 Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C,
  65   const Adaptor3d_Surface& S,
  66   const Standard_Integer NbT,
  67   const Standard_Integer NbU,
  68   const Standard_Integer NbV,
  69   const Standard_Real tmin,
  70   const Standard_Real tsup,
  71   const Standard_Real Umin,
  72   const Standard_Real Usup,
  73   const Standard_Real Vmin,
  74   const Standard_Real Vsup,
  75   const Standard_Real Tol1,
  76   const Standard_Real Tol2)
  77 {
  78   Initialize(S, NbU, NbV, Umin,Usup,Vmin,Vsup,Tol2);
  79   Perform(C, NbT, tmin, tsup, Tol1);
  80 }
  81
  82 //=======================================================================
  83 //function : Initialize
  84 //purpose  :
  85 //=======================================================================
  86
  87 void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
  88   const Standard_Integer NbU,
  89   const Standard_Integer NbV,
  90   const Standard_Real Tol2)
  91 {
  92   myumin = S.FirstUParameter();
  93   myusup = S.LastUParameter();
  94   myvmin = S.FirstVParameter();
  95   myvsup = S.LastVParameter();
  96   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,Tol2);
  97 }
  98
  99 //=======================================================================
 100 //function : Initialize
 101 //purpose  :
 102 //=======================================================================
 103
 104 void Extrema_GenExtCS::Initialize(const Adaptor3d_Surface& S,
 105   const Standard_Integer NbU,
 106   const Standard_Integer NbV,
 107   const Standard_Real Umin,
 108   const Standard_Real Usup,
 109   const Standard_Real Vmin,
 110   const Standard_Real Vsup,
 111   const Standard_Real Tol2)
 112 {
 113   myS = (Adaptor3d_SurfacePtr)&S;
 114   myusample = NbU;
 115   myvsample = NbV;
 116   myumin = Umin;
 117   myusup = Usup;
 118   myvmin = Vmin;
 119   myvsup = Vsup;
 120   mytol2 = Tol2;
 121 }
 122
 123 //=======================================================================
 124 //function : Perform
 125 //purpose  :
 126 //=======================================================================
 127
 128 void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
 129   const Standard_Integer NbT,
 130   const Standard_Real Tol1)
 131 {
 132   mytmin = C.FirstParameter();
 133   mytsup = C.LastParameter();
 134   Perform(C, NbT, mytmin, mytsup,Tol1);
 135 }
 136
 137 //=======================================================================
 138 //function : Perform
 139 //purpose  :
 140 //=======================================================================
 141
 142 void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
 143   const Standard_Integer NbT,
 144   const Standard_Real tmin,
 145   const Standard_Real tsup,
 146   const Standard_Real Tol1)
 147 {
 148   myDone = Standard_False;
 149   myF.Initialize(C,*myS);
 150   mytmin = tmin;
 151   mytsup = tsup;
 152   mytol1 = Tol1;
 153   mytsample = NbT;
 154   // Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
 155
 156   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
 157   if (Precision::IsInfinite(trimusup)){
 158     trimusup = 1.0e+10;
 159   }
 160   if (Precision::IsInfinite(trimvsup)){
 161     trimvsup = 1.0e+10;
 162   }
 163   if (Precision::IsInfinite(trimumin)){
 164     trimumin = -1.0e+10;
 165   }
 166   if (Precision::IsInfinite(trimvmin)){
 167     trimvmin = -1.0e+10;
 168   }
 169   //
 170   math_Vector Tol(1, 3);
 171   Tol(1) = mytol1;
 172   Tol(2) = mytol2;
 173   Tol(3) = mytol2;
 174   math_Vector UV(1,3), UVinf(1,3), UVsup(1,3);
 175   UVinf(1) = mytmin;
 176   UVinf(2) = trimumin;
 177   UVinf(3) = trimvmin;
 178
 179   UVsup(1) = mytsup;
 180   UVsup(2) = trimusup;
 181   UVsup(3) = trimvsup;
 182   // 18/02/02 akm vvv : (OCC163) bad extremas - extrusion surface versus the line.
 183   //                    Try to preset the initial solution as extrema between
 184   //                    extrusion direction and the curve.
 185   if (myS->GetType() == GeomAbs_SurfaceOfExtrusion)
 186   {
 187     gp_Dir aDir = myS->Direction();
 188     Handle(Adaptor3d_HCurve) aCurve = myS->BasisCurve();
 189     Standard_Real dfUFirst = aCurve->FirstParameter();
 190     // Create iso line of U=U0
 191     GeomAdaptor_Curve anAx(new Geom_Line(aCurve->Value(dfUFirst), aDir),
 192       trimvmin, trimvsup);
 193     Extrema_ExtCC aLocator(C, anAx);
 194     if (aLocator.IsDone() && aLocator.NbExt()>0)
 195     {
 196       Standard_Integer iExt;
 197       // Try to find all extremas
 198       Extrema_POnCurv aP1, aP2;
 199       for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
 200       {
 201         aLocator.Points (iExt, aP1, aP2);
 202         // Parameter on curve
 203         UV(1) = aP1.Parameter();
 204         // To find parameters on surf, try ExtPS
 205         Extrema_ExtPS aPreciser (aP1.Value(), *myS, mytol2, mytol2);
 206         if (aPreciser.IsDone())
 207         {
 208           // Managed to find extremas between point and surface
 209           Standard_Integer iPExt;
 210           for (iPExt=1; iPExt<=aPreciser.NbExt(); iPExt++)
 211           {
 212             aPreciser.Point(iPExt).Parameter(UV(2),UV(3));
 213             math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
 214           }
 215         }
 216         else
 217         {
 218           // Failed... try the point on iso line
 219           UV(2) = dfUFirst;
 220           UV(3) = aP2.Parameter();
 221           math_FunctionSetRoot S1 (myF,UV,Tol,UVinf,UVsup);
 222         }
 223       } // for (iExt=1; iExt<=aLocator.NbExt(); iExt++)
 224     } // if (aLocator.IsDone() && aLocator.NbExt()>0)
 225   } // if (myS.Type() == GeomAbs_ExtrusionSurface)
 226   else
 227   {
 228     Standard_Real aCUAdd = (mytsup - mytmin) / mytsample;
 229     Standard_Real aSUAdd = (myusup - myumin) / myusample;
 230     Standard_Real aSVAdd = (myvsup - myvmin) / myvsample;
 231     Standard_Real tres = C.Resolution(1.);
 232     Standard_Real ures = myS->UResolution(1.);
 233     Standard_Real vres = myS->VResolution(1.);
 234     tres = aCUAdd / tres;
 235     ures = aSUAdd / ures;
 236     vres = aSVAdd / vres;
 237     Standard_Real minres = Min(tres, Min(ures, vres));
 238     Standard_Real factor = 5.;
 239     Standard_Integer maxnbs = 50;
 240     minres *= factor;
 241     if(minres > Epsilon(1.))
 242     {
 243       if(tres > minres)
 244       {
 245         Standard_Real rsample = mytsample * tres / minres;
 246         if(rsample > maxnbs)
 247         {
 248           mytsample = maxnbs;
 249         }
 250         else
 251         {
 252           mytsample = RealToInt(rsample);
 253         }
 254         aCUAdd = (mytsup - mytmin) / mytsample;
 255       }
 256       if(ures > minres)
 257       {
 258         Standard_Real rsample = myusample * ures / minres;
 259         if(rsample > maxnbs)
 260         {
 261           myusample = maxnbs;
 262         }
 263         else
 264         {
 265           myusample = RealToInt(rsample);
 266         }
 267         aSUAdd = (myusup - myumin) / myusample;
 268       }
 269       if(vres > minres)
 270       {
 271         Standard_Real rsample = myvsample * vres / minres;
 272         if(rsample > maxnbs)
 273         {
 274           myvsample = maxnbs;
 275         }
 276         else
 277         {
 278           myvsample = RealToInt(rsample);
 279         }
 280         aSVAdd = (myvsup - myvmin) / myvsample;
 281       }
 282     }
 283
 284     TColgp_HArray1OfPnt aCPs(1, mytsample);
 285     TColgp_HArray2OfPnt aSPs(1, myusample, 1, myvsample);
 286     Standard_Integer aRestIterCount = 3;
 287     // The value is calculated by the bug CR23830.
 288     Standard_Integer aCUDen = 2, aSUDen = 2, aSVDen = 2;
 289     Standard_Boolean anAreAvSqsInited = Standard_False;
 290     Standard_Real aCUSq = 0, aSUSq = 0, aSVSq = 0;
 291     while (aRestIterCount--)
 292     {
 293       Standard_Real aMinCU = 0., aMinSU = 0., aMinSV = 0., aMaxCU = 0., aMaxSU = 0., aMaxSV = 0.;
 294       Standard_Real aMinSqDist = DBL_MAX, aMaxSqDist = -DBL_MAX;
 295       for (Standard_Integer aSUNom = 1; aSUNom < aSUDen; aSUNom += 2)
 296       {
 297         Standard_Real aSU0 = myumin + (aSUNom * aSUAdd) / aSUDen;
 298         for (Standard_Integer aSVNom = 1; aSVNom < aSVDen; aSVNom += 2)
 299         {
 300           Standard_Real aSV0 = myvmin + (aSVNom * aSVAdd) / aSVDen;
 301           for (Standard_Integer aCUNom = 1; aCUNom < aCUDen; aCUNom += 2)
 302           {
 303             Standard_Real aCU0 = mytmin + (aCUNom * aCUAdd) / aCUDen;
 304             Standard_Real aCU = aCU0;
 305             for (Standard_Integer aCUI = 1; aCUI <= mytsample;
 306               aCUI++, aCU += aCUAdd)
 307             {
 308               aCPs.SetValue(aCUI, C.Value(aCU));
 309             }
 310             //
 311             aCU = aCU0;
 312             Standard_Real aSU = aSU0;
 313             for (Standard_Integer aSUI = 1; aSUI <= myusample;
 314               aSUI++, aSU += aSUAdd)
 315             {
 316               Standard_Real aSV = aSV0;
 317               for (Standard_Integer aSVI = 1; aSVI <= myvsample;
 318                 aSVI++, aSV += aSVAdd)
 319               {
 320                 gp_Pnt aSP = myS->Value(aSU, aSV);
 321                 aSPs.ChangeValue(aSUI, aSVI) = aSP;
 322                 Standard_Real aCU = aCU0;
 323                 for (Standard_Integer aCUI = 1; aCUI <= mytsample;
 324                   aCUI++, aCU += aCUAdd)
 325                 {
 326                   Standard_Real aSqDist = aSP.SquareDistance(aCPs.Value(aCUI));
 327                   if (aSqDist < aMinSqDist)
 328                   {
 329                     aMinSqDist = aSqDist;
 330                     aMinCU = aCU;
 331                     aMinSU = aSU;
 332                     aMinSV = aSV;
 333                   }
 334                   if (aMaxSqDist < aSqDist)
 335                   {
 336                     aMaxSqDist = aSqDist;
 337                     aMaxCU = aCU;
 338                     aMaxSU = aSU;
 339                     aMaxSV = aSV;
 340                   }
 341                 }
 342               }
 343             }
 344           }
 345         }
 346       }
 347       // Find min approximation.
 348       UV(1) = aMinCU;
 349       UV(2) = aMinSU;
 350       UV(3) = aMinSV;
 351       math_FunctionSetRoot anA(myF, UV, Tol, UVinf, UVsup, 100, aRestIterCount);
 352       // Find max approximation.
 353       if (!anA.IsDivergent() || !aRestIterCount)
 354       {
 355         UV(1) = aMaxCU;
 356         UV(2) = aMaxSU;
 357         UV(3) = aMaxSV;
 358         math_FunctionSetRoot aFunc(myF, UV, Tol, UVinf, UVsup);
 359         break;
 360       }
 361       //
 362       if (!anAreAvSqsInited)
 363       {
 364         anAreAvSqsInited = Standard_True;
 365         //
 366         const gp_Pnt * aCP1 = &aCPs.Value(1);
 367         for (Standard_Integer aCUI = 2; aCUI <= mytsample; aCUI++)
 368         {
 369           const gp_Pnt & aCP2 = aCPs.Value(aCUI);
 370           aCUSq += aCP1->SquareDistance(aCP2);
 371           aCP1 = &aCP2;
 372         }
 373         aCUSq /= mytsample - 1;
 374         //
 375         for (Standard_Integer aSVI = 1; aSVI <= myvsample; aSVI++)
 376         {
 377           const gp_Pnt * aSP1 = &aSPs.Value(1, aSVI);
 378           for (Standard_Integer aSUI = 2; aSUI <= myusample; aSUI++)
 379           {
 380             const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
 381             aSUSq += aSP1->SquareDistance(aSP2);
 382             aSP1 = &aSP2;
 383           }
 384         }
 385         aSUSq /= myvsample * (myusample - 1);
 386         //
 387         for (Standard_Integer aSUI = 1; aSUI <= myusample; aSUI++)
 388         {
 389           const gp_Pnt * aSP1 = &aSPs.Value(aSUI, 1);
 390           for (Standard_Integer aSVI = 2; aSVI <= myvsample; aSVI++)
 391           {
 392             const gp_Pnt & aSP2 = aSPs.Value(aSUI, aSVI);
 393             aSVSq += aSP1->SquareDistance(aSP2);
 394             aSP1 = &aSP2;
 395           }
 396         }
 397         aSVSq /= myusample * (myvsample - 1);
 398       }
 399       //
 400       if ((aSUSq <= aCUSq) && (aSVSq <= aCUSq))
 401       {
 402         aCUDen += aCUDen;
 403         aCUSq /= 4;
 404       }
 405       else if ((aCUSq <= aSUSq) && (aSVSq <= aSUSq))
 406       {
 407         aSUDen += aSUDen;
 408         aSUSq /= 4;
 409       }
 410       else
 411       {
 412         aSVDen += aSVDen;
 413         aSVSq /= 4;
 414       }
 415     }
 416   }
 417
 418   myDone = Standard_True;
 419 }
 420
 421 //=======================================================================
 422 //function : IsDone
 423 //purpose  :
 424 //=======================================================================
 425
 426 Standard_Boolean Extrema_GenExtCS::IsDone() const
 427 {
 428   return myDone;
 429 }
 430
 431 //=======================================================================
 432 //function : NbExt
 433 //purpose  :
 434 //=======================================================================
 435
 436 Standard_Integer Extrema_GenExtCS::NbExt() const
 437 {
 438   if (!IsDone()) { StdFail_NotDone::Raise(); }
 439   return myF.NbExt();
 440 }
 441
 442 //=======================================================================
 443 //function : Value
 444 //purpose  :
 445 //=======================================================================
 446
 447 Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const
 448 {
 449   if (!IsDone()) { StdFail_NotDone::Raise(); }
 450   return myF.SquareDistance(N);
 451 }
 452
 453 //=======================================================================
 454 //function : PointOnCurve
 455 //purpose  :
 456 //=======================================================================
 457
 458 const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const
 459 {
 460   if (!IsDone()) { StdFail_NotDone::Raise(); }
 461   return myF.PointOnCurve(N);
 462 }
 463
 464 //=======================================================================
 465 //function : PointOnSurface
 466 //purpose  :
 467 //=======================================================================
 468
 469 const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const
 470 {
 471   if (!IsDone()) { StdFail_NotDone::Raise(); }
 472   return myF.PointOnSurface(N);
 473 }
 474
 475 //=======================================================================
 476 //function : BidonSurface
 477 //purpose  :
 478 //=======================================================================
 479
 480 Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const
 481 {
 482   return (Adaptor3d_SurfacePtr)0L;
 483 }
 484
 485 //=======================================================================
 486 //function : BidonCurve
 487 //purpose  :
 488 //=======================================================================
 489
 490 Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const
 491 {
 492   return (Adaptor3d_CurvePtr)0L;
 493 }
 494