Well, the first problem you define is slightly more problematic, but the actual task you mention is easier.
A Voronoi tessellation is simply an L2 distance boundary. This means if you generate your random centers, you can voxelize it by assigning each voxel to the closest center.
It would look something like:
[X,Y,Z] = meshgrid(xgv,ygv,zgv)
Space = [X(:),Y(:),Z(:)];
IDX = knnsearch(CENTERS,Space);
Voxelized = reshape(IDX,sizeofVolume);