from queue import PriorityQueue
import math

import random
from PIL import Image
from PIL import ImageDraw
import itertools
##
#   This file holds methods related to the generation of
# Voronoi Diagrams, and their Centroidal counterpart. 
# This approach is covered in A. Secord's 2002 paper:
# https://www.cs.ubc.ca/labs/imager/tr/2002/secord2002b/secord.2002b.pdf
##


##
#   This function returns a 2D array of the Discrete Voronoi Diagram
# generated by the genPts list.  The size of the diagram is described
# by the tuple sz.
##
def getVoronoi(genPts, sz):
  # defines a matrix to hold the cluster indeces
  matrix = [[-1 for x in range(sz[0])] for y in range(sz[1])]

  q = PriorityQueue()
  for k in range(len(genPts)):
    if (isValid(genPts[k], sz)):
      q.put((0, genPts[k], k))

  # I'll approximate the Voronoi diagram by a multisource BFS
  while (not q.empty()):
    val = q.get() # get closest point to *some* cluster source
    pt = val[1]
    k = val[2]

    # if it's already colored, don't recolor
    if (matrix[pt[1]][pt[0]] != -1):
      continue

    matrix[pt[1]][pt[0]] = k # color with cluster source

    for n in neighbors(pt, sz):
      if isValid(n, sz) and matrix[n[1]][n[0]] == -1:
        newD = dist(genPts[k], n)
        q.put((newD, n, k))

  return matrix
    
def dist(pt1, pt2):
  return math.sqrt((pt1[0] - pt2[0])**2 + (pt1[1] - pt2[1])**2)

def neighbors(pt, sz):
  dx = [0, 1,  0, -1]
  dy = [1, 0, -1,  0]

  retVal = []
  for i in range(len(dx)):
    retVal.append((pt[0] + dx[i], pt[1] + dy[i]))

  return retVal

def isValid(pt, sz):
  return (pt[0] >= 0) and (pt[0] < sz[0]) and (pt[1] >= 0) and (pt[1] < sz[1])

#######################################################
# Iteration and Centroid tools (Floyd's algorithm)

##
# Precompute: precomputes density integration
##

def findCentroids(A, sz, numGen, rho):
  mX = [0]*numGen
  mY = [0]*numGen
  m  = [0]*numGen
  
  for y in range(sz[1]):
    for x in range(sz[0]):
      k = A[y][x]
      m[k]  +=   rho(x, y)
      mY[k] += y*rho(x, y)
      mX[k] += x*rho(x, y)

  centroids = [(mX[k]/m[k], mY[k]/m[k]) for k in range(numGen) if m[k] > 0]
  return centroids

#######################################################
# Testing framework:

def drawCirc(draw, pt, r):
  pt0 = (pt[0]-r, pt[1]-r)
  pt1 = (pt[0]+r, pt[1]+r)
  draw.ellipse([pt0, pt1], fill=(0))

if (__name__ == '__main__'):
  # Generate some arbitrary Voronoi Diagram
  (xSz, ySz) = (100, 100)
  w = 25
  genPts = []
  """
  for x in itertools.product(range(0, xSz, 10), range(0, ySz, 10)):
    genPts.append(x)
  """
  for i in range(40):
    genPts.append((int(xSz*random.random()), int(ySz*random.random())))
  colors = []
  for i in range(len(genPts)):
    colors.append((int(256*random.random()), int(256*random.random()), int(256*random.random())))

  for i in range(20):
    m = getVoronoi(genPts, (xSz, ySz))

    im = Image.new('RGB', (xSz*w, ySz*w), (255, 255, 255))
    draw = ImageDraw.Draw(im)
    for pt in itertools.product(range(xSz), range(ySz)):
      pt0 = (pt[0]*w, pt[1]*w)
      pt1 = (w*pt[0]+w-1, w*pt[1]+w-1)
      draw.rectangle([pt0, pt1], colors[m[pt[1]][pt[0]]], (0, 0, 0))

    for pt in genPts:
      drawCirc(draw, (pt[0]*w + w/2, pt[1]*w + w/2), 10)
    
    centroids = findCentroids(m, (xSz, ySz), len(genPts), lambda x, y : 1)
    for pt in centroids:
      drawCirc(draw, (pt[0]*w + w/2, pt[1]*w + w/2), 5)

    genPts = [(round(pt[0]), round(pt[1])) for pt in centroids]
    
    im.save('./Img' + str(i) + '.jpg')

  im.show()
