# /**********************************************************************
# *                                                                     *
# * Copyright (c) 2021 Javier Braña <javier.branagutierrez@gmail.com>  *
# *                                                                     *
# * This program is free software; you can redistribute it and/or modify*
# * it under the terms of the GNU Lesser General Public License (LGPL)  *
# * as published by the Free Software Foundation; either version 2 of   *
# * the License, or (at your option) any later version.                 *
# * for detail see the LICENCE text file.                               *
# *                                                                     *
# * This program is distributed in the hope that it will be useful,     *
# * but WITHOUT ANY WARRANTY; without even the implied warranty of      *
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
# * GNU Library General Public License for more details.                *
# *                                                                     *
# * You should have received a copy of the GNU Library General Public   *
# * License along with this program; if not, write to the Free Software *
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307*
# * USA                                                                 *
# *                                                                     *
# ***********************************************************************

import FreeCAD, FreeCADGui, Draft
from PySide import QtCore, QtGui, QtSvg
from PySide.QtCore import QT_TRANSLATE_NOOP

try:
    _fromUtf8 = QtCore.QString.fromUtf8
except AttributeError:
    def _fromUtf8(s):
        return s

import os, math
from PVPlantResources import DirIcons as DirIcons

class _TaskPanel:
    def __init__(self, obj = None):
        self.obj = None
        self.select = 0

        self.form = FreeCADGui.PySideUic.loadUi(os.path.dirname(__file__) + "/PVPlantCreateTerrainMesh.ui")
        self.form.buttonAdd.clicked.connect(self.add)

    def add(self):
        sel = FreeCADGui.Selection.getSelection()
        if len(sel) > 0:
            self.obj = sel[0]
            self.form.editCloud.setText(self.obj.Label)

    def MaxLength(self, P1, P2, P3, MaxlengthLE):
        """
        Calculation of the 2D length between triangle edges
        """

        p1 = FreeCAD.Vector(P1[0], P1[1], 0)
        p2 = FreeCAD.Vector(P2[0], P2[1], 0)
        p3 = FreeCAD.Vector(P3[0], P3[1], 0)

        List = [[p1, p2], [p2, p3], [p3, p1]]
        for i, j in List:
            vec = i.sub(j)
            if vec.Length > MaxlengthLE:
                return False

        return True

    def MaxAngle(self, P1, P2, P3, MaxAngleLE):
        """
        Calculation of the 2D angle between triangle edges
        """

        p1 = FreeCAD.Vector(P1[0], P1[1], 0)
        p2 = FreeCAD.Vector(P2[0], P2[1], 0)
        p3 = FreeCAD.Vector(P3[0], P3[1], 0)

        List = [[p1, p2, p3], [p2, p3, p1], [p3, p1, p2]]
        for j, k, l in List:
            vec1 = j.sub(k)
            vec2 = l.sub(k)
            radian = vec1.getAngle(vec2)
            if radian > MaxAngleLE:
                return False

        return True

    def accept(self):
        from datetime import datetime
        starttime = datetime.now()       

        import Part
        import numpy as np
        from scipy.spatial import Delaunay

        bnd = FreeCAD.ActiveDocument.Site.Terrain.CuttingBoundary.Shape
        if len(bnd.Faces) == 0:
            bnd = Part.Face(bnd)

        # Get user input
        MaxlengthLE = int(self.form.MaxlengthLE.text()) * 1000
        MaxAngleLE = math.radians(int(self.form.MaxAngleLE.text()))

        firstPoint = self.obj.Points.Points[0]
        nbase = FreeCAD.Vector(firstPoint.x, firstPoint.y, firstPoint.z)
        data = []
        for point in self.obj.Points.Points:
            tmp = FreeCAD.Vector(0, 0, 0).add(point)
            tmp.z = 0
            if bnd.isInside(tmp, 0, True):
                p = point - nbase
                data.append([float(p.x), float(p.y), float(p.z)])

        Data = np.array(data)
        data.clear()

        ''' not working:
        import multiprocessing
        rows = Data.shape[0]

        cpus = multiprocessing.cpu_count()
        steps = math.ceil(rows / cpus)

        for cpu in range(cpus - 1):
            start = steps * cpu
            end = steps * (cpu + 1)
            if end > rows:
                end = rows - start
            tmp = Data[start : end, :]
            p = multiprocessing.Process(target = Triangulate, args = (tmp,))
            p.start()
            p.join()
        return
        '''

        # TODO: si es muy grande, dividir el cálculo de la maya en varias etapas
        # Create delaunay triangulation
        tri = Delaunay(Data[:, :2])
        print("tiempo delaunay:", datetime.now() - starttime)

        faces = tri.simplices
        from stl import mesh
        wireframe = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
        for i, f in enumerate(faces):
            if self.MaxLength(Data[f[0]], Data[f[1]], Data[f[2]], MaxlengthLE) and \
               self.MaxAngle(Data[f[0]], Data[f[1]], Data[f[2]], MaxAngleLE):
                for j in range(3):
                    wireframe.vectors[i][j] = Data[f[j], :]

        import Mesh
        MeshObject = Mesh.Mesh(wireframe.vectors.tolist())
        MeshObject.Placement.move(nbase)
        MeshObject.harmonizeNormals()
        Surface = FreeCAD.ActiveDocument.addObject("Mesh::Feature", self.form.SurfaceNameLE.text())
        Surface.Mesh = MeshObject
        Surface.Label = self.form.SurfaceNameLE.text()

        shape = MeshToShape(MeshObject)
        import Part
        Part.show(shape)

        FreeCAD.ActiveDocument.recompute()
        FreeCADGui.Control.closeDialog()
        print(" --- Tiempo tardado:", datetime.now() - starttime)

    def reject(self):
        FreeCADGui.Control.closeDialog()


#----------------------------------------------------------------------------------------------
#from PySide.QtWidgets import QVBoxLayout, QLabel, QPushButton, QWidget, QMainWindow, QApplication

def Triangulate(Points, MaxlengthLE = 8000, MaxAngleLE = math.pi, use3d=True):
    import numpy as np
    from scipy.spatial import Delaunay
    from stl import mesh as stlmesh
    import Mesh

    if Points.__class__ is list:
        Points = np.array(Points)
    tri = Delaunay(Points[:, :2])
    faces = tri.simplices
    wireframe = stlmesh.Mesh(np.zeros(faces.shape[0], dtype=stlmesh.Mesh.dtype))

    for i, f in enumerate(faces):
        if MaxLength(Points[f[0]], Points[f[1]], Points[f[2]], MaxlengthLE) and \
           MaxAngle(Points[f[0]], Points[f[1]], Points[f[2]], MaxAngleLE):
            for j in range(3):
                wireframe.vectors[i][j] = Points[f[j], :]
                if not use3d:
                    wireframe.vectors[i][j][2] = 0

    MeshObject = Mesh.Mesh(wireframe.vectors.tolist())
    if len(MeshObject.Facets) == 0:
        return None
    MeshObject.harmonizeNormals()
    if MeshObject.Facets[0].Normal.z < 0:
        MeshObject.flipNormals()
    return MeshObject

def MaxLength(P1, P2, P3, MaxlengthLE):
    """ Calculation of the 2D length between triangle edges """
    p1 = FreeCAD.Vector(P1[0], P1[1], 0)
    p2 = FreeCAD.Vector(P2[0], P2[1], 0)
    p3 = FreeCAD.Vector(P3[0], P3[1], 0)

    List = [[p1, p2], [p2, p3], [p3, p1]]
    for i, j in List:
        vec = i.sub(j)
        if vec.Length > MaxlengthLE:
            return False
    return True

def MaxAngle(P1, P2, P3, MaxAngleLE):
    """ Calculation of the 2D angle between triangle edges """
    p1 = FreeCAD.Vector(P1[0], P1[1], 0)
    p2 = FreeCAD.Vector(P2[0], P2[1], 0)
    p3 = FreeCAD.Vector(P3[0], P3[1], 0)
    List = [[p1, p2, p3], [p2, p3, p1], [p3, p1, p2]]
    for j, k, l in List:
        vec1 = j.sub(k)
        vec2 = l.sub(k)
        radian = vec1.getAngle(vec2)
        if radian > MaxAngleLE:
            return False
    return True

def Open3DTriangle(point_cloud):
    import numpy as np
    import open3d as o3d

    '''
    input_path = "your_path_to_file/"
    output_path = "your_path_to_output_folder/"
    dataname = "sample.xyz"
    point_cloud = np.loadtxt(input_path + dataname, skiprows=1)
    '''
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector(point_cloud)

    pcd.normals = o3d.utility.Vector3dVector(np.zeros((1, 3)))
    pcd.estimate_normals()
    pcd.orient_normals_consistent_tangent_plane(100)

    mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(pcd,
                                                                             depth=8,
                                                                             width=0,
                                                                             scale=1.1,
                                                                             linear_fit=False,
                                                                             n_threads=8)
    o3d.visualization.draw_geometries([mesh])
    #bbox = pcd.get_axis_aligned_bounding_box()
    #p_mesh_crop = mesh.crop(bbox)
    return mesh

class _PVPlantCreateTerrainMesh:

    def GetResources(self):
        return {'Pixmap': str(os.path.join(DirIcons, "surface.svg")),
                'MenuText': QT_TRANSLATE_NOOP("PVPlant", "Create Surface"),
                'Accel': "C, S",
                'ToolTip': QT_TRANSLATE_NOOP("PVPlant", "Creates a surface form a cloud of points.")}

    def IsActive(self):
        return not FreeCAD.ActiveDocument is None

    def Activated(self):
        self.TaskPanel = _TaskPanel()
        FreeCADGui.Control.showDialog(self.TaskPanel)


if FreeCAD.GuiUp:
    FreeCADGui.addCommand('PVPlantCreateTerrainMesh', _PVPlantCreateTerrainMesh())
