#
#---------------------------------------------------
#
# Simulainvasion v.0.1 - 1 / 08 / 2013
# Authors:
# Alessio Papini - Department of Plant Biology University of Florence Italy, Via La Pira, 4 Firenze, mail alpapiniATunifi.it
# Ugo Santosuosso - Department of Anatomny, Istology and Forensic Medicine, Largo Brambilla, 1 Firenze, mail ugoATunifi.it
# 
# Implements:
# generates a set of points on a map representing a simulation of the spreading of an invasive species
# a second version generates a spreading from a central point to 4 (or else) secondary points and, from these secondary
# points a short distance spreading: it should be useful to test the efficiency of the kmeans/Silhouette methods
# input is a starting bitmap map; a starting point (origin of the invasion); output: akmeans list of coordinates corresponding to
# simulated invasion location, in csv format. On these last a kmeans approach can be applied
# The version 0_3 uses as starting point for clusters ALSO the general origin of the invasion
# Version 0_31 uses the square measure of the euclidean distance for the distance decay function, rather than the simple distance
# moreover it shows the results on a graphic and reduces the exclusion zone
#version 0_34 implements a control to avoid that points are too close to image boundaries
#---------------------------------------------------
#
# importo le librerie 
#
import os, sys
import math
import numpy
import csv
from pylab import plot,show
from numpy import vstack,array
from numpy.random import rand
from scipy.cluster.vq import kmeans,vq,whiten
from numpy import *
import Image
import ImageDraw
import ImageFont
import random
#
# -----------------------------------------------------
#
print "Execute the program as follows: python Simulainvasion0_31.py namemap.bmp originofinvasioncoordx originofinvasioncoordy numberofclustersofsecondaryspreads"
#NUM_ARGS=3
imagemap=sys.argv[1]
xcoordstart=int(sys.argv[2])
ycoordstart=int(sys.argv[3])
clusters=int(sys.argv[4])

# i2 - name of the image containing the map where the coordinates have been calculated, 
#  with its path to directory if different from the current working directory
#

#
# input the number of clusters as an integer
#
nclust=int(sys.argv[4])

im2 = Image.open(imagemap)
WIDTH, HEIGHT = im2.size
#image dimension

rossmodist = ((math.sqrt((math.pow(WIDTH,2))+(math.pow(HEIGHT,2))))/12)
#It takes as limit of the distance decay funciton 1/3 of the image diagonal
rossmoexcl = ((math.sqrt((math.pow(WIDTH,2))+(math.pow(HEIGHT, 2))))/200)                   
#It takes as exclusion zone limit 1/100 of the image diagonal
rossmodistlong = ((math.sqrt((math.pow(WIDTH,2))+(math.pow(HEIGHT,2))))/3)
#It takes as limit of the distance decay funciton 1/3 of the image diagonal
rossmoexcllong = ((math.sqrt((math.pow(WIDTH,2))+(math.pow(HEIGHT, 2))))/10)                   
#It takes as exclusion zone limit 1/20 of the image diagonal
listaxlong=[]
listaylong=[]
listax=[]
listay=[]
listaxlong.append(xcoordstart)
listaylong.append(ycoordstart)
listax.append(xcoordstart)
listay.append(ycoordstart)  

samplelong=0
while samplelong < (nclust-1):
#long distance dispersal producing the clusters
# (nclust-1) since the starting point provided by the user is also one of the centroids of the clusters
# rossmoexcl is the size of the exclusion zone while rossmodist is the limit data deriving from the distance decay function: points with dist higher than that are excluded
 rand = random.random()
 m = int(WIDTH*rand)
 rand = random.random()
 n = int(HEIGHT*rand)
 distlong=(math.sqrt((math.pow((xcoordstart-m),2))+(math.pow((ycoordstart-n),2))))
#distance between the random generated point and the user provided starting point
 randdistdeclong = distlong*(random.random())
 if(distlong>rossmoexcllong and randdistdeclong<rossmodistlong and m<WIDTH-5 and n<HEIGHT-5):
#distance (decay distance function) of the new point with respect to the origin moltiplicated for a random number between 0 and 1
  listaxlong.append(m)
  listaylong.append(n)
  listax.append(m)
  listay.append(n)   
  samplelong=samplelong+1


salvacent=open("simulcentroids.csv","w")
for k in range (len(listaxlong)):
 salvacent.write((str(listaxlong[k]))+","+(str(listaylong[k]))+"\n")
#
salvacent.close()

for l in range (nclust):                     
 sample = 0
 while sample < (int(100/nclust)):

#100 since the program produces only 100 data elements and these will be divided for each cluster
#an alternative may be to have clusters with random number of elements

# rossmoexcl is the size of the exclusion zone while rossmodist is the limit data deriving from the distance decay function: points with dist higher than that are excluded
  rand = random.random()
  x = int(WIDTH*rand)
  rand = random.random()
  y = int(HEIGHT*rand)
  dist=(math.sqrt((math.pow((listaxlong[l]-x),2))+(math.pow((listaylong[l]-y),2))))
#distance between the random generated point and the user provided starting point

  squadist=((math.pow((listaxlong[l]-x),2))+(math.pow((listaylong[l]-y),2)))
#the above measure is the square distance of the random point from the origin
#in the old version the distance decay was proportional to distance

  randdistdec = dist*(random.random())
  if ((dist>rossmoexcl) and (dist<rossmodist)  and m<WIDTH-5 and n<HEIGHT-5):
#distance of the new point with respect to the origin moltiplicated for a random number between 0 and 1
   listax.append(x)
   listay.append(y)   
   sample=sample+1
   print "ok", randdistdec
  else:
   print randdistdec

print "centroidi", listaxlong, listaylong  
print "datisimulati", sample, listax, listay

salva=open("simuldata.csv","w")
for i in range (len(listax)):
 salva.write((str(listax[i]))+","+(str(listay[i]))+"\n")
#
salva.close()

#with open('simuldata.csv', 'wb') as csvfile:
#    spamwriter = csv.writer(csvfile)
#    spamwriter.writerow(listax[1])
#    for i in range (len(listax)):
#     spamwriter.writerow([listax[i] for x in (len(listax))])
#csvfile.close()





# listcoordx=
# listcoordy=.append
#values inserted in the lists

#if the two conditions are not respected it goes back with the cycle (i don't know if it works)
#------------------------output an image containing the simulated points----------------------------------------- 
print" Drawing points.. ",
draw = ImageDraw.Draw(im2)
for cross in range (len(listax)):
  draw.ellipse((listax[cross]-1, listay[cross]-1, listax[cross]+1, listay[cross]+1), fill=(0,0,0))
#  testo="n."+str(int(cross+1))+" ("+str(weight[cross])+"%)"
#  lung_testo=len(testo)
#  draw.text((listax[cross]-int(5*lung_testo/2), listay[cross]-13),testo,fill=(255,255,0))

# with the above instruction it says that 255-result is not an integer but a double and i agree
# may be i should normalize tha value but in the c instruction ai don't see anything like that
# im4.putpixel((i,j),(255, 255, 0)) if i use this instruction i obtain an image that is completely yellow
#
print" Saving image and deleting temporary files"
im2.save("simulpoints.bmp")