import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import os import time import random import numpy as np from skimage import io,exposure from numpy import median,mean,cov,linalg,transpose,sum class Zsj_Cluster: def __init__(self,file_name): self.file_name=file_name self.img_value = np.int64(io.imread(file_name)) self.weith=self.img_value.shape[1] self.heigh=self.img_value.shape[0] self.img_blue = self.img_value.transpose((2, 0, 1))[0] self.img_green = self.img_value.transpose((2, 0, 1))[1] self.img_red = self.img_value.transpose((2, 0, 1))[2] self.img_nir = self.img_value.transpose((2, 0, 1))[3] self.img_ndwi=(self.img_green - self.img_nir) / (self.img_green self.img_nir) self.img_nndwi = (self.img_red self.img_nir) / (self.img_green * 2) self.color=[[0,0，0]，[0，255，255]，[255，0，255]，[255，0]，[255，0]，[255，0，0]，[0，0，0，255]， [0，128，128]，[128，0，128]，[128，128，0，0，0，0，0，0，0，0，0，0，0，0，0，0，0，0，0，128], [255，128，128]，[128，255，128]，[128，128，255，255]，[255，255]，[255，255，255，255，255，255，255，255，255，255]] self.ndwi=0.09 #设定的水体指数ndwi阈值 def zsj_fenlei_ndwi(self,yuzhi,outwenjian): #fazhi1.确定水体指数ndwi分类阀值；outwenjian:保存的文件地址和名称 star=time.time() 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 >= yuzhi))[0]): new_img1[i]=[0, 0, 0] new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3)) io.imsave(outwenjian, new_img) end = time.time() print('{:.2}秒'.format(end-star)) def zsj_fenlei_nndwi(self,yuzhi,outwenjian): #fazhi1.确定水体指数ndwi分类阀值；outwenjian:保存的文件地址和名称 star=time.time() ndwi = self.img_nndwi 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 <= yuzhi)[0]): new_img1[i]=[0, 0, 0] new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3)) io.imsave(outwenjian, new_img) end = time.time() print('{:.2}秒'.format(end-star)) def zsj_fenleik_distance1(self,e1): #根据图像展开二维数组的索引值，获得距离矩阵 blue_dis = self.img_blue - (self.img_blue.reshape(self.heigh * self.weith)[e1]) green_dis = self.img_green - (self.img_green.reshape(self.heigh * self.weith)[e1]) red_dis = self.img_red - (self.img_red.reshape(self.heigh * self.weith)[e1]) nir_dis = self.img_nir - (self.img_nir.reshape(self.heigh * self.weith)[e1]) dis= np.sqrt(blue_dis**2 green_dis**2 red_dis**2 nir_dis**2) return dis def zsj_fenleik_distance2(self,e1): #根据获得的维度坐标，获得距离矩阵 blue_dis = self.img_blue - e1[0] green_dis = self.img_green - e1[1] red_dis = self.img_red - e1[2] nir_dis = self.img_nir - e1[3] dis= np.sqrt(blue_dis**2 green_dis**2 red_dis**2 nir_dis**2) return dis def ndwi_variation(self): ndwi1 = [1 for x in range(self.heigh * self.weith)] for i in np.array(np.where(self.img_ndwi.reshape(self.heigh * self.weith) >= self.ndwi)[0]): ndwi1[i] = 3600 ndwi_dis = np.array(ndwi1).reshape(self.heigh, self.weith) return ndwi_dis def zsj_fenleik_distance3(self,e1): #根据图像展开二维数组的索引值，获得距离矩阵 blue_dis = self.img_blue - (self.img_blue.reshape(self.heigh * self.weith)[e1]) green_dis = self.img_green - (self.img_green.reshape(self.heigh * self.weith)[e1]) red_dis = self.img_red - (self.img_red.reshape(self.heigh * self.weith)[e1]) nir_dis = self.img_nir - (self.img_nir.reshape(self.heigh * self.weith)[e1]) ndwi_dis=self.ndwi_variation() -self.ndwi_variation().reshape(self.heigh*self.weith)[e1] dis= np.sqrt(blue_dis**2 green_dis**2 red_dis**2 nir_dis**2 ndwi_dis**2) # dis = np.sqrt(blue_dis ** 2 green_dis ** 2 red_dis ** 2 nir_dis ** 2) ndwi_dis return dis def zsj_fenleik_distance4(self,e1): #根据获得的纬度坐标，获得距离矩阵 blue_dis = self.img_blue - e1[0] green_dis = self.img_green - e1[1] red_dis = self.img_red - e1[2] nir_dis = self.img_nir - e1[3] ndwi_dis=self.ndwi_variation() - e1[4] dis= np.sqrt(blue_dis**2 green_dis**2 red_dis**2 nir_dis**2 ndwi_dis**2) return dis def zsj_ms_dis1(self,e1): #根据获得的图像展开二维数组的索引值，获得距离矩阵 img2=self.img_value img1 = img2.reshape(self.heigh * self.weith, 4) img = np.append(img1.transpose(1, 0), [self.ndwi_variation().reshape(self.heigh * self.weith)],axis=0).transpose(1, 0) dis=[] imgT=img.T #转置 D = np.cov(mgT) #协方差
invD = np.linalg.pinv(D) # 计算逆矩阵
tp = img - img[e1]
for i in range(self.heigh*self.weith):
dis.append(np.sqrt(np.dot(np.dot(tp[i], invD), tp[i].T)))
return np.array(dis)
def zsj_ms_dis2(self,e1): #按照得到的维度坐标，获得距离矩阵
img2 = self.img_value
img1 = img2.reshape(self.heigh * self.weith, 4)
img = np.append(img1.transpose(1, 0), [self.ndwi_variation().reshape(self.heigh * self.weith)],
axis=0).transpose(1, 0)
dis = []
imgT = img.T # 转置
D = np.cov(imgT) # 协方差
invD = np.linalg.pinv(D) # 计算逆矩阵
tp = img - e1
for i in range(self.heigh * self.weith):
dis.append(np.sqrt(np.dot(np.dot(tp[i], invD), tp[i].T)))
return np.array(dis)
def zsj_ms_dis3(self,e1): #按照得到的图像展开为二维数组的索引值，获得距离矩阵
img=self.img_value.reshape(self.heigh*self.weith,4)
dis=[]
imgT=img.T #转置
D = np.cov(imgT) #协方差
invD = np.linalg.pinv(D) # 计算逆矩阵
tp = img - img[e1]
for i in range(self.heigh*self.weith):
dis.append(np.sqrt(np.dot(np.dot(tp[i], invD), tp[i].T)))
return np.array(dis)
def zsj_ms_dis4(self,e1): #按照得到的维度坐标，获得距离矩阵
img = self.img_value.reshape(self.heigh * self.weith, 4)
dis = []
imgT = img.T # 转置
D = np.cov(imgT) # 协方差
invD = np.linalg.pinv(D) # 计算逆矩阵
tp = img - e1
for i in range(self.heigh * self.weith):
dis.append(np.sqrt(np.dot(np.dot(tp[i], invD), tp[i].T)))
return np.array(dis)
def zsj_Kcenter(self,K): #k-means++算法选取K个中心点
center,dis1,i= [random.randint(0, self.heigh * self.weith - 1)],[],0
while True:
dis1.append(self.zsj_fenleik_distance1(center[i]).reshape(self.heigh * self.weith))
dis2=np.array(dis1).min(0)
center.append(np.argmax(dis2 ** 2 / sum(dis2 ** 2)))
center1=center[-K:]
i += 1
if i-2**3+1 == K: break
return center1
def zsj_Kcenter2(self,K): #k-means++算法选取K个中心点
center,dis1,i= [random.randint(0, self.heigh * self.weith - 1)],[],0
while True:
dis1.append(self.zsj_fenleik_distance3(center[i]).reshape(self.heigh * self.weith))
dis2=np.array(dis1).min(0)
center.append(np.argmax(dis2 ** 2 / sum(dis2 ** 2)))
center1=center[-K:]
i += 1
if i-2**3+1 == K: break
return center1
def zsj_Kmeansplus(self,K,jingdu,outwenjian):
img=self.img_value
img3 = img.reshape(self.heigh * self.weith, 4)
dis,center,kk=[],self.zsj_Kcenter(K),np.arange(0,K,1)
for i in range(K):
dis.append(self.zsj_fenleik_distance1(center[i]).reshape(self.heigh * self.weith))
iii,KK=0,[2**16,2**16,2**16,2**16]
while True:
dis1, ii,new_img1,dis3,fen,fen11,fen12= np.argmin(dis, axis=0), 0,[],np.copy(dis),\
[[] for i in range(K)],[[] for i in range(K)],[[] for i in range(K)]
for i1 in dis1:
new_img1.append(self.color[-1])
fen[np.where(kk==i1)[0][0]].append(img3[ii])
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
new_center,mm=[],[]
for i2 in range(K):
new_center.append(median(fen[i2],axis=0))
dis[i2] = self.zsj_fenleik_distance2(new_center[i2]).reshape(self.heigh * self.weith)
mm.append(median(fen11[i2]))
print(mean(abs(mean(dis,axis=1)-mean(dis3,axis=1))))
if np.argmax(mm) >= 0.1:
for i3 in fen12[np.argmax(mm)]:
new_img1[i3] = self.color[0]
else:
print('非水体')
new_img = np.uint8(np.array(new_img1).reshape(self.heigh, self.weith, 3))
io.imsave(outwenjian, new_img)
KK.append(mean(abs(np.mean(dis, axis=1) - mean(dis3, axis=1))))
# if mean(abs(np.mean(dis,axis=1)-mean(dis3,axis=1)))<jingdu:
if KK[-1]-KK[-4] > 0 :break
iii+=1
def zsj_Kmeansplus2(self,K,jingdu,outwenjian):
img=self.img_value
img1 = img.reshape(self.heigh * self.weith, 4)
img2 = np.append(img1.transpose(1, 0),[self.ndwi_variation().reshape(self.heigh * self.weith)],axis=0).transpose(1, 0)
dis,center,kk=[],self.zsj_Kcenter2(K),np.arange(0,K,1)
for i in range(K):
dis.append(self.zsj_fenleik_distance3(center[i]).reshape(self.heigh * self.weith))
iii=0
while True:
star=time.time()
dis1, ii,new_img1,dis3,fen,fen11,fen12= np.argmin(dis, axis=0), 0,[],np.copy(dis),\
[[] for i in range(K)],[[] for i in range(K)],[[] for i in range(K)]
for i1 in dis1:
new_img1.append(self.color[-1])
fen[np.where(kk==i1)[0][0]].append(img2[ii])
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
new_center,mm=[],[]
for i2 in range(K):
new_center.append(np.mean(fen[i2],axis=0))
dis[i2] = self.zsj_fenleik_distance4(new_center[i2]).reshape(self.heigh * self.weith)
mm.append(median(fen11[i2]))
end = time.time()
if np.argmax(mm) >= 0.05:
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)
print(mean(abs(np.mean(dis,axis=1)-mean(dis3,axis=1))),'{:.3}秒'.format(end-star))
if np.mean(abs(np.mean(dis,axis=1)-np.mean(dis3,axis=1)))<jingdu:
break
iii+=1
def zsj_Kmeansplus3(self,K,outwenjian):
img=self.img_value
img2 = img.reshape(self.heigh * self.weith, 4)
dis,center,kk=[],self.zsj_Kcenter(K),np.arange(0,K,1)
for i in range(K):
dis.append(self.zsj_ms_dis3(center[i]).reshape(self.heigh * self.weith))
iii,KK=0,[2**16,2**16,2**16,2**16]
while True:
dis1, ii,new_img1,dis3,fen,fen11,fen12= np.argmin(dis, axis=0), 0,[],np.copy(dis),\
[[] for i in range(K)],[[] for i in range(K)],[[] for i in range(K)]
for i1 in dis1:
new_img1.append(self.color[-1])
fen[np.where(kk==i1)[0][0]].append(img2[ii])
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
new_center,mm=[],[]
for i2 in range(K):
new_center.append(mean(fen[i2],axis=0))
dis[i2] = self.zsj_ms_dis4(new_center[i2]).reshape(self.heigh * self.weith)
mm.append(median(fen11[i2]))
if np.argmax(mm)>=0.09:
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)
else:
print('非水体')
break
print(mean(abs(np.mean(dis,axis=1)-mean(dis3,axis=1))))
KK.append(mean(abs(np.mean(dis, axis=1) - mean(dis3, axis=1))))
if KK[-1]-KK[-4] >= 0:break
iii+=1
# data=csv.reader("D:/BS/xuexi/1/tezhen.xlsx")
# print(data)
# for i in fen[np.argmax(mm)]:
# data.append((7, 3.14, 'good'))
# df = pd.DataFrame('tezhen.xlsx')
img=Zsj_Cluster("D:/BS/xuexi/1/cj006.tif")
# img.zsj_Kmeansplus(50,10,"D:/BS/xuexi/test/ts051.tif")
# img.zsj_Kmeansplus2(10,10,"D:/BS/xuexi/1/cj00121.tif")
# img.zsj_Kmeansplus3(6,"D:/BS/xuexi/1/cj00131.tif")
# img.zsj_fenlei_ndwi(0.02,"D:/BS/xuexi/1/cj00601.tif")
# print(np.max(img.img_red),np.max(img.img_blue),np.max(img.img_green),np.max(img.img_nir))