射影変換のお勉強

はじめに

射影変換はある平面を別の平面に射影することができる変換です。斜めから見たものを、もし正面から見たらどうなるかを計算できます。

変換式

変換前の座標(x,y)を(x',y')に変換するための行列Hを求めることが目的です。
（参考：CiNii 論文 -  反復を要しない射影変換の高精度解法）

x'とHxは平行になるので、x' × Hx=0と書けます。

とすると、内積

に変形できます。

最小二乗法

以下の式のJを最小とするhを求めることが最小二乗法になります。

これは２次形式なので、行列Mの最小固有値に対する固有ベクトルがhとなります。
実行結果が以下の図です。４点マウスで選択して左の画像を変換すると下の画像になります。

void Homography::Process(List<Point^>^ list1, List<Point^>^ list2)
{
  if (list1->Count != list2->Count)
    throw gcnew Exception("Data count unmatched");

  const double f0 = 1.0;
  if (is_first)
  {
     h = CalcH(list1, list2, f0); // 最小二乗法でHを計算

     h.resize(3, 3);
     MatrixXd temp = h.transpose();
     h = temp;	
     is_first = false;
  }

  cv::Mat_<cv::Vec3b> camera_image = usbcamra_->GrabImage();

  cv::Mat_<double> cvh;	
  cv::eigen2cv(h, cvh);

  cv::Mat_<cv::Vec3b> image;
  cv::warpPerspective(camera_image, image, cvh, image.size());

  cv::imshow("Result", image);
}

MatrixXd Homography::CalcH(List<Point^>^ list1, List<Point^>^ list2, double f0)
{
  const int N = list1->Count;
  MatrixXd mls = MatrixXd::Zero(9, 9);
  for (int i = 0; i < N; i++)
  {
    const double x1 = list1[i]->X;
    const double y1 = list1[i]->Y;
    const double x2 = list2[i]->X;
    const double y2 = list2[i]->Y;

    VectorXd xi[3] = {VectorXd::Zero(9), VectorXd::Zero(9), VectorXd::Zero(9)};
    xi[0] << 0, 0, 0, -f0*x1, -f0*y1, -f0*f0, x1*y2, y1*y2, f0*y2;
    xi[1] << f0*x1, f0*y1, f0*f0, 0, 0, 0, -x1*x2, -y1*x2, -f0*x2;
    xi[2] << -x1*y2, -y1*y2, -f0*y2, x1*x2, y1*x2, f0*x2, 0, 0, 0; 

    for (int j = 0; j < 3; j++)
      mls += xi[j] * xi[j].transpose();
  }
  mls /= N;

  // 最小固有ベクトルを求めて、最小２乗法をとく
  int min_eigen_index = 0;
  EigenSolver<MatrixXd> es(mls);
  es.eigenvalues().real().minCoeff(&min_eigen_index);
  return es.eigenvectors().col(min_eigen_index).real();
}