用numpy做图像处理（上）

Image Processing with Numpy —— github

在python中进行图像处理，我们有三个好工具：OpenCV, SciKit-Image 和 Pillow。但是在本文中，为了理解一些简单图像处理技术的基础，我们将使用numpy。所以这也是练习numpy的良好教程。

• 涵盖知识：
  图像读取与裁剪，通道分离，颜色变换（爱因斯坦求和约定），动态gif生成，灰阶转换
  图像卷积， 图像分割（大津二值，KMeans像素聚类），图像矢量化（轮廓提取，内边缘判断，轮廓填充）

---

导入库并加载图像

import numpy as np
import matplotlib.pylab as plt
# 加载图像
im = plt.imread("BTD.jpg") # 加载当前文件夹中名为BTD.jpg的图片
print(im.shape) # 输出图像尺寸
# (4608, 2592, 3)即(y轴像素点数， x轴像素点数，图像通道数)
# 这里用的是RGB三通道图像，通道数为3

裁剪图像

def plti(im, **kwargs):
    """
    画图的辅助函数
    """
    plt.imshow(im, interpolation="none", **kwargs)
    plt.axis('off') # 去掉坐标轴
    plt.show() # 弹窗显示图像

im = im[400:3800,:2000,:]  # 直接切片对图像进行裁剪
plti(im)

分离各通道的图像

fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(15,5)) 
# 将一张图分为1x3个子图，axs为各子图对象构成的列表。figsize为显示窗口的横纵比。

for c, ax in zip(range(3), axs): # 使用zip来同时循环3通道和3个子图对象
    tmp_im = np.zeros(im.shape) # 初始化一个和原图像大小相同的三维数组
    # 注意 tmp_im 仍然是三通道
    tmp_im[:,:,c] = im[:,:,c] # 只复制某一通道
    one_channel = im[:,:,c].flatten() # 索引该通道并展平至一维
    print("channel", c, " max = ", max(one_channel), "min = ", min(one_channel)) # 输出该通道最大最小的像素值
    ax.imshow(tmp_im) # 在子图上绘制
    ax.set_axis_off() # 去掉子图坐标轴
# 注意以上 tmp_im 采用的是切片复制
plt.show()

#输出：
#channel 0  max =  220 min =  11
#channel 1  max =  203 min =  10
#channel 2  max =  185 min =  0

颜色变换

• 一个RGB图像，其中每个点由RGB三通道的值组合为最终颜色，如果将RGB三个通道的值分别对应到XYZ轴上，一个RGB颜色就表示为一个三维空间中的一个点，而对三维空间中的点我们可以做空间变换（即乘一个三维矩阵），通过空间变换，一个颜色就变为另一个颜色（这之间需要做一些归一化和反归一的变换）。让我们尝试下对一张图上的所有点（颜色）都做相同的变换，看看会有什么效果。
  （以下作为代码看不明白时的补充）
• Sigmoid函数
  
• 爱因斯坦求和约定（Einstein notation）
• Numpy中的爱因斯坦求和约定（np.einsum）-Einstein Summation in Numpy(非常好的解释)

# assert( np.log(np.e) == 1.0)
# np.log 即 ln() - 以e为底的对数函数

def do_normalise(im):
    return -np.log(1/((1 + im)/257) - 1)
# 预处理函数
# im中的像素值为 [0, 255] 闭区间， 则 (1+im) 为 [1, 256]
# 先做 (1+im)/257 操作将值归一化到 (0, 1) 开区间内
# 再使用 sigmoid函数 的反函数，效果见sigmod函数图像
# -np.log(1/((1 + 0)/257) - 1) = -5.5451774444795623
# -np.log(1/((1 + 255)/257) - 1) = 5.5451774444795623

def undo_normalise(im):
    return (1/(np.exp(-im) + 1) * 257 - 1).astype("uint8")
# 预处理函数的反函数
# 即先使用sigmod函数，再将值变换到(0, 257)区间再减1，通过astype保证值位于[0, 255]
# 关于 astype("uint8") ：
# np.array([-1]).astype("uint8") = array([255], dtype=uint8)
# np.array([256]).astype("uint8") = array([0], dtype=uint8)

def rotation_matrix(theta):
    """
    3D 旋转矩阵，围绕X轴旋转theta角
    """
    return np.c_[
        [1,0,0],
        [0,np.cos(theta),-np.sin(theta)],
        [0,np.sin(theta),np.cos(theta)]
    ]
# np.c_[ ] 将列表中的元素在第二维上拼接起来
# np.c_[[1,2],[3,4],[5,6]] =
# array([[1, 3, 5],
#        [2, 4, 6]])


im_normed = do_normalise(im)
im_rotated = np.einsum("ijk,lk->ijl", im_normed, rotation_matrix(np.pi))
# 利用爱因斯坦求和约定做矩阵乘法，实际上是将每个RGB像素点表示的三维空间点绕X轴（即红色通道轴）旋转180°。
im2 = undo_normalise(im_rotated)

plti(im2)

尝试动态效果

• 现在让我们连续地旋转这些像素点，看看效果如何。
• 我们将利用matplotlib库中的 FuncAnimation工具
• 介绍FuncAnimation的优秀博客
• 为了将图片保存为GIF（matplotlib必须要外部支持才能将动图保存为gif），我们需要下载imagemagick工具（我用的是win7平台，下载的是ImageMagick-7.0.3-4-Q16-x64-static.exe版本）并安装，我直接默认安装在了C盘。
• 为了顺利使用，我们还要做点配置。先找到matplotlib配置文件路径：

import matplotlib
print(matplotlib.matplotlib_fname())
# 我的输出 C:\Anaconda3\lib\site-packages\matplotlib\mpl-data\matplotlibrc

得到文件路径后编辑该文件，在末尾添加一行（冒号后面为你的magick.exe工具的安装路径）

animation.convert_path: C:\Program Files\ImageMagick-7.0.3-Q16\magick.exe

这样就OK了。

from matplotlib.animation import FuncAnimation

fig, ax = plt.subplots(figsize=(5,8))

def update(i):
    im_normed = do_normalise(im)
    im_rotated = np.einsum("ijk,lk->ijl", im_normed, rotation_matrix(i * np.pi/10))
    im2 = undo_normalise(im_rotated)
    # 在更新函数里根据i来改变旋转的角度
    ax.imshow(im2)
    ax.set_title("Angle: {}*pi/10".format(i), fontsize=20)
    ax.set_axis_off()
    # 将旋转后的图绘出
# 以上其余注释见前文

anim = FuncAnimation(fig, update, frames=np.arange(0, 20), interval=50)
# fig是图像句柄
# update是更新函数
# frames为帧数列表，将值依次提供给更新函数
# interval表示每帧间隔ms数
anim.save('colour_rotation.gif', dpi=80, writer='imagemagick')
plt.close()

哈哈，是不是有点迷幻呢。

灰阶

• np.tile(A, reps)：按照reps指定的次数在相应的维度上重复A矩阵来构建一个新的矩阵。
• np.array的*操作符不像其他大多数语言，是elementwise的，即两个矩阵中对应元素逐一互乘，而矩阵乘法要使用dot()方法。

def to_grayscale(im, weights = np.c_[0.2989, 0.5870, 0.1140]):
    """
    取原始图像的RGB值的加权平均来将图片转换为灰阶，权重矩阵为tile
    """
    # 默认的 weights = array([[ 0.2989,  0.587 ,  0.114 ]])
    tile = np.tile(weights, reps=(im.shape[0],im.shape[1],1))
    # assert( tile.shape == im.shape )
    return np.sum(tile * im, axis=2)
    # np.sum意味着沿某一轴求和，axis=2为第三维（0为第一维）
    # 整个乘法意味着由图像每个像素点的RGB 得到 (R*0.2989+ G*0.5870+ B*0.1140)灰阶值，图像的二维尺寸不变，而减为单通道。
img = to_grayscale(im)
plti(img, cmap='Greys') # 注意要以灰度形式画出

1.png

卷积

• 卷积是图像处理的基本操作，它的公式是
  C(x,y)=∫dx′dy′I(x+x′,y+y′)W(x′,y′)
  
• C是卷积后的图像，I是原始图像，W是一个窗口函数。本质上就是将每个像素点的值替换为它与其相邻元素点的值的加权和。理解公式的话，考虑W是一个固定大小的矩阵，对于每个确定的(x,y)，将W中的每个元素与I中相应位置的元素相乘并求和就得到该(x,y)点处卷积后的值。更直观的理解参看卷积神经网络CNN基本概念笔记。
• 由于卷积操作消耗较大，我们先将图片缩小再对图像运用一个均匀窗口，它能够模糊图像，其实就是以每个像素点和其相邻像素点的均值来替代该点原来的值。

from scipy.ndimage.interpolation import zoom
im_small = zoom(im, (0.2,0.2,1))
# zoom 将图片每一维以相应系数缩小
# im.shape = (3400, 2000, 3)
# im_small.shape = (680, 400, 3)

from scipy.signal import convolve2d
# 引入二维卷积函数
def convolve_all_colours(im, window):
    """
    用窗口window卷积图像，依次对图像的每个通道卷积
    """
    ims = []
    # 用ims作为每个通道转换结果的暂存列表
    for d in range(3):
    # 对图像的三个通道循环处理
        im_conv_d = convolve2d(im[:,:,d], window, mode="same", boundary="symm")
        # mode决定输出尺寸，boundary决定边界条件，这里输出尺寸与原图相同，采用对称边界条件
        ims.append(im_conv_d)
        # 将单通道转换结果添加到列表

    im_conv = np.stack(ims, axis=2).astype("uint8")
    # 在第三维上堆叠ims列表中的每个元素，并通过astype保证值在0-255
    return im_conv

n=50
window = np.ones((n,n))
# 构建50x50的全1矩阵
window /= np.sum(window)
# 矩阵每个元素除以矩阵所有元素的和，使矩阵所有元素的和为1
plti(convolve_all_colours(im_small, window))

• 实际上，模糊图像有许多不同的方法。最常用的就是均匀窗口、高斯窗口、中值滤波。为了对它们的处理效果有个直观的感受，让我们以不同的窗口大小使用它们。

from scipy.ndimage import median_filter

def make_guassian_window(n, sigma=1):
    """
    使用高斯分布的权重创建一个n*n的方形窗口
    """
    nn = int((n-1)/2)
    a = np.asarray([[x**2 + y**2 for x in range(-nn,nn+1)] for y in range(-nn,nn+1)])
    # np.asarray可以将输入转化为np.array, 这里输入为一个列表推导式
    return np.exp(-a/(2*sigma**2))

def median_filter_all_colours(im_small, window_size):
    """
    对图像所有通道运用中值滤波
    """
    ims = []
    for d in range(3):
        im_conv_d = median_filter(im_small[:,:,d], size=(window_size,window_size))
        ims.append(im_conv_d)

    im_conv = np.stack(ims, axis=2).astype("uint8")
    
    return im_conv

window_sizes = [9,17,33,65]
fig, axs = plt.subplots(nrows=3, ncols=len(window_sizes), figsize=(15,15));

# 均值滤波 - 均匀窗口
for w, ax in zip(window_sizes, axs[0]):
    window = np.ones((w,w))
    window /= np.sum(window)
    ax.imshow(convolve_all_colours(im_small, window));
    ax.set_title("Mean Filter: window size: {}".format(w));
    ax.set_axis_off();
    
# 高斯滤波 - 高斯窗口
for w, ax in zip(window_sizes, axs[1]):
    window = make_guassian_window(w,sigma=w)
    window /= np.sum(window)
    ax.imshow(convolve_all_colours(im_small, window));
    ax.set_title("Guassian Filter: window size: {}".format(w));
    ax.set_axis_off();
    
# 中值滤波
for w, ax in zip(window_sizes, axs[2]):
    ax.imshow(median_filter_all_colours(im_small, w));
    ax.set_title("Median Filter: window size: {}".format(w));
    ax.set_axis_off();

用numpy做图像处理（下）