Wavelet equalizer inprovemts and fixes

-now works in Lab space (on all three components)
-now works in preview (but only at 100% scale, as sharpening and denoise)
-several small fixes and refactorings

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);

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++;
     }
 

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);

rtengine/ipequalizer.cc

@@ -22,18 +22,55 @@
 
 #include <wavelet_dec.h>
 
+#include <iostream>
+
 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);
+    if (!params->equalizer.enabled) {
-    r.reconstruct(image->r, params.c);
+        return;
+    }
 
-    wavelet_decomposition g(image->g, fw, fh);
+    limiter<wavelet_decomposition::internal_type> l(0, 65535);
-    g.reconstruct(image->g, params.c);
 
-    wavelet_decomposition b(image->b, fw, fh);
+    wavelet_decomposition r(image->r, image->getW(), image->getH());
-    b.reconstruct(image->b, params.c);
+    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<wavelet_decomposition::internal_type> l1(0, 65535);
+
+        wavelet_decomposition L(image->L, image->W, image->H);
+        L.reconstruct(image->L, params->equalizer.c, l1);
+    }
+
+    if (chromaticity) {
+        limiter<wavelet_decomposition::internal_type> 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;
 
 }
 

rtengine/refreshmap.cc

@@ -103,6 +103,7 @@ ALL,    //    EvResizeEnabled
 ALL,    //    EvProfileChangeNotification
 RETINEX,    //    EvShrHighQuality
 TRANSFORM,  //    EvPerspCorr
-ALL     // ????? EvEqualizer
+EQUALIZER,  // EvEqualizer
+EQUALIZER   // EvEqlEnabled
 };
 

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

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; i<fh; i++)
@@ -143,6 +139,7 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
     if (pl) 
         pl->setProgress (0.5);
 
+
     // luminance histogram update
     memset (hist16, 0, 65536*sizeof(int));
     for (int i=0; i<fh; i++)
@@ -162,6 +159,9 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
     ipf.colorCurve (labView, labView);
     ipf.colordenoise (labView, buffer);
 
+    // wavelet equalizer
+    ipf.waveletEqualizer (labView, true, true);
+
     for (int i=0; i<fh; i++)
         delete [] buffer[i];
     delete [] buffer;

rtengine/wavelet.h

@@ -21,7 +21,13 @@
 #ifndef WAVELET_H_INCLUDED
 #define WAVELET_H_INCLUDED
 
-template<class T>
+/////////////////////////////////////////////////////////////////////////////
+//                        Haar wavelets                                    //
+/////////////////////////////////////////////////////////////////////////////
+
+// Bad, strong block effect
+
+template<typename T>
 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<class T>
+template<typename T>
 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<class T>
+template<typename T>
 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<class T>
+template<typename T>
 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 <raananf@cs.huji.ac.il>,
+//     Hebrew University of Jerusalem, Israel
+//
+// WCDF variant from this paper is used here
+// T must be one of floating-point types 
+/////////////////////////////////////////////////////////////////////////////
+
+template<typename T>
+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<typename T>
+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

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<internal_type>(src, m_w1, m_h1);
-    nlevels = 1;
-    
-    while(nlevels < maxlevels)
-    {
-        m_c[nlevels] = new wavelet_level<internal_type>(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<internal_type>::CLIP);
-}
-
 };
 

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<internal_type> * m_c[maxlevels];
     
 public:
-    wavelet_decomposition(unsigned short ** src, size_t w, size_t h);
+
+    template<typename E>
+    wavelet_decomposition(E ** src, size_t w, size_t h);
     
     ~wavelet_decomposition();
     
-    void reconstruct(unsigned short ** dst, const int * c);
+    template<typename E, typename L>
+    void reconstruct(E ** dst, const int * c, L & limiter);
 };
 
 //////////////////////////////////////////////////////////////////////////////
 
+template<typename E>
+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<internal_type>(src, m_w1, m_h1);
+    nlevels = 1;
+    
+    while(nlevels < maxlevels)
+    {
+        m_c[nlevels] = new wavelet_level<internal_type>(m_c[nlevels - 1]->lowfreq(), m_c[nlevels-1]->lfw(), m_c[nlevels-1]->lfh());
+        nlevels ++;
+    }
+}
+
+
+template<typename E, typename L>
+void wavelet_decomposition::reconstruct(E ** dst, const int * c, L & l)
+{
+    noop<internal_type> 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

rtengine/wavelet_level.h

@@ -28,7 +28,36 @@
 namespace rtengine
 {
 
-template<class T>
+template<typename T>
+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<typename T>
+class noop
+{
+public:
+    T operator()(T x)
+    {
+        return x;
+    }
+};
+
+template<typename T>
 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 <class A, class B>
+template <typename A, typename B, typename L>
-void plane_copy(A ** a, B ** b, size_t w, size_t h)
+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<B>(a[j][i]);
+            b[j][i] = static_cast<B>(l(a[j][i]));
 }
 
 //////////////////////////////////////////////////////////////////////////////
 
-template<class T>
+template<typename T>
 class wavelet_level
 {
     // full size
@@ -73,9 +102,7 @@ class wavelet_level
 
 public:
 
-    static const bool CLIP = true;
+    template<typename E>
-    
-    template<class E>
     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<class E>
+    template<typename E>
     void decompose(E ** src);
 
-    template<class E>
+    template<typename E, typename L>
-    void reconstruct(E ** dst, int alpha, bool do_clip = false);
+    void reconstruct(E ** dst, int alpha, L & limiter);
 };
 
 //////////////////////////////////////////////////////////////////////////////
 
-template<class T>
+template<typename T>
 void wavelet_level<T>::dwt_2d(size_t w, size_t h)
 {
     T * buffer = new T[std::max(w, h) + 4];
@@ -133,7 +160,7 @@ void wavelet_level<T>::dwt_2d(size_t w, size_t h)
     delete buffer;
 }
 
-template<class T>
+template<typename T>
 void wavelet_level<T>::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<T>::idwt_2d(size_t w, size_t h, int alpha)
 
 }
 
-template<class T>
+template<typename T>
 void wavelet_level<T>::create()
 {
     // 16-byte alignment: no effect
@@ -168,7 +195,7 @@ void wavelet_level<T>::create()
     }
 }
 
-template<class T>
+template<typename T>
 void wavelet_level<T>::destroy()
 {
     if(m_coeffs)
@@ -179,28 +206,22 @@ void wavelet_level<T>::destroy()
     }
 }
 
-template<class T> template<class E>
+template<typename T> template<typename E>
 void wavelet_level<T>::decompose(E ** src)
 {
-    plane_copy(src, m_coeffs, m_w, m_h);
+    noop<T> l;
+
+    plane_copy(src, m_coeffs, m_w, m_h, l);
 
     dwt_2d(m_w, m_h);
 }
 
-template<class T> template<class E>
+template<typename T> template<typename E, typename L>
-void wavelet_level<T>::reconstruct(E ** dst, int alpha, bool do_clip)
+void wavelet_level<T>::reconstruct(E ** dst, int alpha, L & l)
 {
     idwt_2d(m_w, m_h, alpha);
 
-    if(do_clip)
+    plane_copy(m_coeffs, dst, m_w, m_h, l);
-    {
-        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);
-
 }
 
 };