// This file is part of libigl, a simple c++ geometry processing library.
// 
// Copyright (C) 2026 Alec Jacobson <alecjacobson@gmail.com>
// 
// This Source Code Form is subject to the terms of the Mozilla Public License 
// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
// obtain one at http://mozilla.org/MPL/2.0/.
#include "box_cubic.h"
#include "../cubic.h"
#include "../parallel_for.h"
#include <cyPolynomial.h>

template < 
  typename DerivedC, 
  typename DerivedB>
IGL_INLINE void igl::cycodebase::box_cubic(
  const Eigen::MatrixBase<DerivedC>& C,
  Eigen::PlainObjectBase<DerivedB>& B1,
  Eigen::PlainObjectBase<DerivedB>& B2)
{
  using Scalar = typename DerivedC::Scalar;
  typedef Eigen::Matrix<Scalar,3,1> VectorS3;
  typedef Eigen::Matrix<Scalar,DerivedC::RowsAtCompileTime,1> VectorSC;

  // Using the control points is a simple (but not tight) bound
  B1 = C({0,3},Eigen::all).colwise().minCoeff();
  B2 = C({0,3},Eigen::all).colwise().maxCoeff();

  // Better to find each of the t values where dC/dt = 0 and evaluate C there
  for(int d = 0;d<C.cols();d++)
  {
    // [3*C1 - 3*C0, 6*C0 - 12*C1 + 6*C2, 9*C1 - 3*C0 - 9*C2 + 3*C3]
    VectorS3 coef(
      Scalar(3) * (C(1,d) - C(0,d)),
      Scalar(6) * (C(0,d) - Scalar(2) * C(1,d) + C(2,d)),
      Scalar(3) * (C(3,d) - C(0,d) + Scalar(3) * (C(1,d) - C(2,d)))
    );
    VectorS3 r;
    int nr = cy::PolynomialRoots<2>(r.data(), coef.data(),0.0,1.0);
    for(int i = 0;i<nr;i++)
    {
      const Scalar t = r[i];
      // Evaluate cubic at t
      Eigen::Matrix<Scalar,1,1> Cdt;
      igl::cubic(C.col(d),t,Cdt);
      B1(d) = std::min(B1(d),Cdt(0));
      B2(d) = std::max(B2(d),Cdt(0));
    }
  }
}

template < 
  typename DerivedP, 
  typename DerivedC, 
  typename DerivedB>
IGL_INLINE void igl::cycodebase::box_cubic(
  const Eigen::MatrixBase<DerivedP>& P,
  const Eigen::MatrixBase<DerivedC>& C,
  Eigen::PlainObjectBase<DerivedB>& B1,
  Eigen::PlainObjectBase<DerivedB>& B2)
{
  B1.resize(C.rows(),P.cols());
  B2.resize(C.rows(),P.cols());
  typedef Eigen::Matrix<typename DerivedP::Scalar,1,DerivedP::ColsAtCompileTime> RowVectorP;

  igl::parallel_for(C.rows(),[&](const int c)
  {
    RowVectorP B1_c, B2_c;
    // Eval copies, but is it really good to make a template for the non copy?
    box_cubic(P(C.row(c),Eigen::all).eval(),B1_c,B2_c);
    B1.row(c) = B1_c;
    B2.row(c) = B2_c;
  },1000);
}

#ifdef IGL_STATIC_LIBRARY
// Explicit template instantiation
// generated by autoexplicit.sh
template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 1, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 1, -1, -1>>&);
template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 1, -1, 1, 1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1>>&);
template void igl::cycodebase::box_cubic<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>>&);
#endif