基于模糊C--均值聚类彩色图像分割

import numpy as np
import matplotlib.pyplot as plt
import cv2 as cv
a=[]
b=[]
c=[]
# 计算N个样本，到C个中心的距离
#  X ： [N,D],D=2
#  Centers :  [C,D],D=2
# 返回 [N,C] 两两之间的距离
def FCM_dist(X, Centers):
    N, D = np.shape(X)#N个样本点（N，2）
    C, D = np.shape(Centers)#中心点（2，2）第一行，第一个中心点（x，y）

    tile_x = np.tile(np.expand_dims(X, axis=1), [1, C, 1])#ND--N1D--NCD
    tile_centers = np.tile(np.expand_dims(Centers, axis=0), [N, 1, 1])#CD--1CD--NCD

    dist = np.sum((tile_x - tile_centers) ** 2, axis=-1)#NCD-NCD,变成NC


    return np.sqrt(dist)


# 获取新的聚类中心
# U 关系矩阵 [N,C]
# X 输入数据 [N,D]
# 返回新的聚类中心 [C,D]

def FCM_getCenters(U, X, m):
    N, D = np.shape(X)
    N, C = np.shape(U)

    um = U ** m

    tile_X = np.tile(np.expand_dims(X, 1), [1, C, 1])#N1D--NCD
    tile_um = np.tile(np.expand_dims(um, -1), [1, 1, D])#NC1--NCD
    temp = tile_X * tile_um
    new_C = np.sum(temp, axis=0) / np.expand_dims(np.sum(um, axis=0), axis=-1)
    #np.sum(temp, axis=0)--CD
    #np.expand_dims(np.sum(um, axis=0), axis=-1)--N,C--c--C,1;系统除法时，会出现广播效应

    return new_C


# 更新关系矩阵
# X : [N,D]
# Centers : [C,D]
# 返回
# U ：[N,C]

def FCM_getU(X, Centers, m):

    N, D = np.shape(X)
    C, D = np.shape(Centers)

    temp = FCM_dist(X, Centers) ** float(2 / (m - 1))  # NC

    tile_temp = np.tile(np.expand_dims(temp, 1), [1, C, 1])  # NC--NCC

    denominator_ = np.expand_dims(temp, -1) / tile_temp  # NC--NC1/NCC=NCC

    return 1 / np.sum(denominator_, axis=-1)  # NC


def FCM_train(X, n_centers, m, max_iter=100, theta=1e-5, seed=0):
    rng = np.random.RandomState(seed)
    N, D = np.shape(X)
    # 随机初始化关系矩阵
    U = rng.uniform(size=(N, n_centers))
    # 保证每行和为1
    U = U / np.sum(U, axis=1, keepdims=True)

    # 开始迭代
    for i in range(max_iter):
        # print(i)
        U_old = U.copy()
        centers = FCM_getCenters(U, X, m)
        U = FCM_getU(X, centers, m)

        # 两次关系矩阵距离过小，结束训练
        #np.linalg.norm:表示求整体的矩阵平方和，再开根号(降维)。加keepdims，保持维度
        if np.linalg.norm(U - U_old) < theta:
            break

    return centers, U

#聚类标签U每行最大
def FCM_getClass(U):
    return np.argmax(U, axis=-1)





if __name__=='__main__':
    n_centers = 3
    m = 2
    img=cv.imread('../data/1.jpg')
#cv中默认BGR，在matplotlib中需要转换为RGB进行显示图像
    img=cv.cvtColor(img,cv.COLOR_BGR2RGB)
    data=img.reshape(-1,3)#将三维降成二维（这个-1代表的是什么意思，我不
    # 知道可以分成多少行，但是我需要分成多少列）
    data1=data[:,0].reshape(-1,1)
    data2=data[:,1].reshape(-1,1)
    data3=data[:,2].reshape(-1,1)
    centers1,U1 = FCM_train(data1, n_centers, m, max_iter=100, theta=1e-5, seed=0)
    centers2, U2 = FCM_train(data2, n_centers, m, max_iter=100, theta=1e-5, seed=0)
    centers3, U3 = FCM_train(data3, n_centers, m, max_iter=100, theta=1e-5, seed=0)
    labels1 = FCM_getClass(U1)
    labels2 = FCM_getClass(U2)
    labels3 = FCM_getClass(U3)
    for i in range(len(centers1)):
        d=centers1[i][0]
        a.append(np.uint8(d))
    for i in range(len(centers2)):
        d=centers2[i][0]
        b.append(np.uint8(d))
    for i in range(len(centers3)):
        d=centers3[i][0]
        c.append(np.uint8(d))
    for i in range(len(labels1)):
        if labels1[i]==0:
            labels1[i]=a[0]
        elif labels1[i]==1:
            labels1[i]=a[1]
        elif labels1[i]==2:
            labels1[i]=a[2]
        elif labels1[i]==3:
            labels1[i]=a[3]
    for i in range(len(labels2)):
        if labels2[i]==0:
            labels2[i]=b[0]
        elif labels2[i]==1:
            labels2[i]=b[1]
        elif labels2[i]==2:
            labels2[i]=b[2]
        elif labels2[i]==3:
            labels2[i]=b[3]

    for i in range(len(labels3)):
        if labels3[i]==0:
            labels3[i]=c[0]
        elif labels3[i]==1:
            labels3[i]=c[1]
        elif labels3[i]==2:
            labels3[i]=c[2]
        elif labels3[i]==3:
            labels3[i]=c[3]

    data=np.c_[labels1,labels2,labels3]
    new_img=np.reshape(data,img.shape)
    titles = ['原始图像', 'Fuzzy-c-means聚类图像']
    images = [img, new_img]
    # 用来正常显示中文标签
    plt.rcParams['font.sans-serif'] = ['SimHei']
    for i in range(2):
        plt.subplot(1, 2, i + 1), plt.imshow(images[i]),
        plt.title(titles[i])
        plt.xticks([]), plt.yticks([])
    plt.show()