local standard deviation of an image

Just checked the file SIGMA_FILTER at the NASA site. I really must spend more time browsing these great sites.

The code is similar, however it does not calculate the true variance under the mask. They calculate for a box width of n, (ignoring centre pixel removal):

   mean_im=(smooth(image, n) )
   dev_im = (image - mean_im)^2
   var_im = smooth(dev_im, n)/(n-1)

This is not the true variance of the pixels under the box mask, as each pixel in the mask is having a different mean subtracted. For example (read this as a formula if you can!):

   Pseudo_Variance = SUM ij ( ( I(x+i,y+j) - MEAN(x+i,y+j)^2) /(n-1)

instead of true variance:

   Variance = SUM ij ( ( I(x+i,y+j) - MEANxy)^2) /(n-1)

which can be reduced to:

   {(SUM ij ( ( I(x+i,y+j)^2 ) - (SUM ij I(x+i,y+j) ) ^2)/n }/(n-1)

hence the non-loop method we use below:

     ; calc box mean
   mean_im = smooth(image, n)
     ; calc box mean of squares
   msq_im = smooth(image^2, n)
     ; hence variance
    var_im = ( msq_im - mean_im^2) * (n/(n-1.0))