lib: remove bezier

Signed-off-by: Silvan <silvan@auterion.com>
This commit is contained in:
Silvan
2025-01-10 11:35:35 +01:00
committed by Silvan Fuhrer
parent b34a5eb6f7
commit b7b6d45e18
9 changed files with 1 additions and 1039 deletions

View File

@@ -5,7 +5,6 @@
# tests command arguments
set(tests
atomic_bitset
bezier
bitset
bson
dataman

View File

@@ -35,7 +35,6 @@ add_subdirectory(adsb EXCLUDE_FROM_ALL)
add_subdirectory(airspeed EXCLUDE_FROM_ALL)
add_subdirectory(atmosphere EXCLUDE_FROM_ALL)
add_subdirectory(battery EXCLUDE_FROM_ALL)
add_subdirectory(bezier EXCLUDE_FROM_ALL)
add_subdirectory(button EXCLUDE_FROM_ALL)
add_subdirectory(cdev EXCLUDE_FROM_ALL)
add_subdirectory(cdrstream EXCLUDE_FROM_ALL)

View File

@@ -1,163 +0,0 @@
/****************************************************************************
*
* Copyright (c) 2020 PX4 Development Team. 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 PX4 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 OWNER 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.
*
****************************************************************************/
/**
* @file BezierN.cpp
* Bezier function
*
* @author Julian Kent <julian@auterion.com>
*/
#include <bezier/BezierN.hpp>
#include <matrix/Dual.hpp>
namespace
{
/*
* Generic in-place bezier implementation. Leaves result in first element.
*
*/
template <typename Scalar, size_t D>
void calculateBezier(matrix::Vector<Scalar, D> *positions, int N, Scalar t, Scalar one_minus_t)
{
for (int bezier_order = 1; bezier_order < N; bezier_order++) {
for (int i = 0; i < N - bezier_order; i++) {
positions[i] = positions[i] * one_minus_t + positions[i + 1] * t;
}
}
}
}
namespace bezier
{
bool calculateBezierPosVel(const matrix::Vector3f *positions, int N, float t,
matrix::Vector3f &position, matrix::Vector3f &velocity)
{
if (positions == nullptr || N == 0 || t < 0 || t > 1) {
return false;
}
using Df = matrix::Dual<float, 1>;
using Vector3Df = matrix::Vector3<Df>;
Vector3Df intermediates[N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < 3; j++) {
intermediates[i](j) = positions[i](j);
}
}
Df dual_t(t, 0); // derivative with respect to time
calculateBezier(intermediates, N, dual_t, Df(1) - dual_t);
position = matrix::collectReals(intermediates[0]);
velocity = matrix::collectDerivatives(intermediates[0]);
return true;
}
bool calculateBezierPosVelAcc(const matrix::Vector3f *positions, int N, float t,
matrix::Vector3f &position, matrix::Vector3f &velocity, matrix::Vector3f &acceleration)
{
if (positions == nullptr || N == 0 || t < 0 || t > 1) {
return false;
}
using Df = matrix::Dual<float, 1>;
using DDf = matrix::Dual<Df, 1>;
using Vector3DDf = matrix::Vector3<DDf>;
Vector3DDf intermediates[N];
for (int i = 0; i < N; i++) {
for (int j = 0; j < 3; j++) {
intermediates[i](j) = Df(positions[i](j));
}
}
DDf dual_t(Df(t, 0), 0); // 1st and 2nd derivative with respect to time
calculateBezier(intermediates, N, dual_t, Df(1) - dual_t);
position = matrix::collectReals(matrix::collectReals(intermediates[0]));
velocity = matrix::collectReals(matrix::collectDerivatives(intermediates[0]));
acceleration = matrix::collectDerivatives(matrix::collectDerivatives(intermediates[0]));
return true;
}
bool calculateBezierYaw(const float *setpoints, int N, float t, float &yaw_setpoint, float &yaw_vel_setpoint)
{
if (setpoints == nullptr || N == 0 || t < 0 || t > 1) {
return false;
}
using Df = matrix::Dual<float, 1>;
using Vector1Df = matrix::Vector<Df, 1>;
Vector1Df intermediates[N];
// all yaw setpoints are wrapped relative to the starting yaw
const float offset = setpoints[0];
for (int i = 0; i < N; i++) {
intermediates[i](0) = matrix::wrap_pi(setpoints[i] - offset);
}
Df dual_t (t, 0); // derivative with respect to time
calculateBezier(intermediates, N, dual_t, Df(1) - dual_t);
Df result = intermediates[0](0);
yaw_setpoint = matrix::wrap_pi(result.value + offset);
yaw_vel_setpoint = result.derivative(0);
return true;
}
bool calculateT(int64_t start_time, int64_t end_time, int64_t now, float &T)
{
if (now < start_time || end_time < now) {
return false;
}
int64_t total_duration = end_time - start_time;
int64_t elapsed_duration = now - start_time;
T = (float) elapsed_duration / (float) total_duration;
return true;
}
}

View File

@@ -1,76 +0,0 @@
/****************************************************************************
*
* Copyright (C) 2020 PX4 Development Team. 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 PX4 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 OWNER 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.
*
****************************************************************************/
/**
* @file BezierN.hpp
*
* @author Julian Kent <julian@auterion.com>
*
* N-order Bezier library designed for time-aware trajectory tracking
*/
#pragma once
#include <matrix/math.hpp>
namespace bezier
{
/*
* Calculates the location and velocity with respect to T on a given bezier curve of any order.
*
*/
bool calculateBezierPosVel(const matrix::Vector3f *positions, int N, float t,
matrix::Vector3f &position, matrix::Vector3f &velocity);
/*
* Calculates the position, velocity and acceleration with respect to T on a given bezier curve of any order.
*
*/
bool calculateBezierPosVelAcc(const matrix::Vector3f *positions, int N, float t,
matrix::Vector3f &position, matrix::Vector3f &velocity, matrix::Vector3f &acceleration);
/*
* Calculates the position and velocity of yaw with respect to t on a bezier curve.
* All yaw setpoints are wrapped relative to the starting yaw.
*
*/
bool calculateBezierYaw(const float *setpoints, int N, float t, float &yaw_setpoint, float &yaw_vel_setpoint);
/*
* Calculates the fraction between the begin and end time which can be used for fast bezier curve lookups
*/
bool calculateT(int64_t start_time, int64_t end_time, int64_t now, float &T);
}

View File

@@ -1,325 +0,0 @@
/****************************************************************************
*
* Copyright (C) 2019 PX4 Development Team. 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 PX4 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 OWNER 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.
*
****************************************************************************/
/**
* Test code for the Velocity Smoothing library
* Run this test only using make tests TESTFILTER=BezierN
*
* @author Julian Kent <julian@auterion.com>
*/
#include <gtest/gtest.h>
#include <matrix/matrix/math.hpp>
#include "BezierN.hpp"
TEST(BezierN_calculateBezier, checks_validity)
{
matrix::Vector3f points[10];
matrix::Vector3f a, b;
EXPECT_FALSE(bezier::calculateBezierPosVel(nullptr, 10, 0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierPosVel(points, 0, 0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierPosVel(points, 10, -0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierPosVel(points, 10, 1.5f, a, b));
}
TEST(BezierN_calculateBezier, checks_validity_accel)
{
matrix::Vector3f points[10];
matrix::Vector3f a, b, c;
EXPECT_FALSE(bezier::calculateBezierPosVelAcc(nullptr, 10, 0.5f, a, b, c));
EXPECT_FALSE(bezier::calculateBezierPosVelAcc(points, 0, 0.5f, a, b, c));
EXPECT_FALSE(bezier::calculateBezierPosVelAcc(points, 10, -0.5f, a, b, c));
EXPECT_FALSE(bezier::calculateBezierPosVelAcc(points, 10, 1.5f, a, b, c));
}
TEST(BezierN_calculateBezier, work_1_point)
{
// GIVEN: a single point bezier curve
matrix::Vector3f points[2] = {matrix::Vector3f(1, 2, 3), matrix::Vector3f(NAN, NAN, NAN)};
matrix::Vector3f pos, vel;
pos *= NAN;
vel *= NAN;
// WHEN: we get the half-way point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 1, 0.5f, pos, vel));
// THEN: it should be the same as the point, and the velocity should be 0
EXPECT_EQ((pos - points[0]).norm(), 0.f);
EXPECT_EQ(vel.norm(), 0.f);
}
TEST(BezierN_calculateBezier, works_2_points)
{
// GIVEN: a 2-point bezier curve
matrix::Vector3f points[3] = {matrix::Vector3f(1, 2, 3), matrix::Vector3f(5, 0, 1), matrix::Vector3f(NAN, NAN, NAN)};
matrix::Vector3f pos, vel;
pos *= NAN;
vel *= NAN;
// WHEN: we get the half-way point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 2, 0.5f, pos, vel));
// THEN: the position should be the mid-point between the start and end, and velocity should be the length
EXPECT_EQ((pos - matrix::Vector3f(3, 1, 2)).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[1] - points[0])).norm(), 0.f);
// WHEN: we get the beginning point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 2, 0.f, pos, vel));
// THEN: the position should be the first point, and the velocity should still be the length
EXPECT_EQ((pos - points[0]).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[1] - points[0])).norm(), 0.f);
// WHEN: we get the end point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 2, 1.f, pos, vel));
// THEN: the position should be the first point, and the velocity should still be the length
EXPECT_EQ((pos - points[1]).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[1] - points[0])).norm(), 0.f);
}
TEST(BezierN_calculateBezier, works_3_points_zero_accel)
{
// GIVEN: 3 points bezier, evenly spaced in a straight line
matrix::Vector3f points[4] = {matrix::Vector3f(1, 2, 3), matrix::Vector3f(5, 0, 1), matrix::Vector3f(9, -2, -1), matrix::Vector3f(NAN, NAN, NAN)};
matrix::Vector3f pos, vel;
pos *= NAN;
vel *= NAN;
// WHEN: we get the half-way point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 0.5f, pos, vel));
// THEN: it should be the middle point, with velocity of 1st to last
EXPECT_FLOAT_EQ((pos - points[1]).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[2] - points[0])).norm(), 0.f);
matrix::Vector3f pos2, vel2, accel2;
// WHEN: we use the accel interface
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 0.5f, pos2, vel2, accel2));
// THEN: it should give same position, velocity as the non-accel interface, and zero accel (since this curve is 0 accel)
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel2 - vel).norm(), 0.f);
EXPECT_FLOAT_EQ(accel2.norm(), 0.f);
// WHEN: we check at the beginning
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 0.f, pos, vel));
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 0.f, pos2, vel2, accel2));
// THEN: it should be the starting point and same velocity
EXPECT_FLOAT_EQ((pos - points[0]).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[2] - points[0])).norm(), 0.f);
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel2 - vel).norm(), 0.f);
EXPECT_FLOAT_EQ(accel2.norm(), 0.f);
// WHEN: we check at the end
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 1.f, pos, vel));
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 1.f, pos2, vel2, accel2));
// THEN: it should be the ending point and same velocity
EXPECT_FLOAT_EQ((pos - points[2]).norm(), 0.f);
EXPECT_FLOAT_EQ((vel - (points[2] - points[0])).norm(), 0.f);
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel2 - vel).norm(), 0.f);
EXPECT_FLOAT_EQ(accel2.norm(), 0.f);
}
TEST(BezierN_calculateBezier, works_3_points_accel)
{
// GIVEN: 3 points bezier, in a curve
matrix::Vector3f points[4] = {matrix::Vector3f(1, 2, 3), matrix::Vector3f(5, 0, 1), matrix::Vector3f(19, -8, 1), matrix::Vector3f(NAN, NAN, NAN)};
matrix::Vector3f pos, vel;
pos *= NAN;
vel *= NAN;
matrix::Vector3f pos2;
pos2 *= NAN;
matrix::Vector3f accel_start, accel_mid, accel_end;
matrix::Vector3f vel_start, vel_mid, vel_end;
// WHEN: we check at the beginning
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 0.f, pos, vel));
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 0.f, pos2, vel_start, accel_start));
// THEN: it should give same position, velocity as the non-accel interface, and non-zero accel
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel_start - vel).norm(), 0.f);
EXPECT_GT(accel_start.norm(), 0.f);
// WHEN: we use the accel interface to get the half-way point
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 0.5f, pos, vel));
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 0.5f, pos2, vel_mid, accel_mid));
// THEN: the values should matche between accel and non-accel version
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel_mid - vel).norm(), 0.f);
// AND: the accel should be the same as the start
EXPECT_FLOAT_EQ((accel_mid - accel_start).norm(), 0.f);
// WHEN: we check at the end
EXPECT_TRUE(bezier::calculateBezierPosVel(points, 3, 1.f, pos, vel));
EXPECT_TRUE(bezier::calculateBezierPosVelAcc(points, 3, 1.f, pos2, vel_end, accel_end));
// THEN: it should be the ending point, and accel should match
EXPECT_FLOAT_EQ((pos - points[2]).norm(), 0.f);
EXPECT_FLOAT_EQ((pos2 - pos).norm(), 0.f);
EXPECT_FLOAT_EQ((vel_end - vel).norm(), 0.f);
EXPECT_FLOAT_EQ((accel_end - accel_start).norm(), 0.f);
// FINALLY: mid point velocity should be average of start and end velocity
EXPECT_FLOAT_EQ((vel_mid - 0.5f * (vel_start + vel_end)).norm(), 0.f);
}
TEST(BezierN_calculateBezierYaw, checks_validity)
{
float points[10];
float a, b;
EXPECT_FALSE(bezier::calculateBezierYaw(nullptr, 10, 0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierYaw(points, 0, 0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierYaw(points, 10, -0.5f, a, b));
EXPECT_FALSE(bezier::calculateBezierYaw(points, 10, 1.5f, a, b));
}
TEST(BezierN_calculateBezierYaw, work_1_point)
{
// GIVEN: a single yaw point
float points[2] = {M_PI / 2, NAN};
float yaw, yaw_speed;
// WHEN: we use it as a 1-point bezier curve
EXPECT_TRUE(bezier::calculateBezierYaw(points, 1, 0.5f, yaw, yaw_speed));
// THEN: it should have that same value, and the velocity should be 0
EXPECT_FLOAT_EQ(yaw, M_PI / 2);
EXPECT_FLOAT_EQ(yaw_speed, 0);
}
TEST(BezierN_calculateBezierYaw, work_2_points)
{
// GIVEN: a single yaw point
float points[3] = {0, M_PI / 2, NAN};
float yaw, yaw_speed;
// WHEN: we get the beginning
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 0.f, yaw, yaw_speed));
// THEN: it should have the beginning value, and the velocity should be the difference between first and last
EXPECT_FLOAT_EQ(yaw, 0);
EXPECT_FLOAT_EQ(yaw_speed, M_PI / 2);
// WHEN: we get the middle
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 0.5f, yaw, yaw_speed));
// THEN: it should have the beginning value, and the velocity should be the difference between first and last
EXPECT_FLOAT_EQ(yaw, M_PI / 4);
EXPECT_FLOAT_EQ(yaw_speed, M_PI / 2);
// WHEN: we get the end
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 1.f, yaw, yaw_speed));
// THEN: it should have the beginning value, and the velocity should be the difference between first and last
EXPECT_FLOAT_EQ(yaw, M_PI / 2);
EXPECT_FLOAT_EQ(yaw_speed, M_PI / 2);
}
TEST(BezierN_calculateBezierYaw, work_2_points_wrap)
{
// GIVEN: 2 yaw points on either side of the +- PI wrap line
float points[3] = {-M_PI + 0.1, M_PI - 0.1, NAN};
float yaw, yaw_speed;
// WHEN: we get the beginning
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 0.f, yaw, yaw_speed));
// THEN: it should have the beginning value, and the velocity should be the wrapped distance between first and last
EXPECT_FLOAT_EQ(yaw, -M_PI + 0.1);
EXPECT_NEAR(yaw_speed, -0.2, 1e-6f);
// WHEN: we get the middle
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 0.5f, yaw, yaw_speed));
// THEN: it should have the wrapped middle value, and the velocity should be the wrapped distance between first and last
EXPECT_FLOAT_EQ(matrix::wrap_pi(yaw - float(M_PI)), 0);
EXPECT_NEAR(yaw_speed, -0.2, 1e-6f);
// WHEN: we get the end
EXPECT_TRUE(bezier::calculateBezierYaw(points, 2, 1.f, yaw, yaw_speed));
// THEN: it should have the end value, and the velocity should be the wrapped distance between first and last
EXPECT_FLOAT_EQ(yaw, M_PI - 0.1);
EXPECT_NEAR(yaw_speed, -0.2, 1e-6f);
}
TEST(BezierN_calculateT, rejects_bad_timestamps)
{
float f = NAN;
EXPECT_FALSE(bezier::calculateT(100, 1000, 99, f));
EXPECT_FALSE(bezier::calculateT(100, 1000, 1001, f));
EXPECT_FALSE(bezier::calculateT(1001, 1000, 1001, f));
}
TEST(BezierN_calculateT, begin_middle_end)
{
float f = NAN;
EXPECT_TRUE(bezier::calculateT(100, 1000, 100, f));
EXPECT_FLOAT_EQ(f, 0.f);
EXPECT_TRUE(bezier::calculateT(100, 1000, 550, f));
EXPECT_FLOAT_EQ(f, 0.5f);
EXPECT_TRUE(bezier::calculateT(100, 1000, 1000, f));
EXPECT_FLOAT_EQ(f, 1.f);
}
TEST(BezierN_calculateT, giant_offset)
{
int64_t offset = 0xFFFFFFFFFFFF; // 48 bit max
float f = NAN;
EXPECT_TRUE(bezier::calculateT(offset + 100, offset + 1000, offset + 100, f));
EXPECT_FLOAT_EQ(f, 0.f);
EXPECT_TRUE(bezier::calculateT(offset + 100, offset + 1000, offset + 550, f));
EXPECT_FLOAT_EQ(f, 0.5f);
EXPECT_TRUE(bezier::calculateT(offset + 100, offset + 1000, offset + 1000, f));
EXPECT_FLOAT_EQ(f, 1.f);
}

View File

@@ -1,223 +0,0 @@
/****************************************************************************
*
* Copyright (c) 2018 PX4 Development Team. 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 PX4 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 OWNER 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.
*
****************************************************************************/
/**
* @file BezierQuad.cpp
* Bezier function
*
* @author Dennis Mannhart <dennis.mannhart@gmail.com>
*/
// The header includes this implementation, so in order to have this library shared with other
// .cpp files and avoid duplication of constants, we need to prevent including it twice
#ifndef BEZIER_QUAD_CPP
#define BEZIER_QUAD_CPP
#include "BezierQuad.hpp"
namespace bezier
{
static constexpr double GOLDEN_RATIO = 1.6180339887; //(sqrt(5)+1)/2
static constexpr double RESOLUTION = 0.0001; //End criterion for golden section search
template<typename Tp>
void BezierQuad<Tp>::setBezier(const Vector3_t &pt0, const Vector3_t &ctrl, const Vector3_t &pt1,
Tp duration)
{
_pt0 = pt0;
_ctrl = ctrl;
_pt1 = pt1;
_duration = duration;
_cached_resolution = (Tp)(-1);
}
template<typename Tp>
void BezierQuad<Tp>::getBezier(Vector3_t &pt0, Vector3_t &ctrl, Vector3_t &pt1)
{
pt0 = _pt0;
ctrl = _ctrl;
pt1 = _pt1;
}
template<typename Tp>
matrix::Vector<Tp, 3> BezierQuad<Tp>::getPoint(const Tp t)
{
return (_pt0 * ((Tp)1 - t / _duration) * ((Tp)1 - t / _duration) + _ctrl * (Tp)2 * ((
Tp)1 - t / _duration) * t / _duration + _pt1 *
t / _duration * t / _duration);
}
template<typename Tp>
matrix::Vector<Tp, 3> BezierQuad<Tp>::getVelocity(const Tp t)
{
return (((_ctrl - _pt0) * _duration + (_pt0 - _ctrl * (Tp)2 + _pt1) * t) * (Tp)2 / (_duration * _duration));
}
template<typename Tp>
matrix::Vector<Tp, 3> BezierQuad<Tp>::getAcceleration()
{
return ((_pt0 - _ctrl * (Tp)2 + _pt1) * (Tp)2 / (_duration * _duration));
}
template<typename Tp>
void BezierQuad<Tp>::getStates(Vector3_t &point, Vector3_t &vel, Vector3_t &acc, const Tp time)
{
point = getPoint(time);
vel = getVelocity(time);
acc = getAcceleration();
}
template<typename Tp>
void BezierQuad<Tp>::getStatesClosest(Vector3_t &point, Vector3_t &vel, Vector3_t &acc,
const Vector3_t pose)
{
// get t that corresponds to point closest on bezier point
Tp t = _goldenSectionSearch(pose);
// get states corresponding to t
getStates(point, vel, acc, t);
}
template<typename Tp>
void BezierQuad<Tp>::setBezFromVel(const Vector3_t &ctrl, const Vector3_t &vel0, const Vector3_t &vel1,
const Tp duration)
{
// update bezier points
_ctrl = ctrl;
_duration = duration;
_pt0 = _ctrl - vel0 * _duration / (Tp)2;
_pt1 = _ctrl + vel1 * _duration / (Tp)2;
_cached_resolution = (Tp)(-1);
}
template<typename Tp>
Tp BezierQuad<Tp>::getArcLength(const Tp resolution)
{
// we don't need to recompute arc length if:
// 1. _cached_resolution is up to date; 2. _cached_resolution is smaller than desired resolution (= more accurate)
if ((_cached_resolution > (Tp)0) && (_cached_resolution <= resolution)) {
return _cached_arc_length;
}
// get number of elements
int n = (int)(roundf(_duration / resolution));
// check if n is even
if (n % 2 == 1) {
n += 1;
}
// step size
Tp h = (_duration) / n;
// get integration
_cached_arc_length = (Tp)0;
Vector3_t y;
for (int i = 1; i < n; i++) {
y = getVelocity(h * i);
if (i % 2 == 1) {
_cached_arc_length += (Tp)4 * y.length();
} else {
_cached_arc_length += (Tp)2 * y.length();
}
}
// velocity length at start and end points
Tp y0 = getVelocity((Tp)0).length();
Tp yn = getVelocity(_duration).length();
// 1/3 simpsons rule
_cached_arc_length = h / (Tp)3 * (y0 + yn + _cached_arc_length);
// update cache
_cached_resolution = resolution;
return _cached_arc_length;
}
template<typename Tp>
Tp BezierQuad<Tp>::getDistToClosestPoint(const Vector3_t &pose)
{
// get t that corresponds to point closest on bezier point
Tp t = _goldenSectionSearch(pose);
// get closest point
Vector3_t point = getPoint(t);
return (pose - point).length();
}
template<typename Tp>
Tp BezierQuad<Tp>::_goldenSectionSearch(const Vector3_t &pose)
{
Tp a, b, c, d;
a = (Tp)0; // represents most left point
b = _duration * (Tp)1; // represents most right point
c = b - (b - a) / GOLDEN_RATIO;
d = a + (b - a) / GOLDEN_RATIO;
while (fabsf(c - d) > RESOLUTION) {
if (_getDistanceSquared(c, pose) < _getDistanceSquared(d, pose)) {
b = d;
} else {
a = c;
}
c = b - (b - a) / GOLDEN_RATIO;
d = a + (b - a) / GOLDEN_RATIO;
}
return (b + a) / (Tp)2;
}
template<typename Tp>
Tp BezierQuad<Tp>::_getDistanceSquared(const Tp t, const Vector3_t &pose)
{
// get point on bezier
Vector3_t vec = getPoint(t);
// get vector from point to pose
vec = vec - pose;
// norm squared
return (vec * vec);
}
}
#endif

View File

@@ -1,210 +0,0 @@
/****************************************************************************
*
* Copyright (C) 2018 PX4 Development Team. 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 PX4 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 OWNER 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.
*
****************************************************************************/
/**
* @file BezierQuad.hpp
*
* Quadratic bezier lib
*
* A quadratic bezier function/spline is completely defined by three 3D points in space and a time scaling factor.
* pt0 and pt1 define the start and end points of the spline. ctrl point is a point in space that effects the curvature
* of the spline. The time scaling factor (= duration) defines the time it takes to travel along the spline from pt0 to
* pt1.
* A bezier spline is a continuous function from which position, velocity and acceleration can be extracted. For a given spline,
* acceleration stays constant.
*/
#pragma once
#include <matrix/math.hpp>
namespace bezier
{
template<typename Tp>
class BezierQuad
{
public:
using Vector3_t = matrix::Vector<Tp, 3>;
/**
* Empty constructor
*/
BezierQuad() :
_pt0(Vector3_t()), _ctrl(Vector3_t()), _pt1(Vector3_t()), _duration(1.0f) {}
/**
* Constructor from array
*/
BezierQuad(const Tp pt0[3], const Tp ctrl[3], const Tp pt1[3], Tp duration = 1.0f) :
_pt0(Vector3_t(pt0)), _ctrl(Vector3_t(ctrl)), _pt1(Vector3_t(pt1)), _duration(duration) {}
/**
* Constructor from vector
*/
BezierQuad(const Vector3_t &pt0, const Vector3_t &ctrl, const Vector3_t &pt1,
Tp duration = 1.0f):
_pt0(pt0), _ctrl(ctrl), _pt1(pt1), _duration(duration) {}
/*
* Get bezier points
*/
void getBezier(Vector3_t &pt0, Vector3_t &ctrl, Vector3_t &pt1);
/*
* Return pt0
*/
Vector3_t getPt0() {return _pt0;}
/*
* Return ctrl
*/
Vector3_t getCtrl() {return _ctrl;}
/*
* Return pt1
*/
Vector3_t getPt1() {return _pt1;}
/**
* Set new bezier points and duration
*/
void setBezier(const Vector3_t &pt0, const Vector3_t &ctrl, const Vector3_t &pt1,
Tp duration = (Tp)1);
/*
* Set duration
*
* @param time is the total time it takes to travel along the bezier spline.
*/
void setDuration(const Tp time) {_duration = time;}
/**
* Return point on bezier point corresponding to time t
*
* @param t is a time in seconds in between [0, duration]
* @return a point on bezier
*/
Vector3_t getPoint(const Tp t);
/*
* Distance to closest point given a position
*
* @param pose is a position in 3D space from which distance to bezier is computed.
* @return distance to closest point on bezier
*/
Tp getDistToClosestPoint(const Vector3_t &pose);
/*
* Return velocity on bezier corresponding to time t
*
* @param t is a time in seconds in between [0, duration]
* @return velocity vector at time t
*/
Vector3_t getVelocity(const Tp t);
/*
* Return acceleration on bezier corresponding to time t
*
* @return constant acceleration of bezier
*/
Vector3_t getAcceleration();
/*
* Get all states on bezier corresponding to time t
*/
void getStates(Vector3_t &point, Vector3_t &vel, Vector3_t &acc, const Tp t);
/*
* Get states on bezier which are closest to pose in space
*
* @param point is a posiiton on the spline that is closest to a given pose
* @param vel is the velocity at that given point
* @param acc is the acceleration for that spline
* @param pose represent a position in space from which closest point is computed
*/
void getStatesClosest(Vector3_t &point, Vector3_t &vel, Vector3_t &acc,
const Vector3_t pose);
/*
* Compute bezier from velocity at bezier end points and ctrl point
*
* The bezier end points are fully defined by a given control point ctrl, the duration and
* the desired velocity vectors at the end points.
*/
void setBezFromVel(const Vector3_t &ctrl, const Vector3_t &vel0, const Vector3_t &vel1,
const Tp duration = (Tp)1);
/*
* Return the arc length of a bezier spline
*
* The arc length is computed with simpsons integration.
* @param resolution in meters.
*/
Tp getArcLength(const Tp resolution);
private:
Vector3_t _pt0; /**< Bezier starting point */
Vector3_t _ctrl; /**< Bezier control point */
Vector3_t _pt1; /**< bezier end point */
Tp _duration = (Tp)1; /**< Total time to travle along spline */
Tp _cached_arc_length = (Tp)0; /**< The saved arc length of the spline */
Tp _cached_resolution = (Tp)(-1); /**< The resolution used to compute the arc length.
Negative number means that cache is not up to date. */
/*
* Golden section search
*/
Tp _goldenSectionSearch(const Vector3_t &pose);
/*
* Get squared distance from 3D pose in space and a point on bezier.
*
* @param t is the time in between [0, duration] that defines a point on the bezier.
* @param pose is a 3D pose in space.
*/
Tp _getDistanceSquared(const Tp t, const Vector3_t &pose);
};
using BezierQuad_f = BezierQuad<float>;
using BezierQuad_d = BezierQuad<double>;
}
// include implementation
#include "BezierQuad.cpp"

View File

@@ -1,39 +0,0 @@
############################################################################
#
# Copyright (c) 2018 PX4 Development Team. 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 PX4 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 OWNER 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.
#
############################################################################
px4_add_library(bezier
BezierQuad.cpp
BezierN.cpp
)
px4_add_unit_gtest(SRC BezierNTest.cpp LINKLIBS bezier)

View File

@@ -38,7 +38,7 @@ px4_add_library(FlightTaskUtility
StickYaw.cpp
)
target_link_libraries(FlightTaskUtility PUBLIC FlightTask hysteresis bezier SlewRate motion_planning mathlib)
target_link_libraries(FlightTaskUtility PUBLIC FlightTask hysteresis SlewRate motion_planning mathlib)
target_include_directories(FlightTaskUtility PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
px4_add_functional_gtest(SRC StickTiltXYTest.cpp LINKLIBS FlightTaskUtility)