;+
; Computes the local moments for a given window size for an array.
; 
; @returns the local mean
; @param image {in}{required}{type=2D numeric array} original image; will be 
;        converted to float or double
; @param width {in}{required} size of window
; @keyword double {in}{optional}{type=boolean} set to do computations as doubles
; @keyword sdev {out}{optional}{type=2D float/double array} local standard 
;          deviation
; @keyword variance {out}{optional}{type=2D float/double array} local variance
; @keyword skewness {out}{optional}{type=2D float/double array} local skewness
; @keyword kurtosis {out}{optional}{type=2D float/double array} local kurtosis
; @keyword edge_truncate {in}{optional}{type=boolean} compute edge values by 
;          repeating
; @keyword edge_wrap {in}{optional}{type=boolean} compute edge values by padding
;          array with zeros
; @keyword edge_zero {in}{optional}{type=boolean} compute edge values by 
;          wrapping
; @requires IDL 6.2
;-
function mg_local_moment, image, width, double=double, $
                          sdev=sdev, variance=variance, skewness=skewness, $ 
                          kurtosis=kurtosis, edge_truncate=edge_truncate, $
                          edge_wrap=edge_wrap, edge_zero=edge_zero
    compile_opt strictarr
    on_error, 2

    ; sanity checking on arguments
    if (n_params() ne 2) then message, 'incorrect number of arguments' 
    if (size(image, /n_dimensions) ne 2) then begin
        message, 'IMAGE parameter must be a 2D array'
    endif
    
    im = keyword_set(double) ? double(image) : float(image) 
    kernel = keyword_set(double) $
             ? dblarr(width, width) + 1.0D $
             : fltarr(width, width) + 1.0
    n = width^2

    ; mean
    local_mean = convol(im, kernel, edge_truncate=edge_truncate, $
                        edge_wrap=edge_wrap, edge_zero=edge_zero) / n

    ; variance or standard deviation
    if (arg_present(sdev) $
        || arg_present(variance) $
        || arg_present(skewness) $
        || arg_present(kurtosis)) then begin
        squared = convol(im^2, kernel, edge_truncate=edge_truncate, $
                         edge_wrap=edge_wrap, edge_zero=edge_zero)
        variance = (squared  - n * local_mean^2) / (n - 1.0)
        sdev = sqrt(variance)
    endif

    ; skewness
    if (arg_present(skewness) || arg_present(kurtosis)) then begin
        cubed = convol(im^3, kernel, edge_truncate=edge_truncate, $
                       edge_wrap=edge_wrap, edge_zero=edge_zero)  
        skewness = (cubed / n - 3.0 * local_mean * squared / n + 2.0 * local_mean^3) / sdev ^ 3
    endif

    ; kurtosis
    if (arg_present(kurtosis)) then begin
        fourth = convol(im^4, kernel, edge_truncate=edge_truncate, $
                        edge_wrap=edge_wrap, edge_zero=edge_zero)
        kurtosis = (fourth / n - 4.0 * local_mean * cubed / n + 6.0 * local_mean^2 * squared / n - 3.0 * local_mean^4) / variance^2 - 3.0
    endif
    
    return, local_mean
end