[occt.git] / src / Extrema / Extrema_GExtPC.gxx

// Created on: 1992-10-19
// Created by: Laurent PAINNOT
// Copyright (c) 1992-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.

// Modified by cma, Fri Dec 10 18:11:56 1993

#include Extrema_EPC_hxx
#include ThePOnC_hxx
#include TheSequenceOfPOnC_hxx
#include ThePoint_hxx
#include TheVector_hxx
#include TheExtPElC_hxx
#include <StdFail_NotDone.hxx>
#include <Standard_Failure.hxx>
#include <GeomAbs_CurveType.hxx>
#include <Precision.hxx>
#include <ElCLib.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_HArray1OfReal.hxx>
#include <NCollection_Array1.hxx>


//=======================================================================
//function : Perform
//purpose  : 
//=======================================================================
void Extrema_GExtPC::Perform(const ThePoint& P)
{
  mySqDist.Clear();
  mypoint.Clear();
  myismin.Clear();
  Standard_Integer i, NbExt, n;                     
  Standard_Real U;
  mysample = 17;
  Standard_Real t3d = Precision::Confusion();
  if (Precision::IsInfinite(myuinf)) mydist1 = RealLast();
  else {
    Pf = TheCurveTool::Value(*((TheCurve*)myC), myuinf);
    mydist1 = P.SquareDistance(Pf);
  }
  
  if (Precision::IsInfinite(myusup)) mydist2 = RealLast();
  else {
    Pl = TheCurveTool::Value(*((TheCurve*)myC), myusup);
    mydist2 = P.SquareDistance(Pl);
  }

  TheCurve & aCurve = *((TheCurve*)myC);

  switch(type) {
  case GeomAbs_Circle: 
    {
      myExtPElC.Perform(P, TheCurveTool::Circle(aCurve), t3d, myuinf, myusup);
    }
    break;
  case GeomAbs_Ellipse: 
    {
      myExtPElC.Perform(P, TheCurveTool::Ellipse(aCurve), t3d, myuinf, myusup);
    }
    break;
  case GeomAbs_Parabola: 
    {
      myExtPElC.Perform(P, TheCurveTool::Parabola(aCurve), t3d,myuinf,myusup);
    }
    break;
  case GeomAbs_Hyperbola: 
    {
      myExtPElC.Perform(P,TheCurveTool::Hyperbola(aCurve),t3d, myuinf, myusup);
    }
    break;
  case GeomAbs_Line: 
    {
      myExtPElC.Perform(P, TheCurveTool::Line(aCurve), t3d, myuinf, myusup);
    }
    break;
  case GeomAbs_BezierCurve:
    {
      myintuinf = myuinf;
      myintusup = myusup;
      mysample = (TheCurveTool::Bezier(aCurve))->NbPoles() * 2;
      myExtPC.Initialize(aCurve);
      IntervalPerform(P);
      return;
    }
  case GeomAbs_BSplineCurve: 
    {
      const Standard_Integer 
        aFirstIdx = TheCurveTool::BSpline(aCurve)->FirstUKnotIndex(),
        aLastIdx  = TheCurveTool::BSpline(aCurve)->LastUKnotIndex();
      // const reference can not be used due to implementation of BRep_Adaptor.
      TColStd_Array1OfReal aKnots(aFirstIdx, aLastIdx);
      TheCurveTool::BSpline(aCurve)->Knots(aKnots); 

      // Workaround to work with:
      // blend, where knots may be moved from parameter space.
      Standard_Real aPeriodJump = 0.0;
      // Avoid problem with too close knots.
      const Standard_Real aTolCoeff = (myusup - myuinf) * Precision::PConfusion();
      if (TheCurveTool::IsPeriodic(aCurve))
      {
        Standard_Integer aPeriodShift =
          Standard_Integer ((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
        
        if (myuinf < aKnots(aFirstIdx) - aTolCoeff)
          aPeriodShift--;

        aPeriodJump = TheCurveTool::Period(aCurve) * aPeriodShift;
      }

      Standard_Integer anIdx;

      // Find first and last used knot.
      Standard_Integer aFirstUsedKnot = aFirstIdx,
                        aLastUsedKnot = aLastIdx;
      for(anIdx = aFirstIdx; anIdx <= aLastIdx; anIdx++)
      {
        Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
        if (myuinf >= aKnot - aTolCoeff)
          aFirstUsedKnot = anIdx;
        else
          break;

      }
      for(anIdx = aLastIdx; anIdx >= aFirstIdx; anIdx--)
      {
        Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
        if (myusup <= aKnot + aTolCoeff)
          aLastUsedKnot = anIdx;
        else
          break;
      }

      if (aFirstUsedKnot == aLastUsedKnot)
      {
        // Degenerated case:
        // Some bounds lies out of curve param space.
        // In this case build one interval with [myuinf, myusup].
        // Parameters of these indexes will be redefined.
        aFirstUsedKnot = aFirstIdx;
        aLastUsedKnot  = aFirstIdx + 1;
      }

      mysample = (TheCurveTool::BSpline(aCurve))->Degree() + 1;

      // Fill sample points.
      Standard_Integer aValIdx = 1;
      NCollection_Array1<Standard_Real> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
      NCollection_Array1<Standard_Real> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
      for(anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
      {
        Standard_Real aF = aKnots(anIdx) + aPeriodJump,
                      aL = aKnots(anIdx + 1) + aPeriodJump;

        if (anIdx == aFirstUsedKnot)
          aF = myuinf;
        if (anIdx == aLastUsedKnot - 1)
          aL = myusup;

        Standard_Real aStep = (aL - aF) / mysample;
        for(Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
        {
          Standard_Real aCurrentParam = aF + aStep * aPntIdx;
          aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
          aParam(aValIdx) = aCurrentParam;
          aValIdx++;
        }
      }
      // Fill last point.
      aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
      aParam(aValIdx) = myusup;

       myExtPC.Initialize(aCurve);

      // Find extremas.
      for(anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
      {
        if (aVal(anIdx) <= Precision::SquareConfusion())
        {
          mySqDist.Append(aVal(anIdx));
          myismin.Append(Standard_True);
          mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
        }
        if ((aVal(anIdx) >= aVal(anIdx + 1) &&
             aVal(anIdx) >= aVal(anIdx - 1)) ||
            (aVal(anIdx) <= aVal(anIdx + 1) &&
             aVal(anIdx) <= aVal(anIdx - 1)) )
        {
          myintuinf = aParam(anIdx - 1);
          myintusup = aParam(anIdx + 1);

          IntervalPerform(P);
        }
      }

      // Solve on the first and last intervals.
      if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
      {
        ThePoint aP1, aP2;
        TheVector aV1, aV2;
        TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()),     aP1, aV1);
        TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
        TheVector aBase1(P, aP1), aBase2(P, aP2);
        Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
        Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
        if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
        {
          // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
          // Necessary condition - when point lies on curve.
          // Necessary condition - when derivative of point is too small.
          if(aVal1 * aVal2 <= 0.0 ||
            aBase1.Dot(aBase2) <= 0.0 ||
            2.0 * Abs(aVal1) < Precision::Confusion() )
          {
            myintuinf = aParam(aVal.Lower());
            myintusup = aParam(aVal.Lower() + 1);
            IntervalPerform(P);
          }
        }
      }

      if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
      {
        ThePoint aP1, aP2;
        TheVector aV1, aV2;
        TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
        TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()),     aP2, aV2);
        TheVector aBase1(P, aP1), aBase2(P, aP2);
        Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
        Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2

        if(!(Precision::IsInfinite(aVal1) || Precision::IsInfinite(aVal2)))
        {
          // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
          // Necessary condition - when point lies on curve.
          // Necessary condition - when derivative of point is too small.
          if(aVal1 * aVal2 <= 0.0 ||
            aBase1.Dot(aBase2) <= 0.0 ||
            2.0 * Abs(aVal2) < Precision::Confusion() )
          {
            myintuinf = aParam(aVal.Upper() - 1);
            myintusup = aParam(aVal.Upper());
            IntervalPerform(P);
          }
        }
      }

      mydone = Standard_True;
      break;
    }
  default:
    {
      const Standard_Integer aMaxSample = 17;
      Standard_Boolean IntExtIsDone = Standard_False;
      Standard_Boolean IntIsNotValid;
      Handle(TColStd_HArray1OfReal) theHInter;
      n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
      if (n > 1)
      {
        theHInter = new TColStd_HArray1OfReal(1, n + 1);
        TheCurveTool::Intervals(aCurve, theHInter->ChangeArray1(), GeomAbs_C2);
      }
      else
      {
        theHInter = TheCurveTool::DeflCurvIntervals(aCurve);
        n = theHInter->Length() - 1;
      }
      mysample = Max(mysample / n, aMaxSample);
      Standard_Real maxint = 0.;
      for (i = 1; i <= n; ++i)
      {
        Standard_Real dt = theHInter->Value(i + 1) - theHInter->Value(i);
        if (maxint < dt)
        {
          maxint = dt;
        }
      }
      Standard_Boolean isPeriodic = TheCurveTool::IsPeriodic(aCurve);
      TheVector V1;
      ThePoint PP;
      Standard_Real s1 = 0.0 ;
      Standard_Real s2 = 0.0;
      myExtPC.Initialize(aCurve);
      for (i = 1; i <= n; i++)
      {
        myintuinf = theHInter->Value(i);
        myintusup = theHInter->Value(i+1);
        mysample = Max(RealToInt(aMaxSample*(myintusup - myintuinf) / maxint), 3);
        
        Standard_Real anInfToCheck = myintuinf;
        Standard_Real aSupToCheck = myintusup;
        
        if (isPeriodic) {
          Standard_Real aPeriod = TheCurveTool::Period(aCurve);
          anInfToCheck = ElCLib::InPeriod(myintuinf, myuinf, myuinf+aPeriod);
          aSupToCheck = myintusup+(anInfToCheck-myintuinf);
        }    
        IntIsNotValid = (myuinf > aSupToCheck) || (myusup < anInfToCheck);

        if(IntIsNotValid) continue;

        if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
        if (myusup <= aSupToCheck) aSupToCheck = myusup;
        if((aSupToCheck - anInfToCheck) <= mytolu) continue;
        
        if (i != 1)
          {
          TheCurveTool::D1(aCurve, myintuinf, PP, V1);
          s1 = (TheVector(P, PP))*V1;
          if (s1*s2 < 0.0) {
            mySqDist.Append(PP.SquareDistance(P));
            myismin.Append((s1 < 0.0));
            mypoint.Append(ThePOnC(myintuinf, PP));
          }
        }
        if (i != n) {
          TheCurveTool::D1(aCurve, myintusup, PP, V1);
          s2 = (TheVector(P, PP))*V1;
        }

        IntervalPerform(P);
        IntExtIsDone = IntExtIsDone || mydone;
      }

      mydone = IntExtIsDone;
      break;
    }
  }

  // Postprocessing.
  if (type == GeomAbs_BSplineCurve ||
      type == GeomAbs_OffsetCurve ||
      type == GeomAbs_OtherCurve)
  {
    // Additional checking if the point is on the first or last point of the curve
    // and does not added yet.
    if (mydist1 < Precision::SquareConfusion() || 
        mydist2 < Precision::SquareConfusion())
    {
      Standard_Boolean isFirstAdded = Standard_False;
      Standard_Boolean isLastAdded  = Standard_False;
      Standard_Integer aNbPoints = mypoint.Length();
      for (i = 1; i <= aNbPoints; i++)
      {
        U = mypoint.Value(i).Parameter();
        if (Abs(U - myuinf) < mytolu)
          isFirstAdded = Standard_True;
        else if (Abs(myusup - U) < mytolu)
          isLastAdded = Standard_True;
      }
      if (!isFirstAdded && mydist1 < Precision::SquareConfusion())
      {
        mySqDist.Prepend(mydist1);
        myismin.Prepend(Standard_True);
        mypoint.Prepend(ThePOnC(myuinf, Pf));
      }
      if (!isLastAdded && mydist2 < Precision::SquareConfusion())
      {
        mySqDist.Append(mydist2);
        myismin.Append(Standard_True);
        mypoint.Append(ThePOnC(myusup, Pl));
      }
      mydone = Standard_True;
    }
  }
  else
  {
    // In analytical case
    mydone = myExtPElC.IsDone();
    if (mydone)
    {
      NbExt = myExtPElC.NbExt();
      for (i = 1; i <= NbExt; i++)
      {
        // Verification de la validite des parametres:
        ThePOnC PC = myExtPElC.Point(i);
        U = PC.Parameter();
        if (TheCurveTool::IsPeriodic(aCurve))
        {
          U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(aCurve));
        }
        if ((U >= myuinf-mytolu) && (U <= myusup+mytolu))
        {
          PC.SetValues(U, myExtPElC.Point(i).Value());
          mySqDist.Append(myExtPElC.SquareDistance(i));
          myismin.Append(myExtPElC.IsMin(i));
          mypoint.Append(PC);
        }
      }
    }
  }
}


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

void Extrema_GExtPC::Initialize(const TheCurve&     C,
				const Standard_Real Uinf,
				const Standard_Real Usup,
				const Standard_Real TolF) 
{
  myC = (Standard_Address)&C;
  myintuinf = myuinf = Uinf;
  myintusup = myusup = Usup;
  mytolf = TolF;
  mytolu = TheCurveTool::Resolution(*((TheCurve*)myC), Precision::Confusion());
  type = TheCurveTool::GetType(C);
  mydone = Standard_False;
  mydist1 = RealLast();
  mydist2 = RealLast();
  mysample = 17;
}


//=======================================================================
//function : IntervalPerform
//purpose  : 
//=======================================================================

void Extrema_GExtPC::IntervalPerform(const ThePoint& P)
{
  Standard_Integer i;
  Standard_Real U;
  myExtPC.Initialize(mysample, myintuinf, myintusup, mytolu, mytolf);
  myExtPC.Perform(P);
  mydone = myExtPC.IsDone();
  if (mydone)
  {
    Standard_Integer NbExt = myExtPC.NbExt();
    for (i = 1; i <= NbExt; i++)
    {
      // Verification de la validite des parametres pour le cas trimme:
      ThePOnC PC = myExtPC.Point(i);
      U = PC.Parameter();
      if (TheCurveTool::IsPeriodic(*((TheCurve*)myC)))
      {
        U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(*((TheCurve*)myC)));
      }
      if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
      {
        PC.SetValues(U, PC.Value());
        mySqDist.Append(myExtPC.SquareDistance(i));
        myismin.Append(myExtPC.IsMin(i));
        mypoint.Append(PC);
      }
    }
  }
}


//=======================================================================
//function : Extrema_GExtPC
//purpose  : 
//=======================================================================

Extrema_GExtPC::Extrema_GExtPC()
{
  myC = 0;
  mydone = Standard_False;
  mydist1 = RealLast();
  mydist2 = RealLast();
  mytolu = 0.0;
  mytolf = 0.0;
  mysample = 17;
  myintuinf = myintusup = myuinf = myusup = Precision::Infinite();
  type = GeomAbs_OtherCurve;
}

//=======================================================================
//function : Extrema_GExtPC
//purpose  : 
//=======================================================================

Extrema_GExtPC::Extrema_GExtPC(const ThePoint&           P, 
			       const TheCurve&           C,
			       const Standard_Real       Uinf,
			       const Standard_Real       Usup,
			       const Standard_Real       TolF) 
{
  Initialize(C, Uinf, Usup, TolF);
  Perform(P);
}

//=======================================================================
//function : Extrema_GExtPC
//purpose  : 
//=======================================================================

Extrema_GExtPC::Extrema_GExtPC(const ThePoint&     P, 
			       const TheCurve&     C,
			       const Standard_Real TolF) 
{
  Initialize(C, TheCurveTool::FirstParameter(C), 
	     TheCurveTool::LastParameter(C), TolF);
  Perform(P);
}


//=======================================================================
//function : IsDone
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_GExtPC::IsDone() const
{
  return mydone;
}


//=======================================================================
//function : Value
//purpose  : 
//=======================================================================

Standard_Real Extrema_GExtPC::SquareDistance(const Standard_Integer N) const 
{
  if ((N < 1) || (N > NbExt())) throw Standard_OutOfRange();
  return mySqDist.Value(N);
}


//=======================================================================
//function : NbExt
//purpose  : 
//=======================================================================

Standard_Integer Extrema_GExtPC::NbExt() const
{
  if (!IsDone()) throw StdFail_NotDone();
  return mySqDist.Length();
}


//=======================================================================
//function : IsMin
//purpose  : 
//=======================================================================

Standard_Boolean Extrema_GExtPC::IsMin(const Standard_Integer N) const
{
  if ((N < 1) || (N > NbExt())) throw Standard_OutOfRange();
  return myismin.Value(N);
}


//=======================================================================
//function : Point
//purpose  : 
//=======================================================================

const ThePOnC & Extrema_GExtPC::Point(const Standard_Integer N) const
{
  if ((N < 1) || (N > NbExt())) throw Standard_OutOfRange();
  return mypoint.Value(N);
}


//=======================================================================
//function : TrimmedDistances
//purpose  : 
//=======================================================================

void Extrema_GExtPC::TrimmedSquareDistances(Standard_Real& dist1, 
				      Standard_Real& dist2,
				      ThePoint& P1,
				      ThePoint& P2) const
{
  dist1 = mydist1;
  dist2 = mydist2;
  P1 = Pf;
  P2 = Pl;
}
Commit	Line	Data
b311480e	1	// Created on: 1992-10-19
	2	// Created by: Laurent PAINNOT
	3	// Copyright (c) 1992-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.
b311480e	16
7fd59977	17	// Modified by cma, Fri Dec 10 18:11:56 1993
7fd59977	18
	19	#include Extrema_EPC_hxx
	20	#include ThePOnC_hxx
	21	#include TheSequenceOfPOnC_hxx
	22	#include ThePoint_hxx
	23	#include TheVector_hxx
	24	#include TheExtPElC_hxx
	25	#include <StdFail_NotDone.hxx>
	26	#include <Standard_Failure.hxx>
	27	#include <GeomAbs_CurveType.hxx>
	28	#include <Precision.hxx>
	29	#include <ElCLib.hxx>
	30	#include <TColStd_Array1OfReal.hxx>
f6b08ecf	31	#include <TColStd_HArray1OfReal.hxx>
c8bf1eb7	32	#include <NCollection_Array1.hxx>
7fd59977	33
	34
	35	//=======================================================================
	36	//function : Perform
	37	//purpose :
	38	//=======================================================================
7fd59977	39	void Extrema_GExtPC::Perform(const ThePoint& P)
	40	{
	41	mySqDist.Clear();
	42	mypoint.Clear();
	43	myismin.Clear();
	44	Standard_Integer i, NbExt, n;
	45	Standard_Real U;
	46	mysample = 17;
	47	Standard_Real t3d = Precision::Confusion();
	48	if (Precision::IsInfinite(myuinf)) mydist1 = RealLast();
	49	else {
	50	Pf = TheCurveTool::Value(((TheCurve)myC), myuinf);
	51	mydist1 = P.SquareDistance(Pf);
	52	}
	53
	54	if (Precision::IsInfinite(myusup)) mydist2 = RealLast();
	55	else {
	56	Pl = TheCurveTool::Value(((TheCurve)myC), myusup);
	57	mydist2 = P.SquareDistance(Pl);
	58	}
	59
c8bf1eb7	60	TheCurve & aCurve = ((TheCurve)myC);
c8bf1eb7	61
7fd59977	62	switch(type) {
	63	case GeomAbs_Circle:
	64	{
c8bf1eb7	65	myExtPElC.Perform(P, TheCurveTool::Circle(aCurve), t3d, myuinf, myusup);
7fd59977	66	}
	67	break;
	68	case GeomAbs_Ellipse:
	69	{
c8bf1eb7	70	myExtPElC.Perform(P, TheCurveTool::Ellipse(aCurve), t3d, myuinf, myusup);
7fd59977	71	}
	72	break;
	73	case GeomAbs_Parabola:
	74	{
c8bf1eb7	75	myExtPElC.Perform(P, TheCurveTool::Parabola(aCurve), t3d,myuinf,myusup);
7fd59977	76	}
	77	break;
	78	case GeomAbs_Hyperbola:
	79	{
c8bf1eb7	80	myExtPElC.Perform(P,TheCurveTool::Hyperbola(aCurve),t3d, myuinf, myusup);
7fd59977	81	}
	82	break;
	83	case GeomAbs_Line:
	84	{
c8bf1eb7	85	myExtPElC.Perform(P, TheCurveTool::Line(aCurve), t3d, myuinf, myusup);
7fd59977	86	}
	87	break;
	88	case GeomAbs_BezierCurve:
	89	{
	90	myintuinf = myuinf;
	91	myintusup = myusup;
c8bf1eb7	92	mysample = (TheCurveTool::Bezier(aCurve))->NbPoles() * 2;
c8bf1eb7	93	myExtPC.Initialize(aCurve);
7fd59977	94	IntervalPerform(P);
	95	return;
	96	}
	97	case GeomAbs_BSplineCurve:
	98	{
c8bf1eb7	99	const Standard_Integer
	100	aFirstIdx = TheCurveTool::BSpline(aCurve)->FirstUKnotIndex(),
	101	aLastIdx = TheCurveTool::BSpline(aCurve)->LastUKnotIndex();
	102	// const reference can not be used due to implementation of BRep_Adaptor.
	103	TColStd_Array1OfReal aKnots(aFirstIdx, aLastIdx);
	104	TheCurveTool::BSpline(aCurve)->Knots(aKnots);
	105
	106	// Workaround to work with:
c5a10cf7	107	// blend, where knots may be moved from parameter space.
c8bf1eb7	108	Standard_Real aPeriodJump = 0.0;
c5a10cf7	109	// Avoid problem with too close knots.
0cbfb9f1	110	const Standard_Real aTolCoeff = (myusup - myuinf) * Precision::PConfusion();
c8bf1eb7	111	if (TheCurveTool::IsPeriodic(aCurve))
	112	{
	113	Standard_Integer aPeriodShift =
	114	Standard_Integer ((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
	115
0cbfb9f1	116	if (myuinf < aKnots(aFirstIdx) - aTolCoeff)
c8bf1eb7	117	aPeriodShift--;
	118
	119	aPeriodJump = TheCurveTool::Period(aCurve) * aPeriodShift;
	120	}
	121
	122	Standard_Integer anIdx;
	123
1581d651	124	// Find first and last used knot.
c8bf1eb7	125	Standard_Integer aFirstUsedKnot = aFirstIdx,
	126	aLastUsedKnot = aLastIdx;
	127	for(anIdx = aFirstIdx; anIdx <= aLastIdx; anIdx++)
	128	{
	129	Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
0cbfb9f1	130	if (myuinf >= aKnot - aTolCoeff)
c8bf1eb7	131	aFirstUsedKnot = anIdx;
	132	else
	133	break;
	134
	135	}
	136	for(anIdx = aLastIdx; anIdx >= aFirstIdx; anIdx--)
	137	{
	138	Standard_Real aKnot = aKnots(anIdx) + aPeriodJump;
0cbfb9f1	139	if (myusup <= aKnot + aTolCoeff)
c8bf1eb7	140	aLastUsedKnot = anIdx;
	141	else
	142	break;
	143	}
	144
1581d651	145	if (aFirstUsedKnot == aLastUsedKnot)
	146	{
	147	// Degenerated case:
	148	// Some bounds lies out of curve param space.
	149	// In this case build one interval with [myuinf, myusup].
	150	// Parameters of these indexes will be redefined.
	151	aFirstUsedKnot = aFirstIdx;
	152	aLastUsedKnot = aFirstIdx + 1;
	153	}
	154
c8bf1eb7	155	mysample = (TheCurveTool::BSpline(aCurve))->Degree() + 1;
	156
	157	// Fill sample points.
	158	Standard_Integer aValIdx = 1;
	159	NCollection_Array1<Standard_Real> aVal(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
	160	NCollection_Array1<Standard_Real> aParam(1, (mysample) * (aLastUsedKnot - aFirstUsedKnot) + 1);
	161	for(anIdx = aFirstUsedKnot; anIdx < aLastUsedKnot; anIdx++)
	162	{
	163	Standard_Real aF = aKnots(anIdx) + aPeriodJump,
	164	aL = aKnots(anIdx + 1) + aPeriodJump;
	165
	166	if (anIdx == aFirstUsedKnot)
	167	aF = myuinf;
	168	if (anIdx == aLastUsedKnot - 1)
	169	aL = myusup;
	170
	171	Standard_Real aStep = (aL - aF) / mysample;
	172	for(Standard_Integer aPntIdx = 0; aPntIdx < mysample; aPntIdx++)
	173	{
	174	Standard_Real aCurrentParam = aF + aStep * aPntIdx;
	175	aVal(aValIdx) = TheCurveTool::Value(aCurve, aCurrentParam).SquareDistance(P);
	176	aParam(aValIdx) = aCurrentParam;
	177	aValIdx++;
	178	}
	179	}
	180	// Fill last point.
	181	aVal(aValIdx) = TheCurveTool::Value(aCurve, myusup).SquareDistance(P);
	182	aParam(aValIdx) = myusup;
	183
	184	myExtPC.Initialize(aCurve);
	185
	186	// Find extremas.
	187	for(anIdx = aVal.Lower() + 1; anIdx < aVal.Upper(); anIdx++)
	188	{
	189	if (aVal(anIdx) <= Precision::SquareConfusion())
	190	{
	191	mySqDist.Append(aVal(anIdx));
	192	myismin.Append(Standard_True);
	193	mypoint.Append(ThePOnC(aParam(anIdx), TheCurveTool::Value(aCurve, aParam(anIdx))));
	194	}
	195	if ((aVal(anIdx) >= aVal(anIdx + 1) &&
	196	aVal(anIdx) >= aVal(anIdx - 1)) \|\|
	197	(aVal(anIdx) <= aVal(anIdx + 1) &&
	198	aVal(anIdx) <= aVal(anIdx - 1)) )
	199	{
	200	myintuinf = aParam(anIdx - 1);
	201	myintusup = aParam(anIdx + 1);
	202
	203	IntervalPerform(P);
	204	}
	205	}
	206
c5a10cf7	207	// Solve on the first and last intervals.
f4dee9bb	208	if (mydist1 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist1))
c8bf1eb7	209	{
	210	ThePoint aP1, aP2;
	211	TheVector aV1, aV2;
	212	TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower()), aP1, aV1);
	213	TheCurveTool::D1(aCurve, aParam.Value(aParam.Lower() + 1), aP2, aV2);
	214	TheVector aBase1(P, aP1), aBase2(P, aP2);
	215	Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
	216	Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
f4dee9bb	217	if(!(Precision::IsInfinite(aVal1) \|\| Precision::IsInfinite(aVal2)))
c8bf1eb7	218	{
f4dee9bb	219	// Derivatives have opposite signs - min or max inside of interval (sufficient condition).
	220	// Necessary condition - when point lies on curve.
	221	// Necessary condition - when derivative of point is too small.
	222	if(aVal1 * aVal2 <= 0.0 \|\|
	223	aBase1.Dot(aBase2) <= 0.0 \|\|
	224	2.0 * Abs(aVal1) < Precision::Confusion() )
	225	{
	226	myintuinf = aParam(aVal.Lower());
	227	myintusup = aParam(aVal.Lower() + 1);
	228	IntervalPerform(P);
	229	}
c8bf1eb7	230	}
	231	}
	232
f4dee9bb	233	if (mydist2 > Precision::SquareConfusion() && !Precision::IsPositiveInfinite(mydist2))
c8bf1eb7	234	{
	235	ThePoint aP1, aP2;
	236	TheVector aV1, aV2;
	237	TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper() - 1), aP1, aV1);
	238	TheCurveTool::D1(aCurve, aParam.Value(aParam.Upper()), aP2, aV2);
	239	TheVector aBase1(P, aP1), aBase2(P, aP2);
	240	Standard_Real aVal1 = aV1.Dot(aBase1); // Derivative of (C(u) - P)^2
	241	Standard_Real aVal2 = aV2.Dot(aBase2); // Derivative of (C(u) - P)^2
	242
f4dee9bb	243	if(!(Precision::IsInfinite(aVal1) \|\| Precision::IsInfinite(aVal2)))
c8bf1eb7	244	{
f4dee9bb	245	// Derivatives have opposite signs - min or max inside of interval (sufficient condition).
	246	// Necessary condition - when point lies on curve.
	247	// Necessary condition - when derivative of point is too small.
	248	if(aVal1 * aVal2 <= 0.0 \|\|
	249	aBase1.Dot(aBase2) <= 0.0 \|\|
	250	2.0 * Abs(aVal2) < Precision::Confusion() )
	251	{
	252	myintuinf = aParam(aVal.Upper() - 1);
	253	myintusup = aParam(aVal.Upper());
	254	IntervalPerform(P);
	255	}
c8bf1eb7	256	}
	257	}
	258
	259	mydone = Standard_True;
	260	break;
7fd59977	261	}
1aec3320	262	default:
7fd59977	263	{
f6b08ecf	264	const Standard_Integer aMaxSample = 17;
7fd59977	265	Standard_Boolean IntExtIsDone = Standard_False;
7fd59977	266	Standard_Boolean IntIsNotValid;
f6b08ecf	267	Handle(TColStd_HArray1OfReal) theHInter;
c8bf1eb7	268	n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
f6b08ecf	269	if (n > 1)
	270	{
	271	theHInter = new TColStd_HArray1OfReal(1, n + 1);
	272	TheCurveTool::Intervals(aCurve, theHInter->ChangeArray1(), GeomAbs_C2);
	273	}
	274	else
	275	{
	276	theHInter = TheCurveTool::DeflCurvIntervals(aCurve);
	277	n = theHInter->Length() - 1;
	278	}
	279	mysample = Max(mysample / n, aMaxSample);
	280	Standard_Real maxint = 0.;
	281	for (i = 1; i <= n; ++i)
	282	{
	283	Standard_Real dt = theHInter->Value(i + 1) - theHInter->Value(i);
	284	if (maxint < dt)
	285	{
	286	maxint = dt;
	287	}
	288	}
c8bf1eb7	289	Standard_Boolean isPeriodic = TheCurveTool::IsPeriodic(aCurve);
7fd59977	290	TheVector V1;
	291	ThePoint PP;
	292	Standard_Real s1 = 0.0 ;
c8bf1eb7	293	Standard_Real s2 = 0.0;
	294	myExtPC.Initialize(aCurve);
	295	for (i = 1; i <= n; i++)
	296	{
f6b08ecf	297	myintuinf = theHInter->Value(i);
	298	myintusup = theHInter->Value(i+1);
	299	mysample = Max(RealToInt(aMaxSample*(myintusup - myintuinf) / maxint), 3);
253881cf	300
	301	Standard_Real anInfToCheck = myintuinf;
	302	Standard_Real aSupToCheck = myintusup;
	303
	304	if (isPeriodic) {
c8bf1eb7	305	Standard_Real aPeriod = TheCurveTool::Period(aCurve);
253881cf	306	anInfToCheck = ElCLib::InPeriod(myintuinf, myuinf, myuinf+aPeriod);
	307	aSupToCheck = myintusup+(anInfToCheck-myintuinf);
	308	}
	309	IntIsNotValid = (myuinf > aSupToCheck) \|\| (myusup < anInfToCheck);
	310
	311	if(IntIsNotValid) continue;
	312
	313	if (myuinf >= anInfToCheck) anInfToCheck = myuinf;
	314	if (myusup <= aSupToCheck) aSupToCheck = myusup;
	315	if((aSupToCheck - anInfToCheck) <= mytolu) continue;
32ca7a51	316
	317	if (i != 1)
	318	{
c8bf1eb7	319	TheCurveTool::D1(aCurve, myintuinf, PP, V1);
253881cf	320	s1 = (TheVector(P, PP))*V1;
	321	if (s1*s2 < 0.0) {
	322	mySqDist.Append(PP.SquareDistance(P));
	323	myismin.Append((s1 < 0.0));
	324	mypoint.Append(ThePOnC(myintuinf, PP));
	325	}
	326	}
	327	if (i != n) {
c8bf1eb7	328	TheCurveTool::D1(aCurve, myintusup, PP, V1);
253881cf	329	s2 = (TheVector(P, PP))*V1;
253881cf	330	}
32ca7a51	331
253881cf	332	IntervalPerform(P);
253881cf	333	IntExtIsDone = IntExtIsDone \|\| mydone;
7fd59977	334	}
c8bf1eb7	335
7fd59977	336	mydone = IntExtIsDone;
c8bf1eb7	337	break;
	338	}
	339	}
94f71cad	340
c8bf1eb7	341	// Postprocessing.
c8bf1eb7	342	if (type == GeomAbs_BSplineCurve \|\|
1aec3320	343	type == GeomAbs_OffsetCurve \|\|
c8bf1eb7	344	type == GeomAbs_OtherCurve)
	345	{
	346	// Additional checking if the point is on the first or last point of the curve
	347	// and does not added yet.
	348	if (mydist1 < Precision::SquareConfusion() \|\|
	349	mydist2 < Precision::SquareConfusion())
	350	{
	351	Standard_Boolean isFirstAdded = Standard_False;
	352	Standard_Boolean isLastAdded = Standard_False;
	353	Standard_Integer aNbPoints = mypoint.Length();
	354	for (i = 1; i <= aNbPoints; i++)
94f71cad	355	{
c8bf1eb7	356	U = mypoint.Value(i).Parameter();
	357	if (Abs(U - myuinf) < mytolu)
	358	isFirstAdded = Standard_True;
	359	else if (Abs(myusup - U) < mytolu)
	360	isLastAdded = Standard_True;
	361	}
	362	if (!isFirstAdded && mydist1 < Precision::SquareConfusion())
	363	{
	364	mySqDist.Prepend(mydist1);
	365	myismin.Prepend(Standard_True);
	366	mypoint.Prepend(ThePOnC(myuinf, Pf));
	367	}
	368	if (!isLastAdded && mydist2 < Precision::SquareConfusion())
	369	{
	370	mySqDist.Append(mydist2);
	371	myismin.Append(Standard_True);
	372	mypoint.Append(ThePOnC(myusup, Pl));
	373	}
	374	mydone = Standard_True;
	375	}
	376	}
	377	else
	378	{
	379	// In analytical case
	380	mydone = myExtPElC.IsDone();
	381	if (mydone)
	382	{
	383	NbExt = myExtPElC.NbExt();
	384	for (i = 1; i <= NbExt; i++)
	385	{
	386	// Verification de la validite des parametres:
	387	ThePOnC PC = myExtPElC.Point(i);
	388	U = PC.Parameter();
	389	if (TheCurveTool::IsPeriodic(aCurve))
94f71cad	390	{
c8bf1eb7	391	U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(aCurve));
94f71cad	392	}
c8bf1eb7	393	if ((U >= myuinf-mytolu) && (U <= myusup+mytolu))
94f71cad	394	{
c8bf1eb7	395	PC.SetValues(U, myExtPElC.Point(i).Value());
	396	mySqDist.Append(myExtPElC.SquareDistance(i));
	397	myismin.Append(myExtPElC.IsMin(i));
	398	mypoint.Append(PC);
94f71cad	399	}
94f71cad	400	}
7fd59977	401	}
7fd59977	402	}
7fd59977	403	}
	404
	405
	406	//=======================================================================
	407	//function : Initialize
	408	//purpose :
	409	//=======================================================================
	410
	411	void Extrema_GExtPC::Initialize(const TheCurve& C,
	412	const Standard_Real Uinf,
	413	const Standard_Real Usup,
	414	const Standard_Real TolF)
	415	{
	416	myC = (Standard_Address)&C;
	417	myintuinf = myuinf = Uinf;
	418	myintusup = myusup = Usup;
	419	mytolf = TolF;
	420	mytolu = TheCurveTool::Resolution(((TheCurve)myC), Precision::Confusion());
	421	type = TheCurveTool::GetType(C);
	422	mydone = Standard_False;
	423	mydist1 = RealLast();
	424	mydist2 = RealLast();
	425	mysample = 17;
	426	}
	427
	428
	429	//=======================================================================
	430	//function : IntervalPerform
	431	//purpose :
	432	//=======================================================================
	433
	434	void Extrema_GExtPC::IntervalPerform(const ThePoint& P)
	435	{
	436	Standard_Integer i;
	437	Standard_Real U;
c8bf1eb7	438	myExtPC.Initialize(mysample, myintuinf, myintusup, mytolu, mytolf);
7fd59977	439	myExtPC.Perform(P);
7fd59977	440	mydone = myExtPC.IsDone();
c8bf1eb7	441	if (mydone)
c8bf1eb7	442	{
7fd59977	443	Standard_Integer NbExt = myExtPC.NbExt();
c8bf1eb7	444	for (i = 1; i <= NbExt; i++)
c8bf1eb7	445	{
7fd59977	446	// Verification de la validite des parametres pour le cas trimme:
	447	ThePOnC PC = myExtPC.Point(i);
	448	U = PC.Parameter();
c8bf1eb7	449	if (TheCurveTool::IsPeriodic(((TheCurve)myC)))
c8bf1eb7	450	{
253881cf	451	U = ElCLib::InPeriod(U, myuinf, myuinf+TheCurveTool::Period(((TheCurve)myC)));
7fd59977	452	}
c8bf1eb7	453	if ((U >= myuinf - mytolu) && (U <= myusup + mytolu))
c8bf1eb7	454	{
253881cf	455	PC.SetValues(U, PC.Value());
	456	mySqDist.Append(myExtPC.SquareDistance(i));
	457	myismin.Append(myExtPC.IsMin(i));
	458	mypoint.Append(PC);
7fd59977	459	}
	460	}
	461	}
	462	}
	463
	464
	465
	466
	467	//=======================================================================
	468	//function : Extrema_GExtPC
	469	//purpose :
	470	//=======================================================================
	471
	472	Extrema_GExtPC::Extrema_GExtPC()
	473	{
	474	myC = 0;
	475	mydone = Standard_False;
	476	mydist1 = RealLast();
	477	mydist2 = RealLast();
	478	mytolu = 0.0;
	479	mytolf = 0.0;
	480	mysample = 17;
	481	myintuinf = myintusup = myuinf = myusup = Precision::Infinite();
	482	type = GeomAbs_OtherCurve;
	483	}
	484
	485	//=======================================================================
	486	//function : Extrema_GExtPC
	487	//purpose :
	488	//=======================================================================
	489
	490	Extrema_GExtPC::Extrema_GExtPC(const ThePoint& P,
	491	const TheCurve& C,
	492	const Standard_Real Uinf,
	493	const Standard_Real Usup,
	494	const Standard_Real TolF)
	495	{
	496	Initialize(C, Uinf, Usup, TolF);
	497	Perform(P);
	498	}
	499
	500	//=======================================================================
	501	//function : Extrema_GExtPC
	502	//purpose :
	503	//=======================================================================
	504
	505	Extrema_GExtPC::Extrema_GExtPC(const ThePoint& P,
	506	const TheCurve& C,
	507	const Standard_Real TolF)
	508	{
	509	Initialize(C, TheCurveTool::FirstParameter(C),
	510	TheCurveTool::LastParameter(C), TolF);
	511	Perform(P);
	512	}
	513
	514
	515	//=======================================================================
	516	//function : IsDone
	517	//purpose :
	518	//=======================================================================
	519
	520	Standard_Boolean Extrema_GExtPC::IsDone() const
	521	{
	522	return mydone;
523	}
524
525
526	//=======================================================================
527	//function : Value
528	//purpose :
529	//=======================================================================
530
531	Standard_Real Extrema_GExtPC::SquareDistance(const Standard_Integer N) const
532	{
638ad7f3	533	if ((N < 1) \|\| (N > NbExt())) throw Standard_OutOfRange();
7fd59977	534	return mySqDist.Value(N);
	535	}
	536
	537
	538	//=======================================================================
	539	//function : NbExt
	540	//purpose :
	541	//=======================================================================
	542
	543	Standard_Integer Extrema_GExtPC::NbExt() const
	544	{
638ad7f3	545	if (!IsDone()) throw StdFail_NotDone();
7fd59977	546	return mySqDist.Length();
	547	}
	548
	549
	550	//=======================================================================
	551	//function : IsMin
	552	//purpose :
	553	//=======================================================================
	554
	555	Standard_Boolean Extrema_GExtPC::IsMin(const Standard_Integer N) const
	556	{
638ad7f3	557	if ((N < 1) \|\| (N > NbExt())) throw Standard_OutOfRange();
7fd59977	558	return myismin.Value(N);
	559	}
	560
	561
	562
	563	//=======================================================================
	564	//function : Point
	565	//purpose :
	566	//=======================================================================
	567
5d99f2c8	568	const ThePOnC & Extrema_GExtPC::Point(const Standard_Integer N) const
7fd59977	569	{
638ad7f3	570	if ((N < 1) \|\| (N > NbExt())) throw Standard_OutOfRange();
7fd59977	571	return mypoint.Value(N);
	572	}
	573
	574
	575	//=======================================================================
	576	//function : TrimmedDistances
	577	//purpose :
	578	//=======================================================================
	579
	580	void Extrema_GExtPC::TrimmedSquareDistances(Standard_Real& dist1,
	581	Standard_Real& dist2,
	582	ThePoint& P1,
	583	ThePoint& P2) const
	584	{
	585	dist1 = mydist1;
	586	dist2 = mydist2;
	587	P1 = Pf;
	588	P2 = Pl;
	589	}