Haar Discrete Wavelet Transform

This was a project for an image processing class, and is useful for learning the method for lossless compression of JPEG-2000 images. Other source code examples exist online, but this one might be easier than others to understand. It is written in Python as well.


# Author: Ethan Webster
# Date: 3/12/2018
# Purpose: Performs forward and reverse Haar Discrete Wavelet Transform
# on an image

from skimage import io
import numpy as np

#constants for layer transformation formula
s0 = 0.5
s1 = 0.5
w0 = 0.5
w1 = -0.5

#Discrete 1D Haar Wavelet Transform
def Discrete1DWaveletTransform(arr):

    #declare empty array for current row/column
    tArr = np.empty(arr.size, np.float64)

    #find bound
    h = arr.size // 2

    #iterate through arrSize/2
    for i in range(0, h):

        #apply difference formula
        k = int(i * 2)
        tArr[i] = float(arr[k] * s0 + arr[k + 1] * s1)
        tArr[i + h] = float(arr[k] * w0 + arr[k + 1] * w1)

    #deep copy the temp to original array
    arr[0:] = np.copy(tArr)

    return

#Discrete 2D Haar Wavelet Transform
def Discrete2DWaveletTransform(imgArr):

    #get image dimensions
    numRows = imgArr.shape[0]
    numColumns = imgArr.shape[1]

    #iterate though rows
    for i in range(0, numRows):

        #copy current row into 1D numpy array
        currRow = np.copy(imgArr[i, :])

        #perform the transform on the current row
        Discrete1DWaveletTransform(currRow)

        #copy temp array to current row; shallow copy causes mayhem
        imgArr[i, :] = np.copy(currRow)

    #iterate though columns
    for j in range(0, numColumns):

        #copy current column into 1D numpy array
        currColumn = np.copy(imgArr[:, j])

        #perform the transform on the current column
        Discrete1DWaveletTransform(currColumn)

        #copy temp array to current column; shallow copy causes mayhem
        imgArr[:, j] = np.copy(currColumn)

    return

#Reconstruction of Discrete 1D Haar Wavelet Transform
def ReverseDiscrete1DWaveletTransform(arr):

    #declare empty array for current row/column
    tArr = np.empty(arr.size, np.float64)

    #find bound
    h = arr.size // 2

    #iterate through arrSize/2
    for i in range(0,h):

        #apply difference formula
        k = 2 * i
        tArr[k] = float((arr[i] * s0 + arr[i + h] * w0) / w0)
        tArr[k + 1] = float((arr[i] * s1 + arr[i + h] * w1) / s0)

    #copy values back
    arr[0:] = np.copy(tArr)

    return

#Reconstruction of 2D Haar Wavelet Transform
def ReverseDiscrete2DWaveletTransform(imgArr):

    #get image dimensions
    numRows = imgArr.shape[0]
    numColumns = imgArr.shape[1]

    #iterate through size of current column
    for j in range(0,numColumns):

    #copy current column into temp array
        currCol = np.copy(imgArr[:, j])

    #perform reverse transform on current column
        ReverseDiscrete1DWaveletTransform(currCol)

    #copy temp array back into column
        imgArr[:, j] = np.copy(currCol)

    #iterate through size of current column
    for i in range(0, numRows):

    #copy current row into temp array
        currRow = np.copy(imgArr[i, :])

    #perform reverse transform on current row
        ReverseDiscrete1DWaveletTransform(currRow)

        #copy temp array back into row
        imgArr[i, :] = np.copy(currRow)

    return

#function for normalizing values of an image to new scale
def normalize(img, newMax):

    #get min/max/range
    minV = img.min()
    maxV = img.max()
    rang = maxV - minV

    #normalize all values within [0, newMax]
    img = ((img - minV)  / rang) * 255

    return img

############### executable section:

def main():

	#read in image
	origImg = io.imread('lenag.png')

	#get dimensions of image
	height = origImg.shape[0]
	width = origImg.shape[1]

	#source for compressed image
	scaledImg = origImg.astype(np.float64) / 255

	#perform discrete transform on image for 8 iterations
	for i in range(9):

		#image number for save name
		strImageNumber = str(i + 1)

		#extract upper left 1/4 of previous region
		region = scaledImg[:int(height/pow(2,i)), :int(width/pow(2,i))]

		#perform the i-th layer transform
		Discrete2DWaveletTransform(region)

		#rescale image into [0,255] to enable visualization
		imgVisible = normalize(scaledImg, 255).astype(np.uint8)

		#save the image
		io.imsave('HDWT/layer_' + strImageNumber + '.png', imgVisible)

	#reconstruct the image from the 8th layer transformed image
	for k in range(9):
		cl = 8 - k
		region = scaledImg[:int(height/pow(2,cl)), :int(width/pow(2,cl))]
		ReverseDiscrete2DWaveletTransform(region)

	#variable for reconstructed image
	reconImg = np.copy(scaledImg)

	#rescale pixel intensities to 0,255
	reconImg = (reconImg * 255 + 0.5).astype(np.uint8)

	#find error to confirm that it should be zero
	err = np.abs(origImg-reconImg)
	print('Error is: ' + str(np.max(err)))

	#show reconstructed image
	viewer = io.imshow(reconImg)
	viewer.show()

	#save the image
	io.imsave('HDWT/reconImage.png', reconImg)

if __name__ == '__main__':
    main()

%d bloggers like this: