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

Merge pull request #4277 from Beep6581/perceptual_curve_speedup

Browse Source 
Perceptual tone curve speedup
This commit is contained in:
Ingo Weyrich 2018-01-06 00:09:54 +01:00 committed by GitHub
commit 2d37fd6676
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 282 additions and 246 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

72
rtengine/ciecam02.cc
View File

@ -608,31 +608,31 @@ void Ciecam02::calculate_ab ( double &aa, double &bb, double h, double e, double
} }
void Ciecam02::calculate_abfloat ( float &aa, float &bb, float h, float e, float t, float nbb, float a ) void Ciecam02::calculate_abfloat ( float &aa, float &bb, float h, float e, float t, float nbb, float a )
{ {
float2 sincosval = xsincosf ((h * rtengine::RT_PI) / 180.0f); float2 sincosval = xsincosf(h * rtengine::RT_PI_F_180);
float sinh = sincosval.x; float sinh = sincosval.x;
float cosh = sincosval.y; float cosh = sincosval.y;
float x = (a / nbb) + 0.305f; float x = (a / nbb) + 0.305f;
float p3 = 1.05f; constexpr float p3 = 1.05f;
bool swapValues = fabs ( sinh ) > fabs ( cosh ); const bool swapValues = fabs(sinh) > fabs(cosh);
if (swapValues) { if (swapValues) {
std::swap (sinh, cosh); std::swap(sinh, cosh);
} }
float c1 = 1.f; float c1 = 1.f;
float c2 = sinh / cosh; float c2 = sinh / cosh;
if (swapValues) { if (swapValues) {
std::swap (c1, c2); std::swap(c1, c2);
} }
float div = ((e / (t * cosh)) - (-0.31362f - (p3 * 0.15681f)) * c1 - ((0.01924f - (p3 * 4.49038f)) * (c2))); float div = ((e / (t * cosh)) - (-0.31362f - (p3 * 0.15681f)) * c1 - ((0.01924f - (p3 * 4.49038f)) * c2));
// for large values of t the above calculation can change its sign which results in a hue shift of 180 degree // for large values of t the above calculation can change its sign which results in a hue shift of 180 degree
// so we have to check the sign to avoid this shift. // so we have to check the sign to avoid this shift.
// Additionally it seems useful to limit the minimum value of div // Additionally it seems useful to limit the minimum value of div
// I limited it, but I'm sure the actual limit is not the best one // I limited it, but I'm sure the actual limit is not the best one
if (signf (div) != signf (cosh) || fabsf (div) <= fabsf (cosh) * 2.f) { if (signf(div) != signf(cosh) || fabsf(div) <= fabsf(cosh) * 2.f) {
div = cosh * 2.f; div = cosh * 2.f;
} }
@ -640,7 +640,7 @@ void Ciecam02::calculate_abfloat ( float &aa, float &bb, float h, float e, float
bb = (aa * sinh) / cosh; bb = (aa * sinh) / cosh;
if (swapValues) { if (swapValues) {
std::swap (aa, bb); std::swap(aa, bb);
} }
} }
#ifdef __SSE2__ #ifdef __SSE2__
@ -1007,9 +1007,18 @@ void Ciecam02::xyz2jch_ciecam02float ( float &J, float &C, float &h, float aw, f
bp = MAXR (bp, 0.0f); bp = MAXR (bp, 0.0f);
} }
rpa = nonlinear_adaptationfloat ( rp, fl ); #ifdef __SSE2__
gpa = nonlinear_adaptationfloat ( gp, fl ); vfloat pv = _mm_setr_ps(rp, gp, bp, 1.f);
bpa = nonlinear_adaptationfloat ( bp, fl ); vfloat fv = F2V(fl);
vfloat outv = nonlinear_adaptationfloat(pv, fv);
rpa = outv[0];
gpa = outv[1];
bpa = outv[2];
#else
rpa = nonlinear_adaptationfloat(rp, fl);
gpa = nonlinear_adaptationfloat(gp, fl);
bpa = nonlinear_adaptationfloat(bp, fl);
#endif
ca = rpa - ((12.0f * gpa) - bpa) / 11.0f; ca = rpa - ((12.0f * gpa) - bpa) / 11.0f;
cb = (0.11111111f) * (rpa + gpa - (2.0f * bpa)); cb = (0.11111111f) * (rpa + gpa - (2.0f * bpa));
@ -1084,26 +1093,43 @@ void Ciecam02::jch2xyz_ciecam02float ( float &x, float &y, float &z, float J, fl
float a, ca, cb; float a, ca, cb;
float e, t; float e, t;
gamu = 1; gamu = 1;
xyz_to_cat02float ( rw, gw, bw, xw, yw, zw, gamu ); xyz_to_cat02float(rw, gw, bw, xw, yw, zw, gamu);
e = ((961.53846f) * nc * ncb) * (xcosf ( ((h * rtengine::RT_PI) / 180.0f) + 2.0f ) + 3.8f); e = ((961.53846f) * nc * ncb) * (xcosf(h * rtengine::RT_PI_F_180 + 2.0f) + 3.8f);
a = pow_F ( J / 100.0f, 1.0f / (c * cz) ) * aw;
t = pow_F ( 10.f * C / (sqrtf ( J ) * pow1), 1.1111111f );
calculate_abfloat ( ca, cb, h, e, t, nbb, a ); #ifdef __SSE2__
Aab_to_rgbfloat ( rpa, gpa, bpa, a, ca, cb, nbb ); vfloat powinv1 = _mm_setr_ps(J / 100.0f, 10.f * C / (sqrtf(J) * pow1), 1.f, 1.f);
vfloat powinv2 = _mm_setr_ps(1.0f / (c * cz), 1.1111111f, 1.f, 1.f);
vfloat powoutv = pow_F(powinv1, powinv2);
a = powoutv[0] * aw;
t = powoutv[1];
#else
a = pow_F(J / 100.0f, 1.0f / (c * cz)) * aw;
t = pow_F(10.f * C / (sqrtf(J) * pow1), 1.1111111f);
#endif
rp = inverse_nonlinear_adaptationfloat ( rpa, fl ); calculate_abfloat(ca, cb, h, e, t, nbb, a);
gp = inverse_nonlinear_adaptationfloat ( gpa, fl ); Aab_to_rgbfloat(rpa, gpa, bpa, a, ca, cb, nbb);
bp = inverse_nonlinear_adaptationfloat ( bpa, fl );
hpe_to_xyzfloat ( x, y, z, rp, gp, bp ); #ifdef __SSE2__
xyz_to_cat02float ( rc, gc, bc, x, y, z, gamu ); vfloat pav = _mm_setr_ps(rpa, gpa, bpa, 1.f);
vfloat fv = F2V(fl);
vfloat outv = inverse_nonlinear_adaptationfloat(pav, fv);
rp = outv[0];
gp = outv[1];
bp = outv[2];
#else
rp = inverse_nonlinear_adaptationfloat(rpa, fl);
gp = inverse_nonlinear_adaptationfloat(gpa, fl);
bp = inverse_nonlinear_adaptationfloat(bpa, fl);
#endif
hpe_to_xyzfloat(x, y, z, rp, gp, bp);
xyz_to_cat02float(rc, gc, bc, x, y, z, gamu);
r = rc / (((yw * d) / rw) + (1.0f - d)); r = rc / (((yw * d) / rw) + (1.0f - d));
g = gc / (((yw * d) / gw) + (1.0f - d)); g = gc / (((yw * d) / gw) + (1.0f - d));
b = bc / (((yw * d) / bw) + (1.0f - d)); b = bc / (((yw * d) / bw) + (1.0f - d));
cat02_to_xyzfloat ( x, y, z, r, g, b, gamu ); cat02_to_xyzfloat(x, y, z, r, g, b, gamu);
} }
#ifdef __SSE2__ #ifdef __SSE2__

82
rtengine/curves.cc
View File

@ -1822,9 +1822,14 @@ float PerceptualToneCurve::calculateToneCurveContrastValue() const
return maxslope; return maxslope;
} }
void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurveState & state) const void PerceptualToneCurve::BatchApply(const size_t start, const size_t end, float *rc, float *gc, float *bc, const PerceptualToneCurveState &state) const
{ {
float x, y, z; const AdobeToneCurve& adobeTC = static_cast<const AdobeToneCurve&>((const ToneCurve&) * this);
for (size_t i = start; i < end; ++i) {
float r = CLIP(rc[i]);
float g = CLIP(gc[i]);
float b = CLIP(bc[i]);
if (!state.isProphoto) { if (!state.isProphoto) {
// convert to prophoto space to make sure the same result is had regardless of working color space // convert to prophoto space to make sure the same result is had regardless of working color space
@ -1836,38 +1841,38 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
b = newb; b = newb;
} }
const AdobeToneCurve& adobeTC = static_cast<const AdobeToneCurve&>((const ToneCurve&) * this);
float ar = r; float ar = r;
float ag = g; float ag = g;
float ab = b; float ab = b;
adobeTC.Apply(ar, ag, ab); adobeTC.Apply(ar, ag, ab);
if (ar >= 65535.f && ag >= 65535.f && ab >= 65535.f) { if (ar >= 65535.f && ag >= 65535.f && ab >= 65535.f) {
// clip fast path, will also avoid strange colors of clipped highlights // clip fast path, will also avoid strange colours of clipped highlights
r = g = b = 65535.f; rc[i] = gc[i] = bc[i] = 65535.f;
return; continue;
} }
if (ar <= 0.f && ag <= 0.f && ab <= 0.f) { if (ar <= 0.f && ag <= 0.f && ab <= 0.f) {
r = g = b = 0; rc[i] = gc[i] = bc[i] = 0;
return; continue;
} }
// ProPhoto constants for luminance, that is xyz_prophoto[1][] // ProPhoto constants for luminance, that is xyz_prophoto[1][]
const float Yr = 0.2880402f; constexpr float Yr = 0.2880402f;
const float Yg = 0.7118741f; constexpr float Yg = 0.7118741f;
const float Yb = 0.0000857f; constexpr float Yb = 0.0000857f;
// we use the Adobe (RGB-HSV hue-stabilized) curve to decide luminance, which generally leads to a less contrasty result // we use the Adobe (RGB-HSV hue-stabilized) curve to decide luminance, which generally leads to a less contrasty result
// compared to a pure luminance curve. We do this to be more compatible with the most popular curves. // compared to a pure luminance curve. We do this to be more compatible with the most popular curves.
float oldLuminance = r * Yr + g * Yg + b * Yb; const float oldLuminance = r * Yr + g * Yg + b * Yb;
float newLuminance = ar * Yr + ag * Yg + ab * Yb; const float newLuminance = ar * Yr + ag * Yg + ab * Yb;
float Lcoef = newLuminance / oldLuminance; const float Lcoef = newLuminance / oldLuminance;
r = LIM<float>(r * Lcoef, 0.f, 65535.f); r = LIM<float>(r * Lcoef, 0.f, 65535.f);
g = LIM<float>(g * Lcoef, 0.f, 65535.f); g = LIM<float>(g * Lcoef, 0.f, 65535.f);
b = LIM<float>(b * Lcoef, 0.f, 65535.f); b = LIM<float>(b * Lcoef, 0.f, 65535.f);
// move to JCh so we can modulate chroma based on the global contrast-related chroma scaling factor // move to JCh so we can modulate chroma based on the global contrast-related chroma scaling factor
float x, y, z;
Color::Prophotoxyz(r, g, b, x, y, z); Color::Prophotoxyz(r, g, b, x, y, z);
float J, C, h; float J, C, h;
@ -1879,7 +1884,7 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
if (!isfinite(J) || !isfinite(C) || !isfinite(h)) { if (!isfinite(J) || !isfinite(C) || !isfinite(h)) {
// this can happen for dark noise colors or colors outside human gamut. Then we just return the curve's result. // this can happen for dark noise colours or colours outside human gamut. Then we just return the curve's result.
if (!state.isProphoto) { if (!state.isProphoto) {
float newr = state.Prophoto2Working[0][0] * r + state.Prophoto2Working[0][1] * g + state.Prophoto2Working[0][2] * b; float newr = state.Prophoto2Working[0][0] * r + state.Prophoto2Working[0][1] * g + state.Prophoto2Working[0][2] * b;
float newg = state.Prophoto2Working[1][0] * r + state.Prophoto2Working[1][1] * g + state.Prophoto2Working[1][2] * b; float newg = state.Prophoto2Working[1][0] * r + state.Prophoto2Working[1][1] * g + state.Prophoto2Working[1][2] * b;
@ -1888,8 +1893,11 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
g = newg; g = newg;
b = newb; b = newb;
} }
rc[i] = r;
gc[i] = g;
bc[i] = b;
return; continue;
} }
float cmul = state.cmul_contrast; // chroma scaling factor float cmul = state.cmul_contrast; // chroma scaling factor
@ -1899,8 +1907,8 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
{ {
// decrease chroma scaling sligthly of extremely saturated colors // decrease chroma scaling sligthly of extremely saturated colors
float saturated_scale_factor = 0.95f; float saturated_scale_factor = 0.95f;
const float lolim = 35.f; // lower limit, below this chroma all colors will keep original chroma scaling factor constexpr float lolim = 35.f; // lower limit, below this chroma all colors will keep original chroma scaling factor
const float hilim = 60.f; // high limit, above this chroma the chroma scaling factor is multiplied with the saturated scale factor value above constexpr float hilim = 60.f; // high limit, above this chroma the chroma scaling factor is multiplied with the saturated scale factor value above
if (C < lolim) { if (C < lolim) {
// chroma is low enough, don't scale // chroma is low enough, don't scale
@ -1928,8 +1936,8 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
float nL = Color::gamma2curve[newLuminance]; // apply gamma so we make comparison and transition with a more perceptual lightness scale float nL = Color::gamma2curve[newLuminance]; // apply gamma so we make comparison and transition with a more perceptual lightness scale
float dark_scale_factor = 1.20f; float dark_scale_factor = 1.20f;
//float dark_scale_factor = 1.0 + state.debug.p2 / 100.0f; //float dark_scale_factor = 1.0 + state.debug.p2 / 100.0f;
const float lolim = 0.15f; constexpr float lolim = 0.15f;
const float hilim = 0.50f; constexpr float hilim = 0.50f;
if (nL < lolim) { if (nL < lolim) {
// do nothing, keep scale factor // do nothing, keep scale factor
@ -1954,8 +1962,8 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
{ {
// to avoid strange CIECAM02 chroma errors on close-to-shadow-clipping colors we reduce chroma scaling towards 1.0 for black colors // to avoid strange CIECAM02 chroma errors on close-to-shadow-clipping colors we reduce chroma scaling towards 1.0 for black colors
float dark_scale_factor = 1.f / cmul; float dark_scale_factor = 1.f / cmul;
const float lolim = 4.f; constexpr float lolim = 4.f;
const float hilim = 7.f; constexpr float hilim = 7.f;
if (J < lolim) { if (J < lolim) {
// do nothing, keep scale factor // do nothing, keep scale factor
@ -1985,7 +1993,7 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
c, nc, 1, pow1, nbb, ncb, fl, cz, d, aw ); c, nc, 1, pow1, nbb, ncb, fl, cz, d, aw );
if (!isfinite(x) || !isfinite(y) || !isfinite(z)) { if (!isfinite(x) || !isfinite(y) || !isfinite(z)) {
// can happen for colors on the rim of being outside gamut, that worked without chroma scaling but not with. Then we return only the curve's result. // can happen for colours on the rim of being outside gamut, that worked without chroma scaling but not with. Then we return only the curve's result.
if (!state.isProphoto) { if (!state.isProphoto) {
float newr = state.Prophoto2Working[0][0] * r + state.Prophoto2Working[0][1] * g + state.Prophoto2Working[0][2] * b; float newr = state.Prophoto2Working[0][0] * r + state.Prophoto2Working[0][1] * g + state.Prophoto2Working[0][2] * b;
float newg = state.Prophoto2Working[1][0] * r + state.Prophoto2Working[1][1] * g + state.Prophoto2Working[1][2] * b; float newg = state.Prophoto2Working[1][0] * r + state.Prophoto2Working[1][1] * g + state.Prophoto2Working[1][2] * b;
@ -1995,7 +2003,11 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
b = newb; b = newb;
} }
return; rc[i] = r;
gc[i] = g;
bc[i] = b;
continue;
} }
Color::xyz2Prophoto(x, y, z, r, g, b); Color::xyz2Prophoto(x, y, z, r, g, b);
@ -2012,14 +2024,13 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
// we use the RGB-HSV hue-stable "Adobe" curve as reference. For S-curve contrast it increases // we use the RGB-HSV hue-stable "Adobe" curve as reference. For S-curve contrast it increases
// saturation greatly, but desaturates extreme highlights and thus provide a smooth transition to // saturation greatly, but desaturates extreme highlights and thus provide a smooth transition to
// the white point. However the desaturation effect is quite strong so we make a weighting // the white point. However the desaturation effect is quite strong so we make a weighting
float ah, as, av, h, s, v; const float as = Color::rgb2s(ar, ag, ab);
Color::rgb2hsv(ar, ag, ab, ah, as, av); const float s = Color::rgb2s(r, g, b);
Color::rgb2hsv(r, g, b, h, s, v);
float sat_scale = as <= 0.f ? 1.f : s / as; // saturation scale compared to Adobe curve const float sat_scale = as <= 0.f ? 1.f : s / as; // saturation scale compared to Adobe curve
float keep = 0.2f; float keep = 0.2f;
const float lolim = 1.00f; // only mix in the Adobe curve if we have increased saturation compared to it constexpr float lolim = 1.00f; // only mix in the Adobe curve if we have increased saturation compared to it
const float hilim = 1.20f; constexpr float hilim = 1.20f;
if (sat_scale < lolim) { if (sat_scale < lolim) {
// saturation is low enough, don't desaturate // saturation is low enough, don't desaturate
@ -2041,9 +2052,9 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
if (keep < 1.f) { if (keep < 1.f) {
// mix in some of the Adobe curve result // mix in some of the Adobe curve result
r = r * keep + (1.f - keep) * ar; r = intp(keep, r, ar);
g = g * keep + (1.f - keep) * ag; g = intp(keep, g, ag);
b = b * keep + (1.f - keep) * ab; b = intp(keep, b, ab);
} }
} }
@ -2055,8 +2066,11 @@ void PerceptualToneCurve::Apply(float &r, float &g, float &b, PerceptualToneCurv
g = newg; g = newg;
b = newb; b = newb;
} }
rc[i] = r;
gc[i] = g;
bc[i] = b;
}
} }
float PerceptualToneCurve::cf_range[2]; float PerceptualToneCurve::cf_range[2];
float PerceptualToneCurve::cf[1000]; float PerceptualToneCurve::cf[1000];
float PerceptualToneCurve::f, PerceptualToneCurve::c, PerceptualToneCurve::nc, PerceptualToneCurve::yb, PerceptualToneCurve::la, PerceptualToneCurve::xw, PerceptualToneCurve::yw, PerceptualToneCurve::zw, PerceptualToneCurve::gamut; float PerceptualToneCurve::f, PerceptualToneCurve::c, PerceptualToneCurve::nc, PerceptualToneCurve::yb, PerceptualToneCurve::la, PerceptualToneCurve::xw, PerceptualToneCurve::yw, PerceptualToneCurve::zw, PerceptualToneCurve::gamut;

2
rtengine/curves.h
View File

@ -872,7 +872,7 @@ private:
public: public:
static void init(); static void init();
void initApplyState(PerceptualToneCurveState & state, Glib::ustring workingSpace) const; void initApplyState(PerceptualToneCurveState & state, Glib::ustring workingSpace) const;
void Apply(float& r, float& g, float& b, PerceptualToneCurveState & state) const; void BatchApply(const size_t start, const size_t end, float *r, float *g, float *b, const PerceptualToneCurveState &state) const;
}; };
// Standard tone curve // Standard tone curve

8
rtengine/improcfun.cc
View File

@ -233,14 +233,8 @@ void customToneCurve(const ToneCurve &customToneCurve, ToneCurveParams::TcMode c
} }
} else if (curveMode == ToneCurveParams::TcMode::PERCEPTUAL) { // apply curve while keeping color appearance constant } else if (curveMode == ToneCurveParams::TcMode::PERCEPTUAL) { // apply curve while keeping color appearance constant
const PerceptualToneCurve& userToneCurve = static_cast<const PerceptualToneCurve&> (customToneCurve); const PerceptualToneCurve& userToneCurve = static_cast<const PerceptualToneCurve&> (customToneCurve);
for (int i = istart, ti = 0; i < tH; i++, ti++) { for (int i = istart, ti = 0; i < tH; i++, ti++) {
for (int j = jstart, tj = 0; j < tW; j++, tj++) { userToneCurve.BatchApply(0, tW - jstart, &rtemp[ti * tileSize], &gtemp[ti * tileSize], &btemp[ti * tileSize], ptcApplyState);
rtemp[ti * tileSize + tj] = CLIP<float> (rtemp[ti * tileSize + tj]);
gtemp[ti * tileSize + tj] = CLIP<float> (gtemp[ti * tileSize + tj]);
btemp[ti * tileSize + tj] = CLIP<float> (btemp[ti * tileSize + tj]);
userToneCurve.Apply(rtemp[ti * tileSize + tj], gtemp[ti * tileSize + tj], btemp[ti * tileSize + tj], ptcApplyState);
}
} }
} }
} }

2
rtengine/rt_math.h
View File

@ -14,6 +14,7 @@ constexpr double MAXVALD = static_cast<double>(MAXVAL); // double version of MAX
constexpr double RT_PI = 3.14159265358979323846; // pi constexpr double RT_PI = 3.14159265358979323846; // pi
constexpr double RT_PI_2 = 1.57079632679489661923; // pi/2 constexpr double RT_PI_2 = 1.57079632679489661923; // pi/2
constexpr double RT_PI_180 = 0.017453292519943295769; // pi/180
constexpr double RT_1_PI = 0.31830988618379067154; // 1/pi constexpr double RT_1_PI = 0.31830988618379067154; // 1/pi
constexpr double RT_2_PI = 0.63661977236758134308; // 2/pi constexpr double RT_2_PI = 0.63661977236758134308; // 2/pi
constexpr double RT_SQRT1_2 = 0.70710678118654752440; // 1/sqrt(2) constexpr double RT_SQRT1_2 = 0.70710678118654752440; // 1/sqrt(2)
@ -23,6 +24,7 @@ constexpr double RT_NAN = std::numeric_limits<double>::quiet_NaN();
constexpr float RT_PI_F = RT_PI; constexpr float RT_PI_F = RT_PI;
constexpr float RT_PI_F_2 = RT_PI_2; constexpr float RT_PI_F_2 = RT_PI_2;
constexpr float RT_PI_F_180 = RT_PI_180;
constexpr float RT_INFINITY_F = std::numeric_limits<float>::infinity(); constexpr float RT_INFINITY_F = std::numeric_limits<float>::infinity();
constexpr float RT_NAN_F = std::numeric_limits<float>::quiet_NaN(); constexpr float RT_NAN_F = std::numeric_limits<float>::quiet_NaN();