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

'''
Ajout de bruit gaussien dans les donnees. A chaque iteration, un bruit poivre 
et sel est ajoute, puis un lissage gaussien autour de ce point est ajoute.
On fait un nombre d'iteration = 1024*complexity/25 ce qui equivaud
a complexity/25 des points qui recoivent le centre du noyau gaussien.
Il y en a beaucoup moins que le bruit poivre et sel, car la transformation
est plutôt aggressive et touche beaucoup de pixels autour du centre 

La grandeur de la gaussienne ainsi que son ecart type sont definit par complexity 
et par une composante aleatoire normale.

On a 25 % de chances d'effectuer le bruitage

Ce fichier prend pour acquis que les images sont donnees une a la fois
sous forme de numpy.array de 1024 (32 x 32) valeurs entre 0 et 1.

Sylvain Pannetier Lebeuf dans le cadre de IFT6266, hiver 2010

'''

import numpy
import random
import scipy
from scipy import ndimage

class BruitGauss():
    
    def __init__(self):
        self.proportion_bruit=0.1 #Le pourcentage des pixels qui seront bruites
        self.nb_chng=10 #Le nombre de pixels changes. Seulement pour fin de calcul
        self.sigma_gauss=3.0  #L'ecart type du noyau gaussien
        self.grandeur=7 #Largeur de la fenetre gaussienne
        self.effectuer=1   #1=effectuer et 0=rien faire
        
    def get_settings_names(self):
        return ['proportion_bruit','sigma_gauss','grandeur','effectuer']

    def regenerate_parameters(self, complexity):
        self.proportion_bruit = float(complexity)/25 
        self.nb_chng=int(1024*self.proportion_bruit)
        
        if float(complexity) > 0:
            self.sigma_gauss=max(0,numpy.random.normal(complexity*5,complexity))
            self.grandeur=int(min(31,max(1,8*complexity*numpy.random.normal(1,float(complexity)/2))))
            self.effectuer =numpy.random.binomial(1,0.25)    ##### On a 25% de faire un bruit #####
        else:
            self.effectuer = 0
            self.sigma_gauss = 0
            self.grandeur=1
        #Un peu de paranoia ici, mais on ne sait jamais
        
        if self.grandeur%2 == 0:
            self.grandeur+=1    #Toujours un nombre impair, plus simple plus tard
        return self._get_current_parameters()

    def _get_current_parameters(self):
        return [self.proportion_bruit,self.sigma_gauss,self.grandeur,self.effectuer]

    
    def transform_image(self, image):
        if self.effectuer == 0:
            return image
        image=image.reshape((32,32))

        #creation du noyau gaussien
        gauss=numpy.zeros((self.grandeur,self.grandeur))
        x0 = y0 = self.grandeur/2
        for i in xrange(0,self.grandeur):
            for j in xrange(0,self.grandeur):
                gauss[i,j]=numpy.exp(-4*numpy.log(2) * ((i-x0)**2 + (j-y0)**2) / self.sigma_gauss**2)
        #pylab.contour(gauss)
        #pylab.show()   #Pour voir si la gaussienne est bien comme desiree
        
        #Chaque tour dans la boucle ajoute un pointpoivre et sel, puis
        #y ajoute un bruit gaussien autour afin d'avoir de la correlation dans
        #les points
        
        for i in xrange(0,self.nb_chng):
            x_bruit=int(numpy.random.randint(0,32))
            y_bruit=int(numpy.random.randint(0,32))
            #l'overlap n'est pas dramatique. Pas besoin de mecanisme pour s'en proteger
            image[x_bruit,y_bruit]=max(0,min(1,numpy.random.normal(0.4,self.proportion_bruit*20)))
            
            bord = int((self.grandeur-1)/2)
            #Faire le "smooting"
            for x in xrange(0,self.grandeur):
                for y in xrange(0,self.grandeur):
                    #pour etre certain de ne pas changer le vide
                    if x_bruit-bord+x < 0:
                        continue
                    if y_bruit-bord+y < 0:
                        continue
                    if x_bruit-bord+x > 31:
                        continue
                    if y_bruit-bord+y > 31:
                        continue
                    image[x_bruit-bord+x,y_bruit-bord+y]=max(image[x_bruit-bord+x,y_bruit-bord+y],gauss[x,y]*image[x_bruit,y_bruit])
                    #image[x_bruit-bord+x,y_bruit-bord+y]=min(1,image[x_bruit-bord+x,y_bruit-bord+y]*(1+gauss[x,y]))
                    #Cette derniere ligne n'est pas très interessante. Elle ajoute le bruit
                    #plutot que de prendre le max entre la valeur presente et le bruit. Ca rend l'image un peu 
                    #chaostique, pas une bonne idee
                    
        return image

#---TESTS---

def _load_image():
    f = open('/home/sylvain/Dropbox/Msc/IFT6266/donnees/lower_test_data.ft')  #Le jeu de donnees est en local. 
    d = ft.read(f)
    w=numpy.asarray(d[0])
    return (w/255.0).astype('float')

def _test(complexite):
    img=_load_image()
    transfo = BruitGauss()
    pylab.imshow(img.reshape((32,32)))
    pylab.show()
    print transfo.get_settings_names()
    print transfo.regenerate_parameters(complexite)
    
    img_trans=transfo.transform_image(img)
    
    pylab.imshow(img_trans.reshape((32,32)))
    pylab.show()
    

if __name__ == '__main__':
    from pylearn.io import filetensor as ft
    import pylab
    _test(0.5)