Image Processing with Numpy —— github

請先看用numpy做圖像處理（上）

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

卷積（下）

• 讓我們再使用一個新的窗口濾波器：Sobel濾波器，它被用來近似圖像在某一方向上的梯度（也就是邊緣檢測），它的窗口形式是這樣的：

  
  
  
  

  讓我們對X和Y方向都使用一次該窗口，然后對每個像素點的兩次處理結(jié)果求平方和再開根號作為該點的值

n=100
sobel_x = np.c_[
    [-1,0,1],
    [-2,0,2],
    [-1,0,1]
]

sobel_y = np.c_[
    [1,2,1],
    [0,0,0],
    [-1,-2,-1]
]
# 或者使用sobel_y = sobel_x.transpose()

ims = []
for d in range(3):
    sx = convolve2d(im_small[:,:,d], sobel_x, mode="same", boundary="symm")
    sy = convolve2d(im_small[:,:,d], sobel_y, mode="same", boundary="symm")
    ims.append(np.sqrt(sx*sx + sy*sy)) # 平方和再開根號，疊加兩次處理效果

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

plti(im_conv)

• 我們再試試組合中值濾波和sobel濾波的效果

im_smoothed = median_filter_all_colours(im_small, 71)
# 先利用之前設(shè)置的函數(shù)進行中值濾波，再對圖像進行sobel濾波
sobel_x = np.c_[
    [-1,0,1],
    [-2,0,2],
    [-1,0,1]
]

sobel_y = np.c_[
    [1,2,1],
    [0,0,0],
    [-1,-2,-1]
]

ims = []
for d in range(3):
    sx = convolve2d(im_smoothed[:,:,d], sobel_x, mode="same", boundary="symm")
    sy = convolve2d(im_smoothed[:,:,d], sobel_y, mode="same", boundary="symm")
    ims.append(np.sqrt(sx*sx + sy*sy))

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

plti(im_conv)

• 到目前為止，我們都是對圖像的所有通道一起操作，讓我們看看如果一次只對圖像的某一通道進行模糊，會有什么效果

fig, axs = plt.subplots(nrows=1, ncols=4, figsize=(15,5))
# 得到1x4的畫圖結(jié)構(gòu)
ax= axs[0]
ax.imshow(im_small)
ax.set_title("Original", fontsize=20)
ax.set_axis_off()
# 畫原圖
w=75
window = np.ones((w,w))
window /= np.sum(window)
# 設(shè)置均值模糊窗口
ax= axs[1]
ims = []
for d in range(3):
    if d == 0:
        im_conv_d = convolve2d(im_small[:,:,d],window, mode="same", boundary="symm")
    else:
        im_conv_d = im_small[:,:,d]
    ims.append(im_conv_d)
ax.imshow(np.stack(ims, axis=2).astype("uint8"))
ax.set_title("Red Blur", fontsize=20)
ax.set_axis_off()
# 只作用于R通道
ax= axs[2]
ims = []
for d in range(3):
    if d == 1:
        im_conv_d = convolve2d(im_small[:,:,d], window, mode="same", boundary="symm")
    else:
        im_conv_d = im_small[:,:,d]
    ims.append(im_conv_d)
ax.imshow(np.stack(ims, axis=2).astype("uint8"))
ax.set_title("Blue Blur", fontsize=20)
ax.set_axis_off()
# 只作用于G通道
ax= axs[3]
ims = []
for d in range(3):
    if d == 2:
        im_conv_d = convolve2d(im_small[:,:,d], window, mode="same", boundary="symm")
    else:
        im_conv_d = im_small[:,:,d]
    ims.append(im_conv_d)
ax.imshow(np.stack(ims, axis=2).astype("uint8"))
ax.set_title("Green Blur", fontsize=20)
ax.set_axis_off()
# 只作用于B通道

分割

• 圖像處理的另一個主要領(lǐng)域分割圖像成不同的區(qū)域，比如前景和背景。
• 最簡單的方式是將圖片轉(zhuǎn)換為灰階，然后找到一個閾值，以像素值與閾值的比較判斷其為前景還是背景。我們先隨便設(shè)一組閾值，看看不同大小的閾值的效果區(qū)別

def simple_threshold(im, threshold=128):
    return ((im > threshold) * 255).astype("uint8")
# 以閾值對灰度圖像做二值處理
# np.array([1,2,3,4,5,6]) > 3 得到
# array([False, False, False,  True,  True,  True], dtype=bool)
# assert( True * 255 == 255 )

thresholds = [100,120,128,138,150]
# 預(yù)設(shè)的閾值列表
fig, axs = plt.subplots(nrows=1, ncols=len(thresholds), figsize=(20,5));
gray_im = to_grayscale(im)
# 先將圖片轉(zhuǎn)換為灰階
for t, ax in zip(thresholds, axs):
    # 循環(huán)閾值列表，依次處理畫圖
    ax.imshow(simple_threshold(gray_im, t), cmap='Greys');
    ax.set_title("Threshold: {}".format(t), fontsize=20);
    ax.set_axis_off();

• 隨機選取閾值具有不確定性，而我們會希望在分割圖片時，使得類內(nèi)像素點值的差距最小，而類間差距最大，一種方式就是使用大津二值化算法

def otsu_threshold(im):

    pixel_counts = [np.sum(im == i) for i in range(256)]
    # 得到圖片的以0-255索引的像素值個數(shù)列表
    s_max = (0,-10)
    ss = []
    for threshold in range(256):
    # 遍歷所有閾值，根據(jù)公式挑選出最好的

        # 更新
        w_0 = sum(pixel_counts[:threshold]) # 得到閾值以下像素個數(shù)
        w_1 = sum(pixel_counts[threshold:]) # 得到閾值以上像素個數(shù)

        mu_0 = sum([i * pixel_counts[i] for i in range(0,threshold)]) / w_0 if w_0 > 0 else 0
        # 得到閾值下所有像素的均值
        # 注意 if else 用法意義： 如果 w_0 > 0 則 mu_0 = sum/w_0 否則 mu_0 = 0
        mu_1 = sum([i * pixel_counts[i] for i in range(threshold, 256)]) / w_1 if w_1 > 0 else 0
        # 得到閾值上所有像素的均值

        # 根據(jù)公式計算
        s = 1.0 * w_0 * w_1 * (mu_0 - mu_1) ** 2
        # 直接使用w_0 * w_1可能會造成整數(shù)相乘溢出，所以先乘一個1.0變?yōu)楦↑c數(shù)
        ss.append(s)

        # 取最大的
        if s > s_max[1]:
            s_max = (threshold, s)

    return s_max[0]

t = otsu_threshold(gray_im)
# 先根據(jù)大津算法算出閾值t，再使用閾值t來分割圖像
plti(simple_threshold(gray_im, t), cmap='Greys')

• 但是處理效果并不理想，很可能是在將圖片轉(zhuǎn)化為灰階的時候我們丟失了一些圖像信息，讓我們對圖像的每個通道分別處理試試

fig, axs = plt.subplots(nrows=1, ncols=3, figsize=(15,5))

c_ims = []
# c_ims收集各通道處理結(jié)果以待進一步處理

for c, ax in zip(range(3), axs): # 依次處理每個通道，畫到相應(yīng)的子圖上
    tmp_im = im[:,:,c]
    t = otsu_threshold(tmp_im)
    tmp_im = simple_threshold(tmp_im, t)
    ax.imshow(tmp_im, cmap='Greys')
    c_ims.append(tmp_im)
    ax.set_axis_off()
    

• 將三個通道的處理結(jié)果合并

# 黑色值為255，白色值為0
# 255 & 255 = 255; 0 & 0 = 0; 0 & 255 = 0
plti(c_ims[0] & c_ims[1] & c_ims[2], cmap='Greys')

• 這樣就好多了
• 如果不是尋找一個門檻，而是發(fā)現(xiàn)色彩的集群，我們可以使用K-means聚類技術(shù)，可以結(jié)合scikit.KMeans的最簡例子理解：
  

from sklearn.cluster import KMeans

h,w = im_small.shape[:2]
# 獲取圖像的高和寬
im_small_long = im_small.reshape((h * w, 3))
# 將圖像轉(zhuǎn)化為 RGB三值表示的三維空間坐標(biāo)點 的列表
im_small_wide = im_small_long.reshape((h,w,3))
# 可以將im_small_long恢復(fù)為原尺寸

km = KMeans(n_clusters=3)
# 分3個聚類
km.fit(im_small_long)
# 對所有RGB三維空間點分類，每個點獲得一個分類標(biāo)簽

cc = km.cluster_centers_.astype(np.uint8)
# 得到每個聚類的中心點坐標(biāo)，并保證在np.uint8范圍內(nèi)
out = np.asarray([cc[i] for i in km.labels_]).reshape((h,w,3))
# 以每一類的中心點坐標(biāo)所映射的RGB色彩 替代該類中所有點的色彩，并轉(zhuǎn)換為原始尺寸
# km.labels_是所有空間點的標(biāo)簽列表，列表中元素的值（即標(biāo)簽）為0,1或2，表示3個不同類

plti(out)

• 看上去還不錯，我們試試看如果隨機給三個類賦予色彩，會不會有意外的藝術(shù)效果：

rng = range(4)

fig, axs = plt.subplots(nrows=1, ncols=len(rng), figsize=(15,5))
                        
for t, ax in zip(rng, axs):
    rnd_cc = np.random.randint(0,256, size = (3,3))
    # 隨機色彩，RGB為1x3, 3個色彩就是3x3的矩陣
    out = np.asarray([rnd_cc[i] for i in km.labels_]).reshape((h,w,3))
    ax.imshow(out)
    ax.set_axis_off()

矢量化

• 到現(xiàn)在為止，我們都是將圖片作位圖處理。最后，我們將把狗的部分從圖像中取出，將它轉(zhuǎn)換為矢量圖
• 我們先利用聚類的結(jié)果將狗的部分提取出來

dog = np.asarray([cc[i] if i == 1 else [0,0,0] for i in km.labels_]).reshape((h,w,3))

plti(dog)

• 現(xiàn)在，我們可以使用來自Scikit-Image的輪廓跟蹤算法（measure.find_contours）來勾勒它

from skimage import measure

seg = np.asarray([(1 if i == 1 else 0)
                  for i in km.labels_]).reshape((h,w))
# 由聚類結(jié)果得到0,1二值單通道圖像

contours = measure.find_contours(seg, 0.5, fully_connected="high")
# 找到輪廓，返回每個連續(xù)輪廓的列表，每個連續(xù)輪廓是一組坐標(biāo)點

simplified_contours = [measure.approximate_polygon(c, tolerance=5) for c in contours]
# 簡化輪廓，近似為多邊形，tolerance決定精度

plt.figure(figsize=(5,10))

for n, contour in enumerate(simplified_contours):
    plt.plot(contour[:, 1], contour[:, 0], linewidth=2)
    # 連點成線，循環(huán)畫出每個輪廓

plt.ylim(h,0) # 設(shè)置y軸顯示范圍
plt.axes().set_aspect('equal') # 設(shè)置x,y軸對數(shù)據(jù)的比例相同

• 要畫出它，我們將輪廓線變?yōu)閳D塊補丁并填充

from matplotlib.patches import PathPatch
from matplotlib.path import Path

# 先翻轉(zhuǎn)輪廓的坐標(biāo)（交換x,y的位置），使其和matplotlib一致
paths = [[[v[1], v[0]] for v in vs] for vs in simplified_contours]

def get_path(vs):
    codes = [Path.MOVETO] + [Path.LINETO] * (len(vs)-2) + [Path.CLOSEPOLY]
    return PathPatch(Path(vs, codes))
    # 返回當(dāng)前輪廓繞出的圖塊補丁對象
    # [Path.MOVETO] + [Path.LINETO] + [Path.CLOSEPOLY] = [1, 2, 79]
    # [Path.LINETO] *3 = [2, 2, 2]

plt.figure(figsize=(5,10))

ax = plt.gca()
# gca獲取當(dāng)前坐標(biāo)軸

for n, path in enumerate(paths):
    ax.add_patch(get_path(path))

plt.ylim(h,0)
plt.xlim(0,w)

plt.axes().set_aspect('equal')

• 我們碰到一個問題，有些輪廓不是形狀的外邊緣，而是內(nèi)邊緣。我們來解決這個問題。
• 我們使用ray casting method來考察一個點是在多邊形內(nèi)還是多邊形外，并且基于我們現(xiàn)在處理的例子中，多邊形之間沒有交疊的狀況：一個多邊形要么是在另一個多邊形里面，要么是在它外面。最終subsume函數(shù)將所有多邊形組合進我們的多邊形對象中，它描述了外邊緣和可能的內(nèi)邊。

class Polygon:
    def __init__(self, outer, inners):
        self.outer = outer
        self.inners = inners
        
    def clone(self):
        return Polygon(self.outer[:], [c.clone() for c in self.inners])

def crosses_line(point, line):
    """
    Checks if the line project in the positive x direction from point
    crosses the line between the two points in line
    """
    p1,p2 = line
    
    x,y = point
    x1,y1 = p1
    x2,y2 = p2
    
    if x <= min(x1,x2) or x >= max(x1,x2):
        return False
    
    dx = x1 - x2
    dy = y1 - y2
    
    if dx == 0:
        return True
    
    m = dy / dx

    if y1 + m * (x - x1) <= y:
        return False
    
    return True

def point_within_polygon(point, polygon):
    """
    Returns true if point within polygon
    """
    crossings = 0
    
    for i in range(len(polygon.outer)-1):
        line = (polygon.outer[i], polygon.outer[i+1])
        crossings += crosses_line(point, line)
        
    if crossings % 2 == 0:
        return False
    
    return True

def polygon_inside_polygon(polygon1, polygon2):
    """
    Returns true if the first point of polygon1 is inside 
    polygon 2
    """
    return point_within_polygon(polygon1.outer[0], polygon2)


    
def subsume(original_list_of_polygons):
    """
    Takes a list of polygons and returns polygons with insides.
    
    This function makes the unreasonable assumption that polygons can
    only be complete contained within other polygons (e.g. not overlaps).
    
    In the case where polygons are nested more then three levels, 
    """
    list_of_polygons = [p.clone() for p in original_list_of_polygons]
    polygons_outside = [p for p in list_of_polygons
                        if all(not polygon_inside_polygon(p, p2) 
                               for p2 in list_of_polygons 
                               if p2 != p)]
    
    polygons_inside = [p for p in list_of_polygons
                        if any(polygon_inside_polygon(p, p2) 
                               for p2 in list_of_polygons
                               if p2 != p)]
    
    for outer_polygon in polygons_outside:
        for inner_polygon in polygons_inside:
            if polygon_inside_polygon(inner_polygon, outer_polygon):
                outer_polygon.inners.append(inner_polygon)
                
    return polygons_outside
                
    for p in polygons_outside:
        p.inners = subsume(p.inners)
        
    return polygons_outside
   
   
polygons = [Polygon(p,[]) for p in paths]
subsumed_polygons = subsume(polygons)


def build_path_patch(polygon, **kwargs):
    """
    Builds a matplotlb patch to plot a complex polygon described
    by the Polygon class.
    """
    verts = polygon.outer[:]
    codes = [Path.MOVETO] + [Path.LINETO] * (len(verts)-2) + [Path.CLOSEPOLY]

    for inside_path in polygon.inners:
        verts_tmp = inside_path.outer
        codes_tmp = [Path.MOVETO] + [Path.LINETO] * (len(verts_tmp)-2) + [Path.CLOSEPOLY]

        verts += verts_tmp
        codes += codes_tmp

    drawn_path = Path(verts, codes)

    return PathPatch(drawn_path, **kwargs)

    
def rnd_color():
    """Returns a randon (R,G,B) tuple"""
    return tuple(np.random.random(size=3))

plt.figure(figsize=(5,10))
ax = plt.gca()  
c = rnd_color()

for n, polygon in enumerate(subsumed_polygons):
    ax.add_patch(build_path_patch(polygon, lw=0, facecolor=c))
    
plt.ylim(h,0)
plt.xlim(0,w)
ax.set_axis_bgcolor('k')
plt.axes().set_aspect('equal')
plt.show()

• 到現(xiàn)在，我們只是做了圖像處理的最基本的一些東西，但這已經(jīng)不少了。如果要進一步學(xué)習(xí)，可以參考一些圖片庫所提供的例子。