From 6ae287b72b8072293de19437f71e7a5039321e19 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Thu, 26 Jul 2018 17:27:57 +0200 Subject: [PATCH 01/10] diagonalcurves: replace cubic splines with (centripetal) Catmull-Rom splines See https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline Fixes #4343 --- rtengine/curves.h | 1 + rtengine/diagonalcurves.cc | 165 ++++++++++++++++++++++++++++++++++--- 2 files changed, 155 insertions(+), 11 deletions(-) diff --git a/rtengine/curves.h b/rtengine/curves.h index 98934f13a..a9f4b4ace 100644 --- a/rtengine/curves.h +++ b/rtengine/curves.h @@ -444,6 +444,7 @@ protected: DiagonalCurveType kind; void spline_cubic_set (); + void catmull_rom_set(); void NURBS_set (); public: diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index a3b00b165..1e9c730d5 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -86,7 +86,8 @@ DiagonalCurve::DiagonalCurve (const std::vector& p, int poly_pn) if (!identity) { if (kind == DCT_Spline && N > 2) { - spline_cubic_set (); + //spline_cubic_set (); + catmull_rom_set(); } else if (kind == DCT_NURBS && N > 2) { NURBS_set (); fillHash(); @@ -270,6 +271,134 @@ void DiagonalCurve::NURBS_set () fillDyByDx(); } + +/***************************************************************************** + * Catmull Rom Spline + * (https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline) + *****************************************************************************/ + +namespace { + +inline double pow2(double x) +{ + return x*x; +} + + +inline double catmull_rom_tj(double ti, + double xi, double yi, + double xj, double yj) +{ + return sqrt(sqrt(pow2(xj-xi) + pow2(yj-yi))) + ti; +} + + +inline void catmull_rom_spline(int n_points, + double p0_x, double p0_y, + double p1_x, double p1_y, + double p2_x, double p2_y, + double p3_x, double p3_y, + std::vector &res_x, + std::vector &res_y) +{ + res_x.reserve(n_points); + res_y.reserve(n_points); + + double t0 = 0; + double t1 = catmull_rom_tj(t0, p0_x, p0_y, p1_x, p1_y); + double t2 = catmull_rom_tj(t1, p1_x, p1_y, p2_x, p2_y); + double t3 = catmull_rom_tj(t2, p2_x, p2_y, p3_x, p3_y); + + double space = (t2-t1) / n_points; + + double t; + int i; + double c, d, A1_x, A1_y, A2_x, A2_y, A3_x, A3_y; + double B1_x, B1_y, B2_x, B2_y, C_x, C_y; + + for (i = 0; i < n_points; ++i) { + t = t1 + space * i; + + c = (t1 - t)/(t1 - t0); + d = (t - t0)/(t1 - t0); + A1_x = c * p0_x + d * p1_x; + A1_y = c * p0_y + d * p1_y; + + c = (t2 - t)/(t2 - t1); + d = (t - t1)/(t2 - t1); + A2_x = c * p1_x + d * p2_x; + A2_y = c * p1_y + d * p2_y; + + c = (t3 - t)/(t3 - t2); + d = (t - t2)/(t3 - t2); + A3_x = c * p2_x + d * p3_x; + A3_y = c * p2_y + d * p3_y; + + c = (t2 - t)/(t2 - t0); + d = (t - t0)/(t2 - t0); + B1_x = c * A1_x + d * A2_x; + B1_y = c * A1_y + d * A2_y; + + c = (t3 - t)/(t3 - t1); + d = (t - t1)/(t3 - t1); + B2_x = c * A2_x + d * A3_x; + B2_y = c * A2_y + d * A3_y; + + c = (t2 - t)/(t2 - t1); + d = (t - t1)/(t2 - t1); + C_x = c * B1_x + d * B2_x; + C_y = c * B1_y + d * B2_y; + + res_x.push_back(C_x); + res_y.push_back(C_y); + } +} + + +void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, + std::vector &res_x, std::vector &res_y) +{ + static const double epsilon = 1e-5; + double xr = x[1] - x[0]; + double yr = y[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]; + 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]; + + int segments = n_cp - 1; + int points_segments = n_points / segments; + + res_x.reserve(n_points); + res_y.reserve(n_points); + + for (int i = 0; i < segments; ++i) { + int n = points_segments + (i == 0 ? n_points % segments : 0); + catmull_rom_spline( + n, i == 0 ? x_first : x[i-1], i == 0 ? y_first : y[i-1], + x[i], y[i], x[i+1], y[i+1], + i == segments-1 ? x_last : x[i+2], + i == segments-1 ? y_last : y[i+2], + res_x, res_y); + } +} + +} // namespace + + +void DiagonalCurve::catmull_rom_set() +{ + int n_points = max(ppn * 65, 65000); + poly_x.clear(); + poly_y.clear(); + catmull_rom_chain(n_points, N, x, y, poly_x, poly_y); +} + +/*****************************************************************************/ + + double DiagonalCurve::getVal (double t) const { @@ -300,7 +429,8 @@ double DiagonalCurve::getVal (double t) const } case DCT_Linear : - case DCT_Spline : { + // case DCT_Spline : + { // values under and over the first and last point if (t > x[N - 1]) { return y[N - 1]; @@ -324,20 +454,33 @@ double DiagonalCurve::getVal (double t) const double h = x[k_hi] - x[k_lo]; // linear - if (kind == DCT_Linear) { + // if (kind == DCT_Linear) { return y[k_lo] + (t - x[k_lo]) * ( y[k_hi] - y[k_lo] ) / h; - } - // spline curve - else { // if (kind==Spline) { - 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 CLIPD(r); - } + // } + // // spline curve + // else { // if (kind==Spline) { + // 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 CLIPD(r); + // } break; } + case DCT_Spline: { + auto it = std::lower_bound(poly_x.begin(), poly_x.end(), t); + if (it == poly_x.end()) { + return poly_y.back(); + } + auto d = it - poly_x.begin(); + if (it+1 < poly_x.end() && t - *it > *(it+1) - t) { + ++d; + } + return *(poly_y.begin() + d); + break; + } + case DCT_NURBS : { // get the hash table entry by rounding the value (previously multiplied by "hashSize") unsigned short int i = (unsigned short int)(t * hashSize); From 2bb27661a2802f9dfefe12e94afc5c8f5e5103b0 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Fri, 10 Aug 2018 08:05:41 +0200 Subject: [PATCH 02/10] various tweaks to catmull-rom --- rtengine/diagonalcurves.cc | 12 +++--- rtengine/histmatching.cc | 79 ++++++++++++++++++++++++++++++++++---- 2 files changed, 79 insertions(+), 12 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index 1e9c730d5..23c1fba47 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -359,12 +359,14 @@ void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, std::vector &res_x, std::vector &res_y) { static const double epsilon = 1e-5; - double xr = x[1] - x[0]; - double yr = y[1] - y[0]; + // double xr = x[1] - x[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]; - yr = y[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]; 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; } diff --git a/rtengine/histmatching.cc b/rtengine/histmatching.cc index 2dfcc6d3b..184152ac8 100644 --- a/rtengine/histmatching.cc +++ b/rtengine/histmatching.cc @@ -97,6 +97,55 @@ int findMatch(int val, const std::vector &cdf, int j) } +class CubicSplineCurve: public DiagonalCurve { +public: + CubicSplineCurve(const std::vector &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 &mapping, std::vector &curve) { curve.clear(); @@ -117,8 +166,9 @@ void mappingToCurve(const std::vector &mapping, std::vector &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 &mapping, std::vector &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,23 +202,38 @@ void mappingToCurve(const std::vector &mapping, std::vector &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(); } - + curve.push_back(1.0); curve.push_back(1.0); 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 source; { From 0d7d1cfc8c0baf30ea6a148730a435615ba18026 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 11 Aug 2018 23:40:21 +0200 Subject: [PATCH 03/10] more tweaks to catmull-rom and histogram matching after the feedback by DrSlony --- rtengine/diagonalcurves.cc | 19 +++++++++---------- rtengine/histmatching.cc | 24 +++++++++++++++--------- 2 files changed, 24 insertions(+), 19 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index 23c1fba47..b2ac7414e 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -289,7 +289,8 @@ inline double catmull_rom_tj(double ti, double xi, double yi, double xj, double yj) { - return sqrt(sqrt(pow2(xj-xi) + pow2(yj-yi))) + ti; + static constexpr double alpha = 0.25; + return pow(sqrt(pow2(xj-xi) + pow2(yj-yi)), alpha) + ti; } @@ -303,7 +304,7 @@ inline void catmull_rom_spline(int n_points, { res_x.reserve(n_points); res_y.reserve(n_points); - + double t0 = 0; double t1 = catmull_rom_tj(t0, p0_x, p0_y, p1_x, p1_y); double t2 = catmull_rom_tj(t1, p1_x, p1_y, p2_x, p2_y); @@ -359,15 +360,13 @@ void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, std::vector &res_x, std::vector &res_y) { static const double epsilon = 1e-5; - // double xr = x[1] - x[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 xr = x[1] - x[0]; + double yr = y[1] - y[0]; + double x_first = x[0] - xr * 0.01; double y_first = xr > epsilon ? (yr / xr) * (x_first - x[0]) + y[0] : y[0]; - // xr = x[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; + xr = x[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.01; double y_last = xr > epsilon ? (yr / xr) * (x_last - x[0]) + y[0] : y[0]; int segments = n_cp - 1; diff --git a/rtengine/histmatching.cc b/rtengine/histmatching.cc index 184152ac8..d4f236d0e 100644 --- a/rtengine/histmatching.cc +++ b/rtengine/histmatching.cc @@ -221,20 +221,26 @@ void mappingToCurve(const std::vector &mapping, std::vector &curve) 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; - } + double mid = coord(idx); + double x = 0.0; + constexpr double shgap = 0.075; + curve = { DCT_Spline }; + while (mid - x > shgap / 2) { curve.push_back(x); curve.push_back(c.getVal(x)); + x += shgap; + } + curve.push_back(mid); + curve.push_back(c.getVal(mid)); + constexpr double hlgap = 0.2; + x = mid + hlgap; + while (1 - x > hlgap / 2) { + curve.push_back(x); + curve.push_back(c.getVal(x)); + x += hlgap; } curve.push_back(1.0); curve.push_back(1.0); - curve.insert(curve.begin(), DCT_Spline); } } From 2f2065cf9a45573fec974d678159c8cadfa04925 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sun, 12 Aug 2018 23:20:31 +0200 Subject: [PATCH 04/10] fixed typo leading to artifacts in catmull-rom splines in extreme cases --- rtengine/diagonalcurves.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index b2ac7414e..3955c8e01 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -367,7 +367,7 @@ void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, xr = x[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.01; - double y_last = xr > epsilon ? (yr / xr) * (x_last - x[0]) + y[0] : y[0]; + double y_last = xr > epsilon ? (yr / xr) * (x_last - x[n_cp-1]) + y[n_cp-1] : y[n_cp-1]; int segments = n_cp - 1; int points_segments = n_points / segments; From 97e73457ee07a2bfa93fb925adda1774f779e90d Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Mon, 13 Aug 2018 14:18:59 +0200 Subject: [PATCH 05/10] histmatching: use better spaced points for the generated curve --- rtengine/histmatching.cc | 19 +++++-------------- 1 file changed, 5 insertions(+), 14 deletions(-) diff --git a/rtengine/histmatching.cc b/rtengine/histmatching.cc index d4f236d0e..b2fe436fd 100644 --- a/rtengine/histmatching.cc +++ b/rtengine/histmatching.cc @@ -221,26 +221,17 @@ void mappingToCurve(const std::vector &mapping, std::vector &curve) curve = { DCT_Linear }; // not enough points, fall back to linear } else { CubicSplineCurve c(curve); - double mid = coord(idx); + double gap = 0.05; double x = 0.0; - constexpr double shgap = 0.075; curve = { DCT_Spline }; - while (mid - x > shgap / 2) { + while (x < 1.0) { curve.push_back(x); curve.push_back(c.getVal(x)); - x += shgap; - } - curve.push_back(mid); - curve.push_back(c.getVal(mid)); - constexpr double hlgap = 0.2; - x = mid + hlgap; - while (1 - x > hlgap / 2) { - curve.push_back(x); - curve.push_back(c.getVal(x)); - x += hlgap; + x += gap; + gap *= 1.4; } curve.push_back(1.0); - curve.push_back(1.0); + curve.push_back(c.getVal(1.0)); } } From d6ca3d65aada2cdf1ba076e65c7fe59e98a25924 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Fri, 17 Aug 2018 15:30:36 +0200 Subject: [PATCH 06/10] catmull-rom: ensure that the curve evaluation is exact at the control points --- rtengine/diagonalcurves.cc | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index 3955c8e01..fc65951fa 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -317,7 +317,10 @@ inline void catmull_rom_spline(int n_points, double c, d, A1_x, A1_y, A2_x, A2_y, A3_x, A3_y; double B1_x, B1_y, B2_x, B2_y, C_x, C_y; - for (i = 0; i < n_points; ++i) { + res_x.push_back(p1_x); + res_y.push_back(p1_y); + + for (i = 1; i < n_points-1; ++i) { t = t1 + space * i; c = (t1 - t)/(t1 - t0); @@ -353,6 +356,9 @@ inline void catmull_rom_spline(int n_points, res_x.push_back(C_x); res_y.push_back(C_y); } + + res_x.push_back(p2_x); + res_y.push_back(p2_y); } From ef57c5da004d75d74ee93efcfef7ac869baf9167 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Fri, 17 Aug 2018 15:52:55 +0200 Subject: [PATCH 07/10] catmull-rom: add special case for evaluating straight segments at 0 or 1 --- rtengine/diagonalcurves.cc | 65 ++++++++++++++++++++++---------------- 1 file changed, 37 insertions(+), 28 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index fc65951fa..e3308d56c 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -320,41 +320,50 @@ inline void catmull_rom_spline(int n_points, res_x.push_back(p1_x); res_y.push_back(p1_y); - for (i = 1; i < n_points-1; ++i) { - t = t1 + space * i; + // special case, a segment at 0 or 1 is computed exactly + if (p1_y == p2_y && (p1_y == 0 || p1_y == 1)) { + for (i = 1; i < n_points-1; ++i) { + t = p1_x + space * i; + res_x.push_back(t); + res_y.push_back(p1_y); + } + } else { + for (i = 1; i < n_points-1; ++i) { + t = t1 + space * i; - c = (t1 - t)/(t1 - t0); - d = (t - t0)/(t1 - t0); - A1_x = c * p0_x + d * p1_x; - A1_y = c * p0_y + d * p1_y; + c = (t1 - t)/(t1 - t0); + d = (t - t0)/(t1 - t0); + A1_x = c * p0_x + d * p1_x; + A1_y = c * p0_y + d * p1_y; - c = (t2 - t)/(t2 - t1); - d = (t - t1)/(t2 - t1); - A2_x = c * p1_x + d * p2_x; - A2_y = c * p1_y + d * p2_y; + c = (t2 - t)/(t2 - t1); + d = (t - t1)/(t2 - t1); + A2_x = c * p1_x + d * p2_x; + A2_y = c * p1_y + d * p2_y; - c = (t3 - t)/(t3 - t2); - d = (t - t2)/(t3 - t2); - A3_x = c * p2_x + d * p3_x; - A3_y = c * p2_y + d * p3_y; + c = (t3 - t)/(t3 - t2); + d = (t - t2)/(t3 - t2); + A3_x = c * p2_x + d * p3_x; + A3_y = c * p2_y + d * p3_y; - c = (t2 - t)/(t2 - t0); - d = (t - t0)/(t2 - t0); - B1_x = c * A1_x + d * A2_x; - B1_y = c * A1_y + d * A2_y; + c = (t2 - t)/(t2 - t0); + d = (t - t0)/(t2 - t0); + B1_x = c * A1_x + d * A2_x; + B1_y = c * A1_y + d * A2_y; - c = (t3 - t)/(t3 - t1); - d = (t - t1)/(t3 - t1); - B2_x = c * A2_x + d * A3_x; - B2_y = c * A2_y + d * A3_y; + c = (t3 - t)/(t3 - t1); + d = (t - t1)/(t3 - t1); + B2_x = c * A2_x + d * A3_x; + B2_y = c * A2_y + d * A3_y; - c = (t2 - t)/(t2 - t1); - d = (t - t1)/(t2 - t1); - C_x = c * B1_x + d * B2_x; - C_y = c * B1_y + d * B2_y; + c = (t2 - t)/(t2 - t1); + d = (t - t1)/(t2 - t1); + C_x = c * B1_x + d * B2_x; + C_y = c * B1_y + d * B2_y; - res_x.push_back(C_x); - res_y.push_back(C_y); + res_x.push_back(C_x); + res_y.push_back(C_y); + } } res_x.push_back(p2_x); From 681aabd0e238f52e2211b9fe007b8f1e01ae5f69 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 18 Aug 2018 12:42:21 +0200 Subject: [PATCH 08/10] catmull-rom: use reflection to calculate the boundary control points --- rtengine/diagonalcurves.cc | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index e3308d56c..d543dd8f5 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -289,7 +289,7 @@ inline double catmull_rom_tj(double ti, double xi, double yi, double xj, double yj) { - static constexpr double alpha = 0.25; + static constexpr double alpha = 0.5; return pow(sqrt(pow2(xj-xi) + pow2(yj-yi)), alpha) + ti; } @@ -371,18 +371,23 @@ inline void catmull_rom_spline(int n_points, } +inline void catmull_rom_reflect(double px, double py, double cx, double cy, + double &rx, double &ry) +{ + double dx = px - cx; + double dy = py - cy; + rx = cx - dx; + ry = cy - dy; +} + + void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, std::vector &res_x, std::vector &res_y) { - static const double epsilon = 1e-5; - double xr = x[1] - x[0]; - double yr = y[1] - y[0]; - double x_first = x[0] - xr * 0.01; - double y_first = xr > epsilon ? (yr / xr) * (x_first - x[0]) + y[0] : y[0]; - xr = x[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.01; - double y_last = xr > epsilon ? (yr / xr) * (x_last - x[n_cp-1]) + y[n_cp-1] : y[n_cp-1]; + double x_first, y_first; + double x_last, y_last; + catmull_rom_reflect(x[1], y[1], x[0], y[0], x_first, y_first); + catmull_rom_reflect(x[n_cp-2], y[n_cp-2], x[n_cp-1], y[n_cp-1], x_last, y_last); int segments = n_cp - 1; int points_segments = n_points / segments; From e3ea0926c22792edda2e4dffee3d59c9a7fb7c30 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 18 Aug 2018 12:43:14 +0200 Subject: [PATCH 09/10] catmull-rom: use uniform spacing of curve evaluation points --- rtengine/diagonalcurves.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index d543dd8f5..1af086925 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -390,13 +390,12 @@ void catmull_rom_chain(int n_points, int n_cp, double *x, double *y, catmull_rom_reflect(x[n_cp-2], y[n_cp-2], x[n_cp-1], y[n_cp-1], x_last, y_last); int segments = n_cp - 1; - int points_segments = n_points / segments; res_x.reserve(n_points); res_y.reserve(n_points); for (int i = 0; i < segments; ++i) { - int n = points_segments + (i == 0 ? n_points % segments : 0); + int n = max(int(n_points * (x[i+1] - x[i]) + 0.5), 2); catmull_rom_spline( n, i == 0 ? x_first : x[i-1], i == 0 ? y_first : y[i-1], x[i], y[i], x[i+1], y[i+1], From 4b2392e44ad08d6a923ecc66bf8d6c13f81f8481 Mon Sep 17 00:00:00 2001 From: Alberto Griggio Date: Sat, 18 Aug 2018 15:11:54 +0200 Subject: [PATCH 10/10] catmull-rom: make curves more rounded See https://github.com/Beep6581/RawTherapee/pull/4701#issuecomment-414054187 --- rtengine/diagonalcurves.cc | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/rtengine/diagonalcurves.cc b/rtengine/diagonalcurves.cc index 1af086925..f478ba719 100644 --- a/rtengine/diagonalcurves.cc +++ b/rtengine/diagonalcurves.cc @@ -289,7 +289,8 @@ inline double catmull_rom_tj(double ti, double xi, double yi, double xj, double yj) { - static constexpr double alpha = 0.5; + // see https://github.com/Beep6581/RawTherapee/pull/4701#issuecomment-414054187 + static constexpr double alpha = 0.375; return pow(sqrt(pow2(xj-xi) + pow2(yj-yi)), alpha) + ti; } @@ -374,10 +375,19 @@ inline void catmull_rom_spline(int n_points, inline void catmull_rom_reflect(double px, double py, double cx, double cy, double &rx, double &ry) { +#if 0 double dx = px - cx; double dy = py - cy; rx = cx - dx; ry = cy - dy; +#else + // see https://github.com/Beep6581/RawTherapee/pull/4701#issuecomment-414054187 + static constexpr double epsilon = 1e-5; + double dx = px - cx; + double dy = py - cy; + rx = cx - dx * 0.01; + ry = dx > epsilon ? (dy / dx) * (rx - cx) + cy : cy; +#endif }