FUNCTION coherence, array, width

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

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

size = Size(array, /Dimensions)
IF N_Elements(size) NE 2 THEN $
   Message, 'Input array must be 2-dimensional.'
nx = size[0]
ny = size[1]
IF nx LT 3 OR ny LT 3 THEN $
   Message, 'Input array must have at least 3 rows and 3 columns.'

IF N_Elements(width) EQ 0 THEN width = 3

;    Make Sum-of-X array (sx) where each element contains the sum of
;    values from the neighborhood of the original array.

kernel = Replicate(1, width, width)     ; kernel for convolution to simply
                                        ; sum all values with equal weight
sx = Convol(Float(array), kernel)

;    Make Sum-of-X-squared array (sx2) where each element contains the sum of
;    *squared* values from the neighborhood of the original array.

sx2 = Convol(array ^ 2.0, kernel)

;    Calculate, for each pixel, the standard deviation of the set of
;    the neighborhood values:
;       SD = Sqrt((Sum(X^2) - (Sum(X)^2 / n)) / (n-1))

n = width ^ 2.0

Return, Sqrt((sx2 - (sx^2 / n)) / (n-1))

END ;; coherence