SmartStereoProjectionFactorPP.h

Go to the documentation of this file.

/* ----------------------------------------------------------------------------
 
 * GTSAM Copyright 2010, Georgia Tech Research Corporation,
 * Atlanta, Georgia 30332-0415
 * All Rights Reserved
 * Authors: Frank Dellaert, et al. (see THANKS for the full author list)
 
 * See LICENSE for the license information
 
 * -------------------------------------------------------------------------- */

 
#pragma once
 
#include <gtsam_unstable/slam/SmartStereoProjectionFactor.h>
 
namespace gtsam {

class GTSAM_UNSTABLE_EXPORT SmartStereoProjectionFactorPP
    : public SmartStereoProjectionFactor {
 protected:
  std::vector<std::shared_ptr<Cal3_S2Stereo>> K_all_;

  KeyVector world_P_body_keys_;

  KeyVector body_P_cam_keys_;
 
 public:
  typedef SmartStereoProjectionFactor Base;

  typedef SmartStereoProjectionFactorPP This;

  typedef std::shared_ptr<This> shared_ptr;
 
  static const int DimBlock = 12;  
  static const int DimPose = 6;  
  static const int ZDim = 3;  
  typedef Eigen::Matrix<double, ZDim, DimBlock> MatrixZD;  // F blocks (derivatives wrt camera)
  typedef std::vector<MatrixZD, Eigen::aligned_allocator<MatrixZD> > FBlocks;  // vector of F blocks

  SmartStereoProjectionFactorPP(const SharedNoiseModel& sharedNoiseModel,
                                const SmartStereoProjectionParams& params =
                                    SmartStereoProjectionParams());

  void add(const StereoPoint2& measured, const Key& world_P_body_key,
           const Key& body_P_cam_key,
           const std::shared_ptr<Cal3_S2Stereo>& K);

  void add(const std::vector<StereoPoint2>& measurements,
           const KeyVector& w_P_body_keys, const KeyVector& body_P_cam_keys,
           const std::vector<std::shared_ptr<Cal3_S2Stereo>>& Ks);

  void add(const std::vector<StereoPoint2>& measurements,
           const KeyVector& w_P_body_keys, const KeyVector& body_P_cam_keys,
           const std::shared_ptr<Cal3_S2Stereo>& K);

  void print(const std::string& s = "", const KeyFormatter& keyFormatter =
                 DefaultKeyFormatter) const override;

  bool equals(const NonlinearFactor& p, double tol = 1e-9) const override;

  const KeyVector& getExtrinsicPoseKeys() const {
    return body_P_cam_keys_;
  }

  double error(const Values& values) const override;

  inline std::vector<std::shared_ptr<Cal3_S2Stereo>> calibration() const {
    return K_all_;
  }

  Base::Cameras cameras(const Values& values) const override;

  void computeJacobiansAndCorrectForMissingMeasurements(
      FBlocks& Fs, Matrix& E, Vector& b, const Values& values) const {
    if (!result_) {
      throw("computeJacobiansWithTriangulatedPoint");
    } else {  // valid result: compute jacobians
      size_t numViews = measured_.size();
      E = Matrix::Zero(3 * numViews, 3); // a StereoPoint2 for each view (point jacobian)
      b = Vector::Zero(3 * numViews);  // a StereoPoint2 for each view
      Matrix dPoseCam_dPoseBody_i, dPoseCam_dPoseExt_i, dProject_dPoseCam_i, Ei;
 
      for (size_t i = 0; i < numViews; i++) {  // for each camera/measurement
        Pose3 w_P_body = values.at<Pose3>(world_P_body_keys_.at(i));
        Pose3 body_P_cam = values.at<Pose3>(body_P_cam_keys_.at(i));
        StereoCamera camera(
            w_P_body.compose(body_P_cam, dPoseCam_dPoseBody_i, dPoseCam_dPoseExt_i),
            K_all_[i]);
        // get jacobians and error vector for current measurement
        StereoPoint2 reprojectionError_i = StereoPoint2(
            camera.project(*result_, dProject_dPoseCam_i, Ei) - measured_.at(i));
        Eigen::Matrix<double, ZDim, DimBlock> J;  // 3 x 12
        J.block<ZDim, 6>(0, 0) = dProject_dPoseCam_i * dPoseCam_dPoseBody_i;  // (3x6) * (6x6)
        J.block<ZDim, 6>(0, 6) = dProject_dPoseCam_i * dPoseCam_dPoseExt_i;  // (3x6) * (6x6)
        // if the right pixel is invalid, fix jacobians
        if (std::isnan(measured_.at(i).uR()))
        {
          J.block<1, 12>(1, 0) = Matrix::Zero(1, 12);
          Ei.block<1, 3>(1, 0) = Matrix::Zero(1, 3);
          reprojectionError_i = StereoPoint2(reprojectionError_i.uL(), 0.0,
                                           reprojectionError_i.v());
        }
        // fit into the output structures
        Fs.push_back(J);
        size_t row = 3 * i;
        b.segment<ZDim>(row) = -reprojectionError_i.vector();
        E.block<3, 3>(row, 0) = Ei;
      }
    }
  }

  std::shared_ptr<RegularHessianFactor<DimPose>> createHessianFactor(
      const Values& values, const double lambda = 0.0,
      bool diagonalDamping = false) const {
    // we may have multiple cameras sharing the same extrinsic cals, hence the number
    // of keys may be smaller than 2 * nrMeasurements (which is the upper bound where we
    // have a body key and an extrinsic calibration key for each measurement)
    size_t nrUniqueKeys = keys_.size();
 
    // Create structures for Hessian Factors
    KeyVector js;
    std::vector<Matrix> Gs(nrUniqueKeys * (nrUniqueKeys + 1) / 2);
    std::vector<Vector> gs(nrUniqueKeys);
 
    if (this->measured_.size() != cameras(values).size())
      throw std::runtime_error("SmartStereoProjectionHessianFactor: this->"
          "measured_.size() inconsistent with input");
 
    // triangulate 3D point at given linearization point
    triangulateSafe(cameras(values));
 
    // failed: return "empty/zero" Hessian
    if (!result_) {
      for (Matrix& m : Gs) m = Matrix::Zero(DimPose, DimPose);
      for (Vector& v : gs) v = Vector::Zero(DimPose);
      return std::make_shared<RegularHessianFactor<DimPose>>(keys_, Gs, gs,
                                                             0.0);
    }
 
    // compute Jacobian given triangulated 3D Point
    FBlocks Fs;
    Matrix F, E;
    Vector b;
    computeJacobiansAndCorrectForMissingMeasurements(Fs, E, b, values);
 
    // Whiten using noise model
    noiseModel_->WhitenSystem(E, b);
    for (size_t i = 0; i < Fs.size(); i++) {
      Fs[i] = noiseModel_->Whiten(Fs[i]);
    }
 
    // build augmented Hessian (with last row/column being the information vector)
    Matrix3 P;
    Cameras::ComputePointCovariance<3>(P, E, lambda, diagonalDamping);
 
    // these are the keys that correspond to the blocks in augmentedHessian (output of SchurComplement)
    KeyVector nonuniqueKeys;
    for (size_t i = 0; i < world_P_body_keys_.size(); i++) {
      nonuniqueKeys.push_back(world_P_body_keys_.at(i));
      nonuniqueKeys.push_back(body_P_cam_keys_.at(i));
    }
    // but we need to get the augumented hessian wrt the unique keys in key_
    SymmetricBlockMatrix augmentedHessianUniqueKeys =
        Base::Cameras::template SchurComplementAndRearrangeBlocks<3, DimBlock,
                                                                  DimPose>(
            Fs, E, P, b, nonuniqueKeys, keys_);
 
    return std::make_shared<RegularHessianFactor<DimPose>>(
        keys_, augmentedHessianUniqueKeys);
  }

  std::shared_ptr<GaussianFactor> linearizeDamped(
      const Values& values, const double lambda = 0.0) const {
    // depending on flag set on construction we may linearize to different linear factors
    switch (params_.linearizationMode) {
      case HESSIAN:
        return createHessianFactor(values, lambda);
      default:
        throw std::runtime_error(
            "SmartStereoProjectionFactorPP: unknown linearization mode");
    }
  }

  std::shared_ptr<GaussianFactor> linearize(const Values& values) const
      override {
    return linearizeDamped(values);
  }
 
 private:
#if GTSAM_ENABLE_BOOST_SERIALIZATION  
  friend class boost::serialization::access;
  template<class ARCHIVE>
  void serialize(ARCHIVE& ar, const unsigned int /*version*/) {
    ar& BOOST_SERIALIZATION_BASE_OBJECT_NVP(Base);
    ar & BOOST_SERIALIZATION_NVP(K_all_);
  }
#endif
};
// end of class declaration

template<>
struct traits<SmartStereoProjectionFactorPP> : public Testable<
    SmartStereoProjectionFactorPP> {
};
 
}  // namespace gtsam