import numpy as np
import math
from joblib import Parallel, delayed

def internal_raycast (equipment_data, xe, ye, ze, x1, y1, z1):
    x2 = xe + 0.5
    y2 = ye + 0.5
    z2 = ze + 0.5

    i = int(math.floor(x1))
    j = int(math.floor(y1))
    k = int(math.floor(z1))

    iend = int(math.floor(x2))
    jend = int(math.floor(y2))
    kend = int(math.floor(z2))

    di = 1 if x1 < x2 else -1 if x1 > x2 else 0
    dj = 1 if y1 < y2 else -1 if y1 > y2 else 0
    dk = 1 if z1 < z2 else -1 if z1 > z2 else 0

    deltatx = float('inf') if di == 0 else 1.0 / abs(x2 - x1)
    deltaty = float('inf') if dj == 0 else 1.0 / abs(y2 - y1)
    deltatz = float('inf') if dk == 0 else 1.0 / abs(z2 - z1)

    minx = math.floor(x1)
    miny = math.floor(y1)
    minz = math.floor(z1)

    maxx = minx + 1
    maxy = miny + 1
    maxz = minz + 1

    tx = deltatx * ((x1 - minx) if x1 > x2 else (maxx - x1))
    ty = deltaty * ((y1 - miny) if y1 > y2 else (maxy - y1))
    tz = deltatz * ((z1 - minz) if z1 > z2 else (maxz - z1))

    result = []
    result_full = []

    stop = False
    while True:
        if (tx <= ty and tx <= tz):
            if (i == iend):
                break
            tx = tx + deltatx
            i = i + di
        elif (ty <= tz):
            if (j == jend):
                break
            ty = ty + deltaty
            j = j + dj
        else:
            if (k == kend):
                break
            tz = tz + deltatz
            k = k + dk

        result_full.append((i,j,k,True))
        if (not stop):
            result.append((i, j, k, True))
            if (equipment_data[i][j][k]):
                stop = True

    return [result, result_full]

def raycast (vertex, shell, equipment):

    dim = shell.dims[0]

    result = np.zeros(shape=(dim,dim,dim), dtype=np.bool_)
    result_full = np.zeros(shape=(dim,dim,dim), dtype=np.bool_)

    result = result.tolist()
    result_full = result_full.tolist()
    equipment_data = equipment.data.tolist()

    xs, ys, zs = vertex
    x1,y1,z1 = [e+0.5 for e in vertex]

    result[xs][ys][zs] = True
    result_full[xs][ys][zs] = True

    r = Parallel(n_jobs=len(shell.data[0]))(delayed(internal_raycast)(equipment_data, shell.data[0][index], shell.data[1][index], shell.data[2][index], x1, y1, z1) for index in range (len(shell.data[0]))
    results_multi, results_full_multi = zip(*r)

        for res in results_multi:
            #for res in result_multi:
            result[res[0]][res[1]][res[2]] = res[3]

        for res_f in results_full_multi:
            #for res_f in result_full_multi:
            result_full[res_f[0]][res_f[2]][res_f[2]] = res_f[3]

    return [np.asarray(result), np.asarray(result_full)]

def voxelCoordinate (vertex, voxel_length, translate):
    return [int((v-t)/voxel_length - 0.5) for v,t in zip(vertex,translate)]

def is_inside (voxel, dim):
    return (voxel[0]>=0 and voxel[1]>=0 and voxel[2]>=0 and voxel[0]<dim and voxel[1]<dim and voxel[2]<dim)

def bresenham3d (vertex, shell, equipment):

    dim = shell.dims[0]

    result = np.zeros(shape=(dim,dim,dim), dtype=np.bool_)
    result_full = np.zeros(shape=(dim,dim,dim), dtype=np.bool_)

    result = result.tolist()
    result_full = result_full.tolist()
    equipment_data = equipment.data.tolist()

    for index in range (len(shell.data[0])):

        #point = vertex[:]
        point = [vertex[0],vertex[1],vertex[2]]
        dx = shell.data[0][index] - point[0]
        dy = shell.data[1][index] - point[1]
        dz = shell.data[2][index] - point[2]

        x_inc = -1 if dx<0 else 1
        l = abs(dx)
        y_inc = -1 if dy<0 else 1
        m = abs(dy)
        z_inc = -1 if dz<0 else 1
        n = abs(dz)
        dx2 = l << 1
        dy2 = m << 1
        dz2 = n << 1

        stop = False
        if ((l >= m) and (l >= n)):
            err_1 = dy2 - l
            err_2 = dz2 - l
            for i in range(l):
                if (not stop):
                    result[point[0]][point[1]][point[2]] = True
                    if (equipment_data[point[0]][point[1]][point[2]]):
                        stop = True
                result_full[point[0]][point[1]][point[2]] = True

                if (err_1 > 0):
                    point[1] = point[1] + y_inc
                    err_1 = err_1 - dx2

                if (err_2 > 0):
                    point[2] = point[2] + z_inc
                    err_2 = err_2 - dx2

                err_1 = err_1 + dy2
                err_2 = err_2 + dz2
                point[0] = point[0] + x_inc

        elif ((m >= l) and (m >= n)):
            err_1 = dx2 - m
            err_2 = dz2 - m
            for i in range(m):
                if (not stop):
                    result[point[0]][point[1]][point[2]] = True
                    if (equipment_data[point[0]][point[1]][point[2]]):
                        stop = True
                result_full[point[0]][point[1]][point[2]] = True

                if (err_1 > 0):
                    point[0] = point[0] + x_inc
                    err_1 = err_1 - dy2

                if (err_2 > 0):
                    point[2] = point[2] + z_inc
                    err_2 = err_2 - dy2

                err_1 = err_1 + dx2
                err_2 = err_2 + dz2
                point[1] = point[1] + y_inc

        else:
            err_1 = dy2 - n
            err_2 = dx2 - n
            for i in range(n):
                if (not stop):
                    result[point[0]][point[1]][point[2]] = True
                    if (equipment_data[point[0]][point[1]][point[2]]):
                        stop = True
                result_full[point[0]][point[1]][point[2]] = True

                if (err_1 > 0):
                    point[1] = point[1] + y_inc
                    err_1 = err_1 - dz2

                if (err_2 > 0):
                    point[0] = point[0] + x_inc
                    err_2 = err_2 - dz2

                err_1 = err_1 + dy2
                err_2 = err_2 + dx2
                point[2] = point[2] + z_inc

        if (not stop):
            result[point[0]][point[1]][point[2]] = True
            if (equipment.data[point[0]][point[1]][point[2]]):
                stop = True
        result_full[point[0]][point[1]][point[2]] = True

    return [np.asarray(result), np.asarray(result_full)]