From 874bc35d2208b4aad289e20fe8958cfcbb67f1bf Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sat, 10 Oct 2026 17:02:13 -0700 Subject: [PATCH 1/7] Imu motion predictor --- CMakeLists.txt | 2 +- .../include/rtabmap/core/ImuMotionPredictor.h | 187 +++++++++ corelib/include/rtabmap/core/Odometry.h | 3 + corelib/include/rtabmap/core/Parameters.h | 4 +- corelib/src/CMakeLists.txt | 1 + corelib/src/ImuMotionPredictor.cpp | 389 ++++++++++++++++++ corelib/src/Odometry.cpp | 46 +++ corelib/test/CMakeLists.txt | 1 + corelib/test/test_imumotionpredictor.cpp | 354 ++++++++++++++++ guilib/src/PreferencesDialog.cpp | 2 + guilib/src/ui/preferencesDialog.ui | 80 +++- package.xml | 2 +- 12 files changed, 1054 insertions(+), 17 deletions(-) create mode 100644 corelib/include/rtabmap/core/ImuMotionPredictor.h create mode 100644 corelib/src/ImuMotionPredictor.cpp create mode 100644 corelib/test/test_imumotionpredictor.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7992c03a..177f8448 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -22,7 +22,7 @@ SET(CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake_modules") ####################### SET(RTABMAP_MAJOR_VERSION 0) SET(RTABMAP_MINOR_VERSION 24) -SET(RTABMAP_PATCH_VERSION 0) +SET(RTABMAP_PATCH_VERSION 1) SET(RTABMAP_VERSION ${RTABMAP_MAJOR_VERSION}.${RTABMAP_MINOR_VERSION}.${RTABMAP_PATCH_VERSION}) diff --git a/corelib/include/rtabmap/core/ImuMotionPredictor.h b/corelib/include/rtabmap/core/ImuMotionPredictor.h new file mode 100644 index 00000000..4b5e9d19 --- /dev/null +++ b/corelib/include/rtabmap/core/ImuMotionPredictor.h @@ -0,0 +1,187 @@ +/* +Copyright (c) 2010-2026, Mathieu Labbe - IntRoLab - Universite de Sherbrooke +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the Universite de Sherbrooke nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#ifndef RTABMAP_CORE_IMUMOTIONPREDICTOR_H_ +#define RTABMAP_CORE_IMUMOTIONPREDICTOR_H_ + +#include +#include +#include + +#include + +#include + +namespace rtabmap { + +/** + * @brief Predicts the pose of the base frame from the last odometry pose and the IMU. + * + * Between two odometry updates, the pose is propagated with the IMU: the orientation is + * the IMU's own, re-expressed in the odometry frame, and the position is integrated from + * the velocity at the last odometry update and the gravity-compensated acceleration. That + * velocity is re-estimated at every odometry update, from the displacement since an + * odometry pose a short window back (see velocityWindow) corrected by the acceleration + * measured in between, so the position never drifts for long: only the motion since the + * last update is predicted. + * + * This is what lidar deskewing needs: the pose at every point's time, during a sweep that + * started after the last pose odometry estimated. + * + * The IMU samples are given as they are measured (see addImu()): + * the orientation of the IMU in its world frame and its specific force, from which gravity + * is removed here. That world frame must be gravity aligned with +z up, as in ROS (REP-103, + * e.g. ENU); its yaw doesn't matter. An orientation given in a frame with z down (NED) + * must be converted first, otherwise gravity is added instead of removed. The lever arm between the IMU and the base origin is + * ignored: its centripetal and tangential accelerations are small over the fraction of a + * second this predicts. + * + * Not thread-safe: a caller sharing it between threads must lock around every call. + */ +class RTABMAP_CORE_EXPORT ImuMotionPredictor +{ +public: + /** + * Without acceleration in the IMU samples, the position follows the last velocity + * (constant velocity model). + * + * @param maxPoseInterval odometry poses older (s) than this are not used to estimate + * the velocity, which is null without one + * @param velocityWindow the velocity is estimated from the displacement since the + * newest pose at least this old (s). Over a single frame + * interval, the noise of the odometry poses would be of the order + * of the velocity itself; with the acceleration, a longer window + * still gives the velocity at the last pose, not an average. + * @param gravity magnitude (m/s^2) of the gravity removed from the specific force + * given to addImu(), standard gravity by default + * + * These are fixed for the life of the predictor (there is no setter): changing them + * while it estimates would mix poses and samples taken under different settings. + */ + explicit ImuMotionPredictor( + double maxPoseInterval = 1.0, + double velocityWindow = 0.5, + double gravity = 9.80665); + + double maxPoseInterval() const {return maxPoseInterval_;} + double velocityWindow() const {return velocityWindow_;} + double gravity() const {return gravity_;} + + /** + * @brief Adds an IMU measurement. + * @param stamp time of the measurement (s) + * @param imu orientation of the IMU in its world frame (gravity aligned, +z up), linear + * acceleration as measured (the specific force, which includes the + * reaction to gravity) and the transform from the base frame to the IMU. + * Without orientation, the measurement is ignored; without linear + * acceleration (all zeros, or a covariance of -1), only its orientation + * is used. + */ + void addImu(double stamp, const IMU & imu); + + /** + * @brief Adds an odometry pose, from which the next poses are predicted. + * @param stamp time of the pose (s) + * @param pose pose of the base frame in the odometry frame; a null pose (odometry + * lost) resets the prediction until the next valid pose + */ + void addPose(double stamp, const rtabmap::Transform & pose); + + /// Forgets the poses and the IMU samples. + void reset(); + + /** + * @brief Predicts the pose of the base frame in the odometry frame. + * @param stamp time (s) of the prediction, normally after the last pose + * @return the predicted pose, null if there is no IMU sample yet. Without a pose (none + * yet, or the last one was null), the orientation alone is predicted, in the + * IMU's world frame and at the origin: still the relative rotation between two + * stamps, which is what deskewing needs most. + */ + rtabmap::Transform predict(double stamp) const; + + /// Whether there is a pose to predict from (none yet, or the last one was null). + bool hasPose() const; + /// Stamp of the last pose added, 0 if there is none. + double lastPoseStamp() const; + /// Velocity (m/s) of the base in the odometry frame at the last pose. + Eigen::Vector3d velocity() const; + /// Number of IMU samples kept. + size_t samples() const; + +private: + // orientation: of the base frame in the IMU's world frame; acceleration: of the base in + // that frame, gravity removed + void addSample(double stamp, const Eigen::Quaterniond & orientation, const Eigen::Vector3d & acceleration); + + struct Sample + { + Eigen::Quaterniond orientation; + Eigen::Vector3d acceleration; + }; + Eigen::Quaterniond orientationAt(double stamp) const; + Eigen::Vector3d accelerationAt(double stamp) const; + void integrate(double from, double to, const Eigen::Quaterniond & rotation, + Eigen::Vector3d & velocity, Eigen::Vector3d & position) const; + void updateIntegration() const; + + // The acceleration integrated from the last pose up to a sample's stamp, in the + // odometry frame: predicting a stamp then only integrates from the sample before it. + struct Integrated + { + Eigen::Vector3d acceleration; // at that stamp + Eigen::Vector3d velocity; // change since the last pose + Eigen::Vector3d position; // change since the last pose (without its velocity) + }; + +private: + double maxPoseInterval_; + double velocityWindow_; + double gravity_; + std::map samples_; + // Recent odometry poses, to estimate the velocity from (see velocityWindow) + std::map poses_; + + // The last odometry pose and the state predictions start from. + double poseStamp_; + rtabmap::Transform pose_; + Eigen::Vector3d velocity_; + // Rotation from the IMU's world frame to the odometry frame, at the last pose: the two + // are both gravity aligned, but their yaw differ. + Eigen::Quaterniond worldToOdom_; + + // Built lazily by predict(), from the last pose to the newest sample; cleared when the + // pose or the samples it was built from change. + mutable std::map integrated_; + // Whether the acceleration at the last pose's stamp (the first entry) is final: it is + // held constant from the newest sample until a sample after that stamp is received. + mutable bool integratedStartIsFinal_; +}; + +} + +#endif /* RTABMAP_CORE_IMUMOTIONPREDICTOR_H_ */ diff --git a/corelib/include/rtabmap/core/Odometry.h b/corelib/include/rtabmap/core/Odometry.h index 21f94b05..524e589a 100644 --- a/corelib/include/rtabmap/core/Odometry.h +++ b/corelib/include/rtabmap/core/Odometry.h @@ -33,6 +33,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include #include +#include namespace rtabmap { @@ -159,6 +160,7 @@ private: bool _force3DoF; bool _holonomic; bool guessFromMotion_; + bool guessImuAcceleration_; float guessSmoothingDelay_; int _filteringStrategy; int _particleSize; @@ -189,6 +191,7 @@ private: std::vector stereoModels_; std::vector models_; std::map imus_; + ImuMotionPredictor imuMotionPredictor_; // used with Odom/GuessImuAcceleration only }; } /* namespace rtabmap */ diff --git a/corelib/include/rtabmap/core/Parameters.h b/corelib/include/rtabmap/core/Parameters.h index 8d076f71..07a854a6 100644 --- a/corelib/include/rtabmap/core/Parameters.h +++ b/corelib/include/rtabmap/core/Parameters.h @@ -512,7 +512,9 @@ class RTABMAP_CORE_EXPORT Parameters RTABMAP_PARAM(Odom, KalmanProcessNoise, float, 0.001, "Process noise covariance value."); RTABMAP_PARAM(Odom, KalmanMeasurementNoise, float, 0.01, "Process measurement covariance value."); RTABMAP_PARAM(Odom, GuessMotion, bool, true, "Guess next transformation from the last motion computed."); - RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate.", kOdomFilteringStrategy().c_str())); + RTABMAP_PARAM(Odom, GuessImuAcceleration, bool, false, uFormat("With an IMU giving orientation and linear acceleration, predict the translation of the motion guess (with \"%s\") by integrating the IMU acceleration from the velocity at the previous frame, instead of with a constant velocity. That velocity is the average over \"%s\" (from the displacement over that delay), carried to the previous frame with the acceleration measured since: it is not delayed like the average alone. Set \"%s\" to about 0.5 s: over a single frame, the noise of the poses is of the order of the velocity itself. Ignored by odometry approaches processing the IMU themselves.", kOdomGuessMotion().c_str(), kOdomGuessSmoothingDelay().c_str(), kOdomGuessSmoothingDelay().c_str())); + RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration with \"%s\". Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's.", kOdomGuessImuAcceleration().c_str())); + RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate. Recommended (~0.5 s) with \"%s\": the translational velocity is then the displacement over this delay corrected by the IMU acceleration, so that it is the velocity at the last frame rather than a delayed average. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing (\"%s\"), where its noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomGuessImuAcceleration().c_str(), kOdomDeskewing().c_str())); RTABMAP_PARAM(Odom, KeyFrameThr, float, 0.3, "[Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, VisKeyFrameThr, int, 150, "[Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, ScanKeyFrameThr, float, 0.9, "[Geometry] Create a new keyframe when the number of ICP inliers drops under this ratio of points in last frame's scan. Setting the value to 0 means that a keyframe is created for each processed frame."); diff --git a/corelib/src/CMakeLists.txt b/corelib/src/CMakeLists.txt index 892279c9..d693707f 100644 --- a/corelib/src/CMakeLists.txt +++ b/corelib/src/CMakeLists.txt @@ -88,6 +88,7 @@ SET(SRC_FILES RegistrationVis.cpp Odometry.cpp + ImuMotionPredictor.cpp OdometryThread.cpp OdometryInfo.cpp odometry/OdometryF2M.cpp diff --git a/corelib/src/ImuMotionPredictor.cpp b/corelib/src/ImuMotionPredictor.cpp new file mode 100644 index 00000000..f00a1959 --- /dev/null +++ b/corelib/src/ImuMotionPredictor.cpp @@ -0,0 +1,389 @@ +/* +Copyright (c) 2010-2026, Mathieu Labbe - IntRoLab - Universite de Sherbrooke +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the Universite de Sherbrooke nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +#include + +#include + +namespace rtabmap { + +// Samples older than this (s) behind the newest one are dropped, so the buffer stays +// bounded while no odometry pose comes to trim it. +static const double kMaxBufferDuration = 10.0; + +ImuMotionPredictor::ImuMotionPredictor(double maxPoseInterval, double velocityWindow, double gravity) : + maxPoseInterval_(maxPoseInterval), + velocityWindow_(velocityWindow), + gravity_(gravity), + integratedStartIsFinal_(false), + poseStamp_(0.0), + velocity_(Eigen::Vector3d::Zero()), + worldToOdom_(Eigen::Quaterniond::Identity()) +{ +} + +void ImuMotionPredictor::addImu(double stamp, const IMU & imu) +{ + const cv::Vec4d & o = imu.orientation(); + const Eigen::Quaterniond imuOrientation(o[3], o[0], o[1], o[2]); + if(imu.empty() || + imuOrientation.norm() < 0.5 || + (!imu.orientationCovariance().empty() && imu.orientationCovariance().at(0,0) == -1.0)) + { + // No orientation + return; + } + const Eigen::Quaterniond worldToImu = imuOrientation.normalized(); + const Eigen::Quaterniond baseToImu = imu.localTransform().isNull()? + Eigen::Quaterniond::Identity(): + imu.localTransform().getQuaterniond(); + + const cv::Vec3d & f = imu.linearAcceleration(); + Eigen::Vector3d acceleration = Eigen::Vector3d::Zero(); + if((f[0] != 0.0 || f[1] != 0.0 || f[2] != 0.0) && + (imu.linearAccelerationCovariance().empty() || imu.linearAccelerationCovariance().at(0,0) != -1.0)) + { + // The specific force includes the reaction to gravity, up in the world frame + acceleration = worldToImu * Eigen::Vector3d(f[0], f[1], f[2]) - Eigen::Vector3d(0, 0, gravity_); + } + addSample(stamp, worldToImu * baseToImu.inverse(), acceleration); +} + +void ImuMotionPredictor::addSample(double stamp, const Eigen::Quaterniond & orientation, const Eigen::Vector3d & acceleration) +{ + Sample sample; + sample.orientation = orientation.normalized(); + sample.acceleration = acceleration; + samples_[stamp] = sample; + if(!integrated_.empty() && stamp <= integrated_.rbegin()->first) + { + // Out of order: what was integrated after it changes + integrated_.clear(); + } + while(samples_.size() > 2 && samples_.begin()->first < stamp - kMaxBufferDuration) + { + samples_.erase(samples_.begin()); + integrated_.clear(); + } +} + +void ImuMotionPredictor::addPose(double stamp, const rtabmap::Transform & pose) +{ + integrated_.clear(); + if(pose.isNull()) + { + pose_.setNull(); + poseStamp_ = 0.0; + velocity_.setZero(); + poses_.clear(); + return; + } + + Eigen::Vector3d velocity = Eigen::Vector3d::Zero(); + Eigen::Quaterniond worldToOdom = Eigen::Quaterniond::Identity(); + if(!samples_.empty()) + { + // The odometry and the IMU agree on the orientation of the base at that time, + // whatever the yaw of their frames. + worldToOdom = (pose.getQuaterniond() * orientationAt(stamp).inverse()).normalized(); + } + + // The velocity is estimated from the displacement since a previous pose. Not the + // last one: odometry's own noise, divided by a frame interval, would be of the order of + // the velocity itself, and since the prediction deskews the next scans, that noise would + // feed back into the next poses. The newest pose at least velocityWindow_ old is used + // instead (or the oldest kept if there is none yet), not older than maxPoseInterval_. + poses_.erase(poses_.lower_bound(stamp), poses_.end()); + poses_.erase(poses_.begin(), poses_.lower_bound(stamp - maxPoseInterval_)); + std::map::const_iterator reference = poses_.begin(); + for(std::map::const_iterator iter = poses_.begin(); + iter != poses_.end() && iter->first <= stamp - velocityWindow_; ++iter) + { + reference = iter; + } + if(reference != poses_.end()) + { + const double interval = stamp - reference->first; + const Eigen::Vector3d displacement( + pose.x() - reference->second.x(), + pose.y() - reference->second.y(), + pose.z() - reference->second.z()); + + if(!samples_.empty() && + samples_.begin()->first <= reference->first && + samples_.rbegin()->first >= stamp) + { + // displacement = v0*T + D, with D the double integral of the acceleration + // over the interval: solve for v0, the velocity at the reference pose, then + // carry it to this pose with the single integral V. + Eigen::Vector3d deltaVelocity; + Eigen::Vector3d deltaPosition; + integrate(reference->first, stamp, worldToOdom, deltaVelocity, deltaPosition); + velocity = (displacement - deltaPosition) / interval + deltaVelocity; + } + else + { + // Average velocity over the interval (the samples don't cover it) + velocity = displacement / interval; + } + } + + pose_ = pose; + poseStamp_ = stamp; + velocity_ = velocity; + worldToOdom_ = worldToOdom; + poses_[stamp] = pose; + + // The next velocities are estimated from the poses kept: keep the samples from the + // oldest one, including the last sample before it to interpolate at its stamp. + std::map::iterator iter = samples_.upper_bound(poses_.begin()->first); + if(iter != samples_.begin()) + { + --iter; + samples_.erase(samples_.begin(), iter); + } +} + +void ImuMotionPredictor::reset() +{ + samples_.clear(); + poses_.clear(); + integrated_.clear(); + pose_.setNull(); + poseStamp_ = 0.0; + velocity_.setZero(); + worldToOdom_.setIdentity(); +} + +rtabmap::Transform ImuMotionPredictor::predict(double stamp) const +{ + if(samples_.empty()) + { + return rtabmap::Transform(); + } + if(pose_.isNull()) + { + // No pose yet (or lost): only the orientation is known, in the IMU's world frame. + const Eigen::Quaterniond orientation = orientationAt(stamp); + return rtabmap::Transform(0, 0, 0, orientation.x(), orientation.y(), orientation.z(), orientation.w()); + } + + const Eigen::Quaterniond orientation = (worldToOdom_ * orientationAt(stamp)).normalized(); + + // With t0 the stamp of the last pose, p0 its position and v0 the velocity there, and + // a(u) the acceleration (gravity removed, in the odometry frame), the position at t is: + // + // p(t) = p0 + v0*(t-t0) + D(t), with D(t) = integral_t0^t integral_t0^s a(u) du ds + // + // D(t) is the displacement due to the change of velocity since t0, V(s) = integral_t0^s a(u) du. + Eigen::Vector3d position(pose_.x(), pose_.y(), pose_.z()); // p0 + position += velocity_ * (stamp - poseStamp_); // + v0*(t-t0) + if(stamp >= poseStamp_) + { + // + D(t). D and V are kept at every sample stamp since t0 (integrated_), so only + // the part from the last sample ti <= t is integrated here, with dt = t-ti: + // + // D(t) = D(ti) + V(ti)*dt + integral_ti^t integral_ti^s a(u) du ds + // + // Between samples the acceleration is linear, from a(ti) to a(t), for which that + // last double integral is exactly dt^2*(2*a(ti) + a(t))/6. Those are the same + // segments integrate() would go through, without redoing all those before ti. + updateIntegration(); + std::map::const_iterator from = integrated_.upper_bound(stamp); + --from; // ti: the first entry is at t0, so there is one + const double dt = stamp - from->first; + const Eigen::Vector3d accelerationB = worldToOdom_ * accelerationAt(stamp); // a(t) + position += from->second.position + // D(ti) + from->second.velocity * dt + // V(ti)*dt + dt * dt * (2.0 * from->second.acceleration + accelerationB) / 6.0; // a(ti) to a(t) + } + else + { + // + D(t), backward from t0: only for points stamped before the last pose, rare + Eigen::Vector3d deltaVelocity; + Eigen::Vector3d deltaPosition; + integrate(poseStamp_, stamp, worldToOdom_, deltaVelocity, deltaPosition); + position += deltaPosition; + } + + return rtabmap::Transform(position.x(), position.y(), position.z(), + orientation.x(), orientation.y(), orientation.z(), orientation.w()); +} + +bool ImuMotionPredictor::hasPose() const +{ + return !pose_.isNull(); +} + +double ImuMotionPredictor::lastPoseStamp() const +{ + return poseStamp_; +} + +Eigen::Vector3d ImuMotionPredictor::velocity() const +{ + return velocity_; +} + +size_t ImuMotionPredictor::samples() const +{ + return samples_.size(); +} + +// samples_ must not be empty for the four below. + +void ImuMotionPredictor::updateIntegration() const +{ + if(!integrated_.empty() && !integratedStartIsFinal_ && samples_.rbegin()->first >= poseStamp_) + { + // The acceleration at the pose was held from the newest sample, which is no + // longer the newest + integrated_.clear(); + } + if(integrated_.empty()) + { + Integrated start; + start.acceleration = worldToOdom_ * accelerationAt(poseStamp_); + start.velocity.setZero(); + start.position.setZero(); + integrated_[poseStamp_] = start; + integratedStartIsFinal_ = samples_.rbegin()->first >= poseStamp_; + } + // Extend to the samples received since + for(std::map::const_iterator iter = samples_.upper_bound(integrated_.rbegin()->first); + iter != samples_.end(); ++iter) + { + const std::pair & previous = *integrated_.rbegin(); + const double dt = iter->first - previous.first; + Integrated next; + // Same segment as in integrate(), from the previous sample ti to this one: + // D(ti+1) = D(ti) + V(ti)*dt + dt^2*(2*a(ti) + a(ti+1))/6, V(ti+1) = V(ti) + dt*(a(ti) + a(ti+1))/2 + next.acceleration = worldToOdom_ * iter->second.acceleration; + next.position = previous.second.position + previous.second.velocity * dt + + dt * dt * (2.0 * previous.second.acceleration + next.acceleration) / 6.0; + next.velocity = previous.second.velocity + dt * (previous.second.acceleration + next.acceleration) / 2.0; + integrated_.insert(integrated_.end(), std::make_pair(iter->first, next)); + } +} + +Eigen::Quaterniond ImuMotionPredictor::orientationAt(double stamp) const +{ + std::map::const_iterator after = samples_.lower_bound(stamp); + if(after == samples_.end()) + { + return samples_.rbegin()->second.orientation; + } + if(after == samples_.begin() || after->first == stamp) + { + return after->second.orientation; + } + std::map::const_iterator before = std::prev(after); + const double ratio = (stamp - before->first) / (after->first - before->first); + return before->second.orientation.slerp(ratio, after->second.orientation); +} + +Eigen::Vector3d ImuMotionPredictor::accelerationAt(double stamp) const +{ + std::map::const_iterator after = samples_.lower_bound(stamp); + if(after == samples_.end()) + { + return samples_.rbegin()->second.acceleration; + } + if(after == samples_.begin() || after->first == stamp) + { + return after->second.acceleration; + } + std::map::const_iterator before = std::prev(after); + const double ratio = (stamp - before->first) / (after->first - before->first); + return before->second.acceleration + ratio * (after->second.acceleration - before->second.acceleration); +} + +void ImuMotionPredictor::integrate(double from, double to, const Eigen::Quaterniond & rotation, + Eigen::Vector3d & velocity, Eigen::Vector3d & position) const +{ + velocity.setZero(); + position.setZero(); + if(from == to) + { + return; + } + + // Breakpoints: the bounds and every sample in between, in the direction of + // integration (backward if "to" is before "from"). Between two of them the + // acceleration is linear, which the segment update below integrates exactly; before + // the first sample and after the last one, it is held constant. + std::vector stamps; + stamps.push_back(from); + if(from < to) + { + for(std::map::const_iterator iter = samples_.upper_bound(from); + iter != samples_.end() && iter->first < to; ++iter) + { + stamps.push_back(iter->first); + } + } + else + { + std::map::const_iterator iter = samples_.lower_bound(from); + while(iter != samples_.begin()) + { + --iter; + if(iter->first <= to) + { + break; + } + stamps.push_back(iter->first); + } + } + stamps.push_back(to); + + // Returned: with a(u) the acceleration rotated by "rotation", from "from" (t0) to "to" (t), + // + // velocity = V(t) = integral_t0^t a(u) du + // position = D(t) = integral_t0^t integral_t0^s a(u) du ds + // + // that is, the change of velocity and the displacement it causes (from a velocity null + // at t0: the caller adds v0*(t-t0)). They are accumulated breakpoint by breakpoint: from + // ti to the next one ti+1, with dt = ti+1 - ti and a(u) linear from a(ti) to a(ti+1), + // + // D(ti+1) = D(ti) + V(ti)*dt + dt^2*(2*a(ti) + a(ti+1))/6 + // V(ti+1) = V(ti) + dt*(a(ti) + a(ti+1))/2 + // + // both exact for a linear acceleration (dt is negative backward, the same formulas hold). + Eigen::Vector3d accelerationA = rotation * accelerationAt(stamps[0]); // a(t0) + for(size_t i=1; i + +#include +#include +#include + +#include +#include + +using rtabmap::ImuMotionPredictor; + +namespace { + +Eigen::Quaterniond yaw(double angle) +{ + return Eigen::Quaterniond(Eigen::AngleAxisd(angle, Eigen::Vector3d::UnitZ())); +} + +rtabmap::Transform pose(const Eigen::Vector3d & position, double angle) +{ + return rtabmap::Transform(position.x(), position.y(), position.z(), 0, 0, angle); +} + +/// An IMU at rest or accelerating, as measured: orientation of the IMU in its world +/// frame, and the specific force (acceleration minus gravity) in the IMU frame. +rtabmap::IMU measuredImu(const Eigen::Quaterniond & worldToImu, + const Eigen::Vector3d & accelerationInWorld, double gravity, + const rtabmap::Transform & baseToImu) +{ + const Eigen::Vector3d f = worldToImu.inverse() * (accelerationInWorld + Eigen::Vector3d(0, 0, gravity)); + const Eigen::Quaterniond q = worldToImu.normalized(); + return rtabmap::IMU(cv::Vec4d(q.x(), q.y(), q.z(), q.w()), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(0,0,0), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(f.x(), f.y(), f.z()), cv::Mat::eye(3,3,CV_64FC1), + baseToImu); +} + +/// IMU measurements at 200 Hz over [from, to] of an IMU at the base origin, oriented and +/// accelerating (in its world frame) as given. Without accelerometer, it measures no +/// linear acceleration at all. +void addSamples(ImuMotionPredictor & predictor, double from, double to, + const std::function & orientation, + const std::function & acceleration, + bool withAccelerometer = true) +{ + for(int i=0; from + i*0.005 <= to + 1e-9; ++i) + { + const double t = from + i*0.005; + rtabmap::IMU imu = measuredImu(orientation(t), acceleration(t), predictor.gravity(), rtabmap::Transform::getIdentity()); + if(!withAccelerometer) + { + imu = rtabmap::IMU(imu.orientation(), imu.orientationCovariance(), + imu.angularVelocity(), imu.angularVelocityCovariance(), + cv::Vec3d(0,0,0), imu.linearAccelerationCovariance(), + imu.localTransform()); + } + predictor.addImu(t, imu); + } +} + +} // namespace + +TEST(ImuMotionPredictor, predicts_only_the_orientation_without_a_pose) +{ + ImuMotionPredictor predictor; + EXPECT_TRUE(predictor.predict(1.0).isNull()); + + predictor.addImu(1.0, measuredImu(yaw(0.3), Eigen::Vector3d(1, 0, 0), 9.80665, rtabmap::Transform::getIdentity())); + predictor.addImu(1.1, measuredImu(yaw(0.5), Eigen::Vector3d(1, 0, 0), 9.80665, rtabmap::Transform::getIdentity())); + const rtabmap::Transform predicted = predictor.predict(1.05); + ASSERT_FALSE(predicted.isNull()) << "the orientation is known without a pose"; + EXPECT_NEAR(predicted.theta(), 0.4, 1e-5); + EXPECT_NEAR(predicted.x(), 0.0, 1e-9) << "no position without a pose"; + + ImuMotionPredictor withoutImu; + withoutImu.addPose(1.0, rtabmap::Transform::getIdentity()); + EXPECT_TRUE(withoutImu.predict(1.0).isNull()) << "no imu yet"; +} + +TEST(ImuMotionPredictor, follows_a_constant_velocity) +{ + ImuMotionPredictor predictor; + const Eigen::Vector3d velocity(1.0, -0.5, 0.2); + addSamples(predictor, 0.0, 0.3, + [](double) { return Eigen::Quaterniond::Identity(); }, + [](double) { return Eigen::Vector3d::Zero(); }); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), 0)); + EXPECT_TRUE(predictor.velocity().isZero()) << "a single pose has no velocity"; + predictor.addPose(0.1, pose(velocity * 0.1, 0)); + + EXPECT_TRUE(predictor.velocity().isApprox(velocity, 1e-6)); + const rtabmap::Transform predicted = predictor.predict(0.25); + ASSERT_FALSE(predicted.isNull()); + EXPECT_NEAR(predicted.x(), velocity.x() * 0.25, 1e-6); + EXPECT_NEAR(predicted.y(), velocity.y() * 0.25, 1e-6); + EXPECT_NEAR(predicted.z(), velocity.z() * 0.25, 1e-6); +} + +TEST(ImuMotionPredictor, integrates_the_acceleration) +{ + // From rest at t=0 with a constant 2 m/s^2: p = t^2, v = 2t. The velocity at the + // second pose is the instantaneous one, not the average over the interval, and the + // prediction keeps accelerating. An IMU without accelerometer gives a constant + // velocity instead. + const double a = 2.0; + for(bool withAccelerometer : {true, false}) + { + ImuMotionPredictor predictor; + addSamples(predictor, -0.05, 0.3, + [](double) { return Eigen::Quaterniond::Identity(); }, + [&](double t) { return Eigen::Vector3d(t < 0.0 ? 0.0 : a, 0, 0); }, + withAccelerometer); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), 0)); + predictor.addPose(0.1, pose(Eigen::Vector3d(0.5*a*0.01, 0, 0), 0)); + + const double t = 0.2; + const rtabmap::Transform predicted = predictor.predict(t); + ASSERT_FALSE(predicted.isNull()); + if(withAccelerometer) + { + EXPECT_NEAR(predictor.velocity().x(), a*0.1, 1e-6); + EXPECT_NEAR(predicted.x(), 0.5*a*t*t, 1e-6); + } + else + { + // Constant velocity model: the average velocity over the last interval. + EXPECT_NEAR(predictor.velocity().x(), 0.5*a*0.1, 1e-6); + EXPECT_NEAR(predicted.x(), 0.5*a*0.01 + 0.5*a*0.1*(t-0.1), 1e-6); + } + } +} + +TEST(ImuMotionPredictor, expresses_the_imu_in_the_odometry_frame) +{ + // The IMU's world frame and the odometry frame differ by 90 degrees of yaw. The base + // turns at 1 rad/s and accelerates along the IMU world's x, which is the odometry's y. + const double rate = 1.0; + const double offset = M_PI/2.0; + const double a = 3.0; + ImuMotionPredictor predictor; + addSamples(predictor, 0.0, 0.3, + [&](double t) { return yaw(rate*t); }, + [&](double) { return Eigen::Vector3d(a, 0, 0); }); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), offset)); + predictor.addPose(0.1, pose(Eigen::Vector3d(0, 0.5*a*0.01, 0), offset + rate*0.1)); + + const double t = 0.25; + const rtabmap::Transform predicted = predictor.predict(t); + ASSERT_FALSE(predicted.isNull()); + EXPECT_NEAR(predicted.x(), 0.0, 1e-6); + EXPECT_NEAR(predicted.y(), 0.5*a*t*t, 1e-6); + EXPECT_NEAR(predicted.theta(), offset + rate*t, 1e-5); +} + +TEST(ImuMotionPredictor, a_lost_pose_resets_the_prediction) +{ + ImuMotionPredictor predictor; + addSamples(predictor, 0.0, 0.5, + [](double) { return Eigen::Quaterniond::Identity(); }, + [](double) { return Eigen::Vector3d::Zero(); }); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), 0)); + predictor.addPose(0.1, pose(Eigen::Vector3d(0.1, 0, 0), 0)); + ASSERT_FALSE(predictor.velocity().isZero()); + + predictor.addPose(0.2, rtabmap::Transform()); + EXPECT_TRUE(predictor.predict(0.25).isIdentity()) << "orientation only, which is constant here"; + + // After a reset of the odometry, the pose jumps: no velocity across it. + predictor.addPose(0.3, pose(Eigen::Vector3d(10, 0, 0), 0)); + EXPECT_TRUE(predictor.velocity().isZero()); + EXPECT_NEAR(predictor.predict(0.4).x(), 10.0, 1e-6); +} + +TEST(ImuMotionPredictor, poses_too_far_apart_give_no_velocity) +{ + ImuMotionPredictor predictor(0.5); + addSamples(predictor, 0.0, 1.5, + [](double) { return Eigen::Quaterniond::Identity(); }, + [](double) { return Eigen::Vector3d::Zero(); }); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), 0)); + predictor.addPose(1.0, pose(Eigen::Vector3d(1, 0, 0), 0)); + EXPECT_TRUE(predictor.velocity().isZero()); +} + +TEST(ImuMotionPredictor, a_longer_window_averages_out_the_pose_noise) +{ + // 1 m/s along x, with odometry poses alternating 1 cm on each side of the truth: the + // worst case for a velocity differenced over one frame, which sees 0.2 m/s of noise. + for(double window : {0.0, 0.5}) + { + ImuMotionPredictor predictor(1.0, window); + addSamples(predictor, 0.0, 2.0, + [](double) { return Eigen::Quaterniond::Identity(); }, + [](double) { return Eigen::Vector3d::Zero(); }); + double maxError = 0.0; + for(int i=0; i<=15; ++i) + { + const double t = i*0.1; + predictor.addPose(t, pose(Eigen::Vector3d(t + (i%2?0.01:-0.01), 0, 0), 0)); + if(i >= 10) + { + maxError = std::max(maxError, std::fabs(predictor.velocity().x() - 1.0)); + } + } + if(window == 0.0) + { + EXPECT_NEAR(maxError, 0.2, 1e-6) << "differenced over one frame"; + } + else + { + EXPECT_LT(maxError, 0.05) << "differenced over half a second"; + } + } +} + +TEST(ImuMotionPredictor, keeps_only_the_samples_since_the_last_pose) +{ + ImuMotionPredictor predictor; + addSamples(predictor, 0.0, 0.2, + [](double) { return Eigen::Quaterniond::Identity(); }, + [](double) { return Eigen::Vector3d::Zero(); }); + ASSERT_EQ(predictor.samples(), 41u); + predictor.addPose(0.1025, pose(Eigen::Vector3d::Zero(), 0)); + // 0.100 (the last one before the pose, to interpolate at its stamp) .. 0.200 + EXPECT_EQ(predictor.samples(), 21u); + + predictor.reset(); + EXPECT_EQ(predictor.samples(), 0u); + EXPECT_TRUE(predictor.predict(0.2).isNull()); +} + + +TEST(ImuMotionPredictor, removes_gravity_from_what_the_imu_measures) +{ + // The IMU is mounted rolled by 90 degrees on a level base at rest: it measures gravity + // along its own y. Once removed, nothing moves, and the base stays level. + const rtabmap::Transform baseToImu(0, 0, 0, M_PI/2.0, 0, 0); + const Eigen::Quaterniond worldToImu = baseToImu.getQuaterniond(); + ImuMotionPredictor predictor; + for(int i=0; i<=60; ++i) + { + predictor.addImu(i*0.005, measuredImu(worldToImu, Eigen::Vector3d::Zero(), 9.80665, baseToImu)); + } + predictor.addPose(0.0, rtabmap::Transform::getIdentity()); + const rtabmap::Transform predicted = predictor.predict(0.3); + ASSERT_FALSE(predicted.isNull()); + EXPECT_NEAR(predicted.getNorm(), 0.0, 1e-6) << "gravity was not removed"; + EXPECT_TRUE(predicted.getQuaterniond().isApprox(Eigen::Quaterniond::Identity(), 1e-6)) << "the base orientation is the imu's, unmounted"; +} + +TEST(ImuMotionPredictor, removes_the_gravity_it_is_given) +{ + // On the Moon, at rest: the IMU measures 1.62 m/s^2 up. + ImuMotionPredictor predictor(1.0, 0.5, 1.62); + for(int i=0; i<=60; ++i) + { + predictor.addImu(i*0.005, measuredImu(Eigen::Quaterniond::Identity(), Eigen::Vector3d::Zero(), 1.62, + rtabmap::Transform::getIdentity())); + } + predictor.addPose(0.0, rtabmap::Transform::getIdentity()); + EXPECT_NEAR(predictor.predict(0.3).z(), 0.0, 1e-6); + + // Earth's gravity removed from the same measurement: it looks like falling. + ImuMotionPredictor earth; + for(int i=0; i<=60; ++i) + { + earth.addImu(i*0.005, measuredImu(Eigen::Quaterniond::Identity(), Eigen::Vector3d::Zero(), 1.62, + rtabmap::Transform::getIdentity())); + } + earth.addPose(0.0, rtabmap::Transform::getIdentity()); + EXPECT_NEAR(earth.predict(0.2).z(), 0.5*(1.62-9.80665)*0.04, 1e-6); +} + +TEST(ImuMotionPredictor, an_imu_without_acceleration_is_not_a_free_fall) +{ + ImuMotionPredictor predictor; + const Eigen::Quaterniond q(Eigen::AngleAxisd(0.3, Eigen::Vector3d::UnitZ())); + for(int i=0; i<=60; ++i) + { + predictor.addImu(i*0.005, rtabmap::IMU(cv::Vec4d(q.x(), q.y(), q.z(), q.w()), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(0,0,0), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(0,0,0), cv::Mat::eye(3,3,CV_64FC1), + rtabmap::Transform::getIdentity())); + } + predictor.addPose(0.0, rtabmap::Transform::getIdentity()); + EXPECT_NEAR(predictor.predict(0.3).getNorm(), 0.0, 1e-6); + + // And one without orientation is ignored altogether. + ImuMotionPredictor noOrientation; + noOrientation.addImu(0.0, rtabmap::IMU(cv::Vec4d(0,0,0,0), cv::Mat(), + cv::Vec3d(0,0,0), cv::Mat(), cv::Vec3d(0,0,9.8), cv::Mat(), rtabmap::Transform::getIdentity())); + EXPECT_EQ(noOrientation.samples(), 0u); +} + +TEST(ImuMotionPredictor, predicts_the_same_whatever_the_order_samples_and_predictions_come_in) +{ + // The integration since the last pose is kept between predictions: it must not go + // stale when samples arrive after a prediction, or out of order. Compared with a + // predictor given every sample before predicting anything. + auto acceleration = [](double t) { return Eigen::Vector3d(std::sin(20*t), std::cos(15*t), 0.3*t); }; + auto orientation = [](double t) { return yaw(0.5*t); }; + auto imu = [&](double t) { return measuredImu(orientation(t), acceleration(t), 9.80665, rtabmap::Transform::getIdentity()); }; + + ImuMotionPredictor reference; + for(int i=0; i<=80; ++i) reference.addImu(i*0.005, imu(i*0.005)); + reference.addPose(0.0, rtabmap::Transform::getIdentity()); + reference.addPose(0.1, pose(Eigen::Vector3d(0.05, 0, 0), 0.05)); + + ImuMotionPredictor incremental; + for(int i=0; i<=40; ++i) if(i != 30) incremental.addImu(i*0.005, imu(i*0.005)); + incremental.addPose(0.0, rtabmap::Transform::getIdentity()); + incremental.addPose(0.1, pose(Eigen::Vector3d(0.05, 0, 0), 0.05)); // velocity needs up to 0.1: covered + incremental.predict(0.12); // integrates without the sample at 0.15 + incremental.predict(0.3); // beyond the newest sample (0.2) + incremental.addImu(0.15, imu(0.15)); // out of order + for(int i=41; i<=80; ++i) + { + incremental.addImu(i*0.005, imu(i*0.005)); + if(i % 7 == 0) incremental.predict(i*0.005 - 0.0012); + } + + for(double t : {0.1, 0.1013, 0.15, 0.2337, 0.4, 0.45}) + { + const rtabmap::Transform a = reference.predict(t); + const rtabmap::Transform b = incremental.predict(t); + ASSERT_FALSE(a.isNull()); + ASSERT_FALSE(b.isNull()); + EXPECT_NEAR(a.x(), b.x(), 1e-6) << "t=" << t; + EXPECT_NEAR(a.y(), b.y(), 1e-6) << "t=" << t; + EXPECT_NEAR(a.z(), b.z(), 1e-6) << "t=" << t; + } +} + +TEST(ImuMotionPredictor, holds_the_acceleration_at_the_pose_until_a_newer_sample) +{ + // The pose comes after the newest sample: the acceleration there is held from that + // sample, until a newer one says otherwise. + ImuMotionPredictor reference; + ImuMotionPredictor incremental; + auto imu = [](double a) { return measuredImu(Eigen::Quaterniond::Identity(), Eigen::Vector3d(a, 0, 0), 9.80665, rtabmap::Transform::getIdentity()); }; + for(ImuMotionPredictor * p : {&reference, &incremental}) + { + p->addImu(0.0, imu(0.0)); + p->addImu(0.1, imu(0.0)); + p->addPose(0.15, rtabmap::Transform::getIdentity()); + } + incremental.predict(0.2); + for(ImuMotionPredictor * p : {&reference, &incremental}) + { + p->addImu(0.2, imu(4.0)); + p->addImu(0.3, imu(4.0)); + } + EXPECT_NEAR(reference.predict(0.3).x(), incremental.predict(0.3).x(), 1e-9); +} diff --git a/guilib/src/PreferencesDialog.cpp b/guilib/src/PreferencesDialog.cpp index 31345853..3a2d9d98 100644 --- a/guilib/src/PreferencesDialog.cpp +++ b/guilib/src/PreferencesDialog.cpp @@ -1504,6 +1504,8 @@ PreferencesDialog::PreferencesDialog(QWidget * parent) : _ui->odom_flow_scanKeyframeThr->setObjectName(Parameters::kOdomScanKeyFrameThr().c_str()); _ui->odom_flow_guessMotion->setObjectName(Parameters::kOdomGuessMotion().c_str()); _ui->odom_guess_smoothing_delay->setObjectName(Parameters::kOdomGuessSmoothingDelay().c_str()); + _ui->odom_guess_imu_acceleration->setObjectName(Parameters::kOdomGuessImuAcceleration().c_str()); + _ui->odom_imu_gravity->setObjectName(Parameters::kOdomImuGravity().c_str()); _ui->odom_imageDecimation->setObjectName(Parameters::kOdomImageDecimation().c_str()); _ui->odom_alignWithGround->setObjectName(Parameters::kOdomAlignWithGround().c_str()); _ui->odom_lidar_deskewing->setObjectName(Parameters::kOdomDeskewing().c_str()); diff --git a/guilib/src/ui/preferencesDialog.ui b/guilib/src/ui/preferencesDialog.ui index d45752b9..bbb6e2bd 100644 --- a/guilib/src/ui/preferencesDialog.ui +++ b/guilib/src/ui/preferencesDialog.ui @@ -16427,7 +16427,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Data buffer size (0 means inf). @@ -16570,7 +16570,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1.000000000000000 @@ -16590,7 +16590,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 9999 @@ -16600,7 +16600,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + [Geometry] Create a new keyframe when the number of inliers drops under this threshold. Setting value to 0 means that a keyframe is created for each processed frame. @@ -16613,7 +16613,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Decimation of the RGB image before registration. If depth size is larger than decimated RGB size, depth is decimated to be always at most equal to RGB size. If Visual Registration -> Visual Feature -> Depth as Mask is true and if depth is smaller than decimated RGB, depth may be interpolated to match RGB size for feature detection. @@ -16709,7 +16709,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Test odometry @@ -16729,14 +16729,14 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 999999 - + [Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame. @@ -16749,7 +16749,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1 @@ -16778,7 +16778,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1.000000000000000 @@ -16794,7 +16794,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if a filtering strategy is set or the delay is below the odometry rate. + Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if a filtering strategy is set or the delay is below the odometry rate. Recommended (~0.5 s) with the IMU acceleration below: the velocity is then corrected by the IMU acceleration, so that it is the velocity at the last frame rather than a delayed average. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses. true @@ -16831,7 +16831,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + [Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame. @@ -16844,7 +16844,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Lidar deskewing. If input lidar has time channel, it will be deskewed with a constant motion model (with IMU orientation and/or guess if provided). @@ -16857,13 +16857,65 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + + + + + + + + + + + + With an IMU giving orientation and linear acceleration, predict the translation of the motion guess by integrating the IMU acceleration from the velocity at the previous frame, instead of with a constant velocity. That velocity is the average over the guess smoothing delay above (set it to about 0.5 s), carried to the previous frame with the acceleration measured since: it is not delayed like the average alone. Ignored by odometry approaches processing the IMU themselves. + + + true + + + Qt::LinksAccessibleByMouse|Qt::TextSelectableByMouse + + + + + + + m/s² + + + 5 + + + 100.000000000000000 + + + 0.010000000000000 + + + 9.806650000000000 + + + + + + + Gravity magnitude removed from the IMU linear acceleration with the IMU acceleration above. Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. + + + true + + + Qt::LinksAccessibleByMouse|Qt::TextSelectableByMouse + + + diff --git a/package.xml b/package.xml index fc53b1f8..ade42f62 100644 --- a/package.xml +++ b/package.xml @@ -1,7 +1,7 @@ rtabmap - 0.24.0 + 0.24.1 RTAB-Map's standalone library. RTAB-Map is a RGB-D SLAM approach with real-time constraints. Mathieu Labbe Mathieu Labbe From de1761bd6d9d14116ee19692800744aaebe2d4bf Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sat, 10 Oct 2026 18:26:20 -0700 Subject: [PATCH 2/7] simplified parameters --- .../include/rtabmap/core/ImuMotionPredictor.h | 14 ++- corelib/include/rtabmap/core/Odometry.h | 3 +- corelib/include/rtabmap/core/Parameters.h | 7 +- corelib/include/rtabmap/core/util3d.h | 22 ++++ corelib/src/ImuMotionPredictor.cpp | 11 +- corelib/src/Odometry.cpp | 95 ++++++++-------- corelib/src/OdometryThread.cpp | 17 ++- corelib/src/util3d.cpp | 104 ++++++++++++------ corelib/test/test_imumotionpredictor.cpp | 15 +++ corelib/test/test_util3d.cpp | 38 +++++++ guilib/src/PreferencesDialog.cpp | 1 - guilib/src/ui/preferencesDialog.ui | 66 ++++------- 12 files changed, 250 insertions(+), 143 deletions(-) diff --git a/corelib/include/rtabmap/core/ImuMotionPredictor.h b/corelib/include/rtabmap/core/ImuMotionPredictor.h index 4b5e9d19..7517f261 100644 --- a/corelib/include/rtabmap/core/ImuMotionPredictor.h +++ b/corelib/include/rtabmap/core/ImuMotionPredictor.h @@ -66,16 +66,20 @@ class RTABMAP_CORE_EXPORT ImuMotionPredictor { public: /** - * Without acceleration in the IMU samples, the position follows the last velocity + * The IMU acceleration is used only with a velocity window (> 0): over a single frame, + * the velocity is too noisy to be carried forward with it. Without it (window of 0, or + * no acceleration in the IMU samples), the position follows the last velocity * (constant velocity model). * * @param maxPoseInterval odometry poses older (s) than this are not used to estimate * the velocity, which is null without one * @param velocityWindow the velocity is estimated from the displacement since the - * newest pose at least this old (s). Over a single frame - * interval, the noise of the odometry poses would be of the order - * of the velocity itself; with the acceleration, a longer window - * still gives the velocity at the last pose, not an average. + * newest pose at least this old (s), corrected by the + * acceleration measured since, so that it is the velocity at the + * last pose, not an average. Over a single frame interval, the + * noise of the odometry poses would be of the order of the + * velocity itself. 0: the displacement since the previous pose, + * without acceleration. * @param gravity magnitude (m/s^2) of the gravity removed from the specific force * given to addImu(), standard gravity by default * diff --git a/corelib/include/rtabmap/core/Odometry.h b/corelib/include/rtabmap/core/Odometry.h index 524e589a..59b3129d 100644 --- a/corelib/include/rtabmap/core/Odometry.h +++ b/corelib/include/rtabmap/core/Odometry.h @@ -160,7 +160,6 @@ private: bool _force3DoF; bool _holonomic; bool guessFromMotion_; - bool guessImuAcceleration_; float guessSmoothingDelay_; int _filteringStrategy; int _particleSize; @@ -191,7 +190,7 @@ private: std::vector stereoModels_; std::vector models_; std::map imus_; - ImuMotionPredictor imuMotionPredictor_; // used with Odom/GuessImuAcceleration only + ImuMotionPredictor imuMotionPredictor_; // fed when IMU is received (motion guess and deskewing) }; } /* namespace rtabmap */ diff --git a/corelib/include/rtabmap/core/Parameters.h b/corelib/include/rtabmap/core/Parameters.h index 07a854a6..cdd4d5fe 100644 --- a/corelib/include/rtabmap/core/Parameters.h +++ b/corelib/include/rtabmap/core/Parameters.h @@ -512,15 +512,14 @@ class RTABMAP_CORE_EXPORT Parameters RTABMAP_PARAM(Odom, KalmanProcessNoise, float, 0.001, "Process noise covariance value."); RTABMAP_PARAM(Odom, KalmanMeasurementNoise, float, 0.01, "Process measurement covariance value."); RTABMAP_PARAM(Odom, GuessMotion, bool, true, "Guess next transformation from the last motion computed."); - RTABMAP_PARAM(Odom, GuessImuAcceleration, bool, false, uFormat("With an IMU giving orientation and linear acceleration, predict the translation of the motion guess (with \"%s\") by integrating the IMU acceleration from the velocity at the previous frame, instead of with a constant velocity. That velocity is the average over \"%s\" (from the displacement over that delay), carried to the previous frame with the acceleration measured since: it is not delayed like the average alone. Set \"%s\" to about 0.5 s: over a single frame, the noise of the poses is of the order of the velocity itself. Ignored by odometry approaches processing the IMU themselves.", kOdomGuessMotion().c_str(), kOdomGuessSmoothingDelay().c_str(), kOdomGuessSmoothingDelay().c_str())); - RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration with \"%s\". Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's.", kOdomGuessImuAcceleration().c_str())); - RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate. Recommended (~0.5 s) with \"%s\": the translational velocity is then the displacement over this delay corrected by the IMU acceleration, so that it is the velocity at the last frame rather than a delayed average. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing (\"%s\"), where its noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomGuessImuAcceleration().c_str(), kOdomDeskewing().c_str())); + RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration (used with \"%s\" > 0). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's.", kOdomGuessSmoothingDelay().c_str())); + RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is then the displacement over this delay corrected by the acceleration measured since, so that it is the velocity at the last frame rather than a delayed average, and the motion guess and lidar deskewing (\"%s\") integrate the acceleration from it. Recommended (~0.5 s) with an IMU. With 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomDeskewing().c_str())); RTABMAP_PARAM(Odom, KeyFrameThr, float, 0.3, "[Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, VisKeyFrameThr, int, 150, "[Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, ScanKeyFrameThr, float, 0.9, "[Geometry] Create a new keyframe when the number of ICP inliers drops under this ratio of points in last frame's scan. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, ImageDecimation, unsigned int, 1, uFormat("Decimation of the RGB image before registration. If depth size is larger than decimated RGB size, depth is decimated to be always at most equal to RGB size. If %s is true and if depth is smaller than decimated RGB, depth may be interpolated to match RGB size for feature detection.", kVisDepthAsMask().c_str())); RTABMAP_PARAM(Odom, AlignWithGround, bool, false, "Align odometry with the ground on initialization."); - RTABMAP_PARAM(Odom, Deskewing, bool, true, "Lidar deskewing. If input lidar has time channel, it will be deskewed with a constant motion model (with IMU orientation and/or guess if provided)."); + RTABMAP_PARAM(Odom, Deskewing, bool, true, uFormat("Lidar deskewing. If input lidar has time channel, it will be deskewed. With an IMU, the pose of every point is predicted from the previous frame: orientation from the IMU, translation from the velocity (with the IMU acceleration if \"%s\" > 0). Without IMU, with a constant motion model (or the guess if provided).", kOdomGuessSmoothingDelay().c_str())); // Odometry Frame-to-Map RTABMAP_PARAM(OdomF2M, MaxSize, int, 2000, "[Visual] Local map size: If > 0 (example 5000), the odometry will maintain a local map of X maximum words."); diff --git a/corelib/include/rtabmap/core/util3d.h b/corelib/include/rtabmap/core/util3d.h index 70264a66..b3136117 100644 --- a/corelib/include/rtabmap/core/util3d.h +++ b/corelib/include/rtabmap/core/util3d.h @@ -28,6 +28,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #ifndef UTIL3D_H_ #define UTIL3D_H_ +#include #include "rtabmap/core/rtabmap_core_export.h" #include @@ -1295,6 +1296,27 @@ LaserScan RTABMAP_CORE_EXPORT deskew( double inputStamp, const rtabmap::Transform & velocity); +/** + * @brief Deskews a lidar scan with a motion given by the caller. + * + * Same as the velocity overload, but the motion during the sweep comes from @p motion, + * for instance a pose predicted from an IMU. + * + * @param input scan with a time channel (`kXYZIT` or `kXYZIRT`) + * @param inputStamp stamp of the scan, which the time channel is relative to + * @param motion for a stamp (s) in the sweep, the pose of the scan's base frame at that + * stamp relative to the base frame at @p inputStamp; a null transform + * aborts deskewing + * @param slerp call @p motion only for the first and last points and interpolate in + * between, instead of calling it for every time of the sweep + * @return the deskewed scan, empty on error + */ +LaserScan RTABMAP_CORE_EXPORT deskew( + const LaserScan & input, + double inputStamp, + const std::function & motion, + bool slerp = false); + } // namespace util3d } // namespace rtabmap diff --git a/corelib/src/ImuMotionPredictor.cpp b/corelib/src/ImuMotionPredictor.cpp index f00a1959..ff12698d 100644 --- a/corelib/src/ImuMotionPredictor.cpp +++ b/corelib/src/ImuMotionPredictor.cpp @@ -133,7 +133,8 @@ void ImuMotionPredictor::addPose(double stamp, const rtabmap::Transform & pose) pose.y() - reference->second.y(), pose.z() - reference->second.z()); - if(!samples_.empty() && + if(velocityWindow_ > 0.0 && + !samples_.empty() && samples_.begin()->first <= reference->first && samples_.rbegin()->first >= stamp) { @@ -147,7 +148,7 @@ void ImuMotionPredictor::addPose(double stamp, const rtabmap::Transform & pose) } else { - // Average velocity over the interval (the samples don't cover it) + // Average velocity over the interval (no window, or the samples don't cover it) velocity = displacement / interval; } } @@ -202,7 +203,11 @@ rtabmap::Transform ImuMotionPredictor::predict(double stamp) const // D(t) is the displacement due to the change of velocity since t0, V(s) = integral_t0^s a(u) du. Eigen::Vector3d position(pose_.x(), pose_.y(), pose_.z()); // p0 position += velocity_ * (stamp - poseStamp_); // + v0*(t-t0) - if(stamp >= poseStamp_) + if(velocityWindow_ <= 0.0) + { + // No acceleration without a velocity window: constant velocity + } + else if(stamp >= poseStamp_) { // + D(t). D and V are kept at every sample stamp since t0 (integrated_), so only // the part from the last sample ti <= t is integrated here, with dt = t-ti: diff --git a/corelib/src/Odometry.cpp b/corelib/src/Odometry.cpp index 51c4971b..ba02ee78 100644 --- a/corelib/src/Odometry.cpp +++ b/corelib/src/Odometry.cpp @@ -134,7 +134,6 @@ Odometry::Odometry(const rtabmap::ParametersMap & parameters) : _force3DoF(Parameters::defaultRegForce3DoF()), _holonomic(Parameters::defaultOdomHolonomic()), guessFromMotion_(Parameters::defaultOdomGuessMotion()), - guessImuAcceleration_(Parameters::defaultOdomGuessImuAcceleration()), guessSmoothingDelay_(Parameters::defaultOdomGuessSmoothingDelay()), _filteringStrategy(Parameters::defaultOdomFilteringStrategy()), _particleSize(Parameters::defaultOdomParticleSize()), @@ -161,13 +160,12 @@ Odometry::Odometry(const rtabmap::ParametersMap & parameters) : Parameters::parse(parameters, Parameters::kRegForce3DoF(), _force3DoF); Parameters::parse(parameters, Parameters::kOdomHolonomic(), _holonomic); Parameters::parse(parameters, Parameters::kOdomGuessMotion(), guessFromMotion_); - Parameters::parse(parameters, Parameters::kOdomGuessImuAcceleration(), guessImuAcceleration_); Parameters::parse(parameters, Parameters::kOdomGuessSmoothingDelay(), guessSmoothingDelay_); - if(guessImuAcceleration_) { float imuGravity = Parameters::defaultOdomImuGravity(); Parameters::parse(parameters, Parameters::kOdomImuGravity(), imuGravity); - // The velocity is estimated over the smoothing delay + // The velocity is estimated over the smoothing delay, and the IMU acceleration used + // only with one (> 0) imuMotionPredictor_ = ImuMotionPredictor(1.0, guessSmoothingDelay_, imuGravity); } Parameters::parse(parameters, Parameters::kOdomFillInfoData(), _fillInfoData); @@ -351,10 +349,7 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet imus_.erase(imus_.begin()); } - if(guessImuAcceleration_) - { - imuMotionPredictor_.addImu(data.stamp(), data.imu()); - } + imuMotionPredictor_.addImu(data.stamp(), data.imu()); } else { @@ -661,10 +656,11 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet orientation.r11(), orientation.r12(), orientation.r13(), guess.x(), orientation.r21(), orientation.r22(), orientation.r23(), guess.y(), orientation.r31(), orientation.r32(), orientation.r33(), guess.z()); - if(guessFromMotion_ && guessImuAcceleration_ && imuMotionPredictor_.hasPose()) + if(guessFromMotion_ && guessSmoothingDelay_ > 0.0f && imuMotionPredictor_.hasPose()) { // Translation (and orientation) predicted from the previous pose with the - // IMU acceleration, instead of a constant velocity. + // IMU acceleration, instead of a constant velocity: the velocity is the one + // over the smoothing delay, carried to the previous frame with the IMU. Transform predicted = imuMotionPredictor_.predict(data.stamp()); if(!predicted.isNull()) { @@ -688,15 +684,48 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet UTimer time; - // Deskewing lidar + // Deskewing lidar, if the scan has a time spread (not already deskewed: deskewing zeroes + // the time channel) + const bool scanHasTimeSpread = + !data.laserScanRaw().empty() && + data.laserScanRaw().hasTime() && + data.laserScanRaw().data().ptr(0, data.laserScanRaw().size()-1)[data.laserScanRaw().getTimeOffset()] != + data.laserScanRaw().data().ptr(0, 0)[data.laserScanRaw().getTimeOffset()]; if( _deskewing && - !data.laserScanRaw().empty() && - data.laserScanRaw().hasTime() && + scanHasTimeSpread && + !imus_.empty() && + !imuMotionPredictor_.predict(data.stamp()).isNull()) + { + UDEBUG("Deskewing with IMU begin"); + // Every point's pose predicted with the IMU since the previous frame: orientation + // from the IMU, translation from the velocity (carried with the IMU acceleration + // with a smoothing delay). Before the first pose, only the orientation. + const Transform referenceInverse = imuMotionPredictor_.predict(data.stamp()).inverse(); + auto motion = [&](double stamp) + { + Transform pose = imuMotionPredictor_.predict(stamp); + if(pose.isNull()) + { + return pose; + } + pose = referenceInverse * pose; + return _force3DoF?pose.to3DoF():pose; + }; + LaserScan scanDeskewed = util3d::deskew(data.laserScanRaw(), data.stamp(), motion); + if(!scanDeskewed.isEmpty()) + { + data.setLaserScan(scanDeskewed); + } + info->timeDeskewing = time.ticks(); + UDEBUG("Deskewing end"); + } + else if( _deskewing && + scanHasTimeSpread && dt > 0 && !guess.isNull()) { UDEBUG("Deskewing begin"); - // Recompute velocity + // Constant velocity float vx,vy,vz, vroll,vpitch,vyaw; guess.getTranslationAndEulerAngles(vx,vy,vz, vroll,vpitch,vyaw); @@ -708,38 +737,6 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet vpitch /= dt; vyaw /= dt; - if(!imus_.empty()) - { - float scanTime = - data.laserScanRaw().data().ptr(0, data.laserScanRaw().size()-1)[data.laserScanRaw().getTimeOffset()] - - data.laserScanRaw().data().ptr(0, 0)[data.laserScanRaw().getTimeOffset()]; - - // replace orientation velocity based on IMU (if available) - Transform imuFirstScan = Transform::getTransform(imus_, - data.stamp() + - data.laserScanRaw().data().ptr(0, 0)[data.laserScanRaw().getTimeOffset()]); - Transform imuLastScan = Transform::getTransform(imus_, - data.stamp() + - data.laserScanRaw().data().ptr(0, data.laserScanRaw().size()-1)[data.laserScanRaw().getTimeOffset()]); - if(!imuFirstScan.isNull() && !imuLastScan.isNull()) - { - Transform orientation = imuFirstScan.inverse() * imuLastScan; - orientation.getEulerAngles(vroll, vpitch, vyaw); - if(_force3DoF) - { - vroll=0; - vpitch=0; - vyaw /= scanTime; - } - else - { - vroll /= scanTime; - vpitch /= scanTime; - vyaw /= scanTime; - } - } - } - Transform velocity(vx,vy,vz,vroll,vpitch,vyaw); LaserScan scanDeskewed = util3d::deskew(data.laserScanRaw(), data.stamp(), velocity); if(!scanDeskewed.isEmpty()) @@ -902,7 +899,7 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet } } - if(t.isNull() && guessImuAcceleration_) + if(t.isNull()) { // Lost: no velocity can be estimated across the reset that follows imuMotionPredictor_.addPose(data.stamp(), Transform()); @@ -1063,11 +1060,11 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet velocityGuess_.setNull(); } - if(guessImuAcceleration_) { const Transform newPose = _pose * t; imuMotionPredictor_.addPose(data.stamp(), newPose); - if(!velocityGuess_.isNull() && _filteringStrategy != 1 && particleFilters_.empty()) + if(guessSmoothingDelay_ > 0.0f && !imus_.empty() && + !velocityGuess_.isNull() && _filteringStrategy != 1 && particleFilters_.empty()) { // The translational velocity over the smoothing delay, carried to this // frame with the IMU acceleration (see ImuMotionPredictor), in this frame. diff --git a/corelib/src/OdometryThread.cpp b/corelib/src/OdometryThread.cpp index 3c7b079a..82af2eda 100644 --- a/corelib/src/OdometryThread.cpp +++ b/corelib/src/OdometryThread.cpp @@ -34,6 +34,8 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include "rtabmap/utilite/ULogger.h" #include "rtabmap/utilite/UTimer.h" +#include + namespace rtabmap { OdometryThread::OdometryThread(Odometry * odometry, unsigned int dataBufferMaxSize) : @@ -236,13 +238,24 @@ bool OdometryThread::getData(SensorEvent & event) { if(!_dataBuffer.empty()) { - // Send IMU up to stamp greater than image (OpenVINS needs this). + // Send IMU up to stamp greater than image (OpenVINS needs this). For a lidar + // scan with a time channel, up to the end of its sweep: deskewing predicts the + // pose of every point with the IMU (approaches processing the IMU themselves + // get it as before). + double imuUntil = _dataBuffer.front().data().stamp(); + const LaserScan & scan = _dataBuffer.front().data().laserScanRaw(); + if(!_odometry->canProcessAsyncIMU() && !scan.isEmpty() && scan.hasTime()) + { + imuUntil += std::max(0.0f, std::max( + scan.data().ptr(0, 0)[scan.getTimeOffset()], + scan.data().ptr(0, scan.size()-1)[scan.getTimeOffset()])); + } while(!_imuBuffer.empty()) { _odometry->process(_imuBuffer.front()); double stamp =_imuBuffer.front().stamp(); _imuBuffer.pop_front(); - if(stamp > _dataBuffer.front().data().stamp()) { + if(stamp > imuUntil) { break; } } diff --git a/corelib/src/util3d.cpp b/corelib/src/util3d.cpp index 895e158a..9c7ecffe 100644 --- a/corelib/src/util3d.cpp +++ b/corelib/src/util3d.cpp @@ -3822,17 +3822,12 @@ pcl::PointCloud::Ptr loadCloud( return util3d::transformPointCloud(cloud, transform); } -LaserScan deskew( +static LaserScan deskewImpl( const LaserScan & input, double inputStamp, - const rtabmap::Transform & velocity) + const std::function & motion, + bool slerp) { - if(velocity.isNull()) - { - UERROR("velocity should be valid!"); - return LaserScan(); - } - if(!input.hasTime()) { UERROR("input scan doesn't have a \"time\" channel! Supported formats: \"%s\", \"%s\".", @@ -3861,33 +3856,28 @@ LaserScan deskew( return LaserScan(); } + // With slerp, the poses of the base frame at the first and last stamps (relative to + // the base frame at inputStamp), interpolated in between rtabmap::Transform firstPose; rtabmap::Transform lastPose; - - float vx,vy,vz, vroll,vpitch,vyaw; - velocity.getTranslationAndEulerAngles(vx,vy,vz, vroll,vpitch,vyaw); - - // 1- The pose of base frame in odom frame at first stamp - // 2- The pose of base frame in odom frame at last stamp - double dt1 = firstStamp - inputStamp; - double dt2 = lastStamp - inputStamp; - - firstPose = rtabmap::Transform(vx*dt1, vy*dt1, vz*dt1, vroll*dt1, vpitch*dt1, vyaw*dt1); - lastPose = rtabmap::Transform(vx*dt2, vy*dt2, vz*dt2, vroll*dt2, vpitch*dt2, vyaw*dt2); - - if(firstPose.isNull()) + if(slerp) { - UERROR("Could not get transform between stamps %f and %f!", - firstStamp, - inputStamp); - return LaserScan(); - } - if(lastPose.isNull()) - { - UERROR("Could not get transform between stamps %f and %f!", - lastStamp, - inputStamp); - return LaserScan(); + firstPose = motion(firstStamp); + lastPose = motion(lastStamp); + if(firstPose.isNull()) + { + UERROR("Could not get transform between stamps %f and %f!", + firstStamp, + inputStamp); + return LaserScan(); + } + if(lastPose.isNull()) + { + UERROR("Could not get transform between stamps %f and %f!", + lastStamp, + inputStamp); + return LaserScan(); + } } double stamp; @@ -3918,7 +3908,12 @@ LaserScan deskew( { const float * inputPtr = input.data().ptr(0, u); stamp = inputStamp + inputPtr[offsetTime]; - rtabmap::Transform transform = firstPose.interpolate((stamp-firstStamp) / scanTime, lastPose); + rtabmap::Transform transform = slerp?firstPose.interpolate((stamp-firstStamp) / scanTime, lastPose):motion(stamp); + if(transform.isNull()) + { + UERROR("Could not get transform between stamps %f and %f!", stamp, inputStamp); + return LaserScan(); + } for(int v=0; v(v, 0); stamp = inputStamp + inputPtr[offsetTime]; - rtabmap::Transform transform = firstPose.interpolate((stamp-firstStamp) / scanTime, lastPose); + rtabmap::Transform transform = slerp?firstPose.interpolate((stamp-firstStamp) / scanTime, lastPose):motion(stamp); + if(transform.isNull()) + { + UERROR("Could not get transform between stamps %f and %f!", stamp, inputStamp); + return LaserScan(); + } for(int u=0; u & motion, + bool slerp) +{ + if(!motion) + { + UERROR("motion should be set!"); + return LaserScan(); + } + return deskewImpl(input, inputStamp, motion, slerp); +} + } diff --git a/corelib/test/test_imumotionpredictor.cpp b/corelib/test/test_imumotionpredictor.cpp index f5a8f9ba..5aff410d 100644 --- a/corelib/test/test_imumotionpredictor.cpp +++ b/corelib/test/test_imumotionpredictor.cpp @@ -352,3 +352,18 @@ TEST(ImuMotionPredictor, holds_the_acceleration_at_the_pose_until_a_newer_sample } EXPECT_NEAR(reference.predict(0.3).x(), incremental.predict(0.3).x(), 1e-9); } + +TEST(ImuMotionPredictor, ignores_the_acceleration_without_a_velocity_window) +{ + // From rest with a constant 2 m/s^2: with a window of 0, the velocity is the one of the + // last interval and is kept constant, as without IMU. + const double a = 2.0; + ImuMotionPredictor predictor(1.0, 0.0); + addSamples(predictor, -0.05, 0.3, + [](double) { return Eigen::Quaterniond::Identity(); }, + [&](double t) { return Eigen::Vector3d(t < 0.0 ? 0.0 : a, 0, 0); }); + predictor.addPose(0.0, pose(Eigen::Vector3d::Zero(), 0)); + predictor.addPose(0.1, pose(Eigen::Vector3d(0.5*a*0.01, 0, 0), 0)); + EXPECT_NEAR(predictor.velocity().x(), 0.5*a*0.1, 1e-6); + EXPECT_NEAR(predictor.predict(0.2).x(), 0.5*a*0.01 + 0.5*a*0.1*0.1, 1e-6); +} diff --git a/corelib/test/test_util3d.cpp b/corelib/test/test_util3d.cpp index de36a4e7..fcedb10e 100644 --- a/corelib/test/test_util3d.cpp +++ b/corelib/test/test_util3d.cpp @@ -2114,3 +2114,41 @@ TEST(Util3dTest, DeskewValidScan) { // empty scan EXPECT_TRUE(util3d::deskew(LaserScan(), inputStamp, velocity).empty()); } + +TEST(Util3dTest, DeskewWithMotion) { + // Three points measured at y=0, 1 s before, at and 1 s after the scan stamp, while the + // base moves along y as (t-stamp)^2: a motion no constant velocity describes. + cv::Mat data = cv::Mat::zeros(1, 3, CV_32FC(5)); + float * dataPtr = (float*)data.data; + dataPtr[0] = 1; dataPtr[4] = -1; + dataPtr[5] = 1; dataPtr[9] = 0; + dataPtr[10] = 1; dataPtr[14] = 1; + LaserScan scan(data, 3, 10.0f, LaserScan::kXYZIT); + const double inputStamp = 1000.0; + + std::vector stamps; + auto motion = [&](double stamp) { + stamps.push_back(stamp); + const double dt = stamp - inputStamp; + return Transform(0.0, dt*dt, 0.0, 0.0, 0.0, 0.0); + }; + + // Asked for every point time, each point is moved by its own pose + LaserScan result = util3d::deskew(scan, inputStamp, motion); + ASSERT_EQ(result.size(), 3); + EXPECT_EQ(stamps.size(), 3u); + EXPECT_FLOAT_EQ(result.field(0, 1), 1.0f); + EXPECT_FLOAT_EQ(result.field(1, 1), 0.0f); + EXPECT_FLOAT_EQ(result.field(2, 1), 1.0f); + + // With slerp, only the ends are asked for, the middle point is interpolated between them + stamps.clear(); + result = util3d::deskew(scan, inputStamp, motion, true); + ASSERT_EQ(result.size(), 3); + EXPECT_EQ(stamps.size(), 2u); + EXPECT_FLOAT_EQ(result.field(1, 1), 1.0f); + + // A failing motion, or none + EXPECT_TRUE(util3d::deskew(scan, inputStamp, [](double) { return Transform(); }).empty()); + EXPECT_TRUE(util3d::deskew(scan, inputStamp, std::function()).empty()); +} diff --git a/guilib/src/PreferencesDialog.cpp b/guilib/src/PreferencesDialog.cpp index 3a2d9d98..e7d09a41 100644 --- a/guilib/src/PreferencesDialog.cpp +++ b/guilib/src/PreferencesDialog.cpp @@ -1504,7 +1504,6 @@ PreferencesDialog::PreferencesDialog(QWidget * parent) : _ui->odom_flow_scanKeyframeThr->setObjectName(Parameters::kOdomScanKeyFrameThr().c_str()); _ui->odom_flow_guessMotion->setObjectName(Parameters::kOdomGuessMotion().c_str()); _ui->odom_guess_smoothing_delay->setObjectName(Parameters::kOdomGuessSmoothingDelay().c_str()); - _ui->odom_guess_imu_acceleration->setObjectName(Parameters::kOdomGuessImuAcceleration().c_str()); _ui->odom_imu_gravity->setObjectName(Parameters::kOdomImuGravity().c_str()); _ui->odom_imageDecimation->setObjectName(Parameters::kOdomImageDecimation().c_str()); _ui->odom_alignWithGround->setObjectName(Parameters::kOdomAlignWithGround().c_str()); diff --git a/guilib/src/ui/preferencesDialog.ui b/guilib/src/ui/preferencesDialog.ui index bbb6e2bd..6a30578e 100644 --- a/guilib/src/ui/preferencesDialog.ui +++ b/guilib/src/ui/preferencesDialog.ui @@ -16427,7 +16427,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Data buffer size (0 means inf). @@ -16570,7 +16570,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1.000000000000000 @@ -16590,7 +16590,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 9999 @@ -16600,7 +16600,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + [Geometry] Create a new keyframe when the number of inliers drops under this threshold. Setting value to 0 means that a keyframe is created for each processed frame. @@ -16613,7 +16613,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Decimation of the RGB image before registration. If depth size is larger than decimated RGB size, depth is decimated to be always at most equal to RGB size. If Visual Registration -> Visual Feature -> Depth as Mask is true and if depth is smaller than decimated RGB, depth may be interpolated to match RGB size for feature detection. @@ -16709,7 +16709,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + Test odometry @@ -16729,14 +16729,14 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 999999 - + [Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame. @@ -16749,7 +16749,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1 @@ -16778,7 +16778,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + 1.000000000000000 @@ -16794,7 +16794,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if a filtering strategy is set or the delay is below the odometry rate. Recommended (~0.5 s) with the IMU acceleration below: the velocity is then corrected by the IMU acceleration, so that it is the velocity at the last frame rather than a delayed average. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses. + Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if a filtering strategy is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is then corrected by the acceleration measured since, so that it is the velocity at the last frame rather than a delayed average, and the motion guess and lidar deskewing integrate the acceleration from it. Recommended (~0.5 s) with an IMU; with 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses. true @@ -16831,7 +16831,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + [Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame. @@ -16844,37 +16844,10 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + - Lidar deskewing. If input lidar has time channel, it will be deskewed with a constant motion model (with IMU orientation and/or guess if provided). - - - true - - - Qt::LinksAccessibleByMouse|Qt::TextSelectableByMouse - - - - - - - - - - - - - - - - - - - - - With an IMU giving orientation and linear acceleration, predict the translation of the motion guess by integrating the IMU acceleration from the velocity at the previous frame, instead of with a constant velocity. That velocity is the average over the guess smoothing delay above (set it to about 0.5 s), carried to the previous frame with the acceleration measured since: it is not delayed like the average alone. Ignored by odometry approaches processing the IMU themselves. + Lidar deskewing. If input lidar has time channel, it will be deskewed. With an IMU, the pose of every point is predicted from the previous frame: orientation from the IMU, translation from the velocity (with the IMU acceleration if the guess smoothing delay is > 0). Without IMU, with a constant motion model (or the guess if provided). true @@ -16885,6 +16858,13 @@ With <0, the length is estimated once for each unique marker, then re-used fo + + + + + + + m/s² @@ -16903,10 +16883,10 @@ With <0, the length is estimated once for each unique marker, then re-used fo - + - Gravity magnitude removed from the IMU linear acceleration with the IMU acceleration above. Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. + Gravity magnitude removed from the IMU linear acceleration (used with a guess smoothing delay > 0 above). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. true From 9e27446236b4fff21a9dc77a2b9e34c2c8a95b94 Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sat, 10 Oct 2026 22:12:50 -0700 Subject: [PATCH 3/7] refactored odom's imu init orientation --- corelib/src/Odometry.cpp | 105 ++++++++++------- corelib/test/test_odometry.cpp | 209 +++++++++++++++++++++++++++++++++ 2 files changed, 273 insertions(+), 41 deletions(-) diff --git a/corelib/src/Odometry.cpp b/corelib/src/Odometry.cpp index ba02ee78..338f58ce 100644 --- a/corelib/src/Odometry.cpp +++ b/corelib/src/Odometry.cpp @@ -320,40 +320,47 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet { UASSERT_MSG(data.id() >= 0, uFormat("Input data should have ID greater or equal than 0 (id=%d)!", data.id()).c_str()); - // cache imu data - if(!data.imu().empty() && !this->canProcessAsyncIMU()) + if(!data.imu().empty()) { - if(!(data.imu().orientation()[0] == 0.0 && data.imu().orientation()[1] == 0.0 && data.imu().orientation()[2] == 0.0)) + if(!this->canProcessAsyncIMU()) { - Transform orientation(0,0,0, data.imu().orientation()[0], data.imu().orientation()[1], data.imu().orientation()[2], data.imu().orientation()[3]); - // orientation includes roll and pitch but not yaw in local transform - Transform imuT = Transform(data.imu().localTransform().x(),data.imu().localTransform().y(),data.imu().localTransform().z(), 0,0,data.imu().localTransform().theta()) * - orientation* - data.imu().localTransform().rotation().inverse(); - - if( this->getPose().r11() == 1.0f && this->getPose().r22() == 1.0f && this->getPose().r33() == 1.0f && - this->framesProcessed() == 0) + // cache imu data + if(!(data.imu().orientation()[0] == 0.0 && data.imu().orientation()[1] == 0.0 && data.imu().orientation()[2] == 0.0)) { - Eigen::Quaterniond imuQuat = imuT.getQuaterniond(); - Transform previous = this->getPose(); - Transform newFramePose = Transform(previous.x(), previous.y(), previous.z(), imuQuat.x(), imuQuat.y(), imuQuat.z(), imuQuat.w()); - UWARN("Updated initial pose from %s to %s with IMU orientation", previous.prettyPrint().c_str(), newFramePose.prettyPrint().c_str()); - std::map imus = imus_; - this->reset(newFramePose); - imus_ = imus; - } + Transform orientation(0,0,0, data.imu().orientation()[0], data.imu().orientation()[1], data.imu().orientation()[2], data.imu().orientation()[3]); + // orientation includes roll and pitch but not yaw in local transform + Transform imuT = Transform(data.imu().localTransform().x(),data.imu().localTransform().y(),data.imu().localTransform().z(), 0,0,data.imu().localTransform().theta()) * + orientation* + data.imu().localTransform().rotation().inverse(); - imus_.insert(std::make_pair(data.stamp(), imuT)); - if(imus_.size() > 1000) + imus_.insert(std::make_pair(data.stamp(), imuT)); + if(imus_.size() > 1000) + { + imus_.erase(imus_.begin()); + } + + imuMotionPredictor_.addImu(data.stamp(), data.imu()); + } + else { - imus_.erase(imus_.begin()); + UWARN("Received IMU doesn't have orientation set! It is ignored."); } - - imuMotionPredictor_.addImu(data.stamp(), data.imu()); } - else + + // IMU-only update: nothing more to do once the IMU is cached, except for approaches + // processing it themselves. A frame that brings its own features carries no image, + // and a frame whose scene was empty carries no feature either, so neither says + // whether there is a frame at all. The calibration does: it is there when a camera + // produced this data. + if(data.imageRaw().empty() && data.imageCompressed().empty() && + data.laserScanRaw().isEmpty() && data.laserScanCompressed().isEmpty() && + data.cameraModels().empty() && data.stereoCameraModels().empty()) { - UWARN("Received IMU doesn't have orientation set! It is ignored."); + if(this->canProcessAsyncIMU()) + { + this->computeTransform(data, Transform(), info); + } + return Transform(); // Return null on IMU-only updates } } @@ -598,12 +605,40 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet } } + // Initial orientation from the IMU at the first frame's own stamp (unless an initial pose + // with a rotation was given). Not from the first IMU sample received, which can be much + // older (e.g., IMU buffered while waiting for the first frame), nor from the newest one, + // which can be after the stamp (lidar deskewing needs IMU up to the end of the sweep). + if(this->framesProcessed() == 0 && !imus_.empty() && + this->getPose().r11() == 1.0f && this->getPose().r22() == 1.0f && this->getPose().r33() == 1.0f) + { + // Interpolated at the stamp, or the closest sample if the IMU doesn't cover it (e.g., + // the only sample received is just after the frame) + Transform imuAtStamp = Transform::getTransform(imus_, data.stamp()); + if(imuAtStamp.isNull()) + { + imuAtStamp = data.stamp() < imus_.begin()->first?imus_.begin()->second:imus_.rbegin()->second; + } + if(!imuAtStamp.isNull()) + { + const Eigen::Quaterniond q = imuAtStamp.getQuaterniond(); + const Transform previous = this->getPose(); + const Transform initialPose(previous.x(), previous.y(), previous.z(), q.x(), q.y(), q.z(), q.w()); + UWARN("Updated initial pose from %s to %s with IMU orientation", previous.prettyPrint().c_str(), initialPose.prettyPrint().c_str()); + std::map imus = imus_; + ImuMotionPredictor imuMotionPredictor = imuMotionPredictor_; + this->reset(initialPose); + imus_ = imus; + imuMotionPredictor_ = imuMotionPredictor; + } + } + // KITTI datasets start with stamp=0 double dt = previousStamp_>0.0f || (previousStamp_==0.0f && framesProcessed()==1)?data.stamp() - previousStamp_:0.0; Transform guess = dt>0.0 && guessFromMotion_ && !velocityGuess_.isNull()?Transform::getIdentity():Transform(); if(!(dt>0.0 || (dt == 0.0 && velocityGuess_.isNull()))) { - if(guessFromMotion_ && (!data.imageRaw().empty() || !data.laserScanRaw().isEmpty())) + if(guessFromMotion_) { UERROR("Guess from motion is set but dt is invalid! Odometry is then computed without guess. (dt=%f previous transform=%s)", dt, velocityGuess_.prettyPrint().c_str()); } @@ -678,7 +713,7 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet UWARN("Could not find imu transform at %f", data.stamp()); } } - else if(!guess.isNull() && (!data.imageRaw().empty() || !data.laserScanRaw().isEmpty())) { + else if(!guess.isNull()) { UDEBUG("Using guess from motion %s", guess.prettyPrint().c_str()); } @@ -859,23 +894,11 @@ Transform Odometry::process(SensorData & data, const Transform & guessIn, Odomet } } } - // A frame that brings its own features carries no image, and a frame whose scene was - // empty carries no feature either, so neither says whether there is a frame at all. - // The calibration does: it is there when a camera produced this data. - else if(!data.imageRaw().empty() || - !data.cameraModels().empty() || - !data.stereoCameraModels().empty() || - !data.laserScanRaw().isEmpty() || - (this->canProcessAsyncIMU() && !data.imu().empty())) + else { t = this->computeTransform(data, guess, info); } - if(data.imageRaw().empty() && data.laserScanRaw().isEmpty() && !data.imu().empty()) - { - return Transform(); // Return null on IMU-only updates - } - if(info) { info->timeEstimation = time.ticks(); diff --git a/corelib/test/test_odometry.cpp b/corelib/test/test_odometry.cpp index d1efe483..307974cf 100644 --- a/corelib/test/test_odometry.cpp +++ b/corelib/test/test_odometry.cpp @@ -11,6 +11,10 @@ #include #include #include +#include +#include +#include +#include #include using namespace rtabmap; @@ -584,3 +588,208 @@ TEST(OdometryTest, RefusesAFirstScanTooSmallForTheCorrespondenceRatio) EXPECT_FALSE(unchecked->process(uncheckedData).isNull()) << "a scan of unknown sweep size was refused"; } + +// --------------------------------------------------------------------------- +// Lidar deskewing with an IMU (Odom/Deskewing): a sensor moving in a box-shaped room, +// whose scans are generated point by point from where the sensor was at each point's +// time, as a spinning lidar measures them. +// --------------------------------------------------------------------------- + +namespace { + +const Eigen::Vector3d kRoomMin(-6.0, -4.0, -1.5); +const Eigen::Vector3d kRoomMax(7.0, 5.0, 2.5); +const double kSweep = 0.1; // s, first column to last +const int kRings = 16; +const int kColumns = 512; +const double kGravity = 9.80665; + +/// Pose of the sensor at time t, and its acceleration (world frame). +struct Trajectory +{ + double yaw = 0.0; // rad, heading at t=0 + double yawRate = 0.0; // rad/s + double acceleration = 0.0; // m/s^2 along x, from rest at t=0 + Transform pose(double t) const + { + return Transform(float(0.5*acceleration*t*t), 0, 0, 0, 0, float(yaw + yawRate*t)); + } + Eigen::Vector3d linearAcceleration() const { return Eigen::Vector3d(acceleration, 0, 0); } +}; + +/// The scan of a sweep starting at @p stamp: organized (rings x columns), time on columns. +LaserScan makeSweep(const Trajectory & trajectory, double stamp) +{ + cv::Mat data(kRings, kColumns, CV_32FC(5)); + for(int u=0; u::max(); + for(int i=0; i<3; ++i) + { + if(world[i] > 1e-9) range = std::min(range, (kRoomMax[i] - origin[i]) / world[i]); + else if(world[i] < -1e-9) range = std::min(range, (kRoomMin[i] - origin[i]) / world[i]); + } + float * p = data.ptr(v, u); + p[0] = float(range*direction.x()); + p[1] = float(range*direction.y()); + p[2] = float(range*direction.z()); + p[3] = 1.0f; + p[4] = float(dt); + } + } + return LaserScan(data, kRings*kColumns, 0.0f, LaserScan::kXYZIT, Transform::getIdentity()); +} + +/// What the IMU, at the sensor's origin, measures at time t. +IMU makeImu(const Trajectory & trajectory, double t) +{ + const Eigen::Quaterniond q = trajectory.pose(t).getQuaterniond(); + const Eigen::Vector3d f = q.inverse() * (trajectory.linearAcceleration() + Eigen::Vector3d(0, 0, kGravity)); + return IMU(cv::Vec4d(q.x(), q.y(), q.z(), q.w()), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(0, 0, trajectory.yawRate), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(f.x(), f.y(), f.z()), cv::Mat::eye(3,3,CV_64FC1), + Transform::getIdentity()); +} + +/// RMS distance (m) of the scan's points to the nearest wall, with the sensor at @p pose. +double wallDistance(const LaserScan & scan, const Transform & pose) +{ + double sum = 0.0; + int n = 0; + for(int i=0; i(0, i); + const Transform world = pose * Transform(p[0], p[1], p[2], 0, 0, 0); + const Eigen::Vector3d w(world.x(), world.y(), world.z()); + double d = std::numeric_limits::max(); + for(int k=0; k<3; ++k) + { + d = std::min(d, std::fabs(w[k] - kRoomMin[k])); + d = std::min(d, std::fabs(w[k] - kRoomMax[k])); + } + sum += d*d; + ++n; + } + return n?std::sqrt(sum/n):0.0; +} + +/** + * Feeds @p frames sweeps, one every 0.1 s from t=0.1, with the IMU at 200 Hz up to the end + * of each sweep (as OdometryThread does), and returns the RMS wall distance of the last + * sweep as odometry left it in the frame (deskewed or not). + */ +double deskewedWallDistance(const Trajectory & trajectory, int frames, const ParametersMap & extra) +{ + ParametersMap parameters; + parameters.insert(ParametersPair(Parameters::kRegStrategy(), "1")); + parameters.insert(ParametersPair(Parameters::kIcpVoxelSize(), "0")); + parameters.insert(ParametersPair(Parameters::kIcpPointToPlane(), "true")); + parameters.insert(ParametersPair(Parameters::kIcpPointToPlaneK(), "10")); + parameters.insert(ParametersPair(Parameters::kIcpMaxTranslation(), "0")); + parameters.insert(ParametersPair(Parameters::kIcpMaxCorrespondenceDistance(), "0.5")); + for(const ParametersPair & p : extra) uInsert(parameters, p); + std::unique_ptr odometry(Odometry::create(parameters)); + + double imuStamp = 0.0; + double distance = -1.0; + for(int i=1; i<=frames; ++i) + { + const double stamp = 0.1*i; + for(; imuStamp <= stamp + kSweep + 0.005; imuStamp += 0.005) + { + SensorData imu(makeImu(trajectory, imuStamp), 0, imuStamp); + odometry->process(imu); + } + SensorData data(makeSweep(trajectory, stamp), cv::Mat(), cv::Mat(), CameraModel(), i, stamp); + OdometryInfo info; + const Transform pose = odometry->process(data, &info); + EXPECT_FALSE(pose.isNull()) << "frame " << i << " not registered"; + distance = wallDistance(data.laserScanRaw(), trajectory.pose(stamp)); + } + return distance; +} + +} // namespace + +TEST(OdometryTest, DeskewsWithTheImuOrientation) +{ + // Turning at 2 rad/s: 11 degrees over a sweep, so the far walls are smeared by about a + // metre. On the first frame, before any pose, only the IMU orientation is used. + Trajectory turning; + turning.yawRate = 2.0; + ParametersMap off; + off.insert(ParametersPair(Parameters::kOdomDeskewing(), "false")); + const double skewed = deskewedWallDistance(turning, 1, off); + const double deskewed = deskewedWallDistance(turning, 1, ParametersMap()); + EXPECT_GT(skewed, 0.1) << "the scan should be smeared without deskewing"; + EXPECT_LT(deskewed, 0.01) << "the points should be back on the walls"; +} + +TEST(OdometryTest, DeskewsWithTheImuAccelerationOverTheSmoothingDelay) +{ + // Accelerating at 2 m/s^2 from rest. With Odom/GuessSmoothingDelay, the velocity is + // carried to each frame with the IMU acceleration, and the prediction integrates it + // over the sweep. At 0, the acceleration is not used: the velocity of the last + // interval lags behind, and the sweep is still skewed. + Trajectory accelerating; + accelerating.acceleration = 2.0; + // Not facing exactly along x: an IMU orientation whose x, y and z are all zero (the + // identity) is read by RTAB-Map's odometry as "not set", and the IMU ignored. + accelerating.yaw = 0.3; + ParametersMap noDelay; + noDelay.insert(ParametersPair(Parameters::kOdomGuessSmoothingDelay(), "0")); + ParametersMap delay; + delay.insert(ParametersPair(Parameters::kOdomGuessSmoothingDelay(), "0.5")); + const double withoutAcceleration = deskewedWallDistance(accelerating, 12, noDelay); + const double withAcceleration = deskewedWallDistance(accelerating, 12, delay); + EXPECT_LT(withAcceleration, 0.005) << "the points should be back on the walls"; + EXPECT_GT(withoutAcceleration, withAcceleration * 3.0); +} + +TEST(OdometryTest, InitialOrientationIsTheImuOrientationAtTheFirstFrame) +{ + // IMU received for 2 s before the first frame, while the sensor turns at 1 rad/s: the + // first frame must start from the orientation at its own stamp, not at the first IMU + // sample (2 rad earlier). + Trajectory turning; + turning.yaw = 0.3; + turning.yawRate = 1.0; + ParametersMap parameters; + parameters.insert(ParametersPair(Parameters::kRegStrategy(), "1")); + parameters.insert(ParametersPair(Parameters::kIcpVoxelSize(), "0")); + std::unique_ptr odometry(Odometry::create(parameters)); + const double stamp = 2.0; + for(double t = 0.0; t <= stamp + kSweep + 0.005; t += 0.005) + { + SensorData imu(makeImu(turning, t), 0, t); + odometry->process(imu); + } + SensorData data(makeSweep(turning, stamp), cv::Mat(), cv::Mat(), CameraModel(), 1, stamp); + OdometryInfo info; + const Transform pose = odometry->process(data, &info); + ASSERT_FALSE(pose.isNull()); + // The IMU was fed up to the end of the sweep (as OdometryThread does), but the orientation + // is interpolated at the frame's stamp + EXPECT_NEAR(pose.theta(), turning.yaw + turning.yawRate*stamp, 0.002); + + // An initial pose given with a rotation is kept + std::unique_ptr given(Odometry::create(parameters)); + given->reset(Transform(0, 0, 0, 0, 0, 1.0f)); + for(double t = 0.0; t <= stamp; t += 0.005) + { + SensorData imu(makeImu(turning, t), 0, t); + given->process(imu); + } + EXPECT_NEAR(given->getPose().theta(), 1.0, 1e-6); +} From 3357599379839052647412b24e21e58f27b99b5a Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sun, 11 Oct 2026 10:45:21 -0700 Subject: [PATCH 4/7] fixing sign issue with x86-64-v3 --- corelib/test/test_imumotionpredictor.cpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/corelib/test/test_imumotionpredictor.cpp b/corelib/test/test_imumotionpredictor.cpp index 5aff410d..b60b42a6 100644 --- a/corelib/test/test_imumotionpredictor.cpp +++ b/corelib/test/test_imumotionpredictor.cpp @@ -37,15 +37,19 @@ rtabmap::IMU measuredImu(const Eigen::Quaterniond & worldToImu, /// IMU measurements at 200 Hz over [from, to] of an IMU at the base origin, oriented and /// accelerating (in its world frame) as given. Without accelerometer, it measures no -/// linear acceleration at all. +/// linear acceleration at all. The stamps are on the 200 Hz grid (from must be too). void addSamples(ImuMotionPredictor & predictor, double from, double to, const std::function & orientation, const std::function & acceleration, bool withAccelerometer = true) { - for(int i=0; from + i*0.005 <= to + 1e-9; ++i) + // An integer over 200 rather than from + i*0.005: with FMA (e.g., -march=x86-64-v3), + // -0.05 + 10*0.005 gives -1.7e-18, not 0, putting that sample on the wrong side of + // a step at t=0. + const double first = std::round(from * 200.0); + for(int i=0; (first + i) / 200.0 <= to + 1e-9; ++i) { - const double t = from + i*0.005; + const double t = (first + i) / 200.0; rtabmap::IMU imu = measuredImu(orientation(t), acceleration(t), predictor.gravity(), rtabmap::Transform::getIdentity()); if(!withAccelerometer) { From 5d447d9bb0a242da61350400aceaba123778786d Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sun, 11 Oct 2026 10:49:07 -0700 Subject: [PATCH 5/7] fixing M_PI not defined in tests --- corelib/test/test_imumotionpredictor.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/corelib/test/test_imumotionpredictor.cpp b/corelib/test/test_imumotionpredictor.cpp index b60b42a6..48c61cc0 100644 --- a/corelib/test/test_imumotionpredictor.cpp +++ b/corelib/test/test_imumotionpredictor.cpp @@ -1,3 +1,7 @@ +// M_PI on MSVC (must come before any header including ) +#ifndef _USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#endif #include #include From 526c4226436279d641005f2565ec39631be1126b Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sun, 11 Oct 2026 11:36:44 -0700 Subject: [PATCH 6/7] Added auto gravity estimation mode --- .../include/rtabmap/core/ImuMotionPredictor.h | 50 ++++++- corelib/include/rtabmap/core/Parameters.h | 2 +- corelib/src/ImuMotionPredictor.cpp | 130 ++++++++++++++++-- corelib/test/test_imumotionpredictor.cpp | 83 +++++++++++ guilib/src/ui/preferencesDialog.ui | 5 +- 5 files changed, 249 insertions(+), 21 deletions(-) diff --git a/corelib/include/rtabmap/core/ImuMotionPredictor.h b/corelib/include/rtabmap/core/ImuMotionPredictor.h index 7517f261..20f355b9 100644 --- a/corelib/include/rtabmap/core/ImuMotionPredictor.h +++ b/corelib/include/rtabmap/core/ImuMotionPredictor.h @@ -35,6 +35,7 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. #include #include +#include namespace rtabmap { @@ -60,6 +61,17 @@ namespace rtabmap { * ignored: its centripetal and tangential accelerations are small over the fraction of a * second this predicts. * + * The gravity removed can be estimated from the IMU itself (see the constructor): an + * accelerometer can read a few percent off at rest (scale or bias), and that error, + * integrated as an acceleration, would bias the predicted velocity along gravity. The + * accelerometer is averaged over the windows where it is quiet: with a gyro and an + * accelerometer that barely change, the IMU is not accelerating (at rest or at constant + * velocity), so it only measures gravity. The IMU must then be held still (or at constant + * velocity) for a moment, e.g., before moving: until then, standard gravity is used, and + * a steady acceleration without rotation would be taken for gravity. Only the error along + * gravity is corrected, which is all of it for an IMU that stays about level; one changing + * attitude would need its axes calibrated. + * * Not thread-safe: a caller sharing it between threads must lock around every call. */ class RTABMAP_CORE_EXPORT ImuMotionPredictor @@ -81,7 +93,9 @@ public: * velocity itself. 0: the displacement since the previous pose, * without acceleration. * @param gravity magnitude (m/s^2) of the gravity removed from the specific force - * given to addImu(), standard gravity by default + * given to addImu(), standard gravity by default. <= 0: estimated + * from the IMU when it is still (see the class description), with + * standard gravity until then. * * These are fixed for the life of the predictor (there is no setter): changing them * while it estimates would mix poses and samples taken under different settings. @@ -93,7 +107,12 @@ public: double maxPoseInterval() const {return maxPoseInterval_;} double velocityWindow() const {return velocityWindow_;} + /// Gravity (m/s^2) removed from the specific force: the one given, or the estimate double gravity() const {return gravity_;} + /// Whether the gravity is estimated from the IMU (see the constructor) + bool isGravityEstimated() const {return gravityEstimated_;} + /// Number of quiet IMU windows the gravity estimate is from (0: not estimated yet) + size_t gravityWindows() const {return gravityWindowsCount_;} /** * @brief Adds an IMU measurement. @@ -115,7 +134,8 @@ public: */ void addPose(double stamp, const rtabmap::Transform & pose); - /// Forgets the poses and the IMU samples. + /// Forgets the poses and the IMU samples. The gravity estimate is kept: it is the + /// sensor's. void reset(); /** @@ -140,13 +160,22 @@ public: private: // orientation: of the base frame in the IMU's world frame; acceleration: of the base in // that frame, gravity removed - void addSample(double stamp, const Eigen::Quaterniond & orientation, const Eigen::Vector3d & acceleration); + // specificForce: in that frame too, including the reaction to gravity (see + // hasAcceleration) + void addSample(double stamp, const Eigen::Quaterniond & orientation, + const Eigen::Vector3d & specificForce, bool hasAcceleration); + // Gravity estimation, from the IMU's own specific force and angular velocity norm + void estimateGravity(double stamp, const Eigen::Vector3d & specificForce, double angularVelocity); struct Sample { Eigen::Quaterniond orientation; - Eigen::Vector3d acceleration; + // Kept with gravity, which is removed when used: the estimate can change + Eigen::Vector3d specificForce; + bool hasAcceleration; }; + // Acceleration of a sample, gravity removed (null without acceleration) + Eigen::Vector3d accelerationOf(const Sample & sample) const; Eigen::Quaterniond orientationAt(double stamp) const; Eigen::Vector3d accelerationAt(double stamp) const; void integrate(double from, double to, const Eigen::Quaterniond & rotation, @@ -166,7 +195,20 @@ private: double maxPoseInterval_; double velocityWindow_; double gravity_; + bool gravityEstimated_; std::map samples_; + + // Gravity estimation: the samples of the current window, and the gravity measured over + // the last quiet windows + struct StillSample + { + double stamp; + Eigen::Vector3d specificForce; // in the IMU frame + double angularVelocity; // norm + }; + std::vector stillWindow_; + std::vector gravityWindows_; + size_t gravityWindowsCount_; // Recent odometry poses, to estimate the velocity from (see velocityWindow) std::map poses_; diff --git a/corelib/include/rtabmap/core/Parameters.h b/corelib/include/rtabmap/core/Parameters.h index cdd4d5fe..6a749928 100644 --- a/corelib/include/rtabmap/core/Parameters.h +++ b/corelib/include/rtabmap/core/Parameters.h @@ -512,7 +512,7 @@ class RTABMAP_CORE_EXPORT Parameters RTABMAP_PARAM(Odom, KalmanProcessNoise, float, 0.001, "Process noise covariance value."); RTABMAP_PARAM(Odom, KalmanMeasurementNoise, float, 0.01, "Process measurement covariance value."); RTABMAP_PARAM(Odom, GuessMotion, bool, true, "Guess next transformation from the last motion computed."); - RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration (used with \"%s\" > 0). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's.", kOdomGuessSmoothingDelay().c_str())); + RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration (used with \"%s\" > 0). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. 0: estimated from the IMU while it is still (gyro and accelerometer quiet over 0.5 s), with standard gravity until then: the robot should then be perfectly still for a moment (e.g., at start), and a steady acceleration without rotation would be taken for gravity.", kOdomGuessSmoothingDelay().c_str())); RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is then the displacement over this delay corrected by the acceleration measured since, so that it is the velocity at the last frame rather than a delayed average, and the motion guess and lidar deskewing (\"%s\") integrate the acceleration from it. Recommended (~0.5 s) with an IMU. With 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomDeskewing().c_str())); RTABMAP_PARAM(Odom, KeyFrameThr, float, 0.3, "[Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, VisKeyFrameThr, int, 150, "[Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame."); diff --git a/corelib/src/ImuMotionPredictor.cpp b/corelib/src/ImuMotionPredictor.cpp index ff12698d..15e293ab 100644 --- a/corelib/src/ImuMotionPredictor.cpp +++ b/corelib/src/ImuMotionPredictor.cpp @@ -26,7 +26,9 @@ SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ #include +#include +#include #include namespace rtabmap { @@ -35,14 +37,29 @@ namespace rtabmap { // bounded while no odometry pose comes to trim it. static const double kMaxBufferDuration = 10.0; +// Gravity estimation (see the class description). The IMU is still over a window if its +// angular velocity stays under kStillMaxAngularVelocity and the standard deviation of its +// specific force, on every axis, under kStillMaxAccelerationStd. Those are above the noise +// and vibrations of an IMU on a robot at rest (e.g., ~0.1 m/s^2 on a drone on the ground, +// against ~1 m/s^2 in flight), and well under what any acceleration that matters would +// cause. The gravity is the median over the last kGravityWindows quiet windows. +static const double kStillWindow = 0.5; // s +static const double kStillMaxAngularVelocity = 0.05; // rad/s +static const double kStillMaxAccelerationStd = 0.15; // m/s^2 +static const size_t kStillMinSamples = 5; +static const size_t kGravityWindows = 100; +static const double kStandardGravity = 9.80665; + ImuMotionPredictor::ImuMotionPredictor(double maxPoseInterval, double velocityWindow, double gravity) : maxPoseInterval_(maxPoseInterval), velocityWindow_(velocityWindow), - gravity_(gravity), - integratedStartIsFinal_(false), + gravity_(gravity > 0.0 ? gravity : kStandardGravity), + gravityEstimated_(gravity <= 0.0), + gravityWindowsCount_(0), poseStamp_(0.0), velocity_(Eigen::Vector3d::Zero()), - worldToOdom_(Eigen::Quaterniond::Identity()) + worldToOdom_(Eigen::Quaterniond::Identity()), + integratedStartIsFinal_(false) { } @@ -63,21 +80,93 @@ void ImuMotionPredictor::addImu(double stamp, const IMU & imu) imu.localTransform().getQuaterniond(); const cv::Vec3d & f = imu.linearAcceleration(); - Eigen::Vector3d acceleration = Eigen::Vector3d::Zero(); - if((f[0] != 0.0 || f[1] != 0.0 || f[2] != 0.0) && - (imu.linearAccelerationCovariance().empty() || imu.linearAccelerationCovariance().at(0,0) != -1.0)) + const bool hasAcceleration = (f[0] != 0.0 || f[1] != 0.0 || f[2] != 0.0) && + (imu.linearAccelerationCovariance().empty() || imu.linearAccelerationCovariance().at(0,0) != -1.0); + if(hasAcceleration && gravityEstimated_) { - // The specific force includes the reaction to gravity, up in the world frame - acceleration = worldToImu * Eigen::Vector3d(f[0], f[1], f[2]) - Eigen::Vector3d(0, 0, gravity_); + const cv::Vec3d & w = imu.angularVelocity(); + estimateGravity(stamp, Eigen::Vector3d(f[0], f[1], f[2]), Eigen::Vector3d(w[0], w[1], w[2]).norm()); } - addSample(stamp, worldToImu * baseToImu.inverse(), acceleration); + // The specific force includes the reaction to gravity, up in the world frame: it is + // removed when used (accelerationOf()) + addSample(stamp, worldToImu * baseToImu.inverse(), + hasAcceleration ? Eigen::Vector3d(worldToImu * Eigen::Vector3d(f[0], f[1], f[2])) : Eigen::Vector3d::Zero(), + hasAcceleration); } -void ImuMotionPredictor::addSample(double stamp, const Eigen::Quaterniond & orientation, const Eigen::Vector3d & acceleration) +void ImuMotionPredictor::estimateGravity(double stamp, const Eigen::Vector3d & specificForce, double angularVelocity) +{ + if(!stillWindow_.empty() && stamp < stillWindow_.back().stamp) + { + // Out of order + stillWindow_.clear(); + } + if(!stillWindow_.empty() && stamp - stillWindow_.front().stamp >= kStillWindow) + { + // The window is complete (windows don't overlap: each sample counts once) + if(stillWindow_.size() >= kStillMinSamples) + { + Eigen::Vector3d mean = Eigen::Vector3d::Zero(); + double maxAngularVelocity = 0.0; + for(size_t i=0; i kGravityWindows) + { + gravityWindows_.erase(gravityWindows_.begin()); + } + ++gravityWindowsCount_; + std::vector sorted = gravityWindows_; + std::nth_element(sorted.begin(), sorted.begin() + sorted.size()/2, sorted.end()); + const double gravity = sorted[sorted.size()/2]; + if(gravityWindowsCount_ == 1) + { + UINFO("IMU gravity estimated at %f m/s^2 (the IMU is still)", gravity); + } + else + { + UDEBUG("IMU gravity estimated at %f m/s^2 (%d still windows)", gravity, (int)gravityWindowsCount_); + } + if(gravity != gravity_) + { + gravity_ = gravity; + // The accelerations integrated changed + integrated_.clear(); + } + } + } + stillWindow_.clear(); + } + StillSample sample; + sample.stamp = stamp; + sample.specificForce = specificForce; + sample.angularVelocity = angularVelocity; + stillWindow_.push_back(sample); +} + +void ImuMotionPredictor::addSample(double stamp, const Eigen::Quaterniond & orientation, + const Eigen::Vector3d & specificForce, bool hasAcceleration) { Sample sample; sample.orientation = orientation.normalized(); - sample.acceleration = acceleration; + sample.specificForce = specificForce; + sample.hasAcceleration = hasAcceleration; samples_[stamp] = sample; if(!integrated_.empty() && stamp <= integrated_.rbegin()->first) { @@ -172,6 +261,7 @@ void ImuMotionPredictor::addPose(double stamp, const rtabmap::Transform & pose) void ImuMotionPredictor::reset() { samples_.clear(); + stillWindow_.clear(); poses_.clear(); integrated_.clear(); pose_.setNull(); @@ -287,7 +377,7 @@ void ImuMotionPredictor::updateIntegration() const Integrated next; // Same segment as in integrate(), from the previous sample ti to this one: // D(ti+1) = D(ti) + V(ti)*dt + dt^2*(2*a(ti) + a(ti+1))/6, V(ti+1) = V(ti) + dt*(a(ti) + a(ti+1))/2 - next.acceleration = worldToOdom_ * iter->second.acceleration; + next.acceleration = worldToOdom_ * accelerationOf(iter->second); next.position = previous.second.position + previous.second.velocity * dt + dt * dt * (2.0 * previous.second.acceleration + next.acceleration) / 6.0; next.velocity = previous.second.velocity + dt * (previous.second.acceleration + next.acceleration) / 2.0; @@ -311,20 +401,30 @@ Eigen::Quaterniond ImuMotionPredictor::orientationAt(double stamp) const return before->second.orientation.slerp(ratio, after->second.orientation); } +Eigen::Vector3d ImuMotionPredictor::accelerationOf(const Sample & sample) const +{ + if(!sample.hasAcceleration) + { + return Eigen::Vector3d::Zero(); + } + return sample.specificForce - Eigen::Vector3d(0, 0, gravity_); +} + Eigen::Vector3d ImuMotionPredictor::accelerationAt(double stamp) const { std::map::const_iterator after = samples_.lower_bound(stamp); if(after == samples_.end()) { - return samples_.rbegin()->second.acceleration; + return accelerationOf(samples_.rbegin()->second); } if(after == samples_.begin() || after->first == stamp) { - return after->second.acceleration; + return accelerationOf(after->second); } std::map::const_iterator before = std::prev(after); const double ratio = (stamp - before->first) / (after->first - before->first); - return before->second.acceleration + ratio * (after->second.acceleration - before->second.acceleration); + const Eigen::Vector3d accelerationBefore = accelerationOf(before->second); + return accelerationBefore + ratio * (accelerationOf(after->second) - accelerationBefore); } void ImuMotionPredictor::integrate(double from, double to, const Eigen::Quaterniond & rotation, diff --git a/corelib/test/test_imumotionpredictor.cpp b/corelib/test/test_imumotionpredictor.cpp index 48c61cc0..ba65ed9f 100644 --- a/corelib/test/test_imumotionpredictor.cpp +++ b/corelib/test/test_imumotionpredictor.cpp @@ -375,3 +375,86 @@ TEST(ImuMotionPredictor, ignores_the_acceleration_without_a_velocity_window) EXPECT_NEAR(predictor.velocity().x(), 0.5*a*0.1, 1e-6); EXPECT_NEAR(predictor.predict(0.2).x(), 0.5*a*0.01 + 0.5*a*0.1*0.1, 1e-6); } + +namespace { + +/// An IMU tilted by 0.1 rad in pitch whose accelerometer reads `measuredGravity` at rest, +/// with a vibration of the given amplitude (m/s^2) on every axis and the given angular +/// velocity (rad/s, about z). +rtabmap::IMU tiltedImu(double t, double measuredGravity, double vibration = 0.05, double angularVelocity = 0.0) +{ + const Eigen::Quaterniond worldToImu(Eigen::AngleAxisd(0.1, Eigen::Vector3d::UnitY())); + const Eigen::Vector3d f = worldToImu.inverse() * Eigen::Vector3d(0, 0, measuredGravity) + + vibration * Eigen::Vector3d(std::sin(t*331.0), std::sin(t*457.0), std::sin(t*563.0)); + return rtabmap::IMU(cv::Vec4d(worldToImu.x(), worldToImu.y(), worldToImu.z(), worldToImu.w()), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(0, 0, angularVelocity), cv::Mat::eye(3,3,CV_64FC1), + cv::Vec3d(f.x(), f.y(), f.z()), cv::Mat::eye(3,3,CV_64FC1), + rtabmap::Transform::getIdentity()); +} + +} // namespace + +TEST(ImuMotionPredictor, estimates_the_gravity_while_the_imu_is_still) +{ + // An accelerometer reading 9.55 m/s^2 at rest, as the netherdrone's does. + ImuMotionPredictor predictor(1.0, 0.5, 0.0); + EXPECT_TRUE(predictor.isGravityEstimated()); + EXPECT_DOUBLE_EQ(predictor.gravity(), 9.80665) << "standard gravity until estimated"; + + for(int i=0; i<100; ++i) // 0.5 s at 200 Hz: the first window is not complete yet + { + predictor.addImu(i*0.005, tiltedImu(i*0.005, 9.55)); + } + EXPECT_EQ(predictor.gravityWindows(), 0u); + EXPECT_DOUBLE_EQ(predictor.gravity(), 9.80665); + + for(int i=100; i<=400; ++i) + { + predictor.addImu(i*0.005, tiltedImu(i*0.005, 9.55)); + } + EXPECT_GE(predictor.gravityWindows(), 3u); + EXPECT_NEAR(predictor.gravity(), 9.55, 0.005) << "whatever the attitude"; + + // At rest, it doesn't look like falling anymore. + predictor.addPose(1.5, rtabmap::Transform::getIdentity()); + predictor.addPose(2.0, rtabmap::Transform::getIdentity()); + EXPECT_NEAR(predictor.predict(2.0).z(), 0.0, 1e-3); + for(int i=401; i<=460; ++i) + { + predictor.addImu(i*0.005, tiltedImu(i*0.005, 9.55)); + } + EXPECT_NEAR(predictor.predict(2.3).z(), 0.0, 2e-3); + + // The estimate is the sensor's: a reset keeps it. + predictor.reset(); + EXPECT_NEAR(predictor.gravity(), 9.55, 0.005); +} + +TEST(ImuMotionPredictor, does_not_take_a_moving_imu_for_a_still_one) +{ + // Turning (centripetal acceleration, gravity moving between the axes) or vibrating + // as in flight: not still, standard gravity is kept. + ImuMotionPredictor turning(1.0, 0.5, 0.0); + ImuMotionPredictor vibrating(1.0, 0.5, 0.0); + for(int i=0; i<=400; ++i) + { + turning.addImu(i*0.005, tiltedImu(i*0.005, 9.55, 0.05, 0.5)); + vibrating.addImu(i*0.005, tiltedImu(i*0.005, 9.55, 1.0)); + } + EXPECT_EQ(turning.gravityWindows(), 0u); + EXPECT_DOUBLE_EQ(turning.gravity(), 9.80665); + EXPECT_EQ(vibrating.gravityWindows(), 0u); + EXPECT_DOUBLE_EQ(vibrating.gravity(), 9.80665); +} + +TEST(ImuMotionPredictor, keeps_the_gravity_it_is_given) +{ + ImuMotionPredictor predictor; + for(int i=0; i<=400; ++i) + { + predictor.addImu(i*0.005, tiltedImu(i*0.005, 9.55)); + } + EXPECT_FALSE(predictor.isGravityEstimated()); + EXPECT_EQ(predictor.gravityWindows(), 0u); + EXPECT_DOUBLE_EQ(predictor.gravity(), 9.80665); +} diff --git a/guilib/src/ui/preferencesDialog.ui b/guilib/src/ui/preferencesDialog.ui index 6a30578e..8406e50d 100644 --- a/guilib/src/ui/preferencesDialog.ui +++ b/guilib/src/ui/preferencesDialog.ui @@ -16866,6 +16866,9 @@ With <0, the length is estimated once for each unique marker, then re-used fo + + Auto + m/s² @@ -16886,7 +16889,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - Gravity magnitude removed from the IMU linear acceleration (used with a guess smoothing delay > 0 above). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. + Gravity magnitude removed from the IMU linear acceleration (used with a guess smoothing delay > 0 above). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. Auto (0): estimated from the IMU while it is still, with standard gravity until then: the robot should then be perfectly still for a moment (e.g., at start). true From 05b5661e6b7c37e5d8921cd61b16df1cd02d6789 Mon Sep 17 00:00:00 2001 From: matlabbe Date: Sun, 11 Oct 2026 11:45:52 -0700 Subject: [PATCH 7/7] Fixing crash on windows when uFormat() is to long --- corelib/include/rtabmap/core/Parameters.h | 4 ++-- corelib/test/CMakeLists.txt | 1 + corelib/test/test_parameters.cpp | 20 ++++++++++++++++++++ guilib/src/ui/preferencesDialog.ui | 2 +- utilite/src/UConversion.cpp | 3 ++- utilite/test/test_uconversion.cpp | 11 +++++++++++ 6 files changed, 37 insertions(+), 4 deletions(-) create mode 100644 corelib/test/test_parameters.cpp diff --git a/corelib/include/rtabmap/core/Parameters.h b/corelib/include/rtabmap/core/Parameters.h index 6a749928..56328910 100644 --- a/corelib/include/rtabmap/core/Parameters.h +++ b/corelib/include/rtabmap/core/Parameters.h @@ -489,7 +489,7 @@ class RTABMAP_CORE_EXPORT Parameters RTABMAP_PARAM(Optimizer, Baseline, double, 0.075, "When doing bundle adjustment with RGB-D data (mono camera + depth), set a fake baseline (m) so the BA backend treats depth as stereo disparity. Applies to all BA-capable backends (g2o, GTSAM, Ceres). Set to 0 to keep the problem mono (depth observations are ignored). For real stereo data the baseline in the calibration (Tx) is used directly."); RTABMAP_PARAM(Optimizer, PixelVariance, double, 1.0, "Pixel variance used on the u/v axes of every bundle adjustment reprojection edge. Applies to all BA-capable backends (g2o, GTSAM, Ceres). Should approximate the squared 1-sigma keypoint localization error in pixels. Set higher (e.g. 4-9) if features are noisy (low texture, motion blur, low light, or large detector scale). Set lower (e.g. 0.01-0.1) if features are sub-pixel refined (Lucas-Kanade tracking, parabolic peak interpolation). Intuition: the lower the pixel variance, the more the optimizer trusts the keypoint positions."); - RTABMAP_PARAM(Optimizer, DisparityVariance, double, 1.0, "Disparity variance used on the disparity axis (u - u_right) of stereo / RGB-D bundle adjustment edges. Applies to all BA-capable backends (g2o, GTSAM, Ceres). Defaults to the same value as PixelVariance for backward compatibility. Set higher (e.g. 2-4) if your depth source is noisier than your feature detector's u/v precision (typical for stereo block matchers / SGM at long range). Set lower (e.g. 0.01-0.1) if your depth source is more accurate than the u/v detector (typical for ToF / LiDAR-fused depth where range is measured directly rather than triangulated). Intuition: the lower the disparity variance, the more the optimizer trusts the depth measurements. Geometric note: wider baseline and/or higher image resolution improve a block matcher's effective disparity precision (larger disparity magnitudes and finer sub-pixel refinement), so wide-baseline high-resolution stereo pairs can usually afford a lower disparity variance (e.g. 0.1-0.5); narrow-baseline low-resolution pairs should keep it higher (e.g. 1-4)."); + RTABMAP_PARAM(Optimizer, DisparityVariance, double, 1.0, "Disparity variance used on the disparity axis (u - u_right) of stereo / RGB-D bundle adjustment edges (g2o, GTSAM, Ceres). Defaults to the same value as PixelVariance for backward compatibility. The lower it is, the more the optimizer trusts the depth. Set higher (e.g. 2-4) if the depth is noisier than the features' u/v precision (stereo block matchers / SGM at long range, narrow-baseline low-resolution stereo). Set lower (e.g. 0.01-0.1) if it is more accurate (ToF / LiDAR-fused depth, measured rather than triangulated), or around 0.1-0.5 for wide-baseline high-resolution stereo."); RTABMAP_PARAM(Optimizer, RobustKernelDelta, double, 8, "Robust kernel delta used for bundle adjustment (0 means don't use robust kernel). Applies to all BA-capable backends (g2o, GTSAM, Ceres). Observations with chi2 over this threshold will be ignored in the second optimization pass."); RTABMAP_PARAM(GTSAM, Optimizer, int, 1, "0=Levenberg 1=GaussNewton 2=Dogleg"); @@ -513,7 +513,7 @@ class RTABMAP_CORE_EXPORT Parameters RTABMAP_PARAM(Odom, KalmanMeasurementNoise, float, 0.01, "Process measurement covariance value."); RTABMAP_PARAM(Odom, GuessMotion, bool, true, "Guess next transformation from the last motion computed."); RTABMAP_PARAM(Odom, ImuGravity, float, 9.80665, uFormat("Gravity magnitude (m/s^2) removed from the IMU linear acceleration (used with \"%s\" > 0). Standard gravity by default. Set it to what the accelerometer reads at rest to compensate for its scale error, or for another gravity than Earth's. 0: estimated from the IMU while it is still (gyro and accelerometer quiet over 0.5 s), with standard gravity until then: the robot should then be perfectly still for a moment (e.g., at start), and a steady acceleration without rotation would be taken for gravity.", kOdomGuessSmoothingDelay().c_str())); - RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if \"%s\" is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is then the displacement over this delay corrected by the acceleration measured since, so that it is the velocity at the last frame rather than a delayed average, and the motion guess and lidar deskewing (\"%s\") integrate the acceleration from it. Recommended (~0.5 s) with an IMU. With 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomDeskewing().c_str())); + RTABMAP_PARAM(Odom, GuessSmoothingDelay, float, 0, uFormat("Guess smoothing delay (s). The velocity is averaged over the last transforms up to this delay, for a smoother velocity prediction. The last velocity is used directly if \"%s\" is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is the displacement over this delay corrected by the acceleration measured since, so it is not delayed, and the motion guess and lidar deskewing (\"%s\") integrate the acceleration from it (~0.5 s recommended). With 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: useful for platforms with inertia, high frame rates, or lidar deskewing, where the velocity noise would feed back into the next poses.", kOdomFilteringStrategy().c_str(), kOdomDeskewing().c_str())); RTABMAP_PARAM(Odom, KeyFrameThr, float, 0.3, "[Visual] Create a new keyframe when the number of inliers drops under this ratio of features in last frame. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, VisKeyFrameThr, int, 150, "[Visual] Create a new keyframe when the number of inliers drops under this threshold. Setting the value to 0 means that a keyframe is created for each processed frame."); RTABMAP_PARAM(Odom, ScanKeyFrameThr, float, 0.9, "[Geometry] Create a new keyframe when the number of ICP inliers drops under this ratio of points in last frame's scan. Setting the value to 0 means that a keyframe is created for each processed frame."); diff --git a/corelib/test/CMakeLists.txt b/corelib/test/CMakeLists.txt index 049b61c6..a1af10df 100644 --- a/corelib/test/CMakeLists.txt +++ b/corelib/test/CMakeLists.txt @@ -42,6 +42,7 @@ set(corelib_test_sources test_link.cpp #Link.h test_optimizer.cpp #Optimizer.h test_gps.cpp #GPS.h + test_parameters.cpp #Parameters.h test_imu.cpp #IMU.h test_imufilter.cpp #IMUFilter.h test_imumotionpredictor.cpp #ImuMotionPredictor.h diff --git a/corelib/test/test_parameters.cpp b/corelib/test/test_parameters.cpp new file mode 100644 index 00000000..20701090 --- /dev/null +++ b/corelib/test/test_parameters.cpp @@ -0,0 +1,20 @@ +#include +#include + +using rtabmap::Parameters; +using rtabmap::ParametersMap; + +TEST(Parameters, descriptions_are_shorter_than_1024_characters) +{ + // Long descriptions don't fit in the GUI's tooltips nor in a ROS parameter listing: + // keep them short, the details belong in the documentation. 1024 was also the size + // beyond which uFormat(), which formats most of them, aborted on Windows. + const ParametersMap & parameters = Parameters::getDefaultParameters(); + ASSERT_FALSE(parameters.empty()); + for(ParametersMap::const_iterator iter = parameters.begin(); iter != parameters.end(); ++iter) + { + const std::string description = Parameters::getDescription(iter->first); + EXPECT_LT(description.size(), 1024u) << iter->first << " has a description of " + << description.size() << " characters: shorten it"; + } +} diff --git a/guilib/src/ui/preferencesDialog.ui b/guilib/src/ui/preferencesDialog.ui index 8406e50d..30a16282 100644 --- a/guilib/src/ui/preferencesDialog.ui +++ b/guilib/src/ui/preferencesDialog.ui @@ -16794,7 +16794,7 @@ With <0, the length is estimated once for each unique marker, then re-used fo - Guess smoothing delay (s). Estimated velocity is averaged based on last transforms up to this maximum delay. This can help to get smoother velocity prediction. Last velocity computed is used directly if a filtering strategy is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is then corrected by the acceleration measured since, so that it is the velocity at the last frame rather than a delayed average, and the motion guess and lidar deskewing integrate the acceleration from it. Recommended (~0.5 s) with an IMU; with 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: it can make sense for platforms with inertia (e.g., cars at high speed, where one bad frame would otherwise change the predicted velocity a lot), for high frame rates, or when the velocity is used for lidar deskewing, where its noise would feed back into the next poses. + Guess smoothing delay (s). The velocity is averaged over the last transforms up to this delay, for a smoother velocity prediction. The last velocity is used directly if a filtering strategy is set or the delay is below the odometry rate. With an IMU giving orientation and linear acceleration, a delay > 0 also enables the IMU acceleration: the velocity is the displacement over this delay corrected by the acceleration measured since, so it is not delayed, and the motion guess and lidar deskewing integrate the acceleration from it (~0.5 s recommended). With 0, the IMU only gives the orientation. Without IMU, the average is delayed by half this delay: useful for platforms with inertia, high frame rates, or lidar deskewing, where the velocity noise would feed back into the next poses. true diff --git a/utilite/src/UConversion.cpp b/utilite/src/UConversion.cpp index 200024c0..c78d5632 100644 --- a/utilite/src/UConversion.cpp +++ b/utilite/src/UConversion.cpp @@ -338,7 +338,8 @@ std::string uFormatv (const char *fmt, va_list args) // Try to vsnprintf into our buffer. #ifdef _MSC_VER - int needed = vsnprintf_s(buf, size, size, fmt, argsTmp); + // _TRUNCATE: return -1 if it doesn't fit (count=size would call the invalid parameter handler and abort) + int needed = vsnprintf_s(buf, size, _TRUNCATE, fmt, argsTmp); #else int needed = vsnprintf (buf, size, fmt, argsTmp); #endif diff --git a/utilite/test/test_uconversion.cpp b/utilite/test/test_uconversion.cpp index ca5a5c03..8218651f 100644 --- a/utilite/test/test_uconversion.cpp +++ b/utilite/test/test_uconversion.cpp @@ -201,3 +201,14 @@ TEST(UConversionTest, UFormat) EXPECT_NE(result.find("42"), std::string::npos); } +TEST(UConversionTest, UFormatLongerThanItsFirstBuffer) +{ + // Over the 1024 bytes uFormat first tries: formatted whole, not truncated (MSVC's + // checked vsnprintf_s used to abort the process instead). + const std::string longText(3000, 'a'); + std::string result = uFormat("%s %d", longText.c_str(), 42); + EXPECT_EQ(result, longText + " 42"); + result = uFormat("%s", std::string(1023, 'b').c_str()); + EXPECT_EQ(result, std::string(1023, 'b')); +} +