# -*- coding: utf-8 -*-
from .geomaticape_algorithm import GeomaticapeAlgorithm
"""
suavizar_poligono.py
====================
Suaviza los contornos de poligonos vectoriales utilizando el algoritmo PAEK
(Polynomial Approximation with Exponential Kernel), equivalente al
Smooth Polygon de ArcGIS.

Parametros:
    - Capa vectorial de poligonos (entrada)
    - Tolerancia de suavizado (unidades del mapa)
    - Correccion topologica (NO_CHECK / FLAG_ERRORS / RESOLVE_ERRORS)
    - Capa de salida suavizada

Autor : Geomatica Ambiental - https://www.geomatica.pe
Plugin: GeoRemote Tools v1.15
Grupo : PostProcesamiento
"""

import numpy as np

from qgis.core import (
    QgsProcessingParameterFeatureSource,
    QgsProcessingParameterNumber,
    QgsProcessingParameterEnum,
    QgsProcessingParameterFeatureSink,
    QgsProcessingParameterBoolean,
    QgsProcessingException,
    QgsProcessing,
    QgsFeatureSink,
    QgsWkbTypes,
    QgsGeometry,
    QgsPointXY,
    QgsFields,
    QgsField,
    QgsFeature,
)
from qgis.PyQt.QtCore import QVariant


# ---------------------------------------------------------------------------
# Constantes
# ---------------------------------------------------------------------------
TOPO_OPTIONS = [
    "NO_CHECK - No verificar errores topologicos",
    "FLAG_ERRORS - Detectar y reportar errores topologicos",
    "RESOLVE_ERRORS - Detectar y corregir errores topologicos",
]
TOPO_NO_CHECK    = 0
TOPO_FLAG_ERRORS = 1
TOPO_RESOLVE     = 2
MIN_VERTICES   = 4
MIN_AREA_RATIO = 1e-8


# ---------------------------------------------------------------------------
# PAEK - nucleo matematico
# ---------------------------------------------------------------------------

def _resample_ring(pts, step):
    """Remuestrea un anillo cerrado a intervalos regulares de longitud `step`."""
    pts = np.asarray(pts, dtype=np.float64)
    if len(pts) < 2:
        return pts
    diffs = np.diff(pts, axis=0)
    seg_len = np.hypot(diffs[:, 0], diffs[:, 1])
    cum = np.concatenate([[0.0], np.cumsum(seg_len)])
    total = cum[-1]
    if total < 1e-10:
        return pts
    n_pts = max(int(total / step), MIN_VERTICES)
    sample_d = np.linspace(0.0, total, n_pts, endpoint=False)
    rx = np.interp(sample_d, cum, pts[:, 0])
    ry = np.interp(sample_d, cum, pts[:, 1])
    rx = np.append(rx, rx[0])
    ry = np.append(ry, ry[0])
    return np.column_stack([rx, ry])


def _paek_smooth_ring(pts, tolerance):
    """
    Aplica el kernel exponencial PAEK a un anillo ya remuestrado.

    Cada punto suavizado es la media ponderada de todos los puntos del anillo,
    con peso  w_i = exp(-d_i / tolerance)  donde d_i es la distancia geodesica
    (arco) entre el punto objetivo y el punto i. La distancia de arco se calcula
    considerando el camino mas corto en ambos sentidos (wrap-around).

    Referencia teorica:
      McMaster, R.B. (1987). Automated line generalization.
      Cartographica 24(2):74-111.
    """
    pts = np.asarray(pts, dtype=np.float64)
    n = len(pts)
    if n < 3:
        return pts
    diffs = np.diff(pts, axis=0)
    seg_len = np.hypot(diffs[:, 0], diffs[:, 1])
    n_unique = n - 1
    cum = np.concatenate([[0.0], np.cumsum(seg_len[:-1])])
    total = cum[-1] + (seg_len[-1] if len(seg_len) > 0 else 0.0)
    if total < 1e-10:
        return pts
    sx = np.zeros(n_unique)
    sy = np.zeros(n_unique)
    for i in range(n_unique):
        raw = np.abs(cum - cum[i])
        d = np.minimum(raw, total - raw)
        w = np.exp(-d / tolerance)
        w_sum = w.sum()
        if w_sum < 1e-15:
            sx[i] = pts[i, 0]
            sy[i] = pts[i, 1]
        else:
            sx[i] = (w * pts[:n_unique, 0]).sum() / w_sum
            sy[i] = (w * pts[:n_unique, 1]).sum() / w_sum
    sx = np.append(sx, sx[0])
    sy = np.append(sy, sy[0])
    return np.column_stack([sx, sy])


def _smooth_ring(ring_pts, tolerance):
    """Pipeline: extraccion de coordenadas -> remuestreo -> PAEK -> lista QgsPointXY."""
    pts = np.asarray([[p.x(), p.y()] for p in ring_pts], dtype=np.float64)
    if len(pts) < MIN_VERTICES:
        return ring_pts
    step = max(tolerance / 4.0, 1e-6)
    resampled = _resample_ring(pts, step)
    if len(resampled) < MIN_VERTICES:
        return ring_pts
    smoothed = _paek_smooth_ring(resampled, tolerance)
    return [QgsPointXY(float(x), float(y)) for x, y in smoothed]


def _smooth_geometry(geom, tolerance):
    """Suaviza todos los anillos de una geometria Polygon o MultiPolygon."""
    geom_type = geom.wkbType()
    flat = QgsWkbTypes.flatType(geom_type)
    if flat == QgsWkbTypes.Type.Polygon:
        polygons = [geom.asPolygon()]
    elif flat == QgsWkbTypes.Type.MultiPolygon:
        polygons = geom.asMultiPolygon()
    else:
        return geom
    smoothed_polys = []
    for polygon in polygons:
        smoothed_rings = [_smooth_ring(ring, tolerance) for ring in polygon]
        smoothed_polys.append(smoothed_rings)
    if flat == QgsWkbTypes.Type.Polygon:
        return QgsGeometry.fromPolygonXY(smoothed_polys[0])
    return QgsGeometry.fromMultiPolygonXY(smoothed_polys)


# ---------------------------------------------------------------------------
# Correccion topologica
# ---------------------------------------------------------------------------

def _check_topology(geom):
    """Verifica si la geometria tiene errores topologicos."""
    if geom is None or geom.isEmpty():
        return False, "Geometria vacia"
    if geom.isGeosValid():
        return True, ""
    err = geom.isGeosValidReason() if hasattr(geom, "isGeosValidReason") else "Error topologico"
    return False, err


def _resolve_topology(geom):
    """Intenta reparar errores topologicos usando buffer(0) y makeValid."""
    fixed = geom.buffer(0, 5)
    if fixed and not fixed.isEmpty() and fixed.isGeosValid():
        return fixed
    if hasattr(geom, "makeValid"):
        valid = geom.makeValid()
        if valid and not valid.isEmpty() and valid.isGeosValid():
            return valid
    return None


# ---------------------------------------------------------------------------
# Clase principal
# ---------------------------------------------------------------------------

class SuavizarPoligono(GeomaticapeAlgorithm):
    _algorithm_name = "suavizar_poligono"
    _icon_name      = "poligonos_tabla.png"

    INPUT          = "INPUT"
    TOLERANCE      = "TOLERANCE"
    TOPO_METHOD    = "TOPO_METHOD"
    PRESERVE_SMALL = "PRESERVE_SMALL"
    OUTPUT         = "OUTPUT"
    OUTPUT_ERRORS  = "OUTPUT_ERRORS"

    def displayName(self):
        return self.tr("Suavizar poligono (PAEK)")

    def group(self):
        return self.tr("PostProcesamiento")

    def groupId(self):
        return "geomaticape_postprocesamiento"

    def shortHelpString(self):
        return """
<h3>Suavizar Poligono - PAEK</h3>
<b>Autor:</b> Geomatica Ambiental<br>
<b>Plugin:</b> GeoRemote Tools<br>
<b>Version:</b> 1.15<br><br>

<b>Descripcion:</b><br>
Suaviza los contornos de poligonos vectoriales utilizando el algoritmo
<b>PAEK</b> (Polynomial Approximation with Exponential Kernel), equivalente
al Smooth Polygon de ArcGIS.<br><br>

El algoritmo reemplaza cada vertice del anillo por la media ponderada de los
vertices vecinos, donde el peso decrece exponencialmente con la distancia de
arco. La <b>tolerancia de suavizado</b> controla el radio de influencia del
kernel: tolerancias mayores producen contornos mas suaves pero mas alejados
del original.<br><br>

<b>Parametros:</b>
<ul>
  <li><b>Capa de poligonos:</b> capa vectorial de entrada (Polygon / MultiPolygon).</li>
  <li><b>Tolerancia de suavizado</b> (unidades del mapa): radio de influencia del kernel.
      Equivale al Smoothing Tolerance de ArcGIS. Usar unidades del CRS
      (metros para proyectados, grados para geograficos).</li>
  <li><b>Correccion topologica:</b>
    <ul>
      <li><b>NO_CHECK:</b> no se verifica. Salida mas rapida.</li>
      <li><b>FLAG_ERRORS:</b> detecta y reporta poligonos invalidos en una capa adicional.</li>
      <li><b>RESOLVE_ERRORS:</b> detecta y corrige errores aplicando buffer(0) + makeValid.</li>
    </ul>
  </li>
  <li><b>Conservar poligonos que colapsan:</b> si el area suavizada es menor al 1%
      de la original, se mantiene la geometria original sin cambios.</li>
</ul>

<b>Capa de errores (FLAG_ERRORS):</b><br>
Se genera una capa de puntos con la ubicacion del centroide de cada poligono
invalido y el mensaje de error GEOS.<br><br>

<b>Recomendaciones de tolerancia:</b>
<ul>
  <li>Datos regionales 1:250000 : 500 - 2000 m</li>
  <li>Datos departamentales 1:100000: 100 - 500 m</li>
  <li>Datos locales 1:25000: 20 - 100 m</li>
  <li>Datos catastrales 1:5000: 1 - 10 m</li>
</ul>

<b>Web:</b> https://www.geomatica.pe/
"""

    def initAlgorithm(self, config=None):
        self.addParameter(
            QgsProcessingParameterFeatureSource(
                self.INPUT,
                self.tr("Capa de poligonos de entrada"),
                types=[QgsProcessing.SourceType.TypeVectorPolygon]
            )
        )
        self.addParameter(
            QgsProcessingParameterNumber(
                self.TOLERANCE,
                self.tr("Tolerancia de suavizado (unidades del mapa)"),
                type=QgsProcessingParameterNumber.Type.Double,
                defaultValue=100.0,
                minValue=0.001
            )
        )
        self.addParameter(
            QgsProcessingParameterEnum(
                self.TOPO_METHOD,
                self.tr("Correccion topologica"),
                options=TOPO_OPTIONS,
                defaultValue=TOPO_NO_CHECK
            )
        )
        self.addParameter(
            QgsProcessingParameterBoolean(
                self.PRESERVE_SMALL,
                self.tr("Conservar poligonos que colapsan tras el suavizado"),
                defaultValue=True
            )
        )
        self.addParameter(
            QgsProcessingParameterFeatureSink(
                self.OUTPUT,
                self.tr("Poligonos suavizados")
            )
        )
        param_err = QgsProcessingParameterFeatureSink(
            self.OUTPUT_ERRORS,
            self.tr("Capa de errores topologicos (solo FLAG_ERRORS)"),
            optional=True,
            createByDefault=False
        )
        self.addParameter(param_err)

    def processAlgorithm(self, parameters, context, feedback):
        source = self.parameterAsSource(parameters, self.INPUT, context)
        if source is None:
            raise QgsProcessingException(self.tr("No se pudo cargar la capa de entrada."))

        tolerance   = self.parameterAsDouble(parameters, self.TOLERANCE, context)
        topo_method = self.parameterAsEnum(parameters, self.TOPO_METHOD, context)
        preserve    = self.parameterAsBool(parameters, self.PRESERVE_SMALL, context)

        if tolerance <= 0:
            raise QgsProcessingException(self.tr("La tolerancia debe ser mayor que 0."))

        (sink, dest_id) = self.parameterAsSink(
            parameters, self.OUTPUT, context,
            source.fields(), source.wkbType(), source.sourceCrs()
        )
        if sink is None:
            raise QgsProcessingException(self.tr("No se pudo crear la capa de salida."))

        # Capa de errores (solo FLAG_ERRORS)
        sink_errors    = None
        dest_errors_id = None
        fields_err     = QgsFields()
        if topo_method == TOPO_FLAG_ERRORS:
            try:
                fields_err.append(QgsField("fid_origen", QVariant.LongLong))
                fields_err.append(QgsField("error_msg",  QVariant.String))
                (sink_errors, dest_errors_id) = self.parameterAsSink(
                    parameters, self.OUTPUT_ERRORS, context,
                    fields_err, QgsWkbTypes.Type.Point, source.sourceCrs()
                )
            except Exception:
                sink_errors = None

        total = source.featureCount()
        step  = 100.0 / total if total > 0 else 0
        n_smoothed = n_errors = n_resolved = n_preserved = n_skipped = 0

        for current, feat in enumerate(source.getFeatures()):
            if feedback.isCanceled():
                break
            feedback.setProgress(int(current * step))

            geom = feat.geometry()
            if geom is None or geom.isEmpty():
                sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                n_skipped += 1
                continue

            flat = QgsWkbTypes.flatType(geom.wkbType())
            if flat not in (QgsWkbTypes.Type.Polygon, QgsWkbTypes.Type.MultiPolygon):
                sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                n_skipped += 1
                continue

            area_original = geom.area()

            # --- Suavizado ---
            try:
                smoothed_geom = _smooth_geometry(geom, tolerance)
            except Exception as e:
                feedback.pushWarning(
                    self.tr(f"Feature {feat.id()}: error en suavizado ({e}). Se conserva original."))
                sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                n_skipped += 1
                continue

            if smoothed_geom is None or smoothed_geom.isEmpty():
                sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                n_preserved += 1
                continue

            # --- Proteger poligonos pequenos ---
            if preserve and area_original > 0:
                if smoothed_geom.area() < area_original * MIN_AREA_RATIO:
                    feedback.pushInfo(
                        self.tr(f"Feature {feat.id()}: area colapsada, se conserva original."))
                    sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                    n_preserved += 1
                    continue

            # --- Correccion topologica ---
            is_valid, error_msg = _check_topology(smoothed_geom)

            if not is_valid:
                n_errors += 1
                if topo_method == TOPO_FLAG_ERRORS:
                    feedback.pushWarning(
                        self.tr(f"Feature {feat.id()}: error topologico - {error_msg}"))
                    if sink_errors is not None:
                        ef = QgsFeature(fields_err)
                        ef.setGeometry(smoothed_geom.centroid())
                        ef.setAttribute("fid_origen", feat.id())
                        ef.setAttribute("error_msg",  error_msg)
                        sink_errors.addFeature(ef, QgsFeatureSink.Flag.FastInsert)
                elif topo_method == TOPO_RESOLVE:
                    fixed = _resolve_topology(smoothed_geom)
                    if fixed is not None:
                        smoothed_geom = fixed
                        n_resolved += 1
                        feedback.pushInfo(
                            self.tr(f"Feature {feat.id()}: error topologico corregido."))
                    else:
                        feedback.pushWarning(
                            self.tr(f"Feature {feat.id()}: no se pudo corregir ({error_msg}). Original conservado."))
                        sink.addFeature(feat, QgsFeatureSink.Flag.FastInsert)
                        n_preserved += 1
                        continue

            # --- Escribir feature suavizado ---
            out_feat = QgsFeature()
            out_feat.setGeometry(smoothed_geom)
            out_feat.setAttributes(feat.attributes())
            sink.addFeature(out_feat, QgsFeatureSink.Flag.FastInsert)
            n_smoothed += 1

        # --- Resumen ---
        feedback.pushInfo("-" * 55)
        feedback.pushInfo(self.tr("RESUMEN - Suavizar Poligono (PAEK)"))
        feedback.pushInfo(self.tr(f"  Tolerancia         : {tolerance} u.m."))
        feedback.pushInfo(self.tr(f"  Metodo topologico  : {TOPO_OPTIONS[topo_method].split(' - ')[0]}"))
        feedback.pushInfo(self.tr(f"  Total features     : {total}"))
        feedback.pushInfo(self.tr(f"  Suavizados OK      : {n_smoothed}"))
        if n_errors    > 0:
            feedback.pushInfo(self.tr(f"  Errores detectados : {n_errors}"))
        if n_resolved  > 0:
            feedback.pushInfo(self.tr(f"  Errores corregidos : {n_resolved}"))
        if n_preserved > 0:
            feedback.pushInfo(self.tr(f"  Conservados        : {n_preserved}"))
        if n_skipped   > 0:
            feedback.pushInfo(self.tr(f"  Omitidos           : {n_skipped}"))
        feedback.pushInfo("-" * 55)

        result = {self.OUTPUT: dest_id}
        if dest_errors_id:
            result[self.OUTPUT_ERRORS] = dest_errors_id
        return result

    def run(self):
        from qgis import processing
        processing.execAlgorithmDialog(self)
