Image Processing with Numpy —— github

在python中進行圖像處理侵浸，我們有三個好工具：OpenCV, SciKit-Image 和 Pillow。但是在本文中氛谜，為了理解一些簡單圖像處理技術(shù)的基礎(chǔ)掏觉，我們將使用numpy。所以這也是練習(xí)numpy的良好教程值漫。

• 涵蓋知識：
  圖像讀取與裁剪履腋，通道分離，顏色變換（愛因斯坦求和約定）惭嚣，動態(tài)gif生成，灰階轉(zhuǎn)換
  圖像卷積悔政， 圖像分割（大津二值晚吞，KMeans像素聚類），圖像矢量化（輪廓提取谋国，內(nèi)邊緣判斷槽地，輪廓填充）

導(dǎo)入庫并加載圖像

import numpy as np
import matplotlib.pylab as plt
# 加載圖像
im = plt.imread("BTD.jpg") # 加載當(dāng)前文件夾中名為BTD.jpg的圖片
print(im.shape) # 輸出圖像尺寸
# (4608, 2592, 3)即(y軸像素點數(shù)， x軸像素點數(shù)，圖像通道數(shù))
# 這里用的是RGB三通道圖像捌蚊，通道數(shù)為3

裁剪圖像

def plti(im, **kwargs):
    """
    畫圖的輔助函數(shù)
    """
    plt.imshow(im, interpolation="none", **kwargs)
    plt.axis('off') # 去掉坐標(biāo)軸
    plt.show() # 彈窗顯示圖像

im = im[400:3800,:2000,:]  # 直接切片對圖像進行裁剪
plti(im)

分離各通道的圖像

fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(15,5)) 
# 將一張圖分為1x3個子圖集畅，axs為各子圖對象構(gòu)成的列表。figsize為顯示窗口的橫縱比缅糟。

for c, ax in zip(range(3), axs): # 使用zip來同時循環(huán)3通道和3個子圖對象
    tmp_im = np.zeros(im.shape) # 初始化一個和原圖像大小相同的三維數(shù)組
    # 注意 tmp_im 仍然是三通道
    tmp_im[:,:,c] = im[:,:,c] # 只復(fù)制某一通道
    one_channel = im[:,:,c].flatten() # 索引該通道并展平至一維
    print("channel", c, " max = ", max(one_channel), "min = ", min(one_channel)) # 輸出該通道最大最小的像素值
    ax.imshow(tmp_im) # 在子圖上繪制
    ax.set_axis_off() # 去掉子圖坐標(biāo)軸
# 注意以上 tmp_im 采用的是切片復(fù)制
plt.show()

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

顏色變換

• 一個RGB圖像挺智，其中每個點由RGB三通道的值組合為最終顏色，如果將RGB三個通道的值分別對應(yīng)到XYZ軸上窗宦，一個RGB顏色就表示為一個三維空間中的一個點赦颇，而對三維空間中的點我們可以做空間變換（即乘一個三維矩陣），通過空間變換赴涵，一個顏色就變?yōu)?strong>另一個顏色（這之間需要做一些歸一化和反歸一的變換）媒怯。讓我們嘗試下對一張圖上的所有點（顏色）都做相同的變換，看看會有什么效果髓窜。
  （以下作為代碼看不明白時的補充）
• Sigmoid函數(shù)
  
• 愛因斯坦求和約定（Einstein notation）
• Numpy中的愛因斯坦求和約定（np.einsum）-Einstein Summation in Numpy(非常好的解釋)

# assert( np.log(np.e) == 1.0)
# np.log 即 ln() - 以e為底的對數(shù)函數(shù)

def do_normalise(im):
    return -np.log(1/((1 + im)/257) - 1)
# 預(yù)處理函數(shù)
# im中的像素值為 [0, 255] 閉區(qū)間扇苞， 則 (1+im) 為 [1, 256]
# 先做 (1+im)/257 操作將值歸一化到 (0, 1) 開區(qū)間內(nèi)
# 再使用 sigmoid函數(shù) 的反函數(shù)，效果見sigmod函數(shù)圖像
# -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")
# 預(yù)處理函數(shù)的反函數(shù)
# 即先使用sigmod函數(shù)寄纵，再將值變換到(0, 257)區(qū)間再減1鳖敷，通過astype保證值位于[0, 255]
# 關(guān)于 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 旋轉(zhuǎn)矩陣，圍繞X軸旋轉(zhuǎn)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軸（即紅色通道軸）旋轉(zhuǎn)180°哄陶。
im2 = undo_normalise(im_rotated)

plti(im2)

嘗試動態(tài)效果

• 現(xiàn)在讓我們連續(xù)地旋轉(zhuǎn)這些像素點，看看效果如何哺壶。
• 我們將利用matplotlib庫中的 FuncAnimation工具
• 介紹FuncAnimation的優(yōu)秀博客
• 為了將圖片保存為GIF（matplotlib必須要外部支持才能將動圖保存為gif）屋吨，我們需要下載imagemagick工具（我用的是win7平臺，下載的是ImageMagick-7.0.3-4-Q16-x64-static.exe版本）并安裝山宾，我直接默認(rèn)安裝在了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)
    # 在更新函數(shù)里根據(jù)i來改變旋轉(zhuǎn)的角度
    ax.imshow(im2)
    ax.set_title("Angle: {}*pi/10".format(i), fontsize=20)
    ax.set_axis_off()
    # 將旋轉(zhuǎn)后的圖繪出
# 以上其余注釋見前文

anim = FuncAnimation(fig, update, frames=np.arange(0, 20), interval=50)
# fig是圖像句柄
# update是更新函數(shù)
# frames為幀數(shù)列表绷杜，將值依次提供給更新函數(shù)
# interval表示每幀間隔ms數(shù)
anim.save('colour_rotation.gif', dpi=80, writer='imagemagick')
plt.close()

哈哈直秆，是不是有點迷幻呢。

灰階

• np.tile(A, reps)：按照reps指定的次數(shù)在相應(yīng)的維度上重復(fù)A矩陣來構(gòu)建一個新的矩陣鞭盟。
• np.array的*操作符不像其他大多數(shù)語言圾结，是elementwise的，即兩個矩陣中對應(yīng)元素逐一互乘齿诉，而矩陣乘法要使用dot()方法筝野。

def to_grayscale(im, weights = np.c_[0.2989, 0.5870, 0.1140]):
    """
    取原始圖像的RGB值的加權(quán)平均來將圖片轉(zhuǎn)換為灰階晌姚，權(quán)重矩陣為tile
    """
    # 默認(rèn)的 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是一個窗口函數(shù)号坡。本質(zhì)上就是將每個像素點的值替換為它與其相鄰元素點的值的加權(quán)和懊烤。理解公式的話，考慮W是一個固定大小的矩陣宽堆，對于每個確定的(x,y)腌紧，將W中的每個元素與I中相應(yīng)位置的元素相乘并求和就得到該(x,y)點處卷積后的值。更直觀的理解參看卷積神經(jīng)網(wǎng)絡(luò)CNN基本概念筆記畜隶。
• 由于卷積操作消耗較大壁肋，我們先將圖片縮小再對圖像運用一個均勻窗口，它能夠模糊圖像籽慢，其實就是以每個像素點和其相鄰像素點的均值來替代該點原來的值浸遗。

from scipy.ndimage.interpolation import zoom
im_small = zoom(im, (0.2,0.2,1))
# zoom 將圖片每一維以相應(yīng)系數(shù)縮小
# im.shape = (3400, 2000, 3)
# im_small.shape = (680, 400, 3)

from scipy.signal import convolve2d
# 引入二維卷積函數(shù)
def convolve_all_colours(im, window):
    """
    用窗口window卷積圖像，依次對圖像的每個通道卷積
    """
    ims = []
    # 用ims作為每個通道轉(zhuǎn)換結(jié)果的暫存列表
    for d in range(3):
    # 對圖像的三個通道循環(huán)處理
        im_conv_d = convolve2d(im[:,:,d], window, mode="same", boundary="symm")
        # mode決定輸出尺寸箱亿，boundary決定邊界條件跛锌，這里輸出尺寸與原圖相同，采用對稱邊界條件
        ims.append(im_conv_d)
        # 將單通道轉(zhuǎn)換結(jié)果添加到列表

    im_conv = np.stack(ims, axis=2).astype("uint8")
    # 在第三維上堆疊ims列表中的每個元素届惋，并通過astype保證值在0-255
    return im_conv

n=50
window = np.ones((n,n))
# 構(gòu)建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):
    """
    使用高斯分布的權(quán)重創(chuàng)建一個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可以將輸入轉(zhuǎn)化為np.array, 這里輸入為一個列表推導(dǎo)式
    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做圖像處理（下）