[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 <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 param space.
      Standard_Real aPeriodJump = 0.0;
      if (TheCurveTool::IsPeriodic(aCurve))
      {
        Standard_Integer aPeriodShift =
          Standard_Integer ((myuinf - aKnots(aFirstIdx)) / TheCurveTool::Period(aCurve));
        
        if (myuinf < aKnots(aFirstIdx))
          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)
          aFirstUsedKnot = anIdx;
        else
          break;

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

      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 first and last interval.
      if (mydist1 > Precision::SquareConfusion())
      {
        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

        // Derivatives have opposite signs - min or max inside of interval (sufficient condition).
        // Necessary condition - when point lies on curve.
        if(aVal1 * aVal2 <= 0.0 ||
           aBase1.Dot(aBase2) <= 0.0)
        {
          myintuinf = aParam(aVal.Lower());
          myintusup = aParam(aVal.Lower() + 1);
          IntervalPerform(P);
        }
      }

      if (mydist2 > Precision::SquareConfusion())
      {
        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

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

      mydone = Standard_True;
      break;
    }
  case GeomAbs_OtherCurve:
    {
      Standard_Boolean IntExtIsDone = Standard_False;
      Standard_Boolean IntIsNotValid;
      n = TheCurveTool::NbIntervals(aCurve, GeomAbs_C2);
      TColStd_Array1OfReal theInter(1, n+1);
      Standard_Boolean isPeriodic = TheCurveTool::IsPeriodic(aCurve);
      TheCurveTool::Intervals(aCurve, theInter, GeomAbs_C2);
      mysample = Max(mysample/n, 17);
      TheVector V1;
      ThePoint PP;
      Standard_Real s1 = 0.0 ;
      Standard_Real s2 = 0.0;
      myExtPC.Initialize(aCurve);
      for (i = 1; i <= n; i++)
      {
        myintuinf = theInter(i);
        myintusup = theInter(i+1);
        
        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_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(!mydone) StdFail_NotDone::Raise();
  if ((N < 1) || (N > mySqDist.Length())) Standard_OutOfRange::Raise();
  return mySqDist.Value(N);
}


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

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


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

Standard_Boolean Extrema_GExtPC::IsMin(const Standard_Integer N) const
{
  if(!mydone) StdFail_NotDone::Raise();
  if ((N < 1) || (N > mySqDist.Length())) Standard_OutOfRange::Raise();
  return myismin.Value(N);
}


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

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