cmake, more tests, unique header and thread-safe

CMakeLists.txt

@@ -0,0 +1,26 @@
+cmake_minimum_required(VERSION 3.1.0)
+
+set(CMAKE_CXX_STANDARD 14)
+
+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,50 @@
 # 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.hpp` which you need. Either
+use `make install` to install it on your system or just copy
+the header into your project directory.
+
+The class is thread-safe.
+
+### Installing the header on your system
+```
+cmake .
+make
+make test
+sudo make install
+```
+or
+
+### Copy [cxx-spline.hpp](cxx-spline.hpp) into your include directory.
+
+### Use `cubic_spline` class like this
 
 ```c++
 #include <vector>
-#include <spline.hpp>
+#include <cxx-spline.hpp>
 
 int main()
 {
@@ -30,21 +53,24 @@ int main()
     std::vector<double> y = { 1, 2, 3, 2, 1, 2 };
 
     // Spline interpolation (and extrapolation) of the points
-    cubic_spline spline(t, y);
+    cubic_spline spline;
+    spline.calc(t, y);
 
     spline(0.5); // 1.44976
     spline(1.5); // 2.65072
     spline(6.0); // 3
 }
 ```
+You can call `calc` as much as you like to update the splines. See the
+demo directory for this file in action.
 
 To interpolate a 2D (or higher dimensional) curve, just create splines for each
 coordinate values. Example:
 
 ```c++
 #include <iostream>
 #include <vector>
-#include <spline.hpp>
+#include <cxx-spline.hpp>
 
 int main()
 {
@@ -54,8 +80,10 @@ int main()
     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);
+    cubic_spline spline_x;
+    spline_x.calc(t, x);
+    cubic_spline spline_y;
+    spline_y.calc(t, y);
 
     for (int i = 0; i <= 100; i++) {
         double sx = spline_x(i / 100.0);
@@ -64,6 +92,7 @@ int main()
     }
 }
 ```
+See again the demo directory for this example in action.
 
 ## Boundary conditions
 
@@ -72,8 +101,10 @@ Pass `cubic_spline::natural` or `cubic_spline::not_a_knot` to the constructor
 to choose the boundary conditions. Natural is the default.
 
 ```c++
-cubic_spline natural_spline(t, x, cubic_spline::natural);
-cubic_spline notaknot_spline(t, x, cubic_spline::not_a_knot);
+cubic_spline natural_spline;
+natural_spline.calc(t, x, cubic_spline::natural);
+cubic_spline notaknot_spline;
+notaknot_spline.calc(t, x, cubic_spline::not_a_knot);
 ```
 
 The *natural* (or free-end) spline gives the minimum-energy curve in a certain
@@ -84,9 +115,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.hpp

@@ -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
@@ -33,7 +34,7 @@
 #include <cstdlib>
 #include <stdexcept>
 #include <vector>
-
+#include <mutex>
 
 namespace detail_cubic_spline
 {
@@ -138,32 +139,52 @@ class cubic_spline
     /*
      * Constructs a cubic spline function that passes through given knots.
      */
-    cubic_spline(
+    void calc(
         std::vector<double> const& t,
         std::vector<double> const& x,
         boundary_conditions bc = natural
     )
     {
-        make_spline(t, x, bc);
-        make_bins();
+	    make_spline(t, x, bc);
+        _splinemutex.unlock();
     }
 
     /*
      * Evaluates the cubic splines at `t`.
      */
-    double operator()(double t) const
+    double operator()(double t)
     {
+        _splinemutex.lock();
         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];
         }
-
+        _splinemutex.unlock();
         return value;
     }
 
+    /**
+     * Gets the lowest timestamp or x value
+     */
+    double getLowerBound() {
+        _splinemutex.lock();
+	const double lower = _lower_bound;
+        _splinemutex.unlock();
+        return lower;
+    }
+
+    /**
+     * Gets the highest timestamp or x value
+     */
+    double getUpperBound() {
+        _splinemutex.lock();
+	const double upper = _upper_bound;
+        _splinemutex.unlock();
+        return upper;
+    }
+
 private:
     /*
      * Constructs spline segments `_splines` for given set of knot points.
@@ -268,6 +289,9 @@ class cubic_spline
         detail_cubic_spline::solve_tridiagonal_system(L, D, U, Y);
         auto const& M = Y;
 
+	_splinemutex.lock();
+	_bins.clear();
+
         // Derive the polynomial coefficients of each spline segment from the
         // second derivatives `M`.
         _splines.clear();
@@ -287,14 +311,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
@@ -319,6 +335,7 @@ class cubic_spline
         }
 
         _bins.shrink_to_fit();
+	_splinemutex.unlock();
     }
 
     /*
@@ -339,8 +356,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;
@@ -353,9 +368,10 @@ class cubic_spline
 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;
+    std::mutex _splinemutex;
 };
 
 #endif

demo/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.1.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,26 @@
+#include <vector>
+#include <cxx-spline.hpp>
+#include <stdio.h>
+
+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;
+    spline.calc(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.getLowerBound(), spline.getUpperBound()
+	    );
+}

demo/series2D.cpp

@@ -0,0 +1,27 @@
+#include <vector>
+#include <cxx-spline.hpp>
+#include <iostream>
+
+#include <iostream>
+#include <vector>
+#include <cxx-spline.hpp>
+
+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;
+    spline_x.calc(t, x);
+    cubic_spline spline_y;
+    spline_y.calc(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';
+    }
+}