<commit_before>/*
 * Copyright 2016, Simula Research Laboratory
 *
 * 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 <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/iterator/indirect_iterator.hpp>
#include <cctag/EdgePoint.hpp>
#include <cctag/Fitting.hpp>
#include <cctag/utils/Defines.hpp>
#include <Eigen/SVD>
#include <Eigen/LU>
#include <cctag/geometry/Ellipse.hpp>
#include <cctag/geometry/Distance.hpp>
#include <cctag/geometry/EllipseFromPoints.hpp>
#include <cctag/geometry/Point.hpp>
#include <cctag/Fitting.hpp>
#include <cmath>
#include <cfloat>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <utility>
#include <Eigen/Eigenvalues>

namespace cctag {
namespace numerical {

namespace geometry
{

using Vector6f = Eigen::Matrix<float, 6, 1>;
using Conic = std::tuple<Vector6f, Eigen::Vector2f>; // 6 coefficients + center offset

template<typename It>
static Eigen::Vector2f get_offset(It begin, It end)
{
  Eigen::Vector2f center(0, 0);
  const size_t n = end - begin;
  for (; begin != end; ++begin)
    center += Eigen::Vector2f((*begin)(0), (*begin)(1));
  return center / n;
}

template<typename It>
static std::tuple<Eigen::Matrix3f,Eigen::Matrix3f,Eigen::Matrix3f>
get_scatter_matrix(It begin, It end, const Eigen::Vector2f& offset)
{
  using namespace Eigen;
  
  const auto qf = [&](It it) {
    Vector2f p((*it)(0), (*it)(1));
    auto pc = p - offset;
    return Vector3f(pc(0)*pc(0), pc(0)*pc(1), pc(1)*pc(1));
  };
  const auto lf = [&](It it) {
    Vector2f p((*it)(0), (*it)(1));
    auto pc = p - offset;
    return Vector3f(pc(0), pc(1), 1);
  };
  
  const size_t n = end - begin;
  MatrixX3f D1(n,3), D2(n,3);
  
  // Construct the quadratic and linear parts.  Doing it in two loops has better cache locality.
  // TODO@stian: Make an nx6 matrix and define D1 and D2 as subblocks so that they're as one memory block
  for (size_t i = 0; begin != end; ++begin, ++i) {
    D1.row(i) = qf(begin);
    D2.row(i) = lf(begin);
  }
  
  // Construct the three parts of the symmetric scatter matrix.
  Matrix3f S1 = D1.transpose() * D1;
  Matrix3f S2 = D1.transpose() * D2;
  Matrix3f S3 = D2.transpose() * D2;
  return std::make_tuple(S1, S2, S3);
}

template<typename It>
static Conic fit_solver(It begin, It end)
{
  using namespace Eigen;
  using std::get;
  
  static const struct C1_Initializer {
    Matrix3f matrix;
    Matrix3f inverse;
    C1_Initializer()
    {
      matrix <<
          0,  0, 2,
          0, -1, 0,
          2,  0, 0;
      inverse <<
            0,  0, 0.5,
            0, -1,   0,
          0.5,  0,   0; 
    };
  } C1;
  
  const auto offset = get_offset(begin, end);
  const auto St = get_scatter_matrix(begin, end, offset);
  const auto& S1 = std::get<0>(St);
  const auto& S2 = std::get<1>(St);
  const auto& S3 = std::get<2>(St);
  bool invertible;
  Matrix3f S3Inv;
  S3.computeInverseWithCheck(S3Inv, invertible);
  if(!invertible)
  {
      throw std::domain_error("fit_solver: the input points appear to be linearly dependent");
  }
  const auto T = -S3.inverse() * S2.transpose();
  const auto M = C1.inverse * (S1 + S2*T);
  
  EigenSolver<Matrix3f> M_ev(M);
  Vector3f cond;
  {
    const auto evr = M_ev.eigenvectors().real().array();
    cond = 4*evr.row(0)*evr.row(2) - evr.row(1)*evr.row(1);
  }

    const auto eps = std::numeric_limits<float>::epsilon();
  float min = std::numeric_limits<float>::max();
  int imin = -1;
  for (int i = 0; i < 3; ++i)
  {
      if (cond(i) > eps && cond(i) < min)
      {
          imin = i;
          min = cond(i);
      }
  }
  
  Vector6f ret = Matrix<float, 6, 1>::Zero();
  if (imin == -1)
  {
      throw std::domain_error("fit_solver: degeneracy");
  }
    Vector3f a1 = M_ev.eigenvectors().real().col(imin);
    Vector3f a2 = T * a1;
    ret.block<3, 1>(0, 0) = a1;
    ret.block<3, 1>(3, 0) = a2;
    return std::make_tuple(ret, offset);
}

// Adapted from OpenCV old code; see
// https://github.com/Itseez/opencv/commit/4eda1662aa01a184e0391a2bb2e557454de7eb86#diff-97c8133c3c171e64ea0df0db4abd033c
void to_ellipse(const Conic& conic, Ellipse& ellipse)
{
  using namespace Eigen;
  const auto eps = std::numeric_limits<float>::epsilon();

  auto coef = std::get<0>(conic);

  float idet = coef(0)*coef(2) - coef(1)*coef(1)/4; // ac-b^2/4<commit_msg>[cctag] renaming to avoid conflict with min<commit_after>/*
 * Copyright 2016, Simula Research Laboratory
 *
 * 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 <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/split.hpp>
#include <boost/foreach.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/special_functions/round.hpp>
#include <boost/iterator/indirect_iterator.hpp>
#include <cctag/EdgePoint.hpp>
#include <cctag/Fitting.hpp>
#include <cctag/utils/Defines.hpp>
#include <Eigen/SVD>
#include <Eigen/LU>
#include <cctag/geometry/Ellipse.hpp>
#include <cctag/geometry/Distance.hpp>
#include <cctag/geometry/EllipseFromPoints.hpp>
#include <cctag/geometry/Point.hpp>
#include <cctag/Fitting.hpp>
#include <cmath>
#include <cfloat>
#include <fstream>
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <vector>
#include <utility>
#include <Eigen/Eigenvalues>

namespace cctag {
namespace numerical {

namespace geometry
{

using Vector6f = Eigen::Matrix<float, 6, 1>;
using Conic = std::tuple<Vector6f, Eigen::Vector2f>; // 6 coefficients + center offset

template<typename It>
static Eigen::Vector2f get_offset(It begin, It end)
{
  Eigen::Vector2f center(0, 0);
  const size_t n = end - begin;
  for (; begin != end; ++begin)
    center += Eigen::Vector2f((*begin)(0), (*begin)(1));
  return center / n;
}

template<typename It>
static std::tuple<Eigen::Matrix3f,Eigen::Matrix3f,Eigen::Matrix3f>
get_scatter_matrix(It begin, It end, const Eigen::Vector2f& offset)
{
  using namespace Eigen;
  
  const auto qf = [&](It it) {
    Vector2f p((*it)(0), (*it)(1));
    auto pc = p - offset;
    return Vector3f(pc(0)*pc(0), pc(0)*pc(1), pc(1)*pc(1));
  };
  const auto lf = [&](It it) {
    Vector2f p((*it)(0), (*it)(1));
    auto pc = p - offset;
    return Vector3f(pc(0), pc(1), 1);
  };
  
  const size_t n = end - begin;
  MatrixX3f D1(n,3), D2(n,3);
  
  // Construct the quadratic and linear parts.  Doing it in two loops has better cache locality.
  // TODO@stian: Make an nx6 matrix and define D1 and D2 as subblocks so that they're as one memory block
  for (size_t i = 0; begin != end; ++begin, ++i) {
    D1.row(i) = qf(begin);
    D2.row(i) = lf(begin);
  }
  
  // Construct the three parts of the symmetric scatter matrix.
  Matrix3f S1 = D1.transpose() * D1;
  Matrix3f S2 = D1.transpose() * D2;
  Matrix3f S3 = D2.transpose() * D2;
  return std::make_tuple(S1, S2, S3);
}

template<typename It>
static Conic fit_solver(It begin, It end)
{
  using namespace Eigen;
  using std::get;
  
  static const struct C1_Initializer {
    Matrix3f matrix;
    Matrix3f inverse;
    C1_Initializer()
    {
      matrix <<
          0,  0, 2,
          0, -1, 0,
          2,  0, 0;
      inverse <<
            0,  0, 0.5,
            0, -1,   0,
          0.5,  0,   0; 
    };
  } C1;
  
  const auto offset = get_offset(begin, end);
  const auto St = get_scatter_matrix(begin, end, offset);
  const auto& S1 = std::get<0>(St);
  const auto& S2 = std::get<1>(St);
  const auto& S3 = std::get<2>(St);
  bool invertible;
  Matrix3f S3Inv;
  S3.computeInverseWithCheck(S3Inv, invertible);
  if(!invertible)
  {
      throw std::domain_error("fit_solver: the input points appear to be linearly dependent");
  }
  const auto T = -S3.inverse() * S2.transpose();
  const auto M = C1.inverse * (S1 + S2*T);
  
  EigenSolver<Matrix3f> M_ev(M);
  Vector3f cond;
  {
    const auto evr = M_ev.eigenvectors().real().array();
    cond = 4*evr.row(0)*evr.row(2) - evr.row(1)*evr.row(1);
  }

  const auto eps = std::numeric_limits<float>::epsilon();
  float minValue = std::numeric_limits<float>::max();
  int imin = -1;
  for (int i = 0; i < 3; ++i)
  {
      if (cond(i) > eps && cond(i) < minValue)
      {
          imin = i;
          minValue = cond(i);
      }
  }
  
  Vector6f ret = Matrix<float, 6, 1>::Zero();
  if (imin == -1)
  {
      throw std::domain_error("fit_solver: degeneracy");
  }
    Vector3f a1 = M_ev.eigenvectors().real().col(imin);
    Vector3f a2 = T * a1;
    ret.block<3, 1>(0, 0) = a1;
    ret.block<3, 1>(3, 0) = a2;
    return std::make_tuple(ret, offset);
}

// Adapted from OpenCV old code; see
// https://github.com/Itseez/opencv/commit/4eda1662aa01a184e0391a2bb2e557454de7eb86#diff-97c8133c3c171e64ea0df0db4abd033c
void to_ellipse(const Conic& conic, Ellipse& ellipse)
{
  using namespace Eigen;
  const auto eps = std::numeric_limits<float>::epsilon();

  auto coef = std::get<0>(conic);

  float idet = coef(0)*coef(2) - coef(1)*coef(1)/4; // ac-b^2/4
{'fraction_non_alphanumeric': 0.15588266440749762, 'fraction_numerical': 0.039817854495393414, 'mean_word_length': 3.4864608076009502, 'pattern_counts': {'":': 0, '<': 101, '<?xml version=': 0, '>': 91, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5071962', 'n_tokens_mistral': 3791, 'n_tokens_neox': 3629, 'n_words': 905}