验证【基于pca计算多边形的主轴】是否可行

结论：不可行，待补充。。。

过程：

#include <iostream>
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
#include <vector>
#include <matplot/matplot.h>


// 利用PCA方法计算多边形的主轴
class PrincipalAxesCalc {
 public:
  static void demean(Eigen::Vector2d &mean, Eigen::MatrixXd &centered, const Eigen::MatrixXd &pts_matrix) {
    mean = pts_matrix.colwise().mean();
    centered = pts_matrix.rowwise() - mean.transpose();
  }

  static bool Calc(const std::vector<Eigen::Vector2d>& pts) {
    if (pts.size() < 3) {
      return false;
    }
    Eigen::MatrixXd pts_matrix(pts.size(), 2);
    for (size_t i = 0; i < pts.size(); ++i) {
      pts_matrix(i, 0) = pts.at(i).x();
      pts_matrix(i, 1) = pts.at(i).y();
    }
    Eigen::Vector2d mean;
    Eigen::MatrixXd centered;
    demean(mean, centered, pts_matrix);
    Eigen::MatrixXd cov = (centered.adjoint() * centered) / (double )pts_matrix.rows();
    std::cout << "original pts matrix is:\n" << pts_matrix << "\n";
    std::cout << "col mean is:\n" << mean << "\n";
    std::cout << "centered is:\n" << centered << "\n";
    std::cout << "cov is:\n" << cov << "\n";

    // draw:
    auto fig_pca = matplot::figure()->current_axes();
    fig_pca->clear();
    fig_pca->axis(matplot::equal);
    fig_pca->hold(matplot::on);
    fig_pca->xlim({-10, 10});
    fig_pca->ylim({-10, 10});
    fig_pca->xlabel("x");
    fig_pca->ylabel("y");
    std::vector<double> ticks;
    for (int i = -10; i < 11; ++i) {
      ticks.push_back(i * 1.0);
    }
    fig_pca->xticks(ticks);
    fig_pca->yticks(ticks);
    std::vector<double> xs;
    std::vector<double> ys;
    for (size_t i = 0; i < pts_matrix.rows(); ++i) {
      xs.push_back(pts_matrix(i, 0));
      ys.push_back(pts_matrix(i, 1));
    }
    fig_pca->plot(xs, ys, "sb--");

    Eigen::EigenSolver<Eigen::MatrixXd> es(cov);
    Eigen::MatrixXd D = es.pseudoEigenvalueMatrix().diagonal();
    Eigen::MatrixXd V = es.pseudoEigenvectors();
    std::cout << "The pseudo-eigenvalue matrix D is:" << "\n" << D << "\n";
    std::cout << "The pseudo-eigenvalue matrix V is:" << "\n" << V << "\n";
    std::vector<size_t> sorted_index(D.rows());
    std::iota(sorted_index.begin(), sorted_index.end(), 0);
    std::sort(sorted_index.begin(), sorted_index.end(), [&](size_t i1, size_t i2) { return D(i1) > D(i2);});
    for (size_t i = 0; i < sorted_index.size(); ++i) {
      std::cout << "the " << i << "th value is:" << sorted_index.at(i) << "\n";
    }

    Eigen::MatrixXd W = V.col(sorted_index.at(0)) * V.col(sorted_index.at(0)).transpose();
    std::cout << "The W is:\n" << W << "\n";
    Eigen::MatrixXd pca_matrix = (centered * W).rowwise() + mean.transpose();
    std::cout << "The pca matrix is:\n" << pca_matrix << "\n";

    // draw:
    xs.clear();
    ys.clear();
    for (size_t i = 0; i < pca_matrix.rows(); ++i) {
      xs.push_back(pca_matrix(i, 0));
      ys.push_back(pca_matrix(i, 1));
    }
    fig_pca->plot(xs, ys, "dr--");

    fig_pca->arrow(mean.x(), mean.y(), W(0, 0), W(1, 0))->color("black");

    for (size_t i = 0; i < pca_matrix.rows(); ++i) {
      fig_pca->arrow(pts_matrix(i, 0), pts_matrix(i, 1), pca_matrix(i, 0), pca_matrix(i, 1))->color("green");
    }

    return true;
  }
};

int main(int argc, char **argv) {
  std::vector<Eigen::Vector2d> pts;
  Eigen::Vector2d beg_pt;
  Eigen::Vector2d end_pt;
  int test_case = 1;
  if (test_case == 0) {
    pts.push_back({-2.0, -1.0});
    pts.push_back({ 2.0, -1.0});
    pts.push_back({ 2.0,  1.0});
    pts.push_back({-2.0,  1.0});
  } else if (test_case == 1) {
    pts.push_back({ 0.0, -4.0});
    pts.push_back({ 4.0,  0.0});
    pts.push_back({ 2.0,  2.9});
    pts.push_back({ 0.0,  9.0});
    pts.push_back({-2.0,  0.0});
  } else if (test_case == 2) {
    pts.push_back({2.5, 2.4});
    pts.push_back({0.5, 0.7});
    pts.push_back({2.2, 2.9});
    pts.push_back({1.9, 2.2});
    pts.push_back({3.1, 3.0});
    pts.push_back({2.3, 2.7});
    pts.push_back({2.0, 1.6});
    pts.push_back({1.0, 1.1});
    pts.push_back({1.5, 1.6});
    pts.push_back({1.1, 0.9});
  } else if (test_case == 3) {
    pts.push_back({-1.0, -2.0});
    pts.push_back({-1.0,  2.0});
    pts.push_back({ 1.0,  2.0});
    pts.push_back({ 1.0, -2.0});
  }
  PrincipalAxesCalc::Calc(pts);
  return 0;
}