# -*- coding: utf-8 -*-
"""
/***************************************************************************
 GLCM
                                 A QGIS plugin
 Grey Level Co-occurence Matrices
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2024-11-11
        git sha              : $Format:%H$
        copyright            : (C) 2024 by Tim Le Bas
        email                : tlb@noc.ac.uk
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""
from qgis.PyQt.QtCore import QSettings, QTranslator, QCoreApplication
from qgis.PyQt.QtGui import QIcon
from qgis.PyQt.QtWidgets import QAction,QFileDialog,QMessageBox 
from qgis.core import QgsProject

# Initialize Qt resources from file resources.py
from .resources import *
# Import the code for the dialog
from .GLCM_dialog import GLCMDialog
import os.path

from qgis.core import Qgis,QgsMessageLog
from qgis.gui import QgsMessageBar

from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterRasterLayer
from qgis.core import QgsProcessingParameterRasterDestination
from qgis.core import QgsProcessingParameterDefinition
from qgis.core import QgsVectorLayer
from qgis.core import (QgsSymbol,QgsSimpleFillSymbolLayer,QgsRendererCategory,QgsCategorizedSymbolRenderer)
from qgis.core import QgsRasterLayer
from qgis.core import QgsRasterBandStats

import processing

class LoadingScreenDlg:
    """Loading screen animation."""
    from qgis.PyQt.QtWidgets import QDialog, QLabel 
    from qgis.PyQt.QtGui import QMovie, QPalette, QColor

    def __init__(self, gif_path):
        self.dlg = self.QDialog()
        self.dlg.setWindowTitle("Please Wait")
        self.dlg.setWindowModality(False)
        self.dlg.setFixedSize(200, 100)
        pal = self.QPalette()
        role = self.QPalette.Background
        pal.setColor(role, self.QColor(255, 255, 255))
        self.dlg.setPalette(pal)
        self.label_animation = self.QLabel(self.dlg)
        self.movie = self.QMovie(gif_path)
        self.label_animation.setMovie(self.movie)

    def start_animation(self):
        self.movie.start()
        self.dlg.show()
        return

    def stop_animation(self):
        self.movie.stop()
        self.dlg.done(0)       

class GLCM:
    """QGIS Plugin Implementation."""
    """
    def __init__(self, iface):
        # Save reference to the QGIS interface
        self.iface = iface
        # initialize plugin directory
        self.plugin_dir = os.path.dirname(__file__)
        # initialize locale
        locale = QSettings().value('locale/userLocale')[0:2]
        locale_path = os.path.join(
            self.plugin_dir,
            'i18n',
            'GLCM_{}.qm'.format(locale))

        if os.path.exists(locale_path):
            self.translator = QTranslator()
            self.translator.load(locale_path)
            QCoreApplication.installTranslator(self.translator)

        # Declare instance attributes
        self.actions = []
        self.menu = self.tr(u'&GLCM')

        # Check if plugin was started the first time in current QGIS session
        # Must be set in initGui() to survive plugin reloads
        self.first_start = None

    # noinspection PyMethodMayBeStatic
    def tr(self, message):

        # noinspection PyTypeChecker,PyArgumentList,PyCallByClass
        return QCoreApplication.translate('GLCM', message)


    def add_action(
        self,
        icon_path,
        text,
        callback,
        enabled_flag=True,
        add_to_menu=True,
        add_to_toolbar=True,
        status_tip=None,
        whats_this=None,
        parent=None):

        icon = QIcon(icon_path)
        action = QAction(icon, text, parent)
        action.triggered.connect(callback)
        action.setEnabled(enabled_flag)

        if status_tip is not None:
            action.setStatusTip(status_tip)

        if whats_this is not None:
            action.setWhatsThis(whats_this)

        if add_to_toolbar:
            # Adds plugin icon to Plugins toolbar
            self.iface.addToolBarIcon(action)

        if add_to_menu:
            self.iface.addPluginToMenu(
                self.menu,
                action)

        self.actions.append(action)

        return action

    def initGui(self):
 
        icon_path = ':/plugins/GLCM/icon.png'
        self.add_action(
            icon_path,
            text=self.tr(u'Creates Grey-Level Co-occurrence Matrices '),
            callback=self.run,
            parent=self.iface.mainWindow())

        # will be set False in run()
        self.first_start = True


    def unload(self):
        for action in self.actions:
            self.iface.removePluginMenu(
                self.tr(u'&GLCM'),
                action)
            self.iface.removeToolBarIcon(action)

    """
    def select_output_file(self): 
        filename, _filter = QFileDialog.getSaveFileName(selfMT.dlg, "Select output image file ","", '*.img') 
        selfMT.dlg.outputFile.setText(filename) 
        if os.path.exists(filename):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
         
    def select_input_file(self): 
        filename, _filter = QFileDialog.getOpenFileName(selfMT.dlg, "Select input image file ","", '*.img *.tif') 
        selfMT.dlg.inputCombo.clear() 
        selfMT.dlg.inputCombo.insertItem(0,filename)
        selfMT.dlg.inputCombo.setCurrentIndex(0)
        selectedLayerIndex = selfMT.dlg.GLCM_type.currentIndex()
        names = ["contrast","dissimilarity","homogeneity","energy","entropy","ASM","std","mean"]
        GLCM_type =names[selectedLayerIndex]
        blocksize = selfMT.dlg.blocksize.text()
        if blocksize == "":
            blocksize = "5"
        sampsize = selfMT.dlg.sampsize.text()
        if sampsize == "":
            sampsize = "1"
        greylevels = selfMT.dlg.levels.text()
        if greylevels == "":
            greylevels = "8"
        testangle = selfMT.dlg.angle.text()
        if testangle == "":
            testangle = "0"
        autoPoly = filename[:-4]+"_"+GLCM_type+"_"+blocksize+"_"+sampsize+"_"+greylevels+"_"+testangle+".tif"
        selfMT.dlg.outputFile.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        
    def indexChanged(self): 
        selectedLayerIndex = selfMT.dlg.inputCombo.currentIndex()
        currentText = selfMT.dlg.inputCombo.currentText()
        layers = QgsProject.instance().mapLayers().values()
        a=0
        filename="NULL"
        for layer in (layer1 for layer1 in layers if str(layer1.type())== "1" or str(layer1.type())== "LayerType.Raster"):
            if a == selectedLayerIndex:
                filename = str(layer.source())
            a=a+1
        filename1= selfMT.dlg.outputFile.text()[0:len(currentText[:-4])]
        if filename1[0:3] == "_co" or currentText[:-4] == filename1[0:len(currentText[:-4])]:
            filename = currentText
        #autofill
        currentGLCMIndex = selfMT.dlg.GLCM_type.currentIndex()
        names = ["contrast","dissimilarity","homogeneity","energy","entropy","ASM","std","mean"]
        GLCM_type =names[currentGLCMIndex]
        blocksize = selfMT.dlg.blocksize.text()
        if blocksize == "":
            blocksize = "5"
        sampsize = selfMT.dlg.sampsize.text()
        if sampsize == "":
            sampsize = "1"
        greylevels = selfMT.dlg.levels.text()
        if greylevels == "":
            greylevels = "8"
        testangle = selfMT.dlg.angle.text()
        if testangle == "":
            testangle = "0"
        autoPoly = filename[:-4]+"_"+GLCM_type+"_"+blocksize+"_"+sampsize+"_"+greylevels+"_"+testangle+".tif"
        selfMT.dlg.outputFile.setText(autoPoly)
        if os.path.exists(autoPoly):
            selfMT.dlg.exists1.setText("Existing file will be overwritten")
        else:
            selfMT.dlg.exists1.setText("")
        
    def help(self): 
        import webbrowser
        import marinetools
        MThelp = os.path.dirname(marinetools.__file__) + "\\GLCM\\GLCM.pdf"
        webbrowser.open(MThelp)  
        
    def run(self):
        import tempfile,glob,os
        import random
        import shutil
        import math
        from marinetools.GLCM.GLCM_dialog import GLCMDialog
        global selfMT

        # Create the dialog with elements (after translation) and keep reference
        # Only create GUI ONCE in callback, so that it will only load when the plugin is started
        if self.first_start == True:
            self.first_start = False
            self.dlg = GLCMDialog()
            selfMT = self
            self.dlg.inputFileDir.clicked.connect(GLCM.select_input_file) 
            self.dlg.outputFileDir.clicked.connect(GLCM.select_output_file) 
            self.dlg.GLCM_type.currentIndexChanged.connect(GLCM.indexChanged) 
            self.dlg.blocksize.textChanged.connect(GLCM.indexChanged) 
            self.dlg.sampsize.textChanged.connect(GLCM.indexChanged) 
            self.dlg.levels.textChanged.connect(GLCM.indexChanged) 
            self.dlg.angle.textChanged.connect(GLCM.indexChanged) 
            self.dlg.helpButton.clicked.connect(GLCM.help) 
            self.dlg.inputCombo.currentIndexChanged.connect(GLCM.indexChanged)

        
        layers = QgsProject.instance().mapLayers().values()
        self.dlg.inputCombo.clear() 
        self.dlg.inputCombo.addItems([layer.name() for layer in layers if str(layer.type())== "1" or str(layer.type())== "LayerType.Raster"])
        GLCM.indexChanged(self) 
        names = ["contrast","dissimilarity","homogeneity","energy","entropy","ASM","std","mean"]
        self.dlg.GLCM_type.clear() 
        # Populate the comboBox with names of all GLCM types   
        self.dlg.GLCM_type.addItems(names)

        # show the dialog
        self.dlg.show()
        # Run the dialog event loop
        result = self.dlg.exec_()
        # See if OK was pressed
        if result:
            selectedLayerIndex = self.dlg.inputCombo.currentIndex()
            currentText = selfMT.dlg.inputCombo.currentText()
            layers = QgsProject.instance().mapLayers().values()
            a=0
            filename="NULL"
            for layer in (layer1 for layer1 in layers if str(layer1.type())== "1" or str(layer1.type())== "LayerType.Raster"):
                if a == selectedLayerIndex:
                    filename = str(layer.source())
                a=a+1
            if currentText not in filename:
                filename = currentText
            inputFile = filename

            # Test input file to find if supported 
            from marinetools.fileTest import fileTest
            inputFile = fileTest.main(inputFile,"raster")
            if inputFile[0:5] == "Error":
                print(inputFile)
                return

            outputFile = self.dlg.outputFile.text()
            if outputFile == "":
                self.iface.messageBar().pushMessage("\nNo outputfile was given")
                return
            currentGLCMIndex = self.dlg.GLCM_type.currentIndex()
            names = ["contrast","dissimilarity","homogeneity","energy","entropy","ASM","std","mean"]
            GLCM_type =names[currentGLCMIndex]
 
            blocksize = self.dlg.blocksize.text()
            if blocksize == "":
                blocksize = "5"
            sampsize = self.dlg.sampsize.text()
            if sampsize == "":
                sampsize = "1"
            greylevels = self.dlg.levels.text()
            if greylevels == "":
                greylevels = "8"
            testangle = self.dlg.angle.text()
            if testangle == "":
                testangle = "0"
                
            try:
                import numpy as np
                import cv2
            except:
                import pip
                import subprocess,sys
                subprocess.call('python3 -m pip install opencv-python')
                import numpy as np
                import cv2

            try:
                import skimage
                from skimage.feature import graycomatrix, graycoprops
                from skimage import data
            except:
                import pip
                import subprocess,sys
                subprocess.call('python3 -m pip install scikit-image')
                
            try:
                import skimage
                from skimage.feature import graycomatrix, graycoprops
                from skimage import data
            except:
                import pip
                import subprocess,sys
                subprocess.call('python3 -m pip install skimage')
                import skimage
                from skimage.feature import graycomatrix, graycoprops
                from skimage import data
                
            #Print file report
            txtFile = open(str(OutputFilename[:-4]) + "_Info.txt", "w")
            txtFile.write("Script: Grey Level Co-occurence Matrices" "\n")
            txtFile.write("\n")
            txtFile.write("File Name: " + str(outputFile) + "\n")
            txtFile.write("Input File: " + os.path.basename(inputFile) + "\n")
            txtFile.write("\n")

            txtFile.write("Matrix used: " + GLCM_type + "\n")
            txtFile.write("Sampling Box Size      : " + str(blocksize) + "\n")
            txtFile.write("Sampling Distance      : " + str(sampsize) + "\n")
            txtFile.write("No. of grey-levels     : " + str(greylevels) + "\n")
            txtFile.write("Pixel pair angle       : " + str(testangle) + "\n")
            txtFile.write("\n")

            txtFile.close()


            infile = QgsRasterLayer(inputFile)
            crs = infile.crs()
            extent = infile.extent()

            result = processing.run("native:rescaleraster", {'INPUT':inputFile,'BAND':1,'MINIMUM':1,'MAXIMUM':255,'NODATA':0,'OUTPUT':'TEMPORARY_OUTPUT'})
            mytemp = result['OUTPUT']
            result = processing.run("gdal:translate", {'INPUT':mytemp,'TARGET_CRS':None,'NODATA':None,'COPY_SUBDATASETS':False,'OPTIONS':'','EXTRA':'','DATA_TYPE':3,'OUTPUT':'TEMPORARY_OUTPUT'})
            myRaster = result['OUTPUT']

            ks = int(blocksize)
            distance = int(sampsize)
            levels= int(greylevels)
            angle = int(testangle)
            mi, ma = 0, 255
            
            plugin_dir = os.path.dirname(__file__)
            gif_path = os.path.join(plugin_dir, "loading.gif")
            self.loading_screen = LoadingScreenDlg(gif_path)  # init loading dlg
            self.loading_screen.start_animation()  # start loading dlg

            #img = data.camera()
            import osgeo
            from osgeo import gdal
            ds = gdal.Open(myRaster)
            img = np.array(ds.GetRasterBand(1).ReadAsArray())
            h,w = img.shape
            
            if GLCM_type == "contrast":
                glcm_contrast = GLCM.fast_glcm_contrast(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_contrast, inputFile, outputFile)
            if GLCM_type == "dissimilarity":
                glcm_dissimilarity = GLCM.fast_glcm_dissimilarity(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_dissimilarity, inputFile, outputFile)
            if GLCM_type == "homogeneity":
                glcm_homogeneity= GLCM.fast_glcm_homogeneity(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_homogeneity, inputFile, outputFile)
            if GLCM_type == "energy":
                glcm_energy = GLCM.fast_glcm_energy(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_energy, inputFile, outputFile)
            if GLCM_type == "entropy":
                glcm_entropy = GLCM.fast_glcm_entropy(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_entropy, inputFile, outputFile)
            if GLCM_type == "ASM":
                glcm_ASM = GLCM.fast_glcm_ASM(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_ASM, inputFile, outputFile)
            if GLCM_type == "std":
                glcm_std = GLCM.fast_glcm_std(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_std, inputFile, outputFile)
            if GLCM_type == "mean":
                glcm_mean = GLCM.fast_glcm_mean(img, mi, ma, levels, ks, distance, angle)
                GLCM.NumPyArrayToRaster(glcm_mean, inputFile, outputFile)

            self.loading_screen.stop_animation()
                    
        pass



    def postExecute( parameters):
        """This method takes place after outputs are processed and added to the display."""
        return


    def fast_glcm( img, vmin=0, vmax=255, levels=8, kernel_size=5, distance=1.0, angle=0.0):
        '''
        Parameters
        ----------
        img: array_like, shape=(h,w), dtype=np.uint8
            input image
        vmin: int
            minimum value of input image
        vmax: int
            maximum value of input image
        levels: int
            number of grey-levels of GLCM
        kernel_size: int
            Patch size to calculate GLCM around the target pixel
        distance: float
            pixel pair distance offsets [pixel] (1.0, 2.0, and etc.)
        angle: float
            pixel pair angles [degree] (0.0, 30.0, 45.0, 90.0, and etc.)

        Returns
        -------
        Grey-level co-occurrence matrix for each pixels 
        shape = (levels, levels, h, w)
        '''

        mi, ma = vmin, vmax
        ks = kernel_size
        h,w = img.shape
        try:
            import numpy as np
            import cv2
        except:
            import pip
            pip.main(['install','opencv-python'])
            import numpy as np
            import cv2
        # digitize
        bins = np.linspace(mi, ma+1, levels+1)
        gl1 = np.digitize(img, bins) - 1

        # make shifted image
        dx = distance*np.cos(np.deg2rad(angle))
        dy = distance*np.sin(np.deg2rad(-angle))
        mat = np.array([[1.0,0.0,-dx], [0.0,1.0,-dy]], dtype=np.float32)
        gl2 = cv2.warpAffine(gl1, mat, (w,h), flags=cv2.INTER_NEAREST, borderMode=cv2.BORDER_REPLICATE)

        # make glcm
        glcm = np.zeros((levels, levels, h, w), dtype=np.uint8)
        for i in range(levels):
            for j in range(levels):
                mask = ((gl1==i) & (gl2==j))
                glcm[i,j, mask] = 1

        kernel = np.ones((ks, ks), dtype=np.uint8)
        for i in range(levels):
            for j in range(levels):
                glcm[i,j] = cv2.filter2D(glcm[i,j], -1, kernel)

        glcm = glcm.astype(np.float32)
        return glcm


    def fast_glcm_mean( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm mean
        '''
        
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        mean = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                mean += glcm[i,j] * i / (levels)**2

        return mean


    def fast_glcm_std(img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm std
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        mean = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                mean += glcm[i,j] * i / (levels)**2

        std2 = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                std2 += (glcm[i,j] * i - mean)**2

        std = np.sqrt(std2)
        return std


    def fast_glcm_contrast( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm contrast
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        cont = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                cont += glcm[i,j] * (i-j)**2

        return cont


    def fast_glcm_dissimilarity( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm dissimilarity
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        diss = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                diss += glcm[i,j] * np.abs(i-j)

        return diss


    def fast_glcm_homogeneity( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm homogeneity
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        homo = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                homo += glcm[i,j] / (1.+(i-j)**2)

        return homo


    def fast_glcm_ASM( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm asm, energy
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        asm = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                asm  += glcm[i,j]**2

        ene = np.sqrt(asm)
        #return asm, ene
        return asm

    def fast_glcm_energy( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm asm, energy
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        asm = np.zeros((h,w), dtype=np.float32)
        for i in range(levels):
            for j in range(levels):
                asm  += glcm[i,j]**2

        ene = np.sqrt(asm)
        #return asm, ene
        return ene


    def fast_glcm_max( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm max
        '''
        import numpy as np
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        max_  = np.max(glcm, axis=(0,1))
        return max_


    def fast_glcm_entropy( img, vmin=0, vmax=255, levels=8, ks=5, distance=1.0, angle=0.0):
        '''
        calc glcm entropy
        '''
        import numpy as np
        h,w = img.shape
        glcm = GLCM.fast_glcm(img, vmin, vmax, levels, ks, distance, angle)
        ent = np.zeros((h,w), dtype=np.float32)
        pnorm = glcm / np.sum(glcm, axis=(0,1), dtype=float) + 1./float(str(ks))**2
        ent  = np.sum(-pnorm * np.log(pnorm), axis=(0,1), dtype=float)
        return ent

    def NumPyArrayToRaster(myarray, inputFile, output_file):
        from osgeo import gdal, osr
        import numpy as np

        ds = gdal.Open(inputFile)
        band = ds.GetRasterBand(1)

        geotransform = ds.GetGeoTransform()
        wkt = ds.GetProjection()

        # Create gtif file
        driver = gdal.GetDriverByName("GTiff")

        dst_ds = driver.Create(output_file, band.XSize, band.YSize, 1, gdal.GDT_Float32)

        new_array = np.array(myarray, dtype=np.float32)

        #writting output raster
        dst_ds.GetRasterBand(1).WriteArray( new_array )
        #setting nodata value
        dst_ds.GetRasterBand(1).SetNoDataValue(0)
        #setting extension of output raster
        # top left x, w-e pixel resolution, rotation, top left y, rotation, n-s pixel resolution
        dst_ds.SetGeoTransform(geotransform)
        # setting spatial reference of output raster
        srs = osr.SpatialReference()
        srs.ImportFromWkt(wkt)
        dst_ds.SetProjection( srs.ExportToWkt() )
        #Close output raster dataset

        ds = None
        dst_ds = None
        fname = os.path.dirname(str(output_file))
        vlayer = QgsRasterLayer(str(output_file), str(output_file[len(fname)+1:-4]))
        QgsProject.instance().addMapLayer(vlayer)
        return



































