diff --git a/.gitignore b/.gitignore
index 3ed22f1..fdb80af 100644
--- a/.gitignore
+++ b/.gitignore
@@ -141,3 +141,13 @@ venv.bak/
.history
CLAUDE.md
+
+################################
+########### C++ BUILD ##########
+################################
+cpp/build/
+*.o
+*.a
+*.so
+*.dylib
+compile_commands.json
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 822c3f7..460cff3 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -12,14 +12,14 @@ repos:
- id: check-toml
- repo: https://github.com/astral-sh/ruff-pre-commit
- rev: 'v0.12.8'
+ rev: 'v0.15.7'
hooks:
- id: ruff
types_or: [python, pyi, jupyter]
args: [ --fix, --exit-non-zero-on-fix, --preview ]
- repo: https://github.com/pre-commit/mirrors-mypy
- rev: v1.17.1
+ rev: v1.19.1
hooks:
- id: mypy
args: [--ignore-missing-imports]
diff --git a/README.md b/README.md
index 77325c2..ac69220 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
# InterpolatePy
-
+
[](https://pepy.tech/projects/interpolatepy)
[](https://github.com/GiorgioMedico/InterpolatePy/actions/workflows/pre-commit.yml)
[](https://github.com/GiorgioMedico/InterpolatePy/actions/workflows/test.yml)
@@ -10,9 +10,9 @@
InterpolatePy provides 20+ algorithms for smooth trajectory generation with precise control over position, velocity, acceleration, and jerk. From cubic splines and B-curves to quaternion interpolation and S-curve motion profiles โ everything you need for professional motion control.
-**โก Fast:** Vectorized NumPy operations, ~1ms for 1000-point cubic splines
-**๐ฏ Precise:** Research-grade algorithms with Cยฒ continuity and bounded derivatives
-**๐ Visual:** Built-in plotting for every algorithm
+**โก Fast:** Optional C++ backend with pybind11; pure-Python fallback uses vectorized NumPy
+**๐ฏ Precise:** Research-grade algorithms with Cยฒ continuity and bounded derivatives
+**๐ Visual:** Built-in plotting for every algorithm
**๐ง Complete:** Splines, motion profiles, quaternions, and path planning in one library
---
@@ -23,7 +23,7 @@ InterpolatePy provides 20+ algorithms for smooth trajectory generation with prec
pip install InterpolatePy
```
-**Requirements:** Python โฅ3.10, NumPy โฅ2.0, SciPy โฅ1.15, Matplotlib โฅ3.10
+**Requirements:** Python โฅ3.11, NumPy โฅ2.3, SciPy โฅ1.16, Matplotlib โฅ3.10.5
Development Installation
@@ -222,12 +222,23 @@ plt.show()
## Performance & Quality
-- **Fast:** Vectorized NumPy operations, optimized algorithms
-- **Reliable:** 85%+ test coverage, continuous integration
-- **Modern:** Python 3.10+, strict typing, dataclass-based APIs
+- **Fast:** Optional C++ backend (Eigen + pybind11) for maximum performance; pure-Python fallback uses vectorized NumPy
+- **Reliable:** 85%+ test coverage, continuous integration, 142 additional C++ unit tests
+- **Modern:** Python 3.11+, strict typing, dataclass-based APIs
- **Research-grade:** Peer-reviewed algorithms from robotics literature
-**Typical Performance:**
+**C++ Backend:**
+
+InterpolatePy includes an optional compiled C++ extension for performance-critical workloads. The API is identical regardless of backend:
+
+```python
+import interpolatepy
+print(f"C++ backend: {interpolatepy.HAS_CPP}") # True if extension is available
+```
+
+Set `INTERPOLATEPY_NO_CPP=1` to force pure-Python mode.
+
+**Typical Performance (pure-Python):**
- Cubic spline (1000 points): ~1ms
- B-spline evaluation (10k points): ~5ms
- S-curve trajectory planning: ~0.5ms
@@ -256,6 +267,8 @@ ruff format interpolatepy/
ruff check interpolatepy/
mypy interpolatepy/
+# Run all pre-commit hooks
+pre-commit run --all-files
```
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
new file mode 100644
index 0000000..06aad54
--- /dev/null
+++ b/cpp/CMakeLists.txt
@@ -0,0 +1,169 @@
+cmake_minimum_required(VERSION 3.21)
+project(InterpolateCpp
+ VERSION 0.1.0
+ LANGUAGES CXX
+ DESCRIPTION "C++ port of InterpolatePy trajectory planning library"
+)
+
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+# Options
+option(INTERPOLATECPP_BUILD_TESTS "Build tests" ON)
+option(INTERPOLATECPP_BUILD_BINDINGS "Build pybind11 bindings" OFF)
+option(INTERPOLATECPP_BUILD_EXAMPLES "Build examples" OFF)
+
+# Version macros
+configure_file(
+ "${CMAKE_CURRENT_SOURCE_DIR}/cmake/version.hpp.in"
+ "${CMAKE_CURRENT_BINARY_DIR}/include/interpolatecpp/version.hpp"
+)
+
+# Dependencies via FetchContent
+include(FetchContent)
+
+FetchContent_Declare(
+ Eigen
+ GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
+ GIT_TAG 3.4.0
+ GIT_SHALLOW TRUE
+ SYSTEM
+)
+
+if(INTERPOLATECPP_BUILD_TESTS)
+ FetchContent_Declare(
+ Catch2
+ GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+ GIT_TAG v3.7.1
+ GIT_SHALLOW TRUE
+ )
+endif()
+
+if(INTERPOLATECPP_BUILD_BINDINGS)
+ FetchContent_Declare(
+ pybind11
+ GIT_REPOSITORY https://github.com/pybind/pybind11.git
+ GIT_TAG v2.13.6
+ GIT_SHALLOW TRUE
+ )
+endif()
+
+# Make Eigen available (suppress its tests and docs)
+set(BUILD_TESTING OFF CACHE BOOL "" FORCE)
+set(EIGEN_BUILD_TESTING OFF CACHE BOOL "" FORCE)
+set(EIGEN_BUILD_DOC OFF CACHE BOOL "" FORCE)
+FetchContent_MakeAvailable(Eigen)
+set(BUILD_TESTING ON CACHE BOOL "" FORCE)
+
+# Library sources
+set(INTERPOLATECPP_SOURCES
+ src/cubic_spline.cpp
+ src/cubic_smoothing_spline.cpp
+ src/cubic_spline_with_acc1.cpp
+ src/cubic_spline_with_acc2.cpp
+ src/smoothing_search.cpp
+ # Phase 2: B-spline
+ src/bspline.cpp
+ src/cubic_bspline_interpolation.cpp
+ src/bspline_interpolator.cpp
+ src/approximation_bspline.cpp
+ src/smoothing_cubic_bspline.cpp
+ # Phase 3: Motion
+ src/polynomial_trajectory.cpp
+ src/double_s_trajectory.cpp
+ src/trapezoidal_trajectory.cpp
+ src/parabolic_blend_trajectory.cpp
+ # Phase 4: Quaternion
+ src/quaternion.cpp
+ src/quaternion_spline.cpp
+ src/squad_c2.cpp
+ src/log_quaternion_interpolation.cpp
+ src/modified_log_quaternion_interpolation.cpp
+ # Phase 5: Path
+ src/linear_path.cpp
+ src/circular_path.cpp
+ src/frenet_frame.cpp
+ src/linear_traj.cpp
+)
+
+add_library(interpolatecpp ${INTERPOLATECPP_SOURCES})
+add_library(interpolatecpp::interpolatecpp ALIAS interpolatecpp)
+
+target_include_directories(interpolatecpp
+ PUBLIC
+ $
+ $
+ $
+)
+
+target_link_libraries(interpolatecpp PUBLIC Eigen3::Eigen)
+
+set_target_properties(interpolatecpp PROPERTIES
+ POSITION_INDEPENDENT_CODE ON
+ CXX_VISIBILITY_PRESET hidden
+ VISIBILITY_INLINES_HIDDEN ON
+)
+
+target_compile_options(interpolatecpp PRIVATE
+ $<$:
+ -Wall -Wextra -Wpedantic -Wconversion -Wshadow
+ >
+ $<$:
+ /W4
+ >
+)
+
+# Tests
+if(INTERPOLATECPP_BUILD_TESTS)
+ enable_testing()
+ FetchContent_MakeAvailable(Catch2)
+ add_subdirectory(tests)
+endif()
+
+# Bindings
+if(INTERPOLATECPP_BUILD_BINDINGS)
+ FetchContent_MakeAvailable(pybind11)
+ add_subdirectory(bindings)
+endif()
+
+# Examples
+if(INTERPOLATECPP_BUILD_EXAMPLES)
+ add_subdirectory(examples)
+endif()
+
+# Install rules
+include(GNUInstallDirs)
+
+install(TARGETS interpolatecpp
+ EXPORT InterpolateCppTargets
+ LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+ ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+ RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+)
+
+install(DIRECTORY include/interpolatecpp
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+)
+
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include/interpolatecpp
+ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+)
+
+install(EXPORT InterpolateCppTargets
+ FILE InterpolateCppTargets.cmake
+ NAMESPACE interpolatecpp::
+ DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/InterpolateCpp
+)
+
+include(CMakePackageConfigHelpers)
+configure_package_config_file(
+ cmake/InterpolateCppConfig.cmake.in
+ "${CMAKE_CURRENT_BINARY_DIR}/InterpolateCppConfig.cmake"
+ INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/InterpolateCpp
+)
+
+install(FILES "${CMAKE_CURRENT_BINARY_DIR}/InterpolateCppConfig.cmake"
+ DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/InterpolateCpp
+)
diff --git a/cpp/bindings/CMakeLists.txt b/cpp/bindings/CMakeLists.txt
new file mode 100644
index 0000000..482d307
--- /dev/null
+++ b/cpp/bindings/CMakeLists.txt
@@ -0,0 +1,26 @@
+pybind11_add_module(interpolatecpp_py
+ module.cpp
+ # Phase 1: Cubic Splines
+ bind_tridiagonal.cpp
+ bind_cubic_spline.cpp
+ bind_smoothing_spline.cpp
+ bind_acc_splines.cpp
+ bind_smoothing_search.cpp
+ # Phase 2: B-Splines
+ bind_bspline.cpp
+ # Phase 3: Motion Profiles
+ bind_motion.cpp
+ # Phase 4: Quaternion Interpolation
+ bind_quaternion.cpp
+ # Phase 5: Geometric Paths
+ bind_paths.cpp
+)
+
+target_link_libraries(interpolatecpp_py
+ PRIVATE
+ interpolatecpp::interpolatecpp
+)
+
+install(TARGETS interpolatecpp_py
+ LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+)
diff --git a/cpp/bindings/bind_acc_splines.cpp b/cpp/bindings/bind_acc_splines.cpp
new file mode 100644
index 0000000..fc8d68f
--- /dev/null
+++ b/cpp/bindings/bind_acc_splines.cpp
@@ -0,0 +1,96 @@
+#include
+#include
+#include
+
+#include
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::spline;
+
+void bind_acc_splines(py::module_& m) {
+ // SplineParameters
+ py::class_(m, "SplineParameters")
+ .def(py::init([](double v0, double vn, std::optional a0, std::optional an,
+ bool debug) {
+ return SplineParameters{v0, vn, a0, an, debug};
+ }),
+ py::arg("v0") = 0.0, py::arg("vn") = 0.0, py::arg("a0") = std::nullopt,
+ py::arg("an") = std::nullopt, py::arg("debug") = false)
+ .def_readwrite("v0", &SplineParameters::v0)
+ .def_readwrite("vn", &SplineParameters::vn)
+ .def_readwrite("a0", &SplineParameters::a0)
+ .def_readwrite("an", &SplineParameters::an)
+ .def_readwrite("debug", &SplineParameters::debug);
+
+ // CubicSplineWithAcceleration1
+ py::class_(m, "CubicSplineWithAcceleration1")
+ .def(py::init([](std::vector t, std::vector q, double v0, double vn,
+ double a0, double an, bool debug) {
+ return CubicSplineWithAcceleration1(t, q, v0, vn, a0, an, debug);
+ }),
+ py::arg("t_points"), py::arg("q_points"), py::arg("v0") = 0.0,
+ py::arg("vn") = 0.0, py::arg("a0") = 0.0, py::arg("an") = 0.0,
+ py::arg("debug") = false)
+ .def("evaluate",
+ py::overload_cast(&CubicSplineWithAcceleration1::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate",
+ py::overload_cast(
+ &CubicSplineWithAcceleration1::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(&CubicSplineWithAcceleration1::evaluate_velocity,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(
+ &CubicSplineWithAcceleration1::evaluate_velocity, py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(&CubicSplineWithAcceleration1::evaluate_acceleration,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(
+ &CubicSplineWithAcceleration1::evaluate_acceleration, py::const_),
+ py::arg("t"))
+ .def_property_readonly("t_points", &CubicSplineWithAcceleration1::t_points)
+ .def_property_readonly("q_points", &CubicSplineWithAcceleration1::q_points)
+ .def_property_readonly("omega", &CubicSplineWithAcceleration1::omega)
+ .def_property_readonly("n_points", &CubicSplineWithAcceleration1::n_points)
+ .def_property_readonly("n_orig", &CubicSplineWithAcceleration1::n_orig);
+
+ // CubicSplineWithAcceleration2
+ py::class_(m, "CubicSplineWithAcceleration2")
+ .def(py::init([](std::vector t, std::vector q, SplineParameters params) {
+ return CubicSplineWithAcceleration2(t, q, params);
+ }),
+ py::arg("t_points"), py::arg("q_points"),
+ py::arg("params") = SplineParameters{})
+ .def("evaluate",
+ py::overload_cast(&CubicSplineWithAcceleration2::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate",
+ py::overload_cast(
+ &CubicSplineWithAcceleration2::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(&CubicSplineWithAcceleration2::evaluate_velocity,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(
+ &CubicSplineWithAcceleration2::evaluate_velocity, py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(&CubicSplineWithAcceleration2::evaluate_acceleration,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(
+ &CubicSplineWithAcceleration2::evaluate_acceleration, py::const_),
+ py::arg("t"))
+ .def_property_readonly("has_quintic_first", &CubicSplineWithAcceleration2::has_quintic_first)
+ .def_property_readonly("has_quintic_last", &CubicSplineWithAcceleration2::has_quintic_last);
+}
diff --git a/cpp/bindings/bind_bspline.cpp b/cpp/bindings/bind_bspline.cpp
new file mode 100644
index 0000000..d178e4c
--- /dev/null
+++ b/cpp/bindings/bind_bspline.cpp
@@ -0,0 +1,116 @@
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::bspline;
+
+void bind_bspline(py::module_& m) {
+ auto bspline_mod = m.def_submodule("bspline", "B-spline interpolation algorithms");
+
+ // Parameterization enum
+ py::enum_(bspline_mod, "Parameterization")
+ .value("EquallySpaced", Parameterization::EquallySpaced)
+ .value("ChordLength", Parameterization::ChordLength)
+ .value("Centripetal", Parameterization::Centripetal);
+
+ // BSpline base class
+ py::class_(bspline_mod, "BSpline")
+ .def(py::init([](int degree, std::vector knots,
+ const Eigen::MatrixXd& control_points) {
+ return BSpline(degree, std::span(knots), control_points);
+ }),
+ py::arg("degree"), py::arg("knots"), py::arg("control_points"))
+ .def("evaluate", &BSpline::evaluate, py::arg("u"))
+ .def("evaluate_derivative", &BSpline::evaluate_derivative, py::arg("u"),
+ py::arg("order") = 1)
+ .def("generate_curve_points", &BSpline::generate_curve_points,
+ py::arg("num_points") = 100)
+ .def("find_knot_span", &BSpline::find_knot_span, py::arg("u"))
+ .def_property_readonly("degree", &BSpline::degree)
+ .def_property_readonly("knots", &BSpline::knots)
+ .def_property_readonly("control_points", &BSpline::control_points)
+ .def_property_readonly("u_min", &BSpline::u_min)
+ .def_property_readonly("u_max", &BSpline::u_max)
+ .def_property_readonly("dimension", &BSpline::dimension)
+ .def_property_readonly("n_control_points", &BSpline::n_control_points)
+ .def_static("create_uniform_knots", &BSpline::create_uniform_knots, py::arg("degree"),
+ py::arg("num_control_points"), py::arg("domain_min") = 0.0,
+ py::arg("domain_max") = 1.0)
+ .def_static("create_periodic_knots", &BSpline::create_periodic_knots,
+ py::arg("degree"), py::arg("num_control_points"),
+ py::arg("domain_min") = 0.0, py::arg("domain_max") = 1.0)
+ .def("basis_functions", &BSpline::basis_functions, py::arg("u"),
+ py::arg("span_index"))
+ .def("basis_function_derivatives", &BSpline::basis_function_derivatives,
+ py::arg("u"), py::arg("span_index"), py::arg("order"));
+
+ // CubicBSplineInterpolation
+ py::class_(bspline_mod, "CubicBSplineInterpolation")
+ .def(py::init&,
+ const std::optional&, Parameterization, bool>(),
+ py::arg("points"), py::arg("v0") = std::nullopt, py::arg("vn") = std::nullopt,
+ py::arg("method") = Parameterization::ChordLength,
+ py::arg("auto_derivatives") = false)
+ .def_property_readonly("interpolation_points",
+ &CubicBSplineInterpolation::interpolation_points)
+ .def_property_readonly("u_bars", &CubicBSplineInterpolation::u_bars);
+
+ // BSplineInterpolator
+ py::class_(bspline_mod, "BSplineInterpolator")
+ .def(py::init&,
+ const std::optional&,
+ const std::optional&,
+ const std::optional&,
+ const std::optional&, bool>(),
+ py::arg("degree"), py::arg("points"), py::arg("times") = std::nullopt,
+ py::arg("initial_velocity") = std::nullopt,
+ py::arg("final_velocity") = std::nullopt,
+ py::arg("initial_acceleration") = std::nullopt,
+ py::arg("final_acceleration") = std::nullopt, py::arg("cyclic") = false)
+ .def_property_readonly("interp_points", &BSplineInterpolator::interp_points)
+ .def_property_readonly("times", &BSplineInterpolator::times);
+
+ // ApproximationBSpline
+ py::class_(bspline_mod, "ApproximationBSpline")
+ .def(py::init&,
+ Parameterization, bool>(),
+ py::arg("points"), py::arg("num_control_points"), py::arg("degree") = 3,
+ py::arg("weights") = std::nullopt,
+ py::arg("method") = Parameterization::ChordLength, py::arg("debug") = false)
+ .def("calculate_approximation_error", &ApproximationBSpline::calculate_approximation_error)
+ .def_property_readonly("original_points", &ApproximationBSpline::original_points)
+ .def_property_readonly("original_parameters",
+ &ApproximationBSpline::original_parameters);
+
+ // BSplineParams (for SmoothingCubicBSpline config)
+ py::class_(bspline_mod, "BSplineParams")
+ .def(py::init<>())
+ .def_readwrite("mu", &BSplineParams::mu)
+ .def_readwrite("weights", &BSplineParams::weights)
+ .def_readwrite("v0", &BSplineParams::v0)
+ .def_readwrite("vn", &BSplineParams::vn)
+ .def_readwrite("method", &BSplineParams::method)
+ .def_readwrite("enforce_endpoints", &BSplineParams::enforce_endpoints)
+ .def_readwrite("auto_derivatives", &BSplineParams::auto_derivatives);
+
+ // SmoothingCubicBSpline
+ py::class_(bspline_mod, "SmoothingCubicBSpline")
+ .def(py::init(), py::arg("points"),
+ py::arg("params") = BSplineParams{})
+ .def("calculate_approximation_error",
+ &SmoothingCubicBSpline::calculate_approximation_error)
+ .def("calculate_total_error", &SmoothingCubicBSpline::calculate_total_error)
+ .def_property_readonly("approximation_points",
+ &SmoothingCubicBSpline::approximation_points)
+ .def_property_readonly("u_bars", &SmoothingCubicBSpline::u_bars)
+ .def_property_readonly("mu", &SmoothingCubicBSpline::mu)
+ .def_property_readonly("lambda_param", &SmoothingCubicBSpline::lambda_param);
+}
diff --git a/cpp/bindings/bind_cubic_spline.cpp b/cpp/bindings/bind_cubic_spline.cpp
new file mode 100644
index 0000000..633fa08
--- /dev/null
+++ b/cpp/bindings/bind_cubic_spline.cpp
@@ -0,0 +1,39 @@
+#include
+#include
+#include
+
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::spline;
+
+void bind_cubic_spline(py::module_& m) {
+ py::class_(m, "CubicSpline")
+ .def(py::init([](std::vector t, std::vector q, double v0, double vn,
+ bool debug) { return CubicSpline(t, q, v0, vn, debug); }),
+ py::arg("t_points"), py::arg("q_points"), py::arg("v0") = 0.0,
+ py::arg("vn") = 0.0, py::arg("debug") = false)
+ .def("evaluate", py::overload_cast(&CubicSpline::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate",
+ py::overload_cast(&CubicSpline::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(&CubicSpline::evaluate_velocity, py::const_), py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(&CubicSpline::evaluate_velocity,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(&CubicSpline::evaluate_acceleration, py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(&CubicSpline::evaluate_acceleration,
+ py::const_),
+ py::arg("t"))
+ .def_property_readonly("t_points", &CubicSpline::t_points)
+ .def_property_readonly("q_points", &CubicSpline::q_points)
+ .def_property_readonly("velocities", &CubicSpline::velocities)
+ .def_property_readonly("coefficients", &CubicSpline::coefficients)
+ .def_property_readonly("n_segments", &CubicSpline::n_segments);
+}
diff --git a/cpp/bindings/bind_motion.cpp b/cpp/bindings/bind_motion.cpp
new file mode 100644
index 0000000..274ecaf
--- /dev/null
+++ b/cpp/bindings/bind_motion.cpp
@@ -0,0 +1,124 @@
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::motion;
+
+void bind_motion(py::module_& m) {
+ auto motion_mod = m.def_submodule("motion", "Motion profile algorithms");
+
+ // Result types
+ py::class_(motion_mod, "TrajectoryResult")
+ .def_readonly("position", &TrajectoryResult::position)
+ .def_readonly("velocity", &TrajectoryResult::velocity)
+ .def_readonly("acceleration", &TrajectoryResult::acceleration);
+
+ py::class_(motion_mod, "FullTrajectoryResult")
+ .def_readonly("position", &FullTrajectoryResult::position)
+ .def_readonly("velocity", &FullTrajectoryResult::velocity)
+ .def_readonly("acceleration", &FullTrajectoryResult::acceleration)
+ .def_readonly("jerk", &FullTrajectoryResult::jerk);
+
+ // BoundaryCondition
+ py::class_(motion_mod, "BoundaryCondition")
+ .def(py::init<>())
+ .def_readwrite("position", &BoundaryCondition::position)
+ .def_readwrite("velocity", &BoundaryCondition::velocity)
+ .def_readwrite("acceleration", &BoundaryCondition::acceleration)
+ .def_readwrite("jerk", &BoundaryCondition::jerk);
+
+ // TimeInterval
+ py::class_(motion_mod, "TimeInterval")
+ .def(py::init<>())
+ .def_readwrite("start", &TimeInterval::start)
+ .def_readwrite("end", &TimeInterval::end)
+ .def("duration", &TimeInterval::duration);
+
+ // StateParams
+ py::class_(motion_mod, "StateParams")
+ .def(py::init([](double q0, double q1, double v0, double v1) {
+ return StateParams{q0, q1, v0, v1};
+ }),
+ py::arg("q_0"), py::arg("q_1"), py::arg("v_0") = 0.0, py::arg("v_1") = 0.0)
+ .def_readonly("q_0", &StateParams::q_0)
+ .def_readonly("q_1", &StateParams::q_1)
+ .def_readonly("v_0", &StateParams::v_0)
+ .def_readonly("v_1", &StateParams::v_1);
+
+ // TrajectoryBounds
+ py::class_(motion_mod, "TrajectoryBounds")
+ .def(py::init(), py::arg("v_bound"), py::arg("a_bound"),
+ py::arg("j_bound"))
+ .def_readonly("v_bound", &TrajectoryBounds::v_bound)
+ .def_readonly("a_bound", &TrajectoryBounds::a_bound)
+ .def_readonly("j_bound", &TrajectoryBounds::j_bound);
+
+ // PolynomialTrajectory
+ py::class_(motion_mod, "PolynomialTrajectory")
+ .def(py::init(),
+ py::arg("bc_start"), py::arg("bc_end"), py::arg("interval"), py::arg("order"))
+ .def("evaluate", &PolynomialTrajectory::evaluate, py::arg("t"))
+ .def_property_readonly("order", &PolynomialTrajectory::order)
+ .def_property_readonly("t_start", &PolynomialTrajectory::t_start)
+ .def_property_readonly("t_end", &PolynomialTrajectory::t_end)
+ .def_property_readonly("duration", &PolynomialTrajectory::duration)
+ .def_property_readonly("coefficients", &PolynomialTrajectory::coefficients)
+ .def_static("heuristic_velocities", &PolynomialTrajectory::heuristic_velocities,
+ py::arg("points"), py::arg("times"))
+ .def_static("multipoint_trajectory", &PolynomialTrajectory::multipoint_trajectory,
+ py::arg("points"), py::arg("times"), py::arg("order") = 3,
+ py::arg("v0") = 0.0, py::arg("vn") = 0.0)
+ .def_static("evaluate_multipoint", &PolynomialTrajectory::evaluate_multipoint,
+ py::arg("segments"), py::arg("t"));
+
+ // DoubleSTrajectory
+ py::class_(motion_mod, "DoubleSTrajectory")
+ .def(py::init(), py::arg("state"),
+ py::arg("bounds"))
+ .def("evaluate", &DoubleSTrajectory::evaluate, py::arg("t"))
+ .def_property_readonly("duration", &DoubleSTrajectory::duration)
+ .def("phase_durations", &DoubleSTrajectory::phase_durations);
+
+ // TrapezoidalTrajectory
+ py::class_(motion_mod, "TrapezoidalTrajectory")
+ .def(py::init(),
+ py::arg("q0"), py::arg("q1"), py::arg("amax"), py::arg("vmax"),
+ py::arg("v0") = 0.0, py::arg("v1") = 0.0, py::arg("t0") = 0.0)
+ .def(py::init([](double q0, double q1, double amax, double v0, double v1,
+ double t0, double duration) {
+ return TrapezoidalTrajectory(TrapezoidalTrajectory::DurationBased{}, q0,
+ q1, amax, v0, v1, t0, duration);
+ }),
+ py::arg("q0"), py::arg("q1"), py::arg("amax"), py::arg("v0"),
+ py::arg("v1"), py::arg("t0"), py::arg("duration"),
+ "Duration-based constructor (computes required acceleration).")
+ .def("evaluate", &TrapezoidalTrajectory::evaluate, py::arg("t"))
+ .def_property_readonly("duration", &TrapezoidalTrajectory::duration)
+ .def_property_readonly("t_start", &TrapezoidalTrajectory::t_start)
+ .def_property_readonly("t_end", &TrapezoidalTrajectory::t_end)
+ .def_static("heuristic_velocities", &TrapezoidalTrajectory::heuristic_velocities,
+ py::arg("points"), py::arg("times"), py::arg("vmax"))
+ .def_static("interpolate_waypoints", &TrapezoidalTrajectory::interpolate_waypoints,
+ py::arg("points"), py::arg("amax"), py::arg("vmax"), py::arg("v0") = 0.0,
+ py::arg("vn") = 0.0, py::arg("times") = std::vector{},
+ py::arg("velocities") = std::vector{})
+ .def_static("evaluate_multipoint", &TrapezoidalTrajectory::evaluate_multipoint,
+ py::arg("segments"), py::arg("t"));
+
+ // ParabolicBlendTrajectory
+ py::class_(motion_mod, "ParabolicBlendTrajectory")
+ .def(py::init&, const std::vector&,
+ const std::vector&>(),
+ py::arg("q"), py::arg("t"), py::arg("dt_blend"))
+ .def("evaluate", &ParabolicBlendTrajectory::evaluate, py::arg("t"))
+ .def_property_readonly("duration", &ParabolicBlendTrajectory::duration)
+ .def_property_readonly("n_waypoints", &ParabolicBlendTrajectory::n_waypoints);
+}
diff --git a/cpp/bindings/bind_paths.cpp b/cpp/bindings/bind_paths.cpp
new file mode 100644
index 0000000..e0f3e5c
--- /dev/null
+++ b/cpp/bindings/bind_paths.cpp
@@ -0,0 +1,73 @@
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::path;
+
+void bind_paths(py::module_& m) {
+ auto path_mod = m.def_submodule("path", "Geometric path algorithms");
+
+ // LinearPath
+ py::class_(path_mod, "LinearPath")
+ .def(py::init(), py::arg("pi"),
+ py::arg("pf"))
+ .def("position",
+ py::overload_cast(&LinearPath::position, py::const_), py::arg("s"))
+ .def("position",
+ py::overload_cast(&LinearPath::position, py::const_),
+ py::arg("s"))
+ .def("velocity", &LinearPath::velocity, py::arg("s"))
+ .def("acceleration", &LinearPath::acceleration, py::arg("s"))
+ .def_property_readonly("length", &LinearPath::length);
+
+ // CircularPath
+ py::class_(path_mod, "CircularPath")
+ .def(py::init(),
+ py::arg("axis"), py::arg("axis_point"), py::arg("circle_point"))
+ .def("position",
+ py::overload_cast(&CircularPath::position, py::const_), py::arg("s"))
+ .def("position",
+ py::overload_cast(&CircularPath::position, py::const_),
+ py::arg("s"))
+ .def("velocity", &CircularPath::velocity, py::arg("s"))
+ .def("acceleration", &CircularPath::acceleration, py::arg("s"))
+ .def_property_readonly("radius", &CircularPath::radius)
+ .def_property_readonly("center", &CircularPath::center);
+
+ // FrenetFrame
+ py::class_(path_mod, "FrenetFrame")
+ .def_readonly("tangent", &FrenetFrame::tangent)
+ .def_readonly("normal", &FrenetFrame::normal)
+ .def_readonly("binormal", &FrenetFrame::binormal)
+ .def_readonly("curvature", &FrenetFrame::curvature)
+ .def_readonly("torsion", &FrenetFrame::torsion);
+
+ // compute_frenet_frames
+ path_mod.def("compute_frenet_frames", &compute_frenet_frames, py::arg("curve"),
+ py::arg("s_values"));
+
+ // Helper trajectory functions
+ path_mod.def("circular_trajectory_with_derivatives",
+ &circular_trajectory_with_derivatives, py::arg("u"), py::arg("r") = 2.0);
+ path_mod.def("helicoidal_trajectory_with_derivatives",
+ &helicoidal_trajectory_with_derivatives, py::arg("u"), py::arg("r") = 2.0,
+ py::arg("d") = 0.5);
+
+ // LinearTrajResult
+ py::class_(path_mod, "LinearTrajResult")
+ .def_readonly("positions", &LinearTrajResult::positions)
+ .def_readonly("velocities", &LinearTrajResult::velocities)
+ .def_readonly("accelerations", &LinearTrajResult::accelerations);
+
+ // linear_traj
+ path_mod.def("linear_traj", &linear_traj, py::arg("p0"), py::arg("p1"), py::arg("t0"),
+ py::arg("t1"), py::arg("num_points") = 100);
+}
diff --git a/cpp/bindings/bind_quaternion.cpp b/cpp/bindings/bind_quaternion.cpp
new file mode 100644
index 0000000..e3a3965
--- /dev/null
+++ b/cpp/bindings/bind_quaternion.cpp
@@ -0,0 +1,142 @@
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::quat;
+
+void bind_quaternion(py::module_& m) {
+ auto quat_mod = m.def_submodule("quat", "Quaternion interpolation algorithms");
+
+ // Quaternion
+ py::class_(quat_mod, "Quaternion")
+ .def(py::init(), py::arg("w") = 1.0,
+ py::arg("x") = 0.0, py::arg("y") = 0.0, py::arg("z") = 0.0)
+ .def_static("identity", &Quaternion::identity)
+ .def_static("from_angle_axis", &Quaternion::from_angle_axis, py::arg("angle"),
+ py::arg("axis"))
+ .def_static("from_euler_angles", &Quaternion::from_euler_angles, py::arg("roll"),
+ py::arg("pitch"), py::arg("yaw"))
+ .def_property_readonly("w", &Quaternion::w)
+ .def_property_readonly("x", &Quaternion::x)
+ .def_property_readonly("y", &Quaternion::y)
+ .def_property_readonly("z", &Quaternion::z)
+ .def_property_readonly("vec", &Quaternion::vec)
+ .def("__mul__",
+ py::overload_cast(&Quaternion::operator*, py::const_))
+ .def("__mul__",
+ py::overload_cast(&Quaternion::operator*, py::const_))
+ .def("__rmul__",
+ [](const Quaternion& q, double s) { return q * s; })
+ .def("__add__",
+ py::overload_cast(&Quaternion::operator+, py::const_))
+ .def("__sub__",
+ py::overload_cast(&Quaternion::operator-, py::const_))
+ .def("__neg__", [](const Quaternion& q) { return -q; })
+ .def("conjugate", &Quaternion::conjugate)
+ .def("inverse", &Quaternion::inverse)
+ .def("unit", &Quaternion::unit)
+ .def("norm", &Quaternion::norm)
+ .def("norm_squared", &Quaternion::norm_squared)
+ .def("dot_product", &Quaternion::dot_product, py::arg("other"))
+ .def("to_rotation_matrix", &Quaternion::to_rotation_matrix)
+ .def("to_transformation_matrix", &Quaternion::to_transformation_matrix)
+ .def("to_axis_angle", &Quaternion::to_axis_angle)
+ .def("to_euler_angles", &Quaternion::to_euler_angles)
+ .def_static("from_rotation_matrix", &Quaternion::from_rotation_matrix,
+ py::arg("rotation_matrix"))
+ // Dynamics
+ .def("E", &Quaternion::E, py::arg("sign"))
+ .def("dot", &Quaternion::dot, py::arg("omega"), py::arg("sign"))
+ .def_static("Omega", &Quaternion::Omega, py::arg("q"), py::arg("q_dot"))
+ .def_static("slerp", &Quaternion::slerp, py::arg("q0"), py::arg("q1"), py::arg("t"))
+ .def_static("slerp_prime", &Quaternion::slerp_prime, py::arg("q0"),
+ py::arg("q1"), py::arg("t"))
+ .def_static("squad", &Quaternion::squad, py::arg("p"), py::arg("a"), py::arg("b"),
+ py::arg("q"), py::arg("t"))
+ .def_static("compute_intermediate_quaternion",
+ &Quaternion::compute_intermediate_quaternion, py::arg("q_prev"),
+ py::arg("q_curr"), py::arg("q_next"))
+ .def_static("exp", &Quaternion::exp, py::arg("q"))
+ .def_static("log", &Quaternion::log, py::arg("q"))
+ .def_static("power", &Quaternion::power, py::arg("q"), py::arg("t"));
+
+ // QuaternionSpline
+ py::enum_(quat_mod, "QuaternionSplineMethod")
+ .value("Slerp", QuaternionSpline::Method::Slerp)
+ .value("Squad", QuaternionSpline::Method::Squad)
+ .value("Auto", QuaternionSpline::Method::Auto);
+
+ py::class_(quat_mod, "QuaternionSpline")
+ .def(py::init&, const std::vector&,
+ QuaternionSpline::Method>(),
+ py::arg("time_points"), py::arg("quaternions"),
+ py::arg("method") = QuaternionSpline::Method::Auto)
+ .def("evaluate", &QuaternionSpline::evaluate, py::arg("t"))
+ .def("evaluate_velocity", &QuaternionSpline::evaluate_velocity, py::arg("t"))
+ .def("evaluate_acceleration", &QuaternionSpline::evaluate_acceleration, py::arg("t"))
+ .def_property_readonly("t_min", &QuaternionSpline::t_min)
+ .def_property_readonly("t_max", &QuaternionSpline::t_max);
+
+ // SquadC2Config
+ py::class_(quat_mod, "SquadC2Config")
+ .def(py::init<>())
+ .def_readwrite("time_points", &SquadC2Config::time_points)
+ .def_readwrite("quaternions", &SquadC2Config::quaternions)
+ .def_readwrite("normalize_quaternions", &SquadC2Config::normalize_quaternions)
+ .def_readwrite("validate_continuity", &SquadC2Config::validate_continuity);
+
+ // SquadC2
+ py::class_(quat_mod, "SquadC2")
+ .def(py::init&, const std::vector&, bool, bool>(),
+ py::arg("time_points"), py::arg("quaternions"),
+ py::arg("normalize_quaternions") = true,
+ py::arg("validate_continuity") = true)
+ .def(py::init(), py::arg("config"))
+ .def("evaluate", &SquadC2::evaluate, py::arg("t"))
+ .def("evaluate_velocity", &SquadC2::evaluate_velocity, py::arg("t"))
+ .def("evaluate_acceleration", &SquadC2::evaluate_acceleration, py::arg("t"))
+ .def_property_readonly("t_min", &SquadC2::t_min)
+ .def_property_readonly("t_max", &SquadC2::t_max)
+ .def_property_readonly("validate_continuity", &SquadC2::validate_continuity);
+
+ // LogQuaternionInterpolation
+ py::class_(quat_mod, "LogQuaternionInterpolation")
+ .def(py::init&, const std::vector&, int,
+ const std::optional&,
+ const std::optional&>(),
+ py::arg("time_points"), py::arg("quaternions"), py::arg("degree") = 3,
+ py::arg("initial_velocity") = std::nullopt,
+ py::arg("final_velocity") = std::nullopt)
+ .def("evaluate", &LogQuaternionInterpolation::evaluate, py::arg("t"))
+ .def("evaluate_velocity", &LogQuaternionInterpolation::evaluate_velocity, py::arg("t"))
+ .def("evaluate_acceleration", &LogQuaternionInterpolation::evaluate_acceleration,
+ py::arg("t"))
+ .def_property_readonly("t_min", &LogQuaternionInterpolation::t_min)
+ .def_property_readonly("t_max", &LogQuaternionInterpolation::t_max);
+
+ // ModifiedLogQuaternionInterpolation
+ py::class_(quat_mod,
+ "ModifiedLogQuaternionInterpolation")
+ .def(py::init&, const std::vector&, int, bool,
+ const std::optional&,
+ const std::optional&>(),
+ py::arg("time_points"), py::arg("quaternions"), py::arg("degree") = 3,
+ py::arg("normalize_axis") = true, py::arg("initial_velocity") = std::nullopt,
+ py::arg("final_velocity") = std::nullopt)
+ .def("evaluate", &ModifiedLogQuaternionInterpolation::evaluate, py::arg("t"))
+ .def("evaluate_velocity", &ModifiedLogQuaternionInterpolation::evaluate_velocity,
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ &ModifiedLogQuaternionInterpolation::evaluate_acceleration, py::arg("t"))
+ .def_property_readonly("t_min", &ModifiedLogQuaternionInterpolation::t_min)
+ .def_property_readonly("t_max", &ModifiedLogQuaternionInterpolation::t_max)
+ .def_property_readonly("normalize_axis",
+ &ModifiedLogQuaternionInterpolation::normalize_axis);
+}
diff --git a/cpp/bindings/bind_smoothing_search.cpp b/cpp/bindings/bind_smoothing_search.cpp
new file mode 100644
index 0000000..b2a1285
--- /dev/null
+++ b/cpp/bindings/bind_smoothing_search.cpp
@@ -0,0 +1,40 @@
+#include
+#include
+#include
+
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::spline;
+
+void bind_smoothing_search(py::module_& m) {
+ // SplineConfig
+ py::class_(m, "SplineConfig")
+ .def(py::init([](std::optional weights, double v0, double vn,
+ int max_iterations, bool debug) {
+ return SplineConfig{weights, v0, vn, max_iterations, debug};
+ }),
+ py::arg("weights") = std::nullopt, py::arg("v0") = 0.0, py::arg("vn") = 0.0,
+ py::arg("max_iterations") = 50, py::arg("debug") = false)
+ .def_readwrite("weights", &SplineConfig::weights)
+ .def_readwrite("v0", &SplineConfig::v0)
+ .def_readwrite("vn", &SplineConfig::vn)
+ .def_readwrite("max_iterations", &SplineConfig::max_iterations)
+ .def_readwrite("debug", &SplineConfig::debug);
+
+ // SmoothingSearchResult
+ py::class_(m, "SmoothingSearchResult")
+ .def_readonly("spline", &SmoothingSearchResult::spline)
+ .def_readonly("mu", &SmoothingSearchResult::mu)
+ .def_readonly("max_error", &SmoothingSearchResult::max_error)
+ .def_readonly("iterations", &SmoothingSearchResult::iterations);
+
+ // Free function
+ m.def(
+ "smoothing_spline_with_tolerance",
+ [](std::vector t, std::vector q, double tolerance,
+ const SplineConfig& config) {
+ return smoothing_spline_with_tolerance(t, q, tolerance, config);
+ },
+ py::arg("t_points"), py::arg("q_points"), py::arg("tolerance"), py::arg("config"));
+}
diff --git a/cpp/bindings/bind_smoothing_spline.cpp b/cpp/bindings/bind_smoothing_spline.cpp
new file mode 100644
index 0000000..8b6ed5f
--- /dev/null
+++ b/cpp/bindings/bind_smoothing_spline.cpp
@@ -0,0 +1,52 @@
+#include
+#include
+#include
+
+#include
+
+#include
+
+namespace py = pybind11;
+using namespace interpolatecpp::spline;
+
+void bind_smoothing_spline(py::module_& m) {
+ py::class_(m, "CubicSmoothingSpline")
+ .def(py::init([](std::vector t, std::vector q, double mu,
+ std::optional> weights, double v0, double vn,
+ bool debug) {
+ std::optional> w_span;
+ if (weights.has_value()) {
+ w_span = std::span(weights->data(), weights->size());
+ }
+ return CubicSmoothingSpline(t, q, mu, w_span, v0, vn, debug);
+ }),
+ py::arg("t_points"), py::arg("q_points"), py::arg("mu") = 0.5,
+ py::arg("weights") = std::nullopt, py::arg("v0") = 0.0, py::arg("vn") = 0.0,
+ py::arg("debug") = false)
+ .def("evaluate",
+ py::overload_cast(&CubicSmoothingSpline::evaluate, py::const_),
+ py::arg("t"))
+ .def("evaluate",
+ py::overload_cast(&CubicSmoothingSpline::evaluate,
+ py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(&CubicSmoothingSpline::evaluate_velocity, py::const_),
+ py::arg("t"))
+ .def("evaluate_velocity",
+ py::overload_cast(
+ &CubicSmoothingSpline::evaluate_velocity, py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(&CubicSmoothingSpline::evaluate_acceleration, py::const_),
+ py::arg("t"))
+ .def("evaluate_acceleration",
+ py::overload_cast(
+ &CubicSmoothingSpline::evaluate_acceleration, py::const_),
+ py::arg("t"))
+ .def_property_readonly("t_points", &CubicSmoothingSpline::t_points)
+ .def_property_readonly("q_points", &CubicSmoothingSpline::q_points)
+ .def_property_readonly("s_points", &CubicSmoothingSpline::s_points)
+ .def_property_readonly("mu", &CubicSmoothingSpline::mu)
+ .def_property_readonly("coefficients", &CubicSmoothingSpline::coefficients);
+}
diff --git a/cpp/bindings/bind_tridiagonal.cpp b/cpp/bindings/bind_tridiagonal.cpp
new file mode 100644
index 0000000..7eae957
--- /dev/null
+++ b/cpp/bindings/bind_tridiagonal.cpp
@@ -0,0 +1,13 @@
+#include
+#include
+
+#include
+
+namespace py = pybind11;
+
+void bind_tridiagonal(py::module_& m) {
+ m.def("solve_tridiagonal", &interpolatecpp::solve_tridiagonal,
+ py::arg("lower_diagonal"), py::arg("main_diagonal"),
+ py::arg("upper_diagonal"), py::arg("right_hand_side"),
+ "Solve a tridiagonal system using the Thomas algorithm");
+}
diff --git a/cpp/bindings/module.cpp b/cpp/bindings/module.cpp
new file mode 100644
index 0000000..e94242d
--- /dev/null
+++ b/cpp/bindings/module.cpp
@@ -0,0 +1,45 @@
+#include
+
+namespace py = pybind11;
+
+// Phase 1: Cubic Splines
+void bind_tridiagonal(py::module_& m);
+void bind_cubic_spline(py::module_& m);
+void bind_smoothing_spline(py::module_& m);
+void bind_acc_splines(py::module_& m);
+void bind_smoothing_search(py::module_& m);
+
+// Phase 2: B-Splines
+void bind_bspline(py::module_& m);
+
+// Phase 3: Motion Profiles
+void bind_motion(py::module_& m);
+
+// Phase 4: Quaternion Interpolation
+void bind_quaternion(py::module_& m);
+
+// Phase 5: Geometric Paths
+void bind_paths(py::module_& m);
+
+PYBIND11_MODULE(interpolatecpp_py, m) {
+ m.doc() = "C++ backend for InterpolatePy trajectory planning library";
+
+ // Phase 1
+ bind_tridiagonal(m);
+ bind_cubic_spline(m);
+ bind_smoothing_spline(m);
+ bind_acc_splines(m);
+ bind_smoothing_search(m);
+
+ // Phase 2
+ bind_bspline(m);
+
+ // Phase 3
+ bind_motion(m);
+
+ // Phase 4
+ bind_quaternion(m);
+
+ // Phase 5
+ bind_paths(m);
+}
diff --git a/cpp/cmake/InterpolateCppConfig.cmake.in b/cpp/cmake/InterpolateCppConfig.cmake.in
new file mode 100644
index 0000000..768a9a4
--- /dev/null
+++ b/cpp/cmake/InterpolateCppConfig.cmake.in
@@ -0,0 +1,5 @@
+@PACKAGE_INIT@
+
+include("${CMAKE_CURRENT_LIST_DIR}/InterpolateCppTargets.cmake")
+
+check_required_components(InterpolateCpp)
diff --git a/cpp/cmake/version.hpp.in b/cpp/cmake/version.hpp.in
new file mode 100644
index 0000000..e0a632f
--- /dev/null
+++ b/cpp/cmake/version.hpp.in
@@ -0,0 +1,6 @@
+#pragma once
+
+#define INTERPOLATECPP_VERSION_MAJOR @PROJECT_VERSION_MAJOR@
+#define INTERPOLATECPP_VERSION_MINOR @PROJECT_VERSION_MINOR@
+#define INTERPOLATECPP_VERSION_PATCH @PROJECT_VERSION_PATCH@
+#define INTERPOLATECPP_VERSION "@PROJECT_VERSION@"
diff --git a/cpp/examples/CMakeLists.txt b/cpp/examples/CMakeLists.txt
new file mode 100644
index 0000000..b4a79eb
--- /dev/null
+++ b/cpp/examples/CMakeLists.txt
@@ -0,0 +1,29 @@
+set(EXAMPLE_SOURCES
+ # Phase 1: Cubic Splines
+ cubic_spline_example.cpp
+ cubic_spline_acc1_example.cpp
+ cubic_spline_acc2_example.cpp
+ cubic_smoothing_example.cpp
+ # Phase 2: B-Splines
+ bspline_example.cpp
+ bspline_cubic_example.cpp
+ bspline_interpolator_example.cpp
+ bspline_approx_smooth_example.cpp
+ # Phase 3: Motion Profiles
+ trapezoidal_example.cpp
+ polynomial_example.cpp
+ double_s_example.cpp
+ parabolic_linear_example.cpp
+ # Phase 4: Quaternion
+ quaternion_example.cpp
+ # Phase 5: Paths & Concepts
+ paths_example.cpp
+ concepts_example.cpp
+)
+
+foreach(source IN LISTS EXAMPLE_SOURCES)
+ get_filename_component(name ${source} NAME_WE)
+ add_executable(${name} ${source})
+ target_link_libraries(${name} PRIVATE interpolatecpp::interpolatecpp)
+ target_include_directories(${name} PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/shared)
+endforeach()
diff --git a/cpp/examples/bspline_approx_smooth_example.cpp b/cpp/examples/bspline_approx_smooth_example.cpp
new file mode 100644
index 0000000..b449cfe
--- /dev/null
+++ b/cpp/examples/bspline_approx_smooth_example.cpp
@@ -0,0 +1,387 @@
+/// B-spline approximation and smoothing example -- C++ port of
+/// examples/b_spline_approx_ex.py AND examples/b_spline_smooth_ex.py
+///
+/// Part 1: ApproximationBSpline -- least-squares fitting with varying control
+/// point counts and degrees.
+/// Part 2: SmoothingCubicBSpline -- smoothing with different lambda/mu values
+/// (Example 8.12).
+
+#include
+#include
+#include
+#include
+
+#include "example_utils.hpp"
+
+#include
+#include
+#include
+#include
+#include
+
+namespace ex = interpolatecpp::examples;
+using namespace interpolatecpp::bspline;
+
+// ============================================================================
+// Part 1: Approximation B-spline
+// ============================================================================
+
+/// Generate sample points from a B-spline defined by the Example 8.10 control points.
+static Eigen::MatrixXd generate_sample_points(int num_samples) {
+ // Control points from Example 8.10
+ Eigen::MatrixXd control_points(10, 2);
+ control_points << 137, 229,
+ 101, 201,
+ 177, 121,
+ 93, 44,
+ 62, 203,
+ 49, 272,
+ 104, 402,
+ 141, 277,
+ 147, 258,
+ 138, 231;
+
+ const int degree = 3;
+ const Eigen::VectorXd knots = BSpline::create_uniform_knots(degree, 10);
+
+ const BSpline spline(degree,
+ std::span(knots.data(), static_cast(knots.size())),
+ control_points);
+
+ const auto [params, curve_pts] = spline.generate_curve_points(num_samples);
+ return curve_pts;
+}
+
+/// Example from Section 8.5: Approximation with different configurations.
+static void example_approximation() {
+ ex::print_header("Part 1 -- B-spline Approximation (Section 8.5)");
+
+ // Print the source control points
+ std::cout << "Control points from Example 8.10:\n";
+ const std::vector> cp_coords = {
+ {137, 229}, {101, 201}, {177, 121}, {93, 44}, {62, 203},
+ {49, 272}, {104, 402}, {141, 277}, {147, 258}, {138, 231}
+ };
+ for (size_t i = 0; i < cp_coords.size(); ++i) {
+ std::cout << " P" << i << ": (" << cp_coords[i].first
+ << ", " << cp_coords[i].second << ")\n";
+ }
+
+ // Generate sample points
+ const Eigen::MatrixXd sample_points = generate_sample_points(84);
+ std::cout << "\nGenerated " << sample_points.rows() << " sample points.\n";
+
+ // Test cases: {num_cps, degree, title}
+ struct TestCase {
+ int num_cps;
+ int degree;
+ std::string title;
+ };
+ const std::vector test_cases = {
+ {10, 3, "Cubic (p=3) with 10 control points"},
+ {10, 4, "Quartic (p=4) with 10 control points"},
+ {20, 3, "Cubic (p=3) with 20 control points"},
+ };
+
+ for (size_t i = 0; i < test_cases.size(); ++i) {
+ const auto& tc = test_cases[i];
+ ex::print_separator('=');
+ std::cout << "Test case " << (i + 1) << ": " << tc.title << "\n\n";
+
+ const ApproximationBSpline approx(
+ sample_points, tc.num_cps, tc.degree);
+
+ const double error = approx.calculate_approximation_error();
+ ex::print_value("Approximation error", error, 2);
+ ex::print_value("Num control points", static_cast(approx.n_control_points()), 0);
+ ex::print_value("Degree", static_cast(approx.degree()), 0);
+
+ std::cout << "\n";
+ ex::print_vector("Knot vector", approx.knots());
+ ex::print_matrix("Control points", approx.control_points());
+
+ // Generate and print a few curve points
+ const auto [params, curve_pts] = approx.generate_curve_points(10);
+ std::cout << "\nSample curve points:\n";
+ const int w = 14;
+ std::cout << std::right << std::setw(w) << "u"
+ << std::setw(w) << "X" << std::setw(w) << "Y" << "\n";
+ ex::print_separator('-', 3 * w);
+ for (Eigen::Index j = 0; j < curve_pts.rows(); ++j) {
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << params(j)
+ << std::setw(w) << curve_pts(j, 0)
+ << std::setw(w) << curve_pts(j, 1) << "\n";
+ }
+ }
+}
+
+/// Demonstrate degree comparison on a heart-shaped curve.
+static void example_degree_comparison() {
+ ex::print_header("Approximation -- Degree Comparison");
+
+ // Generate heart-shaped points
+ const int n = 100;
+ Eigen::MatrixXd heart_points(n, 2);
+ for (int i = 0; i < n; ++i) {
+ const double t = 2.0 * M_PI * static_cast(i) / n;
+ heart_points(i, 0) = 16.0 * std::pow(std::sin(t), 3) * 10.0 + 150.0;
+ heart_points(i, 1) = (13.0 * std::cos(t) - 5.0 * std::cos(2.0 * t)
+ - 2.0 * std::cos(3.0 * t) - std::cos(4.0 * t)) * 10.0 + 150.0;
+ }
+
+ std::cout << "Heart shape: " << heart_points.rows() << " sample points\n\n";
+
+ const int num_cps = 12;
+ const std::vector degrees = {2, 3, 4};
+
+ for (const int degree : degrees) {
+ const ApproximationBSpline approx(heart_points, num_cps, degree);
+ const double error = approx.calculate_approximation_error();
+
+ std::cout << "Degree " << degree
+ << " (n_cp=" << num_cps << "): error = "
+ << std::fixed << std::setprecision(2) << error << "\n";
+ }
+}
+
+/// Demonstrate control-point-count comparison on a spiral.
+static void example_cp_count_comparison() {
+ ex::print_header("Approximation -- Control Point Count Comparison");
+
+ // Generate spiral points
+ const int n = 100;
+ Eigen::MatrixXd spiral_points(n, 2);
+ for (int i = 0; i < n; ++i) {
+ const double t = 6.0 * M_PI * static_cast(i) / (n - 1);
+ const double r = 5.0 + 15.0 * t;
+ spiral_points(i, 0) = r * std::cos(t) + 150.0;
+ spiral_points(i, 1) = r * std::sin(t) + 150.0;
+ }
+
+ std::cout << "Spiral: " << spiral_points.rows() << " sample points\n\n";
+
+ const int degree = 3;
+ const std::vector cp_counts = {8, 15, 25};
+
+ for (const int num_cps : cp_counts) {
+ const ApproximationBSpline approx(spiral_points, num_cps, degree);
+ const double error = approx.calculate_approximation_error();
+
+ std::cout << "CPs = " << std::setw(3) << num_cps
+ << " (degree=" << degree << "): error = "
+ << std::fixed << std::setprecision(2) << error << "\n";
+ }
+}
+
+/// Demonstrate weighted approximation.
+static void example_weighted_approximation() {
+ ex::print_header("Approximation -- Weighted Fit");
+
+ // Circle with noise
+ const int n = 60;
+ Eigen::MatrixXd circle_points(n, 2);
+ for (int i = 0; i < n; ++i) {
+ const double t = 2.0 * M_PI * static_cast(i) / n;
+ circle_points(i, 0) = 100.0 * std::cos(t) + 150.0;
+ circle_points(i, 1) = 100.0 * std::sin(t) + 150.0;
+ }
+
+ const int num_cps = 10;
+ const int degree = 3;
+
+ // Uniform weights
+ const ApproximationBSpline approx_uniform(circle_points, num_cps, degree);
+ const double error_uniform = approx_uniform.calculate_approximation_error();
+
+ // Custom weights: emphasize first half of the points
+ Eigen::VectorXd weights = Eigen::VectorXd::Ones(n);
+ for (int i = 0; i < n / 2; ++i) {
+ weights(i) = 5.0;
+ }
+
+ const ApproximationBSpline approx_weighted(
+ circle_points, num_cps, degree, weights);
+ const double error_weighted = approx_weighted.calculate_approximation_error();
+
+ std::cout << "Circle approximation (" << n << " points, " << num_cps << " CPs):\n";
+ ex::print_value("Uniform weights error", error_uniform, 2);
+ ex::print_value("Weighted (first half emphasized) error", error_weighted, 2);
+}
+
+// ============================================================================
+// Part 2: Smoothing Cubic B-spline
+// ============================================================================
+
+/// Example 8.12: Smoothing B-spline with different lambda values.
+static void example_8_12() {
+ ex::print_header("Part 2 -- Smoothing Cubic B-spline (Example 8.12)");
+
+ // Points from the example
+ Eigen::MatrixXd points(6, 3);
+ points << 0, 0, 0,
+ 1, 2, 1,
+ 2, 3, 0,
+ 4, 3, 0,
+ 5, 2, 2,
+ 6, 0, 2;
+
+ std::cout << "Approximation points:\n";
+ ex::print_matrix("Points", points);
+
+ // Test different lambda values
+ const std::vector lambda_values = {1e-4, 1e-5, 1e-6};
+
+ for (const double lambda_val : lambda_values) {
+ ex::print_separator('=');
+
+ // Convert lambda to mu: lambda = (1 - mu) / (6 * mu) => mu = 1 / (6*lambda + 1)
+ const double mu = 1.0 / (6.0 * lambda_val + 1.0);
+
+ std::cout << "Lambda = " << std::scientific << std::setprecision(6) << lambda_val
+ << ", Mu = " << std::fixed << std::setprecision(6) << mu << "\n\n";
+
+ BSplineParams params;
+ params.mu = mu;
+ params.method = Parameterization::ChordLength;
+ params.enforce_endpoints = true;
+ params.auto_derivatives = true;
+
+ const SmoothingCubicBSpline spline(points, params);
+
+ ex::print_value("Mu (stored)", spline.mu());
+ ex::print_value("Lambda (stored)", spline.lambda_param());
+ ex::print_value("Degree", static_cast(spline.degree()), 0);
+ ex::print_value("Num control points", static_cast(spline.n_control_points()), 0);
+
+ std::cout << "\n";
+ ex::print_matrix("Control points", spline.control_points());
+
+ // Approximation errors
+ const Eigen::VectorXd errors = spline.calculate_approximation_error();
+ std::cout << "\n";
+ ex::print_vector("Per-point errors", errors);
+
+ const double total_error = spline.calculate_total_error();
+ ex::print_value("Total error", total_error);
+
+ // Evaluate along the curve
+ std::cout << "\nCurve samples:\n";
+ const int w = 14;
+ std::cout << std::right
+ << std::setw(w) << "u"
+ << std::setw(w) << "X"
+ << std::setw(w) << "Y"
+ << std::setw(w) << "Z" << "\n";
+ ex::print_separator('-', 4 * w);
+
+ const int n_samples = 10;
+ const double u_start = spline.u_min();
+ const double u_end = spline.u_max();
+ for (int i = 0; i <= n_samples; ++i) {
+ const double u = u_start + (u_end - u_start) * static_cast(i) / n_samples;
+ const Eigen::VectorXd pt = spline.evaluate(u);
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << u
+ << std::setw(w) << pt(0)
+ << std::setw(w) << pt(1)
+ << std::setw(w) << pt(2) << "\n";
+ }
+ }
+}
+
+/// Demonstrate smoothing with different mu values on 2D data.
+static void example_smoothing_mu_comparison() {
+ ex::print_header("Smoothing -- Mu Comparison (2D)");
+
+ // Simple 2D points
+ Eigen::MatrixXd points(8, 2);
+ points << 0, 0,
+ 1, 3,
+ 2, 1,
+ 3, 4,
+ 4, 2,
+ 5, 5,
+ 6, 1,
+ 7, 0;
+
+ const std::vector mu_values = {0.1, 0.5, 0.9};
+
+ for (const double mu : mu_values) {
+ BSplineParams params;
+ params.mu = mu;
+ params.method = Parameterization::ChordLength;
+ params.enforce_endpoints = true;
+ params.auto_derivatives = true;
+
+ const SmoothingCubicBSpline spline(points, params);
+
+ const double total_error = spline.calculate_total_error();
+ std::cout << "Mu = " << std::fixed << std::setprecision(2) << mu
+ << ": total_error = " << std::setprecision(6) << total_error
+ << ", n_cp = " << spline.n_control_points() << "\n";
+ }
+}
+
+/// Demonstrate smoothing with explicit endpoint derivatives.
+static void example_smoothing_with_derivatives() {
+ ex::print_header("Smoothing -- Explicit Endpoint Derivatives");
+
+ Eigen::MatrixXd points(6, 3);
+ points << 0, 0, 0,
+ 1, 2, 1,
+ 2, 3, 0,
+ 4, 3, 0,
+ 5, 2, 2,
+ 6, 0, 2;
+
+ Eigen::VectorXd v0(3), vn(3);
+ v0 << 4.43, 8.87, 4.43;
+ vn << 4.85, -9.71, 0.0;
+
+ BSplineParams params;
+ params.mu = 0.999;
+ params.v0 = v0;
+ params.vn = vn;
+ params.method = Parameterization::ChordLength;
+ params.enforce_endpoints = true;
+ params.auto_derivatives = false;
+
+ const SmoothingCubicBSpline spline(points, params);
+
+ std::cout << "With explicit endpoint derivatives:\n";
+ ex::print_vector("v0", v0);
+ ex::print_vector("vn", vn);
+ ex::print_value("Mu", spline.mu());
+ ex::print_value("Total error", spline.calculate_total_error());
+ std::cout << "\n";
+ ex::print_matrix("Control points", spline.control_points());
+
+ // Second derivative magnitudes at a few points
+ ex::print_separator();
+ std::cout << "Second derivative magnitude along curve:\n\n";
+ const int n_samples = 8;
+ for (int i = 0; i <= n_samples; ++i) {
+ const double u = spline.u_min()
+ + (spline.u_max() - spline.u_min()) * static_cast(i) / n_samples;
+ const Eigen::VectorXd d2 = spline.evaluate_derivative(u, 2);
+ const double mag = d2.norm();
+ std::cout << " u = " << std::fixed << std::setprecision(4) << u
+ << ": ||s''(u)|| = " << std::setprecision(4) << mag << "\n";
+ }
+}
+
+int main() {
+ // Part 1: Approximation
+ example_approximation();
+ example_degree_comparison();
+ example_cp_count_comparison();
+ example_weighted_approximation();
+
+ // Part 2: Smoothing
+ example_8_12();
+ example_smoothing_mu_comparison();
+ example_smoothing_with_derivatives();
+
+ return 0;
+}
diff --git a/cpp/examples/bspline_cubic_example.cpp b/cpp/examples/bspline_cubic_example.cpp
new file mode 100644
index 0000000..7a12cf0
--- /dev/null
+++ b/cpp/examples/bspline_cubic_example.cpp
@@ -0,0 +1,182 @@
+/// Cubic B-spline interpolation example -- C++ port of examples/b_spline_cubic_ex.py
+///
+/// Demonstrates CubicBSplineInterpolation through 3D points (Example 8.8)
+/// with chord-length parameterization and auto-derivative computation.
+
+#include
+#include
+
+#include "example_utils.hpp"
+
+#include
+
+namespace ex = interpolatecpp::examples;
+using namespace interpolatecpp::bspline;
+
+/// Example 8.8: 3D cubic B-spline interpolation.
+static void example_8_8() {
+ ex::print_header("Example 8.8 -- 3D Cubic B-spline Interpolation");
+
+ // Define the interpolation points from the example
+ Eigen::MatrixXd points(10, 3);
+ points << 83, -54, 119,
+ -64, 10, 124,
+ 42, 79, 226,
+ -98, 23, 222,
+ -13, 125, 102,
+ 140, 81, 92,
+ 43, 32, 92,
+ -65, -17, 134,
+ -45, -89, 182,
+ 71, 90, 192;
+
+ // Create cubic B-spline interpolation with chord-length parameterization
+ const CubicBSplineInterpolation interpolation(
+ points,
+ std::nullopt, // v0 (auto)
+ std::nullopt, // vn (auto)
+ Parameterization::ChordLength,
+ true // auto_derivatives
+ );
+
+ // Print basic info
+ std::cout << "Cubic B-spline Interpolation Information:\n";
+ ex::print_value("Number of interpolation points",
+ static_cast(points.rows()), 0);
+ ex::print_value("Degree", static_cast(interpolation.degree()), 0);
+ ex::print_value("Dimension", static_cast(interpolation.dimension()), 0);
+ ex::print_value("u_min", interpolation.u_min());
+ ex::print_value("u_max", interpolation.u_max());
+
+ // Print parameter values (u-bars)
+ std::cout << "\n";
+ ex::print_vector("Parameter values (u_bars)", interpolation.u_bars());
+
+ // Print the knot vector
+ ex::print_vector("Knot vector", interpolation.knots());
+
+ // Print start/end derivatives
+ std::cout << "\n";
+ ex::print_vector("Start derivative (v0)", interpolation.start_derivative());
+ ex::print_vector("End derivative (vn)", interpolation.end_derivative());
+
+ // Print control points
+ std::cout << "\n";
+ ex::print_matrix("Control points", interpolation.control_points());
+
+ // Print interpolation points
+ std::cout << "\n";
+ ex::print_matrix("Interpolation points", interpolation.interpolation_points());
+
+ // Verify interpolation: evaluate at each parameter value
+ ex::print_separator('=');
+ std::cout << "Interpolation verification (evaluate at each u_bar):\n\n";
+
+ const int w = 12;
+ std::cout << std::right
+ << std::setw(8) << "u_bar"
+ << std::setw(w) << "X_eval"
+ << std::setw(w) << "Y_eval"
+ << std::setw(w) << "Z_eval"
+ << std::setw(w) << "X_orig"
+ << std::setw(w) << "Y_orig"
+ << std::setw(w) << "Z_orig" << "\n";
+ ex::print_separator('-', 8 + 6 * w);
+
+ for (Eigen::Index i = 0; i < points.rows(); ++i) {
+ const double u = interpolation.u_bars()(i);
+ const Eigen::VectorXd evaluated = interpolation.evaluate(u);
+
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(8) << u
+ << std::setw(w) << evaluated(0)
+ << std::setw(w) << evaluated(1)
+ << std::setw(w) << evaluated(2)
+ << std::setw(w) << points(i, 0)
+ << std::setw(w) << points(i, 1)
+ << std::setw(w) << points(i, 2) << "\n";
+ }
+
+ // Evaluate derivatives at a few points along the curve
+ ex::print_separator('=');
+ std::cout << "Derivative evaluation at sample parameter values:\n\n";
+
+ const int n_samples = 5;
+ const double u_start = interpolation.u_min();
+ const double u_end = interpolation.u_max();
+
+ for (int i = 0; i <= n_samples; ++i) {
+ const double u = u_start + (u_end - u_start) * static_cast(i) / n_samples;
+ const Eigen::VectorXd pos = interpolation.evaluate(u);
+ const Eigen::VectorXd vel = interpolation.evaluate_derivative(u, 1);
+ const Eigen::VectorXd acc = interpolation.evaluate_derivative(u, 2);
+
+ std::cout << "u = " << std::fixed << std::setprecision(4) << u << ":\n";
+ ex::print_vector(" Position", pos);
+ ex::print_vector(" Velocity", vel);
+ ex::print_vector(" Acceleration", acc);
+ std::cout << "\n";
+ }
+}
+
+/// Demonstrate with explicit endpoint derivatives.
+static void example_with_derivatives() {
+ ex::print_header("Cubic B-spline with Explicit Derivatives");
+
+ // Simple 2D points
+ Eigen::MatrixXd points(5, 2);
+ points << 0, 0,
+ 1, 2,
+ 3, 3,
+ 5, 1,
+ 6, 0;
+
+ // Specify endpoint derivatives
+ Eigen::VectorXd v0(2);
+ v0 << 1.0, 3.0;
+
+ Eigen::VectorXd vn(2);
+ vn << 1.0, -2.0;
+
+ const CubicBSplineInterpolation interpolation(
+ points, v0, vn,
+ Parameterization::ChordLength,
+ false // not auto_derivatives
+ );
+
+ std::cout << "Interpolation with explicit derivatives:\n";
+ ex::print_value("Number of points", static_cast(points.rows()), 0);
+ ex::print_value("Degree", static_cast(interpolation.degree()), 0);
+ ex::print_vector("Start derivative (v0)", interpolation.start_derivative());
+ ex::print_vector("End derivative (vn)", interpolation.end_derivative());
+
+ std::cout << "\n";
+ ex::print_vector("u_bars", interpolation.u_bars());
+ ex::print_vector("Knots", interpolation.knots());
+
+ std::cout << "\n";
+ ex::print_matrix("Control points", interpolation.control_points());
+
+ // Generate curve points and print a summary
+ const auto [params, curve_pts] = interpolation.generate_curve_points(15);
+ std::cout << "\nCurve samples:\n";
+
+ const int w = 14;
+ std::cout << std::right << std::setw(w) << "u"
+ << std::setw(w) << "X" << std::setw(w) << "Y" << "\n";
+ ex::print_separator('-', 3 * w);
+
+ for (Eigen::Index i = 0; i < curve_pts.rows(); ++i) {
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << params(i)
+ << std::setw(w) << curve_pts(i, 0)
+ << std::setw(w) << curve_pts(i, 1) << "\n";
+ }
+}
+
+int main() {
+ example_8_8();
+ example_with_derivatives();
+
+ return 0;
+}
diff --git a/cpp/examples/bspline_example.cpp b/cpp/examples/bspline_example.cpp
new file mode 100644
index 0000000..0d855c0
--- /dev/null
+++ b/cpp/examples/bspline_example.cpp
@@ -0,0 +1,229 @@
+/// B-spline example -- C++ port of examples/b_spline_ex.py
+///
+/// Demonstrates BSpline construction, basis function evaluation,
+/// derivative computation, curve point generation, and 3D curves.
+
+#include
+
+#include "example_utils.hpp"
+
+#include
+#include
+#include
+
+namespace ex = interpolatecpp::examples;
+using namespace interpolatecpp::bspline;
+
+/// Create the 2D B-spline from the document example.
+static BSpline create_example_bspline() {
+ const int degree = 3;
+
+ Eigen::MatrixXd control_points(7, 2);
+ control_points << 1, 2,
+ 2, 3,
+ 3, -3,
+ 4, 4,
+ 5, 5,
+ 6, -5,
+ 7, -6;
+
+ const std::vector knots = {0, 0, 0, 0, 1, 2, 4, 7, 7, 7, 7};
+
+ return BSpline(degree, std::span(knots), control_points);
+}
+
+/// Demonstrate basic B-spline evaluation and basis functions.
+static void demonstrate_basic_bspline() {
+ ex::print_header("Basic B-spline Evaluation");
+
+ const auto bspline = create_example_bspline();
+
+ std::cout << "B-spline properties:\n";
+ ex::print_value("Degree", static_cast(bspline.degree()), 0);
+ ex::print_value("Number of control points", static_cast(bspline.n_control_points()), 0);
+ ex::print_value("Dimension", static_cast(bspline.dimension()), 0);
+ ex::print_value("u_min", bspline.u_min());
+ ex::print_value("u_max", bspline.u_max());
+
+ // Evaluate at u = 1.5
+ const double u_value = 1.5;
+ const Eigen::VectorXd point = bspline.evaluate(u_value);
+
+ std::cout << "\nPoint at u = " << u_value << ":\n";
+ ex::print_vector("Position", point);
+
+ // Basis functions at u = 1.5
+ const int span = bspline.find_knot_span(u_value);
+ const Eigen::VectorXd basis = bspline.basis_functions(u_value, span);
+
+ std::cout << "\nKnot span index at u = " << u_value << ": " << span << "\n";
+ std::cout << "Non-zero basis functions:\n";
+ for (Eigen::Index i = 0; i < basis.size(); ++i) {
+ std::cout << " B^" << bspline.degree() << "_"
+ << (span - bspline.degree() + static_cast(i))
+ << " = " << std::fixed << std::setprecision(4) << basis(i) << "\n";
+ }
+
+ // Print knot vector
+ std::cout << "\n";
+ ex::print_vector("Knot vector", bspline.knots());
+ ex::print_matrix("Control points", bspline.control_points());
+}
+
+/// Example B.6: Basis function derivatives at u = 4.5.
+static void example_b6() {
+ ex::print_header("Example B.6 -- Basis Function Derivatives");
+
+ const int degree = 3;
+ const std::vector knots = {0, 0, 0, 0, 1, 2, 4, 7, 7, 7, 7};
+
+ // Dummy control points (basis functions don't depend on them)
+ Eigen::MatrixXd control_points = Eigen::MatrixXd::Zero(7, 2);
+
+ const BSpline bspline(degree, std::span(knots), control_points);
+
+ const double u_value = 4.5;
+ const int span = bspline.find_knot_span(u_value);
+ std::cout << "For u = " << u_value << ", the knot span index is: " << span << "\n";
+
+ // Calculate derivatives up to order 3
+ const Eigen::MatrixXd derivatives = bspline.basis_function_derivatives(u_value, span, 3);
+
+ std::cout << "\nBasis function values and derivatives at u = 4.5:\n";
+ ex::print_separator('-', 80);
+
+ for (int k = 0; k < 4; ++k) {
+ std::cout << "Ders[" << k << "]: ";
+ for (int j = 0; j < 4; ++j) {
+ if (j > 0) std::cout << ", ";
+ std::cout << std::fixed << std::setprecision(4) << derivatives(k, j);
+ }
+ std::cout << "\n";
+ }
+
+ std::cout << "\nWhich correspond to:\n";
+ ex::print_separator('-', 80);
+
+ const std::vector labels = {"B_i^3 ", "B_i^3(1) ", "B_i^3(2) ", "B_i^3(3) "};
+ for (int k = 0; k < 4; ++k) {
+ std::cout << labels[static_cast(k)] << ": ";
+ for (int j = 0; j < 4; ++j) {
+ const int idx = span - degree + j;
+ if (j > 0) std::cout << ", ";
+ std::cout << "B_" << idx << " = " << std::fixed << std::setprecision(4)
+ << derivatives(k, j);
+ }
+ std::cout << "\n";
+ }
+
+ std::cout << "\nAll the other terms B_j^3(k) are null.\n";
+}
+
+/// Demonstrate curve point generation.
+static void demonstrate_curve_generation() {
+ ex::print_header("Curve Point Generation");
+
+ const auto bspline = create_example_bspline();
+
+ // Generate curve points
+ const int num_points = 20;
+ const auto [params, curve_points] = bspline.generate_curve_points(num_points);
+
+ std::cout << "Generated " << curve_points.rows() << " curve points:\n\n";
+
+ const int w = 14;
+ std::cout << std::right << std::setw(w) << "u"
+ << std::setw(w) << "X" << std::setw(w) << "Y" << "\n";
+ ex::print_separator('-', 3 * w);
+
+ for (Eigen::Index i = 0; i < curve_points.rows(); ++i) {
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << params(i)
+ << std::setw(w) << curve_points(i, 0)
+ << std::setw(w) << curve_points(i, 1) << "\n";
+ }
+}
+
+/// Demonstrate 3D B-spline with uniform knots.
+static void demonstrate_3d_bspline() {
+ ex::print_header("3D B-spline Curve");
+
+ const int degree = 3;
+
+ Eigen::MatrixXd control_points(6, 3);
+ control_points << 0, 0, 0,
+ 1, 1, 2,
+ 2, -1, 1,
+ 3, 0, 3,
+ 4, 2, 0,
+ 5, 0, 1;
+
+ const Eigen::VectorXd knots = BSpline::create_uniform_knots(degree, 6);
+
+ const BSpline bspline(degree, std::span(knots.data(), static_cast(knots.size())),
+ control_points);
+
+ std::cout << "B-spline properties:\n";
+ ex::print_value("Degree", static_cast(bspline.degree()), 0);
+ ex::print_value("Number of control points", static_cast(bspline.n_control_points()), 0);
+ ex::print_value("Dimension", static_cast(bspline.dimension()), 0);
+ ex::print_value("u_min", bspline.u_min());
+ ex::print_value("u_max", bspline.u_max());
+
+ std::cout << "\n";
+ ex::print_vector("Uniform knot vector", bspline.knots());
+ ex::print_matrix("Control points (3D)", bspline.control_points());
+
+ // Evaluate at midpoint
+ const double u_mid = (bspline.u_min() + bspline.u_max()) / 2.0;
+ const Eigen::VectorXd mid_point = bspline.evaluate(u_mid);
+ std::cout << "\nPoint at u_mid = " << std::fixed << std::setprecision(4) << u_mid << ":\n";
+ ex::print_vector("Position", mid_point);
+
+ // First derivative at midpoint
+ const Eigen::VectorXd deriv1 = bspline.evaluate_derivative(u_mid, 1);
+ ex::print_vector("1st derivative", deriv1);
+
+ // Second derivative at midpoint
+ const Eigen::VectorXd deriv2 = bspline.evaluate_derivative(u_mid, 2);
+ ex::print_vector("2nd derivative", deriv2);
+
+ // Generate 3D curve points
+ const auto [params, curve_pts] = bspline.generate_curve_points(10);
+ std::cout << "\nSample 3D curve points:\n";
+ const int w = 14;
+ std::cout << std::right << std::setw(w) << "u"
+ << std::setw(w) << "X" << std::setw(w) << "Y" << std::setw(w) << "Z" << "\n";
+ ex::print_separator('-', 4 * w);
+ for (Eigen::Index i = 0; i < curve_pts.rows(); ++i) {
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << params(i)
+ << std::setw(w) << curve_pts(i, 0)
+ << std::setw(w) << curve_pts(i, 1)
+ << std::setw(w) << curve_pts(i, 2) << "\n";
+ }
+}
+
+/// Demonstrate periodic knot generation.
+static void demonstrate_periodic_knots() {
+ ex::print_header("Periodic Knots");
+
+ const int degree = 3;
+ const int n_cp = 6;
+
+ const Eigen::VectorXd uniform_knots = BSpline::create_uniform_knots(degree, n_cp, 0.0, 1.0);
+ const Eigen::VectorXd periodic_knots = BSpline::create_periodic_knots(degree, n_cp, 0.0, 1.0);
+
+ ex::print_vector("Uniform knots (p=3, n_cp=6)", uniform_knots);
+ ex::print_vector("Periodic knots (p=3, n_cp=6)", periodic_knots);
+}
+
+int main() {
+ demonstrate_basic_bspline();
+ example_b6();
+ demonstrate_curve_generation();
+ demonstrate_3d_bspline();
+ demonstrate_periodic_knots();
+
+ return 0;
+}
diff --git a/cpp/examples/bspline_interpolator_example.cpp b/cpp/examples/bspline_interpolator_example.cpp
new file mode 100644
index 0000000..fdf091f
--- /dev/null
+++ b/cpp/examples/bspline_interpolator_example.cpp
@@ -0,0 +1,360 @@
+/// B-spline interpolator example -- C++ port of examples/b_spline_interpolate_ex.py
+///
+/// Demonstrates BSplineInterpolator with multiple configurations:
+/// cubic with velocity constraints, degree 4 (jerk-continuous),
+/// degree 5 with acceleration constraints, cyclic, and 3D interpolation.
+
+#include
+
+#include "example_utils.hpp"
+
+#include
+#include
+
+namespace ex = interpolatecpp::examples;
+using namespace interpolatecpp::bspline;
+
+/// Helper: print a scalar trajectory table for BSplineInterpolator.
+static void print_scalar_trajectory(const BSplineInterpolator& interp,
+ double t_start, double t_end,
+ int num_samples = 15) {
+ ex::print_trajectory_table(
+ [&](double t) { return interp.evaluate(t)(0); },
+ [&](double t) { return interp.evaluate_derivative(t, 1)(0); },
+ [&](double t) { return interp.evaluate_derivative(t, 2)(0); },
+ t_start, t_end, num_samples);
+}
+
+/// Example 1: Cubic B-spline with velocity constraints (Example 4.16).
+static void example_cubic_bspline() {
+ ex::print_header("Example 1 -- Cubic B-spline with Velocity Constraints");
+
+ Eigen::VectorXd times(7);
+ times << 0, 5, 7, 8, 10, 15, 18;
+
+ Eigen::MatrixXd points(7, 1);
+ points << 3, -2, -5, 0, 6, 12, 8;
+
+ Eigen::VectorXd v0(1);
+ v0 << 2.0;
+ Eigen::VectorXd vn(1);
+ vn << -3.0;
+
+ const BSplineInterpolator interp(
+ 3, // degree
+ points,
+ times,
+ v0, // initial_velocity
+ vn // final_velocity
+ );
+
+ std::cout << "Properties:\n";
+ ex::print_value("Degree", static_cast(interp.degree()), 0);
+ ex::print_value("Num control points", static_cast(interp.n_control_points()), 0);
+ ex::print_value("u_min", interp.u_min());
+ ex::print_value("u_max", interp.u_max());
+
+ std::cout << "\n";
+ ex::print_vector("Times", interp.times());
+ ex::print_vector("Knots", interp.knots());
+
+ std::cout << "\nTrajectory evaluation:\n\n";
+ print_scalar_trajectory(interp, times(0), times(times.size() - 1));
+
+ // Verify boundary velocities
+ ex::print_separator('=');
+ std::cout << "Boundary verification:\n";
+ ex::print_value("Velocity at t_start", interp.evaluate_derivative(times(0), 1)(0));
+ ex::print_value("Expected v0", 2.0);
+ ex::print_value("Velocity at t_end", interp.evaluate_derivative(times(times.size() - 1), 1)(0));
+ ex::print_value("Expected vn", -3.0);
+}
+
+/// Example 2: Degree 4 B-spline (jerk-continuous).
+static void example_jerk_continuous() {
+ ex::print_header("Example 2 -- Degree 4 B-spline (Jerk Continuous)");
+
+ Eigen::VectorXd times(7);
+ times << 0, 5, 7, 8, 10, 15, 18;
+
+ Eigen::MatrixXd points(7, 1);
+ points << 3, -2, -5, 0, 6, 12, 8;
+
+ Eigen::VectorXd v0(1), vn(1), a0(1), an(1);
+ v0 << 2.0;
+ vn << -3.0;
+ a0 << 0.0;
+ an << 0.0;
+
+ const BSplineInterpolator interp(
+ 4, // degree
+ points,
+ times,
+ v0, // initial_velocity
+ vn, // final_velocity
+ a0, // initial_acceleration
+ an // final_acceleration
+ );
+
+ std::cout << "Properties:\n";
+ ex::print_value("Degree", static_cast(interp.degree()), 0);
+ ex::print_value("Num control points", static_cast(interp.n_control_points()), 0);
+
+ std::cout << "\n";
+ ex::print_vector("Knots", interp.knots());
+
+ // Print trajectory with jerk
+ std::cout << "\nTrajectory evaluation (with jerk):\n\n";
+ const double t_start = times(0);
+ const double t_end = times(times.size() - 1);
+ const int num_samples = 15;
+
+ ex::print_full_trajectory_table(
+ [&](double t) -> std::tuple {
+ return {interp.evaluate(t)(0),
+ interp.evaluate_derivative(t, 1)(0),
+ interp.evaluate_derivative(t, 2)(0),
+ interp.evaluate_derivative(t, 3)(0)};
+ },
+ t_start, t_end, num_samples);
+
+ // Verify boundary conditions
+ ex::print_separator('=');
+ std::cout << "Boundary verification:\n";
+ ex::print_value("Velocity at t_start", interp.evaluate_derivative(t_start, 1)(0));
+ ex::print_value("Acceleration at t_start", interp.evaluate_derivative(t_start, 2)(0));
+ ex::print_value("Velocity at t_end", interp.evaluate_derivative(t_end, 1)(0));
+ ex::print_value("Acceleration at t_end", interp.evaluate_derivative(t_end, 2)(0));
+}
+
+/// Example 3: Degree 5 B-spline with acceleration constraints.
+static void example_degree5() {
+ ex::print_header("Example 3 -- Degree 5 B-spline with Acceleration Constraints");
+
+ Eigen::VectorXd times(7);
+ times << 0, 1, 2, 3, 4, 5, 6;
+
+ Eigen::MatrixXd points(7, 1);
+ points << 0, 2, 1, 3, 2, 4, 3;
+
+ Eigen::VectorXd v0(1), vn(1), a0(1), an(1);
+ v0 << 1.0;
+ vn << -1.0;
+ a0 << 0.0;
+ an << 0.0;
+
+ const BSplineInterpolator interp(
+ 5, // degree
+ points,
+ times,
+ v0, vn,
+ a0, an
+ );
+
+ std::cout << "Properties:\n";
+ ex::print_value("Degree", static_cast(interp.degree()), 0);
+ ex::print_value("Num control points", static_cast(interp.n_control_points()), 0);
+
+ std::cout << "\n";
+ ex::print_vector("Knots", interp.knots());
+
+ std::cout << "\nTrajectory evaluation:\n\n";
+ print_scalar_trajectory(interp, times(0), times(times.size() - 1));
+
+ // Verify boundary conditions
+ ex::print_separator('=');
+ std::cout << "Boundary verification:\n";
+ const double t_start = times(0);
+ const double t_end = times(times.size() - 1);
+ ex::print_value("Velocity at t_start", interp.evaluate_derivative(t_start, 1)(0));
+ ex::print_value("Acceleration at t_start", interp.evaluate_derivative(t_start, 2)(0));
+ ex::print_value("Velocity at t_end", interp.evaluate_derivative(t_end, 1)(0));
+ ex::print_value("Acceleration at t_end", interp.evaluate_derivative(t_end, 2)(0));
+}
+
+/// Example 4: Cyclic B-spline (Example 4.17).
+static void example_cyclic() {
+ ex::print_header("Example 4 -- Cyclic B-spline (Degree 4)");
+
+ Eigen::VectorXd times(7);
+ times << 0, 5, 7, 8, 10, 15, 18;
+
+ // Last point equals first point for cyclic
+ Eigen::MatrixXd points(7, 1);
+ points << 3, -2, -5, 0, 6, 12, 3;
+
+ const BSplineInterpolator interp(
+ 4, // degree
+ points,
+ times,
+ std::nullopt, // initial_velocity
+ std::nullopt, // final_velocity
+ std::nullopt, // initial_acceleration
+ std::nullopt, // final_acceleration
+ true // cyclic
+ );
+
+ std::cout << "Properties:\n";
+ ex::print_value("Degree", static_cast(interp.degree()), 0);
+ ex::print_value("Num control points", static_cast(interp.n_control_points()), 0);
+ ex::print_value("Cyclic", 1.0, 0);
+
+ std::cout << "\n";
+ ex::print_vector("Knots", interp.knots());
+
+ // Trajectory with jerk and snap
+ std::cout << "\nTrajectory evaluation:\n\n";
+ const double t_start = times(0);
+ const double t_end = times(times.size() - 1);
+
+ const int w = 14;
+ const int p = 6;
+ std::cout << std::right
+ << std::setw(w) << "Time"
+ << std::setw(w) << "Position"
+ << std::setw(w) << "Velocity"
+ << std::setw(w) << "Acceleration"
+ << std::setw(w) << "Jerk"
+ << std::setw(w) << "Snap" << "\n";
+ ex::print_separator('-', 6 * w);
+
+ const int num_samples = 15;
+ for (int i = 0; i <= num_samples; ++i) {
+ const double t = t_start + (t_end - t_start) * static_cast(i) / num_samples;
+ std::cout << std::fixed << std::setprecision(p)
+ << std::setw(w) << t
+ << std::setw(w) << interp.evaluate(t)(0)
+ << std::setw(w) << interp.evaluate_derivative(t, 1)(0)
+ << std::setw(w) << interp.evaluate_derivative(t, 2)(0)
+ << std::setw(w) << interp.evaluate_derivative(t, 3)(0)
+ << std::setw(w) << interp.evaluate_derivative(t, 4)(0)
+ << "\n";
+ }
+ std::cout << "\n";
+
+ // Verify cyclic continuity: values at start and end should match
+ ex::print_separator('=');
+ std::cout << "Cyclic continuity verification:\n";
+ ex::print_value("Position at t_start", interp.evaluate(t_start)(0));
+ ex::print_value("Position at t_end", interp.evaluate(t_end)(0));
+ ex::print_value("Velocity at t_start", interp.evaluate_derivative(t_start, 1)(0));
+ ex::print_value("Velocity at t_end", interp.evaluate_derivative(t_end, 1)(0));
+ ex::print_value("Acceleration at t_start", interp.evaluate_derivative(t_start, 2)(0));
+ ex::print_value("Acceleration at t_end", interp.evaluate_derivative(t_end, 2)(0));
+}
+
+/// Example 5: 3D B-spline interpolation.
+static void example_3d() {
+ ex::print_header("Example 5 -- 3D B-spline Interpolation");
+
+ Eigen::VectorXd times(5);
+ times << 0, 1, 2, 3, 4;
+
+ Eigen::MatrixXd points(5, 3);
+ points << 0, 0, 0,
+ 1, 1, 2,
+ 2, 0, 3,
+ 3, -1, 2,
+ 4, 0, 0;
+
+ const int degree = 3;
+ const BSplineInterpolator interp(degree, points, times);
+
+ std::cout << "Properties:\n";
+ ex::print_value("Degree", static_cast(interp.degree()), 0);
+ ex::print_value("Num control points", static_cast(interp.n_control_points()), 0);
+ ex::print_value("Dimension", static_cast(interp.dimension()), 0);
+ std::cout << "Continuity: C^" << (degree - 1) << " (continuous acceleration)\n";
+
+ std::cout << "\n";
+ ex::print_vector("Knots", interp.knots());
+ ex::print_matrix("Control points (3D)", interp.control_points());
+
+ // Print original and interpolated points
+ std::cout << "\nOriginal points to interpolate:\n";
+ for (Eigen::Index i = 0; i < points.rows(); ++i) {
+ std::cout << " Point " << i << " (t=" << times(i) << "): ("
+ << std::fixed << std::setprecision(1)
+ << points(i, 0) << ", " << points(i, 1) << ", " << points(i, 2) << ")\n";
+ }
+
+ // Evaluate at intermediate times
+ std::cout << "\nInterpolated points at intermediate times:\n";
+ const std::vector t_samples = {0.5, 1.5, 2.5, 3.5};
+ for (const double t : t_samples) {
+ const Eigen::VectorXd pt = interp.evaluate(t);
+ std::cout << " t = " << std::fixed << std::setprecision(1) << t << ": ("
+ << std::setprecision(3) << pt(0) << ", " << pt(1) << ", " << pt(2) << ")\n";
+ }
+
+ // Print full 3D trajectory
+ std::cout << "\n3D trajectory table:\n\n";
+ const int w = 14;
+ std::cout << std::right
+ << std::setw(w) << "Time"
+ << std::setw(w) << "X" << std::setw(w) << "Y" << std::setw(w) << "Z" << "\n";
+ ex::print_separator('-', 4 * w);
+
+ const int num_samples = 12;
+ for (int i = 0; i <= num_samples; ++i) {
+ const double t = times(0) + (times(times.size() - 1) - times(0))
+ * static_cast(i) / num_samples;
+ const Eigen::VectorXd pt = interp.evaluate(t);
+ std::cout << std::fixed << std::setprecision(4)
+ << std::setw(w) << t
+ << std::setw(w) << pt(0)
+ << std::setw(w) << pt(1)
+ << std::setw(w) << pt(2) << "\n";
+ }
+
+ // Try higher degrees with more points
+ ex::print_separator('=');
+ std::cout << "Degree sensitivity with 5 points:\n\n";
+
+ // Degree 3 works (already demonstrated). Try degree 4 and 5.
+ for (const int deg : {4, 5}) {
+ std::cout << " Trying degree " << deg << " with " << points.rows() << " points... ";
+ try {
+ const BSplineInterpolator test_interp(deg, points, times);
+ std::cout << "Succeeded! (n_cp=" << test_interp.n_control_points() << ")\n";
+ } catch (const std::exception& e) {
+ std::cout << "Failed: " << e.what() << "\n";
+ }
+ }
+
+ // Extended points for higher degrees
+ ex::print_separator();
+ std::cout << "Extended points (7 points) for higher degrees:\n\n";
+
+ Eigen::VectorXd times_ext(7);
+ times_ext << 0, 1, 2, 3, 4, 5, 6;
+
+ Eigen::MatrixXd points_ext(7, 3);
+ points_ext << 0, 0, 0,
+ 1, 1, 2,
+ 2, 0, 3,
+ 3, -1, 2,
+ 4, 0, 0,
+ 5, 1, -1,
+ 6, 0, -2;
+
+ for (const int deg : {3, 4, 5}) {
+ std::cout << " Degree " << deg << " with 7 points: ";
+ try {
+ const BSplineInterpolator ext_interp(deg, points_ext, times_ext);
+ std::cout << "Succeeded! (n_cp=" << ext_interp.n_control_points() << ")\n";
+ } catch (const std::exception& e) {
+ std::cout << "Failed: " << e.what() << "\n";
+ }
+ }
+}
+
+int main() {
+ example_cubic_bspline();
+ example_jerk_continuous();
+ example_degree5();
+ example_cyclic();
+ example_3d();
+
+ return 0;
+}
diff --git a/cpp/examples/concepts_example.cpp b/cpp/examples/concepts_example.cpp
new file mode 100644
index 0000000..abde367
--- /dev/null
+++ b/cpp/examples/concepts_example.cpp
@@ -0,0 +1,398 @@
+/// C++20 concepts example -- C++ port of examples/protocols_ex.py
+///
+/// Demonstrates how C++20 concepts serve as the compile-time equivalent of
+/// Python protocols (PEP 544). Template functions constrained by each concept
+/// accept any concrete type that satisfies the required interface, with no
+/// inheritance needed.
+
+#include
+
+// Concrete types for ScalarTrajectory
+#include
+#include
+
+// Concrete types for CurveEvaluator
+#include