4
$\begingroup$

I try to understand what following algorithm from this gimp plugin does to denoise an image:

 /* Wavelet denoise GIMP plugin
 * 
 * wavelet.c
 * Copyright 2008 by Marco Rossini
 * 
 * Implements the wavelet denoise code of UFRaw by Udi Fuchs
 * which itself bases on the code by Dave Coffin
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License version 2
 * as published by the Free Software Foundation.
 * 
 * Instructions:
 * compile with gimptool, eg. 'gimptool-2.0 --install wavelet-denoise.c'
 */

void
wavelet_denoise (float *fimg[3], unsigned int width,
         unsigned int height, float threshold, double low, float a,
         float b)
{
  float *temp, thold;
  unsigned int i, lev, lpass, hpass, size, col, row;
  double stdev[5];
  unsigned int samples[5];

  size = width * height;

  /* FIXME: replace by GIMP functions */
  temp = (float *) malloc (MAX2 (width, height) * sizeof (float));

  hpass = 0;
  for (lev = 0; lev < 5; lev++)
    {
      if (b != 0)
    gimp_progress_update (a + b * lev / 5.0);
      lpass = ((lev & 1) + 1);
      for (row = 0; row < height; row++)
    {
      hat_transform (temp, fimg[hpass] + row * width, 1, width, 1 << lev);
      for (col = 0; col < width; col++)
        {
          fimg[lpass][row * width + col] = temp[col] * 0.25;
        }
    }
      if (b != 0)
    gimp_progress_update (a + b * (lev + 0.25) / 5.0);
      for (col = 0; col < width; col++)
    {
      hat_transform (temp, fimg[lpass] + col, width, height, 1 << lev);
      for (row = 0; row < height; row++)
        {
          fimg[lpass][row * width + col] = temp[row] * 0.25;
        }
    }
      if (b != 0)
    gimp_progress_update (a + b * (lev + 0.5) / 5.0);

      thold =
    5.0 / (1 << 6) * exp (-2.6 * sqrt (lev + 1)) * 0.8002 / exp (-2.6);

      /* initialize stdev values for all intensities */
      stdev[0] = stdev[1] = stdev[2] = stdev[3] = stdev[4] = 0.0;
      samples[0] = samples[1] = samples[2] = samples[3] = samples[4] = 0;

      /* calculate stdevs for all intensities */
      for (i = 0; i < size; i++)
    {
      fimg[hpass][i] -= fimg[lpass][i];
      if (fimg[hpass][i] < thold && fimg[hpass][i] > -thold)
        {
          if (fimg[lpass][i] > 0.8) {
            stdev[4] += fimg[hpass][i] * fimg[hpass][i];
            samples[4]++;
          } else if (fimg[lpass][i] > 0.6) {
            stdev[3] += fimg[hpass][i] * fimg[hpass][i];
            samples[3]++;
          } else if (fimg[lpass][i] > 0.4) {
            stdev[2] += fimg[hpass][i] * fimg[hpass][i];
            samples[2]++;
          } else if (fimg[lpass][i] > 0.2) {
            stdev[1] += fimg[hpass][i] * fimg[hpass][i];
            samples[1]++;
          } else {
            stdev[0] += fimg[hpass][i] * fimg[hpass][i];
            samples[0]++;
          }
        }
    }
      stdev[0] = sqrt (stdev[0] / (samples[0] + 1));
      stdev[1] = sqrt (stdev[1] / (samples[1] + 1));
      stdev[2] = sqrt (stdev[2] / (samples[2] + 1));
      stdev[3] = sqrt (stdev[3] / (samples[3] + 1));
      stdev[4] = sqrt (stdev[4] / (samples[4] + 1));

      if (b != 0)
    gimp_progress_update (a + b * (lev + 0.75) / 5.0);

      /* do thresholding */
      for (i = 0; i < size; i++)
    {
      if (fimg[lpass][i] > 0.8) {
        thold = threshold * stdev[4];
      } else if (fimg[lpass][i] > 0.6) {
        thold = threshold * stdev[3];
      } else if (fimg[lpass][i] > 0.4) {
        thold = threshold * stdev[2];
      } else if (fimg[lpass][i] > 0.2) {
        thold = threshold * stdev[1];
      } else {
        thold = threshold * stdev[0];
      }

      if (fimg[hpass][i] < -thold)
        fimg[hpass][i] += thold - thold * low;
      else if (fimg[hpass][i] > thold)
        fimg[hpass][i] -= thold - thold * low;
      else
        fimg[hpass][i] *= low;

      if (hpass)
        fimg[0][i] += fimg[hpass][i];
    }
      hpass = lpass;
    }

  for (i = 0; i < size; i++)
    fimg[0][i] = fimg[0][i] + fimg[lpass][i];

  /* FIXME: replace by GIMP functions */
  free (temp);
}

Code of UFRaw by Udi Fuchs, which itself bases on the code by Dave Coffin

Why are wavelets key to do that? How are the pixels calculated?

Any help would be greatly appreciated!

$\endgroup$
  • $\begingroup$ Very interesting... $\endgroup$ – Royi Jan 18 '15 at 17:44
  • $\begingroup$ @lennon310 Thanks for the edit :) Have added the denoising tag too. $\endgroup$ – p2or Jan 18 '15 at 19:57
  • 1
    $\begingroup$ Looks like it is using wavelet thresholding to denoise the image. cs.haifa.ac.il/hagit/courses/seminars/wavelets/Presentations/… $\endgroup$ – Aaron Jan 28 '15 at 20:01
  • 1
    $\begingroup$ wavelet is to study nonstationary signals, for nonstationary signals, you often have nonstationary noise so yo can't just use a FIR filter. Wavelet thresholding method such as SURE estimates the noisy wavelet coefficients and eliminates them from the total wavelet coefficients, then transforms the signal back to time domain $\endgroup$ – Aåkon Jan 30 '15 at 2:50
  • $\begingroup$ Thanks @IllegalImmigrant! Sounds like you are very familiar with it. :) Why don't you write an answer? Thanks again for sharing your knowledge. $\endgroup$ – p2or Feb 1 '15 at 14:41
4
+50
$\begingroup$

Wavelets are not key to denoising. There are different ways to denoise an image, for example in the original signal domain or in the transform domain (i.e. Fourier or wavelet). Wavelets work best for additive noise, where the noise is random & not correlated in time.

wavelet_denoise (float *fimg[3], unsigned int width,
     unsigned int height, float threshold, double low, float a, float b)

In the above function, the image is passed in as fimg[3], representing the RGB channels of the image. The image is decomposed, thresholded & re-written back into fimg.

More specifically, the wavelet transform is a fast hierarchical scheme for processing an image using a set of consecutive lowpass and highpass filters, followed by a decimation. This results in a decomposition into different scales which can be regarded as different “frequency bands”, determined by the mother wavelet, see:

hat_transform (temp, fimg[lpass] + col, width, height, 1 << lev);

Lowpass filter outputs are called approximation coefficients & highpass filters are referred to as detail coefficients. The detail coefficients have high values in the noisy parts of the signal. To threshold, each coefficient is compared to a threshold value and attenuated/shrunk by some factor. Finally, the image (fimg) is reconstructed from the thresholded wavelet coefficients.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.