% weiner filter problem

rand('state',0);randn('state',0);
nx = 16; ny = 16;
[xx yy] = ndgrid([-nx/2:nx/2-1],[-ny/2:ny/2-1]);
h1 = ones([3 3]);

p=0.015;
ff_tmp = (rand(nx,ny) < p) - p;
ff = conv2(ff_tmp,h1,'same')/sqrt(p*(1-p));
gg = ff + 2*randn(size(ff));  % what is the noise variance?
g1 = gg(:);  % make image a vector

% build arrays to calculate covariance matrix
x1 = xx(:);  
[x1a x1b] = ndgrid(x1);
% the following array contains the difference between 
% the x-coordinates for the covariance elements
xinddiff = (x1a - x1b);  

y1 = yy(:);
[y1a y1b] = ndgrid(y1);
% the following array contains the difference between 
% the x-coordinates for the covariance elements
yinddiff = (y1a - y1b);

% now calculate the covariance matrices 

% calc estimate of f
fhat = ????*g1;
fhat2 = reshape(fhat,size(ff));

mseor = mean(mean((gg-ff).^2))
msew = mean(mean((fhat2-ff).^2))

%diplay filtered results
clim = [min(min(gg)) max(max(gg))];
subplot(221); imagesc(ff,clim); colormap gray; colorbar; title('Original Image')
subplot(222); imagesc(gg,clim); colormap gray; colorbar; title('Original Image + Noise')
subplot(223); imagesc(fhat2,clim); colormap gray; colorbar; title('Wiener Filtered Image')