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