|
@@ -0,0 +1,551 @@
|
|
|
+/*
|
|
|
+ @author Roberto G. Valenti <robertogl.valenti@gmail.com>
|
|
|
+
|
|
|
+ @section LICENSE
|
|
|
+ Copyright (c) 2015, City University of New York
|
|
|
+ CCNY Robotics Lab <http://robotics.ccny.cuny.edu>
|
|
|
+ All rights reserved.
|
|
|
+
|
|
|
+ Redistribution and use in source and binary forms, with or without
|
|
|
+ modification, are permitted provided that the following conditions are met:
|
|
|
+ 1. Redistributions of source code must retain the above copyright
|
|
|
+ notice, this list of conditions and the following disclaimer.
|
|
|
+ 2. 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.
|
|
|
+ 3. Neither the name of the City College of New York 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 CCNY ROBOTICS LAB 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 "imu_complementary_filter/complementary_filter.h"
|
|
|
+
|
|
|
+#include <cstdio>
|
|
|
+#include <cmath>
|
|
|
+#include <iostream>
|
|
|
+
|
|
|
+namespace imu_tools {
|
|
|
+
|
|
|
+const double ComplementaryFilter::kGravity = 9.81;
|
|
|
+const double ComplementaryFilter::gamma_ = 0.01;
|
|
|
+// Bias estimation steady state thresholds
|
|
|
+const double ComplementaryFilter::kAngularVelocityThreshold = 0.2;
|
|
|
+const double ComplementaryFilter::kAccelerationThreshold = 0.1;
|
|
|
+const double ComplementaryFilter::kDeltaAngularVelocityThreshold = 0.01;
|
|
|
+
|
|
|
+ComplementaryFilter::ComplementaryFilter() :
|
|
|
+ gain_acc_(0.01),
|
|
|
+ gain_mag_(0.01),
|
|
|
+ bias_alpha_(0.01),
|
|
|
+ do_bias_estimation_(true),
|
|
|
+ do_adaptive_gain_(false),
|
|
|
+ initialized_(false),
|
|
|
+ steady_state_(false),
|
|
|
+ q0_(1), q1_(0), q2_(0), q3_(0),
|
|
|
+ wx_prev_(0), wy_prev_(0), wz_prev_(0),
|
|
|
+ wx_bias_(0), wy_bias_(0), wz_bias_(0) { }
|
|
|
+
|
|
|
+ComplementaryFilter::~ComplementaryFilter() { }
|
|
|
+
|
|
|
+void ComplementaryFilter::setDoBiasEstimation(bool do_bias_estimation)
|
|
|
+{
|
|
|
+ do_bias_estimation_ = do_bias_estimation;
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::getDoBiasEstimation() const
|
|
|
+{
|
|
|
+ return do_bias_estimation_;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::setDoAdaptiveGain(bool do_adaptive_gain)
|
|
|
+{
|
|
|
+ do_adaptive_gain_ = do_adaptive_gain;
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::getDoAdaptiveGain() const
|
|
|
+{
|
|
|
+ return do_adaptive_gain_;
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::setGainAcc(double gain)
|
|
|
+{
|
|
|
+ if (gain >= 0 && gain <= 1.0)
|
|
|
+ {
|
|
|
+ gain_acc_ = gain;
|
|
|
+ return true;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ return false;
|
|
|
+}
|
|
|
+bool ComplementaryFilter::setGainMag(double gain)
|
|
|
+{
|
|
|
+ if (gain >= 0 && gain <= 1.0)
|
|
|
+ {
|
|
|
+ gain_mag_ = gain;
|
|
|
+ return true;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ return false;
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getGainAcc() const
|
|
|
+{
|
|
|
+ return gain_acc_;
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getGainMag() const
|
|
|
+{
|
|
|
+ return gain_mag_;
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::getSteadyState() const
|
|
|
+{
|
|
|
+ return steady_state_;
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::setBiasAlpha(double bias_alpha)
|
|
|
+{
|
|
|
+ if (bias_alpha >= 0 && bias_alpha <= 1.0)
|
|
|
+ {
|
|
|
+ bias_alpha_ = bias_alpha;
|
|
|
+ return true;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ return false;
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getBiasAlpha() const
|
|
|
+{
|
|
|
+ return bias_alpha_;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::setOrientation(
|
|
|
+ double q0, double q1, double q2, double q3)
|
|
|
+{
|
|
|
+ // Set the state to inverse (state is fixed wrt body).
|
|
|
+ invertQuaternion(q0, q1, q2, q3, q0_, q1_, q2_, q3_);
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+double ComplementaryFilter::getAngularVelocityBiasX() const
|
|
|
+{
|
|
|
+ return wx_bias_;
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getAngularVelocityBiasY() const
|
|
|
+{
|
|
|
+ return wy_bias_;
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getAngularVelocityBiasZ() const
|
|
|
+{
|
|
|
+ return wz_bias_;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::update(double ax, double ay, double az,
|
|
|
+ double wx, double wy, double wz,
|
|
|
+ double dt)
|
|
|
+{
|
|
|
+ if (!initialized_)
|
|
|
+ {
|
|
|
+ // First time - ignore prediction:
|
|
|
+ getMeasurement(ax, ay, az,
|
|
|
+ q0_, q1_, q2_, q3_);
|
|
|
+ initialized_ = true;
|
|
|
+ return;
|
|
|
+ }
|
|
|
+
|
|
|
+ // Bias estimation.
|
|
|
+ if (do_bias_estimation_)
|
|
|
+ updateBiases(ax, ay, az, wx, wy, wz);
|
|
|
+
|
|
|
+ // Prediction.
|
|
|
+ double q0_pred, q1_pred, q2_pred, q3_pred;
|
|
|
+ getPrediction(wx, wy, wz, dt,
|
|
|
+ q0_pred, q1_pred, q2_pred, q3_pred);
|
|
|
+
|
|
|
+ // Correction (from acc):
|
|
|
+ // q_ = q_pred * [(1-gain) * qI + gain * dq_acc]
|
|
|
+ // where qI = identity quaternion
|
|
|
+ double dq0_acc, dq1_acc, dq2_acc, dq3_acc;
|
|
|
+ getAccCorrection(ax, ay, az,
|
|
|
+ q0_pred, q1_pred, q2_pred, q3_pred,
|
|
|
+ dq0_acc, dq1_acc, dq2_acc, dq3_acc);
|
|
|
+
|
|
|
+ double gain;
|
|
|
+ if (do_adaptive_gain_)
|
|
|
+ {
|
|
|
+ gain = getAdaptiveGain(gain_acc_, ax, ay, az);
|
|
|
+
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ gain = gain_acc_;
|
|
|
+
|
|
|
+ }
|
|
|
+
|
|
|
+ scaleQuaternion(gain, dq0_acc, dq1_acc, dq2_acc, dq3_acc);
|
|
|
+
|
|
|
+ quaternionMultiplication(q0_pred, q1_pred, q2_pred, q3_pred,
|
|
|
+ dq0_acc, dq1_acc, dq2_acc, dq3_acc,
|
|
|
+ q0_, q1_, q2_, q3_);
|
|
|
+
|
|
|
+ normalizeQuaternion(q0_, q1_, q2_, q3_);
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::update(double ax, double ay, double az,
|
|
|
+ double wx, double wy, double wz,
|
|
|
+ double mx, double my, double mz,
|
|
|
+ double dt)
|
|
|
+{
|
|
|
+ if (!initialized_)
|
|
|
+ {
|
|
|
+ // First time - ignore prediction:
|
|
|
+ getMeasurement(ax, ay, az,
|
|
|
+ mx, my, mz,
|
|
|
+ q0_, q1_, q2_, q3_);
|
|
|
+ initialized_ = true;
|
|
|
+ return;
|
|
|
+ }
|
|
|
+
|
|
|
+ // Bias estimation.
|
|
|
+ if (do_bias_estimation_)
|
|
|
+ updateBiases(ax, ay, az, wx, wy, wz);
|
|
|
+
|
|
|
+ // Prediction.
|
|
|
+ double q0_pred, q1_pred, q2_pred, q3_pred;
|
|
|
+ getPrediction(wx, wy, wz, dt,
|
|
|
+ q0_pred, q1_pred, q2_pred, q3_pred);
|
|
|
+
|
|
|
+ // Correction (from acc):
|
|
|
+ // q_temp = q_pred * [(1-gain) * qI + gain * dq_acc]
|
|
|
+ // where qI = identity quaternion
|
|
|
+ double dq0_acc, dq1_acc, dq2_acc, dq3_acc;
|
|
|
+ getAccCorrection(ax, ay, az,
|
|
|
+ q0_pred, q1_pred, q2_pred, q3_pred,
|
|
|
+ dq0_acc, dq1_acc, dq2_acc, dq3_acc);
|
|
|
+ double alpha = gain_acc_;
|
|
|
+ if (do_adaptive_gain_)
|
|
|
+ alpha = getAdaptiveGain(gain_acc_, ax, ay, az);
|
|
|
+ scaleQuaternion(alpha, dq0_acc, dq1_acc, dq2_acc, dq3_acc);
|
|
|
+
|
|
|
+ double q0_temp, q1_temp, q2_temp, q3_temp;
|
|
|
+ quaternionMultiplication(q0_pred, q1_pred, q2_pred, q3_pred,
|
|
|
+ dq0_acc, dq1_acc, dq2_acc, dq3_acc,
|
|
|
+ q0_temp, q1_temp, q2_temp, q3_temp);
|
|
|
+
|
|
|
+ // Correction (from mag):
|
|
|
+ // q_ = q_temp * [(1-gain) * qI + gain * dq_mag]
|
|
|
+ // where qI = identity quaternion
|
|
|
+ double dq0_mag, dq1_mag, dq2_mag, dq3_mag;
|
|
|
+ getMagCorrection(mx, my, mz,
|
|
|
+ q0_temp, q1_temp, q2_temp, q3_temp,
|
|
|
+ dq0_mag, dq1_mag, dq2_mag, dq3_mag);
|
|
|
+
|
|
|
+ scaleQuaternion(gain_mag_, dq0_mag, dq1_mag, dq2_mag, dq3_mag);
|
|
|
+
|
|
|
+ quaternionMultiplication(q0_temp, q1_temp, q2_temp, q3_temp,
|
|
|
+ dq0_mag, dq1_mag, dq2_mag, dq3_mag,
|
|
|
+ q0_, q1_, q2_, q3_);
|
|
|
+
|
|
|
+ normalizeQuaternion(q0_, q1_, q2_, q3_);
|
|
|
+}
|
|
|
+
|
|
|
+bool ComplementaryFilter::checkState(double ax, double ay, double az,
|
|
|
+ double wx, double wy, double wz) const
|
|
|
+{
|
|
|
+ double acc_magnitude = sqrt(ax*ax + ay*ay + az*az);
|
|
|
+ if (fabs(acc_magnitude - kGravity) > kAccelerationThreshold)
|
|
|
+ return false;
|
|
|
+
|
|
|
+ if (fabs(wx - wx_prev_) > kDeltaAngularVelocityThreshold ||
|
|
|
+ fabs(wy - wy_prev_) > kDeltaAngularVelocityThreshold ||
|
|
|
+ fabs(wz - wz_prev_) > kDeltaAngularVelocityThreshold)
|
|
|
+ return false;
|
|
|
+
|
|
|
+ if (fabs(wx - wx_bias_) > kAngularVelocityThreshold ||
|
|
|
+ fabs(wy - wy_bias_) > kAngularVelocityThreshold ||
|
|
|
+ fabs(wz - wz_bias_) > kAngularVelocityThreshold)
|
|
|
+ return false;
|
|
|
+
|
|
|
+ return true;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::updateBiases(double ax, double ay, double az,
|
|
|
+ double wx, double wy, double wz)
|
|
|
+{
|
|
|
+ steady_state_ = checkState(ax, ay, az, wx, wy, wz);
|
|
|
+
|
|
|
+ if (steady_state_)
|
|
|
+ {
|
|
|
+ wx_bias_ += bias_alpha_ * (wx - wx_bias_);
|
|
|
+ wy_bias_ += bias_alpha_ * (wy - wy_bias_);
|
|
|
+ wz_bias_ += bias_alpha_ * (wz - wz_bias_);
|
|
|
+ }
|
|
|
+
|
|
|
+ wx_prev_ = wx;
|
|
|
+ wy_prev_ = wy;
|
|
|
+ wz_prev_ = wz;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::getPrediction(
|
|
|
+ double wx, double wy, double wz, double dt,
|
|
|
+ double& q0_pred, double& q1_pred, double& q2_pred, double& q3_pred) const
|
|
|
+{
|
|
|
+ double wx_unb = wx - wx_bias_;
|
|
|
+ double wy_unb = wy - wy_bias_;
|
|
|
+ double wz_unb = wz - wz_bias_;
|
|
|
+
|
|
|
+ q0_pred = q0_ + 0.5*dt*( wx_unb*q1_ + wy_unb*q2_ + wz_unb*q3_);
|
|
|
+ q1_pred = q1_ + 0.5*dt*(-wx_unb*q0_ - wy_unb*q3_ + wz_unb*q2_);
|
|
|
+ q2_pred = q2_ + 0.5*dt*( wx_unb*q3_ - wy_unb*q0_ - wz_unb*q1_);
|
|
|
+ q3_pred = q3_ + 0.5*dt*(-wx_unb*q2_ + wy_unb*q1_ - wz_unb*q0_);
|
|
|
+
|
|
|
+ normalizeQuaternion(q0_pred, q1_pred, q2_pred, q3_pred);
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::getMeasurement(
|
|
|
+ double ax, double ay, double az,
|
|
|
+ double mx, double my, double mz,
|
|
|
+ double& q0_meas, double& q1_meas, double& q2_meas, double& q3_meas)
|
|
|
+{
|
|
|
+ // q_acc is the quaternion obtained from the acceleration vector representing
|
|
|
+ // the orientation of the Global frame wrt the Local frame with arbitrary yaw
|
|
|
+ // (intermediary frame). q3_acc is defined as 0.
|
|
|
+ double q0_acc, q1_acc, q2_acc, q3_acc;
|
|
|
+
|
|
|
+ // Normalize acceleration vector.
|
|
|
+ normalizeVector(ax, ay, az);
|
|
|
+ if (az >=0)
|
|
|
+ {
|
|
|
+ q0_acc = sqrt((az + 1) * 0.5);
|
|
|
+ q1_acc = -ay/(2.0 * q0_acc);
|
|
|
+ q2_acc = ax/(2.0 * q0_acc);
|
|
|
+ q3_acc = 0;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ double X = sqrt((1 - az) * 0.5);
|
|
|
+ q0_acc = -ay/(2.0 * X);
|
|
|
+ q1_acc = X;
|
|
|
+ q2_acc = 0;
|
|
|
+ q3_acc = ax/(2.0 * X);
|
|
|
+ }
|
|
|
+
|
|
|
+ // [lx, ly, lz] is the magnetic field reading, rotated into the intermediary
|
|
|
+ // frame by the inverse of q_acc.
|
|
|
+ // l = R(q_acc)^-1 m
|
|
|
+ double lx = (q0_acc*q0_acc + q1_acc*q1_acc - q2_acc*q2_acc)*mx +
|
|
|
+ 2.0 * (q1_acc*q2_acc)*my - 2.0 * (q0_acc*q2_acc)*mz;
|
|
|
+ double ly = 2.0 * (q1_acc*q2_acc)*mx + (q0_acc*q0_acc - q1_acc*q1_acc +
|
|
|
+ q2_acc*q2_acc)*my + 2.0 * (q0_acc*q1_acc)*mz;
|
|
|
+
|
|
|
+ // q_mag is the quaternion that rotates the Global frame (North West Up) into
|
|
|
+ // the intermediary frame. q1_mag and q2_mag are defined as 0.
|
|
|
+ double gamma = lx*lx + ly*ly;
|
|
|
+ double beta = sqrt(gamma + lx*sqrt(gamma));
|
|
|
+ double q0_mag = beta / (sqrt(2.0 * gamma));
|
|
|
+ double q3_mag = ly / (sqrt(2.0) * beta);
|
|
|
+
|
|
|
+ // The quaternion multiplication between q_acc and q_mag represents the
|
|
|
+ // quaternion, orientation of the Global frame wrt the local frame.
|
|
|
+ // q = q_acc times q_mag
|
|
|
+ quaternionMultiplication(q0_acc, q1_acc, q2_acc, q3_acc,
|
|
|
+ q0_mag, 0, 0, q3_mag,
|
|
|
+ q0_meas, q1_meas, q2_meas, q3_meas );
|
|
|
+ //q0_meas = q0_acc*q0_mag;
|
|
|
+ //q1_meas = q1_acc*q0_mag + q2_acc*q3_mag;
|
|
|
+ //q2_meas = q2_acc*q0_mag - q1_acc*q3_mag;
|
|
|
+ //q3_meas = q0_acc*q3_mag;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+void ComplementaryFilter::getMeasurement(
|
|
|
+ double ax, double ay, double az,
|
|
|
+ double& q0_meas, double& q1_meas, double& q2_meas, double& q3_meas)
|
|
|
+{
|
|
|
+ // q_acc is the quaternion obtained from the acceleration vector representing
|
|
|
+ // the orientation of the Global frame wrt the Local frame with arbitrary yaw
|
|
|
+ // (intermediary frame). q3_acc is defined as 0.
|
|
|
+
|
|
|
+ // Normalize acceleration vector.
|
|
|
+ normalizeVector(ax, ay, az);
|
|
|
+
|
|
|
+ if (az >=0)
|
|
|
+ {
|
|
|
+ q0_meas = sqrt((az + 1) * 0.5);
|
|
|
+ q1_meas = -ay/(2.0 * q0_meas);
|
|
|
+ q2_meas = ax/(2.0 * q0_meas);
|
|
|
+ q3_meas = 0;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ double X = sqrt((1 - az) * 0.5);
|
|
|
+ q0_meas = -ay/(2.0 * X);
|
|
|
+ q1_meas = X;
|
|
|
+ q2_meas = 0;
|
|
|
+ q3_meas = ax/(2.0 * X);
|
|
|
+ }
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::getAccCorrection(
|
|
|
+ double ax, double ay, double az,
|
|
|
+ double p0, double p1, double p2, double p3,
|
|
|
+ double& dq0, double& dq1, double& dq2, double& dq3)
|
|
|
+{
|
|
|
+ // Normalize acceleration vector.
|
|
|
+ normalizeVector(ax, ay, az);
|
|
|
+
|
|
|
+ // Acceleration reading rotated into the world frame by the inverse predicted
|
|
|
+ // quaternion (predicted gravity):
|
|
|
+ double gx, gy, gz;
|
|
|
+ rotateVectorByQuaternion(ax, ay, az,
|
|
|
+ p0, -p1, -p2, -p3,
|
|
|
+ gx, gy, gz);
|
|
|
+
|
|
|
+ // Delta quaternion that rotates the predicted gravity into the real gravity:
|
|
|
+ dq0 = sqrt((gz + 1) * 0.5);
|
|
|
+ dq1 = -gy/(2.0 * dq0);
|
|
|
+ dq2 = gx/(2.0 * dq0);
|
|
|
+ dq3 = 0.0;
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::getMagCorrection(
|
|
|
+ double mx, double my, double mz,
|
|
|
+ double p0, double p1, double p2, double p3,
|
|
|
+ double& dq0, double& dq1, double& dq2, double& dq3)
|
|
|
+{
|
|
|
+ // Magnetic reading rotated into the world frame by the inverse predicted
|
|
|
+ // quaternion:
|
|
|
+ double lx, ly, lz;
|
|
|
+ rotateVectorByQuaternion(mx, my, mz,
|
|
|
+ p0, -p1, -p2, -p3,
|
|
|
+ lx, ly, lz);
|
|
|
+
|
|
|
+ // Delta quaternion that rotates the l so that it lies in the xz-plane
|
|
|
+ // (points north):
|
|
|
+ double gamma = lx*lx + ly*ly;
|
|
|
+ double beta = sqrt(gamma + lx*sqrt(gamma));
|
|
|
+ dq0 = beta / (sqrt(2.0 * gamma));
|
|
|
+ dq1 = 0.0;
|
|
|
+ dq2 = 0.0;
|
|
|
+ dq3 = ly / (sqrt(2.0) * beta);
|
|
|
+}
|
|
|
+
|
|
|
+void ComplementaryFilter::getOrientation(
|
|
|
+ double& q0, double& q1, double& q2, double& q3) const
|
|
|
+{
|
|
|
+ // Return the inverse of the state (state is fixed wrt body).
|
|
|
+ invertQuaternion(q0_, q1_, q2_, q3_, q0, q1, q2, q3);
|
|
|
+}
|
|
|
+
|
|
|
+double ComplementaryFilter::getAdaptiveGain(double alpha, double ax, double ay, double az)
|
|
|
+{
|
|
|
+ double a_mag = sqrt(ax*ax + ay*ay + az*az);
|
|
|
+ double error = fabs(a_mag - kGravity)/kGravity;
|
|
|
+ double factor;
|
|
|
+ double error1 = 0.1;
|
|
|
+ double error2 = 0.2;
|
|
|
+ double m = 1.0/(error1 - error2);
|
|
|
+ double b = 1.0 - m*error1;
|
|
|
+ if (error < error1)
|
|
|
+ factor = 1.0;
|
|
|
+ else if (error < error2)
|
|
|
+ factor = m*error + b;
|
|
|
+ else
|
|
|
+ factor = 0.0;
|
|
|
+ //printf("FACTOR: %f \n", factor);
|
|
|
+ return factor*alpha;
|
|
|
+}
|
|
|
+
|
|
|
+void normalizeVector(double& x, double& y, double& z)
|
|
|
+{
|
|
|
+ double norm = sqrt(x*x + y*y + z*z);
|
|
|
+
|
|
|
+ x /= norm;
|
|
|
+ y /= norm;
|
|
|
+ z /= norm;
|
|
|
+}
|
|
|
+
|
|
|
+void normalizeQuaternion(double& q0, double& q1, double& q2, double& q3)
|
|
|
+{
|
|
|
+ double norm = sqrt(q0*q0 + q1*q1 + q2*q2 + q3*q3);
|
|
|
+ q0 /= norm;
|
|
|
+ q1 /= norm;
|
|
|
+ q2 /= norm;
|
|
|
+ q3 /= norm;
|
|
|
+}
|
|
|
+
|
|
|
+void invertQuaternion(
|
|
|
+ double q0, double q1, double q2, double q3,
|
|
|
+ double& q0_inv, double& q1_inv, double& q2_inv, double& q3_inv)
|
|
|
+{
|
|
|
+ // Assumes quaternion is normalized.
|
|
|
+ q0_inv = q0;
|
|
|
+ q1_inv = -q1;
|
|
|
+ q2_inv = -q2;
|
|
|
+ q3_inv = -q3;
|
|
|
+}
|
|
|
+
|
|
|
+void scaleQuaternion(
|
|
|
+ double gain,
|
|
|
+ double& dq0, double& dq1, double& dq2, double& dq3)
|
|
|
+{
|
|
|
+ if (dq0 < 0.0)//0.9
|
|
|
+ {
|
|
|
+ // Slerp (Spherical linear interpolation):
|
|
|
+ double angle = acos(dq0);
|
|
|
+ double A = sin(angle*(1.0 - gain))/sin(angle);
|
|
|
+ double B = sin(angle * gain)/sin(angle);
|
|
|
+ dq0 = A + B * dq0;
|
|
|
+ dq1 = B * dq1;
|
|
|
+ dq2 = B * dq2;
|
|
|
+ dq3 = B * dq3;
|
|
|
+ }
|
|
|
+ else
|
|
|
+ {
|
|
|
+ // Lerp (Linear interpolation):
|
|
|
+ dq0 = (1.0 - gain) + gain * dq0;
|
|
|
+ dq1 = gain * dq1;
|
|
|
+ dq2 = gain * dq2;
|
|
|
+ dq3 = gain * dq3;
|
|
|
+ }
|
|
|
+
|
|
|
+ normalizeQuaternion(dq0, dq1, dq2, dq3);
|
|
|
+}
|
|
|
+
|
|
|
+void quaternionMultiplication(
|
|
|
+ double p0, double p1, double p2, double p3,
|
|
|
+ double q0, double q1, double q2, double q3,
|
|
|
+ double& r0, double& r1, double& r2, double& r3)
|
|
|
+{
|
|
|
+ // r = p q
|
|
|
+ r0 = p0*q0 - p1*q1 - p2*q2 - p3*q3;
|
|
|
+ r1 = p0*q1 + p1*q0 + p2*q3 - p3*q2;
|
|
|
+ r2 = p0*q2 - p1*q3 + p2*q0 + p3*q1;
|
|
|
+ r3 = p0*q3 + p1*q2 - p2*q1 + p3*q0;
|
|
|
+}
|
|
|
+
|
|
|
+void rotateVectorByQuaternion(
|
|
|
+ double x, double y, double z,
|
|
|
+ double q0, double q1, double q2, double q3,
|
|
|
+ double& vx, double& vy, double& vz)
|
|
|
+{
|
|
|
+ vx = (q0*q0 + q1*q1 - q2*q2 - q3*q3)*x + 2*(q1*q2 - q0*q3)*y + 2*(q1*q3 + q0*q2)*z;
|
|
|
+ vy = 2*(q1*q2 + q0*q3)*x + (q0*q0 - q1*q1 + q2*q2 - q3*q3)*y + 2*(q2*q3 - q0*q1)*z;
|
|
|
+ vz = 2*(q1*q3 - q0*q2)*x + 2*(q2*q3 + q0*q1)*y + (q0*q0 - q1*q1 - q2*q2 + q3*q3)*z;
|
|
|
+}
|
|
|
+
|
|
|
+
|
|
|
+} // namespace imu_tools
|