[occt.git] / src / GCPnts / GCPnts_TangentialDeflection.pxx

// Created on: 1996-11-08
// Created by: Jean Claude VAUTHIER
// Copyright (c) 1996-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

#include <GCPnts_DeflectionType.hxx>
#include <Standard_ConstructionError.hxx>
#include <Precision.hxx>
#include <gp_XYZ.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <gp.hxx>
#include <NCollection_List.hxx>
#include <math_PSO.hxx>
#include <math_BrentMinimum.hxx>

#define Us3 0.3333333333333333333333333333

void GCPnts_TangentialDeflection::EvaluateDu (
  const TheCurve& C,
  const Standard_Real      U,
  gp_Pnt&   P,
  Standard_Real&     Du,
  Standard_Boolean&  NotDone) const {

  gp_Vec T, N;
  D2 (C, U, P, T, N);
  Standard_Real Lt   = T.Magnitude ();
  Standard_Real LTol = Precision::Confusion ();
  if (Lt > LTol && N.Magnitude () > LTol) {
    Standard_Real Lc = N.CrossMagnitude (T);
    Standard_Real Ln = Lc/Lt;
    if (Ln > LTol) {
      Du = sqrt (8.0 * Max(curvatureDeflection, myMinLen) / Ln);
      NotDone = Standard_False;
    }
  }
}


//=======================================================================
//function : GCPnts_TangentialDeflection
//purpose  : 
//=======================================================================

GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
 const TheCurve&        C,
 const Standard_Real    AngularDeflection,
 const Standard_Real    CurvatureDeflection,
 const Standard_Integer MinimumOfPoints,
 const Standard_Real    UTol,
 const Standard_Real    theMinLen)

{ 
  Initialize (C,AngularDeflection,CurvatureDeflection,MinimumOfPoints,UTol,theMinLen); 
}


//=======================================================================
//function : GCPnts_TangentialDeflection
//purpose  : 
//=======================================================================

GCPnts_TangentialDeflection::GCPnts_TangentialDeflection (
 const TheCurve&        C,
 const Standard_Real    FirstParameter,
 const Standard_Real    LastParameter,
 const Standard_Real    AngularDeflection,
 const Standard_Real    CurvatureDeflection,
 const Standard_Integer MinimumOfPoints,
 const Standard_Real    UTol,
 const Standard_Real    theMinLen)

{ 
  Initialize (C, 
        FirstParameter, 
        LastParameter,
        AngularDeflection, 
        CurvatureDeflection, 
        MinimumOfPoints,
        UTol, theMinLen);
}


//=======================================================================
//function : Initialize
//purpose  : 
//=======================================================================

void GCPnts_TangentialDeflection::Initialize (
 const TheCurve&        C, 
 const Standard_Real    AngularDeflection, 
 const Standard_Real    CurvatureDeflection,
 const Standard_Integer MinimumOfPoints,
 const Standard_Real    UTol,
 const Standard_Real    theMinLen)

{
  Initialize (C, 
              C.FirstParameter (), 
              C.LastParameter (),
              AngularDeflection, 
              CurvatureDeflection,
              MinimumOfPoints,
              UTol, theMinLen);
}


//=======================================================================
//function : Initialize
//purpose  : 
//=======================================================================

void GCPnts_TangentialDeflection::Initialize (
 const TheCurve&        C, 
 const Standard_Real    FirstParameter, 
 const Standard_Real    LastParameter,
 const Standard_Real    AngularDeflection, 
 const Standard_Real    CurvatureDeflection,
 const Standard_Integer MinimumOfPoints,
 const Standard_Real    UTol,
 const Standard_Real    theMinLen)

{
  
  Standard_ConstructionError_Raise_if (CurvatureDeflection < Precision::Confusion () ||
                                       AngularDeflection < Precision::Angular (),
                                       "GCPnts_TangentialDeflection::Initialize - Zero Deflection")

 parameters.Clear ();
 points    .Clear ();
 if (FirstParameter < LastParameter) {
   firstu = FirstParameter;
   lastu  = LastParameter;
 }
 else {
   lastu  = FirstParameter;
   firstu = LastParameter;
 }
 uTol                = UTol;
 angularDeflection   = AngularDeflection;
 curvatureDeflection = CurvatureDeflection;
 minNbPnts           = Max (MinimumOfPoints, 2);
 myMinLen            = Max(theMinLen, Precision::Confusion());

 switch (C.GetType()) {

 case GeomAbs_Line:   
   PerformLinear (C);
   break;

 case GeomAbs_Circle: 
   PerformCircular (C);
   break;

 case GeomAbs_BSplineCurve:
   {
     Handle_TheBSplineCurve BS = C.BSpline() ;
     if (BS->NbPoles() == 2 ) PerformLinear (C);
     else                     PerformCurve  (C);
     break;
   }
 case GeomAbs_BezierCurve:
   {
     Handle_TheBezierCurve  BZ = C.Bezier(); 
     if (BZ->NbPoles() == 2) PerformLinear (C);
     else                    PerformCurve  (C);
     break;
   }
 default: PerformCurve (C);
   
 }
}


//=======================================================================
//function : PerformLinear
//purpose  : 
//=======================================================================

void GCPnts_TangentialDeflection::PerformLinear (const TheCurve& C) {

  gp_Pnt P;
  D0 (C, firstu, P);
  parameters.Append (firstu);
  points    .Append (P);
  if (minNbPnts > 2)
  {
    Standard_Real Du = (lastu - firstu) / minNbPnts;
    Standard_Real U = firstu + Du;
    for (Standard_Integer i = 2; i < minNbPnts; i++)
    {
      D0 (C, U, P);
      parameters.Append (U);
      points    .Append (P);
      U += Du;
    }
  }
  D0 (C, lastu, P);
  parameters.Append (lastu);
  points    .Append (P);
}

//=======================================================================
//function : PerformCircular
//purpose  : 
//=======================================================================

void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C) 
{
  // akm 8/01/02 : check the radius before divide by it
  Standard_Real dfR = C.Circle().Radius();
  Standard_Real Du = GCPnts_TangentialDeflection::ArcAngularStep(
    dfR, curvatureDeflection, angularDeflection, myMinLen);
    
  const Standard_Real aDiff = lastu - firstu;
  // Round up number of points to satisfy curvatureDeflection more precisely
  Standard_Integer NbPoints = (Standard_Integer)Min(Ceiling(aDiff / Du), 1.0e+6);
  NbPoints = Max(NbPoints, minNbPnts - 1);
  Du       = aDiff / NbPoints;

  gp_Pnt P;
  Standard_Real U = firstu;
  for (Standard_Integer i = 1; i <= NbPoints; i++)
  {
    D0 (C, U, P);
    parameters.Append (U);
    points    .Append (P);
    U += Du;
  }
  D0 (C, lastu, P);
  parameters.Append (lastu);
  points    .Append (P);
}


//=======================================================================
//function : PerformCurve
//purpose  : On respecte ll'angle et la fleche, on peut imposer un nombre
//           minimum de points sur un element lineaire
//=======================================================================
void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)

{
  Standard_Integer i, j;       
  gp_XYZ V1, V2;
  gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
  Standard_Real Du, Dusave, MiddleU, L1, L2;

  Standard_Real U1       = firstu;   
  Standard_Real LTol     = Precision::Confusion ();  //protection longueur nulle
  Standard_Real ATol     = 1.e-2 * angularDeflection;
  if(ATol > 1.e-2)
    ATol = 1.e-2;
  else if(ATol < 1.e-7)
    ATol = 1.e-7;

  D0 (C, lastu, LastPoint);

  //Initialization du calcul

  Standard_Boolean NotDone = Standard_True;
  Dusave = (lastu - firstu)*Us3;
  Du     = Dusave;
  EvaluateDu (C, U1, CurrentPoint, Du, NotDone);
  parameters.Append (U1);
  points    .Append (CurrentPoint);

  // Used to detect "isLine" current bspline and in Du computation in general handling.
  Standard_Integer NbInterv = C.NbIntervals(GeomAbs_CN);
  TColStd_Array1OfReal Intervs(1, NbInterv+1);
  C.Intervals(Intervs, GeomAbs_CN);

  if (NotDone || Du > 5. * Dusave) {
    //C'est soit une droite, soit une singularite :
    V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
    L1 = V1.Modulus ();
    if (L1 > LTol)
    {
      //Si c'est une droite on verifie en calculant minNbPoints :
      Standard_Boolean IsLine   = Standard_True;
      Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
      switch (C.GetType()) {
      case GeomAbs_BSplineCurve:
        {
          Handle_TheBSplineCurve BS = C.BSpline() ;
          NbPoints = Max(BS->Degree() + 1, NbPoints);
          break;
        }
      case GeomAbs_BezierCurve:
        {
          Handle_TheBezierCurve  BZ = C.Bezier();
          NbPoints = Max(BZ->Degree() + 1, NbPoints);
          break;
        }
      default:
      ;}
      ////
      Standard_Real param = 0.;
      for (i = 1; i <= NbInterv && IsLine; ++i)
      {
        // Avoid usage intervals out of [firstu, lastu].
        if ((Intervs(i+1) < firstu) ||
            (Intervs(i)   > lastu))
        {
          continue;
        }
        // Fix border points in applicable intervals, to avoid be out of target interval.
        if ((Intervs(i)   < firstu) &&
            (Intervs(i+1) > firstu))
        {
          Intervs(i) = firstu;
        }
        if ((Intervs(i)   < lastu) &&
            (Intervs(i+1) > lastu))
        {
          Intervs(i + 1) = lastu;
        }

        Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
        for (j = 1; j <= NbPoints && IsLine; ++j)
        {
          param = Intervs(i) + j*delta;
          D0 (C, param, MiddlePoint);
          V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
          L2 = V2.Modulus ();
          if (L2 > LTol)
          {
            const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
            IsLine = (aAngle < ATol);
          }
        }
      }

      if (IsLine)
      {
        parameters.Clear();
        points    .Clear();

        PerformLinear(C);
        return;
      }
      else
      {
        //c'etait une singularite on continue :
        //Du = Dusave;
        EvaluateDu (C, param, MiddlePoint, Du, NotDone);
      }
    }
    else
    {
      Du = (lastu-firstu)/2.1;
      MiddleU = firstu + Du;
      D0 (C, MiddleU, MiddlePoint);
      V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
      L1 = V1.Modulus ();
      if (L1 < LTol)
      {
        // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
        // on renvoi un segment de 2 points   (protection)
        parameters.Append (lastu);
        points    .Append (LastPoint);
        return;
      }
    }
  }

  if (Du > Dusave) Du = Dusave;
  else             Dusave = Du;

  if (Du < uTol) {
    Du = lastu - firstu;
    if (Du < uTol) {
      parameters.Append (lastu);
      points    .Append (LastPoint);
      return;
    }
  }

  //Traitement normal pour une courbe
  Standard_Boolean MorePoints = Standard_True;
  Standard_Real U2            = firstu;   
  Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
  Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
  Standard_Boolean isNeedToCheck = Standard_False;
  gp_Pnt aPrevPoint = points.Last();

  while (MorePoints) {
    aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
    U2 += Du;

    if (U2 >= lastu) {                       //Bout de courbe
      U2 = lastu;
      CurrentPoint = LastPoint;
      Du = U2-U1;
      Dusave = Du;
    }
    else D0 (C, U2, CurrentPoint);           //Point suivant

    Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
    Standard_Boolean Correction, TooLarge, TooSmall;
    TooLarge   = Standard_False;
    Correction = Standard_True;
    TooSmall = Standard_False;

    while (Correction) {                     //Ajustement Du
      if (isNeedToCheck)
      {
        aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
        if (aIdx[1] > aIdx[0]) // Jump to another polynom.
        {
          if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
          {
            Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
            U2 = U1 + Du;
            if (U2 > lastu)
              U2 = lastu;
            D0 (C, U2, CurrentPoint);
          }
        }
      }
      MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
      D0 (C, MiddleU, MiddlePoint);

      V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
      V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
      L1 = V1.Modulus ();

      FCoef = (L1 > myMinLen) ? 
        V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;

      V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
      L1 = V1.Modulus ();
      L2 = V2.Modulus ();
      if (L1 > myMinLen && L2 > myMinLen)
      {
        Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
        ACoef = angg / AngleMax;
      }
      else
        ACoef = 0.0;

      //On retient le plus penalisant
      Coef = Max(ACoef, FCoef);

      if (isNeedToCheck && Coef < 0.55)
      {
        isNeedToCheck = Standard_False;
        Du = Dusave;
        U2 = U1 + Du;
        if (U2 > lastu)
          U2 = lastu;
        D0 (C, U2, CurrentPoint);
        continue;
      }

      if (Coef <= 1.0) {
        if (Abs (lastu-U2) < uTol) {
          parameters.Append (lastu);
          points    .Append (LastPoint);
          MorePoints = Standard_False;
          Correction = Standard_False;
        }
        else {
          if (Coef >= 0.55 || TooLarge) { 
            parameters.Append (U2);
            points    .Append (CurrentPoint);
            aPrevPoint = CurrentPoint;
            Correction = Standard_False;
            isNeedToCheck = Standard_True;
          }
          else if (TooSmall) {
            Correction = Standard_False;
            aPrevPoint = CurrentPoint;
          }
          else {
            TooSmall = Standard_True;
            //Standard_Real UUU2 = U2;
            Du += Min((U2-U1)*(1.-Coef), Du*Us3);

            U2 = U1 + Du;
            if (U2 > lastu)
              U2 = lastu;
            D0 (C, U2, CurrentPoint);
          }
        }
      }
      else {

        if (Coef >= 1.5) {
          if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion()))
          {
            parameters.Append (U1);
            points    .Append (aPrevPoint);
          }
          U2 = MiddleU;
          Du  = U2-U1;
          CurrentPoint = MiddlePoint;
        }
        else {
          Du*=0.9;
          U2 = U1 + Du;
          D0 (C, U2, CurrentPoint);
          TooLarge = Standard_True;
        }

      }
    }
    
    Du  = U2-U1;

    if (MorePoints) {
      if (U1 > firstu) {
        if (FCoef > ACoef) {
          //La fleche est critere de decoupage
          EvaluateDu (C, U2, CurrentPoint, Du, NotDone);
          if (NotDone) {
            Du += (Du-Dusave)*(Du/Dusave);
            if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
            if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
          }
        }
        else {
          //L'angle est le critere de decoupage
          Du += (Du-Dusave)*(Du/Dusave);
          if (Du > 1.5 * Dusave) Du = 1.5  * Dusave;
          if (Du < 0.75* Dusave) Du = 0.75 * Dusave;
        }
      }
      
      if (Du < uTol) {
        Du = lastu - U2;
        if (Du < uTol) {
          parameters.Append (lastu);
          points    .Append (LastPoint);
          MorePoints = Standard_False;
        }
        else if (Du*Us3 > uTol) Du*=Us3;
      }
      U1 = U2;
      Dusave = Du;
    }
  }
  //Recalage avant dernier point :
  i = points.Length()-1;
//  Real d = points (i).Distance (points (i+1));
// if (Abs(parameters (i) - parameters (i+1))<= 0.000001 || d < Precision::Confusion()) {
//    cout<<"deux points confondus"<<endl;
//    parameters.Remove (i+1);
//    points.Remove (i+1);
//    i--;
//  }
  if (i >= 2) {
    MiddleU = parameters (i-1);
    MiddleU = (lastu + MiddleU)*0.5;
    D0 (C, MiddleU, MiddlePoint);
    parameters.SetValue (i, MiddleU);
    points    .SetValue (i, MiddlePoint);
  }

  //-- On rajoute des points aux milieux des segments si le nombre
  //-- mini de points n'est pas atteint
  //--
  Standard_Integer Nbp = points.Length();

  //std::cout << "GCPnts_TangentialDeflection: Number of Points (" << Nbp << " " << minNbPnts << " )" << std::endl;

  while(Nbp < minNbPnts)
  { 
    for (i = 2; i <= Nbp; i += 2)
    {
      MiddleU = (parameters.Value(i-1)+parameters.Value(i))*0.5;
      D0 (C, MiddleU, MiddlePoint); 
      parameters.InsertBefore(i,MiddleU);
      points.InsertBefore(i,MiddlePoint);
      Nbp++;
    }
  }
  //Additional check for intervals
  Standard_Real MinLen2 = myMinLen * myMinLen;
  Standard_Integer MaxNbp = 10 * Nbp;
  for(i = 1; i < Nbp; ++i)
  {
    U1 = parameters(i);
    U2 = parameters(i + 1);

    if (U2 - U1 <= uTol)
    {
      continue;
    }

    // Check maximal deflection on interval;
    Standard_Real dmax = 0.;
    Standard_Real umax = 0.;
    Standard_Real amax = 0.;
    EstimDefl(C, U1, U2, dmax, umax);
    const gp_Pnt& P1 = points(i);
    const gp_Pnt& P2 = points(i+1);
    D0(C, umax, MiddlePoint);
    amax = EstimAngl(P1, MiddlePoint, P2);
    if(dmax > curvatureDeflection || amax > AngleMax)
    {
      if(umax - U1 > uTol && U2 - umax > uTol)
      {
        if (P1.SquareDistance(MiddlePoint) > MinLen2 && P2.SquareDistance(MiddlePoint) > MinLen2)
        {
          parameters.InsertAfter(i, umax);
          points.InsertAfter(i, MiddlePoint);
          ++Nbp;
          --i; //To compensate ++i in loop header: i must point to first part of splitted interval
          if(Nbp > MaxNbp)
          {
            break;
          }
        }
      }
    }
  }
  // 
}

//=======================================================================
//function : EstimDefl
//purpose  : Estimation of maximal deflection for interval [U1, U2]
//           
//=======================================================================
void GCPnts_TangentialDeflection::EstimDefl (const TheCurve& C,
                            const Standard_Real U1, const Standard_Real U2, 
                            Standard_Real& MaxDefl, Standard_Real& UMax)
{
  Standard_Real Du = (lastu - firstu);
  //
  TheMaxCurvLinDist aFunc(C, U1, U2);
  //
  const Standard_Integer aNbIter = 100;
  Standard_Real reltol = Max(1.e-3, 2.*uTol/((Abs(U1) + Abs(U2))));
  //
  math_BrentMinimum anOptLoc(reltol, aNbIter, uTol);
  anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2);
  if(anOptLoc.IsDone())
  {
    MaxDefl = Sqrt(-anOptLoc.Minimum());
    UMax = anOptLoc.Location();
    return;
  }
  //
  math_Vector aLowBorder(1,1);
  math_Vector aUppBorder(1,1);
  math_Vector aSteps(1,1);
  //
  aSteps(1) = Max(0.1 * Du, 100. * uTol);
  Standard_Integer aNbParticles = Max(8, RealToInt(32 * (U2 - U1) / Du));
  //
  aLowBorder(1) = U1;
  aUppBorder(1) = U2;
  //
  //
  Standard_Real aValue;
  math_Vector aT(1,1);
  TheMaxCurvLinDistMV aFuncMV(aFunc);

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