From 8d22262afd0f91b49cb046d090b186c4305d470b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Nusser?= <andre.nusser@googlemail.com>
Date: Sat, 25 Apr 2020 22:45:47 +0200
Subject: Introduce Steffen monotone cubic splines.

---
 src/powermap.cc | 65 ++++++++++++++++++++++++++++++++++++++++++++-------------
 src/powermap.h  |  6 +++++-
 2 files changed, 56 insertions(+), 15 deletions(-)

diff --git a/src/powermap.cc b/src/powermap.cc
index 2dd5df4..ca17a60 100644
--- a/src/powermap.cc
+++ b/src/powermap.cc
@@ -156,18 +156,43 @@ void Powermap::updateSpline()
 	assert(0. <= fixed[0].out && fixed[0].out <= fixed[1].out &&
 	       fixed[1].out <= fixed[2].out && fixed[2].out <= 1.);
 
-	// TODO: What to do if fixed[0] is (0,0) or fixed[2] is (1,1)??
 	Powers X = shelf ? Powers{fixed[0].in, fixed[1].in, fixed[2].in}
 	                 : Powers{0., fixed[0].in, fixed[1].in, fixed[2].in, 1.};
-	Powers P = shelf ? Powers{fixed[0].out, fixed[1].out, fixed[2].out}
+	Powers Y = shelf ? Powers{fixed[0].out, fixed[1].out, fixed[2].out}
 	                 : Powers{0., fixed[0].out, fixed[1].out, fixed[2].out, 1.};
 
+	auto slopes = calcSlopes(X, Y);
+
+	if (shelf) {
+		assert(slopes.size() == 3);
+		this->m[1] = slopes[0];
+		this->m[2] = slopes[1];
+		this->m[3] = slopes[2];
+	}
+	else {
+		assert(slopes.size() == 5);
+		for (std::size_t i = 0; i < m.size(); ++i) {
+			this->m[i] = slopes[i];
+		}
+	}
+
+	spline_needs_update = false;
+}
+
+std::vector<float> Powermap::calcSlopes(Powers const& X, Powers const& Y)
+{
+	// return calcSlopesFritschCarlson(X, Y);
+	return calcSlopesSteffen(X, Y);
+}
+
+std::vector<float> Powermap::calcSlopesFritschCarlson(Powers const& X, Powers const& Y)
+{
 	Powers deltas(X.size()-1);
 	Powers m(X.size());
 
 	// 1
 	for (std::size_t i = 0; i < deltas.size(); ++i) {
-		deltas[i] = (P[i+1]-P[i])/(X[i+1]-X[i]);
+		deltas[i] = (Y[i+1]-Y[i])/(X[i+1]-X[i]);
 	}
 
 	// 2
@@ -227,20 +252,32 @@ void Powermap::updateSpline()
 		}
 	}
 
-	if (shelf) {
-		assert(m.size() == 3);
-		this->m[1] = m[0];
-		this->m[2] = m[1];
-		this->m[3] = m[2];
+	return m;
+}
+
+std::vector<float> Powermap::calcSlopesSteffen(Powers const& X, Powers const& Y)
+{
+	Powers m(X.size());
+
+	Powers d(X.size()-1);
+	Powers h(X.size()-1);
+	for (std::size_t i = 0; i < d.size(); ++i) {
+		h[i] = X[i+1]-X[i];
+		d[i] = (Y[i+1]-Y[i])/h[i];
 	}
-	else {
-		assert(m.size() == 5);
-		for (std::size_t i = 0; i < m.size(); ++i) {
-			this->m[i] = m[i];
-		}
+
+	m.front() = d.front();
+	for (std::size_t i = 1; i < m.size()-1; ++i) {
+		m[i] = (d[i-1] + d[i])/2.;
 	}
+	m.back() = d.back();
 
-	spline_needs_update = false;
+	for (std::size_t i = 1; i < m.size()-1; ++i) {
+		auto const min_d = 2*std::min(d[i-1], d[i]);
+		m[i] = std::min<float>(min_d, (h[i]*d[i-1] + h[i-1]*d[i])/(h[i-1]+h[i]));
+	}
+
+	return m;
 }
 
 Power Powermap::clamp(Power in, Power min, Power max) const
diff --git a/src/powermap.h b/src/powermap.h
index 29a3c81..b942722 100644
--- a/src/powermap.h
+++ b/src/powermap.h
@@ -64,9 +64,13 @@ private:
 
 	// spline parameters (deterministically computed from the input parameters)
 	std::array<float, 5> m;
-
 	bool spline_needs_update;
+
 	void updateSpline();
+	std::vector<float> calcSlopes(Powers const& X, Powers const& P);
+	std::vector<float> calcSlopesFritschCarlson(Powers const& X, Powers const& P);
+	std::vector<float> calcSlopesSteffen(Powers const& X, Powers const& P);
+
 	Power clamp(Power in, Power min, Power max) const;
 
 	const Power eps = 1e-3;
-- 
cgit v1.2.3