import random import matplotlib.pyplot as plt import numpy as np from skimage import io, exposure class Fcm_s: def __init__(self, file_name): self.file_name = file_name # self.img_value = np.int64(io.imread(self.file_name)) self.img_value = np.int64(io.imread(self.file_name).transpose((2, 0, 1)))[0] # self.img_value = np.int64(io.imread(self.file_name)).transpose((2, 0, 1))[3] # self.img_value = np.delete(np.int64(io.imread(self.file_name).transpose( # (2, 0, )，[1]，axis=0).transpose((1,2,0)) self.weith = self.img_value.shape[1] self.heigh = self.img_value.shape[0] # self.weidu=self.img_value.shape[2] self.weidu = 1 # # self.img_blue = np.int64(io.imread(self.file_name).transpose((2, 0, 1)))[0] # self.img_green = np.int64(io.imread(self.file_name).transpose((2, 0, 1)))[1] # self.img_red = np.int64(io.imread(self.file_name).transpose((2, 0, 1)))[2] # self.img_nir = np.int64(io.imread(self.file_name).transpose((2, 0, 1)))[3] # self.img_ndwi=(self.img_green - self.img_nir) / (self.img_green self.img_nir) self.color = [[0, 0, 0], [255, 255, 255], [128, 128, 128], [0, 255, 255], [255, 0, 255], [255, 255, 0], [0, 255, 0], [255, 0, 0], [0, 0, 255], [0, 128, 128], [128, 0, 128], [128, 128, 0], [128, 0, 0], [0, 128, 0], [0, 0, 128], [255, 128, 128], [128, 255, 128], [128, 128, 255], [128, 255, 255], [255, 128, 255], [255, 255, 128]] self.ndwi = 0.1 # 设置的ndwi值 self.dx = [-1, 0, 1, -1, 1, -1, 0, 1] self.dy = [-1, -1, -1, 0, 0, 1, 1, 1] # self.dx = [-2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2] # self.dy = [-2, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2] def zsj_ndwi(self): # fazhi1.确定水体指数ndwi分类阀值；outwenjian:保存的文件地址和名称 ndwi = self.img_ndwi ndwi_1 = np.array(ndwi).reshape(self.heigh * self.weith) new_img1 = [[255, 255, 255] for x in range(len(ndwi_1))] for i in np.array(np.where(ndwi_1 >= self.ndwi)[0]): new_img1[i] = [0, 0, 0] new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3)) return new_img def produce_center(self, K): a = np.ones((K, self.weidu)) for i in range(K): x, y = random.randint(0, self.heigh - 1), random.randint(0, self.weith - 1) a[i] = self.img_value[x][y] return a def init_weight(self, K): u = [] for i in range(K): u.append(random.random()) sum_u = sum(u) for i in range(len(u)): u[i] = u[i] / sum_u return u def produce_weight(self, K): weight = np.zeros((self.heigh, self.weith, K)) for i in range(self.heigh): for j in range(self.weith): miu = self.init_weight(K) for k in range(K): weight[i][j][k] = miu[k] return weight def space1(self): tj = np.zeros((self.heigh, self.weith)) for i in range(self.heigh): for j in range(self.weith): for r in range(len(self.dx)): index_x = i self.dx[r] index_y = j self.dy[r] if i self.dx[r] < 0 or i self.dx[r] >= self.heigh: index_x = i if j self.dy[r] < 0 or j self.dy[r] >= self.weith: index_y = j tj[i][j] = np.linalg.norm(self.img_value[i][j] - self.img_value[index_x][index_y]) / len(self.dx) pj = 1-(tj - np.min(tj)) / (np.max(tj) - np.min(tj)) return pj def space(self, weight, K, pj): cc = np.zeros((self.heigh, self.weith, K, len(self.dx))) for i in range(self.heigh): for j in range(self.weith): for k in range(K): for r in range(len(self.dx)): index_x = i self.dx[r] index_y = j self.dy[r] if i self.dx[r] < 0 or i self.dx[r] >= self.heigh: index_x = i if j self.dy[r] < 0 or j self.dy[r] >= self.weith: index_y = j for k in range(K): cc[i][j][k][r] = weight[i][j][k] * weight[index_x][index_y][k] * \ pj[i][j] / max([abs(self.dx[r]), abs(self.dy[r])]) return cc def cal_weight(self, K, weight, center, m=2, a=2): # 更新隶属度和聚类中心 distance = np.zeros((self.heigh, self.weith, K)) neighbor_dist = np.zeros((self.heigh, self.eith, K))
miu_fenzi = np.zeros((self.heigh, self.weith, K, self.weidu))
for i in range(self.heigh):
for j in range(self.weith):
for k in range(K):
distance[i][j][k] = pow(np.linalg.norm(self.img_value[i][j] - center[k]), 2)
# aa=[] #中值
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
neighbor_dist[i][j][k] += (1 / len(self.dx)) * pow(np.linalg.norm(
self.img_value[index_x][index_y] - center[k]), 2) # 均值
# np.append(aa,self.img_value[index_x][index_y]) #中值
# neighbor_dist[i][j][k]= pow(np.linalg.norm(np.median(np.array(aa), axis=0) - center[k]),2) #中值
miu_fenzi[i][j][k] = pow(distance[i][j][k] + a * neighbor_dist[i][j][k], 1 / (m - 1))
miu_fenmu = np.sum(miu_fenzi, axis=2)
for i in range(self.heigh):
for j in range(self.weith):
for k in range(K):
weight[i][j][k] = miu_fenzi[i][j][k][0] / (miu_fenmu[i][j][0] + 1)
return weight
def cal_center(self, K, weight, center, m=2, a=2):
b = np.sum(pow(weight, m), axis=0)
c = np.sum(b, axis=0)
weight_sum = (1 + a) * c
sum_neighbor = np.zeros((self.heigh, self.weith, self.weidu))
for i in range(self.heigh):
for j in range(self.weith):
# aa=[] #中值
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
sum_neighbor[i][j] += (self.img_value[index_x][index_y]) / len(self.dx) # 均值
# np.append(aa,self.img_value[index_x][index_y]) #中值
# sum_neighbor[i][j]=np.median(np.array(aa),axis=0) #中值
pixel_information = np.zeros((self.heigh, self.weith, K, self.weidu))
for k in range(K):
for i in range(self.heigh):
for j in range(self.weith):
pixel_information[i][j][k] += pow(weight[i][j][k], m) * \
(self.img_value[i][j] + a * sum_neighbor[i][j])
c2 = np.zeros((K, self.weidu))
for i in range(self.heigh):
for j in range(self.weith):
for k in range(K):
c2[k] += pixel_information[i][j][k]
# print('c2:',c2)
for k in range(K):
center[k] = c2[k] / weight_sum[k]
return center
def cal_weight2(self, K, weight, center, cc, m=2, a=2): # 更新隶属度和聚类中心
distance = np.zeros((self.heigh, self.weith, K))
neighbor_dist = np.zeros((self.heigh, self.weith, K))
miu_fenzi = np.zeros((self.heigh, self.weith, K, self.weidu))
for i in range(self.heigh):
for j in range(self.weith):
for k in range(K):
distance[i][j][k] = pow(np.linalg.norm(self.img_value[i][j] - center[k]), 2)
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
neighbor_dist[i][j][k] += (1 / len(self.dx)) * pow((1-cc[i][j][k][r]),2) * pow(np.linalg.norm(
self.img_value[index_x][index_y] - center[k]), 2)
miu_fenzi[i][j][k] = pow(distance[i][j][k] + a*neighbor_dist[i][j][k], 1 / (m - 1))
miu_fenmu = np.sum(miu_fenzi, axis=2)
for i in range(self.heigh):
for j in range(self.weith):
if miu_fenmu[i][j] == 0: miu_fenmu[i][j] = 1
for k in range(K):
weight[i][j][k] = miu_fenzi[i][j][k] / (miu_fenmu[i][j])
return weight
def cal_center2(self, K, weight, center, cc, m=2, a=2):
bb = np.sum(np.sum(pow(weight, 2), axis=0), axis=0)
cc1 = np.sum((1 - cc) / len(self.dx), 3)
sum_neighbor, sum_neighbor1 = np.zeros((self.heigh, self.weith, K, self.weidu)), \
np.zeros((self.heigh, self.weith, K))
for i in range(self.heigh):
for j in range(self.weith):
for k in range(K):
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
sum_neighbor[i][j][k] += (self.img_value[index_x][index_y]) * pow((1-cc[i][j][k][r]),2) / len(
self.dx)
pixel_information = np.zeros((self.heigh, self.weith, K, self.weidu))
for k in range(K):
for i in range(self.heigh):
for j in range(self.weith):
pixel_information[i][j][k] += pow(weight[i][j][k], m) * \
(self.img_value[i][j] + a * sum_neighbor[i][j][k])
pixel_information1 = pow(weight, m) * cc1
c2 = np.sum(np.sum(pixel_information, 0), 0)
c3 = np.sum(np.sum(pixel_information1, 0), 0) + bb
for k in range(K):
center[k] = c2[k] / c3[k]
return center
def separation(self, weight, K, outwenjian): # 对结果进行划分，对于同一类的像素点赋予相同的值
index, kk, new_img1, fen, fen11, fen12 = 0, np.arange(0, K, 1), [], [[] for i in range(K)], \
[[] for i in range(K)], [[] for i in range(K)]
weight1 = weight.reshape(self.heigh * self.weith, K).transpose(1, 0)
weight2 = np.argmax(weight1, axis=0)
ii = 0
ddd = self.zsj_ndwi()
for i1 in weight2:
new_img1.append(self.color[1])
fen11[np.where(kk == i1)[0][0]].append(self.img_ndwi.reshape(self.heigh * self.weith)[ii])
fen12[np.where(kk == i1)[0][0]].append(ii)
ii += 1
mm = []
for i2 in range(K):
if len(fen11[i2]) == 0:
mm.append(0)
else:
mm.append(len(np.where(np.array(fen11[i2]) > self.ndwi)[0]) / len(fen11[i2]))
ff = np.where(np.array(fen11[i2]) > self.ndwi)[0]
if np.max(mm) > 0.5:
for i3 in fen12[np.argmax(mm)]:
new_img1[i3] = self.color[0]
new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3))
io.imsave(outwenjian, new_img)
# for i in range(self.heigh):
# for j in range(self.weith):
# if(ddd[i][j][0]!=new_img[i][j][0]):
# new_img[i][j]=[255,255,255]
# io.imsave(outwenjian, new_img)
def separation2(self, weight, K, outwenjian): # 对结果进行划分，对于同一类的像素点赋予相同的值
new_img1 = []
weight1 = weight.reshape(self.heigh * self.weith, K).transpose(1, 0)
weight2 = np.argmax(weight1, axis=0)
for i1 in weight2:
new_img1.append(self.color[i1])
new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3))
io.imsave(outwenjian, new_img)
def mubiao(self, K, weight, center, cc):
aa = np.ones((self.heigh, self.weith, K))
dd = np.ones((self.heigh, self.weith, K))
ee = np.ones((K))
for k in range(K):
for i in range(self.heigh):
for j in range(self.weith):
bb = 0
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
bb += (1 - cc[i][j][k][r]) * pow(np.linalg.norm(self.img_value[index_x][index_y] - center[k]),
2) / len(self.dx)
dd[i][j][k] = weight[i][j][k] * bb
aa[j][j][k] = weight[i][j][k] * pow(np.linalg.norm(self.img_value[i][j] - center[k]), 2)
ee = np.sum(np.sum(dd + aa, 0), 0)
EE = np.sum(ee)
return EE
def mubiaoFCMS(self, K, weight, center):
aa = np.ones((self.heigh, self.weith, K))
dd = np.ones((self.heigh, self.weith, K))
ee = np.ones((K))
for k in range(K):
for i in range(self.heigh):
for j in range(self.weith):
bb = 0
for r in range(len(self.dx)):
index_x = i + self.dx[r]
index_y = j + self.dy[r]
if i + self.dx[r] < 0 or i + self.dx[r] >= self.heigh:
index_x = i
if j + self.dy[r] < 0 or j + self.dy[r] >= self.weith:
index_y = j
bb += 3 * pow(np.linalg.norm(self.img_value[index_x][index_y] - center[k]), 2) / len(self.dx)
dd[i][j][k] = weight[i][j][k] * bb
aa[j][j][k] = weight[i][j][k] * pow(np.linalg.norm(self.img_value[i][j] - center[k]), 2)
ee = np.sum(np.sum(dd + aa, 0), 0)
EE = np.sum(ee)
return EE
def main(self, T, K, outwenjian):
# center = self.produce_center(K)
center = np.array([[31], [226]])
weight = self.produce_weight(K)
print("center: \n", center)
ii = 0
while True:
# plt.plot(ii, self.mubiaoFCMS(K, weight, center), marker='+', color='coral')
# print(self.mubiaoFCMS(K, weight, center))
weight1 = np.copy(weight)
weight = self.cal_weight(K, weight, center)
center = self.cal_center(K, weight, center)
weight2 = weight
print("迭代次数:", ii + 1)
# print("center:\n", center)
# print(self.DIR(K, weight))
self.separation2(weight, K, outwenjian)
if ii + 1 == T:
self.separation2(weight, K, outwenjian)
break
ii += 1
# plt.show()
def main2(self, T, K, outwenjian):
# center = self.produce_center(K)
center = np.array([[31], [226]])
print("center: \n", center)
weight = self.produce_weight(K)
cc1 = self.space1()
cc = self.space(weight, K, cc1)
ii = 0
while True:
# plt.plot(ii, self.mubiao(K, weight, center,cc), marker='+', color='coral')
# print(self.mubiao(K, weight, center,cc))
cc = self.space(weight, K, cc1)
weight = self.cal_weight2(K, weight, center, cc)
center = self.cal_center2(K, weight, center, cc)
print("迭代次数:", ii + 1)
print("center:\n", center)
self.separation2(weight, K, outwenjian)
if ii + 1 == T:
self.separation2(weight, K, outwenjian)
break
ii += 1
# plt.show()
img = Fcm_s("D:/BS/xuexi/test/gs02.png")
# print(img.img_value)
# img.main(20,2,"D:/BS/xuexi/test/gs_fcms502.png")
img.main2(20, 2, "D:/BS/xuexi/test/gs_fcmsyh502.png")
def jingdu(data1, data2, data3, data4):
img1 = np.int64(io.imread(data1).transpose((2, 0, 1)))[0]
img2 = np.int64(io.imread(data2).transpose((2, 0, 1)))[0]
img3 = np.int64(io.imread(data3).transpose((2, 0, 1)))[0]
img4 = np.int64(io.imread(data4).transpose((2, 0, 1)))[0]
img12 = abs(img1 - img2).reshape(256 * 256)
img13 = abs(img1 - img3).reshape(256 * 256)
img14 = abs(img1 - img4).reshape(256 * 256)
PA1 = (len(np.where(img12 == 0)[0]) + len(np.where(img12 == 30)[0])) / (256 * 256)
PA2 = (len(np.where(img13 == 0)[0]) + len(np.where(img13 == 30)[0])) / (256 * 256)
PA3 = (len(np.where(img14 == 0)[0]) + len(np.where(img14 == 30)[0])) / (256 * 256)
# jindu1 = 1-len(np.where(img12 > 30)[0]) / (256 * 256)
# jindu2 = 1-len(np.where(img13 > 30)[0]) / (256 * 256)
# jindu3 = 1-len(np.where(img14 > 30)[0]) / (256 * 256)
# print(PA1)
# print(PA2)
# print(PA3)
# jingdu("D:/BS/xuexi/test/4.png","D:/BS/xuexi/test/jy_fcm02.png",
# "D:/BS/xuexi/test/jy_fcms02.png","D:/BS/xuexi/test/jy_fcmsyh02.png")