FUNCTION coherence, img

;;    Returns 'coherence' matrix from 2D input matrix.
;;    Each element in the result is the Std. Dev. of the
;;    nine-pixel neighborhood from the original image.

;;    Dick Jackson, Fanning Software Consulting
;;    djackson@dfanning.com

size = Size(img, /Dimensions)
nx = size[0]
ny = size[1]
IF nx LT 3 OR ny LT 3 THEN $
   Message, 'Input image must have at least 3 rows and 3 columns'

img2 = img ^ 2.0
sx = FltArr(nx, ny)
sx2 = sx

FOR dy=0, 2 DO FOR dx=0, 2 DO BEGIN
   sx = sx + img[dx:nx-3+dx, dy:ny-3+dy]   
   sx2 = sx2 + img2[dx:nx-3+dx, dy:ny-3+dy]
ENDFOR

coh = FltArr(nx, ny)
coh[1:nx-2, 1:ny-2] = Sqrt((sx2 - (sx^2 / 9)) / 8)

Return, coh

END ;; coherence