Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates.

    Parameters:
    -----------
    input :     (1) full HEALPix mask of booleans, or (2) a set of (lon, lat) coordinates for catalog objects that define the occupied pixels.
    nside_pix : nside_pix is meant to be at coarser resolution than the input mask or catalog object positions
    so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. 

    Returns:
    --------
    lon,lat,area : Return the longitude and latitude of the random positions (deg) and the total area (deg^2).

def randomPositions(input, nside_pix, n=1):
    """
    Generate n random positions within a full HEALPix mask of booleans, or a set of (lon, lat) coordinates.

    Parameters:
    -----------
    input :     (1) full HEALPix mask of booleans, or (2) a set of (lon, lat) coordinates for catalog objects that define the occupied pixels.
    nside_pix : nside_pix is meant to be at coarser resolution than the input mask or catalog object positions
    so that gaps from star holes, bleed trails, cosmic rays, etc. are filled in. 

    Returns:
    --------
    lon,lat,area : Return the longitude and latitude of the random positions (deg) and the total area (deg^2).

    """
    input = np.array(input)
    if len(input.shape) == 1:
        if hp.npix2nside(len(input)) < nside_pix:
            logger.warning('Expected coarser resolution nside_pix in skymap.randomPositions')
        subpix = np.nonzero(input)[0] # All the valid pixels in the mask at the NSIDE for the input mask
        lon, lat = pix2ang(hp.npix2nside(len(input)), subpix)
    elif len(input.shape) == 2:
        lon, lat = input[0], input[1] # All catalog object positions
    else:
        logger.warning('Unexpected input dimensions for skymap.randomPositions')
    pix = surveyPixel(lon, lat, nside_pix)

    # Area with which the random points are thrown
    area = len(pix) * hp.nside2pixarea(nside_pix, degrees=True)

    # Create mask at the coarser resolution
    mask = np.tile(False, hp.nside2npix(nside_pix))
    mask[pix] = True

    # Estimate the number of points that need to be thrown based off
    # coverage fraction of the HEALPix mask
    coverage_fraction = float(np.sum(mask)) / len(mask) 
    n_throw = int(n / coverage_fraction)
        
    lon, lat = [], []
    count = 0
    while len(lon) < n:
        lon_throw = np.random.uniform(0., 360., n_throw)
        lat_throw = np.degrees(np.arcsin(np.random.uniform(-1., 1., n_throw)))

        pix_throw = ugali.utils.healpix.angToPix(nside_pix, lon_throw, lat_throw)
        cut = mask[pix_throw].astype(bool)

        lon = np.append(lon, lon_throw[cut])
        lat = np.append(lat, lat_throw[cut])

        count += 1
        if count > 10:
            raise RuntimeError('Too many loops...')

    return lon[0:n], lat[0:n], area