ardupilot/libraries/AP_OpticalFlow/AP_OpticalFlow_Calibrator.cpp

326 lines
13 KiB
C++

/*
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "AP_OpticalFlow_Calibrator.h"
#include <GCS_MAVLink/GCS.h>
#include <AP_Logger/AP_Logger.h>
const uint32_t AP_OPTICALFLOW_CAL_TIMEOUT_SEC = 120; // calibration timesout after 120 seconds
const uint32_t AP_OPTICALFLOW_CAL_STATUSINTERVAL_SEC = 3; // status updates printed at 3 second intervals
const float AP_OPTICALFLOW_CAL_YAW_MAX_RADS = radians(20); // maximum yaw rotation (must be low to ensure good scaling)
const float AP_OPTICALFLOW_CAL_ROLLPITCH_MIN_RADS = radians(20); // minimum acceptable roll or pitch rotation
const float AP_OPTICALFLOW_CAL_SCALE_MIN = 0.20f; // min acceptable scaling value. If resulting scaling is below this then the calibration fails
const float AP_OPTICALFLOW_CAL_SCALE_MAX = 4.0f; // max acceptable scaling value. If resulting scaling is above this then the calibration fails
const float AP_OPTICALFLOW_CAL_FITNESS_THRESH = 0.5f; // min acceptable fitness
const float AP_OPTICALFLOW_CAL_RMS_FAILED = 1.0e30f; // calc_mean_squared_residuals returns this value when it fails to calculate a good value
extern const AP_HAL::HAL& hal;
// start the calibration
void AP_OpticalFlow_Calibrator::start()
{
// exit immediately if already running
if (_cal_state == CalState::RUNNING) {
return;
}
_cal_state = CalState::RUNNING;
_start_time_ms = AP_HAL::millis();
// clear samples buffers
_cal_data[0].num_samples = 0;
_cal_data[1].num_samples = 0;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "FlowCal: Started");
}
void AP_OpticalFlow_Calibrator::stop()
{
// exit immediately if already stopped
if (_cal_state == CalState::NOT_STARTED) {
return;
}
_cal_state = CalState::NOT_STARTED;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "FlowCal: Stopped");
}
// update the state machine and calculate scaling
bool AP_OpticalFlow_Calibrator::update()
{
// prefix for reporting
const char* prefix_str = "FlowCal:";
// while running add samples
if (_cal_state == CalState::RUNNING) {
uint32_t now_ms = AP_HAL::millis();
uint32_t timestamp_ms;
Vector2f flow_rate, body_rate, los_pred;
if (AP::ahrs().getOptFlowSample(timestamp_ms, flow_rate, body_rate, los_pred)) {
add_sample(timestamp_ms, flow_rate, body_rate, los_pred);
// while collecting samples display percentage complete
if (now_ms - _last_report_ms > AP_OPTICALFLOW_CAL_STATUSINTERVAL_SEC * 1000UL) {
_last_report_ms = now_ms;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s x:%d%% y:%d%%",
prefix_str,
(int)((_cal_data[0].num_samples * 100.0 / AP_OPTICALFLOW_CAL_MAX_SAMPLES)),
(int)((_cal_data[1].num_samples * 100.0 / AP_OPTICALFLOW_CAL_MAX_SAMPLES)));
}
// advance state once sample buffers are full
if (sample_buffers_full()) {
_cal_state = CalState::READY_TO_CALIBRATE;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s samples collected", prefix_str);
}
}
// check for timeout
if (now_ms - _start_time_ms > AP_OPTICALFLOW_CAL_TIMEOUT_SEC * 1000UL) {
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s timeout", prefix_str);
_cal_state = CalState::FAILED;
}
}
// start calibration
if (_cal_state == CalState::READY_TO_CALIBRATE) {
// run calibration and mark failure or success
if (run_calibration()) {
_cal_state = CalState::SUCCESS;
return true;
} else {
_cal_state = CalState::FAILED;
}
}
// return indicating calibration is not complete
return false;
}
// get final scaling values
// scaling values used during sample collection should be multiplied by these scalars
Vector2f AP_OpticalFlow_Calibrator::get_scalars()
{
// return best scaling values
return Vector2f{_cal_data[0].best_scalar, _cal_data[1].best_scalar};
}
// add new sample to the calibrator
void AP_OpticalFlow_Calibrator::add_sample(uint32_t timestamp_ms, const Vector2f& flow_rate, const Vector2f& body_rate, const Vector2f& los_pred)
{
// return immediately if not running
if (_cal_state != CalState::RUNNING) {
return;
}
// check for duplicates
if (timestamp_ms == _last_sample_timestamp_ms) {
return;
}
_last_sample_timestamp_ms = timestamp_ms;
// check yaw rotation is low
const Vector3f gyro = AP::ahrs().get_gyro();
if (fabsf(gyro.z) > AP_OPTICALFLOW_CAL_YAW_MAX_RADS) {
return;
}
// check enough roll or pitch movement and record sample
const bool rates_x_sufficient = (fabsf(body_rate.x) >= AP_OPTICALFLOW_CAL_ROLLPITCH_MIN_RADS) && (fabsf(flow_rate.x) >= AP_OPTICALFLOW_CAL_ROLLPITCH_MIN_RADS);
if (rates_x_sufficient && (_cal_data[0].num_samples < ARRAY_SIZE(_cal_data[0].samples))) {
log_sample(0, _cal_data[0].num_samples, flow_rate.x, body_rate.x, los_pred.x);
_cal_data[0].samples[_cal_data[0].num_samples].flow_rate = flow_rate.x;
_cal_data[0].samples[_cal_data[0].num_samples].body_rate = body_rate.x;
_cal_data[0].samples[_cal_data[0].num_samples].los_pred = los_pred.x;
_cal_data[0].num_samples++;
}
const bool rates_y_sufficient = (fabsf(body_rate.y) >= AP_OPTICALFLOW_CAL_ROLLPITCH_MIN_RADS) && (fabsf(flow_rate.y) >= AP_OPTICALFLOW_CAL_ROLLPITCH_MIN_RADS);
if (rates_y_sufficient && (_cal_data[1].num_samples < ARRAY_SIZE(_cal_data[1].samples))) {
log_sample(1, _cal_data[1].num_samples, flow_rate.y, body_rate.y, los_pred.y);
_cal_data[1].samples[_cal_data[1].num_samples].flow_rate = flow_rate.y;
_cal_data[1].samples[_cal_data[1].num_samples].body_rate = body_rate.y;
_cal_data[1].samples[_cal_data[1].num_samples].los_pred = los_pred.y;
_cal_data[1].num_samples++;
}
}
// returns true once the sample buffer is full
bool AP_OpticalFlow_Calibrator::sample_buffers_full() const
{
return ((_cal_data[0].num_samples >= ARRAY_SIZE(_cal_data[0].samples)) && (_cal_data[1].num_samples >= ARRAY_SIZE(_cal_data[1].samples)));
}
// run calibration algorithm for both axis
// returns true on success and updates _cal_data[0,1].best_scale and best_scale_fitness
bool AP_OpticalFlow_Calibrator::run_calibration()
{
// run calibration for x and y axis
const bool calx_res = calc_scalars(0, _cal_data[0].best_scalar, _cal_data[0].best_scalar_fitness);
const bool caly_res = calc_scalars(1, _cal_data[1].best_scalar, _cal_data[1].best_scalar_fitness);
return calx_res && caly_res;
}
// Run Gauss Newton fitting algorithm for all samples of the given axis
// returns a scalar and fitness (lower numbers mean a better result) in the arguments provided
bool AP_OpticalFlow_Calibrator::calc_scalars(uint8_t axis, float& scalar, float& fitness)
{
// prefix for reporting
const char* prefix_str = "FlowCal:";
const char* axis_str = axis == 0 ? "x" : "y";
// check we have samples
// this should never fail because this method should only be called once the sample buffer is full
const uint8_t num_samples = _cal_data[axis].num_samples;
if (num_samples == 0) {
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s failed because no samples", prefix_str);
return false;
}
// calculate total absolute residual from all samples
float total_abs_residual = 0;
for (uint8_t i = 0; i < num_samples; i++) {
const sample_t& samplei = _cal_data[axis].samples[i];
total_abs_residual += fabsf(calc_sample_residual(samplei, 1.0));
}
// if there are no residuals then scaling is perfect
if (is_zero(total_abs_residual)) {
scalar = 1.0;
fitness = 0;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s perfect scalar%s of 1.0", prefix_str, axis_str);
return true;
}
// for each sample calculate the residual and scalar that best reduces the residual
float best_scalar_total = 0;
for (uint8_t i = 0; i < num_samples; i++) {
float sample_best_scalar;
const sample_t& samplei = _cal_data[axis].samples[i];
if (!calc_sample_best_scalar(samplei, sample_best_scalar)) {
// failed to find the best scalar for a single sample
// this should never happen because of checks when capturing samples
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s failed because of zero flow rate", prefix_str);
INTERNAL_ERROR(AP_InternalError::error_t::flow_of_control);
return false;
}
const float sample_residual = calc_sample_residual(samplei, 1.0);
best_scalar_total += sample_best_scalar * fabsf(sample_residual) / total_abs_residual;
}
// check for out of range results
if (best_scalar_total < AP_OPTICALFLOW_CAL_SCALE_MIN) {
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s scalar%s:%4.3f too low (<%3.1f)", prefix_str, axis_str, (double)best_scalar_total, (double)AP_OPTICALFLOW_CAL_SCALE_MIN);
return false;
}
if (best_scalar_total > AP_OPTICALFLOW_CAL_SCALE_MAX) {
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s scalar%s:%4.3f too high (>%3.1f)", prefix_str, axis_str, (double)best_scalar_total, (double)AP_OPTICALFLOW_CAL_SCALE_MAX);
return false;
}
// check for poor fitness
float fitness_new = calc_mean_squared_residuals(axis, best_scalar_total);
if (fitness_new > AP_OPTICALFLOW_CAL_FITNESS_THRESH) {
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s scalar%s:%4.3f fit:%4.3f too high (>%3.1f)", prefix_str, axis_str, (double)scalar, (double)fitness_new, (double)AP_OPTICALFLOW_CAL_FITNESS_THRESH);
return false;
}
// success if fitness has improved
float fitness_orig = calc_mean_squared_residuals(axis, 1.0);
if (fitness_new <= fitness_orig) {
scalar = best_scalar_total;
fitness = fitness_new;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s scalar%s:%4.3f fit:%4.2f", prefix_str, axis_str, (double)scalar, (double)fitness);
return true;
}
// failed to find a better scalar than 1.0
scalar = 1.0;
fitness = fitness_orig;
GCS_SEND_TEXT(MAV_SEVERITY_INFO, "%s no better scalar%s:%4.3f (fit:%4.3f > orig:%4.3f)", prefix_str, axis_str, (double)best_scalar_total, (double)fitness_new, (double)fitness_orig);
return false;
}
// calculate a single sample's residual
float AP_OpticalFlow_Calibrator::calc_sample_residual(const sample_t& sample, float scalar) const
{
return (sample.body_rate + ((sample.flow_rate * scalar) - sample.los_pred));
}
// calculate the scalar that minimises the residual for a single sample
// returns true on success and populates the best_scalar argument
bool AP_OpticalFlow_Calibrator::calc_sample_best_scalar(const sample_t& sample, float& best_scalar) const
{
// if sample's flow_rate is zero scalar has no effect
// this should never happen because samples should have been checked before being added
if (is_zero(sample.flow_rate)) {
return false;
}
best_scalar = (sample.los_pred - sample.body_rate) / sample.flow_rate;
return true;
}
// calculate mean squared residual for all samples of a single axis (0 or 1) given a scalar parameter
float AP_OpticalFlow_Calibrator::calc_mean_squared_residuals(uint8_t axis, float scalar) const
{
// sanity check axis
if (axis >= 2) {
return AP_OPTICALFLOW_CAL_RMS_FAILED;
}
// calculate and sum residuals of each sample
float sum = 0.0f;
uint16_t num_samples = 0;
for (uint8_t i = 0; i < _cal_data[axis].num_samples; i++) {
sum += sq(calc_sample_residual(_cal_data[axis].samples[i], scalar));
num_samples++;
}
// return a huge residual if no samples
if (num_samples == 0) {
return AP_OPTICALFLOW_CAL_RMS_FAILED;
}
sum /= num_samples;
return sum;
}
// log all samples
void AP_OpticalFlow_Calibrator::log_sample(uint8_t axis, uint8_t sample_num, float flow_rate, float body_rate, float los_pred)
{
// @LoggerMessage: OFCA
// @Description: Optical Flow Calibration sample
// @Field: TimeUS: Time since system startup
// @Field: Axis: Axis (X=0 Y=1)
// @Field: Num: Sample number
// @Field: FRate: Flow rate
// @Field: BRate: Body rate
// @Field: LPred: Los pred
AP::logger().Write(
"OFCA",
"TimeUS,Axis,Num,FRate,BRate,LPred",
"s#-EEE",
"F00000",
"QBBfff",
AP_HAL::micros64(),
(unsigned)axis,
(unsigned)sample_num,
(double)flow_rate,
(double)body_rate,
(double)los_pred);
}