# Brouillon

#print(img1[57,42,])
#igne_data=c( 0.2126, 0.2126, 0.2126, 0.7152,0.7152,0.7152,0.0722,0.0722,0.0722)
#A=matrix(data=ligne_data,nrow=3,ncol=3)
#B=A[1,]

#Chargement de l'image
img=readPNG("~/Dropbox/sylvain-rubenthaler/enseignement/m1im/monte-carlo/shot-noise.png")
width<-ncol(img)
height<-nrow(img)
#On transforme l'image en une image noir & blanc
img3<-matrix(rep(0,width*height),height,width)
for (i in 1:height)
{
	for (j in 1:width)
	{
		pixel=img[i,j,]	
		img3[i,j]=floor(min(1,4*pixel[1]))
	}
}
display(img3)


#-------------------------------------------------------------------------------

library("png")
img1=readPNG("~/Dropbox/sylvain-rubenthaler/enseignement/m1im/monte-carlo/simpsons_PNG88.png")

 # On transorme l'image en image niveaux de gris
 
grayscale<-function(src)
{
	width<-ncol(src)
	height<-nrow(src)
	img<-matrix(rep(0,width*height),height,width)
	for (i in 1:height)
	{
		for (j in 1:width)
		{
			pixel=src[i,j,]	
			X=matrix(data=pixel[1:3],nrow=3,ncol=1)
			img[i,j]=B%*%X
		}
	}
	return(img)
}

img2=grayscale(img1)

width<-ncol(img2)
height<-nrow(img2)

#On transforme l'image en une image noir & blanc

img3<-matrix(rep(0,width*height),height,width)
for (i in 1:height)
{
	for (j in 1:width)
	{
		pixel=img1[i,j,]	
		img3[i,j]=floor(min(1,4*pixel[1]))
	}
}
display(img3)

# Bruitage de l'image

seuil=0.1
img4=img3
for (i in 1:height)
{
	for (j in 1:width)
	{
		u=runif(1,0,1)
		if (u<seuil)
		{ img4[i,j]=1-img4[i,j] }
	}
}
display(img4)



#Calcul de la fonction H (on oublie les aretes pres du bord)
H<-function(x)
{
	s=0
	xp=2*x-1 # on veut obtenir une matrice de 1 et -1
	for (i in 2:(height-1))
	{
		for (j in 1:(width/2-1))
		{
			k=2*j+i%%2
			s=s+(xp[i,k]-xp[i-1,k])^2+(xp[i,k]-xp[i,k-1])^2
			s=s+(xp[i,k]-x[i+1,k])^2+(xp[i,k]-xp[i,k+1])^2
		}
	}
	return(0.5*s)
}

a<-function(x,y)
{
	s=0
	for (i in 1:height)
	{
		for (j in 1:width)
		{
			if(x[i,j]==y[i,j])
			{ s=s+1 }
		}
	}
	return(s)
}

d<-function(x,y)
{
	s=0
	for (i in 1:height)
	{
		for (j in 1:width)
		{
			if(abs(x[i,j]-y[i,j])>0)
			{ s=s+1 }
		}
	}
	return(s)
}

# Calcul de la densite
mup<-function(x)
{
	return(exp(-beta*H(x))*p^(a(x,img3))*(1-p)^(d(x,img3)))
}

# Rapport des densités
rapport<-function(y,x)
{
	return(exp(-beta*H(y)+beta*H(x)+log(p)*(a(y,img3)-a(x,img3))+log(1-p)*(d(y,img3)-d(x,img3))))
}
	
# transition +
tplus<-function(x)
{
	y=2*x-1
	z=x
	for (i in 2:(height-1))
	{
		for (j in 1:((width-2)/2))
		{
			k=2*j+i%%2
			M=y[i-1,k]+y[i,k-1]+y[i+1,k]+y[i,k+1]
			if (img3[i,k]==1)
			{ pp=p}
			else
			{ pp=1-p}
			q=exp(2*beta*M)*pp/(exp(2*beta*M)*pp+exp(-2*beta*M)*(1-pp))
			z[i,k]=rbinom(1,1,q)
		}
	}
	return(z)
}

	
	
# transition -
tmoins<-function(x)
{
	y=2*x-1
	z=x
	for (i in 2:(height-1))
	{
		for (j in 1:((width-2)/2))
		{
			k=2*j+1-i%%2
		  M=y[i-1,k]+y[i,k-1]+y[i+1,k]+y[i,k+1]
			if (img3[i,k]==1)
			{ pp=p}
			else
			{ pp=1-p}
			q=exp(2*beta*M)*pp/(exp(2*beta*M)*pp+exp(-2*beta*M)*(1-pp))
			z[i,k]=rbinom(1,1,q)
		}
	}
	return(z)
}

# transition
transit<-function(x)
{
	u=runif(1,0,1)
	if (u<0.5)
	{ z=tplus(x) }
	else
	{ z=tmoins(x) }
	return(z)
}

#Constantes
beta=4
p=0.9


# Chaine de Metropolis

x=img4
width=ncol(img4)
height=nrow(img4)
#png(filename="~/Dropbox/sylvain-rubenthaler/enseignement/m1im/monte-carlo/simpsons-bruitee.png")
display(x)
#dev.off()
n=100
for (i in 1:n)
{
	y=transit(x)
	#display(y)
	print("coucou")
	alpha=min(1,rapport(y,x))
	print(alpha)
	u=runif(1,0,1)
	if (u<alpha)
	{
		print("accepte") 
		x=y 
	}
	display(x) 
}
	