xy<-rbind(rep(0:255,256),rep(0:255,rep(256,256)))
indw<-c(1:12,29:48,73:100,133:168,209:256)
w0<-matrix(rep(0,256*256),ncol=256)
w0[indw,]<-1
w0[,indw]<-!w0[,indw]
w0<-w0-.5
w0[((xy[1,]-129)^2+(xy[2,]-129)^2)<=10000&((xy[1,]-129)^2+(xy[2,]-129)^2)>=4900]<- 0
w0[abs(xy[1,]-xy[2,])<=20&((xy[1,]-129)^2+(xy[2,]-129)^2)<4900]<- 0
w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625]<- 0
w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[2,]>27&xy[2,]<31]<- -.5
w0[((xy[1,]-225)^2+2*(xy[2,]-30)^2)-(xy[1,]-225)*(xy[2,]-30)<=625&xy[1,]>223&xy[1,]<227]<- .5
w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625]<- 0
w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[1,]>27&xy[1,]<31]<- -.5
w0[((xy[2,]-225)^2+2*(xy[1,]-30)^2)+(xy[2,]-225)*(xy[1,]-30)<=625&xy[2,]>223&xy[2,]<227]<- .5
w0[((xy[2,]-225)^2+(xy[1,]-225)^2)+1*(xy[2,]-225)*(xy[1,]-225)<=400]<- 0
w0[((xy[2,]-30)^2+(xy[1,]-30)^2)<=256]<-0
sigma<-.25
y<-w0+rnorm(w0,0,sigma)
# increase rmax for better results
yhat<-awsbi(y,rmax=3)
par(mfrow=c(1,3))
image(y,col=gray((0:255)/255))
title("Noisy image")
image(yhat$yhat,zlim=range(y),col=gray((0:255)/255))
title("AWS reconstruction")
image(w0,zlim=range(y),col=gray((0:255)/255))
title("Original image")
rm(y,w0,yhat,xy)
Run the code above in your browser using DataLab