K-means 聚类实验：三种实现的比较

第一次课程实验内容:

目标: 主要目标是对图像数据进行 K-means 聚类分析，并比较不同版本的 K-means 算法的性能

import numpy as np
import time
import matplotlib.pyplot as plt  # 用于添加标题
import matplotlib.image as mpimg  # 用于读取图像文件
import seg_tools  # 假设是自定义的库，提供 init() 和 draw_labels() 等功能
from sklearn.cluster import KMeans

# 定义 kmeans_v1
def kmeans_v1(data, centers):
    dim_v = np.size(v)
    # vi = np.reshape(v, (dim_v, 1))
    dim_c = np.size(c)
    label = np.zeros((dim_v, 1), dtype=int)
    i = 0
    flag = 1
    while flag:
        i += 1
        for j in range(0, dim_v):
            vij = v[j]
            for k in range(0, dim_c):
                d = abs(vij - c[k])
                if k == 0 or d < d_min:
                    d_min = d
                    k_min = k
            label[j] = k_min
        for k in range(0, 3):
            c[k] = np.mean(v[label == k])
        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:
            flag = 0
        label_mem = np.copy(label)
    return label, i


# 定义 kmeans_v2
def kmeans_v2(v, c):
    dim_v = np.size(v)
    # vi = np.reshape(v, (dim_v, 1))
    dim_c = np.size(c)
    ci = np.reshape(c, (1, dim_c))
    i = 0
    flag = 1
    while flag:
        i += 1
        d = np.abs(v * np.ones((1, dim_c)) - np.ones((dim_v, 1)) * ci)
        label = np.argmin(d, 1)
        for k in range(0, dim_c):
            ci[0, k] = np.mean(v[label == k])
        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:
            flag = 0
        label_mem = np.copy(label)
    return label, i


# 定义 kmeans_v3
def kmeans_v3(v, K):
    data_set = np.dstack(v.flatten())
    model = KMeans(n_clusters=K, n_init=10)
    model.fit(data_set[0].T)
    return model.labels_


# 主程序部分
seg_tools.init()
I = mpimg.imread('implant.bmp')
dim_y, dim_x = np.shape(I)
v = np.reshape(I.astype(float), (dim_y * dim_x, 1))
c = np.array([50, 150, 200])

# 使用 kmeans_v1
t1 = time.perf_counter()
for n in range(0, 100):
    label1, i = kmeans_v1(v, c)
t2 = time.perf_counter()
print('V1 iterations:', i)
print('V1 temps:', (t2 - t1) / 100)

# 使用 kmeans_v2
t1 = time.perf_counter()
for n in range(0, 1000):
    label2, i = kmeans_v2(v, c)
t2 = time.perf_counter()
print('V2 iterations:', i)
print('V2 temps:', (t2 - t1) / 1000)

# 使用 kmeans_v3
t1 = time.perf_counter()
for n in range(0, 10):
    label3 = kmeans_v3(v, 3)
t2 = time.perf_counter()
print('V3 temps:', (t2 - t1) / 10)

# 绘制结果
seg_tools.draw_uint8(I)
seg_tools.draw_labels(np.reshape(label1, (dim_y, dim_x)))
seg_tools.draw_labels(np.reshape(label2, (dim_y, dim_x)))
seg_tools.draw_labels(np.reshape(label3, (dim_y, dim_x)))

以下是代码详细解释

函数定义

kmeans_v1

接下来我们要实现第一个版本的 K-means 算法，是最基础的实现:

使用双重循环逐点计算每个数据点到各个聚类中心的距离，然后选择最近的聚类中心。
对于每个聚类中心，计算对应数据点的均值作为新的中心。

# 定义 kmeans_v1
def kmeans_v1(data, centers):
    dim_v = np.size(v)
    # vi = np.reshape(v, (dim_v, 1))
    dim_c = np.size(c)
    label = np.zeros((dim_v, 1), dtype=int)
    i = 0
    flag = 1
    while flag:
        i += 1
        for j in range(0, dim_v):
            vij = v[j]
            for k in range(0, dim_c):
                d = abs(vij - c[k])
                if k == 0 or d < d_min:
                    d_min = d
                    k_min = k
            label[j] = k_min
        for k in range(0, 3):
            c[k] = np.mean(v[label == k])
        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:
            flag = 0
        label_mem = np.copy(label)
    return label, i

它有两个核心步骤，a) 将数据点归类到最近的中心 b) 更新中心位置

def kmeans_v1(data, centers):

参数 data 是需要聚类的数据（图像转换成的像素值向量）。
参数 centers 是初始的聚类中心。

    dim_v = np.size(v)

计算数据向量的大小，用 np.size(v) 得到数据点的总数。比如，如果图像有 100 \times 100 像素，展开后共有 10,000 个数据点。

    dim_c = np.size(c)

计算初始聚类中心的数量，用 np.size(c) 得到中心的个数。在这里，初始中心是 [50, 150, 200] ，因此聚类中心数是 3 。

    label = np.zeros((dim_v, 1), dtype=int)

初始化标签数组 label，每个数据点的标签会存储在这个数组中。它的作用是表示每个数据点属于哪个聚类中心。

label 的维度是 (\text{dim\_v}, 1) ，初始值全部为 0 。

    i = 0

初始化迭代计数器 i，用来记录算法进行了几轮迭代。

    flag = 1

定义一个标志位 flag，控制 while 循环。只要 flag 为 1 ，算法会继续迭代，直到满足停止条件。

    while flag:

while 循环是 K-means 算法的核心部分，每次循环代表一轮迭代。

        i += 1

每次进入循环，迭代计数器 i 加 1 ，表示当前是第 i 次迭代。

        for j in range(0, dim_v):

用一个 for 循环遍历每个数据点。 j 从 0 到 \text{dim\_v}-1 。

            vij = v[j]

获取第 j 个数据点的值，并存储到变量 vij 中。例如，如果第 j 个像素值是 120 ，则 vij = 120。

            for k in range(0, dim_c):

再用一个 for 循环，遍历每个聚类中心。 k 是聚类中心的索引，从 0 到 \text{dim\_c}-1 。

                d = abs(vij - c[k])

计算当前数据点到第 k 个聚类中心的距离，在这里，使用的是绝对值距离 \lvert vij - c[k] \rvert 。

例如，如果 vij=120，而当前中心 c[k]=150 ，则 d = 30 。

                if k == 0 or d < d_min:

判断条件：

如果这是第一个聚类中心（ k=0 ），直接认为当前距离是最小距离。
如果当前距离 d 小于之前记录的最小距离 d_min，则更新最小距离。

                    d_min = d

将当前距离 d 赋值为新的最小距离 d_min。

                    k_min = k

将当前聚类中心的索引 k 赋值给 k_min，表示这是当前最近的聚类中心。

            label[j] = k_min

在所有聚类中心中找到最近的一个，将其索引 k_min 分配给数据点 j 的标签。

例如，如果数据点 v[10] 距离中心 c[2] 最近，则 label[10] = 2。

        for k in range(0, 3):

遍历所有聚类中心，准备更新每个中心的位置。假设有 3 个聚类中心，因此循环范围是 0, 1, 2 。

            c[k] = np.mean(v[label == k])

更新第 k 个聚类中心的位置。

取所有属于该聚类中心（即标签为 k 的数据点）的均值作为新中心位置。
例如，如果 k=0 ，而标签为 0 的数据点是 [50, 52, 48] ，则 c[0] = (50 + 52 + 48)/3 = 50 。

        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:

停止条件，若满足以下这两个条件，说明聚类已经收敛：

如果当前迭代不是第一次（ i > 1 ）。
所有数据点的标签在本轮迭代中没有变化，即标签数组 label 与上一次的 label_mem 相等。

            flag = 0

如果满足停止条件，将标志位 flag 置为 0 ，停止迭代。

        label_mem = np.copy(label)

保存当前迭代的标签数组 label，用于下一轮比较。

    return label, i

返回结果：

label 是最终的聚类标签数组，每个数据点对应一个聚类中心的索引。
i 是总共的迭代次数，用来评估算法的收敛速度。

kmeans_v2

接下来进行第二版的 K-means 算法，这部分算法主要改进了距离计算和标签更新

距离计算: 采用矩阵运算，用 NumPy 向量化实现批量计算所有数据点到各聚类中心的距离（避免了双重功能循环）
**标签更新：**直接通过 np.argmin(d, axis=1) 找到距离最近的聚类中心，效率比手动循环高得多。

# 定义 kmeans_v2
def kmeans_v2(v, c):
    dim_v = np.size(v)
    # vi = np.reshape(v, (dim_v, 1))
    dim_c = np.size(c)
    ci = np.reshape(c, (1, dim_c))
    i = 0
    flag = 1
    while flag:
        i += 1
        d = np.abs(v * np.ones((1, dim_c)) - np.ones((dim_v, 1)) * ci)
        label = np.argmin(d, 1)
        for k in range(0, dim_c):
            ci[0, k] = np.mean(v[label == k])
        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:
            flag = 0
        label_mem = np.copy(label)
    return label, i

代码详细解释

def kmeans_v2(v, c):

定义函数 kmeans_v2，接收两个参数：v 是展开后的数据向量，c 是初始的聚类中心数组。

    dim_v = np.size(v)

计算数据点的总数 dim_v，表示数据集中有多少个点。

    dim_c = np.size(c)

计算聚类中心的总数 dim_c，即中心的个数。

    ci = np.reshape(c, (1, dim_c))

将聚类中心数组 c 重塑为一个 1 行 dim_c 列的二维数组 ci，方便后续矩阵运算。

    i = 0

初始化迭代计数器 i，记录算法运行的轮数。

    flag = 1

初始化标志位 flag，用于控制 while 循环，表示是否继续迭代。

    while flag:

开始 while 循环，表示每一轮的聚类操作。

        i += 1

每轮迭代开始时将计数器 i 加 1，记录当前是第几轮。

        d = np.abs(v * np.ones((1, dim_c)) - np.ones((dim_v, 1)) * ci)

计算每个数据点到所有聚类中心的距离矩阵 d。通过 NumPy 的矩阵运算：

v * np.ones((1, dim_c)) 将数据点向量扩展为一个与中心数组形状一致的矩阵；
np.ones((dim_v, 1)) * ci 将聚类中心扩展为与数据点匹配的矩阵；
np.abs(...) 逐元素计算绝对值差，得到距离矩阵 d。

        label = np.argmin(d, 1)

对距离矩阵 d 的每一行找出最小值的索引，将每个数据点分配到最近的聚类中心。结果存储在标签数组 label 中。

        for k in range(0, dim_c):

启动一个循环，逐一更新每个聚类中心的位置。

            ci[0, k] = np.mean(v[label == k])

更新第 k 个聚类中心的位置。筛选出标签为 k 的数据点，计算它们的均值，作为新中心位置。

        if i > 1 and np.sum(np.abs(label - label_mem)) == 0:

检查停止条件：如果当前不是第一次迭代，且标签数组 label 与上一轮相比没有变化，说明聚类已经收敛。

            flag = 0

如果满足停止条件，将标志位 flag 置为 0，停止循环。

        label_mem = np.copy(label)

保存当前的标签数组 label，用于下一轮迭代的比较。

    return label, i

返回最终结果：label 是聚类标签数组，i 是总共的迭代次数。

kmeans_v3

接下来进行第三版的 K-means 算法，这部分算法同样改进了距离计算，速度最快

距离计算：使用了优化的底层实现，速度更快，支持多种距离度量方式。

# 定义 kmeans_v3
def kmeans_v3(v, K):
    data_set = np.dstack(v.flatten())
    model = KMeans(n_clusters=K, n_init=10)
    model.fit(data_set[0].T)
    return model.labels_

代码详细解释:

def kmeans_v3(v, K):

定义函数 kmeans_v3，接收两个参数：

v：数据向量，表示展开后的图像像素值数组。
K：聚类中心的数量。

    data_set = np.dstack(v.flatten())

将数据向量 v 转换为三维数组 data_set：

v.flatten() 将多维数组 v 展平为一维数组。
np.dstack() 将一维数组扩展为三维数组，便于后续操作。
结果 data_set 的形状是 (1, 1, \text{dim\_v}) ，其中 \text{dim\_v} 是数据点的数量。

    model = KMeans(n_clusters=K, n_init=10)

创建一个 KMeans 模型实例：

n_clusters=K 指定聚类中心的数量。
n_init=10 设置模型重复运行 10 次，每次随机初始化聚类中心，并选择最优结果。

    model.fit(data_set[0].T)

使用 fit() 方法对数据进行聚类：

data_set[0] 提取三维数组的第一层，变成一个二维数组，形状是 (1, \text{dim\_v}) 。
.T 对数组转置，得到 (\text{dim\_v}, 1) 的形状，与 KMeans 的输入要求匹配。
此步操作会对输入数据执行聚类并计算最终的聚类中心。

    return model.labels_

返回 model.labels_，即每个数据点的聚类标签。这是一个一维数组，长度为 \text{dim\_v} ，每个值对应数据点所属的聚类中心编号（从 0 到 K-1 ）。

主程序部分

I = mpimg.imread('implant.bmp')

使用 matplotlib.image.imread 读取 implant.bmp 图像。

dim_y, dim_x = np.shape(I)

将二维图像数据转换为一维向量 v，以便后续进行聚类。

v = np.reshape(I.astype(float), (dim_y * dim_x, 1))

转换为浮点数以满足 K-means 的输入要求。

c = np.array([50, 150, 200])

设置初始的聚类中心，即 3 个灰度值：50（黑色）、150 （中灰色）和 200（白色）

# 使用 kmeans_v1
t1 = time.perf_counter()
for n in range(0, 100):
    label1, i = kmeans_v1(v, c)
t2 = time.perf_counter()
print('V1 iterations:', i)
print('V1 temps:', (t2 - t1) / 100)

逐行解释：

t1 = time.perf_counter()

使用 time.perf_counter() 记录当前时间，作为计时的起点。

for n in range(0, 100):

循环 1000 次，用于多次运行 kmeans_v1 算法以统计其平均运行时间，减小单次测量中的随机波动，得到更准确的性能评估。

    label1, i = kmeans_v1(v, c)

在每次循环中调用 kmeans_v2 函数，传入数据向量 v 和初始聚类中心 c。label2 是返回的聚类标签数组，表示每个数据点被分配到哪个聚类中心。

t2 = time.perf_counter()

在循环结束后，再次使用 time.perf_counter() 记录当前时间，作为计时的终点。

print('V1 iterations:', i)

输出最后一次运行 kmeans_v1 的迭代次数 i。

print('V1 temps:', (t2 - t1) / 100)

计算并输出 kmeans_v1 的平均运行时间：

(t2 - t1) 是总运行时间，表示 100 次运行的累计耗时（单位是秒），总时间除以 1000 ，得到单次运行的平均时间。

其他 kmeans_v2 kmeans_v3 同理

菜单

分享

K-means 聚类实验：三种实现的比较

函数定义

主程序部分

评论

优化理论（一）：导论、希尔伯特空间与凸性

Hello Halo

延迟系统稳定化的谱方法（研讨会）

置信函数理论在图像处理与多源融合中的应用

信念函数理论导论（研讨会报告）

心血管模型比较与 HRV 相干性分析（Master ATSI 团队项目）

医学信号处理 TP：相干函数性质证明与脑连接性实验

医学信号处理 TP：第 2 部分问题解答（LMS 自适应噪声消除）

医学信号处理 TP：自适应滤波去眼部伪差、ICA 与脑连接性（完整解答）

医学信号处理 TP（ENSEIRB-MATMECA 2023–2024）：实践工作题目