# be lazy and use the modifier system to compute 3D voronoi cells.

import bpy
from mathutils import *
import random


def center_cube(sz):
    name = "cube %f" % sz
    rval = None # bpy.data.meshes.get(name)
    # I need a single-user mesh if I am applying modifiers, so do not try to reuse
    if rval is None:
        rval = bpy.data.meshes.new(name)
        verts = [
            [-sz,-sz,-sz],
            [sz,-sz,-sz],
            [-sz,sz,-sz],
            [sz,sz,-sz],
            [-sz,-sz,sz],
            [sz,-sz,sz],
            [-sz,sz,sz],
            [sz,sz,sz],
        ]
        faces = [
            [0,2,3,1],
            [0,1,5,4],
            [1,3,7,5],
            [3,2,6,7],
            [2,0,4,6],
            [4,5,7,6]
        ]
        rval.from_pydata(verts, [], faces)
    return rval

def get_half_plane(sz):
    name = "half plane %f" % sz
    rval = bpy.data.meshes.get(name)
    if rval is None:
        rval = bpy.data.meshes.new(name)
        verts = [
            [-sz,-sz,0],
            [sz,-sz,0],
            [-sz,sz,0],
            [sz,sz,0],
            [-sz,-sz,sz],
            [sz,-sz,sz],
            [-sz,sz,sz],
            [sz,sz,sz],
        ]
        faces = [
            [0,2,3,1],
            [0,1,5,4],
            [1,3,7,5],
            [3,2,6,7],
            [2,0,4,6],
            [4,5,7,6]
        ]
        rval.from_pydata(verts, [], faces)
    return rval

def matrix_for_voronoi_face(center, other):
    mid = (center+other)/2

    delta = (center-other).normalized()

    z_axis = delta
    if z_axis[2] == max(z_axis):
        x_axis = z_axis.cross([0,1,0]).normalized()
    else:
        x_axis = z_axis.cross([0,0,1]).normalized()
    y_axis = z_axis.cross(x_axis)

    rot = Matrix([x_axis,y_axis, z_axis]).transposed().to_4x4()

    return Matrix.Translation(mid)*rot




def fab_one_cell(coords, idx, scn, horizon=3):

    center = Vector(coords[idx])

    half_plane_mesh = get_half_plane(2*horizon)
    half_planes = []

    for i in range(len(coords)):
        if i==idx:
            continue

        vi = Vector(coords[i])
        if (max(center - vi)>horizon
            or
            min(center-vi) < -horizon):
            continue

        obj = bpy.data.objects.new("half plane", half_plane_mesh)
        obj.matrix_world = matrix_for_voronoi_face(center, vi)

        scn.objects.link(obj)

        half_planes.append(obj)
        obj.layers = [ i==5 for i in range(20)]


    prime = bpy.data.objects.new("prime", center_cube(2*horizon))
    prime.location = center
    scn.objects.link(prime)
    prime.layers = [ i==1 for i in range(20)]

    if False:

        for i in range(len(half_planes)):
            mod =prime.modifiers.new("boolean", 'BOOLEAN')
            mod.object = half_planes[i]
            half_planes[i].layers = [ i==5 for i in range(20)]
            #print(half_planes[i].layers)

        mesh = prime.to_mesh(scene = scn, apply_modifiers=True, settings ='PREVIEW')
        prime.data = mesh
        print(len(mesh.vertices))
        for i in range(len(half_planes)):
            scn.objects.unlink(half_planes[i])

    else:
        scn.objects.active = prime

        for i in range(len(half_planes)):
            mod =prime.modifiers.new("boolean", 'BOOLEAN')
            mod.object = half_planes[i]
            print(mod.name)
            bpy.ops.object.modifier_apply(modifier=mod.name)
            scn.objects.unlink(half_planes[i])

    return prime

def rand2(radius):
    return radius * 2 * (random.random()-0.5);

def displaced_vertex_grid( min_c, max_c, radius):

    rval = []
    for x in range(min_c[0], max_c[0]):
        for y in range(min_c[1], max_c[1]):
            for z in range(min_c[2], max_c[2]):
                xyz = [ x+rand2(radius), y+rand2(radius), z+rand2(radius)]
                rval.append(xyz)
    return rval

#
#

coords = [ [0,0,0] ]

#coords.append( [1,1,1] )
#coords.append( [1,1,0] )
radius = 10
coords.extend(displaced_vertex_grid( [-radius,-radius, -radius], [radius, radius, radius], 3) )

for i in range(min(200,len(coords))):
    center =coords[i]
    q = max( max(center), -min(center))
    if (q+2>radius):
        continue
    cell = fab_one_cell(coords, i, bpy.context.scene)
    cell.scale = [0.8,0.8,0.8]
    print(cell)