Update 2 files

- /Deblur_Algorithm.py
- /F16.tiff

Deblur_Algorithm.py

+"""
+Deblur Algorithm
+
+input image : F16.tiff
+filters : Mean Filter, Gaussian Filter
+paddings : zero, replicate, mirror
+"""
+
+import cv2
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.sparse import csr_matrix
+from scipy.sparse.linalg import spsolve
+from datetime import datetime
+from pytz import timezone
+
+# calculate global matrix M in CSR format
+def M(input_data, filter_data, stride):
+    num_elements= n*m*filter_data.size
+    # allocate memory for data, row indices, column indices
+    data=np.zeros(num_elements)
+    row_indices= np.zeros(num_elements, dtype=int)
+    col_indices= np.zeros(num_elements, dtype=int)
+    # initialize counters for data storage
+    count_i = 0 # for row indices
+    count_data=0 # for the current position
+    # operation for convolution 
+    for i in range(0, input_data.shape[0] - filter_data.shape[0] + 1, stride):
+        for j in range(0, input_data.shape[1] - filter_data.shape[1] + 1, stride):
+            tmp = np.array(input_data[i:i + filter_data.shape[0], j:j + filter_data.shape[1]]).ravel()
+            tmp2 = np.nonzero(tmp)[0]
+            for k in range(len(tmp2)):
+                data[count_data]= filter_data.ravel()[tmp2[k]]
+                row_indices[count_data]= count_i
+                col_indices[count_data]= tmp[tmp2[k]] - 1
+                count_data +=1
+            count_i += 1
+    data= data[:count_data]
+    row_indices=row_indices[:count_data]
+    col_indices=col_indices[:count_data]
+    result_sparse = csr_matrix((data, (row_indices, col_indices)), shape=(n * m, n * m)) #store in CSR format
+    return result_sparse
+
+# solve MX=Y using spsolve 
+def solve(M_sparse, img_blur):
+    try:
+        Y = np.array(img_blur)
+        img_deblur = spsolve(M_sparse, Y.ravel())
+        img_deblur = np.reshape(img_deblur, (n, m))
+    # if inverse matrix doesn't exists, exit the program
+    except (np.linalg.LinAlgError, RuntimeError): 
+        img_deblur = np.zeros((n, m))
+        print("Inverse of M does not exist.")
+        exit()
+    return img_deblur
+
+# calculate errors 
+def InfNorm(A, B):
+    if A.shape != B.shape:
+        raise ValueError("The shapes of two matrixes must be the same")
+    diff = np.abs(A - B)
+    infnorm = np.max((diff))
+    return infnorm
+
+def apply_filter_and_deblur(img, K):
+    # blur using openCV
+    if border == 0:
+        img_blur = cv2.filter2D(img, -1, K, borderType=cv2.BORDER_CONSTANT)
+    elif border == 1:
+        img_blur = cv2.filter2D(img, -1, K, borderType=cv2.BORDER_REPLICATE)
+    elif border == 2:
+        img_blur = cv2.filter2D(img, -1, K, borderType=cv2.BORDER_REFLECT_101)
+    # calculate M
+    M_sparse = M(border_X, K, 1)
+    # solve MX=Y using spsolve
+    img_deblur = solve(M_sparse, img_blur)
+    # calculate errors between original image and deblurred image
+    error = InfNorm(img128, img_deblur)
+    return img_blur, img_deblur, error
+
+# input image
+img = cv2.imread("images/F16.tiff", cv2.IMREAD_GRAYSCALE)
+img128 = cv2.resize(img, (128, 128))
+img128 = img128.astype(np.float64)/255.0 # image normalization
+
+# show results using matplotlib
+fig = plt.figure()
+rows = 1
+cols = 3
+
+# result 1. original image
+ax1 = fig.add_subplot(rows, cols, 1)
+ax1.imshow(img128, cmap='gray')
+ax1.set_title('F16')
+ax1.axis("off")
+
+n, m = 128, 128
+X = (np.arange(1, n * m + 1)).reshape(n, m) 
+
+# choose the padding type
+border_type = ['zero', 'replicate', 'mirror']
+border = int(input(f"{border_type[:]}(Choose from 0-2) : " ))
+if border == 0: # zero padding
+    border_X = np.pad(X, ((1, 1), (1, 1)), mode='constant', constant_values=0) 
+elif border == 1: # replicate padding
+    border_X = np.pad(X, ((1, 1), (1, 1)), mode='edge')
+elif border == 2: # mirror padding
+    border_X = np.pad(X, ((1, 1), (1, 1)), mode='reflect')
+
+# choose the filter type 
+filter_type = ['Mean', 'Gaussian']
+kernel_size = 3
+kernel = int(input(f"{filter_type[:]}(Choose from 0-1) : " ))
+if kernel == 0: # Mean filter
+    K = np.ones((kernel_size,kernel_size),np.float64)/(kernel_size*kernel_size)
+    img_blur, img_deblur, min_error = apply_filter_and_deblur(img128, K)
+elif kernel == 1: # Gaussian filter
+    min_error = float('inf')
+    best_sigma = None
+    sigma_values = np.arange(1, 101, 1)
+    # find best sigma values
+    for sigma in sigma_values:
+        tmp3 = cv2.getGaussianKernel(kernel_size, sigma)
+        K = np.outer(tmp3, tmp3.T)
+        img_blur, img_deblur, error = apply_filter_and_deblur(img128, K)
+        if error < min_error:
+            min_error = error
+            best_sigma = sigma
+            best_deblur = img_deblur
+    img_deblur = best_deblur
+    print("Best sigma value: ", best_sigma)
+
+# result 2. blurred image using openCV
+ax2 = fig.add_subplot(rows, cols, 2)
+ax2.imshow(img_blur, cmap='gray')
+ax2.set_title(f'{filter_type[kernel]}')
+ax2.axis("off")
+
+#result 3. deblurred image using M()
+ax3 = fig.add_subplot(rows, cols, 3)
+ax3.imshow(img_deblur, cmap='gray')
+ax3.set_title('Deblur')
+ax3.axis("off")
+
+# print error
+print("Error:", min_error)
+
+# save results
+date_today = datetime.now(timezone('Asia/Seoul')).strftime("%Y%m%d_%H%M")
+plt.savefig('Results/{}.jpg'.format(date_today), format='jpg')
+plt.close()
\ No newline at end of file