from scipy.spatial import Delaunay
from itertools import combinations
import numpy as np
from im2mesh.utils import voxels


class MultiGridExtractor(object):
    def __init__(self, resolution0, threshold):
        # Attributes
        self.resolution = resolution0
        self.threshold = threshold

        # Voxels are active or inactive,
        # values live on the space between voxels and are either
        # known exactly or guessed by interpolation (unknown)
        shape_voxels = (resolution0,) * 3
        shape_values = (resolution0 + 1,) * 3
        self.values = np.empty(shape_values)
        self.value_known = np.full(shape_values, False)
        self.voxel_active = np.full(shape_voxels, True)

    def query(self):
        # Query locations in grid that are active but unkown
        idx1, idx2, idx3 = np.where(
            ~self.value_known & self.value_active
        )
        points = np.stack([idx1, idx2, idx3], axis=-1)
        return points

    def update(self, points, values):
        # Update locations and set known status to true
        idx0, idx1, idx2 = points.transpose()
        self.values[idx0, idx1, idx2] = values
        self.value_known[idx0, idx1, idx2] = True

        # Update activity status of voxels accordings to new values
        self.voxel_active = ~self.voxel_empty
        # (
        #     # self.voxel_active &
        #     self.voxel_known & ~self.voxel_empty
        # )

    def increase_resolution(self):
        self.resolution = 2 * self.resolution
        shape_values = (self.resolution + 1,) * 3

        value_known = np.full(shape_values, False)
        value_known[::2, ::2, ::2] = self.value_known
        values = upsample3d_nn(self.values)
        values = values[:-1, :-1, :-1]

        self.values = values
        self.value_known = value_known
        self.voxel_active = upsample3d_nn(self.voxel_active)

    @property
    def occupancies(self):
        return (self.values < self.threshold)

    @property
    def value_active(self):
        value_active = np.full(self.values.shape, False)
        # Active if adjacent to active voxel
        value_active[:-1, :-1, :-1] |= self.voxel_active
        value_active[:-1, :-1, 1:] |= self.voxel_active
        value_active[:-1, 1:, :-1] |= self.voxel_active
        value_active[:-1, 1:, 1:] |= self.voxel_active
        value_active[1:, :-1, :-1] |= self.voxel_active
        value_active[1:, :-1, 1:] |= self.voxel_active
        value_active[1:, 1:, :-1] |= self.voxel_active
        value_active[1:, 1:, 1:] |= self.voxel_active

        return value_active

    @property
    def voxel_known(self):
        value_known = self.value_known
        voxel_known = voxels.check_voxel_occupied(value_known)
        return voxel_known

    @property
    def voxel_empty(self):
        occ = self.occupancies
        return ~voxels.check_voxel_boundary(occ)


def upsample3d_nn(x):
    xshape = x.shape
    yshape = (2*xshape[0], 2*xshape[1], 2*xshape[2])

    y = np.zeros(yshape, dtype=x.dtype)
    y[::2, ::2, ::2] = x
    y[::2, ::2, 1::2] = x
    y[::2, 1::2, ::2] = x
    y[::2, 1::2, 1::2] = x
    y[1::2, ::2, ::2] = x
    y[1::2, ::2, 1::2] = x
    y[1::2, 1::2, ::2] = x
    y[1::2, 1::2, 1::2] = x

    return y


class DelauneyMeshExtractor(object):
    """Algorithm for extacting meshes from implicit function using
    delauney triangulation and random sampling."""
    def __init__(self, points, values, threshold=0.):
        self.points = points
        self.values = values
        self.delaunay = Delaunay(self.points)
        self.threshold = threshold

    def update(self, points, values, reduce_to_active=True):
        # Find all active points
        if reduce_to_active:
            active_simplices = self.active_simplices()
            active_point_idx = np.unique(active_simplices.flatten())
            self.points = self.points[active_point_idx]
            self.values = self.values[active_point_idx]

        self.points = np.concatenate([self.points, points], axis=0)
        self.values = np.concatenate([self.values, values], axis=0)
        self.delaunay = Delaunay(self.points)

    def extract_mesh(self):
        threshold = self.threshold
        vertices = []
        triangles = []
        vertex_dict = dict()

        active_simplices = self.active_simplices()
        active_simplices.sort(axis=1)
        for simplex in active_simplices:
            new_vertices = []
            for i1, i2 in combinations(simplex, 2):
                assert(i1 < i2)
                v1 = self.values[i1]
                v2 = self.values[i2]
                if (v1 < threshold) ^ (v2 < threshold):
                    # Subdivide edge
                    vertex_idx = vertex_dict.get((i1, i2), len(vertices))
                    vertex_idx = len(vertices)
                    if vertex_idx == len(vertices):
                        tau = (threshold - v1) / (v2 - v1)
                        assert(0 <= tau <= 1)
                        p = (1 - tau) * self.points[i1] + tau * self.points[i2]
                        vertices.append(p)
                        vertex_dict[i1, i2] = vertex_idx
                    new_vertices.append(vertex_idx)

            assert(len(new_vertices) in (3, 4))
            p0 = self.points[simplex[0]]
            v0 = self.values[simplex[0]]
            if len(new_vertices) == 3:
                i1, i2, i3 = new_vertices
                p1, p2, p3 = vertices[i1], vertices[i2], vertices[i3]
                vol = get_tetrahedon_volume(np.asarray([p0, p1, p2, p3]))
                if vol * (v0 - threshold) <= 0:
                    triangles.append((i1, i2, i3))
                else:
                    triangles.append((i1, i3, i2))
            elif len(new_vertices) == 4:
                i1, i2, i3, i4 = new_vertices
                p1, p2, p3, p4 = \
                    vertices[i1], vertices[i2], vertices[i3], vertices[i4]
                vol = get_tetrahedon_volume(np.asarray([p0, p1, p2, p3]))
                if vol * (v0 - threshold) <= 0:
                    triangles.append((i1, i2, i3))
                else:
                    triangles.append((i1, i3, i2))

                vol = get_tetrahedon_volume(np.asarray([p0, p2, p3, p4]))
                if vol * (v0 - threshold) <= 0:
                    triangles.append((i2, i3, i4))
                else:
                    triangles.append((i2, i4, i3))

        vertices = np.asarray(vertices, dtype=np.float32)
        triangles = np.asarray(triangles, dtype=np.int32)

        return vertices, triangles

    def query(self, size):
        active_simplices = self.active_simplices()
        active_simplices_points = self.points[active_simplices]
        new_points = sample_tetraheda(active_simplices_points, size=size)
        return new_points

    def active_simplices(self):
        occ = (self.values >= self.threshold)
        simplices = self.delaunay.simplices
        simplices_occ = occ[simplices]

        active = (
            np.any(simplices_occ, axis=1) & np.any(~simplices_occ, axis=1)
        )

        simplices = self.delaunay.simplices[active]
        return simplices


def sample_tetraheda(tetraheda_points, size):
    N_tetraheda = tetraheda_points.shape[0]
    volume = np.abs(get_tetrahedon_volume(tetraheda_points))
    probs = volume / volume.sum()

    tetraheda_rnd = np.random.choice(range(N_tetraheda), p=probs, size=size)
    tetraheda_rnd_points = tetraheda_points[tetraheda_rnd]
    weights_rnd = np.random.dirichlet([1, 1, 1, 1], size=size)
    weights_rnd = weights_rnd.reshape(size, 4, 1)
    points_rnd = (weights_rnd * tetraheda_rnd_points).sum(axis=1)
    # points_rnd = tetraheda_rnd_points.mean(1)

    return points_rnd


def get_tetrahedon_volume(points):
    vectors = points[..., :3, :] - points[..., 3:, :]
    volume = 1/6 * np.linalg.det(vectors)
    return volume
