三维点云-聚类

1 概述

聚类是一种机器学习和数据挖掘技术,用于将数据集中的对象划分为具有相似特征的多个组。这些组被称为“簇”,而包含相似对象的簇内对象之间的相似度要高于其他簇内对象之间的相似度。

聚类算法的目标是发现数据中的内在结构,以便能够更好地理解数据并作出预测。常见的聚类算法包括K均值聚类、层次聚类、DBSCAN(基于密度的空间聚类)等。

聚类在各种领域中都有广泛的应用,包括市场营销、社交网络分析、生物信息学、图像分析等。通过聚类,我们可以发现数据中的模式和规律,从而做出更准确的预测和决策。

2 算法

常见的聚类算法包括:

1. K均值聚类(K-means clustering):这是一种基于距离的聚类方法,它将数据点分成K个簇,每个簇都有一个代表性的中心点,然后通过迭代优化来最小化数据点与其所属簇中心点之间的距离。

2 层次聚类(Hierarchical clustering):这种方法根据数据点之间的相似性逐渐构建聚类层次结构,可以是自下而上的凝聚聚类(agglomerative clustering)或自上而下的分裂聚类(divisive clustering)。

3. DSCAN(Density-Based Spatial Clustering of Applications with Noise):这是一种基于密度的聚类算法,能够识别任意形状的簇,并且能够在噪声和离群点存在的情况下有效地工作。

4. 密度峰值聚类(Density Peak Clustering):该算法通过寻找数据点中的密度峰值来进行聚类,能够有效处理具有不同密度的簇。

5. 高斯混合模型(Gaussian Mixture Model, GMM):这是一种基于概率统计的聚类方法,假设数据是由多个高斯分布组合而成,通过最大期望算法(EM算法)来估计参数。

2.1 K_Means

2.1.1 介绍

K均值(K-means)聚类是一种常见的聚类算法,它将数据点划分为K个簇,每个簇具有自己的中心,目标是最小化数据点与其所属簇中心之间的距离平方和。K均值聚类的基本思想是通过迭代优化来找到最优的簇中心位置,使得簇内的数据点相似度高,而不同簇之间的数据点相似度较低。

K均值聚类的算法步骤如下:

1. 随机初始化K个簇中心点。

2. 对于每个数据点,计算其与各个簇中心的距离,并将其划分到距离最近的簇中。

3. 更新每个簇的中心,使其成为该簇所有数据点的平均值。

4. 重复步骤2和步骤3,直到簇中心不再发生变化,或者达到预先设定的迭代次数。

K均值聚类的优点包括简单易实现、计算效率高;而缺点包括对初始簇中心的选择敏感,对异常值敏感,且需要事先知道簇的数量K。

在实际应用中,K均值聚类常被用于图像压缩、客户分群、文本挖掘等领域。

2.1.2 实现

https://blog.csdn.net/qinlele1994/article/details/106180756/

代码实现如下:

# 文件功能: 实现 K-Means 算法

'''
https://blog.csdn.net/qinlele1994/article/details/106180756/
'''

import numpy as np
import random

class K_Means(object):
    # k是分组数;tolerance‘中心点误差’;max_iter是迭代次数
    def __init__(self, n_clusters=2, tolerance=0.0001, max_iter=300):
        self.k_ = n_clusters
        self.tolerance_ = tolerance
        self.max_iter_ = max_iter

    # 计算各个类别中心点的坐标
    def fit(self, data):
        # 作业1
        # 屏蔽开始
        # 随机选取数据中的中心点
        centers = data[random.sample(range(data.shape[0]), self.k_)]  # 从 range(data.shape[0]) 数据中,随机抽取self.k_ 作为一个列表
        old_centers = np.copy(centers)  # 将旧的中心点 保存到old_centers
        labels = [[] for i in range(self.k_)]
        # [[1, 2], [1.5, 1.8], [5, 8], [8, 8], [1, 0.6], [9, 11]
        for iter_ in range(self.max_iter_):  # 循环一定的次数
            for idx, point in enumerate(data):  # enumerate 函数用于一个可遍历的数据对象组合为一个索引序列,同时列出数据和数据的下标
                # 默认的二范数 就是距离  将每一个点计算到每个中心的距离  有2类,0 ,1 就是计算一点到2个中心点的距离
                diff = np.linalg.norm(old_centers - point, axis=1)  # 一个点分别到两个中心点的距离不同,
                diff2 = (np.argmin(diff))  # np.argmin(diff) 表示最小值在数组中的位置  选取距离小的那一点的索引 也就代表了属于哪个类
                labels[diff2].append(idx)  # 选取距离小的那一点的索引 也就代表了属于哪个类

            for i in range(self.k_):
                points = data[labels[i], :]  # 所有在第k类中的所有点
                centers[i] = points.mean(axis=0)  # 均值 作为新的聚类中心
            if np.sum(
                    np.abs(centers - old_centers)) < self.tolerance_ * self.k_:  # 如果前后聚类中心的距离相差小于self.tolerance_ * self.k_ 输出
                break
            old_centers = np.copy(centers)
        self.centers = centers
        self.fitted = True
        # 屏蔽结束

    #计算出各个点的类别
    def predict(self, p_datas):
        result = []
        # 作业2
        # 屏蔽开始
        if not self.fitted:
            print('Unfitter. ')
            return result
        for point in p_datas:
            diff = np.linalg.norm(self.centers - point, axis=1)
            result.append(np.argmin(diff))
        # 屏蔽结束
        return result

if __name__ == '__main__':
    x = np.array([[1, 2], [1.5, 1.8], [5, 8], [8, 8], [1, 0.6], [9, 11]])
    k_means = K_Means(n_clusters=2)
    k_means.fit(x) #计算聚类中心

    cat = k_means.predict(x) #确定了聚类中心以后,计算每个点属于哪个聚类中心
    print(cat)

运行结果:

PS D:\codePy\pointCloud-ljx\3-PointCloud-Homework> python .\KMeans.py
[0, 0, 1, 1, 0, 1]

2.2 Mean Shift

Mean Shift(均值偏移,也叫均值漂移或均值平移)是一种非参数化的密度估计和聚类算法,用于寻找数据集中的局部极大值点(也称为模式)或聚类中心。它的基本思想是通过不断移动数据点的均值来寻找数据分布的峰值

Mean Shift算法的过程如下:

  1. 初始化每个数据点的位置作为当前的均值。
  2. 对于每个数据点,计算它与其他数据点之间的距离,并将距离小于指定半径的点加入到当前点的邻域内。
  3. 计算邻域内点的平均值,得到新的均值。
  4. 将当前点的位置更新为新的均值。
  5. 重复步骤2-4,直到均值不再发生显著变化或达到预定义的收敛条件。

Mean Shift算法的关键是确定邻域的半径大小,它决定了算法的收敛速度和结果。通常可以通过交叉验证或其他启发式方法来选择合适的半径。

Mean Shift算法在计算上比较简单,且对数据分布的形状没有假设限制。它被广泛用于图像分割、目标跟踪、密度估计等领域。然而,由于Mean Shift算法是基于密度的方法,对于高维数据集和大规模数据集,计算复杂度可能较高。

2.3 GMM

2.3.1 介绍

高斯混合模型(Gaussian Mixture Model, GMM)是一种用来描述多维数据分布的统计模型。它假设数据是由有限个高斯分布组成的混合体,每个高斯分布代表了混合模型中的一个簇。换句话说,GMM可以用来对复杂的数据集进行建模,并且可以发现数据中隐藏的潜在结构。

在GMM中,每个高斯分布都由其均值、协方差矩阵和权重所确定。因此,当使用GMM对数据进行建模时,需要估计这些参数。常用的估计方法是通过最大期望算法(Expectation-Maximization, EM算法)来实现。EM算法通过迭代优化来找到最大似然估计下的参数值,从而能够对数据进行合适的聚类。

GMM在实际应用中被广泛运用于聚类、密度估计、异常检测等领域。它可以灵活地适应不同形状和密度的数据分布,因此在处理复杂数据集时具有一定的优势。同时,GMM也可以用于生成新的样本数据,因为它提供了对数据分布的良好拟合。

2.3.2 实现

https://blog.csdn.net/weixin_43384504/article/details/105754032

代码实现如下:

# 文件功能:实现 GMM 算法

'''
https://blog.csdn.net/weixin_43384504/article/details/105754032
'''

import numpy as np
from numpy import *
import pylab
import random,math

import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.stats import multivariate_normal
plt.style.use('seaborn')

class GMM(object):
    def __init__(self, n_clusters, max_iter=50):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
    
    # 屏蔽开始
    # 更新W
    def update_W(self, X, Mu, Var, Pi):
        n_points, n_clusters = len(X), len(Pi)
        pdfs = np.zeros(((n_points, n_clusters)))
        for i in range(n_clusters):
            pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i], np.diag(Var[i]))
        W = pdfs / pdfs.sum(axis=1).reshape(-1, 1)
        return W

    # 更新pi
    def update_Pi(self, W):
        Pi = W.sum(axis=0) / W.sum()
        return Pi

    # 画出聚类图像
    def plot_clusters(self, X, Mu, Var, Mu_true=None, Var_true=None):
        colors = ['b', 'g', 'r']
        n_clusters = len(Mu)
        plt.figure(figsize=(10, 8))
        plt.axis([-10, 15, -5, 15])
        plt.scatter(X[:, 0], X[:, 1], s=5)
        ax = plt.gca()
        for i in range(n_clusters):
            plot_args = {'fc': 'None', 'lw': 2, 'edgecolor': colors[i], 'ls': ':'}
            ellipse = Ellipse(Mu[i], 3 * Var[i][0], 3 * Var[i][1], **plot_args)
            ax.add_patch(ellipse)
        if (Mu_true is not None) & (Var_true is not None):
            for i in range(n_clusters):
                plot_args = {'fc': 'None', 'lw': 2, 'edgecolor': colors[i], 'alpha': 0.5}
                ellipse = Ellipse(Mu_true[i], 3 * Var_true[i][0], 3 * Var_true[i][1], **plot_args)
                ax.add_patch(ellipse)
        plt.show()

    # 更新Mu
    def update_Mu(self, X, W):
        n_clusters = W.shape[1]
        Mu = np.zeros((n_clusters, 2))
        for i in range(n_clusters):
            Mu[i] = np.average(X, axis=0, weights=W[:, i])
        return Mu

    # 更新Var
    def update_Var(self, X, Mu, W):
        n_clusters = W.shape[1]
        Var = np.zeros((n_clusters, 2))
        for i in range(n_clusters):
            Var[i] = np.average((X - Mu[i]) ** 2, axis=0, weights=W[:, i])
        return Var

    # 计算log似然函数(代价函数)
    def logLH(self, X, Pi, Mu, Var):
        n_points, n_clusters = len(X), len(Pi)
        pdfs = np.zeros(((n_points, n_clusters)))
        for i in range(n_clusters):
            pdfs[:, i] = Pi[i] * multivariate_normal.pdf(X, Mu[i], np.diag(Var[i]))
        return np.mean(np.log(pdfs.sum(axis=1)))
    # 屏蔽结束
    # 这里主要是计算E M steps 返回后验概率值
    def fit(self, data):
        # 作业3
        # 屏蔽开始
        n_points = len(data)
        Mu = [[0, 0], [0, 0], [0, 0]]
        Var = [[1, 1], [1, 1], [1, 1]]
        Pi = [1 / self.n_clusters] * 3
        W = np.ones((n_points, self.n_clusters)) / self.n_clusters
        Pi = W.sum(axis=0) / W.sum()
        loglh = []
        # 给定一批数据点,首先要初始化参数,在主函数中初始化了
        # 这里要重新写一下
        for i in range(5):
            # self.plot_clusters(data, Mu, Var, true_Mu, true_Var)
            loglh.append(self.logLH(data, Pi, Mu, Var))
            # 输出后验概率,每一行中都有属于高斯模型的可能性
            W = self.update_W(data, Mu, Var, Pi)
            # 输出参数
            Pi = self.update_Pi(W)
            Mu = self.update_Mu(data, W)
            print('log-likehood:%.3f' % loglh[-1])
            Var = self.update_Var(data, Mu, W)
        return W
        # 屏蔽结束

    # 每个点在每个分类的概率,要利用上面获得的参数
    def predict(self, data):
        # 屏蔽开始
        W = self.fit(data)
        Label = np.zeros((data.shape[0], 1))
        # 遍历W的每一行 去找概率最大的所以的索引
        Label = W.argmax(axis=1)
        # 屏蔽结束

# 生成仿真数据
def generate_X(true_Mu, true_Var):
    # 第一簇的数据
    num1, mu1, var1 = 400, true_Mu[0], true_Var[0]
    X1 = np.random.multivariate_normal(mu1, np.diag(var1), num1)
    # 第二簇的数据
    num2, mu2, var2 = 600, true_Mu[1], true_Var[1]
    X2 = np.random.multivariate_normal(mu2, np.diag(var2), num2)
    # 第三簇的数据
    num3, mu3, var3 = 1000, true_Mu[2], true_Var[2]
    X3 = np.random.multivariate_normal(mu3, np.diag(var3), num3)
    # 合并在一起
    X = np.vstack((X1, X2, X3))
    # 显示数据
    plt.figure(figsize=(10, 8))
    plt.axis([-10, 15, -5, 15])
    plt.scatter(X1[:, 0], X1[:, 1], s=5)
    plt.scatter(X2[:, 0], X2[:, 1], s=5)
    plt.scatter(X3[:, 0], X3[:, 1], s=5)
    plt.show()
    return X

if __name__ == '__main__':
    # 生成数据
    true_Mu = [[0.5, 0.5], [5.5, 2.5], [1, 7]]
    true_Var = [[1, 3], [2, 2], [6, 2]]
    X = generate_X(true_Mu, true_Var)

    gmm = GMM(n_clusters=3)
    gmm.fit(X)
    cat = gmm.predict(X)
    print(cat)
    # 初始化

运行结果:

图1 GMM聚类

2.4 DBSCAN

DBSCAN(Density-Based Spatial Clustering of Applications with Noise)是一种基于密度的聚类算法,能够发现任意形状的聚类,并且能够有效地识别噪声数据点

DBSCAN根据数据点的密度来进行聚类,相比于K-means等传统聚类算法,它不需要预先指定聚类的数量

DBSCAN算法的核心思想是通过两个参数来定义数据点的密度:邻域半径(ε)和最小邻居数(minPts)。

对于一个数据点,如果它的ε-邻域内包含至少minPts个数据点,则将其视为核心点(core point);如果一个点在另一个核心点的ε-邻域内,那么它也被视为核心点的一部分;剩下的既不是核心点也不在核心点的ε-邻域内的点则被视为噪声点(noise point)。基于这些定义,DBSCAN将数据点划分为核心点、边界点(border point)和噪声点。

DBSCAN算法的运行步骤如下:

  1. 选择一个未被访问的数据点;
  2. 检查该点的ε-邻域内是否包含至少minPts个数据点,若是,则将其标记为核心点,并将其密度直接可达的所有点归为同一簇;
  3. 若该点不是核心点,则继续选择下一个未被访问的数据点;
  4. 重复以上步骤,直到所有的数据点都被访问过。

DBSCAN算法具有以下优点:不需要预先指定簇的数量;能够发现任意形状的簇;对噪声点具有鲁棒性。缺点,对于高维数据集,由于维度灾难的影响,DBSCAN的性能可能会下降。

4 比较聚类算法

代码:

# 文件功能:
# 1. 生成模拟聚类数据
# 2. 调用聚类算法,包括自编的K-Means、GMM,以及sklearn中自带的算法,进行聚类提取
# 3. 在同一张图片中显示聚类结果和算法运行时间,对比查看

import time
import warnings

import numpy as np
import matplotlib.pyplot as plt

from sklearn import cluster, datasets, mixture
from sklearn.neighbors import kneighbors_graph
from sklearn.preprocessing import StandardScaler
from itertools import cycle, islice

from KMeans import K_Means
from GMM import GMM

np.random.seed(0)

# ============
# 模拟原始数据
# ============
print('start generate datasets ...')
n_samples = 1500
noisy_circles = datasets.make_circles(n_samples=n_samples, factor=.5,
                                      noise=.05)
noisy_moons = datasets.make_moons(n_samples=n_samples, noise=.05)
blobs = datasets.make_blobs(n_samples=n_samples, random_state=8)
no_structure = np.random.rand(n_samples, 2), None

# Anisotropicly distributed data
random_state = 170
X, y = datasets.make_blobs(n_samples=n_samples, random_state=random_state)
transformation = [[0.6, -0.6], [-0.4, 0.8]]
X_aniso = np.dot(X, transformation)
aniso = (X_aniso, y)

# blobs with varied variances
varied = datasets.make_blobs(n_samples=n_samples,
                             cluster_std=[1.0, 2.5, 0.5],
                             random_state=random_state)
print('datasets generated over')
# ============
# 设置聚类算法参数
# ============
plt.figure(figsize=(9 * 2 + 3, 12.5))
plt.subplots_adjust(left=.02, right=.98, bottom=.001, top=.96, wspace=.05,
                    hspace=.01)

plot_num = 1

default_base = {'quantile': .3,
                'eps': .3,
                'damping': .9,
                'preference': -200,
                'n_neighbors': 10,
                'n_clusters': 3,
                'min_samples': 20,
                'xi': 0.05,
                'min_cluster_size': 0.1}

datasets = [
    (noisy_circles, {'damping': .77, 'preference': -240,
                     'quantile': .2, 'n_clusters': 2,
                     'min_samples': 20, 'xi': 0.25}),
    (noisy_moons, {'damping': .75, 'preference': -220, 'n_clusters': 2}),
    (varied, {'eps': .18, 'n_neighbors': 2,
              'min_samples': 5, 'xi': 0.035, 'min_cluster_size': .2}),
    (aniso, {'eps': .15, 'n_neighbors': 2,
             'min_samples': 20, 'xi': 0.1, 'min_cluster_size': .2}),
    (blobs, {}),
    (no_structure, {})]


# 一共两层循环,实现每个仿真数据被每种聚类算法调用以此
# 此处是外层循环,遍历所有仿真数据
for i_dataset, (dataset, algo_params) in enumerate(datasets):
    # update parameters with dataset-specific values
    params = default_base.copy()
    params.update(algo_params)
    print('start deal dataset_' + str(i_dataset) + '...')
    print('dataset params is:')
    print(params)

    X, y = dataset

    # normalize dataset for easier parameter selection
    X = StandardScaler().fit_transform(X)

    # estimate bandwidth for mean shift
    bandwidth = cluster.estimate_bandwidth(X, quantile=params['quantile'])

    # connectivity matrix for structured Ward
    connectivity = kneighbors_graph(
        X, n_neighbors=params['n_neighbors'], include_self=False)
    # make connectivity symmetric
    connectivity = 0.5 * (connectivity + connectivity.T)

    # ============
    # 初始化所有聚类算法
    # ============
    # 自编的K-Means、GMM算法
    my_kmeans = K_Means(n_clusters=params['n_clusters'])
    my_gmm = GMM(n_clusters=params['n_clusters'])
    # sklearn中自带的算法
    ms = cluster.MeanShift(bandwidth=bandwidth, bin_seeding=True)
    two_means = cluster.MiniBatchKMeans(n_clusters=params['n_clusters'])
    ward = cluster.AgglomerativeClustering(
        n_clusters=params['n_clusters'], linkage='ward',
        connectivity=connectivity)
    spectral = cluster.SpectralClustering(
        n_clusters=params['n_clusters'], eigen_solver='arpack',
        affinity="nearest_neighbors")
    dbscan = cluster.DBSCAN(eps=params['eps'])
    optics = cluster.OPTICS(min_samples=params['min_samples'],
                            xi=params['xi'],
                            min_cluster_size=params['min_cluster_size'])
    affinity_propagation = cluster.AffinityPropagation(
        damping=params['damping'], preference=params['preference'])
    average_linkage = cluster.AgglomerativeClustering(
        linkage="average", affinity="cityblock",
        n_clusters=params['n_clusters'], connectivity=connectivity)
    birch = cluster.Birch(n_clusters=params['n_clusters'])
    gmm = mixture.GaussianMixture(
        n_components=params['n_clusters'], covariance_type='full')
    
    clustering_algorithms = (
        # ('My_KMeans', my_kmeans),
        # ('My_GMM', my_gmm),
        ('MiniBatchKMeans', two_means),
        ('AffinityPropagation', affinity_propagation),
        ('MeanShift', ms),
        ('SpectralClustering', spectral),
        ('Ward', ward),
        ('AgglomerativeClustering', average_linkage),
        ('DBSCAN', dbscan),
        ('OPTICS', optics),
        ('Birch', birch),
        ('GaussianMixture', gmm)
    )

    # 此处是内层循环,遍历每种算法
    for name, algorithm in clustering_algorithms:
        print(name, algorithm)
        t0 = time.time()

        print('clustering with algorithm', name, '...')
        # catch warnings related to kneighbors_graph
        with warnings.catch_warnings():
            warnings.filterwarnings(
                "ignore",
                message="the number of connected components of the " +
                "connectivity matrix is [0-9]{1,2}" +
                " > 1. Completing it to avoid stopping the tree early.",
                category=UserWarning)
            warnings.filterwarnings(
                "ignore",
                message="Graph is not fully connected, spectral embedding" +
                " may not work as expected.",
                category=UserWarning)
            algorithm.fit(X)

        t1 = time.time()
        if hasattr(algorithm, 'labels_'):
            y_pred = algorithm.labels_.astype(np.int)
        else:
            y_pred = algorithm.predict(X)

        plt.subplot(len(datasets), len(clustering_algorithms), plot_num)
        if i_dataset == 0:
            plt.title(name, size=7)

        print(y_pred)
        colors = np.array(list(islice(cycle(['#377eb8', '#ff7f00', '#4daf4a',
                                             '#f781bf', '#a65628', '#984ea3',
                                             '#999999', '#e41a1c', '#dede00']),
                                      int(max(y_pred) + 1))))
        # add black color for outliers (if any)
        colors = np.append(colors, ["#000000"])
        plt.scatter(X[:, 0], X[:, 1], s=10, color=colors[y_pred])

        plt.xlim(-2.5, 2.5)
        plt.ylim(-2.5, 2.5)
        plt.xticks(())
        plt.yticks(())
        plt.text(.99, .01, ('%.2fs' % (t1 - t0)).lstrip('0'),
                 transform=plt.gca().transAxes, size=15,
                 horizontalalignment='right')
        plot_num += 1

plt.show()

运行结果:

图2 聚类算法比较

参考

K-means:

https://blog.csdn.net/qinlele1994/article/details/106180756/

GMM:

https://blog.csdn.net/weixin_43384504/article/details/105754032

  • 0
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值