def bg_thresholds( dark_arr, n_std=3 ):
    """
    Calculate band-wise mean radiance plus 3 standard deviations for pixels in
    `dark_arr`. Lyzenga et al. 2006 says:

    "...the blue and green bands are thresholded at the deep-water mean
    radiance plus three standard deviations."

    This method will calculate the mean + 3 std for all bands. You'll have to
    pick out the blue and green ones later if that's what you're after.

    Parameters
    ----------
    dark_arr : numpy Masked Array
        Typically, this will be the output of `Lyzenga2006.dark_pixels_array`.
    n_std : int (Default value = 3)
        The number of standard deviations to add to the mean. Lyzenga et al.
        2006 uses 3 so that's the default.

    Returns
    -------
    numpy array
        A 1D array with as many elements as there are bands in `dark_arr`. Each
        element corresponds to the threshold for its respective band.
    """
    nbands = dark_arr.shape[-1]
    darkmeans = dark_arr.reshape(-1,nbands).mean(0).data
    darkstds = dark_arr.reshape(-1,nbands).std(0).data
    return darkmeans + n_std * darkstds