// 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;
      // Avoid prolem 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 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;
}