[MVG] 2차원 사영변환(Homography) 추정 (+ CPP 코드)

haeryong·2023년 11월 24일

Ceres-Solver MultipleViewGeometry findHomography

Multiple View Geometry 4장 "2차원 사영변환의 추정"에 대한 내용과 OpenCV, Ceres-Solver 라이브러리를 이용해 이를 직접 구현한 C++ 코드를 포함하고 있다.

1. 2차원 사영변환(Homography)

2차원에서 2차원으로의 점 대응 $x_i\leftrightarrow x_i'$ 이 존재할 때 점 대응 사이의 변환을 H(homography)라 하며 $x_i'=Hx_i$ 로 나타낼 수 있다.

여기서 $x_i=(u_i,v_i,w_i)^T, x_i'=(u_i',v_i',w_i')^T$ 은 2차원 동차 좌표계로 3차원 벡터이며, 여기서는 편의를 위해 모든 점의 마지막 원소 $w= w' = 1$ 으로 선택한다.

Homography를 이용하면 아래와 같이 Perspective correction(사영 정정)을 수행할 수 있다.

그림 출처

H를 추정하기 위해서 다음과 같은 과정을 거쳐야한다.

DLT 또는 RANSAC 알고리즘을 이용해 초기값 $H_0$ 계산.
비선형 최적화를 통해 기하 거리를 최소화하는 H 계산.

$H\in R^{3\times3}$ 는 9개의 원소로 이루어져 있고, 1개의 scale을 제외하면 8의 자유도를 가진다. 그리고 1쌍의 점대응을 통해 2개의 방정식을 얻을 수 있으므로 최소 4쌍의 점 대응을 통해 H를 결정할 수 있다.

따라서 4쌍 이상의 점 대응이 주어진 경우에 우리는 RANSAC 알고리즘을 통해 4쌍의 점 대응을 선택하거나, DLT 알고리즘을 통해 $H_0$ 를 계산해야 한다. 여기서는 DLT 알고리즘을 사용하기로 한다.

2. Direct Linear Transform(DLT) 알고리즘

먼저, 한 쌍의 점대응을 통해 2개의 방정식을 얻는 부분에 대해 설명하겠다.

$x_i'=Hx_i$ 의 관계를 두 벡터의 외적 $x_i'\times Hx_i=0$ 으로 표현할 수 있다.

$H$ 의 i번째 행벡터를 $h_i=(h_{i1},h_{i2},h_{i3})$ 라 하면 두 벡터의 외적은 다음과 같다.

$x_i'\times Hx_i=\begin{pmatrix}v_i'h_3x_i-w_i'h_2x_i\\w_i'h_1x_i-u_i'h_3x_i\\u_i'h_2x_i-v_i'h_1x_i\end{pmatrix}$

행렬 H를 $h= \begin{pmatrix}h_1&h_2&h_3\end{pmatrix}^T \in R^{9\times1}$ 와 같이 벡터로 정의하고,

$h_jx_i=x_i^Th_j^T \in R$ 이므로 $h_j$ 부분을 뒤쪽에 배치하면 아래와 같은 방정식을 얻는다.

$\begin{pmatrix}0_{1\times3}&-w_i'x_i^T&v_i'x_i^T\\w_i'x_i^T&0_{1\times3}&-u_i'x_i^T\\-v_i'x_i^T&u_i'x_i^T&0_{1\times3}\end{pmatrix}h=0$

3개의 선형 방정식 중 2개만 선형 독립이기 때문에 마지막 행를 제외한 2개의 방정식을 얻게 된다.

마지막 행을 제외하고, $w=1$ 을 적용하면 최종적으로 한 쌍의 점 대응으로부터 $A_i \in R^{2\times9}$ 를 얻을 수 있다.

$A_i=\begin{pmatrix}0&0&0&-u_i&-v_i&-1&v_i'u_i&v_i'v_i&v_i'\\u_i&v_i&1&0&0&0&-u_i'u_i&-u_i'v_i&-u_i'\end{pmatrix}$

네 쌍의 점대응이 주어지면 4개의 $A_i$ 를 합쳐 $A \in R^{8\times9}$ 를 만들고, $Ah=0$ 을 풀면 H를 계산할 수 있다.

그러나 네 쌍 이상의 점대응이 주어지고, 이미지 좌표 측정에 노이즈가 존재한다고 할 때, $A \in R^{2n\times9}$ 의 $Ah=0$ 을 만족하는 $h$ 는 존재하지 않는다.

따라서 적절한 비용함수를 최소화하는 h의 근사해를 구해야하며, DLT 알고리즘에서는 $||h||=1$ 의 제약조건을 만족하면서 $||Ah||$ 을 최소화하는 h를 계산하게 된다.

결과적으로 h의 근사해는 $A^TA$ 의 가장 작은 고윳값(eigen value)에 해당하는 고유벡터 또는 SVD를 사용해서 구할 수 있는 A의 가장 작은 특이값(singular value)에 해당하는 특이벡터이다. 여기서는 SVD방법 사용하였다.

데이터 정규화

DLT 알고리즘을 수행할 때 데이터 정규화를 필수적으로 수행해야한다.
SVD를 이용해서 h의 근사해를 구하는 과정은 $A$ 와 차이가 최소가 되는 $\hat{A}$ 에 대한 $\hat{A}h=0$ 의 해를 구하는 것과 같다.

이미지 점의 좌표가 $x_i = (100,100,1)$ 라고 하면 $A_i$ 의 원소의 값은 대략적으로 $u_iu_i'=10000, \quad u_iw_i=100, \quad w_i'w_i=1$ 을 가질 것이다. 어떤 차원의 좌표가 곱해지느냐에 따라 크기 차이가 발생한다. $A와 \hat{A}$ 의 차이를 최소화하는 과정에서 $u_iu_i'$ 와 $w_i'w_i$ 가 곱해진 원소의 값을 동일하게 변화시킬 때, 비용함수는 동일한 반면, w_i쪽이 훨씬 더 이미지에 변화가 크다. 따라서 데이터의 각 좌표의 평균, 표준편차 등의 값을 정규화해야한다.

정규화된 점 $\tilde{x}$ 는 평균이 (0, 0)이고, 원점에서의 평균 거리가 $\sqrt{2}$ 이 되어야 하고, 이러한 변환을 만족하는 scale과 translation으로 이루어진 닮음 변환 T를 계산할 수 있다. MVG의 내용과 달리 OpenCV 라이브러리의 경우 평균이 (0, 0)이고, u좌표의 평균 크기가 1, v좌표의 평균 크기가 1이 되도록 scale_x, scale_y로 나누어서 T를 계산한다. 여기서는 OpenCV의 방법에 따라 구현하였다.

정규화된 데이터를 이용해서 DLT 알고리즘을 수행하면 정규화된 $\tilde{H}$ 를 구할 수 있고, 이것을 다시 denormalize 해야한다.

$x_i$ 를 $\tilde{x_i}$ 로, $x_i'$ 를 $\tilde{x_i}'$ 로 정규화하는 닮음 변환을 각각 $T, T' \in R^{3\times3}$ 라 할 때, $H=T'^{-1}\tilde{H}T$ 를 계산해 최종적인 H를 구한다.

3. 비선형 최적화

DLT 방법을 통해 ||Ah||를 최소화하였지만 이는 대수 오차를 최소화한 것으로, 비선형 최적화를 통해 기하 오차를 최소화 해주어야 한다.

비선형 최적화를 위해 여러가지 오차함수가 존재하지만 가장 간단한 이미지 하나에서의 오차를 사용하였다.

비용함수 : $\Sigma d(x_i',Hx_i)^2$

비선형 최적화를 통해 비용함수를 최소화하는 H를 찾게 된다.

반복 최소화를 위해 Gauss-Newton, Levenberg-Marquatt 방법 등을 사용하며, Ceres-Sovler를 이용해 최적화를 진행하였다.

4. 전체 C++ 코드

두 이미지 사이의 점 대응 구한 뒤 find_homography 함수를 호출해 H를 얻을 수 있다.

#include <opencv2/opencv.hpp>
#include <ceres/ceres.h>
#include <ceres/rotation.h>

// Homography를 추정하기 위한 비용함수가 정의되어 있음. 
// Ceres-Solver 비선형 최적화를 위해 필요하다.
struct homography_reproj_error
{
  homography_reproj_error(const cv::Point2f &src_point, const cv::Point2f &dst_point)
  : src_point(src_point), dst_point(dst_point)
  {
  }
  template<typename T>
  bool operator()(const T *const H, T *residuals) const
  {
    T pred_z = H[6] * T(src_point.x) + H[7] * T(src_point.y) + H[8];
    T dst_point_est[2];
    dst_point_est[0] = (H[0] * T(src_point.x) + H[1] * T(src_point.y) + H[2]) / pred_z;
    dst_point_est[1] = (H[3] * T(src_point.x) + H[4] * T(src_point.y) + H[5]) / pred_z;
    residuals[0] = dst_point_est[0] - T(dst_point.x);
    residuals[1] = dst_point_est[1] - T(dst_point.y);
    return true;
  }
  static ceres::CostFunction *create(const cv::Point2f &src_point, const cv::Point2f &dst_point)
  {
    return (new ceres::AutoDiffCostFunction<homography_reproj_error, 2, 9>(new homography_reproj_error(src_point, dst_point)));
  }
  const cv::Point2f src_point;
  const cv::Point2f dst_point;
};


// 점 대응이 존재할 때 homography 추정
cv::Mat find_homography(
  const std::vector<cv::Point2f> &src_points,
  const std::vector<cv::Point2f> &dst_points)
{
  std::size_t num_points = src_points.size();
  std::vector<cv::Point2f> normalized_src_points(src_points);
  std::vector<cv::Point2f> normalized_dst_points(dst_points);

  cv::Mat T_src = normalize_image_points(normalized_src_points);
  cv::Mat T_dst = normalize_image_points(normalized_dst_points);

  cv::Mat A = cv::Mat::zeros(2 * num_points, 9, CV_64F);

  for (std::size_t i = 0; i < num_points; ++i)
  {
    auto &src = normalized_src_points[i];
    auto &dst = normalized_dst_points[i];
    A.at<double>(2 * i, 0) = src.x;
    A.at<double>(2 * i, 1) = src.y;
    A.at<double>(2 * i, 2) = 1.0;
    A.at<double>(2 * i, 6) = -1.0 * dst.x * src.x;
    A.at<double>(2 * i, 7) = -1.0 * dst.x * src.y;
    A.at<double>(2 * i, 8) = -1.0 * dst.x;
    A.at<double>(2 * i + 1, 3) = src.x;
    A.at<double>(2 * i + 1, 4) = src.y;
    A.at<double>(2 * i + 1, 5) = 1.0;
    A.at<double>(2 * i + 1, 6) = -1.0 * dst.y * src.x;
    A.at<double>(2 * i + 1, 7) = -1.0 * dst.y * src.y;
    A.at<double>(2 * i + 1, 8) = -1.0 * dst.y;
  }
  cv::Mat S, U, Vt;
  cv::SVD::compute(A, S, U, Vt, cv::SVD::Flags::FULL_UV);
  cv::Mat H = cv::Mat(3, 3, CV_64F, Vt.row(Vt.rows - 1).data);
  H = T_dst.inv() * H * T_src;
  H /= H.at<double>(2, 2);

  refine_homography(H, src_points, dst_points);
  H /= H.at<double>(2, 2);

  return H;
}

// 데이터 정규화
cv::Mat normalize_image_points(std::vector<cv::Point2f>& points)
{
  cv::Mat T = cv::Mat::eye(3, 3, CV_64F);
  std::size_t num_points = points.size();

  double x_center = 0.0, y_center = 0.0;
  for(const auto& point : points)
  {
    x_center += static_cast<double>(point.x);
    y_center += static_cast<double>(point.y);
  }
  x_center /= static_cast<double>(num_points);
  y_center /= static_cast<double>(num_points);

  T.at<double>(0, 2) = -1.0 * x_center;
  T.at<double>(1, 2) = -1.0 * y_center;

  for(auto& point : points)
  {
    point.x -= static_cast<float>(x_center);
    point.y -= static_cast<float>(y_center);
  }

  double avg_dist_x = 0.0, avg_dist_y = 0.0;
  for(const auto& point : points)
  {
    avg_dist_x += static_cast<double>(std::fabs(point.x));
    avg_dist_y += static_cast<double>(std::fabs(point.y));
  }
  avg_dist_x /= static_cast<double>(num_points);
  avg_dist_y /= static_cast<double>(num_points);

  double scale_x = 1.0 / avg_dist_x;
  double scale_y = 1.0 / avg_dist_y;
  T.at<double>(0, 0) = scale_x;
  T.at<double>(1, 1) = scale_y;
  T.at<double>(0, 2) *= scale_x;
  T.at<double>(1, 2) *= scale_y;

  for(auto& point : points)
  {
    point.x *= static_cast<double>(scale_x);
    point.y *= static_cast<double>(scale_y);
  }

  return T;
}

// 비선형 최적화
void refine_homography(
  cv::Mat &H,
  const std::vector<cv::Point2f> &src_points,
  const std::vector<cv::Point2f> &dst_points)
{
  std::size_t num_points = src_points.size();
  std::vector<double> homography(9);
  for(std::size_t i = 0; i < 3; ++i)
  {
    for(std::size_t j = 0; j < 3; ++j)
    {
      homography[3 * i + j] = H.at<double>(i, j);
    }
  }

  auto problem = std::make_unique<ceres::Problem>();
  ceres::Solver::Options options;
  ceres::Solver::Summary summary;
  options.linear_solver_type = ceres::ITERATIVE_SCHUR;
  options.num_threads = 8;
  options.minimizer_progress_to_stdout = true;
  options.parameter_tolerance = 1e-12;
  options.gradient_tolerance = 1e-14;
  options.function_tolerance = 1e-10;
  
  for (std::size_t i = 0; i < num_points; ++i)
  {
    ceres::CostFunction *cost_function = homography_reproj_error::create(src_points[i], dst_points[i]);
    problem->AddResidualBlock(cost_function, NULL, homography.data());
  }
  ceres::Solve(options, problem.get(), &summary);
  
  
  for(std::size_t i = 0; i < 3; ++i)
  {
    for(std::size_t j = 0; j < 3; ++j)
    {
      H.at<double>(i, j) = homography[3 * i + j];
    }
  }
}

haeryong

이전 포스트

[Learn OpenGL] Part1. Getting started 메모

다음 포스트