liz/rawTherapee
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #4701 from Beep6581/catmull-rom-spline

Browse Source 
diagonalcurves: replace cubic splines with (centripetal) Catmull-Rom splines
This commit is contained in:
Alberto Griggio
2018-08-18 15:21:15 +02:00
committed by GitHub
parent 7aeef817e8 4b2392e44a
commit a6bc6a0ab5
3 changed files with 255 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
rtengine/curves.h
View File

@@ -444,6 +444,7 @@ protected:
DiagonalCurveType kind;
void spline_cubic_set ();
void catmull_rom_set();
void NURBS_set ();
public:

195
rtengine/diagonalcurves.cc
View File

@@ -86,7 +86,8 @@ DiagonalCurve::DiagonalCurve (const std::vector<double>& 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,164 @@ 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)
{
// 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;
}
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<double> &res_x,
std::vector<double> &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;
res_x.push_back(p1_x);
res_y.push_back(p1_y);
// 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 = (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);
}
}
res_x.push_back(p2_x);
res_y.push_back(p2_y);
}
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
}
void catmull_rom_chain(int n_points, int n_cp, double *x, double *y,
std::vector<double> &res_x, std::vector<double> &res_y)
{
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;
res_x.reserve(n_points);
res_y.reserve(n_points);
for (int i = 0; i < segments; ++i) {
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],
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 +459,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 +484,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 LIM01(*(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);

78
rtengine/histmatching.cc
View File

@@ -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,24 +202,36 @@ 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();
}
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 {
curve.insert(curve.begin(), DCT_Spline);
CubicSplineCurve c(curve);
double gap = 0.05;
double x = 0.0;
curve = { DCT_Spline };
while (x < 1.0) {
curve.push_back(x);
curve.push_back(c.getVal(x));
x += gap;
gap *= 1.4;
}
curve.push_back(1.0);
curve.push_back(c.getVal(1.0));
}
}
@@ -220,7 +282,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;
{