
import json

from geogst.core.mathematics.deformations import *
from geogst.core.geometries.shapes.segments import Segment
from geogst.core.utils.arrays import *


class TriangBeam(object):
    """
    represents a 'fascio', a 2D semi-infinite geometrical object,
    defined by an apex (a Point object) and two semi-infinite segments, originating
    from the apex and defined by two versors (Vect objects, with init name 'versor_1' and 'versor_2').
    Its maximum width
    """

    def __init__(self, apex_pt3d, vector_1, vector_2):
        """
        assert almost_zero(versor_1.length() - 1.0)
        assert almost_zero(versor_2.length() - 1.0)
        """

        self._apex = apex_pt3d
        self._versor_1 = vector_1.versor()
        self._versor_2 = vector_2.versor()

    def fangle_degr(self):
        """
        rot_angle 'sotteso' by the 'fascio'
        in the 0 - 180 degrees range
        """

        return self._versor_1.angle_with(self._versor_2)

    def point_fangles_degr(self, pt_3d):
        """
        angles
        """

        vector_pt = Segment(self._apex, pt_3d).as_vector3d()

        angle_side_1 = self._versor_1.angle_with(vector_pt)
        angle_side_2 = self._versor_2.angle_with(vector_pt)

        return angle_side_1, angle_side_2

    def is_within_fascio(self, pt_3d):

        apertura = self.fangle_degr()

        if apertura >= 180.0:
            raise ValueError("apertura must be < 180.0")

        ang1, ang2 = self.point_fangles_degr(pt_3d)
        angle_sum = ang1 + ang2

        return almost_zero(apertura - angle_sum)


class CartesianTriangle(object):

    def __init__(self, pt_3d_1, pt_3d_2, pt_3d_3):

        self._pt_1 = pt_3d_1
        self._pt_2 = pt_3d_2
        self._pt_3 = pt_3d_3

    def is_pt_within(self, pt_3d):

        def versor3d(pt_1, pt_2):

            return Segment(pt_1, pt_2).as_vector3d().versor()

        def is_pt_in_fascio(pt_1, pt_2, pt_3):

            apex = pt_1
            versor_1 = versor3d(pt_1, pt_2)
            versor_2 = versor3d(pt_1, pt_3)

            fascio = TriangBeam(apex, versor_1, versor_2)
            if not fascio.is_within_fascio(pt_3d):
                return False
            else:
                return True

        if not (is_pt_in_fascio(self._pt_1, self._pt_2, self._pt_3) and
                    is_pt_in_fascio(self._pt_2, self._pt_1, self._pt_3)):
            return False
        else:
            return True


def geographic_scale_matrix(a_range, b_range, grid_height, grid_width):

    if a_range <= 0.0:
        raise ValueError("a_range must be > 0.0")
    if b_range <= 0.0:
        raise ValueError("b_range must be > 0.0")
    if grid_height <= 0.0:
        raise ValueError("grid_height must be > 0.0")
    if grid_width <= 0.0:
        raise ValueError("grid_width must be > 0.0")

    sx = grid_width / a_range
    sy = grid_height / b_range
    sz = 1

    return np.array([(sx, 0.0, 0.0), (0.0, sy, 0.0), (0.0, 0.0, sz)])


def geographic_rotation_matrix(grid_rot_angle_degr):
    grid_rot_angle_rad = radians(grid_rot_angle_degr)
    sin_rot_angle = sin(grid_rot_angle_rad)
    cos_rot_angle = cos(grid_rot_angle_rad)

    return np.array([(cos_rot_angle, -sin_rot_angle, 0.0),
                     (sin_rot_angle, cos_rot_angle, 0.0),
                     (0.0, 0.0, 1.0)])


def geographic_offset(transformation_matrix, llc_point_matr, llc_point_geog):

    return llc_point_geog - np.dot(transformation_matrix, llc_point_matr)


def geosurface_export_vtk(output_filepath, geodata):

    geosurface_XYZ, grid_dims = geodata
    X, Y, Z = geosurface_XYZ

    X_arr = np.array(X, dtype=float)
    Y_arr = np.array(Y, dtype=float)
    Z_arr = np.array(Z, dtype=float)

    n_points = np.size(X_arr)

    n_rows, n_cols = grid_dims

    with open(output_filepath, 'w') as outfile:

        outfile.write('# vtk DataFile Version 2.0\n')
        outfile.write('Geosurface - qgSurf vers. 0.3.0\n')
        outfile.write('ASCII\n')
        outfile.write('\nDATASET POLYDATA\n')

        outfile.write('POINTS %d float\n' % n_points)
        for n in range(n_points):
            outfile.write('%.4f %.4f %.4f\n' % (X_arr[n], Y_arr[n], Z_arr[n]))

        outfile.write('\n')

        outfile.write('TRIANGLE_STRIPS %d %d\n' % (n_cols - 1, (n_cols - 1) * (1 + n_rows * 2)))

        num_points_strip = n_rows * 2
        for l in range(n_cols - 1):
            triangle_strip_string = "%d " % num_points_strip
            for p in range(n_rows):
                triangle_strip_string += "%d %d " % ((l + 1) * n_rows + p, l * n_rows + p)
            triangle_strip_string += "\n"
            outfile.write(triangle_strip_string)


def geosurface_export_grass(output_filepath, geodata):
    # Save in Grass format

    geosurface_XYZ, grid_dims = geodata
    X, Y, Z = geosurface_XYZ

    X_arr = np.array(X, dtype=float)
    Y_arr = np.array(Y, dtype=float)
    Z_arr = np.array(Z, dtype=float)

    n_rows, n_cols = grid_dims

    with open(output_filepath, 'w') as outfile:
        outfile.write('VERTI:\n')
        for l in range(n_cols - 1):
            for p in range(n_rows - 1):
                start_point_ndx = l * n_rows + p
                forward_line_point_ndx = start_point_ndx + n_rows
                outfile.write('F 4\n')
                outfile.write(
                    ' %.4f %.4f %.4f\n' % (X_arr[start_point_ndx], Y_arr[start_point_ndx], Z_arr[start_point_ndx]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[start_point_ndx + 1], Y_arr[start_point_ndx + 1], Z_arr[start_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))
                outfile.write(
                    ' %.4f %.4f %.4f\n' % (X_arr[start_point_ndx], Y_arr[start_point_ndx], Z_arr[start_point_ndx]))
                outfile.write('F 4\n')
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[start_point_ndx + 1], Y_arr[start_point_ndx + 1], Z_arr[start_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx + 1], Y_arr[forward_line_point_ndx + 1],
                    Z_arr[forward_line_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))

'''
def geosurface_export_esri_generate(output_filepath, geodata):
    # Save geosurface (GAS) in Esri generate  format

    geosurface_XYZ, grid_dims = geodata
    X, Y, Z = geosurface_XYZ

    X_arr = np.array(X, dtype=float)
    Y_arr = np.array(Y, dtype=float)
    Z_arr = np.array(Z, dtype=float)

    n_rows, n_cols = grid_dims

    progr_id = 0
    with open(output_filepath, 'w') as outfile:
        outfile.write('VERTI:\n')
        for l in range(n_cols - 1):
            for p in range(n_rows - 1):
                start_point_ndx = l * n_rows + p
                forward_line_point_ndx = start_point_ndx + n_rows
                progr_id += 1
                outfile.write('%d\n' % progr_id)
                outfile.write(
                    ' %.4f %.4f %.4f\n' % (X_arr[start_point_ndx], Y_arr[start_point_ndx], Z_arr[start_point_ndx]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[start_point_ndx + 1], Y_arr[start_point_ndx + 1], Z_arr[start_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))
                outfile.write(
                    ' %.4f %.4f %.4f\n' % (X_arr[start_point_ndx], Y_arr[start_point_ndx], Z_arr[start_point_ndx]))
                outfile.write('END\n')
                progr_id += 1
                outfile.write('%d\n' % progr_id)
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[start_point_ndx + 1], Y_arr[start_point_ndx + 1], Z_arr[start_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx + 1], Y_arr[forward_line_point_ndx + 1],
                    Z_arr[forward_line_point_ndx + 1]))
                outfile.write(' %.4f %.4f %.4f\n' % (
                    X_arr[forward_line_point_ndx], Y_arr[forward_line_point_ndx], Z_arr[forward_line_point_ndx]))
                outfile.write('END\n')
        outfile.write('END\n')
'''

def geosurface_save_gas(output_filepath, geodata):

    with open(output_filepath, 'w') as outfile:
        json.dump(geodata, outfile)


def geosurface_read_gas_input(infile_path):

    try:
        with open(infile_path, 'r') as infile:
            input_geosurface = json.load(infile)
    except:
        raise Exception("Check input file name")

    src_analytical_params = input_geosurface['analytical surface']
    src_geographical_params = input_geosurface['geographical params']
    try:
        src_deformational_params = input_geosurface['deformational params']
    except:
        src_deformational_params = []

    return src_analytical_params, src_geographical_params, src_deformational_params


if __name__ == "__main__":
    import doctest

    doctest.testmod()
