# -*- coding: utf-8 -*-
"""
/***************************************************************************
 DemShading - ambient occlusion
 This algorithm simulates ambiental occlusion effect over an elevation model (DEM)
 Generated by Plugin Builder: http://g-sherman.github.io/Qgis-Plugin-Builder/
                              -------------------
        begin                : 2020-02-20
        copyright            : (C) 2020 by Zoran Čučković
 ***************************************************************************/
/***************************************************************************
 *                                                                         *
 *   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.                                   *
 *                                                                         *
 ***************************************************************************/
"""

__author__ = 'Zoran Čučković'
__date__ = '2020-02-05'
__copyright__ = '(C) 2020 by Zoran Čučković'
# This will get replaced with a git SHA1 when you do a git archive
__revision__ = '$Format:%H$'
from os import sys, path
from PyQt5.QtCore import QCoreApplication
from qgis.core import (QgsProcessing,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterRasterLayer,
                       QgsProcessingParameterRasterDestination,
                        QgsProcessingParameterBoolean,
                      QgsProcessingParameterNumber,
                       QgsProcessingParameterEnum,
                       QgsProcessingUtils,
                        QgsRasterBandStats,
                       QgsSingleBandGrayRenderer,
                       QgsContrastEnhancement
                        )
from processing.core.ProcessingConfig import ProcessingConfig
try:
    from osgeo import gdal
except ImportError:
    import gdal
import numpy as np
from .modules import Raster as rs
from .modules.helpers import view, window_loop, filter3, median_filter
from qgis.core import QgsMessageLog # for testing
class OcclusionAlgorithm(QgsProcessingAlgorithm):
    """
    This algorithm simulates ambient lighting over a raster DEM (in input). 
    """
    # Constants used to refer to parameters and outputs. They will be
    # used when calling the algorithm from another algorithm, or when
    # calling from the QGIS console.
    INPUT = 'INPUT'
    RADIUS= 'RADIUS'
    DENOISE = 'DENOISE'
    INVERT = 'INVERT'
    ANALYSIS_TYPE='ANALYSIS_TYPE'
    SYMMETRIC='SYMMETRIC'
    #RANGE = 'RANGE'
    OUTPUT = 'OUTPUT'
    ANALYSIS_TYPES = ['Sky-view','Openness']
    DENOISE_TYPES= ['None', 'Mean', 'Median', 'Mean and median']
    output_model = None #for post-processing
    def initAlgorithm(self, config):
        """
        Here we define the inputs and output of the algorithm, along
        with some other properties.
        """
        self.addParameter(
            QgsProcessingParameterRasterLayer(
                self.INPUT,
                self.tr('Digital elevation model')
            ) )
        
        self.addParameter(QgsProcessingParameterEnum (
            self.ANALYSIS_TYPE,
            self.tr('Analysis type'),
            self.ANALYSIS_TYPES,
            defaultValue=0))
        
        self.addParameter(QgsProcessingParameterEnum (
            self.ANALYSIS_TYPE,
            self.tr('Analysis type'),
            self.ANALYSIS_TYPES,
            defaultValue=0))
        
        self.addParameter(QgsProcessingParameterBoolean(
            self.SYMMETRIC,
            self.tr('Symmetric angles'),
            False, False)) 
        """
        self.addParameter(QgsProcessingParameterBoolean(
            self.RANGE,
            self.tr('Range (minimum to maximum angle'),
            False, False)) 
        """              
        self.addParameter(QgsProcessingParameterBoolean(
            self.INVERT,
            self.tr('Inverted DEM'),
            False, False)) 
                    
        self.addParameter(QgsProcessingParameterNumber(
            self.RADIUS,
            self.tr('Radius (pixels)'),
            0, # QgsProcessingParameterNumber.Integer = 0
            7, False, 0, 100))
        
        self.addParameter(QgsProcessingParameterEnum(
            self.DENOISE,
            self.tr('Denoise'),
            self.DENOISE_TYPES,
            defaultValue=0)) 
        
        self.addParameter(
            QgsProcessingParameterRasterDestination(
                self.OUTPUT,
            self.tr("Ambient occlusion")))
        
    def processAlgorithm(self, parameters, context, feedback):
        """
        Here is where the processing itself takes place.
        """
        elevation_model= self.parameterAsRasterLayer(parameters,self.INPUT, context)
        
        self.output_model = self.parameterAsOutputLayer(parameters,self.OUTPUT,context)

        #direction = self.parameterAsDouble(parameters,self.DIRECTION, context)
                
        denoise = self.parameterAsInt(parameters,self.DENOISE, context)  
        
        invert = self.parameterAsInt(parameters,self.INVERT, context)  
       
        radius =self.parameterAsInt(parameters,self.RADIUS, context)
        # 0 : sky view        
        openness = self.parameterAsInt(parameters,self.ANALYSIS_TYPE, context)
        symmetric = self.parameterAsInt(parameters,self.SYMMETRIC, context)
        
        # STILL TESTING
        difference = None #self.parameterAsInt(parameters,self.RANGE, context)
        
        dem = rs.Raster(elevation_model)
        
        err, fatal = dem.verify_raster()
        if err: feedback.reportError(err, fatalError = fatal)

        dem.set_output(self.output_model)
                               
        overlap = radius if not denoise else radius +1
             
        chunk_slice = (dem.ysize, dem.chunk_x + 2 * overlap)
                        
        mx_z = np.zeros(chunk_slice)
        mx_a = np.zeros(mx_z.shape)
      
        if symmetric: mx_b = mx_a
        else: mx_b = np.zeros(mx_z.shape)
        
        mx_a = np.zeros(mx_z.shape + ((radius ,) ))
        mx_b = np.zeros(mx_z.shape + ((radius ,) ))
        
        out =  np.zeros(mx_z.shape)
                
        # intialise the count of lines per pixel
        mx_cnt = np.ones(mx_z.shape)
        # set borders first  
        mx_cnt[:]= 5 if not symmetric else 3
        # main area : 8 lines per pixel (or 4 if symmetric algo)
        mx_cnt[1:-1, 1:-1] = 8 if not symmetric else 4
        # corners
        for v in [(0,0),(-1,-1),(0,-1), (-1,0)]: mx_cnt[v] = 3
               
        counter = 0
            
        for mx_view_in, gdal_take, mx_view_out, gdal_put in window_loop ( 
            shape = (dem.xsize, dem.ysize), 
            chunk = dem.chunk_x,
            overlap = overlap) :
            
            if counter : out[:] = 0 # reset
            
            dem.rst.ReadAsArray(*gdal_take, mx_z[mx_view_in]).astype(float)
            # NODATA : TODO !
            # mx_z[mx_z == nodata] = 0
            
            if invert : mx_z *= -1
            
            if denoise in [2, 3] : mx_z = median_filter(mx_z, radius= 3)
            if denoise in [1, 3] : mx_z = filter3(mx_z) # after the median filter
            

            # 8 standard lines, we use symmetry to optimise
            for dy, dx in [(0,1), (1,0), (1, -1), (1,1)]:
                
                if dx * dy : pix = np.sqrt( dem.pix_y**2 + dem.pix_x**2)
                else : pix = dem.pix_y if dx else dem.pix_x #swapped x, y
                
                for r in range (1, radius + 1): 
                    # we could probably sample over radius, not all pixels are needed...
                                                 
                    view_in, view_out = view(r * dx, r * dy, mx_z[mx_view_in].shape)

                    angles = mx_z[view_in] - mx_z[view_out]
                                                   # diagonals         
                    dist = r * pix
                    
                    angles /= dist 
     
                    mx_a[view_out][:,:,r-1] = angles
                    mx_b[view_in][:,:,r-1] = -angles
                    
                    # a patch for irregular pixels : take care of the length of the LOS
                    if dist > min(dem.pix_x * radius, dem.pix_y * radius) : break
               
                if not openness:  # sky view factor - remove negative angles
                    mx_a [mx_a < 0] = 0; mx_b[mx_b < 0] = 0
                 
                max_a = np.max(mx_a, axis = 2)
                max_b = np.max(mx_b, axis = 2)
                
                if difference : 
                    max_a -= np.min(mx_a, axis = 2)
                    max_b -= np.min(mx_b, axis = 2)
                              
                # average of angles: see Kokalj et al. 2011
                # these operations are costly, however ...
                if symmetric : 
                    # find the highest angle for each *pair* of LOS
                    np.maximum(max_a, max_b, out=max_a)
                    out += np.sin(np.arctan(max_a) )
                else: 
                    out += np.sin(np.arctan(max_a)) 
                    out += np.sin(np.arctan(max_b)) 
                
                mx_a[:] = 0 ; mx_b[:] = 0 # clean-up, remove old values 

                counter += 1
                
                feedback.setProgress(100 * dem.chunk_x * (counter/4) /  dem.xsize)
                if feedback.isCanceled(): return {}
            
            # this is a patch : last chunk is often spilling outside raster edge 
            # so, move the edge values to match raster edge
            end = gdal_take[2]           
            if end + gdal_take[0] == dem.xsize : 
                mx_cnt[:, end-1: end] = mx_cnt[:, -1:] 
            
            out /= mx_cnt
            out = 1 - out

            dem.add_to_buffer(out[mx_view_out], gdal_put)
            
        return {self.OUTPUT: self.output_model}

    def postProcessAlgorithm(self, context, feedback):

        output = QgsProcessingUtils.mapLayerFromString(self.output_model, context)
        provider = output.dataProvider()

        stats = provider.bandStatistics(1,QgsRasterBandStats.All,output.extent(),0)
        mean, sd = stats.mean, stats.stdDev
        
        rnd = QgsSingleBandGrayRenderer(provider, 1)
        ce = QgsContrastEnhancement(provider.dataType(1))
        ce.setContrastEnhancementAlgorithm(QgsContrastEnhancement.StretchToMinimumMaximum)
       
        ce.setMinimumValue(mean-3*sd)
        ce.setMaximumValue( mean+2*sd)

        rnd.setContrastEnhancement(ce)

        output.setRenderer(rnd)
        
        output.triggerRepaint()

        return {self.OUTPUT: self.output_model}

    def name(self):
        """
        Returns the algorithm name, used for identifying the algorithm. This
        string should be fixed for the algorithm, and must not be localised.
        The name should be unique within each provider. Names should contain
        lowercase alphanumeric characters only and no spaces or other
        formatting characters.
        """
        return 'Ambient occlusion'

    def displayName(self):
        """
        Returns the translated algorithm name, which should be used for any
        user-visible display of the algorithm name.
        """
        return self.tr(self.name()+ " (sky-view)")

    def tr(self, string):
        return QCoreApplication.translate('Processing', string)

    def shortHelpString(self):
        curr_dir = path.dirname(path.realpath(__file__))
        h = ( """
                Ambient occlusion of a locale is the proportion of ambient light that it recieves. This algorithm assumes equal light intensity from all directions (simple ambient lighting).
                Parameters:
                 - Analysis type: sky-view allows only for light sources situated above the horizontal plane (i.e. above the horizon), while openness takes into account all possible light sources.
                 - Symmetric: For each pair of opposite directions, take the one with higher horizon. Nice visual effect.
                 - Inverted DEM: Invert high and low values (multiply the DEM by -1) 
                 - Radius: The ambient occlusion is calculated within a defined radius for each raster pixel (computation time is directly dependent on the analysis radius).
                 - Denoise: Apply a smoothing filter.
                NB. This algorithm is made for terrain visualisation, it is not appropriate for precise calculation of solar exposition or of incident light.
                For more information, check <a href = "https://landscapearchaeology.org/qgis-terrain-shading/" >the manual</a>.
             
                If you find this tool useful, consider to :
                 
             <a href='https://ko-fi.com/D1D41HYSW' target='_blank'><img height='30' style='border:0px;height:36px;' src='%s/help/kofi2.webp' /></a>
            """) % curr_dir
		
        return self.tr(h)

    def createInstance(self):
        return OcclusionAlgorithm()
