cmake, more tests, unique header and thread-safe

CMakeLists.txt

@@ -0,0 +1,26 @@
+cmake_minimum_required(VERSION 3.8.0)
+
+set(CMAKE_CXX_STANDARD 11)
+
+project(cxx-spline)
+
+add_subdirectory(test)
+enable_testing ()
+
+add_subdirectory(demo)
+
+add_library(cxx-spline INTERFACE)
+
+target_include_directories(cxx-spline INTERFACE
+  $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+  $<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>)
+
+set_target_properties(cxx-spline PROPERTIES
+  PUBLIC_HEADER cxx-spline.hpp)
+
+install(TARGETS cxx-spline
+        EXPORT ${PROJECT_NAME}_Targets
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+	PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR})

README.md

@@ -1,27 +1,48 @@
 # Cubic splines for C++
 
 ![C++11][cxx-badge]
-![Test status][test-badge]
 [![Boost License][license-badge]](LICENSE.txt)
 
+![alt tag](teaser.png)
+
 Header-only cubic spline interpolator for C++11 and above.
 
 [cxx-badge]: https://img.shields.io/badge/C%2B%2B-11-orange.svg
-[test-badge]: https://github.com/snsinfu/cxx-spline/workflows/test/badge.svg
 [license-badge]: https://img.shields.io/badge/license-Boost-blue.svg
 
 - [Usage](#usage)
-- [Testing](#testing)
 - [License](#license)
 
+## Prerequisites
+
+Catch (standard Linux package from the package manager):
+
+```
+sudo apt install catch
+```
+
 ## Usage
 
-Copy [spline.hpp](include/spline.hpp) into your include directory. Use
-`cubic_spline` class like this:
+It's just the header `cxx-spline.h` which you need. Either
+use `make install` to install it on your system or just copy
+the header into your project directory.
+
+### Installing the header on your system
+```
+cmake .
+make
+make test
+sudo make install
+```
+or
+
+### Copy [cxx-spline.h](cxx-spline.h) into your include directory.
+
+### Use `cubic_spline` class like this
 
 ```c++
 #include <vector>
-#include <spline.hpp>
+#include <cxx-spline.h>
 
 int main()
 {
@@ -44,7 +65,7 @@ coordinate values. Example:
 ```c++
 #include <iostream>
 #include <vector>
-#include <spline.hpp>
+#include <cxx-spline.h>
 
 int main()
 {
@@ -64,6 +85,7 @@ int main()
     }
 }
 ```
+See again the demo directory for this example in action.
 
 ## Boundary conditions
 
@@ -84,9 +106,9 @@ used for visualization.
 ## Testing
 
 ```sh
-git clone https://github.com/snsinfu/cxx-spline.git
-cd cxx-spline/test
+cmake .
 make
+make test
 ```
 
 ## License

include/spline.hpp → cxx-spline.h

@@ -1,6 +1,7 @@
 // Cubic spline
 
 // Copyright snsinfu 2020-2021.
+// Copyright Bernd Porr 2022.
 // Distributed under the Boost Software License, Version 1.0.
 //
 // Permission is hereby granted, free of charge, to any person or organization
@@ -34,7 +35,6 @@
 #include <stdexcept>
 #include <vector>
 
-
 namespace detail_cubic_spline
 {
     /*
@@ -106,7 +106,7 @@ namespace detail_cubic_spline
 
 
 /*
- * Cubic spline functor for interpolating one-dimensional series of values.
+ * Cubic spline class for interpolating one-dimensional series of values.
  */
 class cubic_spline
 {
@@ -135,43 +135,53 @@ class cubic_spline
     static constexpr auto natural = boundary_conditions::natural;
     static constexpr auto not_a_knot = boundary_conditions::not_a_knot;
 
-    /*
-     * Constructs a cubic spline function that passes through given knots.
-     */
-    cubic_spline(
-        std::vector<double> const& t,
-        std::vector<double> const& x,
-        boundary_conditions bc = natural
-    )
-    {
-        make_spline(t, x, bc);
-        make_bins();
-    }
-
-    /*
+    /**
      * Evaluates the cubic splines at `t`.
      */
     double operator()(double t) const
     {
+        if (_splines.empty()) throw std::invalid_argument("No spline data");
         spline_data const& spline = find_spline(t);
-
         double value = spline.coefficients[order];
         for (int i = order - 1; i >= 0; i--) {
             value *= t - spline.knot;
             value += spline.coefficients[i];
         }
-
         return value;
     }
 
-private:
+    /**
+     * Checks if we have spline data.
+     */
+    bool hasSpline() const {
+        return !(_splines.empty());
+    }
+
+    /**
+     * Gets the lowest timestamp or x value
+     */
+    double lower_bound() const {
+        if (_splines.empty()) throw std::invalid_argument("No spline data");
+        return _lower_bound;
+    }
+
+    /**
+     * Gets the highest timestamp or x value
+     */
+    double upper_bound() const {
+        if (_splines.empty()) throw std::invalid_argument("No spline data");
+        return _upper_bound;
+    }
+
+    cubic_spline () = default;
+
     /*
-     * Constructs spline segments `_splines` for given set of knot points.
+     * Constructs a cubic spline function that passes through given knots.
      */
-    void make_spline(
+    cubic_spline (
         std::vector<double> const& knots,
         std::vector<double> const& values,
-        boundary_conditions bc
+        boundary_conditions bc = natural
     )
     {
         if (knots.size() != values.size()) {
@@ -268,6 +278,8 @@ class cubic_spline
         detail_cubic_spline::solve_tridiagonal_system(L, D, U, Y);
         auto const& M = Y;
 
+	_bins.clear();
+
         // Derive the polynomial coefficients of each spline segment from the
         // second derivatives `M`.
         _splines.clear();
@@ -287,14 +299,6 @@ class cubic_spline
 
         _lower_bound = knots.front();
         _upper_bound = knots.back();
-    }
-
-    /*
-     * Builds a bin-based index `_bins` that is used to quickly find a spline
-     * segment covering a query point.
-     */
-    void make_bins()
-    {
         _bin_interval = (_upper_bound - _lower_bound) / double(_splines.size());
 
         // We need to map uniformly-spaced bins to spline segments that may
@@ -321,6 +325,8 @@ class cubic_spline
         _bins.shrink_to_fit();
     }
 
+private:
+
     /*
      * Returns the spline segment that covers given point `t`.
      */
@@ -339,8 +345,6 @@ class cubic_spline
         }
         std::size_t index = _bins[bin];
 
-        assert(t >= _splines[index].knot);
-
         for (; index + 1 < _splines.size(); index++) {
             if (t < _splines[index + 1].knot) {
                 break;
@@ -350,12 +354,11 @@ class cubic_spline
         return _splines[index];
     }
 
-private:
     std::vector<spline_data> _splines;
     std::vector<std::size_t> _bins;
-    double _lower_bound;
-    double _upper_bound;
-    double _bin_interval;
+    double _lower_bound = 0;
+    double _upper_bound = 0;
+    double _bin_interval = 0;
 };
 
 #endif

demo/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.8.0)
+
+set(CMAKE_CXX_STANDARD 11)
+
+include_directories(..)
+
+add_executable(series1D series1D.cpp)
+add_executable(series2D series2D.cpp)

demo/README.md

@@ -0,0 +1,14 @@
+# Spline demos
+
+## 1D
+
+Run `series1D` which outputs `1Dspline.dat` of the interpolated data.
+It also shows that it's possible to do prediction into the future
+by adding values of t which are greater than the upper bound of t. 
+
+Plot the data with `plot1Dspline.py`.
+
+## 2D
+
+Run `series2D` which approximates two independent coordinates.
+

demo/plot1Dspline.py

@@ -0,0 +1,18 @@
+#!/usr/bin/python3
+
+import matplotlib.pyplot as plt
+import scipy
+import numpy as np
+
+t = np.array([0, 1, 2, 3, 4, 5 ]);
+y = np.array([1, 2, 3, 2, 1, 2 ]);
+
+
+sp = np.loadtxt("1Dspline.dat")
+
+tp = sp[:,0]
+yp = sp[:,1]
+
+plt.plot(t,y,"x")
+plt.plot(tp,yp)
+plt.show()

demo/series1D.cpp

@@ -0,0 +1,25 @@
+#include <vector>
+#include <cxx-spline.h>
+#include <cstdio>
+
+int main()
+{
+    // Some points (t,y) on a curve y = f(t)
+    const std::vector<double> t = { 0, 1, 2, 3, 4, 5 };
+    const std::vector<double> y = { 1, 2, 3, 2, 1, 2 };
+
+    // Spline interpolation (and extrapolation) of the points
+    cubic_spline spline(t, y);
+
+    FILE* f=fopen("1Dspline.dat","wt");
+    for(double tp=-0.5;tp<6;tp+=0.1) {
+	    const double yp = spline(tp);
+	    printf("f(%f)=%f\n",tp,yp);
+	    fprintf(f,"%f %f\n",tp,yp);
+    }
+    fclose(f);
+
+    printf("Lower bound = %f, Upper bound = %f\n",
+	   spline.lower_bound(), spline.upper_bound()
+	    );
+}

demo/series2D.cpp

@@ -0,0 +1,22 @@
+#include <vector>
+#include <cxx-spline.h>
+#include <iostream>
+
+
+int main()
+{
+    // Samples from a circle.
+    std::vector<double> t = { 0, 0.25, 0.5, 0.75, 1 };
+    std::vector<double> x = { 1, 0, -1, 0, 1 };
+    std::vector<double> y = { 0, 1, 0, -1, 0 };
+
+    // Spline interpolation of each coordinate.
+    cubic_spline spline_x(t, x);
+    cubic_spline spline_y(t, y);
+
+    for (int i = 0; i <= 100; i++) {
+        double sx = spline_x(i / 100.0);
+        double sy = spline_y(i / 100.0);
+        std::cout << sx << '\t' << sy << '\n';
+    }
+}