From 4c616cf4a001f3b4adda87889a98e403c2e8ab8a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Andr=C3=A9=20Nusser?= <andre.nusser@googlemail.com>
Date: Sat, 25 Apr 2020 23:03:57 +0200
Subject: Just leave the Steffen method in and clean up powermap.

---
 src/powermap.cc | 78 ++-------------------------------------------------------
 src/powermap.h  |  7 ++----
 2 files changed, 4 insertions(+), 81 deletions(-)

diff --git a/src/powermap.cc b/src/powermap.cc
index ca17a60..914b5c1 100644
--- a/src/powermap.cc
+++ b/src/powermap.cc
@@ -179,83 +179,9 @@ void Powermap::updateSpline()
 	spline_needs_update = false;
 }
 
+// This follows the monotone cubic spline algorithm of Steffen, from:
+// "A Simple Method for Monotonic Interpolation in One Dimension"
 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] = (Y[i+1]-Y[i])/(X[i+1]-X[i]);
-	}
-
-	// 2
-	m[0] = deltas[0];
-	for (std::size_t i = 1; i < deltas.size(); ++i) {
-		m[i] =  (deltas[i-1] + deltas[i])/2.;
-	}
-	m.back() = deltas.back();
-
-	// 3
-	std::vector<bool> ignore(deltas.size(), false);
-	for (std::size_t i = 0; i < deltas.size(); ++i) {
-		if (deltas[i] == 0) {
-			m[i] = 0;
-			m[i+1] = 0;
-			ignore[i] = true;
-		}
-	}
-
-	for (std::size_t i = 0; i < deltas.size(); ++i) {
-		if (ignore[i]) {
-			continue;
-		}
-
-		// 4
-		auto alpha = m[i]/deltas[i];
-		auto beta = m[i+1]/deltas[i];
-		assert(alpha >= 0.);
-		assert(beta >= 0.);
-
-		// 5
-		// TODO: expose this parameter for testing both
-		bool const option1 = false;
-		if (option1) {
-			if (alpha > 3 || beta > 3) {
-				m[i] = 3*deltas[i];
-			}
-		}
-		else {
-			auto const a2b2 = alpha*alpha + beta*beta;
-
-			// hard change to enforce monotonicity
-			// if (a2b2 > 9) {
-			//     auto const tau = 3./sqrt(a2b2);
-			//     m[i] = tau*alpha*deltas[i];
-			//     m[i+1] = tau*alpha*deltas[i];
-			// }
-
-			// soft change to enforce monotonicity
-			Power threshold = 4.5; // must be >= 0 and < 9.
-			if (a2b2 >= threshold) {
-				auto const l = std::min(1., (a2b2-threshold)/(9.-threshold));
-				auto const tau = 3./sqrt(a2b2);
-				m[i] = (1-l)*m[i] + l*tau*alpha*deltas[i];
-				m[i+1] = (1-l)*m[i+1] + l*tau*alpha*deltas[i];
-			}
-		}
-	}
-
-	return m;
-}
-
-std::vector<float> Powermap::calcSlopesSteffen(Powers const& X, Powers const& Y)
 {
 	Powers m(X.size());
 
diff --git a/src/powermap.h b/src/powermap.h
index b942722..2cb7efb 100644
--- a/src/powermap.h
+++ b/src/powermap.h
@@ -63,15 +63,12 @@ private:
 	bool shelf;
 
 	// spline parameters (deterministically computed from the input parameters)
-	std::array<float, 5> m;
 	bool spline_needs_update;
+	std::array<float, 5> m;
+	const Power eps = 1e-3;
 
 	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