% weiner filter problem
rand('state',0);randn('state',0);
nx = 64; ny = 64;

h1 = ones([9 9]);

p=0.002;
ss_tmp = (rand(nx,ny) < p) - p;
ss = conv2(ss_tmp,h1,'same')/sqrt(p*(1-p));

yy = ss + 5*randn(size(ss));  % what is the noise variance?

%
% determine and display autocorrelation functions here
%

% zero pad to prevent circular filtering
px = 2*nx; py = 2*ny;

% build Wiener filter on higher-res (px by py) array
wx = [-px/2:px/2-1]*2*pi/px;
wy = [-py/2:py/2-1]*2*pi/py;

Hw = 1;  % this is a place holder for the true Wiener filter

% filter in Fourier domain
Y = fftshift(fft2(yy,px,py));
shatpad = real(ifft2(fftshift(Y.*Hw)));
shat = shatpad(1:nx,1:ny);

shatad = anisodiff(yy,4,10,.2,1);

msew = mean(mean((shat-ss).^2))
msead = mean(mean((shatad-ss).^2))

%diplay filtered results
clim = [min(min(yy)) max(max(yy))];
subplot(221); imagesc(ss,clim); colormap gray; colorbar; title('Original Image')
subplot(222); imagesc(yy,clim); colormap gray; colorbar; title('Original Image + Noise')
subplot(223); imagesc(shat,clim); colormap gray; colorbar; title('Wiener Filtered Image')
subplot(224); imagesc(shatad,clim); colormap gray; colorbar; title('Adaptive Filtered Image')