diff --git a/gtsam/basis/Chebyshev2.cpp b/gtsam/basis/Chebyshev2.cpp index df7bf5f19..0e97590e0 100644 --- a/gtsam/basis/Chebyshev2.cpp +++ b/gtsam/basis/Chebyshev2.cpp @@ -21,12 +21,42 @@ namespace gtsam { -double Chebyshev2::Point(size_t N, int j, double a, double b) { +double Chebyshev2::Point(size_t N, int j) { + if (N == 1) return 0.0; assert(j >= 0 && size_t(j) < N); - const double dtheta = M_PI / (N > 1 ? (N - 1) : 1); - // We add -PI so that we get values ordered from -1 to +1 - // sin(-M_PI_2 + dtheta*j); also works - return a + (b - a) * (1. + cos(-M_PI + dtheta * j)) / 2; + const double dtheta = M_PI / (N - 1); + return -cos(dtheta * j); +} + +double Chebyshev2::Point(size_t N, int j, double a, double b) { + if (N == 1) return (a + b) / 2; + return a + (b - a) * (Point(N, j) + 1.0) / 2.0; +} + +Vector Chebyshev2::Points(size_t N) { + Vector points(N); + if (N == 1) { + points(0) = 0.0; + return points; + } + size_t half = N / 2; + const double dtheta = M_PI / (N - 1); + for (size_t j = 0; j < half; ++j) { + double c = cos(j * dtheta); + points(j) = -c; + points(N - 1 - j) = c; + } + if (N % 2 == 1) { + points(half) = 0.0; + } + return points; +} + +Vector Chebyshev2::Points(size_t N, double a, double b) { + Vector points = Points(N); + const double T1 = (a + b) / 2, T2 = (b - a) / 2; + points = T1 + (T2 * points).array(); + return points; } Weights Chebyshev2::CalculateWeights(size_t N, double x, double a, double b) { diff --git a/gtsam/basis/Chebyshev2.h b/gtsam/basis/Chebyshev2.h index 207f4b617..a66de20d5 100644 --- a/gtsam/basis/Chebyshev2.h +++ b/gtsam/basis/Chebyshev2.h @@ -51,9 +51,17 @@ class GTSAM_EXPORT Chebyshev2 : public Basis { using Parameters = Eigen::Matrix; using DiffMatrix = Eigen::Matrix; + /** + * @brief Specific Chebyshev point, within [-1,1] interval. + * + * @param N The degree of the polynomial + * @param j The index of the Chebyshev point + * @return double + */ + static double Point(size_t N, int j); + /** * @brief Specific Chebyshev point, within [a,b] interval. - * Default interval is [-1, 1] * * @param N The degree of the polynomial * @param j The index of the Chebyshev point @@ -61,24 +69,13 @@ class GTSAM_EXPORT Chebyshev2 : public Basis { * @param b Upper bound of interval (default: 1) * @return double */ - static double Point(size_t N, int j, double a = -1, double b = 1); + static double Point(size_t N, int j, double a, double b); /// All Chebyshev points - static Vector Points(size_t N) { - Vector points(N); - for (size_t j = 0; j < N; j++) { - points(j) = Point(N, j); - } - return points; - } + static Vector Points(size_t N); /// All Chebyshev points, within [a,b] interval - static Vector Points(size_t N, double a, double b) { - Vector points = Points(N); - const double T1 = (a + b) / 2, T2 = (b - a) / 2; - points = T1 + (T2 * points).array(); - return points; - } + static Vector Points(size_t N, double a, double b); /** * Evaluate Chebyshev Weights on [-1,1] at any x up to order N-1 (N values)