关于高斯滤波与高斯模糊

上次突然被人问啥是高斯模糊，这个听起来如此熟悉，又如此模糊，忘记了，亏自己还是个数学系的学生，真让人汗颜，下面就回顾一下。

首先先说一下高斯滤波，高斯滤波是一种图像处理方法，能够使图片变得更平滑，去掉图像中的噪声（噪声在图像当中常表现为一引起较强视觉效果的孤立像素点或像素块。简单来说，噪声的出现会给图像带来干扰，让图像变得不清楚）。

高斯滤波类似于深度学习中的卷积，作用于整张图像后得到想要效果的图片。
首先高斯核，即高斯函数一维与二维，函数形式如下，原点即为（0，0）：

注：σ的大小决定了高斯函数的宽度。

高斯模糊即对图片做平均的一个过程，具体的，“中间点”取”周围点”的平均值，如下图，2周围8个点的均值为1，处理完后像素值就会变成1。在数值上，这是一种”平滑化”。在图形上，就相当于产生”模糊”效果，”中间点”失去细节。

显然，计算平均值时，取值范围越大，”模糊效果”越强烈。


上面分别是原图、模糊半径3像素、模糊半径10像素的效果。模糊半径越大，图像就越模糊。从数值角度看，就是数值越平滑。

如果使用简单平均，显然不是很合理，因为图像都是连续的，越靠近的点关系越密切，越远离的点关系越疏远。因此，加权平均更合理，距离越近的点权重越大，距离越远的点权重越小。

而正态分布显然是一种可取的权重分配模式。由于图像是二维的，所以需要使用二维的高斯函数，具体权重矩阵咋算的，如下，
假定中心点的坐标是（0,0），那么距离它最近的8个点的坐标则为：

假定sigma=1.5, 则上面九个点的权重矩阵为

例如（0，0）处的权重，就是把（0，0）带入二维的高斯函数，计算的数值则为对应的权重，f(0,0)=1/(23.14151.5*1.5)=1/14.13675=0.0707

这9个点的权重总和等于0.4787147，如果只计算这9个点的加权平均，还必须让它们的权重之和等于1，因此上面9个值还要分别除以0.4787147，得到最终的权重矩阵：

有了权重矩阵，就可以计算高斯模糊的值了。

中心点以及周边n个点，每个点乘以自己的权重值并将这些值相加，即将上图最后9个值加起来，就是中心点的高斯滤波的值。
对所有点重复这个过程，就得到了高斯模糊后的图像。

如果原图是彩色图片，可以对RGB三个通道分别做高斯模糊。

对于边界点来说，周边没有足够的点，一个变通方法是把已有的点拷贝到另一面的对应位置，模拟出完整的矩阵。

代码部分：

#对上文高斯模糊的方法阐释的相关代码  -xlxw
from PIL import Image as p
import numpy as np
from time import clock
import math

#define
sizepic = [0,0]
timer = [0,0,0,0]
PI = 3.1415926

def getrgb(path):#得到图像中各个点像素的RGB三通道值
    timer[0]=clock()
    pd = p.open(path)
    sizepic[0] = pd.size[0]
    sizepic[1] = pd.size[1]
    nr = np.empty((sizepic[0],sizepic[1]))
    ng = np.empty((sizepic[0],sizepic[1]))
    nb = np.empty((sizepic[0],sizepic[1]))
    for i in range(0,sizepic[0]):
        for j in range(0,sizepic[1]):
            nr[i][j] = pd.getpixel((i,j))[0]
            ng[i][j] = pd.getpixel((i,j))[1]
            nb[i][j] = pd.getpixel((i,j))[2]
    print("已经得到所有像素的R,G,B的值，所花时间为{:.3f}s".format(clock()-timer[0]))
    return nr,ng,nb
    
def Matrixmaker(r):#通过半径和坐标计算高斯函数矩阵
    summat = 0
    timer[1] = clock()
    ma = np.empty((2*r+1,2*r+1))
    for i in range(0,2*r+1):
        for j in range(0,2*r+1):
            gaussp = (1/(2*PI*(r**2))) * math.e**(-((i-r)**2+(j-r)**2)/(2*(r**2))) 
            ma[i][j] = gaussp
            summat += gaussp
    print(ma)
    print(summat)
    for i in range(0,2*r+1):
        for j in range(0,2*r+1):
            ma[i][j] = ma[i][j]/summat
    print("已经计算出高斯函数矩阵，所花时间为{:.3f}s".format(clock()-timer[1]))
    print("矩阵如下:")
    return ma
          
def newrgb(ma,nr,ng,nb,r):#生成新的像素rgb矩阵
    timer[2] = clock()
    newr = np.empty((sizepic[0],sizepic[1]))
    newg = np.empty((sizepic[0],sizepic[1]))
    newb = np.empty((sizepic[0],sizepic[1]))
    for i in range(r+1,sizepic[0]-r):
        for j in range(r+1,sizepic[1]-r):
            o = 0 
            for x in range(i-r,i+r+1):
                p = 0
                for y in range(j-r,j+r+1):
                    #print("x{},y{},o{},p{}".format(x,y,o,p))
                    newr[i][j] += nr[x][y]*ma[o][p]
                    newg[i][j] += ng[x][y]*ma[o][p]
                    newb[i][j] += nb[x][y]*ma[o][p]
                    p += 1
                o += 1            
    print("已经计算出新的三通道矩阵，所花时间为{:.3f}s".format(timer[2]))
    return newr,newg,newb

def cpic(r,g,b,path,rd):
    timer[3] = clock()
    pd = p.open(path)
    for i in range(rd+1,sizepic[0]-rd+1):
        for j in range(rd+1,sizepic[1]-rd+1):
            pd.putpixel((i,j),(int(r[i][j]),int(g[i][j]),int(b[i][j])))
    print("已经完成生成，所花时间为{:.3f}s".format(timer[3]))
    print("正在导出图片..")
    pd.save("blurred.jpg")
    
def main():
    rd = eval(input("请输入模糊的半径："))
    path = input("请输入图片的地址(包括后缀)：")
    nr,ng,nb = getrgb(path)
    matx = Matrixmaker(rd)
    print(matx)
    print("正在转换..")
    newr,newg,newb = newrgb(matx,nr,ng,nb,rd)
    print("正在准备输出..")
    cpic(newr,newg,newb,path,rd)
    print("{} - >> {}".format(path.split('/')[-1],"blurred.png"))
    print("总计耗时:{:.3f}s,感谢您的使用.".format(timer[0]+timer[1]+timer[2]+timer[3]))
main()