機器學習算法實踐—K-Means算法與圖像分割

一、理論準備
1.1、圖像分割
圖像分割是圖像處理中的一種方法，圖像分割是指將一幅圖像分解成若干互不相交區域的集合，其實質可以看成是一種像素的聚類過程。通常使用到的圖像分割的方法可以分為：
基于邊緣的技術
基于區域的技術
基于聚類算法的圖像分割屬于基于區域的技術。
1.2、K-Means算法
K-Means算法是基于距離相似性的聚類算法，通過比較樣本之間的相似性，將形式的樣本劃分到同一個類別中，K-Means算法的基本過程為：
初始化常數 ，隨機初始化k個聚類中心
重復計算以下過程，直到聚類中心不再改變
計算每個樣本與每個聚類中心之間的相似度，將樣本劃分到最相似的類別中
計算劃分到每個類別中的所有樣本特征的均值，并將該均值作為每個類新的聚類中心
輸出最終的聚類中心以及每個樣本所屬的類別
在K-Means算法中，需要隨機初始化k個聚類中心，而K-Means算法對初始聚類中心的選取較為敏感，若選擇的聚類中心不好，則得到的聚類結果會非常差，因此，對K-Means算法提出了很多的改進的方法，如K-Means++算法，在K-Means++算法中，希望初始化的k個聚類中心之間的距離盡可能的大，其具體過程為：

在數據集中隨機選擇一個樣本點作為第一個初始化的聚類中心
選擇出其余的聚類中心：
計算樣本中的每一個樣本點與已經初始化的聚類中心之間的距離，并選擇其中最短的距離
以概率選擇距離最大的樣本作為新的聚類中心，重復上述過程，直到 個聚類中心都被確定
對k個初始化的聚類中心，利用K-Means算法計算最終的聚類中心。
對于K-Means算法的具體過程可以參考博文簡單易學的機器學習算法——kMeans，K-Means++算法的具體過程稍后會補充。

二、實踐準備

實踐中使用Python作為開發語言，使用到的模塊包括numpy和Image。numpy模塊是python中矩陣計算使用最多的模塊。

Image模塊是PIL(Python Imaging Library)中的模塊，對于Image模塊，主要是對圖像的一些操作：

模塊的頭文件
import Image as image

打開圖片
fp = open("003.JPG", "rb")
im = image.open(fp)

首先是以二進制文件的形式打開文件，再利用Image模塊的open方法導入圖片。

對于如下的圖片(圣托里尼)：

圖片的屬性
im.format, im.size, im.mode

得到的結果為：JPEG (1600, 1067) RGB

通道分離：
r,g,b = im.split()

分割成三個通道，此時r,g,b分別為三個圖像對象。

取得像素點的值
im.getpixel((4,4))

由于是RGB三通道的，因此此處的值為：(151, 169, 205)

改變單個像素點的值
im.putpixel(xy, color)

圖像類型轉換：
im=im.convert("L")

由RGB的圖像轉成灰度的圖像，其結果為：

生成新的圖像
Image.new(mode, size)

Image.new(mode, size, color)

如：newImg = Image.new(“GBA”,(640,480),(0,255,0))

保存圖片
im.save("save.gif","GIF")

三、利用K-Means++算法進行圖像分割

3.1、利用K-Means++聚類

在利用K-Means++算法進行圖像分割時，將圖像中的每一個像素點作為一個樣本，對RGB圖像來說，每個樣本包括三維：(151, 169, 205)，通過歸一化，將每個通道的值壓縮到[0,1]區間上。數據的導入和處理如下面程序所示：

#coding:UTF-8
import Image as image
import numpy as np
from KMeanspp import run_kmeanspp

def load_data(file_path):
    '''導入數據
    input:  file_path(string):文件的存儲位置
    output: data(mat):數據
    '''
    f = open(file_path, "rb")  # 以二進制的方式打開圖像文件
    data = []
    im = image.open(f)  # 導入圖片
    m, n = im.size  # 得到圖片的大小
    print m, n
    for i in xrange(m):
        for j in xrange(n):
            tmp = []
            x, y, z = im.getpixel((i, j))
            tmp.append(x / 256.0)
            tmp.append(y / 256.0)
            tmp.append(z / 256.0)
            data.append(tmp)
    f.close()
    return np.mat(data)
最終保存成矩陣的形式，矩陣的行為樣本的個數，列為每一個通道的數值(RGB)。在利用K-Means++算法對樣本進行聚類。主函數如下述代碼所示：

if __name__ == "__main__":
k = 10#聚類中心的個數
# 1、導入數據
print "---------- 1.load data ------------"
data = load_data("001.jpg")
# 2、利用kMeans++聚類
print "---------- 2.run kmeans++ ------------"
run_kmeanspp(data, k)
k表示的是聚類的個數。K-Means++程序的實現如下面程序所示：

# coding:UTF-8
'''
Date:20160923
@author: zhaozhiyong
'''

import numpy as np
from random import random
from KMeans import distance, kmeans, save_result

FLOAT_MAX = 1e100  # 設置一個較大的值作為初始化的最小的距離

def nearest(point, cluster_centers):
    '''計算point和cluster_centers之間的最小距離
    input:  point(mat):當前的樣本點
        cluster_centers(mat):當前已經初始化的聚類中心
    output: min_dist(float):點point和當前的聚類中心之間的最短距離
    '''
    min_dist = FLOAT_MAX
    m = np.shape(cluster_centers)[0]  # 當前已經初始化的聚類中心的個數
    for i in xrange(m):
        # 計算point與每個聚類中心之間的距離
        d = distance(point, cluster_centers[i, ])
        # 選擇最短距離
        if min_dist > d:
            min_dist = d
    return min_dist

def get_centroids(points, k):
    '''KMeans++的初始化聚類中心的方法
    input:  points(mat):樣本
        k(int):聚類中心的個數
    output: cluster_centers(mat):初始化后的聚類中心
    '''
    m, n = np.shape(points)
    cluster_centers = np.mat(np.zeros((k , n)))
    # 1、隨機選擇一個樣本點為第一個聚類中心
    index = np.random.randint(0, m)
    cluster_centers[0, ] = np.copy(points[index, ])
    # 2、初始化一個距離的序列
    d = [0.0 for _ in xrange(m)]

    for i in xrange(1, k):
        sum_all = 0
        for j in xrange(m):
            # 3、對每一個樣本找到最近的聚類中心點
            d[j] = nearest(points[j, ], cluster_centers[0:i, ])
            # 4、將所有的最短距離相加
            sum_all += d[j]
        # 5、取得sum_all之間的隨機值
        sum_all *= random()
        # 6、獲得距離最遠的樣本點作為聚類中心點
        for j, di in enumerate(d):
            sum_all -= di
            if sum_all > 0:
                continue
            cluster_centers[i] = np.copy(points[j, ])
            break
    return cluster_centers

def run_kmeanspp(data, k):
    # 1、KMeans++的聚類中心初始化方法
    print "\t---------- 1.K-Means++ generate centers ------------"
    centroids = get_centroids(data, k)
    # 2、聚類計算
    print "\t---------- 2.kmeans ------------"
    subCenter = kmeans(data, k, centroids)
    # 3、保存所屬的類別文件
    print "\t---------- 3.save subCenter ------------"
    save_result("sub_pp", subCenter)
    # 4、保存聚類中心
    print "\t---------- 4.save centroids ------------"
save_result("center_pp", centroids)
在上述代碼中主要是初始化k個聚類中心，K-Means算法的核心程序如下所示：

# coding:UTF-8
'''
Date:20160923
@author: zhaozhiyong
'''
import numpy as np

def distance(vecA, vecB):
    '''計算vecA與vecB之間的歐式距離的平方
    input:  vecA(mat)A點坐標
        vecB(mat)B點坐標
    output: dist[0, 0](float)A點與B點距離的平方
    '''
    dist = (vecA - vecB) * (vecA - vecB).T
    return dist[0, 0]

def randCent(data, k):
    '''隨機初始化聚類中心
    input:  data(mat):訓練數據
        k(int):類別個數
    output: centroids(mat):聚類中心
    '''
    n = np.shape(data)[1]  # 屬性的個數
    centroids = np.mat(np.zeros((k, n)))  # 初始化k個聚類中心
    for j in xrange(n):  # 初始化聚類中心每一維的坐標
        minJ = np.min(data[:, j])
        rangeJ = np.max(data[:, j]) - minJ
        # 在最大值和最小值之間隨機初始化
        centroids[:, j] = minJ * np.mat(np.ones((k , 1))) + np.random.rand(k, 1) * rangeJ
    return centroids

def kmeans(data, k, centroids):
    '''根據KMeans算法求解聚類中心
    input:  data(mat):訓練數據
        k(int):類別個數
        centroids(mat):隨機初始化的聚類中心
    output: centroids(mat):訓練完成的聚類中心
        subCenter(mat):每一個樣本所屬的類別
    '''
    m, n = np.shape(data)  # m：樣本的個數，n：特征的維度
    subCenter = np.mat(np.zeros((m, 2)))  # 初始化每一個樣本所屬的類別
    change = True  # 判斷是否需要重新計算聚類中心
    while change == True:
        change = False  # 重置
        for i in xrange(m):
            minDist = np.inf  # 設置樣本與聚類中心之間的最小的距離，初始值為爭取窮
            minIndex = 0  # 所屬的類別
            for j in xrange(k):
                # 計算i和每個聚類中心之間的距離
                dist = distance(data[i, ], centroids[j, ])
                if dist < minDist:
                    minDist = dist
                    minIndex = j
            # 判斷是否需要改變
            if subCenter[i, 0] <> minIndex:  # 需要改變
                change = True
                subCenter[i, ] = np.mat([minIndex, minDist])
        # 重新計算聚類中心
        for j in xrange(k):
            sum_all = np.mat(np.zeros((1, n)))
            r = 0  # 每個類別中的樣本的個數
            for i in xrange(m):
                if subCenter[i, 0] == j:  # 計算第j個類別
                sum_all += data[i, ]
                r += 1
            for z in xrange(n):
                try:
                    centroids[j, z] = sum_all[0, z] / r
                    print r
                except:
                    print " r is zero"   
    return subCenter

def save_result(file_name, source):
    '''保存source中的結果到file_name文件中
    input:  file_name(string):文件名
        source(mat):需要保存的數據
    output:
    '''
    m, n = np.shape(source)
    f = open(file_name, "w")
    for i in xrange(m):
        tmp = []
        for j in xrange(n):
            tmp.append(str(source[i, j]))
        f.write("\t".join(tmp) + "\n")
    f.close()
3.2、利用聚類結果生成新的圖片

上述的過程中，對每一個像素點進行了聚類，最終利用聚類中心點的RGB值替換原圖中每一個像素點的值，便得到了最終的分割后的圖片，代碼如下所示：數據分析師培訓

#coding:UTF-8

import Image as image

f_center = open("center_pp")

center = []
for line in f_center.readlines():
    lines = line.strip().split("\t")
    tmp = []
    for x in lines:
        tmp.append(int(float(x) * 256))
    center.append(tuple(tmp))
print center
f_center.close()

fp = open("001.jpg", "rb")
im = image.open(fp)
# 新建一個圖片
m, n = im.size
pic_new = image.new("RGB", (m, n))

f_sub = open("sub_pp")
i = 0
for line in f_sub.readlines():
    index = float((line.strip().split("\t"))[0])
    index_n = int(index)
    pic_new.putpixel(((i/n),(i % n)),center[index_n])
    i = i + 1
f_sub.close()

pic_new.save("result.jpg", "JPEG")       
對于上述的圣托里尼的圖片，取不同的k值，得到如下的一些結果：
原圖

k=3

k=5

k=7

k=10