Wavelet equalizer: changed to edge-avoiding wavelets (as of Fattal)

rtengine/wavelet.h

@@ -192,13 +192,16 @@ void idwt_53(T * data, size_t pitch, T * buffer, size_t n, int alpha)
 template<typename T>
 inline T wcdf_weight(T a, T b)
 {
-    static const T eps = 1e-5;
+    static const T eps = 1;
     static const T one = 1.0;
-    return one / (abs(a - b) + eps); // will use overloaded abs for floating numbers
+    return one / (fabs(a - b) + eps);
 }
 
+// buffer is a temporary storage
+// buffer2 must be preserved between dwt and idwt
+
 template<typename T>
-void dwt_wcdf(T * data, size_t pitch, T * buffer, T * buffer2, size_t n)
+void dwt_wcdf(T * data, size_t pitch, T * buffer, size_t n, T * buffer2)
 {
     size_t n2 = n/2;
     size_t n2a = (n + 1) / 2;
@@ -225,10 +228,11 @@ void dwt_wcdf(T * data, size_t pitch, T * buffer, T * buffer2, size_t n)
     for(ptrdiff_t i = 0; i < (ptrdiff_t)n - 1; i++)
     {
         w[i] = wcdf_weight(tmp[i], tmp[i+1]);
+        //w[i] = 1;
     }
     
     w[-1] = w[-2] = (T)0.0;
-    w[n] = w[n + 1] = (T)0.0;
+    w[n-1] = w[n] = w[n + 1] = (T)0.0;
 
     // calculate coefficients
     
@@ -255,6 +259,54 @@ void dwt_wcdf(T * data, size_t pitch, T * buffer, T * buffer2, size_t n)
     }
 }
 
+template<typename T>
+void idwt_wcdf(T * data, size_t pitch, T * buffer, size_t n, int alpha, T * buffer2)
+{
+    size_t n2 = n/2;
+    size_t n2a = (n + 1) / 2;
+    T * tmp = buffer + 2;
+    T * w = buffer2 + 2;
 
+    // copy with reordering
+    
+    for(size_t i = 0, j = 0; i < n; i+=2, j += pitch)
+    {
+        tmp[i] = data[j];
+    }
+
+    for(size_t i = 1, j = n2a*pitch; i < n; i+=2, j += pitch)
+    {
+        tmp[i] = (alpha * data[j]) / 1024;
+    }
+
+    // extend mirror-like
+    
+    tmp[-1] = tmp[1];
+    tmp[-2] = tmp[2];
+    
+    tmp[n] = tmp[n-2];
+    tmp[n+1] = tmp[n-3];
+
+    // calculate coefficients
+
+    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]);
+    }
+
+    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]);
+    }
+   
+    // copy data
+
+    for(size_t i = 0, j = 0; i < n; i++, j += pitch)
+    {
+        data[j] = tmp[i];
+    }
+}
+
+//////////////////////////////////////////////////////////////////////////////
 
 #endif

rtengine/wavelet_dec.h

@@ -30,7 +30,7 @@ class wavelet_decomposition
 {
 public:
 
-    typedef int internal_type;
+    typedef float internal_type;
 
 private:
 

rtengine/wavelet_level.h

@@ -86,16 +86,17 @@ class wavelet_level
     // size of low frequency part
     size_t m_w2, m_h2;
 
-    // distance between lines in the array of coeffs
-    size_t m_pitch;
-
     // array of pointers to lines of coeffs
     // actually is a single contiguous data array pointed by m_coeffs[0]
     T ** m_coeffs;
-
+    
+    // weights storage
+    T ** m_weights_rows;
+    T ** m_weights_cols;
+    
     // allocation and destruction of data storage
-    void create();
+    T ** create(size_t w, size_t h);
-    void destroy();
+    void destroy(T ** p);
 
     void dwt_2d(size_t w, size_t h);
     void idwt_2d(size_t w, size_t h, int alpha);
@@ -104,16 +105,21 @@ public:
 
     template<typename 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)
+    : m_w(w), m_h(h), m_w2((w+1)/2), m_h2((h+1)/2), 
+      m_coeffs(NULL), m_weights_rows(NULL), m_weights_cols(NULL)
     {
-        create();
+        m_coeffs = create(w, h);
+        m_weights_rows = create(w + 4, h);
+        m_weights_cols = create(h + 4, w);
         
         decompose(src);
     }
     
     ~wavelet_level()
     {
-        destroy();
+        destroy(m_coeffs);
+        destroy(m_weights_rows);
+        destroy(m_weights_cols);
     }
     
     T ** lowfreq() const
@@ -148,13 +154,15 @@ void wavelet_level<T>::dwt_2d(size_t w, size_t h)
     for(size_t j = 0; j < h; j++)
     {
         //dwt_haar(m_coeffs[j], 1, buffer, w);
-        dwt_53(m_coeffs[j], 1, buffer, w);
+        //dwt_53(m_coeffs[j], 1, buffer, w);
+        dwt_wcdf(m_coeffs[j], 1, buffer, w, m_weights_rows[j]);
     }
     
     for(size_t i = 0; i < w; i++)
     {
         //dwt_haar(&m_coeffs[0][i], m_pitch, buffer, h);
-        dwt_53(&m_coeffs[0][i], m_pitch, buffer, h);
+        //dwt_53(&m_coeffs[0][i], w, buffer, h);
+        dwt_wcdf(&m_coeffs[0][i], w, buffer, h, m_weights_cols[i]);
     }
     
     delete[] buffer;
@@ -165,44 +173,45 @@ void wavelet_level<T>::idwt_2d(size_t w, size_t h, int alpha)
 {
     T * buffer = new T[std::max(w, h) + 4];
     
-    for(size_t j = 0; j < h; j++)
-    {
-        //idwt_haar(m_coeffs[j], 1, buffer, w, alpha);
-        idwt_53(m_coeffs[j], 1, buffer, w, alpha);
-    }
-    
     for(size_t i = 0; i < w; i++)
     {
         //idwt_haar(&m_coeffs[0][i], m_pitch, buffer, h, alpha);
-        idwt_53(&m_coeffs[0][i], m_pitch, buffer, h, alpha);
+        //idwt_53(&m_coeffs[0][i], w, buffer, h, alpha);
+        idwt_wcdf(&m_coeffs[0][i], w, buffer, h, alpha, m_weights_cols[i]);
+        //idwt_noop(&m_coeffs[0][i], w, buffer, h, alpha);
     }
     
+    for(size_t j = 0; j < h; j++)
+    {
+        //idwt_haar(m_coeffs[j], 1, buffer, w, alpha);
+        //idwt_53(m_coeffs[j], 1, buffer, w, alpha);
+        idwt_wcdf(m_coeffs[j], 1, buffer, w, alpha, m_weights_rows[j]);
+        //idwt_noop(m_coeffs[j], 1, buffer, w, alpha);
+    }
+
     delete[] buffer;
-
 }
 
 template<typename T>
-void wavelet_level<T>::create()
+T ** wavelet_level<T>::create(size_t w, size_t h)
 {
-    // 16-byte alignment: no effect
+    T * data = new T[w * h];
-    //m_pitch = (((m_w*sizeof(T) + 15) / 16) * 16) / sizeof(T);
+    T ** p = new T*[h];
-    m_pitch = m_w;
+    for(size_t j = 0; j < h; j++)
-    T * data = new T[m_pitch * m_h];
-    m_coeffs = new T*[m_h];
-    for(size_t j = 0; j < m_h; j++)
     {
-        m_coeffs[j] = data + m_pitch * j;
+        p[j] = data + w * j;
     }
+    return p;
 }
 
 template<typename T>
-void wavelet_level<T>::destroy()
+void wavelet_level<T>::destroy(T ** p)
 {
-    if(m_coeffs)
+    if(p)
     {
-        delete[] m_coeffs[0];
+        delete[] p[0];
     
-        delete[] m_coeffs;
+        delete[] p;
     }
 }