import math


def trilinear_interpolation(idx, grid):
    """
    Return trilinear interpolated value for a voxel with index idx within the grid
    """

    int_idx = [math.floor(f) for f in idx]
    frac_idx = [f % 1 for f in idx]

    l1 = [[[0 for x in range(2)] for x in range(2)] for x in range(2)]
    for x in range(0, 2):
        for y in range(0, 2):
            for z in range(0, 2):
                l1[x][y][z] = grid[int_idx[0] + x, int_idx[1] + y, int_idx[2] + z]

    l2 = [[0 for x in range(2)] for x in range(2)]
    for y in range(0, 2):
        for z in range(0, 2):
            l2[y][z] = l1[0][y][z] * (1 - frac_idx[0]) + l1[1][y][z] * frac_idx[0]

    l3 = [0 for x in range(2)]
    for z in range(0, 2):
        l3[z] = l2[0][z] * (1 - frac_idx[1]) + l2[1][z] * frac_idx[1]

    return l3[0] * (1 - frac_idx[2]) + l3[1] * frac_idx[2]