2D CDF 9/7 Wavelet Transform in Python

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.

Sample input image for transforming.

Sample input image for transforming.

Image after a single level CDF 9/7 transform.

Image after a single level CDF 9/7 transform.

2 level CDF 9/7 wavelet transformed image

2 level CDF 9/7 wavelet transformed image

This entry was posted in Coding and tagged , , , , , , . Bookmark the permalink.

10 Responses to 2D CDF 9/7 Wavelet Transform in Python

  1. 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. olhovsky says:

    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. Chris Murphy says:

    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. olhovsky says:

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

  5. 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

  6. olhovsky says:

    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.

  7. Cesar says:

    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

  8. Tobias says:

    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

  9. Krish Vikram S says:

    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?

  10. Vinh says:

    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?

Leave a Reply

Your email address will not be published. Required fields are marked *

*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>