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 internets.
Imagine sending 1/10th the data when world geometry information needs to change. This would use far less CPU and GPU time transferring data to the GPU, as well as freeing up CPU time that would normally be used decompressing the data, in situations where HDD space is the reason for the compression.
HDD reads will be shorter, space used in the GPU RAM will be saved.
Coming soon; a C or C# implementation of the CDF 9/7 transform, along with a quantizer and encoder, to realize a GPU based compression scheme.
The goals of my compression scheme are to reduce both compression and decompression time (while balancing compression ratio). Every other scheme I’ve seen on the internets don’t seem to focus on reducing the compression time. As you’ll soon see, this is because my larger project requires fast compression (to allow quick alterations of compressed vertex data).
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.

Sample input image for transforming.

Image after a single level CDF 9/7 transform.

2 level CDF 9/7 wavelet transformed image
#1 by GERAUD BOUSQUET on -
Quote
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
#2 by olhovsky on -
Quote
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).
#3 by Chris Murphy on -
Quote
Hi,
Could you please add an explicit license for this code, or declare that you release it to the public domain?
Thank you,
Chris
#4 by olhovsky on -
Quote
I’ve edited the code as per your suggestion, thanks Chris.