diff --git a/plugins/In_one/main.py b/plugins/In_one/main.py index c2a0c7c..4d09d7c 100644 --- a/plugins/In_one/main.py +++ b/plugins/In_one/main.py @@ -1,4 +1,6 @@ from asyncio.windows_events import NULL +from concurrent.futures import thread +from email.policy import default import os import pdb from threading import Thread @@ -6,15 +8,18 @@ import numpy as np from plugins.basic_change.main import MyDialog from rscder.gui.actions import ActionManager from rscder.plugins.basic import BasicPlugin -from PyQt5.QtWidgets import QAction, QDialog, QHBoxLayout, QVBoxLayout, QPushButton,QWidget,QLabel,QLineEdit,QPushButton +from PyQt5.QtWidgets import QAction, QDialog, QHBoxLayout, QVBoxLayout, QPushButton,QWidget,QLabel,QLineEdit,QPushButton,QComboBox from PyQt5.QtGui import QIcon,QPixmap from PyQt5.QtCore import Qt -from rscder.gui.layercombox import RasterLayerCombox +from rscder.gui.layercombox import LayerCombox,PairLayerCombox from rscder.utils.icons import IconInstance -from rscder.utils.project import Project, RasterLayer, SingleBandRasterLayer -from threshold.otsu import OTSU +from rscder.utils.project import Project, RasterLayer, SingleBandRasterLayer,ResultPointLayer +from In_one.otsu import OTSU from osgeo import gdal from plugins.In_one import pic +import math +from skimage.filters import rank +from skimage.morphology import disk, rectangle class LockerButton(QPushButton): def __init__(self,parent=NULL): super(LockerButton,self).__init__(parent) @@ -39,29 +44,50 @@ class LockerButton(QPushButton): def SetTextLabel(self, text): self.m_textLabel.setText(text) +class selectCombox(QComboBox): + def __init__(self, parent,list:list,default='--') : + super(selectCombox,self).__init__(parent) + self.choose=None + self.list=list + self.default=default + self.addItem(default, None) + self.addItems(list) + self.currentIndexChanged.connect(self.on_change) + + def on_change(self,index): + if index == 0: + self.choose=self.default + else: + self.choose=self.list[index-1] + # print(self.choose) + class AllInOne(QDialog): - def __init__(self, parent=None): + def __init__(self, pre,cd,threshold,parent=None): super(AllInOne, self).__init__(parent) self.setWindowTitle('变化检测') self.setWindowIcon(IconInstance().LOGO) + self.pre=pre#['均值滤波','test滤波'] + self.cd=cd#['差分法','test法'] + self.threshold=threshold#['OTSU'] self.initUI() def initUI(self): - + #图层 + self.layer_combox = PairLayerCombox(self) + layerbox = QHBoxLayout() + layerbox.addWidget(self.layer_combox) + #预处理 filterWeight=QWidget(self) filterlayout=QVBoxLayout() - filerButton =LockerButton(filterWeight); - filerButton.setObjectName("LockerButton") - filerButton.SetTextLabel("大小") - filerButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头 下.svg')) - filerButton.setStyleSheet("#LockerButton{background-color:transparent;border:none;}" - "#LockerButton:hover{background-color:rgba(195,195,195,0.4);border:none;}") - self.layer_combox = RasterLayerCombox(self) - layer_label = QLabel('图层:') - hbox = QHBoxLayout() - hbox.addWidget(layer_label) - hbox.addWidget(self.layer_combox) + filerButton =LockerButton(filterWeight) + filerButton.setObjectName("filerButton") + filerButton.SetTextLabel("预处理") + filerButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png')) + filerButton.setStyleSheet("#filerButton{background-color:transparent;border:none;}" + "#filerButton:hover{background-color:rgba(195,195,195,0.4);border:none;}") + self.pre_select=selectCombox(self,self.pre) + x_size_input = QLineEdit(self) x_size_input.setText('3') y_size_input = QLineEdit(self) @@ -78,35 +104,101 @@ class AllInOne(QDialog): hlayout1.addWidget(time_label) hlayout1.addWidget(y_size_input) vlayout = QVBoxLayout() - # vlayout.addWidget(filerButton) - vlayout.addLayout(hbox) + vlayout.addWidget(self.pre_select) vlayout.addLayout(hlayout1) filterWeight.setLayout(vlayout) filterlayout.addWidget(filerButton) filterlayout.addWidget(filterWeight) #变化检测 + changelayout=QVBoxLayout() changeWeight=QWidget(self) + changeButton =LockerButton(changeWeight) + changeButton.setObjectName("changeButton") + changeButton.SetTextLabel("变化检测") + changeButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png')) + changeButton.setStyleSheet("#changeButton{background-color:transparent;border:none;}" + "#changeButton:hover{background-color:rgba(195,195,195,0.4);border:none;}") + changeWeightlayout=QVBoxLayout() + self.cd_select=selectCombox(self,self.cd) + changeWeightlayout.addWidget(self.cd_select) + changeWeight.setLayout(changeWeightlayout) + changelayout.addWidget(changeButton) + changelayout.addWidget(changeWeight) + #阈值处理 + thresholdlayout=QVBoxLayout() + thresholdWeight=QWidget(self) + thresholdButton =LockerButton(thresholdWeight) + thresholdButton.setObjectName("thresholdButton") + thresholdButton.SetTextLabel("阈值处理") + thresholdButton.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png')) + thresholdButton.setStyleSheet("#thresholdButton{background-color:transparent;border:none;}" + "#thresholdButton:hover{background-color:rgba(195,195,195,0.4);border:none;}") + self.threshold_select=selectCombox(self,self.threshold,default='手动阈值') + self.threshold_input=QLineEdit(self) + self.threshold_input.setText('0.5') + self.threshold_select.currentIndexChanged.connect(lambda index:self.hide_(self.threshold_input,index==0)) + thresholdWeightlayout=QVBoxLayout() + thresholdWeightlayout.addWidget(self.threshold_select) + thresholdWeightlayout.addWidget(self.threshold_input) + + thresholdWeight.setLayout(thresholdWeightlayout) + thresholdlayout.addWidget(thresholdButton) + thresholdlayout.addWidget(thresholdWeight) + + #确认 + oklayout=QHBoxLayout() + self.ok_button = QPushButton('确定', self) + self.ok_button.setIcon(IconInstance().OK) + self.ok_button.clicked.connect(self.accept) + + self.cancel_button = QPushButton('取消', self) + self.cancel_button.setIcon(IconInstance().CANCEL) + self.cancel_button.clicked.connect(self.reject) + oklayout.addWidget(self.ok_button) + oklayout.addWidget(self.cancel_button) totalvlayout=QVBoxLayout() + totalvlayout.addLayout(layerbox) totalvlayout.addLayout(filterlayout) - + totalvlayout.addLayout(changelayout) + totalvlayout.addLayout(thresholdlayout) + totalvlayout.addLayout(oklayout) totalvlayout.addStretch() self.setLayout(totalvlayout) filerButton.clicked.connect(lambda: self.hide(filerButton,filterWeight)) + changeButton.clicked.connect(lambda: self.hide(changeButton,changeWeight)) + thresholdButton.clicked.connect(lambda: self.hide(thresholdButton,thresholdWeight)) + + def hide(self,button:LockerButton,weight:QWidget): if ((button.hide_)%2)==1: weight.setVisible(False) - button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头 下.svg')) + button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表向右.png')) else: weight.setVisible(True) - button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头 右.svg')) + button.SetImageLabel(QPixmap('plugins/In_one/pic/箭头_列表展开.png')) button.hide_=(button.hide_)%2+1 + def hide_(self,widget:QWidget,h:bool): + if h: + widget.setVisible(True) + else: + widget.setVisible(False) + def hideWidget(self,widget:QWidget): + if widget.isVisible: + widget.setVisible(False) + else: + widget.setVisible(True) class InOnePlugin(BasicPlugin): + pre=['均值滤波'] + cd=['差分法'] + threshold=['OTSU阈值'] + + @staticmethod def info(): return { @@ -118,12 +210,287 @@ class InOnePlugin(BasicPlugin): def set_action(self): - basic_diff_method_in_one = QAction('inone差分法') + basic_diff_method_in_one = QAction('总流程') ActionManager().change_detection_menu.addAction(basic_diff_method_in_one) self.basic_diff_method_in_one = basic_diff_method_in_one basic_diff_method_in_one.triggered.connect(self.run) def run(self): - myDialog=AllInOne(self.mainwindow) - myDialog.show() \ No newline at end of file + myDialog=AllInOne(self.pre,self.cd,self.threshold,self.mainwindow) + myDialog.show() + if myDialog.exec_()==QDialog.Accepted: + t=Thread(target=self.run_alg,args=(myDialog,)) + t.start() + + def run_alg(self,w:AllInOne): + layer1=w.layer_combox.layer1 + layer2=w.layer_combox.layer2 + if not layer1.compare(layer2): + self.send_message.emit('两个图层的尺寸不同') + return + pth1 = w.layer_combox.layer1.path + pth2 = w.layer_combox.layer2.path + name=layer1.layer_parent.name + if w.pre_select.choose==self.pre[0]: + + pth1=Meanfilter(w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer1) + self.send_message.emit('均值滤波图像{}'.format(w.layer_combox.layer1.name)) + pth2=Meanfilter(w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer2) + self.send_message.emit('均值滤波图像{}'.format(w.layer_combox.layer2.name)) + name=name+'_mean_filter' + else: + pass + + cdpth=None + if w.cd_select.choose==self.cd[0]: + cdpth=basic_cd(pth1,pth2,w.layer_combox.layer1.layer_parent,self.send_message) + name += '_basic_cd' + else: + pass + + thpth=None + if w.threshold_select.choose==self.threshold[0]: + thpth=otsu(cdpth,w.layer_combox.layer1.layer_parent.name,self.send_message) + name+='_otsu' + elif w.threshold_select.choose=='手动阈值': + thpth=thresh(cdpth,float(w.threshold_input.text()),w.layer_combox.layer1.layer_parent.name,self.send_message) + else: + pass + + table_layer(thpth,layer1,name,self.send_message) + +def Meanfilter(x_size,y_size,layer:RasterLayer): + x_size = int(x_size) + y_size = int(y_size) + pth = layer.path + if pth is None: + return + + ds = gdal.Open(pth) + band_count = ds.RasterCount + out_path = os.path.join(Project().other_path, '{}_mean_filter.tif'.format(layer.name)) + out_ds = gdal.GetDriverByName('GTiff').Create(out_path, ds.RasterXSize, ds.RasterYSize, band_count, ds.GetRasterBand(1).DataType) + out_ds.SetProjection(ds.GetProjection()) + out_ds.SetGeoTransform(ds.GetGeoTransform()) + + for i in range(band_count): + band = ds.GetRasterBand(i+1) + data = band.ReadAsArray() + + data = rank.mean(data, rectangle(y_size, x_size)) + + out_band = out_ds.GetRasterBand(i+1) + out_band.WriteArray(data) + + out_ds.FlushCache() + del out_ds + del ds + return out_path + +def basic_cd(pth1,pth2,layer_parent,send_message): + + ds1 = gdal.Open(pth1) + ds2 = gdal.Open(pth2) + cell_size = layer_parent.cell_size + xsize = ds1.RasterXSize + ysize = ds1.RasterYSize + band = ds1.RasterCount + yblocks = ysize // cell_size[1] + + driver = gdal.GetDriverByName('GTiff') + out_tif = os.path.join(Project().cmi_path, 'temp.tif') + out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) + out_ds.SetGeoTransform(ds1.GetGeoTransform()) + out_ds.SetProjection(ds1.GetProjection()) + max_diff = 0 + min_diff = math.inf + for j in range(yblocks + 1): + + send_message.emit(f'计算{j}/{yblocks}') + block_xy = (0, j * cell_size[1]) + if block_xy[1] > ysize: + break + block_size = (xsize, cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + block_data1 = ds1.ReadAsArray(*block_xy, *block_size) + block_data2 = ds2.ReadAsArray(*block_xy, *block_size) + + if band == 1: + block_data1 = block_data1[None, ...] + block_data2 = block_data2[None, ...] + # pdb.set_trace() + block_diff = block_data1.sum(0) - block_data2.sum(0) + block_diff = block_diff.astype(np.float32) + block_diff = np.abs(block_diff) + + min_diff = min(min_diff, block_diff[block_diff > 0].min()) + max_diff = max(max_diff, block_diff.max()) + out_ds.GetRasterBand(1).WriteArray(block_diff, *block_xy) + + send_message.emit(f'完成{j}/{yblocks}') + + out_ds.FlushCache() + del out_ds + send_message.emit('归一化概率中...') + temp_in_ds = gdal.Open(out_tif) + + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + out_normal_ds = driver.Create(out_normal_tif, xsize, ysize, 1, gdal.GDT_Byte) + out_normal_ds.SetGeoTransform(ds1.GetGeoTransform()) + out_normal_ds.SetProjection(ds1.GetProjection()) + # hist = np.zeros(256, dtype=np.int32) + for j in range(yblocks+1): + block_xy = (0, j * cell_size[1]) + if block_xy[1] > ysize: + break + block_size = (xsize, cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + block_data = temp_in_ds.ReadAsArray(*block_xy, *block_size) + block_data = (block_data - min_diff) / (max_diff - min_diff) * 255 + block_data = block_data.astype(np.uint8) + out_normal_ds.GetRasterBand(1).WriteArray(block_data, *block_xy) + # hist_t, _ = np.histogram(block_data, bins=256, range=(0, 256)) + # hist += hist_t + # print(hist) + del temp_in_ds + del out_normal_ds + try: + os.remove(out_tif) + except: + pass + + # raster_result_layer = SingleBandRasterLayer(None, True, out_normal_tif, BasicLayer.BOATH_VIEW) + + # layer1.layer_parent.add_result_layer(point_result_lalyer) + # layer1.layer_parent.add_result_layer(raster_result_layer) + + # self.send_message.emit('完成计算变化表格') + + send_message.emit('差分法计算完成') + return out_normal_tif + +def otsu(pth,name,send_message): + ds = gdal.Open(pth) + band = ds.GetRasterBand(1) + # band_count = ds.RasterCount + + hist = np.zeros(256, dtype=np.int) + xsize = ds.RasterXSize + ysize = ds.RasterYSize + + max_pixels = 1e7 + max_rows = max_pixels // xsize + if max_rows < 1: + max_rows = 1 + if max_rows > ysize: + max_rows = ysize + block_count = ysize // max_rows + 1 + for i in range(block_count): + start_row = i * max_rows + end_row = min((i + 1) * max_rows, ysize) + block = band.ReadAsArray(0, start_row, xsize, end_row - start_row) + hist += np.histogram(block.flatten(), bins=256, range=(0, 255))[0] + hist = hist.astype(np.float32) + gap = OTSU(hist) + send_message.emit('阈值为：{}'.format(gap)) + + out_th = os.path.join(Project().bcdm_path, '{}_otsu_bcdm.tif'.format(name)) + out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte) + out_ds.SetGeoTransform(ds.GetGeoTransform()) + out_ds.SetProjection(ds.GetProjection()) + out_band = out_ds.GetRasterBand(1) + + for i in range(block_count): + start_row = i * max_rows + end_row = min((i + 1) * max_rows, ysize) + block = band.ReadAsArray(0, start_row, xsize, end_row - start_row) + out_band.WriteArray(block > gap, 0, start_row) + out_band.FlushCache() + out_ds = None + ds = None + send_message.emit('OTSU阈值完成') + return out_th + #otsu_layer = SingleBandRasterLayer(path = out_th, style_info={}) + #layer.layer_parent.add_result_layer(otsu_layer) + +def thresh(pth,gap,name,send_message): + ds = gdal.Open(pth) + band = ds.GetRasterBand(1) + # band_count = ds.RasterCount + + + xsize = ds.RasterXSize + ysize = ds.RasterYSize + + max_pixels = 1e7 + max_rows = max_pixels // xsize + if max_rows < 1: + max_rows = 1 + if max_rows > ysize: + max_rows = ysize + block_count = ysize // max_rows + 1 + # for i in range(block_count): + # start_row = i * max_rows + # end_row = min((i + 1) * max_rows, ysize) + # block = band.ReadAsArray(0, start_row, xsize, end_row - start_row) + # hist += np.histogram(block.flatten(), bins=256, range=(0, 255))[0] + # hist = hist.astype(np.float32) + send_message.emit('阈值为：{}'.format(gap)) + + out_th = os.path.join(Project().bcdm_path, '{}_thresh{}_bcdm.tif'.format(name,gap)) + out_ds = gdal.GetDriverByName('GTiff').Create(out_th, xsize, ysize, 1, gdal.GDT_Byte) + out_ds.SetGeoTransform(ds.GetGeoTransform()) + out_ds.SetProjection(ds.GetProjection()) + out_band = out_ds.GetRasterBand(1) + + for i in range(block_count): + start_row = i * max_rows + end_row = min((i + 1) * max_rows, ysize) + block = band.ReadAsArray(0, start_row, xsize, end_row - start_row) + out_band.WriteArray(block > gap, 0, start_row) + out_band.FlushCache() + out_ds = None + ds = None + send_message.emit('自定义阈值分割完成') + return out_th + #otsu_layer = SingleBandRasterLayer(path = out_th, style_info={}) + #layer.layer_parent.add_result_layer(otsu_layer) + +def table_layer(pth,layer,name,send_message): + send_message.emit('正在计算表格结果...') + cell_size = layer.layer_parent.cell_size + ds = gdal.Open(pth) + xsize = ds.RasterXSize + ysize = ds.RasterYSize + geo = ds.GetGeoTransform() + + out_csv = os.path.join(Project().other_path, f'{name}_table_result.csv') + yblocks = ysize // cell_size[1] + 1 + xblocks = xsize // cell_size[0] + 1 + with open(out_csv, 'w') as f: + f.write('x,y,diff,status\n') + for j in range(yblocks): + block_xy = (0, j * cell_size[1]) + block_size = (xsize, cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + block_data = ds.ReadAsArray(*block_xy, *block_size) + for i in range(xblocks): + start_x = i * cell_size[0] + end_x = start_x + cell_size[0] + if end_x > xsize: + end_x = xsize + block_data_xy = block_data[:, start_x:end_x] + if block_data_xy.mean() > 0.5: + center_x = start_x + cell_size[0] // 2 + center_y = j * cell_size[1] + cell_size[1] // 2 + center_x = center_x * geo[1] + geo [0] + center_y = center_y * geo[5] + geo [3] + f.write(f'{center_x},{center_y},{block_data_xy.mean() * 100},1\n') + + result_layer = ResultPointLayer(out_csv, enable=True, proj=layer.proj, geo=layer.geo) + layer.layer_parent.add_result_layer(result_layer) + send_message.emit('计算完成') \ No newline at end of file diff --git a/plugins/In_one/otsu.py b/plugins/In_one/otsu.py new file mode 100644 index 0000000..962075f --- /dev/null +++ b/plugins/In_one/otsu.py @@ -0,0 +1,42 @@ +import numpy as np + +def OTSU(hist): + + u1=0.0#背景像素的平均灰度值 + u2=0.0#前景像素的平均灰度值 + th=0.0 + + #总的像素数目 + PixSum= np.sum(hist) + #各灰度值所占总像素数的比例 + PixRate=hist / PixSum + #统计各个灰度值的像素个数 + Max_var = 0 + #确定最大类间方差对应的阈值 + GrayScale = len(hist) + for i in range(1,len(hist)):#从1开始是为了避免w1为0. + u1_tem=0.0 + u2_tem=0.0 + #背景像素的比列 + w1=np.sum(PixRate[:i]) + #前景像素的比例 + w2=1.0-w1 + if w1==0 or w2==0: + pass + else:#背景像素的平均灰度值 + for m in range(i): + u1_tem=u1_tem+PixRate[m]*m + u1 = u1_tem * 1.0 / w1 + #前景像素的平均灰度值 + for n in range(i,GrayScale): + u2_tem = u2_tem + PixRate[n]*n + u2 = u2_tem / w2 + #print(u1) + #类间方差公式：G=w1*w2*(u1-u2)**2 + tem_var=w1*w2*np.power((u1-u2),2) + #print(tem_var) + #判断当前类间方差是否为最大值。 + if Max_var