# >>>>>>>> implementation of Gibbs sampler with binary flip errors
#

setwd('u:/public_html/530/Class17')
source("MCout.R")
colrs <- c("red", "white")

## F <- t(as.matrix(read.csv("image.txt",header=FALSE)))      # has to be a square image

F <- read.csv("im1.txt",header=FALSE)
F <- t(as.matrix(F))
dim(F) <- c(40,40)

S <- F

image(F, main="true image")

n <- dim(F)[1]
p <- 0.25        # proportion of pixels flipped
flip <- sample(c(0,1), n^2, replace=TRUE, prob=c(1-p,p))    # randomly chose sites for flipping
dim(flip) <- c(n,n)

G <- F * (1-flip) + (1-F)*flip   # degraded image

G <- rbind(rep(0,n), G, rep(0,n))
G <- cbind(rep(0,n+2), G, rep(0,n+2))  # padding with an extra frame of 0's 

n <- n+1
n1 <- n+1
X <- sample(c(0,1), (n1)^2, replace=TRUE, prob = c(0.5,0.5) )                           # start with random initial configuration
dim(X) <- c(n1,n1)
X[,1] <- rep(0,n1);  X[,n1] <- rep(0,n1); X[1,] <- rep(0,n1);  X[n1,] <- rep(0,n1); 

image(1:n1,1:n1,G, main="noisy image")
image(1:n1,1:n1,X)

C <- 0*G               # count 1-updates
T <- 0*G               # total number of updates


M <- 400000           # 25 samples per site; approx. 1000 iter per sec.
skip <- 1
k <- 0.8

for (m in 1:M){
    i <- sample((2:n),1)
    j <- sample((2:n),1)   ## random site updating

    SV <- (X[i+1,j] + X[i-1,j] + X[i,j+1] + X[i,j-1])               # sum of values of X at neighbor sites 
    mult <- ifelse(G[i,j] == 0, 1-p, p)          # multiplier in the likhd
    expkSV <- exp(k*(SV-2)*2)
    p0 <- mult/expkSV                                        
    p1 <- (1-mult)*expkSV
    newcol <- sample(c(0,1),1, prob=c(p0,p1) ) 
    rect(i-0.5, j-0.5, i+0.5, j+0.5, col = colrs[newcol+1], border= NA)  # rendering
    X[i,j] <- newcol
    C[i,j] <- C[i,j] + newcol
    T[i,j] <- T[i,j] + 1
}

A <- C/T

par(mfrow=c(2,2))

image(S, main="Original Signal")
image(1:n1,1:n1, A, main="Average GS output")        
  

Fest <- A[2:n,2:n] < 0.5      # estimated image
mean(abs(Fest-F))             # proportion of pixels estimated correctly
 


image(Fest, main = "Estimated image")

image(Fest+F-1, main = "Difference")










par(mfrow=c(2,2), mar=rep(3,4))
  cutoff <- c(0.3,0.5,0.8)

  for (i in 1:3){
     image(1:n1,1:n1,A > cutoff[i])
  }
 image(1:n1,1:n1,G)