一、模糊
“模糊”:一个元素可以不同程度的属于某几个子集,也就是说元素对于集合的隶属度可以在[0,1]上取连续值。
二、步骤
2.1步骤
翻译一下:
S1:初始化参数:加权指数m,簇心数目C,以及迭代停止阈值ε。
S2:随机初始化隶属度矩阵U,注意满足式(2-2)。
S3:式(2-3)更新簇心c。
S4:式(2-4)更新隶属度矩阵U
S5:如果隶属度矩阵U满足式(2-5)则返回U并结束算法,否则转到S2
各种式如下:
2.2流程图
流程图比步骤多的一项是计算目标函数J,这一步可有可无,为了观察目标函数的变化趋势,我们在流程图和代码里都加上了计算目标函数这一步。
三、代码
import numpy as np
import matplotlib.pyplot as plt
import time
star = time.time() # 计时
img = plt.imread('1.jpg') # 读取图片信息,存储在一个三维数组中
row = img.shape[0]
col = img.shape[1]
plt.figure(1)
plt.subplot(221)
plt.imshow(img)
def fcm(data, threshold, k, m):
# 0.初始化
data = data.reshape(-1, 3)
cluster_center = np.zeros([k, 3]) # 簇心
distance = np.zeros([k, row*col]) # 欧氏距离
times = 0 # 迭代次数
goal_j = np.array([]) # 迭代终止条件:目标函数
goal_u = np.array([]) # 迭代终止条件:隶属度矩阵元素最大变化量
# 1.初始化U
u = np.random.dirichlet(np.ones(k), row*col).T # 形状(k, col*rol),任意一列元素和=1
# for s in range(50):
while 1:
times += 1
print('循环:', times)
# 2.簇心update
for i in range(k):
cluster_center[i] = np.sum((np.tile(u[i] ** m, (3, 1))).T * data, axis=0) / np.sum(u[i] ** m)
# 3.U update
# 3.1欧拉距离
for i in range(k):
distance[i] = np.sqrt(np.sum((data - np.tile(cluster_center[i], (row * col, 1))) ** 2, axis=1))
# 3.2目标函数
goal_j = np.append(goal_j, np.sum((u**m)*distance**2))
# 3.3 更新隶属度矩阵
oldu = u.copy() # 记录上一次隶属度矩阵
u = np.zeros([k, row * col])
for i in range(k):
for j in range(k):
u[i] += (distance[i] / distance[j]) ** (2 / (m - 1))
u[i] = 1/u[i]
goal_u = np.append(goal_u, np.max(u - oldu)) # 隶属度元素最大变化量
print('隶属度元素最大变化量', np.max(u - oldu), '目标函数', np.sum((u**m)*distance**2))
# 4.判断:隶属度矩阵元素最大变化量是否小于阈值
if np.max(u - oldu) <= threshold:
break
return u, goal_j, goal_u
if __name__ == '__main__':
img_show, goal1_j, goal2_u = fcm(img, 1e-09, 5, 2)
img_show = np.argmax(img_show, axis=0)
# plt.figure(2)
plt.subplot(223)
plt.plot(goal1_j)
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
plt.title('目标函数变化曲线')
plt.xlabel('迭代次数')
plt.ylabel('目标函数')
# plt.figure(3)
plt.subplot(224)
plt.plot(goal2_u)
plt.title('隶属度矩阵相邻两次迭代的元素最大变化量变化曲线')
plt.xlabel('迭代次数')
plt.ylabel('隶属度矩阵相邻两次迭代的元素最大变化量')
# plt.figure(1)
plt.subplot(222)
plt.imshow(img_show.reshape([row, col]))
end = time.time()
print('用时:', end - star)
plt.show()
四、代码怎么用?
第一步,需要把待分割图片放入项目文件,like this:
第二步 ,代码的第5行:'img = plt.imread()',在括号里加入图片路径:“图片名.文件格式”,比如,1.jpeg、2.bmp等等。
五、运行结果
1.迭代终止条件的阈值:1e-09;簇总数:5;m=2,迭代了114次
2. 迭代终止条件:1e-09;簇总数:2;m=2,迭代了16次
3.迭代终止条件:1e-09;簇总数:5;m=2,迭代了170次,用时2306s
六、终止条件为目标函数变化值小于阈值:
6.1步骤
S1:初始化参数:加权指数m,簇心数目C,以及迭代停止阈值ε 。
S2:随机初始化隶属度矩阵U,注意满足式(2-2)。
S3:式(2-3)更新簇心c。
S4:式(2-1)计算目标函数,如果变化值小于阈值,则返回U并结束算法,否则转到S2
S5:式(2-4)更新隶属度矩阵U
6.2 流程图
七、代码详解
7.1产生和为1的随机数
7.1.1德雷克分布
产生(m ,n)的二维数组,每行元素和=1
import numpy as np
x = np.random.dirichlet(np.ones(n), size=m)
7.1.2返回值
形状:(1, 20),二维数组。
转换为向量:x=x[0]
类型:ndarray
7.2三维数组
7.2.1三维和二维
形状(k, m, n)的三维数组,相对于k个形状(m, n)的二维数组
7.2.2索引二维数组
三维数组arr,形状(k,m,n)。arr第1个二维数组为arr[0]
7.3未知形状数组加入新值
7.3.1应用场景
在某一循环内,在循环的总次数不清楚的情况下,需要记录每次循环的结果。即需要给一个不清楚长度数组加入元素。
比如,在代码中,每次循环都记录目标函数的值,但是循环的总次数每次都不确定。
7.3.2实现
先定义x是一个空数组,再通过x.append()加入新值。
x = np.array([])
for i in range(n):
…
x.append(data)
还有一些关于语法的解释,可以看下面这篇,虽然是kmeans算法但是编程上原理差不多。