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
  18 #include <Extrema_GenExtCS.hxx>
  19 #include <Adaptor3d_Curve.hxx>
  20 #include <Adaptor3d_HCurve.hxx>
  21 #include <Adaptor3d_Surface.hxx>
  22 #include <Adaptor3d_HSurface.hxx>
  23 #include <Geom_OffsetCurve.hxx>
  24 #include <Extrema_GlobOptFuncCS.hxx>
  25 #include <Extrema_GlobOptFuncConicS.hxx>
  26 #include <Extrema_GlobOptFuncCQuadric.hxx>
  27 #include <Extrema_POnCurv.hxx>
  28 #include <Extrema_POnSurf.hxx>
  29 #include <Geom_Hyperbola.hxx>
  30 #include <math_FunctionSetRoot.hxx>
  31 #include <math_PSO.hxx>
  32 #include <math_PSOParticlesPool.hxx>
  33 #include <math_Vector.hxx>
  34 #include <Precision.hxx>
  35 #include <Standard_OutOfRange.hxx>
  36 #include <Standard_TypeMismatch.hxx>
  37 #include <StdFail_NotDone.hxx>
  38 #include <TColgp_Array1OfPnt.hxx>
  39 #include <Geom_TrimmedCurve.hxx>
  40 #include <ElCLib.hxx>
  41
  42 const Standard_Real MaxParamVal = 1.0e+10;
  43 const Standard_Real aBorderDivisor = 1.0e+4;
  44 const Standard_Real HyperbolaLimit = 23.; //ln(MaxParamVal)
  45
  46 static Standard_Boolean IsQuadric(const GeomAbs_SurfaceType theSType)
  47 {
  48   if (theSType == GeomAbs_Plane) return Standard_True;
  49   if (theSType == GeomAbs_Cylinder) return Standard_True;
  50   if (theSType == GeomAbs_Cone) return Standard_True;
  51   if (theSType == GeomAbs_Sphere) return Standard_True;
  52   if (theSType == GeomAbs_Torus) return Standard_True;
  53   return  Standard_False;
  54 }
  55 static Standard_Boolean IsConic(const GeomAbs_CurveType theCType)
  56 {
  57   if (theCType == GeomAbs_Line) return Standard_True;
  58   if (theCType == GeomAbs_Circle) return Standard_True;
  59   if (theCType == GeomAbs_Ellipse) return Standard_True;
  60   if (theCType == GeomAbs_Hyperbola) return Standard_True;
  61   if (theCType == GeomAbs_Parabola) return Standard_True;
  62   return  Standard_False;
  63 }
  64 // restrict maximal parameter on hyperbola to avoid FPE
  65 static Standard_Real GetCurvMaxParamVal (const Adaptor3d_Curve& theC)
  66 {
  67   if (theC.GetType() == GeomAbs_Hyperbola)
  68   {
  69     return HyperbolaLimit;
  70   }
  71   if (theC.GetType() == GeomAbs_OffsetCurve)
  72   {
  73     Handle(Geom_Curve) aBC (theC.OffsetCurve()->BasisCurve());
  74     Handle(Geom_TrimmedCurve) aTC = Handle(Geom_TrimmedCurve)::DownCast (aBC);
  75     if (! aTC.IsNull())
  76     {
  77       aBC = aTC->BasisCurve();
  78     }
  79     if (aBC->IsKind (STANDARD_TYPE(Geom_Hyperbola)))
  80       return HyperbolaLimit;
  81   }
  82   return MaxParamVal;
  83 }
  84
  85 // restrict maximal parameter on surfaces based on hyperbola to avoid FPE
  86 static void GetSurfMaxParamVals (const Adaptor3d_Surface& theS,
  87                                  Standard_Real& theUmax, Standard_Real& theVmax)
  88 {
  89   theUmax = theVmax = MaxParamVal;
  90
  91   if (theS.GetType() == GeomAbs_SurfaceOfExtrusion)
  92   {
  93     theUmax = GetCurvMaxParamVal (theS.BasisCurve()->Curve());
  94   }
  95   else if (theS.GetType() == GeomAbs_SurfaceOfRevolution)
  96   {
  97     theVmax = GetCurvMaxParamVal (theS.BasisCurve()->Curve());
  98   }
  99   else if (theS.GetType() == GeomAbs_OffsetSurface)
 100   {
 101     GetSurfMaxParamVals (theS.BasisSurface()->Surface(), theUmax, theVmax);
 102   }
 103 }
 104
 105 //=======================================================================
 106 //function : Extrema_GenExtCS
 107 //purpose  :
 108 //=======================================================================
 109 Extrema_GenExtCS::Extrema_GenExtCS()
 110 : myDone(Standard_False),
 111   mytmin(0.0),
 112   mytsup(0.0),
 113   myumin(0.0),
 114   myusup(0.0),
 115   myvmin(0.0),
 116   myvsup(0.0),
 117   mytsample(0),
 118   myusample(0),
 119   myvsample(0),
 120   mytol1(0.0),
 121   mytol2(0.0),
 122   myS(NULL)
 123 {
 124 }
 125
 126 //=======================================================================
 127 //function : Extrema_GenExtCS
 128 //purpose  :
 129 //=======================================================================
 130 Extrema_GenExtCS::Extrema_GenExtCS(const Adaptor3d_Curve& C,
 131                                    const Adaptor3d_Surface& S,
 132                                    const Standard_Integer NbT,
 133                                    const Standard_Integer NbU,
 134                                    const Standard_Integer NbV,
 135                                    const Standard_Real Tol1,
 136                                    const Standard_Real Tol2)
 137 {
 138   Initialize(S, NbU, NbV, Tol2);
 139   Perform(C, NbT, Tol1);
 140 }
 141
 142 //=======================================================================
 143 //function : Extrema_GenExtCS
 144 //purpose  :
 145 //=======================================================================
 146
 147 Extrema_GenExtCS::Extrema_GenExtCS (const Adaptor3d_Curve& C,
 148                                     const Adaptor3d_Surface& S,
 149                                     const Standard_Integer NbT,
 150                                     const Standard_Integer NbU,
 151                                     const Standard_Integer NbV,
 152                                     const Standard_Real tmin,
 153                                     const Standard_Real tsup,
 154                                     const Standard_Real Umin,
 155                                     const Standard_Real Usup,
 156                                     const Standard_Real Vmin,
 157                                     const Standard_Real Vsup,
 158                                     const Standard_Real Tol1,
 159                                     const Standard_Real Tol2)
 160 {
 161   Initialize(S, NbU, NbV, Umin,Usup,Vmin,Vsup,Tol2);
 162   Perform(C, NbT, tmin, tsup, Tol1);
 163 }
 164
 165 //=======================================================================
 166 //function : Initialize
 167 //purpose  :
 168 //=======================================================================
 169 void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S,
 170                                    const Standard_Integer NbU,
 171                                    const Standard_Integer NbV,
 172                                    const Standard_Real Tol2)
 173 {
 174   myumin = S.FirstUParameter();
 175   myusup = S.LastUParameter();
 176   myvmin = S.FirstVParameter();
 177   myvsup = S.LastVParameter();
 178   Initialize(S,NbU,NbV,myumin,myusup,myvmin,myvsup,Tol2);
 179 }
 180
 181 //=======================================================================
 182 //function : Initialize
 183 //purpose  :
 184 //=======================================================================
 185 void Extrema_GenExtCS::Initialize (const Adaptor3d_Surface& S,
 186                                    const Standard_Integer NbU,
 187                                    const Standard_Integer NbV,
 188                                    const Standard_Real Umin,
 189                                    const Standard_Real Usup,
 190                                    const Standard_Real Vmin,
 191                                    const Standard_Real Vsup,
 192                                    const Standard_Real Tol2)
 193 {
 194   myS = (Adaptor3d_SurfacePtr)&S;
 195   myusample = NbU;
 196   myvsample = NbV;
 197   myumin = Umin;
 198   myusup = Usup;
 199   myvmin = Vmin;
 200   myvsup = Vsup;
 201   mytol2 = Tol2;
 202
 203   Standard_Real umaxpar, vmaxpar;
 204   GetSurfMaxParamVals(*myS, umaxpar, vmaxpar);
 205
 206   if(Precision::IsInfinite (myusup))
 207   {
 208     myusup = umaxpar;
 209   }
 210   if(Precision::IsInfinite (myumin))
 211   {
 212     myumin = -umaxpar;
 213   }
 214   if(Precision::IsInfinite (myvsup))
 215   {
 216     myvsup =  vmaxpar;
 217   }
 218   if(Precision::IsInfinite (myvmin))
 219   {
 220     myvmin = -vmaxpar;
 221   }
 222
 223   Standard_Real du = (myusup - myumin) / aBorderDivisor;
 224   Standard_Real dv = (myvsup - myvmin) / aBorderDivisor;
 225   const Standard_Real aMinU = myumin + du;
 226   const Standard_Real aMinV = myvmin + dv;
 227   const Standard_Real aMaxU = myusup - du;
 228   const Standard_Real aMaxV = myvsup - dv;
 229
 230   const Standard_Real aStepSU = (aMaxU - aMinU) / myusample;
 231   const Standard_Real aStepSV = (aMaxV - aMinV) / myvsample;
 232
 233   mySurfPnts = new TColgp_HArray2OfPnt (0, myusample, 0, myvsample);
 234
 235   Standard_Real aSU = aMinU;
 236   for (Standard_Integer aSUI = 0; aSUI <= myusample; aSUI++, aSU += aStepSU)
 237   {
 238     Standard_Real aSV = aMinV;
 239     for (Standard_Integer aSVI = 0; aSVI <= myvsample; aSVI++, aSV += aStepSV)
 240     {
 241       mySurfPnts->ChangeValue (aSUI, aSVI) = myS->Value (aSU, aSV);
 242     }
 243   }
 244 }
 245
 246 //=======================================================================
 247 //function : Perform
 248 //purpose  :
 249 //=======================================================================
 250 void Extrema_GenExtCS::Perform(const Adaptor3d_Curve& C,
 251   const Standard_Integer NbT,
 252   const Standard_Real Tol1)
 253 {
 254   mytmin = C.FirstParameter();
 255   mytsup = C.LastParameter();
 256   Perform(C, NbT, mytmin, mytsup,Tol1);
 257 }
 258
 259 //=======================================================================
 260 //function : Perform
 261 //purpose  :
 262 //=======================================================================
 263 void Extrema_GenExtCS::Perform (const Adaptor3d_Curve& C,
 264                                 const Standard_Integer NbT,
 265                                 const Standard_Real tmin,
 266                                 const Standard_Real tsup,
 267                                 const Standard_Real Tol1)
 268 {
 269   myDone = Standard_False;
 270   myF.Initialize(C,*myS);
 271   mytmin = tmin;
 272   mytsup = tsup;
 273   mytol1 = Tol1;
 274   mytsample = NbT;
 275   // Modif de lvt pour trimer la surface non pas aux infinis mais  a +/- 10000
 276
 277   Standard_Real trimusup = myusup, trimumin = myumin,trimvsup = myvsup,trimvmin = myvmin;
 278   Standard_Real aCMaxVal = GetCurvMaxParamVal (C);
 279   if (Precision::IsInfinite(mytsup)){
 280     mytsup = aCMaxVal;
 281   }
 282   if (Precision::IsInfinite(mytmin)){
 283     mytmin = -aCMaxVal;
 284   }
 285   //
 286   Standard_Integer aNbVar = 3;
 287   GeomAbs_SurfaceType aSType = myS->GetType();
 288   if (IsQuadric(aSType))
 289   {
 290     aNbVar = 1;
 291   }
 292   else
 293   {
 294     GeomAbs_CurveType aCType = C.GetType();
 295     if (IsConic(aCType))
 296     {
 297       aNbVar = 2;
 298     }
 299   }
 300
 301   math_Vector Tol(1, 3), TUV(1, 3), TUVinf(1, 3), TUVsup(1, 3);
 302   //
 303   Tol(1) = mytol1;
 304   Tol(2) = mytol2;
 305   Tol(3) = mytol2;
 306   //
 307   TUVinf(1) = mytmin;
 308   TUVinf(2) = trimumin;
 309   TUVinf(3) = trimvmin;
 310   //
 311   TUVsup(1) = mytsup;
 312   TUVsup(2) = trimusup;
 313   TUVsup(3) = trimvsup;
 314   //
 315   // Number of particles used in PSO algorithm (particle swarm optimization).
 316   const Standard_Integer aNbParticles = 48;
 317   //
 318   if (aNbVar == 3)
 319   {
 320     GlobMinGenCS(C, aNbParticles, TUVinf, TUVsup, TUV);
 321   }
 322   else if (aNbVar == 2)
 323   {
 324     GlobMinConicS(C, aNbParticles, TUVinf, TUVsup, TUV);
 325   }
 326   else
 327   {
 328     GlobMinCQuadric(C, aNbParticles, TUVinf, TUVsup, TUV);
 329   }
 330
 331   // Find min approximation
 332   math_FunctionSetRoot anA(myF, Tol);
 333   anA.Perform(myF, TUV, TUVinf, TUVsup);
 334
 335   myDone = Standard_True;
 336 }
 337 //=======================================================================
 338 //function : GlobMinGenCS
 339 //purpose  :
 340 //=======================================================================
 341 void Extrema_GenExtCS::GlobMinGenCS(const Adaptor3d_Curve& theC,
 342                                     const Standard_Integer theNbParticles,
 343                                     const math_Vector& theTUVinf,
 344                                     const math_Vector& theTUVsup,
 345                                     math_Vector& theTUV)
 346 {
 347   math_PSOParticlesPool aParticles(theNbParticles, 3);
 348
 349   math_Vector aMinTUV(1, 3);
 350   aMinTUV = theTUVinf + (theTUVsup - theTUVinf) / aBorderDivisor;
 351
 352   math_Vector aMaxTUV(1, 3);
 353   aMaxTUV = theTUVsup - (theTUVsup - theTUVinf) / aBorderDivisor;
 354
 355   Standard_Real aStepCU = (aMaxTUV(1) - aMinTUV(1)) / mytsample;
 356   Standard_Real aStepSU = (aMaxTUV(2) - aMinTUV(2)) / myusample;
 357   Standard_Real aStepSV = (aMaxTUV(3) - aMinTUV(3)) / myvsample;
 358
 359   // Correct number of curve samples in case of low resolution
 360   Standard_Integer aNewCsample = mytsample;
 361   Standard_Real aScaleFactor = 5.0;
 362   Standard_Real aResolutionCU = aStepCU / theC.Resolution(1.0);
 363
 364   Standard_Real aMinResolution = aScaleFactor * Min(aResolutionCU,
 365     Min(aStepSU / myS->UResolution(1.0), aStepSV / myS->VResolution(1.0)));
 366
 367   if (aMinResolution > Epsilon(1.0))
 368   {
 369     if (aResolutionCU > aMinResolution)
 370     {
 371       const Standard_Integer aMaxNbNodes = 50;
 372
 373       aNewCsample = Min(aMaxNbNodes,
 374         RealToInt(mytsample * aResolutionCU / aMinResolution));
 375
 376       aStepCU = (aMaxTUV(1) - aMinTUV(1)) / aNewCsample;
 377     }
 378   }
 379
 380   // Pre-compute curve sample points.
 381   TColgp_Array1OfPnt aCurvPnts(0, aNewCsample);
 382
 383   Standard_Real aCU1 = aMinTUV(1);
 384   for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCU1 += aStepCU)
 385     aCurvPnts.SetValue(aCUI, theC.Value(aCU1));
 386
 387   PSO_Particle* aParticle = aParticles.GetWorstParticle();
 388   // Select specified number of particles from pre-computed set of samples
 389   Standard_Real aSU = aMinTUV(2);
 390   for (Standard_Integer aSUI = 0; aSUI <= myusample; aSUI++, aSU += aStepSU)
 391   {
 392     Standard_Real aSV = aMinTUV(3);
 393     for (Standard_Integer aSVI = 0; aSVI <= myvsample; aSVI++, aSV += aStepSV)
 394     {
 395       Standard_Real aCU2 = aMinTUV(1);
 396       for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCU2 += aStepCU)
 397       {
 398         Standard_Real aSqDist = mySurfPnts->Value(aSUI, aSVI).SquareDistance(aCurvPnts.Value(aCUI));
 399
 400         if (aSqDist < aParticle->Distance)
 401         {
 402           aParticle->Position[0] = aCU2;
 403           aParticle->Position[1] = aSU;
 404           aParticle->Position[2] = aSV;
 405
 406           aParticle->BestPosition[0] = aCU2;
 407           aParticle->BestPosition[1] = aSU;
 408           aParticle->BestPosition[2] = aSV;
 409
 410           aParticle->Distance = aSqDist;
 411           aParticle->BestDistance = aSqDist;
 412
 413           aParticle = aParticles.GetWorstParticle();
 414         }
 415       }
 416     }
 417   }
 418
 419   math_Vector aStep(1, 3);
 420   aStep(1) = aStepCU;
 421   aStep(2) = aStepSU;
 422   aStep(3) = aStepSV;
 423
 424   // Find min approximation
 425   Standard_Real aValue;
 426   Extrema_GlobOptFuncCS aFunc(&theC, myS);
 427   math_PSO aPSO(&aFunc, theTUVinf, theTUVsup, aStep);
 428   aPSO.Perform(aParticles, theNbParticles, aValue, theTUV);
 429
 430 }
 431 //=======================================================================
 432 //function : GlobMinConicS
 433 //purpose  :
 434 //=======================================================================
 435 void Extrema_GenExtCS::GlobMinConicS(const Adaptor3d_Curve& theC,
 436   const Standard_Integer theNbParticles,
 437   const math_Vector& theTUVinf,
 438   const math_Vector& theTUVsup,
 439   math_Vector& theTUV)
 440 {
 441   Standard_Integer aNbVar = 2;
 442   math_Vector  anUVinf(1, aNbVar), anUVsup(1, aNbVar), anUV(1, aNbVar);
 443   Standard_Integer i;
 444   for (i = 1; i <= aNbVar; ++i)
 445   {
 446     anUVinf(i) = theTUVinf(i + 1);
 447     anUVsup(i) = theTUVsup(i + 1);
 448   }
 449   //
 450   math_PSOParticlesPool aParticles(theNbParticles, aNbVar);
 451
 452   math_Vector aMinUV(1, aNbVar);
 453   aMinUV = anUVinf + (anUVsup - anUVinf) / aBorderDivisor;
 454
 455   math_Vector aMaxUV(1, aNbVar);
 456   aMaxUV = anUVsup - (anUVsup - anUVinf) / aBorderDivisor;
 457
 458   //Increase numbers of UV samples to improve searching global minimum
 459   Standard_Integer anAddsample = Max(mytsample / 2, 3);
 460   Standard_Integer anUsample = myusample + anAddsample;
 461   Standard_Integer aVsample = myvsample + anAddsample;
 462   //
 463   Standard_Real aStepSU = (aMaxUV(1) - aMinUV(1)) / anUsample;
 464   Standard_Real aStepSV = (aMaxUV(2) - aMinUV(2)) / aVsample;
 465   //
 466   Extrema_GlobOptFuncConicS aFunc(myS, anUVinf(1), anUVsup(1), anUVinf(2), anUVsup(2));
 467   aFunc.LoadConic(&theC, theTUVinf(1), theTUVsup(1));
 468
 469
 470   PSO_Particle* aParticle = aParticles.GetWorstParticle();
 471   // Select specified number of particles from pre-computed set of samples
 472   Standard_Real aSU = aMinUV(1);
 473
 474   for (Standard_Integer aSUI = 0; aSUI <= anUsample; aSUI++, aSU += aStepSU)
 475   {
 476     anUV(1) = aSU;
 477     Standard_Real aSV = aMinUV(2);
 478     for (Standard_Integer aSVI = 0; aSVI <= aVsample; aSVI++, aSV += aStepSV)
 479     {
 480       anUV(2) = aSV;
 481       Standard_Real aSqDist;
 482       if (!aFunc.Value(anUV, aSqDist))
 483       {
 484         aSqDist = Precision::Infinite();
 485       }
 486
 487       if (aSqDist < aParticle->Distance)
 488       {
 489         aParticle->Position[0] = aSU;
 490         aParticle->Position[1] = aSV;
 491
 492         aParticle->BestPosition[0] = aSU;
 493         aParticle->BestPosition[1] = aSV;
 494
 495         aParticle->Distance = aSqDist;
 496         aParticle->BestDistance = aSqDist;
 497
 498         aParticle = aParticles.GetWorstParticle();
 499       }
 500     }
 501   }
 502
 503   math_Vector aStep(1, aNbVar);
 504   aStep(1) = aStepSU;
 505   aStep(2) = aStepSV;
 506
 507   // Find min approximation
 508   Standard_Real aValue;
 509   math_PSO aPSO(&aFunc, anUVinf, anUVsup, aStep);
 510   aPSO.Perform(aParticles, theNbParticles, aValue, anUV);
 511   //
 512   Standard_Real aCT = aFunc.ConicParameter(anUV);
 513   if (theC.IsPeriodic())
 514   {
 515     if (aCT <  theTUVinf(1) - Precision::PConfusion() || aCT >  theTUVsup(1) + Precision::PConfusion())
 516     {
 517       aCT = ElCLib::InPeriod(aCT, theTUVinf(1), theTUVinf(1) + 2. * M_PI);
 518     }
 519   }
 520   theTUV(1) = aCT;
 521   theTUV(2) = anUV(1);
 522   theTUV(3) = anUV(2);
 523
 524 }
 525 //=======================================================================
 526 //function : GlobMinCQuadric
 527 //purpose  :
 528 //=======================================================================
 529 void Extrema_GenExtCS::GlobMinCQuadric(const Adaptor3d_Curve& theC,
 530   const Standard_Integer theNbParticles,
 531   const math_Vector& theTUVinf,
 532   const math_Vector& theTUVsup,
 533   math_Vector& theTUV)
 534 {
 535   Standard_Integer aNbVar = 1;
 536   math_Vector aTinf(1, aNbVar), aTsup(1, aNbVar), aT(1, aNbVar);
 537   aTinf(1) = theTUVinf(1);
 538   aTsup(1) = theTUVsup(1);
 539   //
 540   math_PSOParticlesPool aParticles(theNbParticles, aNbVar);
 541
 542   math_Vector aMinT(1, aNbVar);
 543   aMinT = aTinf + (aTsup - aTinf) / aBorderDivisor;
 544
 545   math_Vector aMaxT(1, aNbVar);
 546   aMaxT = aTsup - (aTsup - aTinf) / aBorderDivisor;
 547   //
 548
 549   //Increase numbers of curve samples to improve searching global minimum
 550   //because dimension of optimisation task is redused
 551   const Standard_Integer aMaxNbNodes = 50;
 552   Standard_Integer aNewCsample = mytsample;
 553   Standard_Integer anAddsample = Max(myusample / 2, 3);
 554   aNewCsample += anAddsample;
 555   aNewCsample = Min(aNewCsample, aMaxNbNodes);
 556   //
 557   // Correct number of curve samples in case of low resolution
 558   Standard_Real aStepCT = (aMaxT(1) - aMinT(1)) / aNewCsample;
 559   Standard_Real aStepSU = (theTUVsup(2) - theTUVinf(2)) / myusample;
 560   Standard_Real aStepSV = (theTUVsup(3) - theTUVinf(3)) / myvsample;
 561   Standard_Real aScaleFactor = 5.0;
 562   Standard_Real aResolutionCU = aStepCT / theC.Resolution(1.0);
 563
 564   Standard_Real aMinResolution = aScaleFactor * Min(aResolutionCU,
 565     Min(aStepSU / myS->UResolution(1.0), aStepSV / myS->VResolution(1.0)));
 566
 567   if (aMinResolution > Epsilon(1.0))
 568   {
 569     if (aResolutionCU > aMinResolution)
 570     {
 571
 572       aNewCsample = Min(aMaxNbNodes,
 573         RealToInt(aNewCsample * aResolutionCU / aMinResolution));
 574
 575       aStepCT = (aMaxT(1) - aMinT(1)) / aNewCsample;
 576     }
 577   }
 578
 579   //
 580   Extrema_GlobOptFuncCQuadric aFunc(&theC, aTinf(1), aTsup(1));
 581   aFunc.LoadQuad(myS, theTUVinf(2), theTUVsup(2), theTUVinf(3), theTUVsup(3));
 582
 583   PSO_Particle* aParticle = aParticles.GetWorstParticle();
 584   // Select specified number of particles from pre-computed set of samples
 585   Standard_Real aCT = aMinT(1);
 586   for (Standard_Integer aCUI = 0; aCUI <= aNewCsample; aCUI++, aCT += aStepCT)
 587   {
 588     aT(1) = aCT;
 589     Standard_Real aSqDist;
 590     if (!aFunc.Value(aT, aSqDist))
 591     {
 592       aSqDist = Precision::Infinite();
 593     }
 594
 595     if (aSqDist < aParticle->Distance)
 596     {
 597       aParticle->Position[0] = aCT;
 598
 599       aParticle->BestPosition[0] = aCT;
 600
 601       aParticle->Distance = aSqDist;
 602       aParticle->BestDistance = aSqDist;
 603
 604       aParticle = aParticles.GetWorstParticle();
 605     }
 606   }
 607   //
 608   math_Vector aStep(1, aNbVar);
 609   aStep(1) = aStepCT;
 610
 611   // Find min approximation
 612   Standard_Real aValue;
 613   math_PSO aPSO(&aFunc, aTinf, aTsup, aStep);
 614   aPSO.Perform(aParticles, theNbParticles, aValue, aT);
 615   //
 616   math_Vector anUV(1, 2);
 617   aFunc.QuadricParameters(aT, anUV);
 618   if (myS->IsUPeriodic())
 619   {
 620     if (anUV(1) <  theTUVinf(2) - Precision::PConfusion() || anUV(1)>  theTUVsup(2) + Precision::PConfusion())
 621     {
 622       anUV(1) = ElCLib::InPeriod(anUV(1), theTUVinf(2), theTUVinf(2) + 2. * M_PI);
 623     }
 624   }
 625   //
 626   if (myS->IsVPeriodic())
 627   {
 628     if (anUV(2) <  theTUVinf(3) - Precision::PConfusion() || anUV(2)>  theTUVsup(3) + Precision::PConfusion())
 629     {
 630       anUV(2) = ElCLib::InPeriod(anUV(2), theTUVinf(3), theTUVinf(3) + 2. * M_PI);
 631     }
 632   }
 633   //
 634   theTUV(1) = aT(1);
 635   theTUV(2) = anUV(1);
 636   theTUV(3) = anUV(2);
 637
 638 }
 639
 640 //=======================================================================
 641 //function : IsDone
 642 //purpose  :
 643 //=======================================================================
 644 Standard_Boolean Extrema_GenExtCS::IsDone() const
 645 {
 646   return myDone;
 647 }
 648
 649 //=======================================================================
 650 //function : NbExt
 651 //purpose  :
 652 //=======================================================================
 653 Standard_Integer Extrema_GenExtCS::NbExt() const
 654 {
 655   if (!IsDone()) { throw StdFail_NotDone(); }
 656   return myF.NbExt();
 657 }
 658
 659 //=======================================================================
 660 //function : Value
 661 //purpose  :
 662 //=======================================================================
 663 Standard_Real Extrema_GenExtCS::SquareDistance(const Standard_Integer N) const
 664 {
 665   if (N < 1 || N > NbExt())
 666   {
 667     throw Standard_OutOfRange();
 668   }
 669
 670   return myF.SquareDistance(N);
 671 }
 672
 673 //=======================================================================
 674 //function : PointOnCurve
 675 //purpose  :
 676 //=======================================================================
 677 const Extrema_POnCurv& Extrema_GenExtCS::PointOnCurve(const Standard_Integer N) const
 678 {
 679   if (N < 1 || N > NbExt())
 680   {
 681     throw Standard_OutOfRange();
 682   }
 683
 684   return myF.PointOnCurve(N);
 685 }
 686
 687 //=======================================================================
 688 //function : PointOnSurface
 689 //purpose  :
 690 //=======================================================================
 691 const Extrema_POnSurf& Extrema_GenExtCS::PointOnSurface(const Standard_Integer N) const
 692 {
 693   if (N < 1 || N > NbExt())
 694   {
 695     throw Standard_OutOfRange();
 696   }
 697
 698   return myF.PointOnSurface(N);
 699 }
 700
 701 //=======================================================================
 702 //function : BidonSurface
 703 //purpose  :
 704 //=======================================================================
 705 Adaptor3d_SurfacePtr Extrema_GenExtCS::BidonSurface() const
 706 {
 707   return (Adaptor3d_SurfacePtr)0L;
 708 }
 709
 710 //=======================================================================
 711 //function : BidonCurve
 712 //purpose  :
 713 //=======================================================================
 714 Adaptor3d_CurvePtr Extrema_GenExtCS::BidonCurve() const
 715 {
 716   return (Adaptor3d_CurvePtr)0L;
 717 }
 718