As promised, here is an implementation of the Cohen-Daubechies-Feauveau 9 tap / 7 tap wavelet transform on a 2D signal in Python. This is the same transform used in the JPEG2000 codec.

''' 2D CDF 9/7 Wavelet Forward and Inverse Transform (lifting implementation) This code is provided "as is" and is given for educational purposes. 2008 - Kris Olhovsky - code.inquiries@olhovsky.com ''' from PIL import Image # Part of the standard Python Library ''' Example matrix as a list of lists: ''' mat4x4 = [ [0, 1, 2, 3], # Row 1 [4, 5, 6, 7], # Row 2 [8, 9, 10, 11], # Row 3 [12, 13, 14, 15], # Row 4 ] # We don't do anything with this matrix. # It's just here for clarification. def fwt97_2d(m, nlevels=1): ''' Perform the CDF 9/7 transform on a 2D matrix signal m. nlevel is the desired number of times to recursively transform the signal. ''' w = len(m[0]) h = len(m) for i in range(nlevels): m = fwt97(m, w, h) # cols m = fwt97(m, w, h) # rows w /= 2 h /= 2 return m def iwt97_2d(m, nlevels=1): ''' Inverse CDF 9/7 transform on a 2D matrix signal m. nlevels must be the same as the nlevels used to perform the fwt. ''' w = len(m[0]) h = len(m) # Find starting size of m: for i in range(nlevels-1): h /= 2 w /= 2 for i in range(nlevels): m = iwt97(m, w, h) # rows m = iwt97(m, w, h) # cols h *= 2 w *= 2 return m def fwt97(s, width, height): ''' Forward Cohen-Daubechies-Feauveau 9 tap / 7 tap wavelet transform performed on all columns of the 2D n*n matrix signal s via lifting. The returned result is s, the modified input matrix. The highpass and lowpass results are stored on the left half and right half of s respectively, after the matrix is transposed. ''' # 9/7 Coefficients: a1 = -1.586134342 a2 = -0.05298011854 a3 = 0.8829110762 a4 = 0.4435068522 # Scale coeff: k1 = 0.81289306611596146 # 1/1.230174104914 k2 = 0.61508705245700002 # 1.230174104914/2 # Another k used by P. Getreuer is 1.1496043988602418 for col in range(width): # Do the 1D transform on all cols: ''' Core 1D lifting process in this loop. ''' ''' Lifting is done on the cols. ''' # Predict 1. y1 for row in range(1, height-1, 2): s[row][col] += a1 * (s[row-1][col] + s[row+1][col]) s[height-1][col] += 2 * a1 * s[height-2][col] # Symmetric extension # Update 1. y0 for row in range(2, height, 2): s[row][col] += a2 * (s[row-1][col] + s[row+1][col]) s[0][col] += 2 * a2 * s[1][col] # Symmetric extension # Predict 2. for row in range(1, height-1, 2): s[row][col] += a3 * (s[row-1][col] + s[row+1][col]) s[height-1][col] += 2 * a3 * s[height-2][col] # Update 2. for row in range(2, height, 2): s[row][col] += a4 * (s[row-1][col] + s[row+1][col]) s[0][col] += 2 * a4 * s[1][col] # de-interleave temp_bank = [[0]*width for i in range(height)] for row in range(height): for col in range(width): # k1 and k2 scale the vals # simultaneously transpose the matrix when deinterleaving if row % 2 == 0: # even temp_bank[col][row/2] = k1 * s[row][col] else: # odd temp_bank[col][row/2 + height/2] = k2 * s[row][col] # write temp_bank to s: for row in range(width): for col in range(height): s[row][col] = temp_bank[row][col] return s def iwt97(s, width, height): ''' Inverse CDF 9/7. ''' # 9/7 inverse coefficients: a1 = 1.586134342 a2 = 0.05298011854 a3 = -0.8829110762 a4 = -0.4435068522 # Inverse scale coeffs: k1 = 1.230174104914 k2 = 1.6257861322319229 # Interleave: temp_bank = [[0]*width for i in range(height)] for col in range(width/2): for row in range(height): # k1 and k2 scale the vals # simultaneously transpose the matrix when interleaving temp_bank[col * 2][row] = k1 * s[row][col] temp_bank[col * 2 + 1][row] = k2 * s[row][col + width/2] # write temp_bank to s: for row in range(width): for col in range(height): s[row][col] = temp_bank[row][col] for col in range(width): # Do the 1D transform on all cols: ''' Perform the inverse 1D transform. ''' # Inverse update 2. for row in range(2, height, 2): s[row][col] += a4 * (s[row-1][col] + s[row+1][col]) s[0][col] += 2 * a4 * s[1][col] # Inverse predict 2. for row in range(1, height-1, 2): s[row][col] += a3 * (s[row-1][col] + s[row+1][col]) s[height-1][col] += 2 * a3 * s[height-2][col] # Inverse update 1. for row in range(2, height, 2): s[row][col] += a2 * (s[row-1][col] + s[row+1][col]) s[0][col] += 2 * a2 * s[1][col] # Symmetric extension # Inverse predict 1. for row in range(1, height-1, 2): s[row][col] += a1 * (s[row-1][col] + s[row+1][col]) s[height-1][col] += 2 * a1 * s[height-2][col] # Symmetric extension return s def seq_to_img(m, pix): ''' Copy matrix m to pixel buffer pix. Assumes m has the same number of rows and cols as pix. ''' for row in range(len(m)): for col in range(len(m[row])): pix[col,row] = m[row][col]

And here is an example usage:

if __name__ == "__main__": # Load image. im = Image.open("test1_512.png") # Must be a single band image! (grey) # Create an image buffer object for fast access. pix = im.load() # Convert the 2d image to a 1d sequence: m = list(im.getdata()) # Convert the 1d sequence to a 2d matrix. # Each sublist represents a row. Access is done via m[row][col]. m = [m[i:i+im.size[0]] for i in range(0, len(m), im.size[0])] # Cast every item in the list to a float: for row in range(0, len(m)): for col in range(0, len(m[0])): m[row][col] = float(m[row][col]) # Perform a forward CDF 9/7 transform on the image: m = fwt97_2d(m, 3) seq_to_img(m, pix) # Convert the list of lists matrix to an image. im.save("test1_512_fwt.png") # Save the transformed image. # Perform an inverse transform: m = iwt97_2d(m, 3) seq_to_img(m, pix) # Convert the inverse list of lists matrix to an image. im.save("test1_512_iwt.png") # Save the inverse transformation.

This is a first step in a larger task of compressing a certain type of vertex data, which can be decompressed entirely by a GPU via an implementation written in shader code.

The reason for a Python implementation is to create a working, easy to understand version of the implementation. Not only for myself, but for the interwebs.

The Python code does use “height” and “width” variables, although in order to use rectangular images, some adjustments need to be made. This code also assumes that your images have width and height of the form 2^n. E.g. 1024 x 1024, 512 x 512, etc.

Hi !

I am interested in the implementation of CDF 9/7 in RealBasic to use it for denoising pictures in my softwar SeamlessMaker (http://www.hypatiasoft.fr). I have already writen all subroutines it for Haar’s wavelets but the results are not very good. Can I use the same subroutines for CDF 9/7 just changing the High pass and lowpass filters ? Can I trust in these coefficients : http://en.wikipedia.org/wiki/Cohen-Daubechies-Feauveau_wavelet ?

Thanks in advance.

Géraud Bousquet

You can trust the coefficients on wikipedia at the time of posting this. The coefficients I used in my python code are also correct, and mine have more precision (which won’t matter unless you use 64 bit floats).

Hi,

Could you please add an explicit license for this code, or declare that you release it to the public domain?

Thank you,

Chris

I’ve edited the code as per your suggestion, thanks Chris.

Hi,

Does this algorithm require grey scale images? If I did the same mathematics on each channel would I be able to make this work with colour images?

Thanks,

-Matt

You could simply treat a color image as 3 separate images, one for each channel.

Then, perform the transform on each channel, and undo the transform on each channel later.

Hi Olhovsky,

From the previous comment by Chris, it isn’t yet very clear which license applies to this code. Are you attaching any particular license to this source code? Can I use it under a LGPL application?

Best regards,

César

Hi,

I am testing your code for a student project an I’m receiving this error:

Traceback (most recent call last):

File “wavelet_02.py”, line 200, in

m = fwt97_2d(m, 3)

File “wavelet_02.py”, line 27, in fwt97_2d

m = fwt97(m, w, h) # cols

File “wavelet_02.py”, line 108, in fwt97

temp_bank[col][row/2 + height/2] = k2 * s[row][col]

IndexError: list assignment index out of range

Mac OS 10.7.3

Python 2.7.3

PIL 1.1.7

Maybe you have seen this error before?

Cheers,

Tobi

I tried implementing this code in Vb.Net for 3-Level Multi resolution on Y component of Color image. But PSNR were around 18-22, whereas for Cb, Cr component range is 24-26. I had used 512×512 color image-set.

Hence i Switched over to Matlab with different wavelet.

But my query is: how JPEG manages to be visually lossless with this wavelet?

Hi,

I wonder how can I utilize your symmetric extension technique in FPGA.

Can I use replicate the first pixel as the previous sample and the last pixel as the following pixel?