### 7.Example.lennon.R
############################## Generat Noisy Image
data(lennon)
sdimage <- sd(as.numeric(lennon))
nlennon <- ncol(lennon); nlevel <- log2(ncol(lennon)); optlevel <- c(3:(nlevel-1)) 

errortype <- "Gaussian" # or
#errortype <- "Mixture"

#set.seed(55)
if(errortype == "Mixture") {
    n1 <- trunc(nlennon^2 * 0.9)
    n2 <- nlennon^2 - n1
    lennonnoise <- lennon + matrix(sample(c(rnorm(n1, 0, sdimage/5), 
                                    rnorm(n2, 0, sdimage*2))), nlennon, nlennon)
}
if(errortype == "Gaussian")        
    lennonnoise <- lennon+matrix(rnorm(nlennon^2, 0, sdimage), nlennon, nlennon)

par(mfrow=c(2,2), oma=c(0,0,0,0), mar=0.1+c(1,1,2.0,1))
image(lennonnoise, axes=F, col=gray(100:0/100), main="Noisy Image") 

############################## Reconstruction by Universal Threshold

lennonwd <- imwd(lennonnoise)
lennonuv <- imwr(threshold(lennonwd, policy="universal"))
image(lennonuv, axes=F, col=gray(100:0/100), main="Universal Threshold")

############################## Reconstruction by EBayes

tmp <- c(8,9,10)
slevel <- NULL
for (j in 1:(diff(range(optlevel))+1))
    slevel <- c(slevel, tmp+4*(j-1))
slevel <- slevel[length(slevel):1]

sdev <- mad(c(lennonwd[[8]], lennonwd[[9]], lennonwd[[10]]))

lennonwdth <- lennonwd
for (j in 1:length(slevel)) 
    lennonwdth[[slevel[j]]] <- ebayesthresh(lennonwdth[[slevel[j]]], sdev=sdev) 
   
lennoneb <- imwr(lennonwdth)
image(lennoneb, axes=F, col=gray(100:0/100), main="EBayes") 
   
############################## Reconstruction by cv.level

lennoncv <- cvwavelet.image(images=lennonnoise, imagewd=lennonwd, 
                     cv.optlevel=optlevel, cv.bsize=c(1,1), cv.kfold=10)$imagecv

image(lennoncv, axes=F, col=gray(100:0/100), main="cv.level")    

######### MSE of Universal Threshold, Ebayes and cv.level
c(mean((lennonuv-lennon)^2), mean((lennoneb-lennon)^2), mean((lennoncv-lennon)^2))