数字图像处理（一）之灰度转换和卷积python实现

使用Python实现数字图像处理中如下功能：

彩色图像转成灰度图像
实现图像的相关&卷积操作
实现图像的高斯核卷积

使用的库和python版本如下：

imageio：2.9.0 用于读取磁盘中的图片文件
numpy：1.20.3 用于矩阵等操作
matplotlib：3.4.2 用于画图
python：3.8.11

读取图像

在进行图像处理操作前，首先需要对图像进行读取。这里使用imageio库对图片进行读取，并将其转成numpy数组。

下面定义一个covert_img_to_array函数，用于读取图片。

def covert_img_to_array(self, path:str) -> np.array:

    """[将图片转成Array便于处理]

        Args:

            path (str): [图片保存位置]

        Returns:

            np.array: [返回numpy数组，数组元素uint8]

        """

    return np.array(imageio.imread(path))

展示图片

使用matplotlib库用于展示图片，为了更高的展示如片，定义下show_img函数，当不指定col或者row时尽量以方正的形式去展示图片。

def show_img(self,title:str, imgs:list, cmaps:list,row:int = 0,col:int = 0):

    """展示图片 len(imgs) must equal to the len of cmaps

    Args:

        title (str): [图像标题]

        imgs (list): [图片元组]

        cmaps (list): [mask,plt以何种形式展示图片，可参考官方文档使用：'gray'表示灰度图，None表示彩色图]

        row (int, optional): [指令row]. Defaults to 0.

        col (int, optional): [指令col]. Defaults to 0.

    """

    if len(imgs) != len(cmaps):

        print("图片和mask的len必须相同")

    else:

        if row == 0 and col !=0:

            row = np.ceil(len(imgs)/col).astype("uint8")

        elif row!=0 and col == 0:

            col = np.ceil(len(imgs)/row).astype("uint8")

        elif row*col < len(imgs):

            # 尽量以方正的形式去展示图片

            row = np.ceil(np.sqrt(len(imgs))).astype("uint8")

            col = np.ceil(len(imgs)/row).astype("uint8")

        for index,img in enumerate(imgs):

            plt.subplot(row,col,index+1)

            plt.imshow(img,cmap=cmaps[index])

        plt.suptitle(title)

        plt.show()

彩色图像转成灰度图像

彩色图像一般来说RGB表示的。也就是说，如果有一张64*64大小的图片，那么它在numpy中便是以64*64*3的shape进行保存的。将RGB图片转成灰度图有两种方式：

\(gray=\frac{R+G+B}{3}\)
\(gray=R*0.2989 + G*0.5870 + B*0.1140\) 这种灰度转换称之为NTSC标准，考虑了人类的彩色感知体验。

下面定义covert_rgb_to_gray函数，其中method如果为average，则使用第一种方式灰度转换方式；默认为NTSC，使用第二种方式转换。

def covert_rgb_to_gray(self, image:np.array, method:str = 'NTSC') -> np.array:

    """将RGB图像转成gray图像

    Args:

        image (np.array): [rgb图像]

        method (str, optional): [转换模式]. Defaults to 'NTSC'.

    Returns:

        Array: [返回的灰度图像]

    """

    if method == 'average':

        gray_img = image[:,:,0]/3+image[:,:,1]/3+image[:,:,2]/3

    else:

        gray_img = image[:,:,0]*0.2989 + image[:,:,1]*0.5870 + image[:,:,2]*0.1140

    return gray_img

图像卷积

图像卷积的公式如下所示，\(g\)代表输入的像素矩阵，\(w\)代表的是权重系数矩阵也就是所谓的卷积核kernel。

\[h(x,y) =\sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t)g(x-s,y-t)
\]

这里有一个很需要值得注意的点，那就是相关操作。相关操作和卷积很类似，相关操作的公式如下：

\[h(x,y) =\sum_{s=-a}^{a} \sum_{t=-b}^{b} w(s,t)g(x+s,y+t)
\]

在网络有一些博客文章，在解释卷积的时候，使用的是第一个公式，但是在做计算或者实现代码的时候却用的是第二个公式，这样做是不对的。因为卷积的kernel与相关的kernel相差了\(180^{\circ}\)。

但是值得注意的是，在卷积神经网络中，实际上使用的数学公式是相关相关运算，如下图所示。因为在CNN中，kernel的参数是学习过来的，kernel是否翻转并不会影响结果。

理解卷积

前置知识：

卷积定理指出，函数卷积的傅里叶变换是函数傅里叶变换的乘积。至于推导，可以查一下资料。

\[\mathcal{F}\{f * g\}=\mathcal{F}\{f\} \cdot \mathcal{F}\{g\}
\]

提一下图像卷积的含义。如果一个如下的均值滤波器对图像进行卷积，从人类的直觉进行出发，可以去除噪声和平滑图像。（在图像中，一般图像噪声的频率比较大，图像边缘部分的频率也比较大。因此使用均值滤波器可以去除噪声和平滑图像。）

\[1 / 9\left[\begin{array}{lll}
1 & 1 & 1 \\
1 & 1 & 1 \\
1 & 1 & 1
\end{array}\right]
\]

那么为什么会造成这种现象呢？如何从数学的角度来解释均值滤波器的作用呢？

如下所示，图左边是一个一维均值滤波器的函数图像，图右边是均值函数在频域上面的图像。在右边图像上，可以发现一个很明显的特点：频率越高，\(F(\mu)\)越小。

那么如果将\(F(\mu)\)与某另外一个频域上面的函数（比如图像）相乘，显而易见，如果图像的频率越高，则\(F(\mu)\)与之相乘被拖下水的的程度就越大。也就是说，相乘之后，频率低的就被抬上去了，频率高的被拉下去了。

说的细一点，其实从上图可以看到，随着频率的增大，\(F(\mu)\)并不是严格的下降，中间有一个波浪的起伏，这样会在边缘造成一些不好的现象。但是高斯滤波不会有这种情况。后面会介绍高斯滤波。

均值滤波器的二维频域图如下所示：

矩阵点积

下面定义矩阵点积函数。

def __matrix_dot_product(self,matrix,kernel):

    """矩阵点乘 [1,2,3]*[4,5,6] = 1*4 + 2*5 + 3*6 = 32

    Args:

        matrix ([type]): [部分图像]

        kernel ([type]): [kernel]

    Returns:

        [type]: [点乘结果]

    """

    if len(matrix) != len(kernel):

        print("点积失败，大小不一致")

    else:

        # 速度快

        return (np.multiply(matrix,kernel)).sum()

        # result = 0

        # for i, row_nums in enumerate(matrix):

        #     for j,num in enumerate(row_nums):

        #         result += num * kernel[i][j]

        # return result

图像padding

如果不对图像进行padding的话，会造成一个现象，图像越卷越小。在卷积的时候，我们希望卷积后的图像大小与原图像保持一致（CNN网络可能会越卷越小），因此需要对图像进行padding。padding有两种方式，一种在填充0，一种是填充与其距离最近的元素。下图中图像周围虚线部分就是padding的元素。

下面是实现padding操作的具体函数。实际上，可以直接使用np.pad操作实现。（但是我的作业要求不能使用pad操作，只能自己实现）

    def __padding(self, padding_type:str, image:np.array, padding_w:int, padding_h:int):

        """对图片进行padding

        Args:

            padding_type (str): [padding方式]

            image (np.array): [图片]

            padding_w (int): [宽度pdding]

            padding_h (int): [高度padding，一般来说padding_w = padding_h]

        Returns:

            [type]: [返回padding之后的结果]

        """

        image_w = image.shape[0]

        image_h = image.shape[1]

        padding_image = np.zeros((image_w+padding_w*2,image_h+padding_h*2))

        padding_image[padding_w:padding_w+image_w,padding_h:padding_h+image_h] = image

        if padding_type == 'zero':

            return padding_image

        if padding_type == "replicate":

            # 补充四个角

            padding_image[0:padding_w+1,0:padding_h+1] = image[0,0]

            padding_image[image_w+padding_w-1:,0:padding_h+1] = image[image_w-1,0]

            padding_image[0:padding_w+1,image_h+padding_h-1:] = image[0,image_h-1]

            padding_image[image_w+padding_w-1:,image_h+padding_h-1:] = image[image_w-1,image_h-1]

            # 补充旁边的元素

            for i in range(padding_w+1,image_w+padding_w-1):

                padding_image[i,0:padding_h] = image[i-padding_w,0]

                padding_image[i,image_h+padding_h:] = image[i-padding_w,image_h-1]

            for i in range(padding_h+1,image_h+padding_h-1):

                padding_image[0:padding_w,i] = image[0,i-padding_h]

                padding_image[image_w+padding_w:,i] = image[image_w-1,i-padding_h]

            return padding_image

如果想使得卷积之后的结果与原图像一致，padding_w,padding_h为卷积核大小的一半（向下取整，卷积核大小一般是奇数）。比如核的大小是\(5 \times 5\)，那么padding的长宽便是\(2\)。

图像相关操作

前面说过图像的卷积实际上就是将kernel进行翻转\(180^{\circ}\)，然后进行相关运算，因此可以先定义相关操作函数：

def corr2D(self, image:np.array, kernel:np.array, padding:str = 'zero') -> np.array:

    """对图片进行相关运算。

    Args:

        image (np.array): [(*,*)shape的图片]

        kernel (np.array): [kernel，kernel为奇数]

        padding (str, optional): [zero以零填充，replicate以邻近的填充]. Defaults to 'zero'.

    Returns:

        [type]: [description]

    """

    kernel_size_w = kernel.shape[0]

    kernel_size_h = kernel.shape[1]

    image_w,image_h = image.shape

    padding_w = kernel_size_w // 2

    padding_h = kernel_size_h // 2

    # 将图片padding起来

    padding_image = self.__padding(padding,image,padding_w,padding_h)

    new_image = np.zeros((image_w,image_h))

    for i in range(image_w):

        for j in range(image_h):

            new_image[i][j] = self.__matrix_dot_product(padding_image[i:i+kernel_size_w,j:j+kernel_size_h],kernel)

    return new_image.clip(0,255).astype("uint8")

卷积操作

旋转kernel

旋转kernel的代码很简单，如下所示，通过以下操作可以将行和列翻转（相当于反转了\(180^{\circ}\)）。

def flip_180(self, arr: np.array) -> np.array:

    return arr[::-1,::-1]

卷积

将kernel继续宁翻转，然后进行相关运算便是卷积了。

def conv2D(self, image:np.array, kernel:np.array, padding:str = 'zero') -> np.array:

    """二维卷积

    Args:

        image (np.array): [(*,*)shape的图片]

        kernel (np.array): [kernel，kernel为奇数]

        padding (str, optional): [zero以零填充，replicate以邻近的填充]. Defaults to 'zero'.

    Returns:

        [type]: [卷积好的结果]

    """

    return self.corr2D(image,self.flip_180(kernel),padding)

高斯核

二维高斯核的公式如下所示：

\[G(x, y,\sigma_x,\sigma_y)=\frac{1}{2 \pi \sigma_{x}\sigma_{y}} e^{-\left(\frac{x^{2}}{2{\sigma_x}^2} + \frac{y^{2}}{2{\sigma_y}^2}\right)}
\]

二维高斯核的频域图如下所示。

下面是二维高斯滤波函数的定义，其中\(\sigma_x=\sigma_y=sig\)。并对卷积核进行归一化，使得所有元素加起来和为1。

    def gauss_2d_kernel(self,sig,m=0):

        """产生高斯核

        Args:

            sig ([type]): [高斯核参数 sigx = sigy]

            m (int, optional): [高斯kernel的大小]. Defaults to 0. if m=0，then m = ceil(3*sig)*2 +1

        Returns:

            [type]: [m*m大小的高斯核]

        """

        fit_m = math.ceil(3 * sig)*2+1

        if m == 0:

            m = fit_m

        if m < fit_m:

            print("你的核的size应该大一点")

        # 中心点

        center = m //2

        kernel = np.zeros(shape=(m,m))

        for i in range(m):

            for j in range(m):

                kernel[i][j] = (1/(2*math.pi*sig**2))*math.e**(-((i-center)**2+(j-center)**2)/(2*sig**2))

		# 归一化

        return kernel/(kernel.sum())

巴特西