import sys
import re
import numba as nb
import numpy as np
import scipy.spatial as spat


@nb.jit
def nb_dot(x, y):
    val = 0
    for x_i, y_i in zip(x, y):
        val += x_i * y_i
    return val


@nb.jit
def nb_cross(x, y):
    val = np.array([x[1] * y[2] - x[2] * y[1],
                    x[2] * y[0] - x[0] * y[2],
                    x[0] * y[1] - x[1] * y[0]])
    return val


@nb.jit
def r2_circumsphere_tetrahedron_single(a, b, c, d):
    ad = a - d
    bd = b - d
    cd = c - d

    ad2 = nb_dot(ad, ad)
    bd2 = nb_dot(bd, bd)
    cd2 = nb_dot(cd, cd)

    cross_1 = nb_cross(bd, cd)
    cross_2 = nb_cross(cd, ad)
    cross_3 = nb_cross(ad, bd)

    q = ad2 * cross_1 + bd2 * cross_2 + cd2 * cross_3
    p = 2 * np.abs(nb_dot(ad, cross_1))
    if p < 1e-10:
        return np.infty

    r2 = nb_dot(q, q) / p ** 2

    return r2


@nb.jit(nopython=True)
def r2_circumsphere_tetrahedron(a, b, c, d):
    len_a = len(a)
    r2 = np.zeros((len_a,))
    for i in range(len_a):
        r2[i] = r2_circumsphere_tetrahedron_single(a[i], b[i], c[i], d[i])
    return r2


@nb.jit
def get_faces(tetrahedron):
    faces = np.zeros((4, 3))
    for n, (i1, i2, i3) in enumerate([(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)]):
        faces[n] = tetrahedron[i1], tetrahedron[i2], tetrahedron[i3]
    return faces


def get_single_faces(triangulation):
    num_faces_single = 4
    num_tetrahedrons = triangulation.shape[0]
    num_faces = num_tetrahedrons * num_faces_single
    faces = np.zeros((num_faces, 3), np.int_)  # 3 is the dimension of the model
    mask = np.ones((num_faces,), np.bool_)
    for n in range(num_tetrahedrons):
        faces[num_faces_single * n: num_faces_single * (n + 1)] = get_faces(triangulation[n])
    orderlist = ["x{}".format(i) for i in range(faces.shape[1])]
    dtype_list = [(el, faces.dtype.str) for el in orderlist]
    faces.view(dtype_list).sort(axis=0)
    for k in range(num_faces - 1):
        if mask[k]:
            if np.all(faces[k] == faces[k + 1]):
                mask[k] = False
                mask[k + 1] = False
    single_faces = faces[mask]
    return single_faces


def get_alpha_shape(pointcloud, alpha):
    pointcloud = np.asarray(pointcloud)

    print(pointcloud,'\n',pointcloud.shape[1])

    assert pointcloud.ndim == 2
    assert pointcloud.shape[1] == 3, "for now, only 3-dimensional analysis is implemented"

    triangulation = spat.Delaunay(pointcloud)

    tetrahedrons = pointcloud[triangulation.simplices]  # remove this copy step, could be fatal
    radii2 = r2_circumsphere_tetrahedron(tetrahedrons[:, 0, :], tetrahedrons[:, 1, :], tetrahedrons[:, 2, :],
                                         tetrahedrons[:, 3, :])
    reduced_triangulation = triangulation.simplices[radii2 < alpha ** 2]
    del radii2, triangulation, tetrahedrons

    outer_triangulation = get_single_faces(reduced_triangulation)

    return outer_triangulation



def main():

    a = np.array([ 0.46421546,  0.50670312,  0.76935034,  0.75152631,  0.41167491,
        0.75871446,  0.54695894,  0.36280667,  0.08900743,  0.32331662])
    b = np.array([ 0.57476821,  0.34486742,  0.41292409,  0.95925496,  0.48904496,
        0.69459014,  0.92621067,  0.18750462,  0.28832875,  0.85921044])

    c = np.array([ 0.28563033,  0.32286582,  0.70590688,  0.21250345,  0.16658889,
        0.80191428,  0.22324416,  0.0867027 ,  0.87449543,  0.55770228])



    d=[]
    for i in range(len(a)):
        d.append((a[i],b[i],c[i]))

    e = get_alpha_shape(d,1)
    print(e,'\n',len(e))




if __name__ == '__main__':
    sys.argv[0] = re.sub(r'(-script\.pyw?|\.exe)?$', '', sys.argv[0])
    sys.exit(main())