1. Median filter

1) Implement the median filter function based on the 5x5 template shown in the figure below (It is not allowed to directly or indirectly call third-party functions). Briefly describe the specific steps of median filtering.

2) Read in any image (grayscale or color image), add salt and pepper noise, display the original image, your own median filtering result and the third-party median filtering result, if the results are different, discuss the possible reasons.

2. Two-dimensional Separable Gaussian Smoothing

1) Students are required to implement a Gaussian smoothing function (It is not allowed to directly or indirectly call third-party cross-correlation, Gaussian smoothing, Gaussian blurring and related functions), the function parameters should at least include, image, template size, and Gaussian standard deviation. After smoothing, the size of the image should remain the same as the original image. The input image can be assumed to be a single-channel grayscale image.

2) Based on the separable nature of the Gaussian template, implement the two-dimensional Gaussian smoothing by convolving in the horizontal and vertical directions separately (It is not allowed to directly or indirectly call the third-party cross-correlation, and Gaussian smoothing, Gaussian blur and other functions cannot be called). Function parameters should remain consistent with 1).

3) Call the functions defined in 1) and 2) respectively to display a 3x3 grid, arranged in the following form:

      (Template size 13x13, standard deviation=2) original image, the result image of 1), the result image of 2),

      (Template size 13x13, standard deviation=4) original image, the result image of 1), the result image of 2),

      (Template size 13x13, standard deviation=6) original image, the result image of 1), the result image of 2),

4) Compare the Gaussian smoothing performance of 1) and 2) functions, set standard deviation=4, template size 2n+1, n=2...25, and record the running time of the function. Taking the template size as the horizontal coordinates and the running time as the vertical coordinates, plot the running time graphs of 1) and 2) functions. Analyze the time complexity of 1) and 2) functions (in big-O notation). Note: The image size can be adjusted based on your computer hardware (images larger than 512x512 may cause the program to run for a long time.) and the range of n.

我的做法：

from PIL import Image

import numpy as np

import CV2 as cv

import matplotlib.pyplot as plt

# 函数

    #处理灰度图 默认为3\n

def MedianFilter(img,k = 3, padding = None):

    imarray = np.array(img)

    height, width = imarray.shape

    if not padding:

            edge = int((k-1)/2)

            if height - 1 - edge <= edge or width - 1 - edge <= edge:

                print("没办法使用"),

                return None

            new_arr = np.zeros((height, width), dtype = "uint8")

            for i in range(height):

                for j in range(width):

                    #对于边缘不进行处理，仍然是原图像像素\n",

                    if i <= edge - 1 or i >= height - 1 - edge or j <= edge - 1 or j >= width - edge - 1:

                        new_arr[i, j] = imarray[i, j]

                    #可以用模板进行处理的点，其5*5内都有像素点了  

                    #                     (i,j-2)

                    #         （i-1,j-1）(i,j-1)(i+1,j-1)         

                    #(i-2，j)(i-1,j)（i，j）(i+1,j)(i+2,j)

                    #           (i-1.j+1)(i,j+1)（i+1,j+1）

                    #                       (i,j+2)

                    else:

                        b=[]

                        a=[]

                        a.append(imarray[i,j])

                        a.append(imarray[i,j+1])

                        a.append(imarray[i,j+2])

                        a.append(imarray[i,j-1])

                        a.append(imarray[i,j-2])

                        a.append(imarray[i-1,j])

                        a.append(imarray[i-1,j+1])

                        a.append(imarray[i-1,j-1])

                        a.append(imarray[i-2,j])

                        a.append(imarray[i+1,j-1])

                        a.append(imarray[i+1,j])

                        a.append(imarray[i+1,j+1])

                        a.append(imarray[i+2,j])

#                         a.append(imarray[i,j])

#                         a.append(imarray[i+1,j])

#                         a.append(imarray[i+2,j])

#                         a.append(imarray[i-1,j])

#                         a.append(imarray[i-2,j])

#                         a.append(imarray[i,j-1])

#                         a.append(imarray[i+1,j-1])

#                         a.append(imarray[i-1,j-1])

#                         a.append(imarray[i,j-2])

#                         a.append(imarray[i-1,j+1])

#                         a.append(imarray[i,j+1])

#                         a.append(imarray[i+1,j+1])

#                         a.append(imarray[i,j+2])

                        b.append(a)

#                         print(imarray[i,j])

#                         a=[imarray[i,j],imarray[i+1,j],imarray[i+2,j],imarray[i-1,j],imarray[i-2，j],imarray[i,j-1],imarray[i+1,j-1],imarray[i-1,j-1],imarray[i,j-2],imarray[i-1,j+1],imarray[i,j+1],imarray[i+1,j+1],imarray[i,j+2]]

#                         print(a)

#                         print(np.median(imarray[i - edge:i + edge + 1, j - edge:j + edge + 1]))

#                         print(np.median(b))

    

                        new_arr[i, j] = np.median(b)

            new_im = Image.fromarray(new_arr)

            new_im.save("chuli.png")

    

#处理彩色图像 默认为3

def median_Blur(img, filiter_size = 5):  #当输入的图像为彩色图像

        image_copy = np.array(img, copy = True).astype(np.float32)

#         print(image_copy.shape)

#         (751, 1283, 4)

        processed = np.zeros_like(image_copy)

        middle = int(filiter_size / 2)

        r = np.zeros(filiter_size * filiter_size)

        g = np.zeros(filiter_size * filiter_size)

        b = np.zeros(filiter_size * filiter_size)

        for i in range(middle, image_copy.shape[0] - middle):

            for j in range(middle, image_copy.shape[1] - middle):

                count = 0

                #取一点模板周围的像素值

                for m in range(i - middle, i + middle +1):

                    for n in range(j - middle, j + middle + 1):

#                         print(n)

#每一个点哪些取哪些不取

                        if (m==i-middle and n==j-middle) or (m==i-middle+1 and n==j-middle) or (m==i-middle+3 and n==j-middle) or (m==i-middle+4 and n==j-middle) or (m==i-middle and n==j-middle+1) or (m==i-middle+4 and n==j-middle+1) or (m==i-middle and n==j-middle+3) or (m==i-middle+4 and n==j-middle+3) or (m==i-middle and n==j-middle+4) or (m==i-middle+1 and n==j-middle+4) or (m==i-middle+3 and n==j-middle+4) or (m==i-middle+4 and n==j-middle+4):

                            pass

                        else:

                

                            r[count] = image_copy[m][n][0]

                            g[count] = image_copy[m][n][1]

                            b[count] = image_copy[m][n][2]

                            count += 1

#                 print(r)   

#                 print(g)

#                 print(b) 

                

                

                r.sort()

                g.sort()

                b.sort()

                processed[i][j][0] = r[int(13/2)]

                processed[i][j][1] = g[int(13/2)]

                processed[i][j][2] = b[int(13/2)]

        processed = np.clip(processed, 0, 255).astype(np.uint8)

        return processed

    #展示图片函数

def show(img):

        if img.ndim==2:

            plt.imshow(img,cmap='gray')

        else:

            plt.imshow(cv.cvtColor(img,cv.COLOR_BGR2RGB))

#     show(image)

    

    #噪声

    

    

def AddNoise(img, probility = 0.05, method ="salt_pepper"):#灰度图像加噪声

    

        imarray = np.array(img)

        height, width = imarray.shape

    

        for i in range(height):

            for j in range(width):

                if np.random.random(1) < probility:#随机加盐或者加椒    

                    if np.random.random(1) < 0.5:

                        imarray[i, j] = 0

                    else:

                        imarray[i, j] = 255

        new_im = Image.fromarray(imarray)

        new_im.save("zaosheng_huidu.png")

    

def RGBAddNoise(img,prob): #同时加杂乱(RGB单噪声)RGB图噪声 prob:噪声占比

        imarray = np.array(img)

        height, width, channels = imarray.shape

    

        #prob = 0.05 #噪声占比 已经比较明显了 >0.1 严重影响画质

        NoiseImg = imarray.copy()

        NoiseNum = int(prob * height * width)

        for i in range(NoiseNum):    

            rows = np.random.randint(0, height - 1)    

            cols = np.random.randint(0, width - 1)    

            channel = np.random.randint(0, 3)    

            if np.random.randint(0, 2) == 0:    

                NoiseImg[rows, cols, channel] = 0    

            else:        

                NoiseImg[rows, cols, channel] = 255

        #return NoiseImg

        new_im = Image.fromarray(NoiseImg)

        new_im.save("zaosheng_RGB.png")

# 使用plt显示CV2格式的图片

def plt_show_CV2_image(image,image2,image3):

#     image_0 = cv.cvtColor(image, cv.COLOR_BGR2RGB)

#     image_1 = cv.cvtColor(image, cv.COLOR_BGR2RGB)

#     image_2 = cv.cvtColor(image3, cv.COLOR_BGR2RGB)

    

    plt.figure(figsize=(400,300))  # 打开一个画布

    plt.axis('off')  # 不打开坐标轴

    plt.subplot(1,3,1)

    plt.title("original")

    plt.imshow(image)

    plt.subplot(1,3,2)

    plt.title("me")

    plt.imshow(image2)

    plt.subplot(1,3,3)

    plt.title("thrid-cv")

    plt.imshow(image3)

    plt.show()

  

image = cv.imread('orignal.png', cv.IMREAD_GRAYSCALE)

#保存灰度图后面用

cv.imwrite("huidu.png",image)   

addnoise=Image.open("huidu.png")

AddNoise(addnoise)

RGB_addnoise=Image.open("orignal.png")

RGBAddNoise(RGB_addnoise,0.05)

    

## 解释

## 我的方法是分灰度和彩色两种方式进行处理，当读入的是灰度图像时，先对图像大小进行评估看是否大于5*5，否则的话直接返回错误信息，对满足大小的图像

## 进行遍历像素点，对于边缘像素点直接用传入图像进行处理，像图像最开始的时候，其左上角是没有像素点的，以及左边等，相当于对整个图像的开头结尾左边右边

## 2行不进行中值划分，然后对有作用域的像素点取模板中的周边像素，并进行排序取其中值，然后替换到相应的像素点上，最后输出即可

## 对于彩色图像来说，同样是四个方向两行直接使用初始像素进行替代，对于拥有合格大小的像素区域进行遍历，分离出三通道人，r,g,b,进行储存

## 每次处理完一个像素点周围的时候，对其进行排序，找出中值并进行替换原来点的像素值，然后结束输出即可

#处理灰度图片

zaosheng_huidu = cv.imread('zaosheng_huidu.png')

zaosheng_huidu1=Image.open('zaosheng_huidu.png')

MedianFilter(zaosheng_huidu1,k=5)

#处理好的文件

#调用第三方cv处理图片

img_median = cv.medianBlur(zaosheng_huidu, 5)

chuli = cv.imread('chuli.png')

plt_show_CV2_image(zaosheng_huidu,chuli,img_median)

#处理彩色图像

zaosheng_RGB = cv.imread('zaosheng_RGB.png')

median_Blur=median_Blur(zaosheng_RGB, filiter_size = 5)

median_Blur=cv.cvtColor(median_Blur, cv.COLOR_BGR2RGB)

zaosheng_RGB2 = cv.cvtColor(zaosheng_RGB, cv.COLOR_BGR2RGB)#颜色模式转换

img_median = cv.medianBlur(zaosheng_RGB2,5)

plt.figure(figsize=(400,300))  # 打开一个画布

plt.axis('off')  # 不打开坐标轴

plt.subplot(1,3,1)

plt.title("original")

plt.imshow(zaosheng_RGB2)

plt.subplot(1,3,2)

plt.title("me")

plt.imshow(median_Blur)

plt.subplot(1,3,3)

plt.title("thrid-cv")

plt.imshow(img_median)

plt.show()

# 2. Two-dimensional Separable Gaussian Smoothing

# 1) Students are required to implement a Gaussian smoothing function (It is not allowed to directly or indirectly call third-party cross-correlation, Gaussian smoothing, Gaussian blurring and related functions), the function parameters should at least include, image, template size, and Gaussian standard deviation. After smoothing, the size of the image should remain the same as the original image. The input image can be assumed to be a single-channel grayscale image.

# 2) Based on the separable nature of the Gaussian template, implement the two-dimensional Gaussian smoothing by convolving in the horizontal and vertical directions separately (It is not allowed to directly or indirectly call the third-party cross-correlation, and Gaussian smoothing, Gaussian blur and other functions cannot be called). Function parameters should remain consistent with 1).

# 3) Call the functions defined in 1) and 2) respectively to display a 3x3 grid, arranged in the following form:

#       (Template size 13x13, standard deviation=2) original image, the result image of 1), the result image of 2),

#       (Template size 13x13, standard deviation=4) original image, the result image of 1), the result image of 2),

#       (Template size 13x13, standard deviation=6) original image, the result image of 1), the result image of 2),

# 4) Compare the Gaussian smoothing performance of 1) and 2) functions, set standard deviation=4, template size 2n+1, n=2...25, and record the running time of the function. Taking the template size as the horizontal coordinates and the running time as the vertical coordinates, plot the running time graphs of 1) and 2) functions. Analyze the time complexity of 1) and 2) functions (in big-O notation). Note: The image size can be adjusted based on your computer hardware (images larger than 512x512 may cause the program to run for a long time.) and the range of n.

#1

#高斯平滑函数

from PIL import Image

import numpy as np

from math import exp

import matplotlib.pyplot as plt

import CV2 as cv

 

def GaussianFilter(img, k = 3, sigma = 1):

 

    imarray = np.array(img)

    height, width = imarray.shape

    new_arr = np.zeros((height, width), dtype = "uint8")#新矩阵的替代

 

    filter = np.zeros((k,k))#设置模板大小

    mid = (k-1)/2#设置边缘

    sigma2 = pow(sigma, 2)

#     模板设置

    for i in range(k):

        for j in range(k):

            tmp = -(pow(i-mid,2)+pow(j-mid,2))/2/sigma2

            filter[i,j] = exp(tmp)/2/3.14/sigma2

    filtersum = filter.sum()

 

    for i in range(height):

        for j in range(width):

            total = 0

            #一个像素点，对模板内进行处理

            for n in range(pow(k,2)):

#k = 3, n = 0, 1, 2 ..., 8, a = -1, 0, 1, b = -1, 0, 1

#k = 5, n = 0, 1, 2, 3 ..., 24, a = -2, -1, 0, 1, 2

                a, b = int(n//k - (k-1)/2), int(n%k - 1)

                #filter_value

                aa, bb = int(n//k), int(n%k)

            #把模板这个点位的值取出来

                f_value = filter[aa, bb]

                if i + a <= 0:

                    if j + b <= 0:

                        total += imarray[0, 0]*f_value

                    elif j + b >= width - 1:

                        total += imarray[0, -1]*f_value

                    else:

                        total += imarray[0, j + b]*f_value

                elif i + a >= height - 1:

                    if j + b <= 0:

                        total += imarray[-1, 0]*f_value

                    elif j + b >= width - 1:

                        total += imarray[-1, -1]*f_value

                    else:

                        total += imarray[-1, j + b]*f_value

                else:

                    if j + b <= 0:

                        total += imarray[i + a, 0]*f_value

                    elif j + b >= width - 1:

                        total += imarray[i + a, -1]*f_value

                    else:

                        total += imarray[i + a, j + b]*f_value

            total /= filtersum

            new_arr[i, j] = int(total)

 

    new_im = Image.fromarray(new_arr)

    new_im.save("GaussianFilter.png")

#测试

# img=Image.open("huidu.png")

# # print(np.array(img).shape)

# GaussianFilter(img,k=7,sigma = 1)

#图片看其起来还行

#2

#二维高斯平滑函数 基于高斯模板的可分离性，通过水平方向和垂直方向分别卷积实现二维高斯平滑(不允许直接或间接调用第三方互相关，不能调用高斯平滑、高斯模糊等函数)。函数参数应与1)一致。

# python

#模板大小，图像，高斯标准差

def conv_2d(kernel_size, img,sigma,mode='fill'):

    kernel = np.zeros((kernel_size, kernel_size))

    center = kernel_size//2

    if sigma<=0:

        sigma = ((kernel_size-1)*0.5-1)*0.3+0.8

    

    s = sigma**2

    sum_val =  0

    for i in range(kernel_size):

        for j in range(kernel_size):

            x, y = i-center, j-center

            

            kernel[i, j] = np.exp(-(x**2+y**2)/(2*s))

            sum_val += kernel[i, j]

    

    kernel = kernel/sum_val

    

    

    

    

    if mode =='fill':

        h = kernel.shape[0]//2

        w = kernel.shape[1]//2

        # opencv filter2d 默认采用reflect填充方式

        # 前面测试了constant edge都不符合

        img = np.pad(img, (h,w), 'reflect')

    

    res_h = img.shape[0]-kernel.shape[0]+1

    res_w = img.shape[1]-kernel.shape[1]+1

    

    res = np.zeros((res_h, res_w))

    

    dh = kernel.shape[0]

    dw = kernel.shape[1]

    for i in range(res_h):

        for j in range(res_w):

            res[i,j] = np.sum(img[i:i+dh, j:j+dw]*kernel)

    cv.imwrite("fl_py.png",res)

#测试    结果与cv函数差不多

# img = np.random.random([5,5])

# img=np.array(Image.open("huidu.png"))

# conv_2d(3,img,5)

# #这样有色差问题

# image0 = cv.cvtColor(img, cv.COLOR_BGR2RGB)

# fl_py=cv.imread("fl_py.png")

# image1 = cv.cvtColor(fl_py, cv.COLOR_BGR2RGB)

# plt.figure(figsize=(400,300))

# plt.subplot(1,2,1)

# plt.title("original")

# plt.imshow(image0)

# plt.subplot(1,2,2)

# plt.title("chuli")

# plt.imshow(image1)

# plt.show()

#3)

img=Image.open("huidu.png")

GaussianFilter(img,k=13,sigma = 2)

fangfa1=cv.imread("GaussianFilter.png")

img_2=np.array(Image.open("huidu.png"))

conv_2d(13,img_2,2)

image0 = cv.cvtColor(img_2, cv.COLOR_BGR2RGB)

fangfa2=cv.imread("fl_py.png")

image1 = cv.cvtColor(fangfa2, cv.COLOR_BGR2RGB)

#模板尺寸13x13，标准差=2   原始图像，结果图像1，结果图像2，

plt.figure(figsize=(600,500))

plt.subplot(1,3,1)

plt.title("original")

plt.imshow(image0)

plt.subplot(1,3,2)

plt.title("result1")

plt.imshow(fangfa1)

plt.subplot(1,3,3)

plt.title("result2")

plt.imshow(fangfa2)

plt.show()

img=Image.open("huidu.png")

GaussianFilter(img,k=13,sigma = 4)

fangfa1=cv.imread("GaussianFilter.png")

img_2=np.array(Image.open("huidu.png"))

conv_2d(13,img_2,4)

image0 = cv.cvtColor(img_2, cv.COLOR_BGR2RGB)

fangfa2=cv.imread("fl_py.png")

image1 = cv.cvtColor(fangfa2, cv.COLOR_BGR2RGB)

#模板尺寸13x13，标准差=4  原始图像，结果图像1，结果图像2，

plt.figure(figsize=(600,500))

plt.subplot(1,3,1)

plt.title("original")

plt.imshow(image0)

plt.subplot(1,3,2)

plt.title("result1")

plt.imshow(fangfa1)

plt.subplot(1,3,3)

plt.title("result2")

plt.imshow(fangfa2)

plt.show()

img=Image.open("huidu.png")

GaussianFilter(img,k=13,sigma = 6)

fangfa1=cv.imread("GaussianFilter.png")

img_2=np.array(Image.open("huidu.png"))

conv_2d(13,img_2,6)

image0 = cv.cvtColor(img_2, cv.COLOR_BGR2RGB)

fangfa2=cv.imread("fl_py.png")

image1 = cv.cvtColor(fangfa2, cv.COLOR_BGR2RGB)

#模板尺寸13x13，标准差=6   原始图像，结果图像1，结果图像2，

plt.figure(figsize=(600,500))

plt.subplot(1,3,1)

plt.title("original")

plt.imshow(image0)

plt.subplot(1,3,2)

plt.title("result1")

plt.imshow(fangfa1)

plt.subplot(1,3,3)

plt.title("result2")

plt.imshow(fangfa2)

plt.show()

#4)

#记录各个模板大小的值

#记录各个函数分别运行的时间

import time

import os

import numpy

# 2n+1, n=2...25

randomByteArray = bytearray(os.urandom(120000))

flatNumpyArray = numpy.array(randomByteArray)

 

#数组转换为一个300*400的灰度图像

grayImage = flatNumpyArray.reshape(300,400)

cv.imwrite('RandomGray.png',grayImage)

#统计

fangfa1_time=[]

fangfa2_time=[]

size_moban=[]

for i in range(2,26):

    zhi=2*i+1#模板大小

    print(i)

    #运行第一个方法

    size_moban.append(zhi)

    img=np.array(img)

    old_time = time.time()

    GaussianFilter(img,k=zhi,sigma = 4)

    current_time = time.time()

    time_end=current_time - old_time

    fangfa1_time.append(time_end)

    

    

    

    

    

    #运行第二个方法

    old_time = time.time()

    conv_2d(13,img,2)

    current_time = time.time()

    time_end=current_time - old_time

    fangfa2_time.append(time_end)

#运行的数据已经有了    

#画图   板大小为横坐标，运行时间为纵坐标

plt.figure(figsize=(20, 8), dpi=80)

# 绘制图像 label 设置标签 color设置颜色 linestyle 设置线条 linewidth 设置线条粗细 alpha设置透明度

plt.plot(size_moban,fangfa1_time,label='one_time', color='red', linestyle=':', marker='.', markersize=5)

plt.plot(size_moban,fangfa2_time,label='two_time', color='black',linestyle='-', marker='.', markersize=5)

# 设置X，Y轴标签

plt.xlabel('size')

plt.ylabel('time')

plt.show()

for i in fangfa2_time:

    print(i)

# 方法二运行时间  

输出：

0.8112778663635254 0.7058887481689453 0.7052502632141113 0.7029576301574707 0.7056605815887451 0.749387264251709 0.7981278896331787 0.7647967338562012 0.7728753089904785 0.8714818954467773 0.8129367828369141 0.7645151615142822 0.9232168197631836 0.7525970935821533 0.8027229309082031 1.0895204544067383 0.7871663570404053 1.352036476135254 0.7733709812164307 0.7695446014404297 0.7612559795379639 0.9162437915802002 0.7639052867889404 1.0024442672729492