I am trying to implement a low pass Wiener filter using the mean-
squared method, to provide elimination of noisy areas, background
texture smoothing, and contrast enhancement. The algorithm I'm using
is (to produce output image Q from input image P), to evaluate a 3x3
neighbourhood for each pixel, and

  Q(x,y) =  mu + [  (sigmasquared - nusquared)/sigmasquared  ]  ( P
(x,y) - mu )

Where:
  mu  is the local mean across the neighbourhood
  sigmasquared is the variance of the 3x3 neighbourhood
  nusquared is the estimated local variance

Now, I don't know the source of the putative high frequency noise, so
I have to estimate nusquared, which is where I am running into
problems. All the code I can find on the net estimates nusquared by
using the variance of the whole of the original image. Not a canonical
source I know, but see, e.g., the bottom of:
  http://www.vislab.uq.edu.au/education/sc3/2000/xinghe/project/theory.html

Now this doesn't work for me, as in images to which noise has NOT been
added, the local variance is consistently SMALLER than the image
variance (in one test image it is never larger than image variance.
This seems to me bound to be the case on (for example) an image that
has been enlarged by bicubic interpolation and thus is already
"smooth". The problem here is that the Wiener filter does terrible
things to the image because the term in square brackets goes
substantially negative (as sigmasquared << nusquared), causing the
output to wildly oscillate.

I think the problem here is that all the code examples I have found
have assumed that the image variance is a good estimate for the local
variance. That's true if you are doing a compsci project where the
input image has had a lot of noise added to it, but isn't universally
true. I think the problem is that I need a better estimate for
nusquared.

What would you recommend here? Taking the mean of the 3x3 variances
across the image?

Alex