various tweaks to catmull-rom

rtengine/diagonalcurves.cc

@@ -359,12 +359,14 @@ void catmull_rom_chain(int n_points, int n_cp, double *x, double *y,
                        std::vector<double> &res_x, std::vector<double> &res_y)
 {
     static const double epsilon = 1e-5;
-    double xr = x[1] - x[0];
+    // double xr = x[1] - x[0];
-    double yr = y[1] - y[0];
+    // double yr = y[1] - y[0];
+    double xr = x[n_cp-1] - x[0];
+    double yr = y[n_cp-1] - y[0];
     double x_first = x[0] - xr * 0.1;
     double y_first = xr > epsilon ? (yr / xr) * (x_first - x[0]) + y[0] : y[0];
-    xr = x[n_cp-1] - x[n_cp-2];
+    // xr = x[n_cp-1] - x[n_cp-2];
-    yr = y[n_cp-1] - x[n_cp-2];
+    // yr = y[n_cp-1] - x[n_cp-2];
     double x_last = x[n_cp-1] + xr * 0.1;
     double y_last = xr > epsilon ? (yr / xr) * (x_last - x[0]) + y[0] : y[0];
 
@@ -477,7 +479,7 @@ double DiagonalCurve::getVal (double t) const
         if (it+1 < poly_x.end() && t - *it > *(it+1) - t) {
             ++d;
         }
-        return *(poly_y.begin() + d);
+        return LIM01(*(poly_y.begin() + d));
         break;
     }
 

rtengine/histmatching.cc

@@ -97,6 +97,55 @@ int findMatch(int val, const std::vector<int> &cdf, int j)
 }
 
 
+class CubicSplineCurve: public DiagonalCurve {
+public:
+    CubicSplineCurve(const std::vector<double> &points):
+        DiagonalCurve({DCT_Linear})
+    {
+        N = points.size() / 2;
+        x = new double[N];
+        y = new double[N];
+
+        for (int i = 0; i < N; ++i) {
+            x[i] = points[2*i];
+            y[i] = points[2*i+1];
+        }
+        kind = DCT_Spline;
+        spline_cubic_set();
+    }
+
+    double getVal(double t) const
+    {
+        // values under and over the first and last point
+        if (t > x[N - 1]) {
+            return y[N - 1];
+        } else if (t < x[0]) {
+            return y[0];
+        }
+
+        // do a binary search for the right interval:
+        unsigned int k_lo = 0, k_hi = N - 1;
+
+        while (k_hi > 1 + k_lo) {
+            unsigned int k = (k_hi + k_lo) / 2;
+
+            if (x[k] > t) {
+                k_hi = k;
+            } else {
+                k_lo = k;
+            }
+        }
+
+        double h = x[k_hi] - x[k_lo];
+
+        double a = (x[k_hi] - t) / h;
+        double b = (t - x[k_lo]) / h;
+        double r = a * y[k_lo] + b * y[k_hi] + ((a * a * a - a) * ypp[k_lo] + (b * b * b - b) * ypp[k_hi]) * (h * h) * 0.1666666666666666666666666666666;
+        return LIM01(r);
+    }
+};
+
+
 void mappingToCurve(const std::vector<int> &mapping, std::vector<double> &curve)
 {
     curve.clear();
@@ -117,8 +166,9 @@ void mappingToCurve(const std::vector<int> &mapping, std::vector<double> &curve)
 
     auto coord = [](int v) -> double { return double(v)/255.0; };
     auto doit =
-        [&](int start, int stop, int step, bool addstart) -> void
+        [&](int start, int stop, int step, bool addstart, int maxdelta=0) -> void
         {
+            if (!maxdelta) maxdelta = step * 2;
             int prev = start;
             if (addstart && mapping[start] >= 0) {
                 curve.push_back(coord(start));
@@ -131,7 +181,7 @@ void mappingToCurve(const std::vector<int> &mapping, std::vector<double> &curve)
                 }
                 bool change = i > 0 && v != mapping[i-1];
                 int diff = i - prev;
-                if ((change && std::abs(diff - step) <= 1) || diff > step * 2) {
+                if ((change && std::abs(diff - step) <= 1) || diff > maxdelta) {
                     curve.push_back(coord(i));
                     curve.push_back(coord(v));
                     prev = i;
@@ -152,13 +202,14 @@ void mappingToCurve(const std::vector<int> &mapping, std::vector<double> &curve)
     int end = mapping.size();
     if (idx <= end / 3) {
         doit(start, idx, idx / 2, true);
-        doit(idx, end, (end - idx) / 3, false);
+        step = (end - idx) / 4;
+        doit(idx, end, step, false, step);
     } else {
         doit(start, idx, idx > step ? step : idx / 2, true);
-        doit(idx, int(mapping.size()), step, idx - step > step / 2 && std::abs(curve[curve.size()-2] - coord(idx)) > 0.01);
+        doit(idx, end, step, idx - step > step / 2 && std::abs(curve[curve.size()-2] - coord(idx)) > 0.01);
     }
     
-    if (curve.size() > 2 && (1 - curve[curve.size()-2] <= step / (256.0 * 3))) {
+    if (curve.size() > 2 && (1 - curve[curve.size()-2] <= coord(step) / 3)) {
         curve.pop_back();
         curve.pop_back();
     }
@@ -169,6 +220,20 @@ void mappingToCurve(const std::vector<int> &mapping, std::vector<double> &curve)
     if (curve.size() < 4) {
         curve = { DCT_Linear }; // not enough points, fall back to linear
     } else {
+        CubicSplineCurve c(curve);
+        curve.pop_back();
+        curve.pop_back();
+        double gap = coord(step);
+        while (1 - curve[curve.size()-2] > gap) {
+            double x = curve[curve.size()-2] + gap;
+            if (1 - x <= gap / 3) {
+                break;
+            }
+            curve.push_back(x);
+            curve.push_back(c.getVal(x));
+        }
+        curve.push_back(1.0);
+        curve.push_back(1.0);
         curve.insert(curve.begin(), DCT_Spline);
     }
 }
@@ -220,7 +285,7 @@ void RawImageSource::getAutoMatchedToneCurve(const ColorManagementParams &cp, st
     neutral.icm = cp;
     neutral.raw.bayersensor.method = RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::Method::FAST);
     neutral.raw.xtranssensor.method = RAWParams::XTransSensor::getMethodString(RAWParams::XTransSensor::Method::FAST);
-    neutral.icm.outputProfile = "sRGB";
+    neutral.icm.outputProfile = ColorManagementParams::NoICMString;
     
     std::unique_ptr<IImage8> source;
     {