def list_and_sort_gaussians(img, patch=None, root=None, sort_by='index'):
    """Returns sorted lists of Gaussians and their names and patch names.

    patch - can be "single", "gaussian", "source", or None

    Returns (outlist, outnames, patchnames)
    outlist is [[g1, g2, g3], [g4], ...]
    outnames is [['root_i2_s1_g1', 'root_i2_s1_g2', 'root_i2_s1_g3'], ...]
    patchnames is ['root_patch_s1', 'root_patch_s2', ...]

    The names are root_iXX_sXX_gXX (or wXX_iXX_sXX_gXX for wavelet Gaussians)
    """
    import numpy as N
    from . import functions as func
    if root is None:
        root = img.parentname
    gauslist = []
    gausname = []
    outlist = []
    outnames = []
    patchnames = []
    patchnames_sorted = []
    gausflux = []
    gausindx = []
    patchflux = []
    patchindx = []
    patchnums = []
    use_mask = False
    if patch not in ['single', 'gaussian', 'source', None]:
        mask_file = img.opts.bbs_patches_mask
        patches_mask, hdr = func.read_image_from_file(mask_file, img, img.indir
            )
        use_mask = True
        act_pixels = patches_mask[0, 0]
        rank = len(act_pixels.shape)
        import scipy.ndimage as nd
        connectivity = nd.generate_binary_structure(rank, rank)
        mask_labels, count = nd.label(act_pixels, connectivity)
    src_list = img.sources
    for src in src_list:
        for g in src.gaussians:
            gauslist.append(g)
            gausflux.append(g.total_flux)
            gausindx.append(g.gaus_num)
            jstr = '_w' + str(g.jlevel)
            gausname.append(root + jstr + '_i' + str(src.island_id) + '_s' +
                str(src.source_id) + '_g' + str(g.gaus_num))
            if patch == 'gaussian':
                outlist.append(gauslist)
                outnames.append(gausname)
                patchnames.append(root + '_patch' + jstr + '_g' + str(g.
                    gaus_num))
                patchflux.append(N.sum(gausflux))
                patchindx.append(g.gaus_num)
                gauslist = []
                gausname = []
                gausflux = []
                gausindx = []
            if use_mask:
                patchnums.append(mask_labels[g.centre_pix[0], g.centre_pix[1]])
        if patch == 'source':
            sorted_gauslist = list(gauslist)
            sorted_gausname = list(gausname)
            if sort_by == 'flux':
                indx = N.argsort(N.array(gausflux)).tolist()
                indx.reverse()
            elif sort_by == 'index':
                indx = N.argsort(N.array(gausindx)).tolist()
            else:
                indx = range(len(gausindx))
            for i, si in enumerate(indx):
                sorted_gauslist[i] = gauslist[si]
                sorted_gausname[i] = gausname[si]
            outlist.append(sorted_gauslist)
            outnames.append(sorted_gausname)
            patchnames.append(root + '_patch' + '_s' + str(src.source_id))
            patchflux.append(N.sum(gausflux))
            patchindx.append(src.source_id)
            gauslist = []
            gausname = []
            gausflux = []
    if use_mask:
        unique_patch_ids = set(patchnums)
        if 0 in unique_patch_ids:
            from . import mylogger
            mylog = mylogger.logging.getLogger('PyBDSM.write_gaul')
            mylog.warning(
                'Some sources fall outside of the regions defined in the mask file. These sources are not included in the output sky model.'
                )
        for p in unique_patch_ids:
            if p != 0:
                in_patch = N.where(patchnums == p)
                outlist.append(N.array(gauslist)[in_patch].tolist())
                outnames.append(N.array(gausname)[in_patch].tolist())
                patchnames.append('patch_' + str(p))
                patchflux.append(N.sum(N.array(gausflux)[in_patch]))
                patchindx.append(p)
    if patch == 'single' or patch is None:
        outlist = [list(gauslist)]
        outlist_sorted = [list(gauslist)]
        outnames = [list(gausname)]
        outnames_sorted = [list(gausname)]
        if patch == 'single':
            patchnames = [root + '_patch']
        else:
            patchnames = [None]
        if sort_by == 'flux':
            indx = N.argsort(N.array(gausflux)).tolist()
            indx.reverse()
        elif sort_by == 'index':
            indx = N.argsort(N.array(gausindx)).tolist()
        else:
            indx = list(range(len(gausindx)))
        for i, si in enumerate(indx):
            outlist_sorted[0][i] = outlist[0][si]
            outnames_sorted[0][i] = outnames[0][si]
            patchnames_sorted = list(patchnames)
    else:
        outlist_sorted = list(outlist)
        outnames_sorted = list(outnames)
        patchnames_sorted = list(patchnames)
        if sort_by == 'flux':
            indx = N.argsort(N.array(patchflux)).tolist()
            indx.reverse()
        elif sort_by == 'index':
            indx = N.argsort(N.array(patchindx)).tolist()
        else:
            indx = list(range(len(gausindx)))
        for i, si in enumerate(indx):
            outlist_sorted[i] = outlist[si]
            outnames_sorted[i] = outnames[si]
            patchnames_sorted[i] = patchnames[si]
    return outlist_sorted, outnames_sorted, patchnames_sorted
