// Created on: 1997-05-13
// Created by: Stagiaire Francois DUMONT
// Copyright (c) 1997-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 <BSplSLib.hxx>
#include <Geom_BSplineSurface.hxx>
#include <GeomLib_DenominatorMultiplier.hxx>
#include <gp_Pnt.hxx>
#include <gp_Vec.hxx>
#include <math_Matrix.hxx>
#include <TColgp_Array2OfPnt.hxx>
#include <TColStd_Array1OfInteger.hxx>
#include <TColStd_Array1OfReal.hxx>
#include <TColStd_Array2OfReal.hxx>

//=================================================================================================

GeomLib_DenominatorMultiplier::GeomLib_DenominatorMultiplier(
  const Handle(Geom_BSplineSurface)& Surface,
  const TColStd_Array1OfReal&        KnotVector)
    : mySurface(Surface),
      myKnotFlatVector(1, KnotVector.Length())
{
  Standard_Integer i;
  for (i = 1; i <= KnotVector.Length(); i++)
    myKnotFlatVector.SetValue(i, KnotVector(i));
}

//=======================================================================
// function : value
// purpose  : give the value of a(UParameter,VParameter)
//=======================================================================

Standard_Real GeomLib_DenominatorMultiplier::Value(const Standard_Real UParameter,
                                                   const Standard_Real VParameter) const

{
  Standard_Real Dumaxv, Duminv, dDduumaxv, dDduuminv, Dv, Buv = 0.0;
  // gp_Pnt         HermPnt;
  gp_Pnt                  N;
  gp_Vec                  Nu, Nv;
  TColgp_Array2OfPnt      surface_poles(1, mySurface->NbUPoles(), 1, mySurface->NbVPoles());
  TColStd_Array2OfReal    surface_weights(1, mySurface->NbUPoles(), 1, mySurface->NbVPoles());
  TColStd_Array1OfReal    surface_u_knots(1, mySurface->NbUKnots());
  TColStd_Array1OfInteger surface_u_mults(1, mySurface->NbUKnots());

  TColStd_Array1OfReal    surface_v_knots(1, mySurface->NbVKnots());
  TColStd_Array1OfInteger surface_v_mults(1, mySurface->NbVKnots());
  Standard_Integer        udegree, vdegree;

  mySurface->UKnots(surface_u_knots);
  mySurface->UMultiplicities(surface_u_mults);
  mySurface->Poles(surface_poles);
  mySurface->Weights(surface_weights);
  mySurface->VKnots(surface_v_knots);
  mySurface->VMultiplicities(surface_v_mults);
  udegree = mySurface->UDegree();
  vdegree = mySurface->VDegree();

  BSplSLib::HomogeneousD1(mySurface->UKnot(mySurface->LastUKnotIndex()),
                          VParameter,
                          0,
                          0,
                          surface_poles,
                          &surface_weights,
                          surface_u_knots,
                          surface_v_knots,
                          &surface_u_mults,
                          &surface_v_mults,
                          udegree,
                          vdegree,
                          mySurface->IsURational(),
                          mySurface->IsVRational(),
                          mySurface->IsUPeriodic(),
                          mySurface->IsVPeriodic(),
                          N,
                          Nu,
                          Nv,
                          Dumaxv,
                          dDduumaxv,
                          Dv);
  BSplSLib::HomogeneousD1(mySurface->UKnot(1),
                          VParameter,
                          0,
                          0,
                          surface_poles,
                          &surface_weights,
                          surface_u_knots,
                          surface_v_knots,
                          &surface_u_mults,
                          &surface_v_mults,
                          udegree,
                          vdegree,
                          mySurface->IsURational(),
                          mySurface->IsVRational(),
                          mySurface->IsUPeriodic(),
                          mySurface->IsVPeriodic(),
                          N,
                          Nu,
                          Nv,
                          Duminv,
                          dDduuminv,
                          Dv);

  math_Matrix   BSplineBasisDeriv(1, 2, 1, 4, 0.0);
  Standard_Real B1prim0, Bprelastprim1,
    lambda = (mySurface->Weight(1, 1) / mySurface->Weight(mySurface->NbUPoles(), 1));
  Standard_Integer index, i;

  BSplCLib::EvalBsplineBasis(1, 4, myKnotFlatVector, 0.0, index, BSplineBasisDeriv);
  B1prim0 = BSplineBasisDeriv(2, 2);

  BSplCLib::EvalBsplineBasis(1, 4, myKnotFlatVector, 1.0, index, BSplineBasisDeriv);
  Bprelastprim1 = BSplineBasisDeriv(2, 3);

  math_Matrix BSplineBasisValue(1, 1, 1, 4, 0.0);
  BSplCLib::EvalBsplineBasis(0, 4, myKnotFlatVector, UParameter, index, BSplineBasisValue);

  TColStd_Array1OfReal value(0, 5);
  TColStd_Array1OfReal Polesenv(0, 5); // poles of a(u,v)

  for (i = 0; i <= 5; i++)
    Polesenv(i) = 0.0;
  Polesenv(0) = Duminv;
  Polesenv(1) = Duminv - dDduuminv / B1prim0;
  Polesenv(4) = lambda * lambda * (Dumaxv - dDduumaxv / Bprelastprim1);
  Polesenv(5) = lambda * lambda * Dumaxv;

  if (myKnotFlatVector.Length() == 8)
  {
    value(0) = BSplineBasisValue(1, 1); // values of the basic functions
    value(1) = BSplineBasisValue(1, 2);
    value(2) = 0.0;
    value(3) = 0.0;
    value(4) = BSplineBasisValue(1, 3);
    value(5) = BSplineBasisValue(1, 4);
  }
  if (myKnotFlatVector.Length() == 9)
  {
    if (index == 1)
    {
      value(0) = BSplineBasisValue(1, 1);
      value(1) = BSplineBasisValue(1, 2);
      value(2) = BSplineBasisValue(1, 3);
      value(3) = 0.0;
      value(4) = BSplineBasisValue(1, 4);
      value(5) = 0.0;
    }
    else
    {
      value(0) = 0.0;
      value(1) = BSplineBasisValue(1, 1);
      value(2) = BSplineBasisValue(1, 2);
      value(3) = 0.0;
      value(4) = BSplineBasisValue(1, 3);
      value(5) = BSplineBasisValue(1, 4);
    }
    Polesenv(2) = (0.5 * (Polesenv(0) + Polesenv(5)));
  }
  if (myKnotFlatVector.Length() == 10)
  {
    if (index == 1)
    {
      value(0) = BSplineBasisValue(1, 1);
      value(1) = BSplineBasisValue(1, 2);
      value(2) = BSplineBasisValue(1, 3);
      value(3) = BSplineBasisValue(1, 4);
      value(4) = 0.0;
      value(5) = 0.0;
    }
    if (index == 2)
    {
      value(0) = 0.0;
      value(1) = BSplineBasisValue(1, 1);
      value(2) = BSplineBasisValue(1, 2);
      value(3) = BSplineBasisValue(1, 3);
      value(4) = BSplineBasisValue(1, 4);
      value(5) = 0.0;
    }
    if (index == 3)
    {
      value(0) = 0.0;
      value(1) = 0.0;
      value(2) = BSplineBasisValue(1, 1);
      value(3) = BSplineBasisValue(1, 2);
      value(4) = BSplineBasisValue(1, 3);
      value(5) = BSplineBasisValue(1, 4);
    }
    Polesenv(2) = (0.5 * (Polesenv(0) + Polesenv(5)));
    Polesenv(3) = Polesenv(2);
  }
  for (i = 0; i <= 5; i++)
    Buv += Polesenv(i) * value(i);
  return Buv;
}