#!/usr/bin/python
# -*- coding: utf-8 -*-

# 2018-05-10 MW
# The objective of this small program is to turn states of an assembly of grain coming from 
#   WOO (DEM) into images using the Kalisphera implementation in spam

import numpy
import os, sys, time
import math
import tifffile
import scipy.ndimage

import spam.kalisphera.makeSphere as kali

folder       = "../dem_data/monotonic/"
folderOutput = "monotonic/"

# input
filenames    = ["triax.-1.eps=0.spheres",
                #"triax.0.eps=-0.0025.spheres",
                #"triax.1.eps=-0.005.spheres",
                #"triax.2.eps=-0.01.spheres",
                #"triax.3.eps=-0.02.spheres",
                #"triax.4.eps=-0.03.spheres",
                #"triax.5.eps=-0.04.spheres",
                #"triax.6.eps=-0.05.spheres",
                #"triax.7.eps=-0.06.spheres",
                #"triax.8.eps=-0.07.spheres",
                #"triax.9.eps=-0.08.spheres",
                #"triax.10.eps=-0.09.spheres",
                #"triax.11.eps=-0.1.spheres",
                #"triax.12.eps=-50124.6.spheres",
                ]

for filename in filenames:
    num_lines = sum(1 for line in open(folder+"positions/"+filename+".txt"))

    pixel_size  = 15.e-6 #m/pixel
    gaussian    = 0.8
    std_dev     = 0.03

    # start of the program

    print "\t ------------------------------------------"
    print "\t Program to turn DEM assemblies into images"
    print "\t ------------------------------------------\n"

    # loading the dem file
    box_size_dem = numpy.genfromtxt(folder+"positions/"+filename+".txt",skip_footer=(num_lines-1),comments='%',usecols=(1,2,3))
    centres   = numpy.loadtxt(folder+"positions/"+filename+".txt",skiprows=0,usecols=(4,3,2)) #! x,y,z
    radii     = numpy.loadtxt(folder+"positions/"+filename+".txt",skiprows=0,usecols=(5))

    # get maximum radius to pad our image (periodic boundaries...)
    r_max         = numpy.amax(radii)
    box_size      = box_size_dem + 3*r_max
    centres[:,:]  = centres[:,:] + 1.5*r_max # move the positions to the new center of the image

    print 'new box size', box_size

    # turn the mm measures into pixels
    box_size = numpy.array([math.ceil(box_size[0]/pixel_size), math.ceil(box_size[1]/pixel_size), math.ceil(box_size[2]/pixel_size)])
    box_size = int(numpy.amax(box_size))
    print 'box size in pixels', box_size

    centres = centres/pixel_size
    radii   = radii/pixel_size

    # create the big image
    Box = numpy.zeros( (box_size, box_size, box_size), dtype="<f8")

    start_time = time.time()

    # create the whole assembly
    kali.run(Box, centres, radii)

    print "finished kalisphera calculations after ... ", time.time() - start_time, " s "

    # some voxels at the contact might get a value higher than 1.0 due to the adding of the spheres...
    # reset these values to 1.0
    Box[numpy.where(Box > 1.0)] = 1.0
    Box[numpy.where(Box < 0.0)] = 0.0

    Box = Box * 0.5
    Box = Box + 0.25

    # apply blur
    if gaussian!=0:
        start_time = time.time()
        Box = scipy.ndimage.filters.gaussian_filter(Box,sigma=gaussian)
        print "finished scipy gaussian calculations after ... ", time.time() - start_time, " s "
        

    # apply noise
    if std_dev!=0:
        Box = numpy.random.normal(Box,scale=std_dev)
        print "\nAdding noise with a std deviation of ", std_dev, "\n"

    # convert to 8bit image
    Box = numpy.rint(Box*255)

    tifffile.imsave( folderOutput+"images/"+filename+"_image_gaussian=%02i"%(gaussian*10)+"_noise=%02i.tif"%(std_dev*100) , Box.astype('uint8') )
