1 #include <boost/foreach.hpp> 6 #include <boost/numeric/ublas/matrix.hpp> 11 #include <Eigen/Cholesky> 13 #include <pcl/filters/filter.h> 14 #include <pcl/point_types.h> 25 #define MEAS_DIM_RBF 3 28 #define COV_MEAS_RBF 1.0 31 _estimation_error_threshold(-1.),
33 _velocity( 0., 0., 0., 0., 0., 0.),
36 _min_num_supporting_feats_to_correct(5),
37 _use_predicted_state_from_kh(false),
38 _use_predicted_measurement_from_kh(false),
39 _last_time_interval_ns(0.0)
41 _pose = Eigen::Matrix4d::Identity();
46 const std::vector<Eigen::Matrix4d>& trajectory,
47 const Eigen::Twistd &initial_velocity,
49 double estimation_error_threshold) :
54 _pose(trajectory[trajectory.size()-1]),
62 this->
_trajectory = std::vector<Eigen::Matrix4d>(trajectory.begin()+1, trajectory.end());
94 if(time_interval_ns < 0)
96 std::cout <<
"Time interval negative!" << std::endl;
102 Eigen::Matrix<double, 6, 6> sys_noise_COV = Eigen::Matrix<double, 6, 6>::Zero();
116 Eigen::Twistd belief_delta_pose_ec = (time_interval_ns/1e9)*
_velocity;
120 Eigen::Matrix<double, 6, 6> delta_pose_covariance = (time_interval_ns/1e9)*
_velocity_covariance;
122 Eigen::Twistd pose_ec;
124 Eigen::Matrix<double,6,6> adjoint = pose_ec.exp(1e-12).adjoint();
199 int meas_dimension = 3 * num_supporting_feats;
201 double feature_depth = 0.;
203 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> > G_t((this->
_G_t_memory.data()), 6, meas_dimension);
204 Eigen::Map<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> > R_inv_times_innovation((this->
_R_inv_times_innovation_memory.data()), meas_dimension, 1);
208 R_inv_times_innovation.setZero();
214 Eigen::Matrix3d R_seq = Eigen::Matrix3d::Zero();
218 Eigen::Matrix<double, 6, 3> P_HT_seq = Eigen::Matrix<double, 6, 3>::Zero();
219 Eigen::Matrix3d S_seq = Eigen::Matrix3d::Zero();
220 Eigen::Matrix3d S_inv_seq = Eigen::Matrix3d::Zero();
221 Eigen::Matrix<double, 6, 3> K_seq = Eigen::Matrix<double, 6, 3>::Zero();
223 Eigen::Matrix<double, 6, 6> temp_seq1 = Eigen::Matrix<double, 6, 6>::Zero();
224 Eigen::Matrix<double, 6, 6> temp_seq2 = Eigen::Matrix<double, 6, 6>::Zero();
226 Eigen::Vector4d feature_innovation;
227 Eigen::Vector4d location_at_birth_eig;
229 Eigen::Vector4d predicted_location;
230 Eigen::Vector4d current_location;
238 if(supporting_feat_ptr->getFeatureAge() > 1)
247 feature_innovation = current_location - predicted_location;
249 feature_depth = supporting_feat_ptr->getLastZ();
259 Eigen::Vector4d predicted_location_feat =
_predicted_pose*location_at_birth_eig;
278 K_seq = P_HT_seq * S_inv_seq;
282 temp_seq2 = temp_seq1 * P_minus_seq;
284 P_minus_seq -= temp_seq2;
287 for (
unsigned int i=0; i<P_minus_seq.rows(); i++ )
289 for (
unsigned int j=0; j<=i; j++ )
291 P_minus_seq(j,i) = P_minus_seq(i,j);
300 Eigen::Matrix<double, 6, 1> pose_update = P_minus_seq * G_t * R_inv_times_innovation;
302 Eigen::Twistd pose_update_ec(pose_update[3], pose_update[4], pose_update[5], pose_update[0], pose_update[1], pose_update[2]);
312 Eigen::Twistd delta_pose_ec;
318 bool isPredictionYBetterThanX(boost::tuple<
size_t,
double, Eigen::Matrix<double, 6, 6> > X, boost::tuple<
size_t,
double, Eigen::Matrix<double, 6, 6> > Y)
322 if(boost::tuples::get<0>(X) < boost::tuples::get<0>(Y))
328 else if(boost::tuples::get<0>(X) > boost::tuples::get<0>(Y))
337 if(boost::tuples::get<1>(X) > boost::tuples::get<1>(Y))
343 else if(boost::tuples::get<1>(X) < boost::tuples::get<1>(Y))
351 int num_larger_x_covariance = 0;
352 for(
int i =0; i<6; i++)
354 if(boost::tuples::get<2>(X)(i,i) > boost::tuples::get<2>(Y)(i,i))
356 num_larger_x_covariance +=1;
359 if(num_larger_x_covariance > 3)
376 std::vector<Feature::Id> alive_supporting_features_ids;
381 alive_supporting_features_ids.push_back(supporting_feat_id);
386 std::vector<boost::tuple<size_t, double, Eigen::Matrix<double, 6, 6> > > num_supporting_feats_and_error_and_pose_cov;
390 std::vector<Feature::Id> supporting_feats_based_on_vh;
391 std::vector<Feature::Id> supporting_feats_based_on_bh;
392 std::vector<Feature::Id> supporting_feats_based_on_kh;
393 double error_vh = 0.0;
394 double error_bh = 0.0;
395 double error_kh = 0.0;
396 double acc_error_vh = 0.0;
397 double acc_error_bh = 0.0;
398 double acc_error_kh = 0.0;
404 error_vh =
L2Distance(last_location, predicted_location_vh);
407 supporting_feats_based_on_vh.push_back(supporting_feat_id);
408 acc_error_vh += error_vh;
412 error_bh =
L2Distance(last_location, predicted_location_bh);
413 if (error_bh < this->_estimation_error_threshold)
415 supporting_feats_based_on_bh.push_back(supporting_feat_id);
416 acc_error_bh += error_bh;
423 error_kh =
L2Distance(last_location, predicted_location_kh);
424 if (error_kh < this->_estimation_error_threshold)
426 supporting_feats_based_on_kh.push_back(supporting_feat_id);
427 acc_error_kh += error_kh;
432 num_supporting_feats_and_error_and_pose_cov.push_back(boost::tuple<
size_t,
double, Eigen::Matrix<double, 6, 6> >(supporting_feats_based_on_vh.size(), acc_error_vh,
_predicted_pose_cov_vh));
433 num_supporting_feats_and_error_and_pose_cov.push_back(boost::tuple<
size_t,
double, Eigen::Matrix<double, 6, 6> >(supporting_feats_based_on_bh.size(), acc_error_bh,
_predicted_pose_cov_bh));
438 num_supporting_feats_and_error_and_pose_cov.push_back(boost::tuple<
size_t,
double, Eigen::Matrix<double, 6, 6> >(supporting_feats_based_on_kh.size(), acc_error_kh,
_predicted_pose_cov_kh));
443 std::vector<boost::tuple<size_t, double, Eigen::Matrix<double, 6, 6> > >::iterator hyp_with_max_support_it =
444 std::max_element(num_supporting_feats_and_error_and_pose_cov.begin(),
445 num_supporting_feats_and_error_and_pose_cov.end(),
447 int hyp_with_max_support_idx = std::distance(num_supporting_feats_and_error_and_pose_cov.begin(), hyp_with_max_support_it);
449 switch(hyp_with_max_support_idx)
454 std::cout <<
"RB" <<
_id<<
" Using velocity prediction" << std::endl;
465 std::cout <<
"RB" <<
_id<<
" Using brake prediction" << std::endl;
477 "RB" <<
_id <<
" Using prediction from higher level!");
499 Eigen::Twistd predicted_pose_kh_ec(hypothesis.pose_wc.twist.angular.x,
500 hypothesis.pose_wc.twist.angular.y,
501 hypothesis.pose_wc.twist.angular.z,
502 hypothesis.pose_wc.twist.linear.x,
503 hypothesis.pose_wc.twist.linear.y,
504 hypothesis.pose_wc.twist.linear.z);
508 hypothesis.velocity_wc.twist.angular.y,
509 hypothesis.velocity_wc.twist.angular.z,
510 hypothesis.velocity_wc.twist.linear.x,
511 hypothesis.velocity_wc.twist.linear.y,
512 hypothesis.velocity_wc.twist.linear.z);
513 for(
int i =0; i<6; i++)
515 for(
int j=0; j<6; j++)
525 Eigen::Twistd shape_based_pose_ec = Eigen::Twistd(twist_refinement.twist.angular.x, twist_refinement.twist.angular.y, twist_refinement.twist.angular.z,
526 twist_refinement.twist.linear.x, twist_refinement.twist.linear.y, twist_refinement.twist.linear.z);
530 ROS_INFO_STREAM_NAMED(
"RBFilter.integrateShapeBasedPose",
"Using shape-based pose (time elapsed=" << pose_time_elapsed_ns/1e9 <<
534 this->
_pose = shape_based_pose_ec.exp(1e-12).toHomogeneousMatrix();
541 Eigen::Twistd delta_pose_ec;
545 ROS_INFO_STREAM_NAMED(
"RBFilter.integrateShapeBasedPose",
"Shape-based pose is too old (time elapsed=" << pose_time_elapsed_ns/1e9 <<
546 " s, loop_period="<<
_loop_period_ns/1e9 <<
" s). Try reducing the framerate");
575 const Eigen::Matrix4d& predicted_pose,
576 bool contemporary_last_feat_loc_and_last_pose_in_trajectory,
578 bool publish_locations)
const 586 if(publish_locations)
592 at_birth_feat->x = boost::tuples::get<0>(at_birth_location);
593 at_birth_feat->y = boost::tuples::get<1>(at_birth_location);
594 at_birth_feat->z = boost::tuples::get<2>(at_birth_location);
595 at_birth_feat->label = feature->getId();
600 at_birth_feat->x = boost::tuples::get<0>(location_at_birth_of_rb);
601 at_birth_feat->y = boost::tuples::get<1>(location_at_birth_of_rb);
602 at_birth_feat->z = boost::tuples::get<2>(location_at_birth_of_rb);
603 at_birth_feat->label = feature->getId();
608 predicted_feat->x = boost::tuples::get<0>(predicted_location);
609 predicted_feat->y = boost::tuples::get<1>(predicted_location);
610 predicted_feat->z = boost::tuples::get<2>(predicted_location);
611 predicted_feat->label = feature->getId();
617 bool contemporary_last_feat_loc_and_last_pose_in_trajectory,
623 size_t feature_age = feature->getFeatureAge();
629 if(!contemporary_last_feat_loc_and_last_pose_in_trajectory)
640 if (feature_age < rigid_body_age + 1)
642 Eigen::Matrix4d transform_from_rb_birth_to_feat_birth = this->
_trajectory.at(rigid_body_age - feature_age);
645 Eigen::Matrix4d inverse_transform_from_rb_birth_to_feat_birth = transform_from_rb_birth_to_feat_birth.inverse();
646 TransformLocation(feat_traj.at(0), inverse_transform_from_rb_birth_to_feat_birth, location_at_birth);
650 location_at_birth = feat_traj.at(feature_age - (rigid_body_age + 1));
687 geometry_msgs::PoseWithCovariancePtr pose_wc_ptr = geometry_msgs::PoseWithCovariancePtr(
new geometry_msgs::PoseWithCovariance());
688 pose_wc_ptr->pose.position.x =
_pose(0,3);
689 pose_wc_ptr->pose.position.y =
_pose(1,3);
690 pose_wc_ptr->pose.position.z =
_pose(2,3);
691 Eigen::Quaterniond quaternion(
_pose.block<3,3>(0,0));
692 pose_wc_ptr->pose.orientation.x = quaternion.x();
693 pose_wc_ptr->pose.orientation.y = quaternion.y();
694 pose_wc_ptr->pose.orientation.z = quaternion.z();
695 pose_wc_ptr->pose.orientation.w = quaternion.w();
697 for (
unsigned int i = 0; i < 6; i++)
698 for (
unsigned int j = 0; j < 6; j++)
706 geometry_msgs::TwistWithCovariancePtr pose_wc_ptr = geometry_msgs::TwistWithCovariancePtr(
new geometry_msgs::TwistWithCovariance());
707 Eigen::Twistd pose_ec;
710 pose_wc_ptr->twist.linear.x = pose_ec.vx();
711 pose_wc_ptr->twist.linear.y = pose_ec.vy();
712 pose_wc_ptr->twist.linear.z = pose_ec.vz();
713 pose_wc_ptr->twist.angular.x = pose_ec.rx();
714 pose_wc_ptr->twist.angular.y = pose_ec.ry();
715 pose_wc_ptr->twist.angular.z = pose_ec.rz();
717 for (
unsigned int i = 0; i < 6; i++)
718 for (
unsigned int j = 0; j < 6; j++)
728 geometry_msgs::TwistWithCovariancePtr vel_t_wc_ptr = geometry_msgs::TwistWithCovariancePtr(
new geometry_msgs::TwistWithCovariance());
729 vel_t_wc_ptr->twist.linear.x =
_velocity.vx();
730 vel_t_wc_ptr->twist.linear.y =
_velocity.vy();
731 vel_t_wc_ptr->twist.linear.z =
_velocity.vz();
732 vel_t_wc_ptr->twist.angular.x = velocity_copy.rx();
733 vel_t_wc_ptr->twist.angular.y = velocity_copy.ry();
734 vel_t_wc_ptr->twist.angular.z = velocity_copy.rz();
736 for (
unsigned int i = 0; i < 6; i++)
737 for (
unsigned int j = 0; j < 6; j++)
763 return predicted_measurement;
785 for (
int i = 0; i < 6; i++)
FeatureCloudPCLwc::Ptr _predicted_measurement_pc_bh
virtual void setPredictedState(omip_msgs::RigidBodyPoseAndVelMsg hypothesis)
virtual void estimateBestPredictionAndSupportingFeatures()
std::map< Feature::Id, Feature::Location > _predicted_measurement_map_vh
boost::tuple< double, double, double > Location
FeatureCloudPCLwc::Ptr _predicted_measurement
double _last_time_interval_ns
virtual std::vector< Feature::Id > getSupportingFeatures() const
virtual std::vector< Eigen::Matrix4d > getTrajectory() const
virtual void addSupportingFeature(Feature::Id supporting_feat_id)
virtual Eigen::Matrix< double, 6, 6 > getPoseCovariance() const
#define ROS_ERROR_STREAM_NAMED(name, args)
virtual geometry_msgs::PoseWithCovariancePtr getPoseWithCovariance() const
virtual void setFeaturesDatabase(FeaturesDataBase::Ptr feats_database)
virtual FeatureCloudPCLwc getPredictedMeasurement(PredictionHypothesis hypothesis=BASED_ON_VELOCITY)
FeatureCloudPCLwc::Ptr _predicted_measurement_pc_vh
virtual void correctState()
bool _use_predicted_state_from_kh
virtual geometry_msgs::TwistWithCovariancePtr getPoseECWithCovariance() const
double L2Distance(const Feature::Location &first, const Feature::Location &second)
std::map< Feature::Id, Feature::Location > _predicted_measurement_map
void TransformMatrix2Twist(const Eigen::Matrix4d &transformation_matrix, Eigen::Twistd &twist)
Eigen::Matrix4d _predicted_pose
Eigen::VectorXd _R_inv_times_innovation_memory
Eigen::Matrix< double, 6, 6 > _predicted_velocity_cov_kh
#define ROS_INFO_STREAM_NAMED(name, args)
FeatureCloudPCLwc::Ptr _predicted_measurement_pc_kh
virtual Eigen::Matrix4d getPose() const
FeatureCloudPCLwc::Ptr _features_at_birth
Feature::Location _initial_location_of_centroid
Eigen::Twistd _predicted_velocity_bh
virtual void integrateShapeBasedPose(const geometry_msgs::TwistWithCovariance twist_refinement, double pose_time_elapsed_ns)
static RB_id_t _rb_id_generator
std::map< Feature::Id, Feature::Location > _predicted_measurement_map_bh
virtual void _initializeAuxiliarMatrices()
FeatureCloudPCLwc::Ptr _features_predicted
std::vector< Feature::Id > _supporting_features_ids
TFSIMD_FORCE_INLINE Quaternion inverse(const Quaternion &q)
Eigen::Matrix< double, 3, 6 > _D_T_p0_circle
void invert3x3MatrixEigen2(const Eigen::Matrix3d &to_inv, Eigen::Matrix3d &inverse)
virtual void GetLocationOfFeatureAtBirthOfRB(Feature::Ptr feature, bool contemporary_last_feat_loc_and_last_pose_in_trajectory, Feature::Location &location_at_birth) const
Get the location of the feature at the frame when the RB was detected first. If the feature is "young...
Eigen::Matrix4d _predicted_pose_kh
bool _use_predicted_measurement_from_kh
Eigen::Matrix< double, 6, 6 > _predicted_velocity_cov_bh
std::vector< Location > Trajectory
virtual geometry_msgs::TwistWithCovariancePtr getVelocityWithCovariance() const
Eigen::Matrix< double, 6, 6 > _predicted_pose_covariance
Eigen::Matrix< double, 6, 6 > _pose_covariance
FeaturesDataBase::Ptr _features_database
void LocationOfFeature2EigenVectorHomogeneous(const Feature::Location &lof, Eigen::Vector4d &eig_vec)
virtual FeaturesDataBase::Ptr getFeaturesDatabase() const
virtual Eigen::Matrix< double, 6, 6 > getVelocityCovariance() const
Eigen::Twistd _predicted_velocity_kh
double _estimation_error_threshold
Eigen::Matrix< double, 6, 6 > _predicted_pose_cov_bh
virtual RB_id_t getId() const
virtual void predictState(double time_interval_ns=-1.)
int _min_num_supporting_feats_to_correct
bool isPredictionYBetterThanX(boost::tuple< size_t, double, Eigen::Matrix< double, 6, 6 > > X, boost::tuple< size_t, double, Eigen::Matrix< double, 6, 6 > > Y)
double _meas_depth_factor
std::vector< Eigen::Matrix4d > _trajectory
virtual Eigen::Twistd getVelocity() const
pcl::PointCloud< FeaturePCLwc > FeatureCloudPCLwc
void TransformLocation(const Feature::Location &origin, const Eigen::Matrix4d &transformation, Feature::Location &new_location)
void LocationAndId2FeaturePCLwc(const Feature::Location &feature_location, const Feature::Id &feature_id, omip::FeaturePCLwc &feature_pcl)
Eigen::Matrix4d _predicted_pose_bh
virtual int getNumberSupportingFeatures() const
std::map< Feature::Id, Feature::Location > _predicted_measurement_map_kh
virtual void addPredictedFeatureLocation(const Feature::Location &predicted_location_velocity, const Feature::Location &predicted_location_brake, const Feature::Id feature_id)
Eigen::Matrix4d _predicted_pose_vh
Eigen::MatrixXd _G_t_memory
Eigen::Matrix< double, 6, 6 > _predicted_pose_cov_kh
Eigen::Matrix< double, 6, 6 > _velocity_covariance
virtual void predictMeasurement()
Second step when updating the filter. The next measurement is predicted from the predicted next state...
virtual void PredictFeatureLocation(Feature::Ptr feature, const Eigen::Matrix4d &predicted_pose, bool contemporary_last_feat_loc_and_last_pose_in_trajectory, Feature::Location &predicted_location, bool publish_locations=false) const
Eigen::Matrix< double, 6, 6 > _predicted_pose_cov_vh