"Normalizing" my D4 wavelet transformation at each step reduces final image quality

Question

Original image:

(Images included are .png images, so no additional distortion was added in save/upload for viewing)

I have used the D4 transform from page 20 of "Ripples in mathematics", which is basically these 5 steps:

Forward d4:

c1 = √3 / 4.0 ;
c2 = (√3 - 2) / 4.0 ;
s[ IEVEN ] += √3 * s[ IODD ] ;
s[ IODD ] -= c1*s[ IEVEN ] + c2*s[ IPREVEVEN ] ;
s[ IEVEN ] -= s[ INEXTODD ] ;
s[ IEVEN ] *= ( √3 - 1 ) / √2 ;
s[ IODD ] *= ( √3 + 1 ) / √2 ;

The inverse:

c1 = √3 / 4.0 ;
c2 = (√3 - 2) / 4.0 ;
s[ IODD ] *= ( √3 - 1 ) / √2 ;
s[ IEVEN ] *= ( √3 + 1 ) / √2  ;
s[ IEVEN ] += s[INEXTODD] ;
s[ IODD ] += c1*s[ IEVEN ] + c2*s[IPREVEVEN] ;
s[ IEVEN ] -= √3 * s[ IODD ] ;

I'm compiling and running this using double precision values from C++. I run this on the rows of the image, then the columns. I use a crude filtration algorithm to remove the lowest 90% of the difference coefficients in the image.

The filtration algorithm is:

• Run through entire transformed image (as a set of numbers)
• Find the largest difference coefficient (maxVal) (in entire 2d image)
• Choose minValToSurvive as 1% of maxVal.
• If a difference coefficient has a magnitude less than minValToSurvive, it is zeroed.

Here's my problem is. When I remove only 83% of the lowest difference coefficients from the image (minValToSurvive=0.01*maxVal), you get this:

normalized

If I remove the normalization steps:

s[ IEVEN ] *= ( √3 - 1 ) / √2 ; // REMOVE
s[ IODD ] *= ( √3 + 1 ) / √2 ;

(in both the fwd and reverse transforms), the result after removing 90% of the components is much better (much less noise)

So I can think of 1 of 2 problems:

• Normalizing the image by the ( √3 - 1 ) / √2 factors is killing precision
• I'm not filtering correctly

Or am I wrong? If I'm filtering (removing insignificant components) incorrectly, what is a better way to filter? If it's the floating point precision, then should I not normalize the transform at every step?

Answer 1

The right answer is you have to perform each of the Update/Predict steps on the input signal completely before doing the next Update/Predict. What I was doing was walking through the signal, and performing each Update/Predict as I go.

On page 158 of "Ripples", there is a reference implementation.

// s is the signal
#define IEVEN (2*j)
#define IODD (2*j + i)
for( int i = 1 ; i <= n/2 ; i *= 2 )
{
  for( int j = 0 ; j <= n/2 - i ; j += i ) // Must do this Predict step COMPLETLEY
    s[ IEVEN ] += √3 * s[ IODD ] ;

  for( int j = 0 ; j <= n/2 - i ; j += i ) // Then this one..
  {
    int prevEvenIndex = IPREVEVEN ;
    s[ IODD ] -= d4c1*s[ IEVEN ] + d4c2*SAFE_PREV(s,prevEvenIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    int nextOddIndex = INEXTODD ;
    s[ IEVEN ] -= SAFE_NEXT(s,nextOddIndex) ;
  }

  for( int j = 0 ; j <= n/2 - i ; j += i )
  {
    s[ IEVEN ] *= d4normEvens ;
    s[ IODD ] *= d4normOdds ;
  }
}

The 98% 0's D4 transform: