function REGION_GROW, array, roiPixels, $
    ALL_NEIGHBORS=allNeighbors, $
    NAN=nan, $
    STDDEV_MULTIPLIER=stdDevMult, $
    THRESHOLD=threshold

    ulong = 0
    CATCH, iErr
    if (iErr ne 0) then begin
        if ((!ERROR_STATE.CODE eq -512) && (ulong eq 0)) then begin
            ; The array may have more regions than fit within a UINT.
            ; If this is the case, use the ULONG keyword.
            ulong = 1
            MESSAGE, /RESET
            goto, doLabelRegion
        endif else begin
            CATCH, /CANCEL
            ON_ERROR, 2
            MESSAGE, /REISSUE_LAST
            return, -1L
        endelse
    endif

    arrDims = SIZE(array, /DIMENSIONS)

    ; Extract the image pixel values that fall within the input region.
    roiPixelData = array[roiPixels]
    nROIPixels = N_ELEMENTS(roiPixels)

    if (TOTAL(FINITE(roiPixelData)) eq 0) then $
        MESSAGE, $
            'Unable to grow region that consists entirely of missing data.'

    ; Set the initial threshold range for the grown region.
    if (N_ELEMENTS(threshold) eq 2) then begin
        thresh_lo = threshold[0]
        thresh_hi = threshold[1]
    endif else $
        thresh_lo = MIN(roiPixelData, MAX=thresh_hi, NAN=nan)

    ; If requested, compute a threshold range based on the standard
    ; deviation and mean.
    if (N_ELEMENTS(stdDevMult) ne 0) then begin
        if (N_ELEMENTS(threshold) eq 2) then begin
            errStr = 'STDDEV_MULTIPLIER and THRESHOLD are mutually ' + $
                     'exclusive.  Using the THRESHOLD value.'
            MESSAGE, /CONTINUE, errStr
        endif else begin
            if (nROIPixels gt 1) then begin
                mean = MEAN(roiPixelData, /DOUBLE, NAN=nan)
                sdev = STDDEV(roiPixelData, /DOUBLE, NAN=nan)
            endif else begin
                mean = roiPixelData
                sdev = 0.0d
            endelse
            offset = stdDevMult * sdev
            if (SIZE(array, /TYPE) le 3) then begin
                thresh_lo = ROUND(mean - offset)
                thresh_hi = ROUND(mean + offset)
            endif else begin
                thresh_lo = mean - offset
                thresh_hi = mean + offset
            endelse
        endelse
    endif

    ; Threshold the original image.
    if (KEYWORD_SET(nan)) then begin
        ; Temporarily overwrite array NaN values with 0s.
        iNan = WHERE(~FINITE(array), nNan)
        if (nNan gt 0) then begin
            oldNanVals = array[iNan]
            array[iNan] = 0
        endif
    endif else $
        nNan = 0

    threshArray = (array ge thresh_lo and array le thresh_hi)

    ; Replace NaN values in input array, and set result of threshArray
    ; accordingly.
    if (nNan gt 0) then begin
        array[iNan] = oldNanVals
        threshArray[iNan] = 0b
    endif

    ; Label the regions within the thresholded image.
doLabelRegion:
    labelArray = LABEL_REGION(threshArray, ALL_NEIGHBORS=allNeighbors, $
        ULONG=ulong)

    ; Determine which labels fall within the ROI.
    if (nROIPixels gt 1) then begin
        labels = WHERE(HISTOGRAM(labelArray[roiPixels], MIN=0) ne 0, nLabels)
    endif else begin
        nLabels = 1
        labels = labelArray[roiPixels]
    endelse

    ; Compute a total of all labeled pixels for the labels within the ROI
    ; (excluding pixels with a label of 0).
    labelHist = HISTOGRAM(labelArray, REVERSE=r, MIN=0)
    nPixels = TOTAL(labelHist[labels]) - (labelHist[labels[0]] * $
                    (labels[0] eq 0))

    ; Create a new array of pixel indices for the labeled pixels.
    if (nPixels gt 0) then begin
        growROIPixels = LONARR(nPixels)
        j = 0
        for i=0L, nLabels-1 do begin
            if (r[labels[i]+1] le N_ELEMENTS(r) and labels[i] ne 0) then begin
                growROIPixels[j] = r[r[labels[i]]:r[labels[i]+1]-1]
                j = j + labelHist[labels[i]]
            endif
        endfor
    endif else $
        growROIPixels = -1L

    return, growROIPixels
end