#  The entire blur process:

#   K = (F D F') Image    (NxN) 

img <- matrix(imgvec,96,64)   # reshape back to the matrix shape 

   cutt <- 16
   addside <- matrix(0,96, cutt)
                                   # cut image --- img <- img[(cutt+1):(96-cutt),]
   img <- cbind(addside,img,addside) # augmented image
 
   m <- 5
   m2 <- 2*m + 1

   vim <- img
   N <- dim(img)[1]

   dim(vim) <- c(N^2,1) 
   v1 <- fft(vim, inverse = T)
 
   co <- rep(1,m)
   for (i in 1:m){
     co <- c(co, rep(0,N-m2), rep(1,m2))
   }
   nk <- length(co)     # now just symmetrize the rest

   c1 <- rep(0,N^2)
   c1[2:nk] <- co[2:nk]
   c1[(N^2):(N^2-nk+2)] <- co[2:nk]
 
   d <- fft(c1)
   v2 <- fft(d*v1)

   max(Im(v2))
   d <- Re(d)
   plot(d)  
   d <- d/max(d)    # normalizing

   blim <- Re(v2)
   blim <- blim/mean(blim)    # normalizing

   g <- blim + rnorm(N^2)*0.001    # add noise with st.dev = 1
     gshow <- g
   dim(gshow) <- c(N,N)

   par(mfrow=c(1,2))
   mi(gshow,'+noise')                # blurred image
   mi(img,'original')

   alpha <- 0.05     # it's ok to set alpha near the noise st.dev.
   dinv <- d/(alpha^2 + d^2)

   inv1 <- fft(g)
   h <- fft(dinv*inv1, inverse = T)
   
   debl <- Re(h)
   dim(debl) <- c(N,N) 
   mi(debl,'Deblurred')
   mi(gshow,'Blurred')