OpenMVGをVisual Studio 2017(x64)でビルド(その3)

その2から随分と間があいてしまいましたが、今回その3として、サンプルソースベースにOpenMVGでSfM(Structure from Motion)の概要について追いかけみたいと思います。

その2までの投稿を参考にVisual Studio 2017でOpenMVGのプロジェクトを開いてください。

ソリューションエクスプローラで[openMVG_sample_multiview_robustEssential]を右クリックして[スタートアッププロジェクトに設定]を選択します。同じくソリューションエクスプローラで[openMVG_sample_multiview_robustEssential]のツリーを開いて、さらにその下の[Source Files]ツリーを開いて[robust_essential.cpp]をダブルクリックして、ソースコードを開きます。

それでは、ざっとソースコードを覗いてみます。

  const std::string sInputDir = stlplus::folder_up(string(THIS_SOURCE_DIR))
    + "/imageData/SceauxCastle/";
  const string jpg_filenameL = sInputDir + "100_7101.jpg";
  const string jpg_filenameR = sInputDir + "100_7102.jpg";

  Image<unsigned char> imageL, imageR;
  ReadImage(jpg_filenameL.c_str(), &imageL);
  ReadImage(jpg_filenameR.c_str(), &imageR);

 

「src\openMVG_Samples\imageData\SceauxCastle」フォルダにある2つの画像ファイル100_7101.jpgと100_7102.jpgを開いています。

これら2枚の画像をSfMにより解析して三次元点群データを生成します。大きな流れとしては以下のようになります。

1. 特徴点を抽出
2. 特徴点のマッチングを実行し、特徴点ペアを生成
3. Robust Estimationで間違ったペアを排除をし、カメラ行列(Relative Pose)を取得
4. 算出したカメラ行列を用いて特徴点ペアの三次元座標を算出し、点群データを生成

 それでは上記ポイントをソースコードで追いかけます。

1.特徴点を抽出

  using namespace openMVG::features;
  std::unique_ptr<Image_describer> image_describer(new SIFT_Anatomy_Image_describer);
  std::map<IndexT, std::unique_ptr<features::Regions>> regions_perImage;
  image_describer->Describe(imageL, regions_perImage[0]);
  image_describer->Describe(imageR, regions_perImage[1]);

  const SIFT_Regions* regionsL = dynamic_cast<SIFT_Regions*>(regions_perImage.at(0).get());
  const SIFT_Regions* regionsR = dynamic_cast<SIFT_Regions*>(regions_perImage.at(1).get());

  const PointFeatures
    featsL = regions_perImage.at(0)->GetRegionsPositions(),
    featsR = regions_perImage.at(1)->GetRegionsPositions();

 

抽出した結果は build\openMVG_Samples\multiview_robust_essential\02_features.svg に出力されています。このファイルをダブルクリックしてInternetExplorer等で確認出来ます。

 

 

2.特徴点のマッチングを実行し、特徴点ペアを生成

  std::vector<IndMatch> vec_PutativeMatches;
  //-- Perform matching -> find Nearest neighbor, filtered with Distance ratio
  {
    // Find corresponding points
    matching::DistanceRatioMatch(
      0.8, matching::BRUTE_FORCE_L2,
      *regions_perImage.at(0).get(),
      *regions_perImage.at(1).get(),
      vec_PutativeMatches);

    IndMatchDecorator<float> matchDeduplicator(
            vec_PutativeMatches, featsL, featsR);
    matchDeduplicator.getDeduplicated(vec_PutativeMatches);

 

以下が求められたペアですが、同じ城を指している点で一部上下の写真で斜めになっている線があると思います。これらが間違ったペアです。

 

3.Robust Estimationで間違ったペアを排除をし、カメラ行列を取得

    Mat3 K;
    //read K from file
    if (!readIntrinsic(stlplus::create_filespec(sInputDir,"K","txt"), K))
    {
      std::cerr << "Cannot read intrinsic parameters." << std::endl;
      return EXIT_FAILURE;
    }

Calibration行列Kをファイルから読み込みます。

1452.44 0 708
0 1452.44 532
0 0 1

CCDカメラのCalibration行列はMultiple View Geometry in Computer Visionの(6.9)式によると以下のようになっています。

 axはayはピクセル単位に変換したカメラの焦点距離、x0とy0はPrinciple Point(カメラの中心座標)です。

    const Pinhole_Intrinsic
      camL(imageL.Width(), imageL.Height(), K(0,0), K(0,2), K(1,2)),
      camR(imageR.Width(), imageR.Height(), K(0,0), K(0,2), K(1,2));

カメラ内部行列を生成します。

    //A. prepare the corresponding putatives points
    Mat xL(2, vec_PutativeMatches.size());
    Mat xR(2, vec_PutativeMatches.size());
    for (size_t k = 0; k < vec_PutativeMatches.size(); ++k)  {
      const PointFeature & imaL = featsL[vec_PutativeMatches[k].i_];
      const PointFeature & imaR = featsR[vec_PutativeMatches[k].j_];
      xL.col(k) = imaL.coords().cast<double>();
      xR.col(k) = imaR.coords().cast<double>();
    }

    //B. Compute the relative pose thanks to a essential matrix estimation
    const std::pair<size_t, size_t> size_imaL(imageL.Width(), imageL.Height());
    const std::pair<size_t, size_t> size_imaR(imageR.Width(), imageR.Height());
    sfm::RelativePose_Info relativePose_info;

    if (!sfm::robustRelativePose(&camL, &camR, xL, xR, relativePose_info, size_imaL, size_imaR, 256))

    {
      std::cerr << " /!\\ Robust relative pose estimation failure."
        << std::endl;
      return EXIT_FAILURE;
    }

robustRelativePoseでRobust Estimationを実行します。算出されたカメラ行列はrelativePose_infoに入ります。

下図がRobust Estimationの結果で残ったペア、つまりInlierです。先の結果で間違ったペアが削除されているのが分かります。

 

 

4.算出したカメラ行列を用いて特徴点ペアの三次元座標を算出し、点群データを生成

    //C. Triangulate and export inliers as a PLY scene
    std::vector<Vec3> vec_3DPoints;

    // Setup camera intrinsic and poses
    Pinhole_Intrinsic intrinsic0(imageL.Width(), imageL.Height(), K(0, 0), K(0, 2), K(1, 2));
    Pinhole_Intrinsic intrinsic1(imageR.Width(), imageR.Height(), K(0, 0), K(0, 2), K(1, 2));

    const Pose3 pose0 = Pose3(Mat3::Identity(), Vec3::Zero());
    const Pose3 pose1 = relativePose_info.relativePose;

    // Init structure by inlier triangulation

    const Mat34 P1 = intrinsic0.get_projective_equivalent(pose0);
    const Mat34 P2 = intrinsic1.get_projective_equivalent(pose1);
Robust Estimationの結果からカメラ行列を生成します。

 

    std::vector<double> vec_residuals;
    vec_residuals.reserve(relativePose_info.vec_inliers.size() * 4);
    for (size_t i = 0; i < relativePose_info.vec_inliers.size(); ++i)  {
      const SIOPointFeature & LL = regionsL->Features()[vec_PutativeMatches[relativePose_info.vec_inliers[i]].i_];
      const SIOPointFeature & RR = regionsR->Features()[vec_PutativeMatches[relativePose_info.vec_inliers[i]].j_];
      // Point triangulation
      Vec3 X;
      TriangulateDLT(
        P1, LL.coords().cast<double>().homogeneous(),
        P2, RR.coords().cast<double>().homogeneous(), &X);
      // Reject point that is behind the camera

TriangulateDLTは、Multiple View Geometry in Computer Visionの12.2章に記載されているLinear triangulation methodsという2点の二次元座標から三次元座標を算出する解法です。

 

      if (pose0.depth(X) < 0 && pose1.depth(X) < 0)
        continue;

      const Vec2 residual0 = intrinsic0.residual(pose0(X), LL.coords().cast<double>());
      const Vec2 residual1 = intrinsic1.residual(pose1(X), RR.coords().cast<double>());
      vec_residuals.push_back(std::abs(residual0(0)));
      vec_residuals.push_back(std::abs(residual0(1)));
      vec_residuals.push_back(std::abs(residual1(0)));
      vec_residuals.push_back(std::abs(residual1(1)));
      vec_3DPoints.emplace_back(X);
    }

 

 求めた三次元座標はEssentialGeometry.plyとして保存されています。少し分かりにくいですが、下図がその保存された三次元座標です。