#
### version 1c is counterpart for Deblur1c, to test the influence of noise
#

myim2 <- function(X){
 nrow <- dim(X)[1]
 image(t(X[nrow:1,]), col=gray(seq(0,1,0.01)), xaxt="n", yaxt="n")
}

myim <- function(X){
 nrow <- dim(X)[1]
 image(t(X[nrow:1,]), col=gray(seq(1,0,-0.01)), xaxt="n", yaxt="n")
}


regty <- function(H){   # regularuty of the array H 
  n1 <- dim(H)[1]
  n2 <- dim(H)[2]
  reg <- sum((H[1:(n1-1),] - H[2:n1,])^2)
  reg <- reg + sum((H[,1:(n2-1)] - H[,2:n2])^2)
  sqrt(reg/(n2*(n1-1) + n1*(n2-1)))
} 
 

nrow <- 96;
ncol <- 64;
N <- nrow*ncol;
i <- 1
 inpName <- sprintf("c:/Oleg/deblur/asc%d.csv",i)
 A <- as.matrix(read.csv(inpName, header=FALSE))
  A1 <- A[,1:ncol]
  A2 <- A[,(ncol+1):(2*ncol)]
  dim(A1) <- N
  dim(A2) <- N
 imgvec <- A1

 img <- matrix(imgvec,96,64)   # reshape back to the matrix shape 
 img <- max(img)-img
 Ni <- dim(img)[1]
 Nj <- dim(img)[2]

m <- 2     # blur order
m2 <- 2*m+1 
B <- matrix(1/m2^2, m2,m2)      # uniform blur
 

  blim2 <- 0*img
  for (i in (m+1):(Ni-m))
   for(j in (m+1):(Nj-m)){
     blim2[(i-m):(i+m), (j-m):(j+m)] <-  blim2[(i-m):(i+m), (j-m):(j+m)] + img[i,j] * B
  }

   par(mfrow=c(1,3))
     myim(img)
     myim(blim2)

 # G <- blim2[(1+m):(Ni-m), (1+m):(Nj-m)]           # this is the source for lkhd

 ##+++++++++++++++++++++++++++  START HERE if RUNNING IN COMPARISON MODE
  
  G <- gshow
  
  F <- G[(1+m):(Ni-m), (1+m):(Nj-m)]  # will be the restored image, edges removed
  NiF <-  Ni- 2*m 
  NjF <-  Nj - 2*m 
  

  Gf <- G*0     # this is the smoothing of F, to be iterated to perfection
  for (i in 1:NiF)
   for(j in 1:NjF){
     Gf[i:(i+2*m), j:(j+2*m)] <-  Gf[i:(i+2*m), j:(j+2*m)] + F[i,j] * B
  }

  
  regty(F)   # actually it's estimate of irregularity, that is sigma
  s <- 0.1    #  a bit higher than noise in F 
  se <- 0.04  # some "small" number, also tunable 
  
   # 5 and 1 give an ok result for m = 2

  Mask <- 0*G
  Mask[(1+2*m):(Ni-4*m), (1+2*m):(Nj-4*m)] <- 1    # to avoid edge effects       	

  lkhdif <- function(i,j,fold,fnew){   # gives loglkhd difference, to be used in Metropolis 
     mask <- Mask[i:(i+2*m),j:(j+2*m)] 
     Gw <- G[i:(i+2*m),j:(j+2*m)]   # affected window
     Gfw <- Gf[i:(i+2*m),j:(j+2*m)]
     Gwn <- Gfw + (fnew-fold)*B      # this is how the new window will look like
   
     loglik <- sum(mask*(Gfw  - Gw)^2)  - sum(mask*(Gwn  - Gw)^2)    
                                           # changes in lkhd from change in f
 
     nhb <- c(F[i+1,j], F[i-1,j], F[i,j+1], F[i,j-1])
     logprior <- sum((fold-nhb)^2) - sum((fnew-nhb)^2)     # changes in regularity from change in f
     (loglik/(se^2) + logprior/(s^2))/2
   }

  MC <- 100
  nskip <- 5000
  delf <- 0.15
  Nacc <- 0
lkhds <- rep(0,MC)
 lkhd <- 0
 meanF <- F*0

  for(mc in 1:MC){
   for (skip in 1:nskip){
    i <- sample(2:(NiF-1),1)    # 1-margin taken off to avoid edge effects
    j <- sample(2:(NjF-1),1)
    fold <- F[i,j]
    fnew <- fold + rnorm(1,0,delf)  
    ldif <- lkhdif(i,j,fold,fnew)  

    accprob <- exp(min(ldif,0))
    u <- runif(1)
    if (u < accprob){                  # Metropolis step
     F[i,j] <- fnew
      Nacc <- Nacc + 1
      lkhd <- lkhd + ldif
      Gf[i:(i+2*m),j:(j+2*m)] <-  Gf[i:(i+2*m),j:(j+2*m)] + (fnew-fold)*B  
    }
  } # end skip
  # output the new pic
   mi(F, paste('iter ',mc))
   lkhds[mc] <- lkhd
   meanF <- meanF + F/MC

     # fitting s 
  dfs <- (NiF*(NjF-1) + NjF*(NiF-1))
  ssq <- regty(F)^2*dfs
  s <- sqrt((ssq + 0.001)/rchisq(1,dfs+5))    # includes prior var and prior df


 }
 (acc.pct <- Nacc/(MC*nskip))
plot(lkhds)
Fprep <- pmin(max(G) + 0*F,pmax(meanF,0*F))
mi(Fprep,'Mean posterior')



sqrt(sum(((Gf-G)*Mask)^2)/sum(Mask))

 11.17
  3.24 (!)

 3.0293  (for m=1, with se = 1 and s = 10)
 


  # resampling B

  # F <-  img[(1+2*m):(Ni-2*m), (1+2*m):(Nj-2*m)]    # cheating

  

 
   










  par(mfrow=c(1,3))
     myim2(blim)
     myim2(debl)
     myim2(blim2)