src/GCPnts/GCPnts_TangentialDeflection.gxx

   1 // Created on: 1996-11-08
   2 // Created by: Jean Claude VAUTHIER
   3 // Copyright (c) 1996-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 <GCPnts_DeflectionType.hxx>
  18 #include <Standard_ConstructionError.hxx>
  19 #include <Precision.hxx>
  20 #include <gp_XYZ.hxx>
  21 #include <gp_Pnt.hxx>
  22 #include <gp_Vec.hxx>
  23 #include <gp.hxx>
  24 #include <NCollection_List.hxx>
  25 #include <math_PSO.hxx>
  26 #include <math_BrentMinimum.hxx>
  27
  28 #define Us3 0.3333333333333333333333333333
  29
  30 void GCPnts_TangentialDeflection::EvaluateDu (
  31   const TheCurve& C,
  32   const Standard_Real      U,
  33   gp_Pnt&   P,
  34   Standard_Real&     Du,
  35   Standard_Boolean&  NotDone) const {
  36
  37   gp_Vec T, N;
  38   D2 (C, U, P, T, N);
  39   Standard_Real Lt   = T.Magnitude ();
  40   Standard_Real LTol = Precision::Confusion ();
  41   if (Lt > LTol && N.Magnitude () > LTol) {
  42     Standard_Real Lc = N.CrossMagnitude (T);
  43     Standard_Real Ln = Lc/Lt;
  44     if (Ln > LTol) {
  45       Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
  46       NotDone = Standard_False;
  47     }
  48   }
  49 }
  50
  51
  52 //=======================================================================
  53 //function : GCPnts_TangentialDeflection
  54 //purpose  :
  55 //=======================================================================
  56
  57 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  58  const TheCurve&        C,
  59  const Standard_Real    AngularDeflection,
  60  const Standard_Real    CurvatureDeflection,
  61  const Standard_Integer MinimumOfPoints,
  62  const Standard_Real    UTol,
  63  const Standard_Real    theMinLen)
  64
  65 {
  66   Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen);
  67 }
  68
  69
  70 //=======================================================================
  71 //function : GCPnts_TangentialDeflection
  72 //purpose  :
  73 //=======================================================================
  74
  75 GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
  76  const TheCurve&        C,
  77  const Standard_Real    FirstParameter,
  78  const Standard_Real    LastParameter,
  79  const Standard_Real    AngularDeflection,
  80  const Standard_Real    CurvatureDeflection,
  81  const Standard_Integer MinimumOfPoints,
  82  const Standard_Real    UTol,
  83  const Standard_Real    theMinLen)
  84
  85 {
  86   Initialize (C,
  87         FirstParameter,
  88         LastParameter,
  89         AngularDeflection,
  90         CurvatureDeflection,
  91         MinimumOfPoints,
  92         UTol, theMinLen);
  93 }
  94
  95
  96
  97 //=======================================================================
  98 //function : Initialize
  99 //purpose  :
 100 //=======================================================================
 101
 102 void GCPnts_TangentialDeflection::Initialize (
 103  const TheCurve&        C,
 104  const Standard_Real    AngularDeflection,
 105  const Standard_Real    CurvatureDeflection,
 106  const Standard_Integer MinimumOfPoints,
 107  const Standard_Real    UTol,
 108  const Standard_Real    theMinLen)
 109
 110 {
 111   Initialize (C,
 112               C.FirstParameter (),
 113               C.LastParameter (),
 114               AngularDeflection,
 115               CurvatureDeflection,
 116               MinimumOfPoints,
 117               UTol, theMinLen);
 118 }
 119
 120
 121 //=======================================================================
 122 //function : Initialize
 123 //purpose  :
 124 //=======================================================================
 125
 126 void GCPnts_TangentialDeflection::Initialize (
 127  const TheCurve&        C,
 128  const Standard_Real    FirstParameter,
 129  const Standard_Real    LastParameter,
 130  const Standard_Real    AngularDeflection,
 131  const Standard_Real    CurvatureDeflection,
 132  const Standard_Integer MinimumOfPoints,
 133  const Standard_Real    UTol,
 134  const Standard_Real    theMinLen)
 135
 136 {
 137
 138   Standard_ConstructionError_Raise_if (CurvatureDeflection <= Precision::Confusion () || AngularDeflection   <= Precision::Angular (), "GCPnts_TangentialDeflection::Initialize - Zero Deflection")
 139
 140  parameters.Clear ();
 141  points    .Clear ();
 142  if (FirstParameter < LastParameter) {
 143    firstu = FirstParameter;
 144    lastu  = LastParameter;
 145  }
 146  else {
 147    lastu  = FirstParameter;
 148    firstu = LastParameter;
 149  }
 150  uTol                = UTol;
 151  angularDeflection   = AngularDeflection;
 152  curvatureDeflection = CurvatureDeflection;
 153  minNbPnts           = Max (MinimumOfPoints, 2);
 154  myMinLen            = Max(theMinLen, Precision::Confusion());
 155
 156  switch (C.GetType()) {
 157
 158  case GeomAbs_Line:
 159    PerformLinear (C);
 160    break;
 161
 162  case GeomAbs_Circle:
 163    PerformCircular (C);
 164    break;
 165
 166  case GeomAbs_BSplineCurve:
 167    {
 168      Handle_TheBSplineCurve BS = C.BSpline() ;
 169      if (BS->NbPoles() == 2 ) PerformLinear (C);
 170      else                     PerformCurve  (C);
 171      break;
 172    }
 173  case GeomAbs_BezierCurve:
 174    {
 175      Handle_TheBezierCurve  BZ = C.Bezier();
 176      if (BZ->NbPoles() == 2) PerformLinear (C);
 177      else                    PerformCurve  (C);
 178      break;
 179    }
 180  default: PerformCurve (C);
 181
 182  }
 183 }
 184
 185
 186 //=======================================================================
 187 //function : PerformLinear
 188 //purpose  :
 189 //=======================================================================
 190
 191 void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) {
 192
 193   gp_Pnt P;
 194   D0 (C, firstu, P);
 195   parameters.Append (firstu);
 196   points    .Append (P);
 197   if (minNbPnts > 2) {
 198     Standard_Real Du = (lastu - firstu) / minNbPnts;
 199     Standard_Real U = firstu + Du;
 200     for (Standard_Integer i = 2; i <= minNbPnts; i++) {
 201       D0 (C, U, P);
 202       parameters.Append (U);
 203       points    .Append (P);
 204       U += Du;
 205     }
 206   }
 207   D0 (C, lastu, P);
 208   parameters.Append (lastu);
 209   points    .Append (P);
 210 }
 211
 212 //=======================================================================
 213 //function : PerformCircular
 214 //purpose  :
 215 //=======================================================================
 216
 217 void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C)
 218 {
 219   // akm 8/01/02 : check the radius before divide by it
 220   Standard_Real dfR = C.Circle().Radius();
 221   Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
 222     dfR, curvatureDeflection, angularDeflection, myMinLen);
 223
 224   const Standard_Real aDiff = lastu - firstu;
 225   // Round up number of points to satisfy curvatureDeflection more precisely
 226   Standard_Integer NbPoints = (Standard_Integer)Ceiling(aDiff / Du);
 227   NbPoints = Max(NbPoints, minNbPnts - 1);
 228   Du       = aDiff / NbPoints;
 229
 230   gp_Pnt P;
 231   Standard_Real U = firstu;
 232   for (Standard_Integer i = 1; i <= NbPoints; i++)
 233   {
 234     D0 (C, U, P);
 235     parameters.Append (U);
 236     points    .Append (P);
 237     U += Du;
 238   }
 239   D0 (C, lastu, P);
 240   parameters.Append (lastu);
 241   points    .Append (P);
 242 }
 243
 244
 245 //=======================================================================
 246 //function : PerformCurve
 247 //purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
 248 //           minimum de points sur un element lineaire
 249 //=======================================================================
 250 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
 251
 252 {
 253   Standard_Integer i, j;
 254   gp_XYZ V1, V2;
 255   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;
 256   Standard_Real Du, Dusave, MiddleU, L1, L2;
 257
 258   Standard_Real U1       = firstu;
 259   Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
 260   Standard_Real ATol     = 1.e-2 * angularDeflection;
 261   if(ATol > 1.e-2)
 262     ATol = 1.e-2;
 263   else if(ATol < 1.e-7)
 264     ATol = 1.e-7;
 265
 266   D0 (C, lastu, LastPoint);
 267
 268   //Initialization du calcul
 269
 270   Standard_Boolean NotDone = Standard_True;
 271   Dusave = (lastu - firstu)*Us3;
 272   Du     = Dusave;
 273   EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
 274   parameters.Append (U1);
 275   points    .Append (CurrentPoint);
 276
 277   // Used to detect "isLine" current bspline and in Du computation in general handling.
 278   Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
 279   TColStd_Array1OfReal Intervs(1, NbInterv+1);
 280   C.Intervals(Intervs, GeomAbs_CN);
 281
 282   if (NotDone || Du > 5. * Dusave) {
 283     //C'est soit une droite, soit une singularite :
 284     V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
 285     L1 = V1.Modulus ();
 286     if (L1 > LTol)
 287     {
 288       //Si c'est une droite on verifie en calculant minNbPoints :
 289       Standard_Boolean IsLine   = Standard_True;
 290       Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
 291       switch (C.GetType()) {
 292       case GeomAbs_BSplineCurve:
 293         {
 294           Handle_TheBSplineCurve BS = C.BSpline() ;
 295           NbPoints = Max(BS->Degree() + 1, NbPoints);
 296           break;
 297         }
 298       case GeomAbs_BezierCurve:
 299         {
 300           Handle_TheBezierCurve  BZ = C.Bezier();
 301           NbPoints = Max(BZ->Degree() + 1, NbPoints);
 302           break;
 303         }
 304       default:
 305       ;}
 306       ////
 307       Standard_Real param = 0.;
 308       for (i = 1; i <= NbInterv && IsLine; ++i)
 309       {
 310         // Avoid usage intervals out of [firstu, lastu].
 311         if ((Intervs(i+1) < firstu) ||
 312             (Intervs(i)   > lastu))
 313         {
 314           continue;
 315         }
 316         // Fix border points in applicable intervals, to avoid be out of target interval.
 317         if ((Intervs(i)   < firstu) &&
 318             (Intervs(i+1) > firstu))
 319         {
 320           Intervs(i) = firstu;
 321         }
 322         if ((Intervs(i)   < lastu) &&
 323             (Intervs(i+1) > lastu))
 324         {
 325           Intervs(i + 1) = lastu;
 326         }
 327
 328         Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
 329         for (j = 1; j <= NbPoints && IsLine; ++j)
 330         {
 331           param = Intervs(i) + j*delta;
 332           D0 (C, param, MiddlePoint);
 333           V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 334           L2 = V2.Modulus ();
 335           if (L2 > LTol)
 336           {
 337             const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
 338             IsLine = (aAngle < ATol);
 339           }
 340         }
 341       }
 342
 343       if (IsLine)
 344       {
 345         parameters.Clear();
 346         points    .Clear();
 347
 348         PerformLinear(C);
 349         return;
 350       }
 351       else
 352       {
 353         //c'etait une singularite on continue :
 354         //Du = Dusave;
 355         EvaluateDu (C, param, MiddlePoint, Du, NotDone);
 356       }
 357     }
 358     else
 359     {
 360       Du = (lastu-firstu)/2.1;
 361       MiddleU = firstu + Du;
 362       D0 (C, MiddleU, MiddlePoint);
 363       V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
 364       L1 = V1.Modulus ();
 365       if (L1 < LTol)
 366       {
 367         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
 368         // on renvoi un segment de 2 points   (protection)
 369         parameters.Append (lastu);
 370         points    .Append (LastPoint);
 371         return;
 372       }
 373     }
 374   }
 375
 376   if (Du > Dusave) Du = Dusave;
 377   else             Dusave = Du;
 378
 379   if (Du < uTol) {
 380     Du = lastu - firstu;
 381     if (Du < uTol) {
 382       parameters.Append (lastu);
 383       points    .Append (LastPoint);
 384       return;
 385     }
 386   }
 387
 388   //Traitement normal pour une courbe
 389   Standard_Boolean MorePoints = Standard_True;
 390   Standard_Real U2            = firstu;
 391   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
 392   Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
 393   Standard_Boolean isNeedToCheck = Standard_False;
 394   gp_Pnt aPrevPoint = points.Last();
 395
 396   while (MorePoints) {
 397     aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
 398     U2 += Du;
 399
 400     if (U2 >= lastu) {                       //Bout de courbe
 401       U2 = lastu;
 402       CurrentPoint = LastPoint;
 403       Du = U2-U1;
 404       Dusave = Du;
 405     }
 406     else D0 (C, U2, CurrentPoint);           //Point suivant
 407
 408     Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
 409     Standard_Boolean Correction, TooLarge, TooSmall;
 410     TooLarge   = Standard_False;
 411     Correction = Standard_True;
 412     TooSmall = Standard_False;
 413
 414     while (Correction) {                     //Ajustement Du
 415       if (isNeedToCheck)
 416       {
 417         aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
 418         if (aIdx[1] > aIdx[0]) // Jump to another polynom.
 419         {
 420           if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
 421           {
 422             Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
 423             U2 = U1 + Du;
 424             if (U2 > lastu)
 425               U2 = lastu;
 426             D0 (C, U2, CurrentPoint);
 427           }
 428         }
 429       }
 430       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
 431       D0 (C, MiddleU, MiddlePoint);
 432
 433       V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
 434       V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
 435       L1 = V1.Modulus ();
 436
 437       FCoef = (L1 > myMinLen) ?
 438         V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
 439
 440       V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
 441       L1 = V1.Modulus ();
 442       L2 = V2.Modulus ();
 443       if (L1 > myMinLen && L2 > myMinLen)
 444       {
 445         Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
 446         ACoef = angg / AngleMax;
 447       }
 448       else
 449         ACoef = 0.0;
 450
 451       //On retient le plus penalisant
 452       Coef = Max(ACoef, FCoef);
 453
 454       if (isNeedToCheck && Coef < 0.55)
 455       {
 456         isNeedToCheck = Standard_False;
 457         Du = Dusave;
 458         U2 = U1 + Du;
 459         if (U2 > lastu)
 460           U2 = lastu;
 461         D0 (C, U2, CurrentPoint);
 462         continue;
 463       }
 464
 465       if (Coef <= 1.0) {
 466         if (Abs (lastu-U2) < uTol) {
 467           parameters.Append (lastu);
 468           points    .Append (LastPoint);
 469           MorePoints = Standard_False;
 470           Correction = Standard_False;
 471         }
 472         else {
 473           if (Coef >= 0.55 || TooLarge) {
 474             parameters.Append (U2);
 475             points    .Append (CurrentPoint);
 476             aPrevPoint = CurrentPoint;
 477             Correction = Standard_False;
 478             isNeedToCheck = Standard_True;
 479           }
 480           else if (TooSmall) {
 481             Correction = Standard_False;
 482             aPrevPoint = CurrentPoint;
 483           }
 484           else {
 485             TooSmall = Standard_True;
 486             //Standard_Real UUU2 = U2;
 487             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
 488
 489             U2 = U1 + Du;
 490             if (U2 > lastu)
 491               U2 = lastu;
 492             D0 (C, U2, CurrentPoint);
 493           }
 494         }
 495       }
 496       else {
 497
 498         if (Coef >= 1.5) {
 499           if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion()))
 500           {
 501             parameters.Append (U1);
 502             points    .Append (aPrevPoint);
 503           }
 504           U2 = MiddleU;
 505           Du  = U2-U1;
 506           CurrentPoint = MiddlePoint;
 507         }
 508         else {
 509           Du*=0.9;
 510           U2 = U1 + Du;
 511           D0 (C, U2, CurrentPoint);
 512           TooLarge = Standard_True;
 513         }
 514
 515       }
 516     }
 517
 518     Du  = U2-U1;
 519
 520     if (MorePoints) {
 521       if (U1 > firstu) {
 522         if (FCoef > ACoef) {
 523           //La fleche est critere de decoupage
 524           EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
 525           if (NotDone) {
 526             Du += (Du-Dusave)*(Du/Dusave);
 527             if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 528             if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 529           }
 530         }
 531         else {
 532           //L'angle est le critere de decoupage
 533           Du += (Du-Dusave)*(Du/Dusave);
 534           if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
 535           if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
 536         }
 537       }
 538
 539       if (Du < uTol) {
 540         Du = lastu - U2;
 541         if (Du < uTol) {
 542           parameters.Append (lastu);
 543           points    .Append (LastPoint);
 544           MorePoints = Standard_False;
 545         }
 546         else if (Du*Us3 > uTol) Du*=Us3;
 547       }
 548       U1 = U2;
 549       Dusave = Du;
 550     }
 551   }
 552   //Recalage avant dernier point :
 553   i = points.Length()-1;
 554 //  Real d = points (i).Distance (points (i+1));
 555 // if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
 556 //    cout<<"deux points confondus"<<endl;
 557 //    parameters.Remove (i+1);
 558 //    points.Remove (i+1);
 559 //    i--;
 560 //  }
 561   if (i >= 2) {
 562     MiddleU = parameters (i-1);
 563     MiddleU = (lastu + MiddleU)*0.5;
 564     D0 (C, MiddleU, MiddlePoint);
 565     parameters.SetValue (i, MiddleU);
 566     points    .SetValue (i, MiddlePoint);
 567   }
 568
 569   //-- On rajoute des points aux milieux des segments si le nombre
 570   //-- mini de points n'est pas atteint
 571   //--
 572   Standard_Integer Nbp =  points.Length();
 573   Standard_Integer MinNb= (9*minNbPnts)/10;
 574   //if(MinNb<4) MinNb=4;
 575
 576   //-- if(Nbp <  MinNb) { cout<<"\n*"; } else {  cout<<"\n."; }
 577   while(Nbp < MinNb) {
 578     //-- cout<<" \nGCPnts TangentialDeflection : Ajout de Points ("<<Nbp<<" "<<minNbPnts<<" )"<<endl;
 579     for(i=2; i<=Nbp; i++) {
 580       MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
 581       D0 (C, MiddleU, MiddlePoint);
 582       parameters.InsertBefore(i,MiddleU);
 583       points.InsertBefore(i,MiddlePoint);
 584       Nbp++;
 585       i++;
 586     }
 587   }
 588   //Additional check for intervals
 589   Standard_Real MinLen2 = myMinLen * myMinLen;
 590   Standard_Integer MaxNbp = 10 * Nbp;
 591   for(i = 1; i < Nbp; ++i)
 592   {
 593     U1 = parameters(i);
 594     U2 = parameters(i + 1);
 595     // Check maximal deflection on interval;
 596     Standard_Real dmax = 0.;
 597     Standard_Real umax = 0.;
 598     Standard_Real amax = 0.;
 599     EstimDefl(C, U1, U2, dmax, umax);
 600     const gp_Pnt& P1 = points(i);
 601     const gp_Pnt& P2 = points(i+1);
 602     D0(C, umax, MiddlePoint);
 603     amax = EstimAngl(P1, MiddlePoint, P2);
 604     if(dmax > curvatureDeflection || amax > AngleMax)
 605     {
 606       if(umax - U1 > uTol && U2 - umax > uTol)
 607       {
 608         if (P1.SquareDistance(MiddlePoint) > MinLen2 && P2.SquareDistance(MiddlePoint) > MinLen2)
 609         {
 610           parameters.InsertAfter(i, umax);
 611           points.InsertAfter(i, MiddlePoint);
 612           ++Nbp;
 613           --i; //To compensate ++i in loop header: i must point to first part of splitted interval
 614           if(Nbp > MaxNbp)
 615           {
 616             break;
 617           }
 618         }
 619       }
 620     }
 621   }
 622   //
 623 }
 624
 625 //=======================================================================
 626 //function : EstimDefl
 627 //purpose  : Estimation of maximal deflection for interval [U1, U2]
 628 //
 629 //=======================================================================
 630 void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& C,
 631                             const Standard_Real U1, const Standard_Real U2,
 632                             Standard_Real& MaxDefl, Standard_Real& UMax)
 633 {
 634   Standard_Real Du = (lastu - firstu);
 635   //
 636   TheMaxCurvLinDist aFunc(C, U1, U2);
 637   //
 638   const Standard_Integer aNbIter = 100;
 639   Standard_Real reltol = Max(1.e-3, 2.*uTol/((Abs(U1) + Abs(U2))));
 640   //
 641   math_BrentMinimum anOptLoc(reltol, aNbIter, uTol);
 642   anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2);
 643   if(anOptLoc.IsDone())
 644   {
 645     MaxDefl = Sqrt(-anOptLoc.Minimum());
 646     UMax = anOptLoc.Location();
 647     return;
 648   }
 649   //
 650   math_Vector aLowBorder(1,1);
 651   math_Vector aUppBorder(1,1);
 652   math_Vector aSteps(1,1);
 653   //
 654   aSteps(1) = Max(0.1 * Du, 100. * uTol);
 655   Standard_Integer aNbParticles = Max(8, RealToInt(32 * (U2 - U1) / Du));
 656   //
 657   aLowBorder(1) = U1;
 658   aUppBorder(1) = U2;
 659   //
 660   //
 661   Standard_Real aValue;
 662   math_Vector aT(1,1);
 663   TheMaxCurvLinDistMV aFuncMV(aFunc);
 664
 665   math_PSO aFinder(&aFuncMV, aLowBorder, aUppBorder, aSteps, aNbParticles);
 666   aFinder.Perform(aSteps, aValue, aT);
 667   //
 668   anOptLoc.Perform(aFunc, Max(aT(1) - aSteps(1), U1) , aT(1), Min(aT(1) + aSteps(1),U2));
 669   if(anOptLoc.IsDone())
 670   {
 671     MaxDefl = Sqrt(-anOptLoc.Minimum());
 672     UMax = anOptLoc.Location();
 673     return;
 674   }
 675   MaxDefl = Sqrt(-aValue);
 676   UMax = aT(1);
 677   //
 678 }