diff --git a/ECD.py b/ECD.py index 1dbf198..495d6c8 100644 --- a/ECD.py +++ b/ECD.py @@ -1,13 +1,10 @@ import os import sys -import argparse sys.path.insert(0, os.path.join(os.path.dirname(__file__), 'libs')) os.environ['PROJ_LIB'] = os.path.join(os.path.dirname(__file__), 'share/proj') os.environ['GDAL_DATA'] = os.path.join(os.path.dirname(__file__), 'share') os.environ['ECD_BASEDIR'] = os.path.dirname(__file__) BASE_DIR = os.path.dirname(__file__) -import distutils -import distutils.version from rscder import MulStart import logging diff --git a/plugins/In_one/__init__.py b/plugins/In_one/__init__.py deleted file mode 100644 index 4df2afb..0000000 --- a/plugins/In_one/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from In_one.main import * \ No newline at end of file diff --git a/plugins/In_one/main.py b/plugins/In_one/main.py deleted file mode 100644 index 8afcf9c..0000000 --- a/plugins/In_one/main.py +++ /dev/null @@ -1,555 +0,0 @@ -import os -from threading import Thread -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,QComboBox,QDialogButtonBox -from PyQt5.QtGui import QPixmap - -from rscder.gui.layercombox import PairLayerCombox -from rscder.utils.icons import IconInstance -from rscder.utils.geomath import geo2imageRC -from rscder.utils.project import Project, PairLayer,ResultPointLayer,MultiBandRasterLayer -from In_one.otsu import OTSU -from osgeo import gdal -from In_one import pic -import math -from skimage.filters import rank -from skimage.morphology import disk, rectangle -from In_one.scripts.UnsupervisedCD import LSTS,CVA,acd,aht,ocd,lhba,sh -def Meanfilter(x_size,y_size,layer:MultiBandRasterLayer): - 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:str,pth2:str,layer_parent:PairLayer,send_message): - ds1:gdal.Dataset=gdal.Open(pth1) - ds2:gdal.Dataset=gdal.Open(pth2) - - cell_size = layer_parent.cell_size - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - - 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) - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - out_ds.SetGeoTransform(geo) - out_ds.SetProjection(proj) - - max_diff = 0 - min_diff = math.inf - - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - for j in range(yblocks + 1):#该改这里了 - - send_message.emit(f'计算{j}/{yblocks}') - block_xy1 = (start1x, start1y+j * cell_size[1]) - block_xy2 = (start2x,start2y+j*cell_size[1]) - block_xy=(0,j * cell_size[1]) - if block_xy1[1] > end1y or block_xy2[1] > end2y: - break - block_size=(xsize, cell_size[1]) - block_size1 = (xsize, cell_size[1]) - block_size2 = (xsize,cell_size[1]) - if block_xy[1] + block_size[1] > ysize: - block_size = (xsize, ysize - block_xy[1]) - if block_xy1[1] + block_size1[1] > end1y: - block_size1 = (xsize,end1y - block_xy1[1]) - if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) - block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) - block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) - - 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}') - del ds2 - del ds1 - 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(geo) - out_normal_ds.SetProjection(proj) - # 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,gap - #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,dict): - 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] - - 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},{int(block_data_xy.mean() > 0.5)}\n') - - result_layer = ResultPointLayer(out_csv, enable=True, proj=layer.proj, geo=layer.geo,result_path=dict) - # print(result_layer.result_path) - layer.layer_parent.add_result_layer(result_layer) - send_message.emit('计算完成') - -class LockerButton(QPushButton): - def __init__(self,parent=None): - super(LockerButton,self).__init__(parent) - m_imageLabel = QLabel(self) - m_imageLabel.setFixedWidth(20) - m_imageLabel.setScaledContents(True) - m_imageLabel.setStyleSheet("QLabel{background-color:transparent;}") - m_textLabel = QLabel(self) - m_textLabel.setStyleSheet("QLabel{background-color:transparent;}") - self.m_imageLabel=m_imageLabel - self.m_textLabel=m_textLabel - self.hide_=1 - mainLayout = QHBoxLayout() - - mainLayout.addWidget(self.m_imageLabel) - mainLayout.addWidget(self.m_textLabel) - mainLayout.setContentsMargins(0,0,0,0) - mainLayout.setSpacing(0) - self.setLayout(mainLayout) - def SetImageLabel(self, pixmap:QPixmap): - self.m_imageLabel.setPixmap(pixmap) - 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, 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("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) - y_size_input.setText('3') - size_label = QLabel(self) - size_label.setText('窗口大小:') - time_label = QLabel(self) - time_label.setText('X') - self.x_size_input = x_size_input - self.y_size_input = y_size_input - hlayout1 = QHBoxLayout() - hlayout1.addWidget(size_label) - hlayout1.addWidget(x_size_input) - hlayout1.addWidget(time_label) - hlayout1.addWidget(y_size_input) - vlayout = QVBoxLayout() - 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) - - #确认 - - self.ok_button = QPushButton('确定', self) - self.ok_button.setIcon(IconInstance().OK) - self.ok_button.clicked.connect(self.accept) - self.ok_button.setDefault(True) - - - self.cancel_button = QPushButton('取消', self) - self.cancel_button.setIcon(IconInstance().CANCEL) - self.cancel_button.clicked.connect(self.reject) - self.cancel_button.setDefault(False) - buttonbox=QDialogButtonBox(self) - buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole) - buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole) - buttonbox.setCenterButtons(True) - #buttonbox.setContentsMargins(,) - - - totalvlayout=QVBoxLayout() - totalvlayout.addLayout(layerbox) - totalvlayout.addLayout(filterlayout) - totalvlayout.addLayout(changelayout) - totalvlayout.addLayout(thresholdlayout) - totalvlayout.addWidget(buttonbox) - 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/箭头_列表向右.png')) - else: - weight.setVisible(True) - 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) - @property - def param(self): - class param(object): - pass - p=param() - p.pre=self.pre - p.cd=self.cd - p.threshold=self.threshold - p.layer_combox=self.layer_combox - p.pre_select=self.pre_select - p.x_size_input=self.x_size_input - p.y_size_input=self.y_size_input - p.threshold_select=self.threshold_select - p.threshold_input=self.threshold_input - p.cd_select=self.cd_select - return p -class InOnePlugin(BasicPlugin): - pre={"均值滤波":Meanfilter}#可添加其他方法 - cd={'差分法':basic_cd,'LSTS':LSTS,'CVA':CVA,'ACD':acd,'AHT':aht,'OCD':ocd,'LHBA':lhba,'SH':sh}#可添加其他方法 - threshold={'OTSU阈值':otsu}#可添加其他方法 - - - @staticmethod - def info(): - return { - 'name': 'AllinOne', - 'description': 'AllinOne', - 'author': 'RSCDER', - 'version': '1.0.0', - } - - def set_action(self): - - basic_diff_method_in_one = QAction(IconInstance().UNSUPERVISED, '&无监督变化检测') - ActionManager().change_detection_menu.addAction(basic_diff_method_in_one) - # ActionManager().menubar.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(list(self.pre.keys()),list(self.cd.keys()),list(self.threshold.keys()),self.mainwindow) - myDialog.show() - if myDialog.exec_()==QDialog.Accepted: - w=myDialog.param - t=Thread(target=self.run_alg,args=(w,)) - t.start() - - def run_alg(self,w:AllInOne): - dict={} - layer1=w.layer_combox.layer1 - - - pth1 = w.layer_combox.layer1.path - pth2 = w.layer_combox.layer2.path - name=layer1.layer_parent.name - # 预处理 - # 若添加的预处理函数接口相同，则无需判断是哪种方法 - # if w.pre_select.choose==self.pre.keys()[0]: - # pass - # el - preKey=w.pre_select.choose - pth1=self.pre[preKey](w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer1) - self.send_message.emit('{}图像{}'.format(preKey,w.layer_combox.layer1.name)) - pth2=self.pre[preKey](w.x_size_input.text(),w.y_size_input.text(),w.layer_combox.layer2) - self.send_message.emit('{}图像{}'.format(preKey,w.layer_combox.layer2.name)) - name=name+'_'+preKey - dict['预处理']=[preKey,'|'.format(pth1,pth2)] - - - cdpth=None - #变化检测 - # if w.cd_select.choose==self.cd[0]: - # if w.cd_select.choose=='ACD': - # cdKey='ACD' - # self.send_message.emit('ACD计算中...') - # cdpth=os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(w.layer_combox.layer1.layer_parent.name, int(np.random.rand() * 100000))) - # print(pth1,pth2,cdpth) - # ACD(pth1,pth2,cdpth) - # else: - # pass - cdKey=w.cd_select.choose - cdpth=self.cd[cdKey](pth1,pth2,w.layer_combox.layer1.layer_parent,self.send_message) - name += '_'+cdKey - dict['变化检测算法']=[cdKey,cdpth] - - #阈值处理 - #例如手动阈值和otsu参数不同，则要做区分 - thpth=None - if w.threshold_select.choose=='手动阈值': - thpth=thresh(cdpth,float(w.threshold_input.text()),w.layer_combox.layer1.layer_parent.name,self.send_message) - dict['后处理']=['手动阈值',[float(w.threshold_input.text())],thpth] - else: - thpth,gap=otsu(cdpth,w.layer_combox.layer1.layer_parent.name,self.send_message) - name+='_otsu' - dict['后处理']=['OTSU阈值',gap,cdpth] - - table_layer(thpth,layer1,name,self.send_message,dict) - diff --git a/plugins/In_one/otsu.py b/plugins/In_one/otsu.py deleted file mode 100644 index 962075f..0000000 --- a/plugins/In_one/otsu.py +++ /dev/null @@ -1,42 +0,0 @@ -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 \ No newline at end of file diff --git a/plugins/In_one/pic/箭头 右.svg b/plugins/In_one/pic/箭头 右.svg deleted file mode 100644 index 2060636..0000000 --- a/plugins/In_one/pic/箭头 右.svg +++ /dev/null @@ -1,2 +0,0 @@ - \ No newline at end of file diff --git a/plugins/In_one/pic/箭头_列表向右.png b/plugins/In_one/pic/箭头_列表向右.png deleted file mode 100644 index d35ae74..0000000 Binary files a/plugins/In_one/pic/箭头_列表向右.png and /dev/null differ diff --git a/plugins/In_one/pic/箭头_列表展开.png b/plugins/In_one/pic/箭头_列表展开.png deleted file mode 100644 index 86a5b19..0000000 Binary files a/plugins/In_one/pic/箭头_列表展开.png and /dev/null differ diff --git a/plugins/In_one/scripts/UnsupervisedCD.py b/plugins/In_one/scripts/UnsupervisedCD.py deleted file mode 100644 index 71b99ac..0000000 --- a/plugins/In_one/scripts/UnsupervisedCD.py +++ /dev/null @@ -1,553 +0,0 @@ -from osgeo import gdal -import math,os -import time -from sklearn.cluster import k_means -from rscder.utils.geomath import geo2imageRC, imageRC2geo -from rscder.utils.project import Project, PairLayer -import numpy as np -from .ACD import ACD -from .AHT import AHT -from .OCD import OCD -from .LHBA import LHBA -from .SH import SH -def warp(file,ds:gdal.Dataset,srcWin=[0,0,0,0]): - driver = gdal.GetDriverByName('GTiff') - xsize=ds.RasterXSize - ysize=ds.RasterYSize - geo=ds.GetGeoTransform() - orj=ds.GetProjection() - band=ds.RasterCount - if os.path.exists(file): - os.remove(file) - out_ds:gdal.Dataset=driver.Create(file, xsize, ysize, band, gdal.GDT_Byte) - out_ds.SetGeoTransform(geo) - out_ds.SetProjection(orj) - for b in range(1,band+1): - out_ds.GetRasterBand(b).WriteArray(ds.ReadAsArray(*srcWin,band_list=[b]),*(0,0)) - del out_ds - -def basic_cd(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - ds1:gdal.Dataset=gdal.Open(pth1) - ds2:gdal.Dataset=gdal.Open(pth2) - - cell_size = layer_parent.cell_size - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - - band = ds1.RasterCount - yblocks = ysize // cell_size[1] - - driver = gdal.GetDriverByName('GTiff') - out_tif = os.path.join(Project().other_path, 'temp.tif') - out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - out_ds.SetGeoTransform(geo) - out_ds.SetProjection(proj) - - max_diff = 0 - min_diff = math.inf - - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - for j in range(yblocks + 1):#该改这里了 - - send_message.emit(f'计算{j}/{yblocks}') - block_xy1 = (start1x, start1y+j * cell_size[1]) - block_xy2 = (start2x,start2y+j*cell_size[1]) - block_xy=(0,j * cell_size[1]) - if block_xy1[1] > end1y or block_xy2[1] > end2y: - break - block_size=(xsize, cell_size[1]) - block_size1 = (xsize, cell_size[1]) - block_size2 = (xsize,cell_size[1]) - if block_xy[1] + block_size[1] > ysize: - block_size = (xsize, ysize - block_xy[1]) - if block_xy1[1] + block_size1[1] > end1y: - block_size1 = (xsize,end1y - block_xy1[1]) - if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) - block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) - block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) - - 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}') - del ds2 - del ds1 - 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(geo) - out_normal_ds.SetProjection(proj) - # 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 - send_message.emit('差分法计算完成') - return out_normal_tif - -def LSTS(pth1:str,pth2:str,layer_parent:PairLayer,send_message,n=5,w_size=(3,3)): - ds1:gdal.Dataset=gdal.Open(pth1) - ds2:gdal.Dataset=gdal.Open(pth2) - - cell_size = layer_parent.cell_size - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - - band = ds1.RasterCount - yblocks = ysize // cell_size[1] - - driver = gdal.GetDriverByName('GTiff') - out_tif = os.path.join(Project().other_path, 'temp.tif') - out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - out_ds.SetGeoTransform(geo) - out_ds.SetProjection(proj) - pixnum=w_size[0]*w_size[1] - print('pixnum',pixnum) - max_diff = 0 - min_diff = math.inf - win_h=w_size[0]//2 #half hight of window - win_w=w_size[1]//2 #half width of window - a=[[(i+1)**j for j in range(n+1)] for i in range(pixnum)] - A=np.array(a).astype(np.float64)# - - k_=np.array(range(1,n+1)) - df1=np.zeros(pixnum).astype(np.float64) - df2=np.zeros(pixnum).astype(np.float64) - - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - for j in range(yblocks + 1): - - send_message.emit(f'计算{j}/{yblocks}') - block_xy1 = (start1x, start1y+j * cell_size[1]) - block_xy2 = (start2x,start2y+j*cell_size[1]) - block_xy=(0,j * cell_size[1]) - if block_xy1[1] > end1y or block_xy2[1] > end2y: - break - block_size=(xsize, cell_size[1]) - block_size1 = (xsize, cell_size[1]) - block_size2 = (xsize,cell_size[1]) - if block_xy[1] + block_size[1] > ysize: - block_size = (xsize, ysize - block_xy[1]) - if block_xy1[1] + block_size1[1] > end1y: - block_size1 = (xsize,end1y - block_xy1[1]) - if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) - block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) - block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) - - if band == 1: - block_data1 = block_data1[None, ...] - block_data2 = block_data2[None, ...] - # pdb.set_trace() - else: - block_data1=np.mean(block_data1,0) - block_data2=np.mean(block_data2,0) - block_diff=np.zeros(block_data1.shape).astype(np.float64) - - for i in range(win_h,block_size1[1]-win_h): - for j_ in range(win_w,block_size1[0]-win_w): - pix=0 - - #get b - # b1=block_data[i+win_h:i+win_h] c in range(j_-win_w,j_+win_w+1) - b1=block_data1[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1] - b2=block_data2[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1] - b1=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b1)] - b2=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b2)] - b1=np.expand_dims(np.concatenate(b1,0),1) - b2=np.expand_dims(np.concatenate(b2,0),1) - - #end get b - - x1=np.squeeze(np.linalg.pinv(A).dot(b1)) - x2=np.squeeze(np.linalg.pinv(A).dot(b2)) - #df - k_=range(1,n+1) - for pix in range(1,pixnum+1): - df1[pix-1]=x1[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_])) - df2[pix-1]=x2[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_])) - - #distance 欧式距离 - block_diff[i][j_]=np.dot(df1-df2,df1-df2)**0.5 - - 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}') - del ds2 - del ds1 - 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(geo) - out_normal_ds.SetProjection(proj) - # 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 - - send_message.emit('LSTS法计算完成') - return out_normal_tif - -def CVA(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - ds1:gdal.Dataset=gdal.Open(pth1) - ds2:gdal.Dataset=gdal.Open(pth2) - - cell_size = layer_parent.cell_size - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - - band = ds1.RasterCount - yblocks = ysize // cell_size[1] - - driver = gdal.GetDriverByName('GTiff') - out_tif = os.path.join(Project().other_path, 'temp.tif') - out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - out_ds.SetGeoTransform(geo) - out_ds.SetProjection(proj) - max_diff = 0 - min_diff = math.inf - - - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - - for j in range(yblocks + 1): - - send_message.emit(f'计算{j}/{yblocks}') - block_xy1 = (start1x, start1y+j * cell_size[1]) - block_xy2 = (start2x,start2y+j*cell_size[1]) - block_xy=(0,j * cell_size[1]) - if block_xy1[1] > end1y or block_xy2[1] > end2y: - break - block_size=(xsize, cell_size[1]) - block_size1 = (xsize, cell_size[1]) - block_size2 = (xsize,cell_size[1]) - if block_xy[1] + block_size[1] > ysize: - block_size = (xsize, ysize - block_xy[1]) - if block_xy1[1] + block_size1[1] > end1y: - block_size1 = (xsize,end1y - block_xy1[1]) - if block_xy2[1] + block_size2[1] > end2y: - block_size2 = (xsize, end2y - block_xy2[1]) - block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) - block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) - - if band == 1: - block_data1 = block_data1[None, ...] - block_data2 = block_data2[None, ...] - # pdb.set_trace() - block_diff=np.sum((block_data1-block_data2)**2,0)**0.5 - 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}') - del ds2 - del ds1 - 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(geo) - out_normal_ds.SetProjection(proj) - # 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 - send_message.emit('欧式距离计算完成') - return out_normal_tif - -def acd(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 - send_message.emit('提取重叠区域数据.....') - - ds2:gdal.Dataset=gdal.Open(pth2) - temp_tif2 = os.path.join(Project().other_path,'temp2.tif') - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) - del ds2 - send_message.emit('图像二提取完成') - - - ds1:gdal.Dataset=gdal.Open(pth1) - temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) - del ds1 - send_message.emit('图像一提取完成') - - - - #运算 - send_message.emit('开始ACD计算.....') - time.sleep(0.1) - out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) - ACD(temp_tif1,temp_tif2,out_normal_tif) - #添加投影 - send_message.emit('录入投影信息.....') - time.sleep(0.1) - ds=gdal.Open(out_normal_tif,1) - ds.SetGeoTransform(geo) - ds.SetProjection(proj) - del ds - - return out_normal_tif - -def aht(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 - send_message.emit('提取重叠区域数据.....') - - ds2:gdal.Dataset=gdal.Open(pth2) - temp_tif2 = os.path.join(Project().other_path,'temp2.tif') - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) - del ds2 - send_message.emit('图像二提取完成') - - - ds1:gdal.Dataset=gdal.Open(pth1) - temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) - del ds1 - send_message.emit('图像一提取完成') - - - - #运算 - send_message.emit('开始AHT计算.....') - time.sleep(0.1) - out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) - AHT(temp_tif1,temp_tif2,out_normal_tif) - #添加投影 - send_message.emit('录入投影信息.....') - time.sleep(0.1) - ds=gdal.Open(out_normal_tif,1) - ds.SetGeoTransform(geo) - ds.SetProjection(proj) - del ds - - return out_normal_tif - -def ocd(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 - send_message.emit('提取重叠区域数据.....') - - ds2:gdal.Dataset=gdal.Open(pth2) - temp_tif2 = os.path.join(Project().other_path,'temp2.tif') - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) - del ds2 - send_message.emit('图像二提取完成') - - - ds1:gdal.Dataset=gdal.Open(pth1) - temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) - del ds1 - send_message.emit('图像一提取完成') - - - - #运算 - send_message.emit('开始OCD计算.....') - time.sleep(0.1) - out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) - OCD(temp_tif1,temp_tif2,out_normal_tif,Project().other_path) - #添加投影 - send_message.emit('录入投影信息.....') - time.sleep(0.1) - ds=gdal.Open(out_normal_tif,1) - ds.SetGeoTransform(geo) - ds.SetProjection(proj) - del ds - - return out_normal_tif - -def lhba(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 - send_message.emit('提取重叠区域数据.....') - - ds2:gdal.Dataset=gdal.Open(pth2) - temp_tif2 = os.path.join(Project().other_path,'temp2.tif') - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) - del ds2 - send_message.emit('图像二提取完成') - - - ds1:gdal.Dataset=gdal.Open(pth1) - temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) - del ds1 - send_message.emit('图像一提取完成') - - #运算 - send_message.emit('开始LHBA计算.....') - time.sleep(0.1) - out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) - LHBA(temp_tif1,temp_tif2,out_normal_tif) - #添加投影 - send_message.emit('录入投影信息.....') - time.sleep(0.1) - ds=gdal.Open(out_normal_tif,1) - ds.SetGeoTransform(geo) - ds.SetProjection(proj) - del ds - return out_normal_tif - -def sh(pth1:str,pth2:str,layer_parent:PairLayer,send_message): - xsize = layer_parent.size[0] - ysize = layer_parent.size[1] - geo=layer_parent.grid.geo - proj=layer_parent.grid.proj - #提取公共部分 - send_message.emit('提取重叠区域数据.....') - - ds2:gdal.Dataset=gdal.Open(pth2) - temp_tif2 = os.path.join(Project().other_path,'temp2.tif') - start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) - del ds2 - send_message.emit('图像二提取完成') - - - ds1:gdal.Dataset=gdal.Open(pth1) - temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') - start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) - end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) - warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) - del ds1 - send_message.emit('图像一提取完成') - - #运算 - send_message.emit('开始SH计算.....') - time.sleep(0.1) - out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) - SH(temp_tif1,temp_tif2,out_normal_tif) - #添加投影 - send_message.emit('录入投影信息.....') - time.sleep(0.1) - ds=gdal.Open(out_normal_tif,1) - ds.SetGeoTransform(geo) - ds.SetProjection(proj) - del ds - return out_normal_tif diff --git a/plugins/In_one/scripts/_ACD.pyd b/plugins/In_one/scripts/_ACD.pyd deleted file mode 100644 index e9cbc58..0000000 Binary files a/plugins/In_one/scripts/_ACD.pyd and /dev/null differ diff --git a/plugins/In_one/scripts/_AHT.pyd b/plugins/In_one/scripts/_AHT.pyd deleted file mode 100644 index 5c78396..0000000 Binary files a/plugins/In_one/scripts/_AHT.pyd and /dev/null differ diff --git a/plugins/In_one/scripts/_LHBA.pyd b/plugins/In_one/scripts/_LHBA.pyd deleted file mode 100644 index 3502e2d..0000000 Binary files a/plugins/In_one/scripts/_LHBA.pyd and /dev/null differ diff --git a/plugins/In_one/scripts/_OCD.pyd b/plugins/In_one/scripts/_OCD.pyd deleted file mode 100644 index 00d560a..0000000 Binary files a/plugins/In_one/scripts/_OCD.pyd and /dev/null differ diff --git a/plugins/In_one/scripts/_SH.pyd b/plugins/In_one/scripts/_SH.pyd deleted file mode 100644 index 1384421..0000000 Binary files a/plugins/In_one/scripts/_SH.pyd and /dev/null differ diff --git a/plugins/In_one/scripts/opencv_ffmpeg3415_64.dll b/plugins/In_one/scripts/opencv_ffmpeg3415_64.dll deleted file mode 100644 index 4b73baa..0000000 Binary files a/plugins/In_one/scripts/opencv_ffmpeg3415_64.dll and /dev/null differ diff --git a/plugins/In_one/scripts/opencv_world3415.dll b/plugins/In_one/scripts/opencv_world3415.dll deleted file mode 100644 index fa03239..0000000 Binary files a/plugins/In_one/scripts/opencv_world3415.dll and /dev/null differ diff --git a/plugins/In_one/test.py b/plugins/In_one/test.py deleted file mode 100644 index 8727dd6..0000000 --- a/plugins/In_one/test.py +++ /dev/null @@ -1,3 +0,0 @@ -from In_one.scripts.USCD import ACD -from osgeo import gdal - diff --git a/plugins/filter_collection/__init__.py b/plugins/filter_collection/__init__.py new file mode 100644 index 0000000..4b95114 --- /dev/null +++ b/plugins/filter_collection/__init__.py @@ -0,0 +1,5 @@ +from misc import Register + +FILTER = Register('filter') + +from filter_collection.main import * \ No newline at end of file diff --git a/plugins/some_filter/main.py b/plugins/filter_collection/main.py similarity index 61% rename from plugins/some_filter/main.py rename to plugins/filter_collection/main.py index a6523cd..d5f5eef 100644 --- a/plugins/some_filter/main.py +++ b/plugins/filter_collection/main.py @@ -1,3 +1,4 @@ +from datetime import datetime import os from threading import Thread from PyQt5.QtWidgets import QDialog, QAction @@ -10,19 +11,98 @@ from rscder.plugins.basic import BasicPlugin from rscder.gui.layercombox import RasterLayerCombox from osgeo import gdal, gdal_array from skimage.filters import rank -from skimage.morphology import disk, rectangle +from skimage.morphology import rectangle +from filter_collection import FILTER +from misc import AlgFrontend + +@FILTER.register +class MainFilter(AlgFrontend): + + @staticmethod + def get_name(): + return '均值滤波' + + @staticmethod + def get_widget(parent=None): + widget = QtWidgets.QWidget(parent) + x_size_input = QtWidgets.QLineEdit(widget) + x_size_input.setText('3') + x_size_input.setValidator(QtGui.QIntValidator()) + x_size_input.setObjectName('xinput') + y_size_input = QtWidgets.QLineEdit(widget) + y_size_input.setValidator(QtGui.QIntValidator()) + y_size_input.setObjectName('yinput') + y_size_input.setText('3') + + size_label = QtWidgets.QLabel(widget) + size_label.setText('窗口大小:') + + time_label = QtWidgets.QLabel(widget) + time_label.setText('X') + + hlayout1 = QtWidgets.QHBoxLayout() + + hlayout1.addWidget(size_label) + hlayout1.addWidget(x_size_input) + hlayout1.addWidget(time_label) + hlayout1.addWidget(y_size_input) + + widget.setLayout(hlayout1) + + return widget + + @staticmethod + def get_params(widget:QtWidgets.QWidget=None): + if widget is None: + return dict(x_size=3, y_size=3) + + x_input = widget.findChild(QtWidgets.QLineEdit, 'xinput') + y_input = widget.findChild(QtWidgets.QLineEdit, 'yinput') + + if x_input is None or y_input is None: + return dict(x_size=3, y_size=3) + + x_size = int(x_input.text()) + y_size = int(y_input.text()) + + return dict(x_size=x_size, y_size=y_size) + + @staticmethod + def run_alg(pth, x_size, y_size, *args, **kargs): + 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(int(datetime.now().timestamp() * 1000))) + 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 + class FilterSetting(QDialog): def __init__(self, parent=None): super(FilterSetting, self).__init__(parent) self.setWindowTitle('滤波设置') - # self.setWindowFlags(Qt.WindowStaysOnTopHint) - # self.setFixedSize(300, 200) - # self.setStyleSheet("QDialog{background-color:rgb(255,255,255);}") self.setWindowIcon(IconInstance().FILTER) - # self.setWindowIconText('Filter Setting') - # self.setWindowModality(Qt.ApplicationModal) self.initUI() - # self.show() def initUI(self): self.layer_combox = RasterLayerCombox(self) @@ -78,10 +158,10 @@ class MainPlugin(BasicPlugin): @staticmethod def info(): return { - 'name': 'mean_filter', + 'name': 'FilterCollection', 'author': 'rscder', 'version': '0.0.1', - 'description': 'Mean Filter' + 'description': 'Filter Collections' } def set_action(self): diff --git a/plugins/misc/__init__.py b/plugins/misc/__init__.py new file mode 100644 index 0000000..3a85146 --- /dev/null +++ b/plugins/misc/__init__.py @@ -0,0 +1,3 @@ +from misc.main import * +from misc.utils import * +from misc.table_layer import * \ No newline at end of file diff --git a/plugins/misc/main.py b/plugins/misc/main.py new file mode 100644 index 0000000..6990ebc --- /dev/null +++ b/plugins/misc/main.py @@ -0,0 +1,80 @@ +from threading import Thread +from plugins.misc.utils import Register +from rscder.plugins.basic import BasicPlugin +from PyQt5.QtWidgets import QWidget, QComboBox, QVBoxLayout + + +class MISCPlugin(BasicPlugin): + + @staticmethod + def info(): + return { + 'name': 'MISC Plugin', + 'description': '基本开发工具', + 'author': 'Wangtong', + 'version': '1.0.0', + } + + def set_action(self): + pass + + +class AlgFrontend(object): + + @staticmethod + def get_name(): + return None + + @staticmethod + def get_widget(parent=None): + return QWidget(parent) + + @staticmethod + def get_params(widget=None): + return dict() + + @staticmethod + def run_alg(*args, **kargs): + pass + + +class AlgSelectWidget(QWidget): + + def __init__(self, parent= None, registery:Register = None) -> None: + super().__init__(parent) + self.reg = registery + if registery is None: + return + + self.selectbox = QComboBox(self) + self.selectbox.addItem('--', None) + for key in self.reg.keys(): + alg:AlgFrontend = self.reg[key] + if alg.get_name() is None: + name = key + else: + name = alg.get_name() + self.selectbox.addItem(name, key) + + self.vbox = QVBoxLayout(self) + self.vbox.addWidget(self.selectbox) + + self.params_widget = QWidget(self) + self.vbox.addWidget(self.params_widget) + + self.selectbox.currentIndexChanged.connect(self.on_changed) + + def on_changed(self, index): + if index == 0: + return + + self.alg:AlgFrontend = self.reg[self.selectbox.itemData(index)] + self.vbox.removeWidget(self.params_widget) + self.params_widget = self.alg.get_widget(self) + self.vbox.addWidget(self.params_widget) + + def get_alg_and_params(self): + if self.selectbox.currentIndex() == 0: + return None, None + + return self.alg, self.alg.get_params(self.params_widget) \ No newline at end of file diff --git a/plugins/misc/table_layer.py b/plugins/misc/table_layer.py new file mode 100644 index 0000000..f3d0f70 --- /dev/null +++ b/plugins/misc/table_layer.py @@ -0,0 +1,46 @@ +from datetime import datetime +import math +from random import random +from osgeo import gdal +from rscder.utils.project import Project +import os +from rscder.utils.project import BasicLayer, ResultPointLayer + +def table_layer(pth:str,layer:BasicLayer, name, send_message = None): + if send_message is not None: + 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_{ int(datetime.now().timestamp() * 1000) }.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] + + 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},{int(block_data_xy.mean() > 0.5)}\n') + + result_layer = ResultPointLayer(out_csv, enable=True, proj=layer.proj, geo=layer.geo,result_path={}) + # print(result_layer.result_path) + layer.layer_parent.add_result_layer(result_layer) + if send_message is not None: + send_message.emit('计算完成') \ No newline at end of file diff --git a/plugins/misc/utils.py b/plugins/misc/utils.py new file mode 100644 index 0000000..8b592e4 --- /dev/null +++ b/plugins/misc/utils.py @@ -0,0 +1,42 @@ +import os +import random +import logging +# def generate_temp_file + +class Register: + + def __init__(self, registry_name): + self._dict = {} + self._name = registry_name + + def __setitem__(self, key, value): + if not callable(value): + raise Exception(f"Value of a Registry must be a callable!\nValue: {value}") + if key is None: + key = value.__name__ + if key in self._dict: + logging.warning("Key %s already in registry %s." % (key, self._name)) + self._dict[key] = value + + def register(self, target): + """Decorator to register a function or class.""" + + def add(key, value): + self[key] = value + return value + + if callable(target): + # @reg.register + return add(None, target) + # @reg.register('alias') + return lambda x: add(target, x) + + def __getitem__(self, key): + return self._dict[key] + + def __contains__(self, key): + return key in self._dict + + def keys(self): + """key""" + return self._dict.keys() \ No newline at end of file diff --git a/plugins/plugins.yaml b/plugins/plugins.yaml index 43de036..169900e 100644 --- a/plugins/plugins.yaml +++ b/plugins/plugins.yaml @@ -13,11 +13,11 @@ path: ./plugin\export_to version: 1.0.0 - author: RSCDER - description: MeanFilter + description: FilterCollection enabled: true - module: some_filter - name: MeanFilter - path: ./plugin\some_filter + module: filter_collection + name: FilterCollection + path: ./plugin\filter_collection version: 1.0.0 - author: RSCDER description: Evaluation @@ -27,11 +27,11 @@ path: ./plugin\evaluation version: 1.0.0 - author: RSCDER - description: AllinOne + description: UnsupervisedPlugin enabled: true - module: In_one - name: basic_diff_AllinOne - path: ./plugin\In_one + module: unsupervised_method + name: UnsupervisedPlugin + path: ./plugin\unsupervsied_method version: 1.0.0 - author: RSCDER description: set Change Rate diff --git a/plugins/some_filter/__init__.py b/plugins/some_filter/__init__.py deleted file mode 100644 index 396c649..0000000 --- a/plugins/some_filter/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from some_filter.main import * \ No newline at end of file diff --git a/plugins/thres/__init__.py b/plugins/thres/__init__.py new file mode 100644 index 0000000..8ab0e4c --- /dev/null +++ b/plugins/thres/__init__.py @@ -0,0 +1,5 @@ +from misc import Register + +THRES = Register('Thres') + +from thres.main import * \ No newline at end of file diff --git a/plugins/thres/main.py b/plugins/thres/main.py new file mode 100644 index 0000000..9ae8537 --- /dev/null +++ b/plugins/thres/main.py @@ -0,0 +1,178 @@ +from osgeo import gdal, gdalconst +from rscder.utils.project import Project +import os +from misc import AlgFrontend +from thres import THRES +from PyQt5.QtWidgets import QWidget, QLabel, QHBoxLayout, QLineEdit +from PyQt5.QtGui import QIntValidator, QDoubleValidator + +@THRES.register +class ManulGapAlg(AlgFrontend): + + @staticmethod + def get_name(): + return '手动阈值' + + @staticmethod + def get_params(widget:QWidget=None): + if widget is None: + return dict(gap=125) + lineedit:QLineEdit = widget.layout().findChild(QLineEdit, 'lineedit') + if lineedit is None: + return dict(gap=125) + + gap = int(lineedit.text()) + + return dict(gap=gap) + @staticmethod + def get_widget(parent=None): + + widget = QWidget(parent) + lineedit = QLineEdit(widget) + lineedit.setObjectName('lineedit') + lineedit.setValidator(QIntValidator(1, 254)) + # lineedit. + label = QLabel('阈值：') + + layout = QHBoxLayout() + layout.addWidget(label) + layout.addWidget(lineedit) + + widget.setLayout(layout) + return widget + + + @staticmethod + def run_alg(pth,name = '',gap = 125, send_message = None): + + ds = gdal.Open(pth) + band = ds.GetRasterBand(1) + + 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 + + if send_message is not None: + 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 + if send_message is not None: + send_message.emit('自定义阈值分割完成') + return out_th + + +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 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 + # ManulGapAlg.run_alg \ No newline at end of file diff --git a/plugins/unsupervised_method/__init__.py b/plugins/unsupervised_method/__init__.py new file mode 100644 index 0000000..c052402 --- /dev/null +++ b/plugins/unsupervised_method/__init__.py @@ -0,0 +1 @@ +from unsupervised_method.main import * \ No newline at end of file diff --git a/plugins/unsupervised_method/main.py b/plugins/unsupervised_method/main.py new file mode 100644 index 0000000..a398a2e --- /dev/null +++ b/plugins/unsupervised_method/main.py @@ -0,0 +1,202 @@ +from threading import Thread +from plugins.misc.main import AlgFrontend +from rscder.gui.actions import ActionManager +from rscder.plugins.basic import BasicPlugin +from PyQt5.QtWidgets import QAction, QToolBar, QMenu, QDialog, QHBoxLayout, QVBoxLayout, QPushButton,QWidget,QLabel,QLineEdit,QPushButton,QComboBox,QDialogButtonBox + +from rscder.gui.layercombox import PairLayerCombox +from rscder.utils.icons import IconInstance +from filter_collection import FILTER +from .scripts import UNSUPER_CD +from thres import THRES +from misc import table_layer, AlgSelectWidget + +class UnsupervisedCDMethod(QDialog): + def __init__(self,parent=None, alg:AlgFrontend=None): + super(UnsupervisedCDMethod, self).__init__(parent) + self.alg = alg + self.setWindowTitle('无监督变化检测:{}'.format(alg.get_name())) + self.setWindowIcon(IconInstance().LOGO) + self.initUI() + self.setMinimumWidth(500) + + def initUI(self): + #图层 + self.layer_combox = PairLayerCombox(self) + layerbox = QHBoxLayout() + layerbox.addWidget(self.layer_combox) + + self.filter_select = AlgSelectWidget(self, FILTER) + self.param_widget = self.alg.get_widget(self) + self.unsupervised_menu = self.param_widget + self.thres_select = AlgSelectWidget(self, THRES) + + self.ok_button = QPushButton('确定', self) + self.ok_button.setIcon(IconInstance().OK) + self.ok_button.clicked.connect(self.accept) + self.ok_button.setDefault(True) + + self.cancel_button = QPushButton('取消', self) + self.cancel_button.setIcon(IconInstance().CANCEL) + self.cancel_button.clicked.connect(self.reject) + self.cancel_button.setDefault(False) + buttonbox=QDialogButtonBox(self) + buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole) + buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole) + buttonbox.setCenterButtons(True) + + totalvlayout=QVBoxLayout() + totalvlayout.addLayout(layerbox) + totalvlayout.addWidget(self.filter_select) + if self.param_widget is not None: + totalvlayout.addWidget(self.param_widget) + totalvlayout.addWidget(self.thres_select) + totalvlayout.addWidget(buttonbox) + totalvlayout.addStretch() + + self.setLayout(totalvlayout) + +class UnsupervisedCD(QDialog): + def __init__(self,parent=None): + super(UnsupervisedCD, self).__init__(parent) + self.setWindowTitle('无监督变化检测') + self.setWindowIcon(IconInstance().LOGO) + self.setMinimumWidth(600) + self.initUI() + + def initUI(self): + #图层 + self.layer_combox = PairLayerCombox(self) + layerbox = QHBoxLayout() + layerbox.addWidget(self.layer_combox) + + self.filter_select = AlgSelectWidget(self, FILTER) + self.unsupervised_select = AlgSelectWidget(self, UNSUPER_CD) + self.thres_select = AlgSelectWidget(self, THRES) + + self.ok_button = QPushButton('确定', self) + self.ok_button.setIcon(IconInstance().OK) + self.ok_button.clicked.connect(self.accept) + self.ok_button.setDefault(True) + + self.cancel_button = QPushButton('取消', self) + self.cancel_button.setIcon(IconInstance().CANCEL) + self.cancel_button.clicked.connect(self.reject) + self.cancel_button.setDefault(False) + buttonbox=QDialogButtonBox(self) + buttonbox.addButton(self.ok_button,QDialogButtonBox.NoRole) + buttonbox.addButton(self.cancel_button,QDialogButtonBox.NoRole) + buttonbox.setCenterButtons(True) + + totalvlayout=QVBoxLayout() + totalvlayout.addLayout(layerbox) + totalvlayout.addWidget(self.filter_select) + totalvlayout.addWidget(self.unsupervised_select) + totalvlayout.addWidget(self.thres_select) + totalvlayout.addWidget(buttonbox) + totalvlayout.addStretch() + + self.setLayout(totalvlayout) + + + +class UnsupervisedPlugin(BasicPlugin): + + + @staticmethod + def info(): + return { + 'name': 'UnsupervisedPlugin', + 'description': 'UnsupervisedPlugin', + 'author': 'RSCDER', + 'version': '1.0.0', + } + + def set_action(self): + unsupervised_menu = QMenu('&无监督变化检测', self.mainwindow) + unsupervised_menu.setIcon(IconInstance().UNSUPERVISED) + ActionManager().change_detection_menu.addMenu(unsupervised_menu) + basic_diff_method_in_one = QAction(IconInstance().UNSUPERVISED, '&无监督变化检测', self.mainwindow) + basic_diff_method_in_one.triggered.connect(self.run) + toolbar:QToolBar = ActionManager().add_toolbar('无监督变化检测') + toolbar.addAction(basic_diff_method_in_one) + # self.mainwindow.addToolbar(toolbar) + # print(UNSUPER_CD.keys()) + for key in UNSUPER_CD.keys(): + alg:AlgFrontend = UNSUPER_CD[key] + if alg.get_name() is None: + name = key + else: + name = alg.get_name() + + action = QAction(name, unsupervised_menu) + action.triggered.connect(lambda : self.run_cd(alg)) + + unsupervised_menu.addAction(action) + + def run(self): + myDialog= UnsupervisedCD(self.mainwindow) + myDialog.show() + if myDialog.exec_()==QDialog.Accepted: + w=myDialog + t=Thread(target=self.run_alg,args=(w,)) + t.start() + + def run_cd(self, alg): + dialog = UnsupervisedCDMethod(self.mainwindow, alg) + dialog.show() + + if dialog.exec_() == QDialog.Accepted: + t = Thread(target=self.run_cd_alg, args=(dialog,)) + t.start() + + def run_cd_alg(self, w:UnsupervisedCDMethod): + + layer1=w.layer_combox.layer1 + pth1 = w.layer_combox.layer1.path + pth2 = w.layer_combox.layer2.path + name = layer1.layer_parent.name + + falg, fparams = w.filter_select.get_alg_and_params() + cdalg = w.alg + cdparams = w.alg.get_params() + thalg, thparams = w.thres_select.get_alg_and_params() + + if cdalg is None or thalg is None: + return + + if falg is not None: + pth1 = falg.run_alg(pth1, name=name, send_message=self.send_message, **fparams) + pth2 = falg.run_alg(pth2, name=name, send_message=self.send_message, **fparams) + + + cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, send_message=self.send_message,**cdparams) + thpth = thalg.run_alg(cdpth, name=name, send_message=self.send_message, **thparams) + + table_layer(thpth,layer1,name,self.send_message) + + def run_alg(self,w:UnsupervisedCD): + + layer1=w.layer_combox.layer1 + pth1 = w.layer_combox.layer1.path + pth2 = w.layer_combox.layer2.path + name = layer1.layer_parent.name + + falg, fparams = w.filter_select.get_alg_and_params() + cdalg, cdparams = w.unsupervised_select.get_alg_and_params() + thalg, thparams = w.thres_select.get_alg_and_params() + + if cdalg is None or thalg is None: + return + + if falg is not None: + pth1 = falg.run_alg(pth1, name=name, send_message=self.send_message, **fparams) + pth2 = falg.run_alg(pth2, name=name, send_message=self.send_message, **fparams) + + + cdpth = cdalg.run_alg(pth1, pth2, layer1.layer_parent, send_message=self.send_message,**cdparams) + + + thpth = thalg.run_alg(cdpth, name=name, send_message=self.send_message, **thparams) + + table_layer(thpth,layer1,name,self.send_message) \ No newline at end of file diff --git a/plugins/In_one/pic.py b/plugins/unsupervised_method/pic.py similarity index 100% rename from plugins/In_one/pic.py rename to plugins/unsupervised_method/pic.py diff --git a/plugins/unsupervised_method/readme.md b/plugins/unsupervised_method/readme.md new file mode 100644 index 0000000..e69de29 diff --git a/plugins/In_one/scripts/ACD.py b/plugins/unsupervised_method/scripts/ACD.py similarity index 100% rename from plugins/In_one/scripts/ACD.py rename to plugins/unsupervised_method/scripts/ACD.py diff --git a/plugins/In_one/scripts/AHT.py b/plugins/unsupervised_method/scripts/AHT.py similarity index 100% rename from plugins/In_one/scripts/AHT.py rename to plugins/unsupervised_method/scripts/AHT.py diff --git a/plugins/In_one/scripts/LHBA.py b/plugins/unsupervised_method/scripts/LHBA.py similarity index 100% rename from plugins/In_one/scripts/LHBA.py rename to plugins/unsupervised_method/scripts/LHBA.py diff --git a/plugins/In_one/scripts/OCD.py b/plugins/unsupervised_method/scripts/OCD.py similarity index 100% rename from plugins/In_one/scripts/OCD.py rename to plugins/unsupervised_method/scripts/OCD.py diff --git a/plugins/In_one/scripts/SH.py b/plugins/unsupervised_method/scripts/SH.py similarity index 100% rename from plugins/In_one/scripts/SH.py rename to plugins/unsupervised_method/scripts/SH.py diff --git a/plugins/unsupervised_method/scripts/__init__.py b/plugins/unsupervised_method/scripts/__init__.py new file mode 100644 index 0000000..f85a810 --- /dev/null +++ b/plugins/unsupervised_method/scripts/__init__.py @@ -0,0 +1,688 @@ +from datetime import datetime +from osgeo import gdal +import math,os +import time +from PyQt5 import QtWidgets +from sklearn.cluster import k_means +from rscder.utils.geomath import geo2imageRC, imageRC2geo +from rscder.utils.project import Project, PairLayer +from misc import Register, AlgFrontend + +UNSUPER_CD = Register('unsuper cd') + +import numpy as np +from .ACD import ACD +from .AHT import AHT +from .OCD import OCD +from .LHBA import LHBA +from .SH import SH + +def warp(file,ds:gdal.Dataset,srcWin=[0,0,0,0]): + driver = gdal.GetDriverByName('GTiff') + xsize=ds.RasterXSize + ysize=ds.RasterYSize + geo=ds.GetGeoTransform() + orj=ds.GetProjection() + band=ds.RasterCount + if os.path.exists(file): + os.remove(file) + out_ds:gdal.Dataset=driver.Create(file, xsize, ysize, band, gdal.GDT_Byte) + out_ds.SetGeoTransform(geo) + out_ds.SetProjection(orj) + for b in range(1,band+1): + out_ds.GetRasterBand(b).WriteArray(ds.ReadAsArray(*srcWin,band_list=[b]),*(0,0)) + del out_ds + +@UNSUPER_CD.register +class BasicCD(AlgFrontend): + + @staticmethod + def get_name(): + return '差分法' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None,*args, **kargs): + + ds1:gdal.Dataset=gdal.Open(pth1) + ds2:gdal.Dataset=gdal.Open(pth2) + + cell_size = layer_parent.cell_size + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + + band = ds1.RasterCount + yblocks = ysize // cell_size[1] + + driver = gdal.GetDriverByName('GTiff') + out_tif = os.path.join(Project().other_path, 'temp.tif') + out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + out_ds.SetGeoTransform(geo) + out_ds.SetProjection(proj) + + max_diff = 0 + min_diff = math.inf + + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + for j in range(yblocks + 1):#该改这里了 + if send_message is not None: + send_message.emit(f'计算{j}/{yblocks}') + block_xy1 = (start1x, start1y+j * cell_size[1]) + block_xy2 = (start2x,start2y+j*cell_size[1]) + block_xy=(0,j * cell_size[1]) + if block_xy1[1] > end1y or block_xy2[1] > end2y: + break + block_size=(xsize, cell_size[1]) + block_size1 = (xsize, cell_size[1]) + block_size2 = (xsize,cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + if block_xy1[1] + block_size1[1] > end1y: + block_size1 = (xsize,end1y - block_xy1[1]) + if block_xy2[1] + block_size2[1] > end2y: + block_size2 = (xsize, end2y - block_xy2[1]) + block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) + block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) + + 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) + if send_message is not None: + + send_message.emit(f'完成{j}/{yblocks}') + del ds2 + del ds1 + out_ds.FlushCache() + del out_ds + if send_message is not None: + 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(geo) + out_normal_ds.SetProjection(proj) + # 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 + if send_message is not None: + send_message.emit('差分法计算完成') + return out_normal_tif + +@UNSUPER_CD.register +class LSTS(AlgFrontend): + + @staticmethod + def get_name(): + return 'LSTS' + + @staticmethod + def get_widget(parent=None): + + widget = QtWidgets.QWidget(parent) + + + def get_params(self): + return dict(n=5, w_size=(3,3)) + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message=None,n=5,w_size=(3,3), *args, **kws): + ds1:gdal.Dataset=gdal.Open(pth1) + ds2:gdal.Dataset=gdal.Open(pth2) + + cell_size = layer_parent.cell_size + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + + band = ds1.RasterCount + yblocks = ysize // cell_size[1] + + driver = gdal.GetDriverByName('GTiff') + out_tif = os.path.join(Project().other_path, '%d.tif'%(int(datetime.now().timestamp() * 1000))) + out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + out_ds.SetGeoTransform(geo) + out_ds.SetProjection(proj) + pixnum=w_size[0]*w_size[1] + # send_message.emit('pixnum:'pixnum) + max_diff = 0 + min_diff = math.inf + win_h=w_size[0]//2 #half hight of window + win_w=w_size[1]//2 #half width of window + a=[[(i+1)**j for j in range(n+1)] for i in range(pixnum)] + A=np.array(a).astype(np.float64)# + + k_=np.array(range(1,n+1)) + df1=np.zeros(pixnum).astype(np.float64) + df2=np.zeros(pixnum).astype(np.float64) + + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + for j in range(yblocks + 1): + if send_message is not None: + send_message.emit(f'计算{j}/{yblocks}') + block_xy1 = (start1x, start1y+j * cell_size[1]) + block_xy2 = (start2x,start2y+j*cell_size[1]) + block_xy=(0,j * cell_size[1]) + if block_xy1[1] > end1y or block_xy2[1] > end2y: + break + block_size=(xsize, cell_size[1]) + block_size1 = (xsize, cell_size[1]) + block_size2 = (xsize,cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + if block_xy1[1] + block_size1[1] > end1y: + block_size1 = (xsize,end1y - block_xy1[1]) + if block_xy2[1] + block_size2[1] > end2y: + block_size2 = (xsize, end2y - block_xy2[1]) + block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) + block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) + + if band == 1: + block_data1 = block_data1[None, ...] + block_data2 = block_data2[None, ...] + # pdb.set_trace() + else: + block_data1=np.mean(block_data1,0) + block_data2=np.mean(block_data2,0) + block_diff=np.zeros(block_data1.shape).astype(np.float64) + + for i in range(win_h,block_size1[1]-win_h): + for j_ in range(win_w,block_size1[0]-win_w): + pix=0 + + #get b + # b1=block_data[i+win_h:i+win_h] c in range(j_-win_w,j_+win_w+1) + b1=block_data1[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1] + b2=block_data2[i-win_h:i+win_h+1,j_-win_w:j_+win_w+1] + b1=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b1)] + b2=[b if (r+1)//2 else b[::-1] for r,b in enumerate(b2)] + b1=np.expand_dims(np.concatenate(b1,0),1) + b2=np.expand_dims(np.concatenate(b2,0),1) + + x1=np.squeeze(np.linalg.pinv(A).dot(b1)) + x2=np.squeeze(np.linalg.pinv(A).dot(b2)) + #df + k_=range(1,n+1) + for pix in range(1,pixnum+1): + df1[pix-1]=x1[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_])) + df2[pix-1]=x2[1:n+1].dot(np.array([k*(pix**(k-1)) for k in k_])) + + #distance 欧式距离 + block_diff[i][j_]=np.dot(df1-df2,df1-df2)**0.5 + + 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}') + del ds2 + del ds1 + out_ds.FlushCache() + del out_ds + if send_message is not None: + 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(geo) + out_normal_ds.SetProjection(proj) + # 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 + + if send_message is not None: + send_message.emit('LSTS法计算完成') + return out_normal_tif + + +@UNSUPER_CD.register +class CVAAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'CVA' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + ds1:gdal.Dataset=gdal.Open(pth1) + ds2:gdal.Dataset=gdal.Open(pth2) + + cell_size = layer_parent.cell_size + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + + band = ds1.RasterCount + yblocks = ysize // cell_size[1] + + driver = gdal.GetDriverByName('GTiff') + out_tif = os.path.join(Project().other_path, 'temp.tif') + out_ds = driver.Create(out_tif, xsize, ysize, 1, gdal.GDT_Float32) + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + out_ds.SetGeoTransform(geo) + out_ds.SetProjection(proj) + max_diff = 0 + min_diff = math.inf + + + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + + for j in range(yblocks + 1): + if send_message is not None: + send_message.emit(f'计算{j}/{yblocks}') + block_xy1 = (start1x, start1y+j * cell_size[1]) + block_xy2 = (start2x,start2y+j*cell_size[1]) + block_xy=(0,j * cell_size[1]) + if block_xy1[1] > end1y or block_xy2[1] > end2y: + break + block_size=(xsize, cell_size[1]) + block_size1 = (xsize, cell_size[1]) + block_size2 = (xsize,cell_size[1]) + if block_xy[1] + block_size[1] > ysize: + block_size = (xsize, ysize - block_xy[1]) + if block_xy1[1] + block_size1[1] > end1y: + block_size1 = (xsize,end1y - block_xy1[1]) + if block_xy2[1] + block_size2[1] > end2y: + block_size2 = (xsize, end2y - block_xy2[1]) + block_data1 = ds1.ReadAsArray(*block_xy1, *block_size1) + block_data2 = ds2.ReadAsArray(*block_xy2, *block_size2) + + if band == 1: + block_data1 = block_data1[None, ...] + block_data2 = block_data2[None, ...] + # pdb.set_trace() + block_diff=np.sum((block_data1-block_data2)**2,0)**0.5 + 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) + if send_message is not None: + send_message.emit(f'完成{j}/{yblocks}') + del ds2 + del ds1 + out_ds.FlushCache() + del out_ds + if send_message is not None: + 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(geo) + out_normal_ds.SetProjection(proj) + # 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 + if send_message is not None: + send_message.emit('欧式距离计算完成') + return out_normal_tif + +@UNSUPER_CD.register +class ACDAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'ACD' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + if send_message is None: + class Empty: + + def emit(self, *args, **kws): + print(args) + send_message = Empty() + # send_message.emit = print + + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + #提取公共部分 + send_message.emit('提取重叠区域数据.....') + + ds2:gdal.Dataset=gdal.Open(pth2) + temp_tif2 = os.path.join(Project().other_path,'temp2.tif') + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) + del ds2 + send_message.emit('图像二提取完成') + + + ds1:gdal.Dataset=gdal.Open(pth1) + temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) + del ds1 + send_message.emit('图像一提取完成') + + + + #运算 + send_message.emit('开始ACD计算.....') + time.sleep(0.1) + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + ACD(temp_tif1,temp_tif2,out_normal_tif) + #添加投影 + send_message.emit('录入投影信息.....') + time.sleep(0.1) + ds=gdal.Open(out_normal_tif,1) + ds.SetGeoTransform(geo) + ds.SetProjection(proj) + del ds + + return out_normal_tif + + +@UNSUPER_CD.register +class AHTAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'AHT' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + if send_message is None: + class Empty: + + def emit(self, *args, **kws): + print(args) + send_message = Empty() + + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + #提取公共部分 + send_message.emit('提取重叠区域数据.....') + + ds2:gdal.Dataset=gdal.Open(pth2) + temp_tif2 = os.path.join(Project().other_path,'temp2.tif') + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) + del ds2 + send_message.emit('图像二提取完成') + + + ds1:gdal.Dataset=gdal.Open(pth1) + temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) + del ds1 + send_message.emit('图像一提取完成') + + + + #运算 + send_message.emit('开始AHT计算.....') + time.sleep(0.1) + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + AHT(temp_tif1,temp_tif2,out_normal_tif) + #添加投影 + send_message.emit('录入投影信息.....') + time.sleep(0.1) + ds=gdal.Open(out_normal_tif,1) + ds.SetGeoTransform(geo) + ds.SetProjection(proj) + del ds + + return out_normal_tif + + +@UNSUPER_CD.register +class OCDAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'OCD' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + if send_message is None: + class Empty: + + def emit(self, *args, **kws): + print(args) + send_message = Empty() + + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + #提取公共部分 + send_message.emit('提取重叠区域数据.....') + + ds2:gdal.Dataset=gdal.Open(pth2) + temp_tif2 = os.path.join(Project().other_path,'temp2.tif') + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) + del ds2 + send_message.emit('图像二提取完成') + + + ds1:gdal.Dataset=gdal.Open(pth1) + temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) + del ds1 + send_message.emit('图像一提取完成') + + + + #运算 + send_message.emit('开始OCD计算.....') + time.sleep(0.1) + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + OCD(temp_tif1,temp_tif2,out_normal_tif,Project().other_path) + #添加投影 + send_message.emit('录入投影信息.....') + time.sleep(0.1) + ds=gdal.Open(out_normal_tif,1) + ds.SetGeoTransform(geo) + ds.SetProjection(proj) + del ds + + return out_normal_tif + +@UNSUPER_CD.register +class LHBAAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'LHBA' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + if send_message is None: + class Empty: + + def emit(self, *args, **kws): + print(args) + send_message = Empty() + + + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + #提取公共部分 + send_message.emit('提取重叠区域数据.....') + + ds2:gdal.Dataset=gdal.Open(pth2) + temp_tif2 = os.path.join(Project().other_path,'temp2.tif') + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) + del ds2 + send_message.emit('图像二提取完成') + + + ds1:gdal.Dataset=gdal.Open(pth1) + temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) + del ds1 + send_message.emit('图像一提取完成') + + #运算 + send_message.emit('开始LHBA计算.....') + time.sleep(0.1) + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + LHBA(temp_tif1,temp_tif2,out_normal_tif) + #添加投影 + send_message.emit('录入投影信息.....') + time.sleep(0.1) + ds=gdal.Open(out_normal_tif,1) + ds.SetGeoTransform(geo) + ds.SetProjection(proj) + del ds + return out_normal_tif + + +@UNSUPER_CD.register +class SHAlg(AlgFrontend): + + @staticmethod + def get_name(): + return 'SH' + + @staticmethod + def run_alg(pth1:str,pth2:str,layer_parent:PairLayer,send_message = None, *args, **kargs): + + if send_message is None: + class Empty: + + def emit(self, *args, **kws): + print(args) + send_message = Empty() + + + + xsize = layer_parent.size[0] + ysize = layer_parent.size[1] + geo=layer_parent.grid.geo + proj=layer_parent.grid.proj + #提取公共部分 + send_message.emit('提取重叠区域数据.....') + + ds2:gdal.Dataset=gdal.Open(pth2) + temp_tif2 = os.path.join(Project().other_path,'temp2.tif') + start2x,start2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end2x,end2y=geo2imageRC(ds2.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif2,ds2,srcWin=[start2x,start2y,xsize,ysize]) + del ds2 + send_message.emit('图像二提取完成') + + + ds1:gdal.Dataset=gdal.Open(pth1) + temp_tif1 = os.path.join(Project().other_path, 'temp1.tif') + start1x,start1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[0],layer_parent.mask.xy[1]) + end1x,end1y=geo2imageRC(ds1.GetGeoTransform(),layer_parent.mask.xy[2],layer_parent.mask.xy[3]) + warp(temp_tif1,ds1,srcWin=[start1x,start1y,xsize,ysize]) + del ds1 + send_message.emit('图像一提取完成') + + #运算 + send_message.emit('开始SH计算.....') + time.sleep(0.1) + out_normal_tif = os.path.join(Project().cmi_path, '{}_{}_cmi.tif'.format(layer_parent.name, int(np.random.rand() * 100000))) + SH(temp_tif1,temp_tif2,out_normal_tif) + #添加投影 + send_message.emit('录入投影信息.....') + time.sleep(0.1) + ds=gdal.Open(out_normal_tif,1) + ds.SetGeoTransform(geo) + ds.SetProjection(proj) + del ds + return out_normal_tif diff --git a/rscder/gui/actions.py b/rscder/gui/actions.py index 43e2c73..19a584d 100644 --- a/rscder/gui/actions.py +++ b/rscder/gui/actions.py @@ -2,7 +2,8 @@ import logging import os from pathlib import Path from PyQt5 import QtCore, QtGui, QtWidgets -from PyQt5.QtWidgets import QAction, QActionGroup, QLabel, QFileDialog, QMenuBar +from PyQt5.QtCore import Qt, QSize, QSettings, pyqtSignal +from PyQt5.QtWidgets import QAction, QActionGroup, QLabel, QFileDialog, QMenuBar, QToolBar from rscder.gui import project from rscder.gui.project import Create from rscder.utils.icons import IconInstance @@ -32,6 +33,8 @@ class ActionManager(QtCore.QObject): self.action_groups = {} self.action_group_actions = {} + self.toolbars = {} + self.double_map = double_map self.layer_tree = layer_tree self.follow_box = follow_box @@ -74,9 +77,21 @@ class ActionManager(QtCore.QObject): self.help_menu = menubar.addMenu( '&帮助') def set_toolbar(self, toolbar): - self.toolbar = toolbar + self.toolbar:QToolBar = toolbar self.toolbar.setIconSize(QtCore.QSize(24, 24)) - + + def add_toolbar(self, name=None): + + toolbar = self.w_parent.addToolBar(name) + toolbar.setMovable(True) + toolbar.setFloatable(False) + toolbar.setIconSize(QSize(32, 32)) + toolbar.setToolButtonStyle(Qt.ToolButtonTextBesideIcon) + toolbar.setContextMenuPolicy(Qt.PreventContextMenu) + toolbar.setLayoutDirection(Qt.LeftToRight) + + return toolbar + def set_status_bar(self, status_bar): self.status_bar = status_bar @@ -131,6 +146,9 @@ class ActionManager(QtCore.QObject): zomm_out.setChecked(False) zomm_in.setCheckable(True) zomm_in.setChecked(False) + + toolbar = self.add_toolbar('Map') + toolbar.addActions([pan, zomm_out, zomm_in, locate]) self.double_map.connect_map_tool(pan, zomm_in, zomm_out) # self.double_map.connect_grid_show(grid_line) @@ -148,6 +166,9 @@ class ActionManager(QtCore.QObject): plugin_list = self.add_action(QAction(IconInstance().PLUGINS,'&插件列表', self.w_parent), 'Plugin') plugin_list.triggered.connect(self.plugin_list) + # toolbar = self.add_toolbar('Plugin') + # toolbar.addAction(plugin_list) + self.plugin_menu.addAction(plugin_list) self.message_box.info('菜单初始化完成') diff --git a/rscder/plugins/loader.py b/rscder/plugins/loader.py index 87ceffe..5ce875b 100644 --- a/rscder/plugins/loader.py +++ b/rscder/plugins/loader.py @@ -62,7 +62,10 @@ class PluginLoader(QObject): self.plugins.append(obj(self.ctx)) break except Exception as e: + # import traceback + # traceback.print_exc() self.ctx['message_box'].error(f'{plugin["name"]} load error: {e}') + self.plugin_loaded.emit() \ No newline at end of file