diff --git a/rtengine/dcrop.cc b/rtengine/dcrop.cc index 427b36cef..4a54fdcee 100644 --- a/rtengine/dcrop.cc +++ b/rtengine/dcrop.cc @@ -166,15 +166,19 @@ void Crop::update (int todo, bool internal) { if (skip==1) { parent->ipf.lumadenoise (labnCrop, cbuffer); parent->ipf.sharpening (labnCrop, (unsigned short**)cbuffer); + parent->ipf.waveletEqualizer(labnCrop, true, false); } } // apply color operations if (todo & M_COLOR) { parent->ipf.colorCurve (laboCrop, labnCrop); - if (skip==1) + if (skip==1) { parent->ipf.colordenoise (labnCrop, cbuffer); + parent->ipf.waveletEqualizer(labnCrop, false, true); + } } + // switch back to rgb parent->ipf.lab2rgb (labnCrop, cropImg); diff --git a/rtengine/improccoordinator.cc b/rtengine/improccoordinator.cc index 6633aa0d2..fc66c93e9 100644 --- a/rtengine/improccoordinator.cc +++ b/rtengine/improccoordinator.cc @@ -171,9 +171,14 @@ void ImProcCoordinator::updatePreviewImage (int todo) { progress ("Sharpening...",100*readyphase/numofphases); ipf.sharpening (nprevl, (unsigned short**)buffer); } + if (scale==1) { + progress ("Wavelet...",100*readyphase/numofphases); + ipf.waveletEqualizer (nprevl, true, false); + } readyphase++; } + if (todo & M_COLOR) { progress ("Applying Color Boost...",100*readyphase/numofphases); ipf.colorCurve (oprevl, nprevl); @@ -182,6 +187,10 @@ void ImProcCoordinator::updatePreviewImage (int todo) { progress ("Denoising color...",100*readyphase/numofphases); ipf.colordenoise (nprevl, buffer); } + if (scale==1) { + progress ("Wavelet...",100*readyphase/numofphases); + ipf.waveletEqualizer (nprevl, false, true); + } readyphase++; } diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h index c81e4dec8..d82f5f6d9 100644 --- a/rtengine/improcfun.h +++ b/rtengine/improcfun.h @@ -89,7 +89,8 @@ class ImProcFunctions { void lab2rgb (LabImage* lab, Image8* image); void resize (Image16* src, Image16* dst); void deconvsharpening(LabImage* lab, unsigned short** buffer); - void waveletEqualizer(Image16 * image, int fw, int fh, const EqualizerParams & params); + void waveletEqualizer(Image16 * image); + void waveletEqualizer(LabImage * image, bool luminance, bool chromaticity); Image8* lab2rgb (LabImage* lab, int cx, int cy, int cw, int ch, Glib::ustring profile); Image16* lab2rgb16 (LabImage* lab, int cx, int cy, int cw, int ch, Glib::ustring profile); diff --git a/rtengine/ipequalizer.cc b/rtengine/ipequalizer.cc index 298c82b3e..8371605f3 100644 --- a/rtengine/ipequalizer.cc +++ b/rtengine/ipequalizer.cc @@ -22,18 +22,55 @@ #include +#include + namespace rtengine { -void ImProcFunctions :: waveletEqualizer(Image16 * image, int fw, int fh, const EqualizerParams & params) { +void ImProcFunctions :: waveletEqualizer(Image16 * image) { - wavelet_decomposition r(image->r, fw, fh); - r.reconstruct(image->r, params.c); + if (!params->equalizer.enabled) { + return; + } - wavelet_decomposition g(image->g, fw, fh); - g.reconstruct(image->g, params.c); + limiter l(0, 65535); - wavelet_decomposition b(image->b, fw, fh); - b.reconstruct(image->b, params.c); + wavelet_decomposition r(image->r, image->getW(), image->getH()); + r.reconstruct(image->r, params->equalizer.c, l); + + wavelet_decomposition g(image->g, image->getW(), image->getH()); + g.reconstruct(image->g, params->equalizer.c, l); + + wavelet_decomposition b(image->b, image->getW(), image->getH()); + b.reconstruct(image->b, params->equalizer.c, l); + +} + +void ImProcFunctions :: waveletEqualizer (LabImage * image, bool luminance, bool chromaticity) { + + if (!params->equalizer.enabled) { + return; + } + + clock_t start = clock(); + + if (luminance) { + limiter l1(0, 65535); + + wavelet_decomposition L(image->L, image->W, image->H); + L.reconstruct(image->L, params->equalizer.c, l1); + } + + if (chromaticity) { + limiter l2(-32768, 32767); + + wavelet_decomposition a(image->a, image->W, image->H); + a.reconstruct(image->a, params->equalizer.c, l2); + + wavelet_decomposition b(image->b, image->W, image->H); + b.reconstruct(image->b, params->equalizer.c, l2); + } + + std::cout << "Wavelets done in " << (double)(clock() - start) / CLOCKS_PER_SEC << std::endl; } diff --git a/rtengine/refreshmap.cc b/rtengine/refreshmap.cc index b13674b45..a27ae2bb7 100644 --- a/rtengine/refreshmap.cc +++ b/rtengine/refreshmap.cc @@ -103,6 +103,7 @@ ALL, // EvResizeEnabled ALL, // EvProfileChangeNotification RETINEX, // EvShrHighQuality TRANSFORM, // EvPerspCorr -ALL // ????? EvEqualizer +EQUALIZER, // EvEqualizer +EQUALIZER // EvEqlEnabled }; diff --git a/rtengine/refreshmap.h b/rtengine/refreshmap.h index 6a9f9700f..962bbcf59 100644 --- a/rtengine/refreshmap.h +++ b/rtengine/refreshmap.h @@ -36,6 +36,7 @@ #define CROP 16384 #define EXIF 32768 #define IPTC 32768 +#define EQUALIZER 3 #define NONE 0 #define M_INIT 128 diff --git a/rtengine/simpleprocess.cc b/rtengine/simpleprocess.cc index 57ebd889a..b64238e3b 100644 --- a/rtengine/simpleprocess.cc +++ b/rtengine/simpleprocess.cc @@ -103,10 +103,6 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p baseImg = trImg; } - if (params.equalizer.enabled) { - ipf.waveletEqualizer (baseImg, fw, fh, params.equalizer); - } - // update blurmap int** buffer = new int*[fh]; for (int i=0; isetProgress (0.5); + // luminance histogram update memset (hist16, 0, 65536*sizeof(int)); for (int i=0; i +///////////////////////////////////////////////////////////////////////////// +// Haar wavelets // +///////////////////////////////////////////////////////////////////////////// + +// Bad, strong block effect + +template void dwt_haar(T * data, size_t pitch, T * buffer, size_t n) { size_t n2a = (n + 1) / 2; @@ -46,7 +52,7 @@ void dwt_haar(T * data, size_t pitch, T * buffer, size_t n) } } -template +template void idwt_haar(T * data, size_t pitch, T * buffer, size_t n, int alpha) { size_t n2a = (n + 1) / 2; @@ -71,9 +77,13 @@ void idwt_haar(T * data, size_t pitch, T * buffer, size_t n, int alpha) } } +///////////////////////////////////////////////////////////////////////////// +// CDF 5/3 wavelets // +///////////////////////////////////////////////////////////////////////////// + // buffer must be of length (n + 4) -template +template void dwt_53(T * data, size_t pitch, T * buffer, size_t n) { size_t n2 = n/2; @@ -97,14 +107,14 @@ void dwt_53(T * data, size_t pitch, T * buffer, size_t n) // calculate coefficients - for(int j = -1; j < (int)n + 1; j += 2) + for(ptrdiff_t i = -1; i < (ptrdiff_t)n + 1; i += 2) { - tmp[j] = tmp[j] - (tmp[j-1] + tmp[j+1]) / 2; + tmp[i] = tmp[i] - (tmp[i-1] + tmp[i+1]) / 2; } - for(int i = 0; i < (int)n; i += 2) + for(ptrdiff_t i = 0; i < (ptrdiff_t)n; i += 2) { - tmp[i] = tmp[i] + (tmp[i-1] + tmp[i+1] + 2) / 4; + tmp[i] = tmp[i] + (tmp[i-1] + tmp[i+1] + 2) / 4; } // copy with reordering @@ -120,7 +130,7 @@ void dwt_53(T * data, size_t pitch, T * buffer, size_t n) } } -template +template void idwt_53(T * data, size_t pitch, T * buffer, size_t n, int alpha) { size_t n2 = n/2; @@ -149,14 +159,14 @@ void idwt_53(T * data, size_t pitch, T * buffer, size_t n, int alpha) // calculate coefficients - for(int i = 0; i < (int)n + 1; i += 2) + for(ptrdiff_t i = 0; i < (ptrdiff_t)n + 1; i += 2) { - tmp[i] = tmp[i] - (tmp[i-1] + tmp[i+1] + 2) / 4; + tmp[i] = tmp[i] - (tmp[i-1] + tmp[i+1] + 2) / 4; } - for(int j = 1; j < (int)n; j += 2) + for(ptrdiff_t i = 1; i < (ptrdiff_t)n; i += 2) { - tmp[j] = tmp[j] + (tmp[j-1] + tmp[j+1]) / 2; + tmp[i] = tmp[i] + (tmp[i-1] + tmp[i+1]) / 2; } // copy data @@ -167,4 +177,84 @@ void idwt_53(T * data, size_t pitch, T * buffer, size_t n, int alpha) } } +///////////////////////////////////////////////////////////////////////////// +// Edge-avoiding wavelets // +///////////////////////////////////////////////////////////////////////////// +// based on +// Edge-Avoiding Wavelets and their Applications +// Raanan Fattal , +// Hebrew University of Jerusalem, Israel +// +// WCDF variant from this paper is used here +// T must be one of floating-point types +///////////////////////////////////////////////////////////////////////////// + +template +inline T wcdf_weight(T a, T b) +{ + static const T eps = 1e-5; + static const T one = 1.0; + return one / (abs(a - b) + eps); // will use overloaded abs for floating numbers +} + +template +void dwt_wcdf(T * data, size_t pitch, T * buffer, T * buffer2, size_t n) +{ + size_t n2 = n/2; + size_t n2a = (n + 1) / 2; + T * tmp = buffer + 2; + T * w = buffer2 + 2; + + // copy data + + for(size_t i = 0, j = 0; i < n; i++, j += pitch) + { + tmp[i] = data[j]; + } + + // extend mirror-like + + tmp[-1] = tmp[1]; + tmp[-2] = tmp[2]; + + tmp[n] = tmp[n-2]; + tmp[n+1] = tmp[n-3]; + + // calculate weights + + for(ptrdiff_t i = 0; i < (ptrdiff_t)n - 1; i++) + { + w[i] = wcdf_weight(tmp[i], tmp[i+1]); + } + + w[-1] = w[-2] = (T)0.0; + w[n] = w[n + 1] = (T)0.0; + + // calculate coefficients + + for(ptrdiff_t i = 1; i < (ptrdiff_t)n; i += 2) + { + tmp[i] = tmp[i] - (w[i-1]*tmp[i-1] + w[i]*tmp[i+1]) / (w[i-1] + w[i]); + } + + for(ptrdiff_t i = 0; i < (ptrdiff_t)n; i += 2) + { + tmp[i] = tmp[i] + (T)0.5 * (w[i-1]*tmp[i-1] + w[i]*tmp[i+1]) / (w[i-1] + w[i]); + } + + // copy with reordering + + for(size_t i = 0, j = 0; i < n; i+=2, j += pitch) + { + data[j] = tmp[i]; + } + + for(size_t i = 1, j = n2a*pitch; i < n; i+=2, j += pitch) + { + data[j] = tmp[i]; + } +} + + + #endif diff --git a/rtengine/wavelet_dec.cc b/rtengine/wavelet_dec.cc index a972d7c7c..5ab9e272a 100644 --- a/rtengine/wavelet_dec.cc +++ b/rtengine/wavelet_dec.cc @@ -22,22 +22,6 @@ namespace rtengine { -wavelet_decomposition::wavelet_decomposition(unsigned short ** src, size_t w, size_t h) -: nlevels(0), m_w(w), m_h(h), m_w1(0), m_h1(0) -{ - m_w1 = w; - m_h1 = h; - - m_c[0] = new wavelet_level(src, m_w1, m_h1); - nlevels = 1; - - while(nlevels < maxlevels) - { - m_c[nlevels] = new wavelet_level(m_c[nlevels - 1]->lowfreq(), m_c[nlevels-1]->lfw(), m_c[nlevels-1]->lfh()); - nlevels ++; - } -} - wavelet_decomposition::~wavelet_decomposition() { for(int i = 0; i < nlevels; i++) @@ -46,17 +30,5 @@ wavelet_decomposition::~wavelet_decomposition() } } -void wavelet_decomposition::reconstruct(unsigned short ** dst, const int * c) -{ - for(int level = nlevels - 1; level > 0; level--) - { - int alpha = 1024 + 10 * c[level]; - m_c[level]->reconstruct(m_c[level-1]->lowfreq(), alpha); - } - - int alpha = 1024 + 10 * c[0]; - m_c[0]->reconstruct(dst, alpha, wavelet_level::CLIP); -} - }; diff --git a/rtengine/wavelet_dec.h b/rtengine/wavelet_dec.h index 372849fec..c53b3a117 100644 --- a/rtengine/wavelet_dec.h +++ b/rtengine/wavelet_dec.h @@ -28,8 +28,12 @@ namespace rtengine { class wavelet_decomposition { +public: + typedef int internal_type; +private: + static const int maxlevels = 8; int nlevels; @@ -39,15 +43,52 @@ class wavelet_decomposition wavelet_level * m_c[maxlevels]; public: - wavelet_decomposition(unsigned short ** src, size_t w, size_t h); + + template + wavelet_decomposition(E ** src, size_t w, size_t h); ~wavelet_decomposition(); - void reconstruct(unsigned short ** dst, const int * c); + template + void reconstruct(E ** dst, const int * c, L & limiter); }; ////////////////////////////////////////////////////////////////////////////// +template +wavelet_decomposition::wavelet_decomposition(E ** src, size_t w, size_t h) +: nlevels(0), m_w(w), m_h(h), m_w1(0), m_h1(0) +{ + m_w1 = w; + m_h1 = h; + + m_c[0] = new wavelet_level(src, m_w1, m_h1); + nlevels = 1; + + while(nlevels < maxlevels) + { + m_c[nlevels] = new wavelet_level(m_c[nlevels - 1]->lowfreq(), m_c[nlevels-1]->lfw(), m_c[nlevels-1]->lfh()); + nlevels ++; + } +} + + +template +void wavelet_decomposition::reconstruct(E ** dst, const int * c, L & l) +{ + noop n; + + for(int level = nlevels - 1; level > 0; level--) + { + int alpha = 1024 + 10 * c[level]; + m_c[level]->reconstruct(m_c[level-1]->lowfreq(), alpha, n); + } + + int alpha = 1024 + 10 * c[0]; + m_c[0]->reconstruct(dst, alpha, l); +} + + }; #endif diff --git a/rtengine/wavelet_level.h b/rtengine/wavelet_level.h index ec950d85e..c70bbde1c 100644 --- a/rtengine/wavelet_level.h +++ b/rtengine/wavelet_level.h @@ -28,7 +28,36 @@ namespace rtengine { -template +template +class limiter +{ + T min_value, max_value; +public: + limiter(T min, T max) + : min_value(min), max_value(max) + {} + + T operator()(T x) + { + if(x < min_value) + return min_value; + if(x > max_value) + return max_value; + return x; + } +}; + +template +class noop +{ +public: + T operator()(T x) + { + return x; + } +}; + +template inline T clip(T x, T min_value, T max_value) { if(x < min_value) @@ -38,17 +67,17 @@ inline T clip(T x, T min_value, T max_value) return x; } -template -void plane_copy(A ** a, B ** b, size_t w, size_t h) +template +void plane_copy(A ** a, B ** b, size_t w, size_t h, L & l) { for(size_t j = 0; j < h; j++) for(size_t i = 0; i < w; i++) - b[j][i] = static_cast(a[j][i]); + b[j][i] = static_cast(l(a[j][i])); } ////////////////////////////////////////////////////////////////////////////// -template +template class wavelet_level { // full size @@ -73,9 +102,7 @@ class wavelet_level public: - static const bool CLIP = true; - - template + template wavelet_level(E ** src, size_t w, size_t h) : m_w(w), m_h(h), m_w2((w+1)/2), m_h2((h+1)/2), m_pitch(0), m_coeffs(NULL) { @@ -104,16 +131,16 @@ public: return m_h2; } - template + template void decompose(E ** src); - template - void reconstruct(E ** dst, int alpha, bool do_clip = false); + template + void reconstruct(E ** dst, int alpha, L & limiter); }; ////////////////////////////////////////////////////////////////////////////// -template +template void wavelet_level::dwt_2d(size_t w, size_t h) { T * buffer = new T[std::max(w, h) + 4]; @@ -133,7 +160,7 @@ void wavelet_level::dwt_2d(size_t w, size_t h) delete buffer; } -template +template void wavelet_level::idwt_2d(size_t w, size_t h, int alpha) { T * buffer = new T[std::max(w, h) + 4]; @@ -154,7 +181,7 @@ void wavelet_level::idwt_2d(size_t w, size_t h, int alpha) } -template +template void wavelet_level::create() { // 16-byte alignment: no effect @@ -168,7 +195,7 @@ void wavelet_level::create() } } -template +template void wavelet_level::destroy() { if(m_coeffs) @@ -179,28 +206,22 @@ void wavelet_level::destroy() } } -template template +template template void wavelet_level::decompose(E ** src) { - plane_copy(src, m_coeffs, m_w, m_h); + noop l; + + plane_copy(src, m_coeffs, m_w, m_h, l); dwt_2d(m_w, m_h); } -template template -void wavelet_level::reconstruct(E ** dst, int alpha, bool do_clip) +template template +void wavelet_level::reconstruct(E ** dst, int alpha, L & l) { idwt_2d(m_w, m_h, alpha); - if(do_clip) - { - for(size_t j = 0; j < m_h; j++) - for(size_t i = 0; i < m_w; i++) - m_coeffs[j][i] = clip(m_coeffs[j][i], 0, 65535); - } - - plane_copy(m_coeffs, dst, m_w, m_h); - + plane_copy(m_coeffs, dst, m_w, m_h, l); } };