!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
import cv2
import copy
import warnings
from scipy import fftpack
import utils
Hit:1 http://security.debian.org/debian-security buster/updates InRelease
Hit:2 http://deb.debian.org/debian buster InRelease
Get:3 http://deb.debian.org/debian buster-updates InRelease [51.9 kB]
Fetched 51.9 kB in 0s (151 kB/s)
2 packages can be upgraded. Run 'apt list --upgradable' to see them.
ffmpeg is already the newest version (7:4.1.6-1~deb10u1).
libsm6 is already the newest version (2:1.2.3-1).
libxext6 is already the newest version (2:1.3.3-1+b2).
0 upgraded, 0 newly installed, 0 to remove and 2 not upgraded.
Requirement already satisfied: opencv-python in /root/venv/lib/python3.7/site-packages (4.5.1.48)
Requirement already satisfied: numpy>=1.14.5 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from opencv-python) (1.19.5)
WARNING: You are using pip version 21.0.1; however, version 21.1.1 is available.
You should consider upgrading via the '/root/venv/bin/python -m pip install --upgrade pip' command.
# 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
denoise = utils.filter_image(img, N, M)
return denoise
plt.imshow(denoise(im_noisy), cmap = 'gray')
plt.show()
def compute_MSE(img, img_new):
err = 0
for i in range(N-1):
for j in range(M-1):
err += (img_new[i,j]-img[i,j])**2
MSE = (1/N)*(1/M)*err
return 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 = np.fft.fft2(X)
# Compute the magnitudes of the coefficients
mgs = np.log(np.abs(cfts))
return cfts, 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
inv = np.fft.ifft2(cfs)
# Get the magnitudes
inv_mags = np.log(np.abs(inv))
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 = np.fft.ifft2(cfs)
imags = np.log(np.abs(icfs))
#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):
# Set empty array for dct coefficients
dct = np.zeros((N, N), dtype=np.float64)
# Iterate through and compute values for all coefficients
for k in range(N):
for l in range(N):
coeffs = 0
################
if k == 0:
c1 = 1/np.sqrt(2)
else:
c1 = 1
if l == 0:
c2 = 1/np.sqrt(2)
else:
c2 = 1
for m in range(34):
for n in range(34):
coeffs += X[m,n]*c1*c2*np.cos(((k*np.pi)/(2*N))*(2*m + 1))*np.cos(((l*np.pi)/(2*N))*(2*n + 1))
dct[k,l] = coeffs
################
# 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 = int(N/8)
def partition_img(img):
ppatches = np.zeros((prange,prange, 8, 8))
for i in range(prange):
for j in range(prange):
for k in range(8):
ppatches[i][j][k] = img[8*i+k][8*j : 8*j+8]
return ppatches
def partitioned_dct(ppatches, K):
# Set empty array of sizes K^2 per patch
dct_patches = np.zeros((prange,prange,8,8))
for row in range(prange):
for column in range(prange):
coefficients, mags = get_coeffs(ppatches[row][column])
sortedmags = np.sort(mags)
mags = []
for mag in sortedmags:
for item in mag:
mags.append(item)
mags.sort(reverse = True)
mags_new = mags[0:K**2]
for i in range(8):
for j in range(8):
if sortedmags[i][j] in mags_new:
dct_patches[row][column][i][j] = sortedmags[i][j]
else:
dct_patches[row][column][i][j] = 0
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 = np.shape(dct_patches)
# Initialize an empty matrix to store the quantized coefficients
q_dct = np.zeros((m,n,k,l))
# Compute and store the quantized coefficients
for row in range(m):
for column in range(n):
for i in range(k):
for j in range(l):
q_dct[row][column][i][j] = round((dct_patches[row][column][i][j]/Qmtrx[i][j]),5)
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 = np.shape(qdct_patches)
# Initialize an empty matrix to store the quantized coefficients
iq_dct = np.zeros((m,n,k,l))
# Compute and store the quantized coefficients
for row in range(m):
for column in range(n):
for i in range(k):
for j in range(l):
iq_dct[row][column][i][j] = qdct_patches[row][column][i][j]*Qmtrx[i][j]
return iq_dct
from scipy import fftpack
# Computes iDCTs for patches
def inverse_dct(patches):
# Set empty array for recovered set of patches image
i_patches = np.zeros((prange,prange, 8, 8), dtype=np.float64)
#iterate through each patch
for row in range(prange):
for column in range(prange):
#For each patch iterate through and compute values for all coefficients
for m in range(8):
for n in range(8):
coeffs = 0
if m == 0:
c1 = 1/np.sqrt(2)
else:
c1 = 1
if n == 0:
c2 = 1/np.sqrt(2)
else:
c2 = 1
for k in range(8):
for l in range(8):
coeffs += patches[row][column][k][l]*c1*c2*np.cos((m*np.pi*(2*k+1))/2*8)*np.cos((n*np.pi*(2*l+1))/2*8)
i_patches[row][column][m,n] = coeffs
#normalize the transform
i_patches = 2*i_patches/np.sqrt(8**2)
"""
fftpack.idctn(patches[i,j], type=2, shape=[8,8], norm='ortho')
"""
return i_patches
def stitch(patches):
recovered = np.zeros((N,N), dtype=np.float64)
for i in range(prange):
for j in range(prange):
for k in range(8):
recovered[8*i+k][8*j : 8*j+8] = patches[i][j][k]
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):
err = 0
for i in range(N-1):
for j in range(N-1):
err += (img_new[i,j]-img[i,j])**2
MSE = (1/N)*(1/N)*err
return MSE
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()