Press "Enter" to skip to content

## 聚类

### 算法

k均值聚类(K-means) 这篇文章中举了一个很不错的应用例子，作者用亚洲15支足球队的2005年到2010年的战绩做了一个向量表，然后用K-Means算法(K=3)把球队聚为三个类，得出了下面的结果，非常真实。

```亚洲一流：日本，韩国，伊朗，沙特

## 代码实战部分

```# segmentation.py
# python 3.6
import numpy as np
import random
from scipy.spatial.distance import squareform, pdist
from skimage.util import img_as_float
### Clustering Methods
def kmeans(features, k, num_iters=100):
N, D = features.shape
assert N >= k, 'Number of clusters cannot be greater than number of points'
# Randomly initalize cluster centers
idxs = np.random.choice(N, size=k, replace=False)
centers = features[idxs]        # 1. 随机中心点
assignments = np.zeros(N)
for n in range(num_iters):
### YOUR CODE HERE
# 2. 分类
for i in range(N):
dist = np.linalg.norm(features[i] - centers, axis=1)    # 每个点和中心点的距离
assignments[i] = np.argmin(dist)        # 第i个点属于最近的中心点
pre_centers = centers.copy()
# 3. 重新计算中心点
for j in range(k):
centers[j] = np.mean(features[assignments == j], axis=0)
# 4. 验证中心点是否改变
if np.array_equal(pre_centers, centers):
break
### END YOUR CODE
return assignments
def kmeans_fast(features, k, num_iters=100):
N, D = features.shape
assert N >= k, 'Number of clusters cannot be greater than number of points'
# Randomly initalize cluster centers
idxs = np.random.choice(N, size=k, replace=False)
centers = features[idxs]
assignments = np.zeros(N)
for n in range(num_iters):
### YOUR CODE HERE
# 计算距离
features_tmp = np.tile(features, (k, 1))        # (k*N, ...)
centers_tmp = np.repeat(centers, N, axis=0)     # (N * k, ...)
dist = np.sum((features_tmp - centers_tmp)**2, axis=1).reshape((k, N))      # 每列 即k个中心点
assignments = np.argmin(dist, axis=0)   # 最近
# 计算新的中心点
pre_centers = centers
# 3. 重新计算中心点
for j in range(k):
centers[j] = np.mean(features[assignments == j], axis=0)
# 4. 验证中心点是否改变
if np.array_equal(pre_centers, centers):
break
### END YOUR CODE
return assignments

def hierarchical_clustering(features, k):
N, D = features.shape
assert N >= k, 'Number of clusters cannot be greater than number of points'
# Assign each point to its own cluster
assignments = np.arange(N)
centers = np.copy(features)
n_clusters = N
while n_clusters > k:
### YOUR CODE HERE
dist = pdist(centers)       # 计算相互之间的距离
matrixDist = squareform(dist)   # 将向量形式变化为矩阵形式
matrixDist = np.where(matrixDist != 0.0, matrixDist, 1e10)      # 将0.0的变为1e10,即为了矩阵中相同的点计算的距离去掉
minValue = np.argmin(matrixDist)        # 最小的值的位置
min_i = minValue // n_clusters          # 行号
min_j = minValue - min_i * n_clusters   # 列号
if min_j < min_i:       # 归并到小号的cluster
min_i, min_j = min_j, min_i  # 交换一下
for i in range(N):
if assignments[i] == min_j:
assignments[i] = min_i     # 两者合并
for i in range(N):
if assignments[i] > min_j:
assignments[i] -= 1     # 合并了一个cluster,因此n_clusters减少一位
centers = np.delete(centers, min_j, axis=0)  # 减少一个
centers[min_i] = np.mean(features[assignments == min_i], axis=0)        # 重新计算中心点
n_clusters -= 1     # 减去1
### END YOUR CODE
return assignments

### Pixel-Level Features
def color_features(img):
H, W, C = img.shape
img = img_as_float(img)
features = np.zeros((H*W, C))
### YOUR CODE HERE
features = img.reshape(H * W, C)        # color作为特征
### END YOUR CODE
return features
def color_position_features(img):
H, W, C = img.shape
color = img_as_float(img)
features = np.zeros((H*W, C+2))
### YOUR CODE HERE
# 坐标
cord = np.dstack(np.mgrid[0:H, 0:W]).reshape((H*W, 2))      # mgrid生成坐标，重新格式为（x,y）的二维
features[:, 0:C] = color.reshape((H*W, C))      # r,g,b
features[:, C:C+2] = cord
features = (features - np.mean(features, axis=0)) / np.std(features, axis=0,  ddof = 0)     # 对特征归一化处理
### END YOUR CODE
return features
def my_features(img):
""" Implement your own features
Args:
img - array of shape (H, W, C)
Returns:
features - array of (H * W, C)
"""
features = None
### YOUR CODE HERE
features = color_position_features(img)
### END YOUR CODE
return features

### Quantitative Evaluation
def compute_accuracy(mask_gt, mask):
accuracy = None
### YOUR CODE HERE
mask_end = mask_gt - mask
count = len(mask_end[np.where(mask_end == 0)])
accuracy = count / (mask_gt.shape[0] * mask_gt.shape[1])
### END YOUR CODE
return accuracy
def evaluate_segmentation(mask_gt, segments):
num_segments = np.max(segments) + 1
best_accuracy = 0
# 将分割结果与真实值进行对比
for i in range(num_segments):
mask = (segments == i).astype(int)
accuracy = compute_accuracy(mask_gt, mask)
best_accuracy = max(accuracy, best_accuracy)
return best_accuracy```

```# test.py
import  numpy as np
from scipy.spatial.distance import  pdist, squareform
if __name__ == "__main__":
a = np.array([[1,1,1,1],[1,0,0,0]])
b = np.array([[1,0,0,1],[1,1,0,0]])
c = (a == 1)
print(c)```

```# utils.py
import numpy as np
import matplotlib.pyplot as plt
from skimage.util import img_as_float
from skimage import transform
from skimage import io
from segmentation import *
import os
def visualize_mean_color_image(img, segments):
img = img_as_float(img)
k = np.max(segments) + 1
mean_color_img = np.zeros(img.shape)
for i in range(k):
mean_color = np.mean(img[segments == i], axis=0)
mean_color_img[segments == i] = mean_color
plt.imshow(mean_color_img)
plt.axis('off')
plt.show()
def compute_segmentation(img, k,
clustering_fn=kmeans_fast,
feature_fn=color_position_features,
scale=0):
""" 计算图像分割结果
首先，从图像的每个像素中提取一个特征向量。然后将聚类算法应用于所有特征向量的集合。当且仅当两个像素的特征向量被分配到同一簇时，两个像素会被分配到同一聚簇。
"""
assert scale <= 1 and scale >= 0, \
'Scale should be in the range between 0 and 1'
H, W, C = img.shape
if scale > 0:
# 缩小图像来获得更快的计算速度
img = transform.rescale(img, scale)
features = feature_fn(img)
assignments = clustering_fn(features, k)
segments = assignments.reshape((img.shape[:2]))
if scale > 0:
# 调整大小分割回图像的原始大小
segments = transform.resize(segments, (H, W), preserve_range=True)
# 调整大小会导致像素值不重叠。
# 像素值四舍五入为最接近的整数
segments = np.rint(segments).astype(int)
return segments

def load_dataset(data_dir):
"""
载入数据集
'imgs/aaa.jpg' is 'gt/aaa.png'
"""
imgs = []
gt_masks = []
for fname in sorted(os.listdir(os.path.join(data_dir, 'imgs'))):
if fname.endswith('.jpg'):
# 读入图像
img = io.imread(os.path.join(data_dir, 'imgs', fname))
imgs.append(img)
# 加载相应的分割mask
mask_fname = fname[:-4] + '.png'
gt_mask = io.imread(os.path.join(data_dir, 'gt', mask_fname))
gt_mask = (gt_mask != 0).astype(int) # 将mask进行二值化
gt_masks.append(gt_mask)
return imgs, gt_masks```

## 运行实施例

```# 初始化
from time import time
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
from skimage import io
from __future__ import print_function
%matplotlib内联
plt.rcParams['figure.figsize'] = (15.0, 12.0) # set default size of plots
plt.rcParams['image.interpolation'] = 'nearest'
plt.rcParams['image.cmap'] = 'gray'
# 自动重载外部模块
%load_ext autoreload
%autoreload 2```

```# 为聚类生成随机数据点
# 采用seed保证生成结果的一致
np.random.seed(0)
# 聚簇 1
mean1 = [-1, 0]
cov1 = [[0.1, 0], [0, 0.1]]
X1 = np.random.multivariate_normal(mean1, cov1, 100)
# 聚簇 2
mean2 = [0, 1]
cov2 = [[0.1, 0], [0, 0.1]]
X2 = np.random.multivariate_normal(mean2, cov2, 100)
# 聚簇 3
mean3 = [1, 0]
cov3 = [[0.1, 0], [0, 0.1]]
X3 = np.random.multivariate_normal(mean3, cov3, 100)
# 聚簇 4
mean4 = [0, -1]
cov4 = [[0.1, 0], [0, 0.1]]
X4 = np.random.multivariate_normal(mean4, cov4, 100)
# 合并两组数据点
X = np.concatenate((X1, X2, X3, X4))
# 绘制数据点
plt.scatter(X[:, 0], X[:, 1])
plt.axis('equal')
plt.show()```

```from segmentation import kmeans
np.random.seed(0)
start = time()
assignments = kmeans(X, 4)
end = time()
kmeans_runtime = end - start
print("kmeans running time: %f seconds." % kmeans_runtime)
for i in range(4):
cluster_i = X[assignments==i]
plt.scatter(cluster_i[:, 0], cluster_i[:, 1])
plt.axis('equal')
plt.show()```

kmeans running time: 0.027956 seconds.

```from segmentation import hierarchical_clustering
start = time()
assignments = hierarchical_clustering(X, 4)
end = time()
print("hierarchical_clustering running time: %f seconds." % (end - start))
for i in range(4):
cluster_i = X[assignments==i]
plt.scatter(cluster_i[:, 0], cluster_i[:, 1])
plt.axis('equal')
plt.show()```

hierarchical_clustering running time: 0.793070 seconds.

```# 载入并显示图像
img = io.imread('train.jpg')
H, W, C = img.shape
plt.imshow(img)
plt.axis('off')
plt.show()```

```from segmentation import color_features
np.random.seed(0)
features = color_features(img)
# 结果检测
assert features.shape == (H * W, C),\
"Incorrect shape! Check your implementation."
assert features.dtype == np.float,\
"dtype of color_features should be float."
assignments = kmeans_fast(features, 8)
segments = assignments.reshape((H, W))
# 展示图像分割结果
plt.imshow(segments, cmap='viridis')
plt.axis('off')
plt.show()```

```from utils import visualize_mean_color_image
visualize_mean_color_image(img, segments)```

```from segmentation import color_position_features
np.random.seed(0)
features = color_position_features(img)
# 结果检测
assert features.shape == (H * W, C + 2),\
"Incorrect shape! Check your implementation."
assert features.dtype == np.float,\
"dtype of color_features should be float."
assignments = kmeans_fast(features, 8)
segments = assignments.reshape((H, W))
# 图像分割结果显示
plt.imshow(segments, cmap='viridis')
plt.axis('off')
plt.show()```

`visualize_mean_color_image(img, segments)`

```from utils import load_dataset, compute_segmentation
from segmentation import evaluate_segmentation
# 载入该小型数据集
imgs, gt_masks = load_dataset('./data')
# 设置图像分割的参数
num_segments = 3
clustering_fn = kmeans_fast
feature_fn = color_features
scale = 0.5
mean_accuracy = 0.0
segmentations = []
for i, (img, gt_mask) in enumerate(zip(imgs, gt_masks)):
# Compute a segmentation for this image
segments = compute_segmentation(img, num_segments,
clustering_fn=clustering_fn,
feature_fn=feature_fn,
scale=scale)

segmentations.append(segments)

# 评估图像分割结果
accuracy = evaluate_segmentation(gt_mask, segments)

print('Accuracy for image %d: %0.4f' %(i, accuracy))
mean_accuracy += accuracy

mean_accuracy = mean_accuracy / len(imgs)
print('Mean accuracy: %0.4f' % mean_accuracy)```

```Accuracy for image 0: 0.8612
Accuracy for image 1: 0.9571
Accuracy for image 2: 0.9824
Accuracy for image 3: 0.9206
Accuracy for image 4: 0.7642
Accuracy for image 5: 0.8062
Accuracy for image 6: 0.6617
Accuracy for image 7: 0.4726
Accuracy for image 8: 0.8317
Accuracy for image 9: 0.7580
Accuracy for image 10: 0.6515
Accuracy for image 11: 0.8261
Accuracy for image 12: 0.7105
Accuracy for image 13: 0.6667
Accuracy for image 14: 0.7623
Accuracy for image 15: 0.5223
Mean accuracy: 0.7597```

```# 可视化图像分割结果
N = len(imgs)
plt.figure(figsize=(15,60))
for i in range(N):
plt.subplot(N, 3, (i * 3) + 1)
plt.imshow(imgs[i])
plt.axis('off')
plt.subplot(N, 3, (i * 3) + 2)
plt.imshow(gt_masks[i])
plt.axis('off')
plt.subplot(N, 3, (i * 3) + 3)
plt.imshow(segmentations[i], cmap='viridis')
plt.axis('off')
plt.show()```