一文詳解bundle adjustment

2021-01-10 3D視覺工坊

來源：公眾號|3D視覺工坊（系投稿）作者：李城「3D視覺工坊」技術交流群已經成立，目前大約有12000人，方向主要涉及3D視覺、CV&深度學習、SLAM、三維重建、點雲後處理、自動駕駛、CV入門、三維測量、VR/AR、3D人臉識別、醫療影像、缺陷檢測、行人重識別、目標跟蹤、視覺產品落地、視覺競賽、車牌識別、硬體選型、學術交流、求職交流、ORB-SLAM系列源碼交流、深度估計等。QQ群「3D視覺研習社」，群號：574432628。

bundle adjustment 的歷史發展

bundle adjustment,中文名稱是光束法平差,經典的BA目的是優化相機的pose和landmark,其在SfM和SLAM 領域中扮演者重要角色.目前大多數書籍或者參老文獻將其翻譯成"捆綁調整"是不太嚴謹的做法.bundle adjustment 最早是19世紀由搞大地測量學(測繪學科)的人提出來的,19世紀中期的時候，geodetics的學者就開始研究large scale triangulations（大型三角剖分）。20世紀中期，隨著camera和computer的出現，photogrammetry(攝影測量學)也開始研究adjustment computation，所以他們給起了個名字叫bundle adjustment(隸屬攝影測量學科前輩的功勞)。21世紀前後，robotics領域開始興起SLAM，最早用的recursive bayesian filter（遞歸貝葉斯濾波），後來把問題搞成個graph然後用least squares方法求解,bundle adjusment歷史發展圖如下:

bundle adjustment 其本質還是離不開最小二乘原理(Gauss功勞)(幾乎所有優化問題其本質都是最小二乘),目前bundle adjustment 優化框架最為代表的是ceres solver和g2o(這裡主要介紹ceres solver).據說ceres的命名是天文學家Piazzi閒暇無事的時候觀測一顆沒有觀測到的星星,最後用least squares算出了這個小行星的軌道,故將這個小行星命名為ceres.

Bundle adjustment 的算法理論

觀測值:像點坐標 優化量(平差量):pose 和landmark 因為一旦涉及平差,就必定有如下公式:觀測值+觀測值改正數=近似值+近似值改正數,那麼bundle adjustment 的公式還是從共線條件方程出發:

四種Bundle adjustment 算法代碼

這裡代碼主要從四個方面來介紹:

優化相機內參及畸變係數,相機的pose(6dof)和landmark 代價函數寫法如下:template <typename CameraModel>class BundleAdjustmentCostFunction { public: explicit BundleAdjustmentCostFunction(const Eigen::Vector2d& point2D) : observed_x_(point2D(0)), observed_y_(point2D(1)) {}//構造函數傳入的是觀測值 static ceres::CostFunction* Create(const Eigen::Vector2d& point2D) { return (new ceres::AutoDiffCostFunction< BundleAdjustmentCostFunction<CameraModel>, 2, 4, 3, 3, CameraModel::kNumParams>( new BundleAdjustmentCostFunction(point2D))); }//優化量:2代表誤差方程個數;4代表pose中的姿態信息,用四元數表示;3代表pose中的位置信息;3代表landmark自由度;CameraModel::kNumParams是相機內參和畸變係數,其取決於相機模型是what// opertator 重載函數的參數即是待優化的量 template <typename T> bool operator()(const T* const qvec, const T* const tvec, const T* const point3D, const T* const camera_params, T* residuals) const { // Rotate and translate. T projection[3]; ceres::UnitQuaternionRotatePoint(qvec, point3D, projection); projection[0] += tvec[0]; projection[1] += tvec[1]; projection[2] += tvec[2]; // Project to image plane. projection[0] /= projection[2]; projection[1] /= projection[2]; // Distort and transform to pixel space. CameraModel::WorldToImage(camera_params, projection[0], projection[1], &residuals[0], &residuals[1]); // Re-projection error. residuals[0] -= T(observed_x_); residuals[1] -= T(observed_y_); return true; } private: const double observed_x_; const double observed_y_;};寫好了代價函數,下面就是需要把參數都加入殘差塊,讓ceres自動求導,代碼如下:

ceres::Problem problem;ceres::CostFunction* cost_function = nullptr; ceres::LossFunction * p_LossFunction = ceres_options_.bUse_loss_function_ ? new ceres::HuberLoss(Square(4.0)) : nullptr; // 關於為何使用損失函數,因為現實中並不是所有觀測過程中的噪聲都服從 //gaussian noise的（或者可以說幾乎沒有）， //遇到有outlier的情況，這些方法非常容易掛掉， //這時候就得用到robust statistics裡面的 //robust cost(*cost也可以叫做loss, 統計學那邊喜歡叫risk) function了， //比較常用的有huber, cauchy等等。cost_function = BundleAdjustmentCostFunction<CameraModel>::Create(point2D.XY()); //將優化量加入殘差塊problem_->AddResidualBlock(cost_function, p_LossFunction, qvec_data, tvec_data, point3D.XYZ().data(), camera_params_data); 如上,case1 的bundle adjustment 就搭建完成!

優化相機內參及畸變係數,pose subset parameterization(pose 信息部分參數優化)和3D landmark,當 只優化姿態信息時候,problem需要添加的代碼如下: //這裡姿態又用歐拉角表示 map_poses[indexPose] = {angleAxis[0], angleAxis[1], angleAxis[2], t(0), t(1), t(2)}; double * parameter_block = &map_poses.at(indexPose)[0]; problem.AddParameterBlock(parameter_block, 6); std::vector<int> vec_constant_extrinsic; vec_constant_extrinsic.insert(vec_constant_extrinsic.end(), {3,4,5}); if (!vec_constant_extrinsic.empty()) { // 主要用到ceres的SubsetParameterization函數 ceres::SubsetParameterization *subset_parameterization = new ceres::SubsetParameterization(6, vec_constant_extrinsic); problem.SetParameterization(parameter_block, subset_parameterization); } 優化相機內參及畸變係數,pose subset parameterization(pose 信息部分參數優化)和3D landmark,當 只優化位置信息時候,problem需要添加的代碼如下(同上面代碼,只需修改一行): //這裡姿態又用歐拉角表示 map_poses[indexPose] = {angleAxis[0], angleAxis[1], angleAxis[2], t(0), t(1), t(2)}; double * parameter_block = &map_poses.at(indexPose)[0]; problem.AddParameterBlock(parameter_block, 6); std::vector<int> vec_constant_extrinsic; vec_constant_extrinsic.insert(vec_constant_extrinsic.end(), {1,2,3}); if (!vec_constant_extrinsic.empty()) { ceres::SubsetParameterization *subset_parameterization = new ceres::SubsetParameterization(6, vec_constant_extrinsic); problem.SetParameterization(parameter_block, subset_parameterization); } 優化相機內參及畸變係數,pose 是常量不優化 和3D landmark. 代價函數寫法如下://相機模型template <typename CameraModel> class BundleAdjustmentConstantPoseCostFunction { public: BundleAdjustmentConstantPoseCostFunction(const Eigen::Vector4d& qvec, const Eigen::Vector3d& tvec, const Eigen::Vector2d& point2D) : qw_(qvec(0)), qx_(qvec(1)), qy_(qvec(2)), qz_(qvec(3)), tx_(tvec(0)), ty_(tvec(1)), tz_(tvec(2)), observed_x_(point2D(0)), observed_y_(point2D(1)) {} static ceres::CostFunction* Create(const Eigen::Vector4d& qvec, const Eigen::Vector3d& tvec, const Eigen::Vector2d& point2D) { return (new ceres::AutoDiffCostFunction< BundleAdjustmentConstantPoseCostFunction<CameraModel>, 2, 3, CameraModel::kNumParams>( new BundleAdjustmentConstantPoseCostFunction(qvec, tvec, point2D))); } template <typename T> bool operator()(const T* const point3D, const T* const camera_params, T* residuals) const { const T qvec[4] = {T(qw_), T(qx_), T(qy_), T(qz_)}; // Rotate and translate. T projection[3]; ceres::UnitQuaternionRotatePoint(qvec, point3D, projection); projection[0] += T(tx_); projection[1] += T(ty_); projection[2] += T(tz_); // Project to image plane. projection[0] /= projection[2]; projection[1] /= projection[2]; // Distort and transform to pixel space. CameraModel::WorldToImage(camera_params, projection[0], projection[1], &residuals[0], &residuals[1]); // Re-projection error. residuals[0] -= T(observed_x_); residuals[1] -= T(observed_y_); return true; } private: const double qw_; const double qx_; const double qy_; const double qz_; const double tx_; const double ty_; const double tz_; const double observed_x_; const double observed_y_;};接下來problem 加入殘差塊代碼如下:

ceres::Problem problem;ceres::CostFunction* cost_function = nullptr; ceres::LossFunction * p_LossFunction = ceres_options_.bUse_loss_function_ ? new ceres::HuberLoss(Square(4.0)) : nullptr; // 關於為何使用損失函數,因為現實中並不是所有觀測過程中的噪聲都服從 //gaussian noise的（或者可以說幾乎沒有）， //遇到有outlier的情況，這些方法非常容易掛掉， //這時候就得用到robust statistics裡面的 //robust cost(*cost也可以叫做loss, 統計學那邊喜歡叫risk) function了， //比較常用的有huber, cauchy等等。cost_function = BundleAdjustmentConstantPoseCostFunction<CameraModel>::Create( \ image.Qvec(), image.Tvec(), point2D.XY());//觀測值輸入 //將優化量加入殘差塊 problem_->AddResidualBlock(cost_function, loss_function, \ point3D.XYZ().data(), camera_params_data);//被優化量加入殘差-3D點和相機內參 以上就是四種BA 的case 當然還可以有很多變種,比如gps約束的BA(即是附有限制條件的間接平差),比如 固定3D landmark,優化pose和相機參數和畸變係數

參考資料

colmap openmvg 原始碼,github 地址:https://github.com/openMVG/openMVGhttps://github.com/colmap/colmap單傑. 光束法平差簡史與概要. 武漢大學學報·信息科學版, 2018, 43(12): 1797-1810.備註：作者也是我們「3D視覺從入門到精通」特邀嘉賓：一個超乾貨的3D視覺學習社區

本文僅做學術分享，如有侵權，請聯繫刪文。

相關焦點

  • 玩轉iPhone:iPhone各種文件路徑詳解
    本文導航1iPhone各種文件路徑詳解（一）2iPhone各種文件路徑詳解（二）3iPhone各種文件路徑詳解（三）4iPhone各種文件路徑詳解（四）<<返回分頁閱讀1iPhone各種文件路徑詳解（一）　　【天極網手機頻道
  • 一文詳解EMI的傳播過程
    打開APP 一文詳解EMI的傳播過程 電子工程師筆記 發表於 2020-12-31 15:17:14 　　電磁幹擾是電子電路設計過程中最常見的問題
  • 基因組denovo組裝:一文詳解
    當我們拿到下機原始數據，並進行預處理後【二、三、四代測序下機數據質控詳解】，就可以開始進行基因組denovo組裝了。
  • 一文詳解HDMI製造過程
    打開APP 一文詳解HDMI製造過程 線纜行業朋友圈 發表於 2021-01-09 11:17:48 隨著HDMI2.1認證的開啟
  • 一文詳解負溫度係數熱敏電阻
    打開APP 一文詳解負溫度係數熱敏電阻 工程師a 發表於 2018-05-20 07:03:00
  • 一文詳解射頻信號源的工作原理
    一文詳解射頻信號源的工作原理 測試那些事兒 發表於 2020-11-19 09:44:27 射頻信號源顧名思義就是產線射頻信號的一個源，或者說是一臺儀表。
  • 一文詳解EV SSL證書,EV、OV和DV證書有何區別?
    一文詳解EV SSL證書，EV、OV和DV證書有何區別？ EV SSL證書是SSL數字證書認證等級裡審核最嚴格、安全性最高的證書。
  • 一文詳解 Word2vec 之 Skip-Gram 模型(訓練篇)
    對優化目標採用「negative sampling」方法，這樣每個訓練樣本的訓練只會更新一小部分的模型權重，從而降低計算負擔。事實證明，對常用詞抽樣並且對優化目標採用「negative sampling」不僅降低了訓練過程中的計算負擔，還提高了訓練的詞向量的質量。
  • 一文讀懂電容傳感器
    藍色標題，獲取文章】 10、一文讀懂光纖傳感器 11、一文讀懂溫溼度傳感器 12
  • 【泡泡一分鐘】挑戰性光照條件下的視覺裡程計多模態跟蹤框架
    每天一分鐘，帶你讀遍機器人頂級會議文章標題：Multimodal tracking framework for visual odometry in challenging illumination conditions作者：Axel Beauvisage, Kenan Ahiska, Nabil Aouf來源：2020 IEEE International
  • 一文讀懂MEMS傳感器(必須收藏)
    【點擊藍色標題，獲取文章】 1
  • SIGIR 2020開幕在即,一文詳解過去五年引用量TOP論文
    AMiner 上線了 SIGIR 2020 最新專題，收錄了今年錄用的所有論文，並對過去歷年來 SIGIR 的錄用論文數據進行了分析，後期還將陸續推出論文數據分析與論文精讀等內容，為大家詳解 SIGIR 2020 最新動態。
  • 高中化學:必修一~必修二知識點詳解!學霸得高分的秘訣
    但是從必修一開始就對這門學科放棄，是絕對不可以的，因為這樣是會影響到之後三年的。因為必修一、必修二就是高中化學的基礎知識，只有學好這一部分的知識才能學好後面的內容，所以說每個想要在高考取得好的成績的，都應該要把這一部分的知識就搞懂搞通。
  • 2020屆深圳市高三年級二模英語試題,答案詳解
    【3題詳解】細節理解題。根據Staten Island Skating Pavilion部分中This large area maintains its year-round frosty temperatures for ice-skating fun.可知，這一大片地區保持了全年的霜凍溫度，為滑冰提供樂趣。
  • 一文詳解紅外遙控模塊工作原理
    本文首先介紹了紅外遙控模塊的基本原理，其次詳解闡述了紅外遙控模塊工作原理，最後介紹了紅外遙控的重要環節及應用。
  • [Surviving Sepsis Campaign]: 拯救膿毒症指南更新一小時集束化治療
    The Surviving Sepsis Campaign has updated the Hour-1 bundle emphasizing time of recognition and element completion
  • 一文讀懂磁傳感器(必須收藏)
    【點擊藍色標題，獲取文章】 一文讀懂溫溼度傳感器 12、一文讀懂圖像傳感器 13、一文讀懂生物傳感器
  • 一文詳解MSMS
    批量生產MEMS採用類似集成電路(IC)的生產工藝和加工過程，用矽微加工工藝在一矽片上可同時製造成百上千個微型機電裝置或完整的MEMS。使MEMS有極高的自動化程度，批量生產可大大降低生產成本;而且地球表層矽的含量為2%。幾乎取之不盡，因此MEMS產品在經濟性方面更具競爭力。
  • 2020年高考英語全國卷一完形填空答案詳解
    第41題詳解：考查名詞詞義辨析。根據下文「If we admit it is a door, they』ll want to go outside constantly . It will drive us crazy.（如果我們承認那是一扇門，他們就會不斷地想出去，而這會讓我們發瘋的。）」可知這正是我們把門解釋成窗戶的原因，所以可知此處需要的意思是「原因」。