#!apt update #!apt-get install ffmpeg libsm6 libxext6 -y #!pip install opencv-python from PIL import Image import matplotlib.pyplot as plt import numpy as np from scipy import ndimage from cv2 import cv2 import copy import warnings from scipy import fftpack import utils

# Height and width of our images in parts 1 and 2 M = 512 N = 512

im_original = cv2.imread("images/puppy1.png", 0) im_noisy = utils.noise_image(im_original, 6.0, N, M)

def denoise(img): #Hint: check out the "filter image" function in utils return utils.filter_image(img, N, M)

plt.imshow(denoise(im_noisy), cmap = 'gray') plt.show()

def compute_MSE(img, img_new): M, N = img.shape error = 0 for i in range(N): for j in range(M): error += (img[i][j] - img_new[i][j])**2 return error/(N*M) noises = [0,2,4,6,8,10,12,14,16] mse = [] for noise in noises: noisy_img = utils.noise_image(im_original, noise, N, M) reconstructed_im = denoise(noisy_img) mse.append(compute_MSE(im_original, reconstructed_im)) plt.plot(noises, mse) plt.title("MSE vs Amount of gausian noise") plt.xlabel("Amount of gausian noise") plt.ylabel("Reconstruction mse")

noise_amounts = np.linspace(0,16,8, dtype=np.float64) errors = [] for i in noise_amounts: noisy_im = utils.noise_image(im_original, i, M, N) filtered_im = denoise(noisy_im) errors.append(compute_MSE(im_original, filtered_im)) goal_curve = [90, 95, 100, 120, 170, 250, 375, 550] plt.plot(noise_amounts, errors, 'bo--', label='errors') plt.plot(noise_amounts, goal_curve, 'gs--', label='goal') plt.legend() plt.xlabel('Gaussian Noise Added') plt.ylabel('Reconstruction Error') plt.show()

im1 = cv2.imread("images/kitten1.png",0) plt.imshow(im1, cmap='gray') plt.show()

def get_coeffs(X): # Compute discrete fourier for the image cfts = cv2.dft(np.float32(X),flags = cv2.DFT_COMPLEX_OUTPUT) cfts_shift = np.fft.fftshift(cfts) # Compute the magnitudes of the coefficients mgs = np.log(cv2.magnitude(cfts_shift[:,:,0],cfts_shift[:,:,1])) return cfts_shift, mgs

coefficients, mags = get_coeffs(im1)

# Pick 2 values for alpha to show the difference in quality of reconstruction alpha1 = .99 alpha2 = .01 coeffs1 = utils.get_alpha_coeffs(coefficients, mags, alpha1, N, M) coeffs2 = utils.get_alpha_coeffs(coefficients, mags, alpha2, N, M)

def inverse_ft(cfs): #Get the IDFT f_ishift = np.fft.ifftshift(cfs) inv = cv2.idft(f_ishift) # Get the magnitudes inv_mags = cv2.magnitude(inv[:,:,0],inv[:,:,1]) return inv, inv_mags

inverse1, inverse_mags1 = inverse_ft(coeffs1) inverse2, inverse_mags2 = inverse_ft(coeffs2)

f = plt.figure(figsize=(12,12)) plt.subplot(131),plt.imshow(im1, cmap = 'gray') plt.xlabel('Orginal', fontsize=8), plt.xticks([]), plt.yticks([]) plt.subplot(132),plt.imshow(inverse_mags1, cmap = 'gray') plt.xlabel("alpha = " + str(alpha1), fontsize=8), plt.xticks([]), plt.yticks([]) plt.subplot(133),plt.imshow(inverse_mags2, cmap = 'gray') plt.xlabel("alpha = " + str(alpha2), fontsize=8), plt.xticks([]), plt.yticks([]) plt.show()

# Computes reconstruction error def compute_error(img, img_new): return compute_MSE(img, img_new)

alphas = np.linspace(.01,.99,25, dtype=np.float64) errors = np.zeros(len(alphas)) for a in range(len(alphas)): cfs = utils.get_alpha_coeffs(coefficients, mags, alphas[a], N, M) icfs = cv2.idft(cfs) imags = cv2.magnitude(icfs[:,:,0],icfs[:,:,1])/(M*N) errors[a] = compute_error(im1,imags) plt.plot(alphas, errors, 'o') plt.xlabel('alpha') plt.ylabel('Reconstruction Error') plt.title('FT Compression/Distortion Trade-Off') plt.show()

# Redefine image dimensions for part 3.1 N = 35

def compute_dct(X): N, N = X.shape # Set empty array for dct coefficients dct = np.zeros((N, N), dtype=np.float64) #discrete cosine def cos_discrete(m,n,k,l, M, N): return np.cos(k*np.pi*(2*m+1)/(2*M))*np.cos(l*np.pi*(2*n+1)/(2*N)) # Iterate through and compute values for all coefficients for k in range(N-1): if k == 0: c1 = 1/np.sqrt(2) else: c1 = 1 for l in range(N-1): if l == 0: c2 = 1/np.sqrt(2) else: c2 = 1 pass for m in range(N-1): for n in range(N-1): dct[k][l] += X[m][n]*c1*c2*cos_discrete(m,n,k,l,N,N) # normalize the transform dct = 2*dct/np.sqrt(N**2) return dct

img = cv2.imread("images/dct_sample1.png", 0) dct = compute_dct(img) dct2 = fftpack.dct(fftpack.dct(img.T, type=2, norm='ortho').T, type=2, norm='ortho') # Plot the image next to the transform and the scipy transform f = plt.figure(figsize=(12,12)) plt.subplot(131),plt.imshow(img, cmap = 'gray') plt.xlabel('', fontsize=8), plt.xticks([]), plt.yticks([]) plt.title('Original Image') plt.subplot(132),plt.imshow(dct, cmap = 'gray') plt.xlabel('', fontsize=8), plt.xticks([]), plt.yticks([]) plt.title('Results of DCT Function') plt.subplot(133),plt.imshow(dct2, cmap = 'gray') plt.xlabel('', fontsize=8), plt.xticks([]), plt.yticks([]) plt.title('Scipys DCT Function') plt.show()

# Redefine image dimensions for part 3.1 N = 256

im1 = cv2.imread("images/puppy3_small.png",0) plt.imshow(im1, cmap='gray') plt.show()

# Number of patches prange = N/8

def partition_img(img): m,n = img.shape patches = np.empty((int(m/8), int(n/8), 8, 8)) for x in range(0,img.shape[0],8): for y in range(0,img.shape[1],8): index_x = int(x/8) index_y = int(y/8) patches[index_x][index_y] = img[x:x+8,y:y+8] return patches

def partitioned_dct(ppatches, K): # Set empty array of sizes K^2 per patch m, n, k, l = ppatches.shape dct_patches = np.zeros((m,n,8,8)) for row in range(m): for col in range(n): #get the dct of a patch dct = compute_dct(ppatches[row][col]) #get indices of largest k^2 elements in that patch ind = np.unravel_index(np.argsort(dct, axis=None), dct.shape) index_row_largest = ind[0][-K**2:] index_col_largest = ind[1][-K**2:] indices = list(zip(index_row_largest, index_col_largest)) #store those elements for element in indices: dct_patches[row][col][element[0]][element[1]] = dct[element[0]][element[1]] return dct_patches

patches = partition_img(im1) for i in range(8): for j in range(8): plt.subplot(8,8,int(1+j+i*8)), plt.imshow(patches[i][j], cmap="gray") plt.xticks([]) plt.yticks([]) plt.show()

partitioned_img = partition_img(im1) dct4 = partitioned_dct(partitioned_img, 2) dct16 = partitioned_dct(partitioned_img, 4) dct36 = partitioned_dct(partitioned_img, 6) dct64 = partitioned_dct(partitioned_img, 8) plt.subplot(141),plt.imshow(dct4[0][0], cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(142),plt.imshow(dct16[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(143),plt.imshow(dct36[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(144),plt.imshow(dct64[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.show()

# Compute the quantized version of the specified transform (DCT or DFT) def quantize(dct_patches): # Initialize the quantization matrix Qmtrx = np.zeros((8,8)) Qmtrx[0,:] = [16, 11, 10, 16, 24, 40, 51, 61] Qmtrx[1,:] = [12, 12, 14, 19, 26, 58, 60, 55] Qmtrx[2,:] = [14, 13, 16, 24, 40, 57, 69, 56] Qmtrx[3,:] = [14, 17, 22, 29, 51, 87, 80, 62] Qmtrx[4,:] = [18, 22, 37, 56, 68, 109, 103, 77] Qmtrx[5,:] = [24, 36, 55, 64, 81, 104, 113, 92] Qmtrx[6,:] = [49, 64, 78, 87, 103, 121, 120, 101] Qmtrx[7,:] = [72, 92, 95, 98, 112, 100, 103, 99] # Store dimensions of DCT coefficient matrix patches for iteration m, n, k, l = dct_patches.shape # Initialize an empty matrix to store the quantized coefficients q_dct = np.empty((m,n,k,l)) # Compute and store the quantized coefficients for row in range(m): for col in range(n): for row_patch in range(k): for col_patch in range(l): q_dct[row][col][row_patch][col_patch] = np.round(dct_patches[row][col][row_patch][col_patch]/Qmtrx[row_patch][col_patch]) return q_dct

qdct4 = quantize(dct4) qdct16 = quantize(dct16) qdct36 = quantize(dct36) qdct64 = quantize(dct64) plt.subplot(141),plt.imshow(qdct4[0][0], cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(142),plt.imshow(qdct16[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(143),plt.imshow(qdct36[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.subplot(144),plt.imshow(qdct64[0][0],cmap='gray'), plt.xticks([]),plt.yticks([]) plt.show()

def unquantize(qdct_patches): # Initialize the quantization matrix Qmtrx = np.zeros((8, 8)) Qmtrx[0, :] = [16, 11, 10, 16, 24, 40, 51, 61] Qmtrx[1, :] = [12, 12, 14, 19, 26, 58, 60, 55] Qmtrx[2, :] = [14, 13, 16, 24, 40, 57, 69, 56] Qmtrx[3, :] = [14, 17, 22, 29, 51, 87, 80, 62] Qmtrx[4, :] = [18, 22, 37, 56, 68, 109, 103, 77] Qmtrx[5, :] = [24, 36, 55, 64, 81, 104, 113, 92] Qmtrx[6, :] = [49, 64, 78, 87, 103, 121, 120, 101] Qmtrx[7, :] = [72, 92, 95, 98, 112, 100, 103, 99] # Store dimensions of DCT coefficient matrix patches for iteration m, n, k, l = qdct_patches.shape # Initialize an empty matrix to store the unquantized coefficients iq_dct = np.empty((m, n, k, l)) # Compute and store the quantized coefficients for row in range(m): for col in range(n): for row_patch in range(k): for col_patch in range(l): iq_dct[row][col][row_patch][col_patch] = ( qdct_patches[row][col][row_patch][col_patch]* Qmtrx[row_patch][col_patch] ) return iq_dct

# Computes iDCTs for patches def inverse_dct(patches): # Set empty array for recovered set of patches image i_patches = np.zeros((int(prange), int(prange), 8, 8), dtype=np.float64) M = 8 N = 8 # inverse discrete cosine def icos_discrete(m, n, k, l, M, N): return np.cos(m * np.pi * (2 * k + 1) / (2 * M)) * np.cos( n * np.pi * (2 * l + 1) / (2 * N) ) # Iterate through and compute values for all coefficients for patch_row in range(int(prange)): for patch_col in range(int(prange)): for m in range(8): if m == 0: c1 = 1 / np.sqrt(2) else: c1 = 1 for n in range(8): if n == 0: c2 = 1 / np.sqrt(2) else: c2 = 1 # for a given patch, for a given component in that patch for k in range(8): for l in range(8): i_patches[patch_row][patch_col][m][n] += ( patches[patch_row][patch_col][k][l] * c1 * c2 * icos_discrete(m, n, k, l, N, N) ) # normalize the transform i_patches[patch_row][patch_col] = ( 2 * i_patches[patch_row][patch_col] / np.sqrt(N ** 2) ) return i_patches

def stitch(patches): N = 256 recovered = np.zeros((N, N), dtype=np.float64) for row in range(0,256,8): for col in range(0,256,8): recovered[row:row+8, col:col+8] = patches[int(row/8)][int(col/8)] return recovered

recovered4 = stitch(inverse_dct(unquantize(qdct4))) recovered16 = stitch(inverse_dct(unquantize(qdct16))) recovered36 = stitch(inverse_dct(unquantize(qdct36))) recovered64 = stitch(inverse_dct(unquantize(qdct64))) plt.subplot(141),plt.imshow(recovered4,cmap='gray'), plt.title('K^2=4',fontsize=8), plt.xticks([]), plt.yticks([]) plt.subplot(142),plt.imshow(recovered16, cmap='gray'), plt.title('K^2 = 16',fontsize=8), plt.xticks([]),plt.yticks([]) plt.subplot(143),plt.imshow(recovered36,cmap='gray'), plt.title('K^2 = 36',fontsize=8), plt.xticks([]),plt.yticks([]) plt.subplot(144),plt.imshow(recovered64, cmap='gray'), plt.title('K^2 = 64',fontsize=8), plt.xticks([]),plt.yticks([]) plt.show()

# Computes reconstruction error def compute_error(img, img_new): return compute_MSE(img, img_new)

dcq_errors = [compute_error(im1, q) for q in [recovered4, recovered16, recovered36, recovered64]] K2list = [4,16,36,64] jpeg1 = cv2.imread("images/puppy3_small1.jpg", 0) j1 = compute_error(im1, jpeg1) j1 = [j1, j1, j1, j1] jpeg2 = cv2.imread("images/puppy3_small2.jpg", 0) j2 = compute_error(im1, jpeg2) j2 = [j2, j2, j2, j2] jpeg3 = cv2.imread("images/puppy3_small3.jpg", 0) j3 = compute_error(im1, jpeg3) j3 = [j3, j3, j3, j3] jpeg4 = cv2.imread("images/puppy3_small4.jpg", 0) j4 = compute_error(im1, jpeg4) j4 = [j4, j4, j4, j4] # Plot the MSE Reconstruction Errors against the Jpeg standard plt.plot(K2list, dcq_errors, 'k-s', label='Quantized DCT') plt.plot(K2list, j1, 'r--', label='JPEG Low Standard') plt.plot(K2list, j2, 'y--', label='JPEG Medium Standard') plt.plot(K2list, j3, 'g--', label='JPEG High Standard') plt.plot(K2list, j4, 'b--', label='JPEG Max Standard') plt.legend() plt.xlabel('K^2') plt.ylabel('Reconstruction Error') plt.title('DFT and DCT Reconstruction Error for Varying K^2') plt.savefig('JPEG_standard_plot.png') plt.show()