def cluster_stats(zimg, mask, height_th, height_control='fpr',
                  cluster_th=0, nulls={}):
    """
    Return a list of clusters, each cluster being represented by a
    dictionary. Clusters are sorted by descending size order. Within
    each cluster, local maxima are sorted by descending depth order.

    Parameters
    ----------
    zimg: z-score image
    mask: mask image
    height_th: cluster forming threshold
    height_control: string
            false positive control meaning of cluster forming
            threshold: 'fpr'|'fdr'|'bonferroni'|'none'
    cluster_th: cluster size threshold
    null_s : cluster-level calibration method: None|'rft'|array

    Notes
    -----
    This works only with three dimensional data
    """
    # Masking
    if len(mask.shape) > 3:
        xyz = np.where((mask.get_data() > 0).squeeze())
        zmap = zimg.get_data().squeeze()[xyz]
    else:
        xyz = np.where(mask.get_data() > 0)
        zmap = zimg.get_data()[xyz]

    xyz = np.array(xyz).T
    nvoxels = np.size(xyz, 0)

    # Thresholding
    if height_control == 'fpr':
        zth = sp_stats.norm.isf(height_th)
    elif height_control == 'fdr':
        zth = empirical_pvalue.gaussian_fdr_threshold(zmap, height_th)
    elif height_control == 'bonferroni':
        zth = sp_stats.norm.isf(height_th / nvoxels)
    else: ## Brute-force thresholding
        zth = height_th
    pth = sp_stats.norm.sf(zth)
    above_th = zmap > zth
    if len(np.where(above_th)[0]) == 0:
        return None, None ## FIXME
    zmap_th = zmap[above_th]
    xyz_th = xyz[above_th]

    # Clustering
    ## Extract local maxima and connex components above some threshold
    ff = field_from_graph_and_data(wgraph_from_3d_grid(xyz_th, k=18), zmap_th)
    maxima, depth = ff.get_local_maxima(th=zth)
    labels = ff.cc()
    ## Make list of clusters, each cluster being a dictionary
    clusters = []
    for k in range(labels.max() + 1):
        s = np.sum(labels == k)
        if s >= cluster_th:
            in_cluster = labels[maxima] == k
            m = maxima[in_cluster]
            d = depth[in_cluster]
            sorted = d.argsort()[::-1]
            clusters.append({'size': s,
                             'maxima': m[sorted],
                             'depth': d[sorted]})

    ## Sort clusters by descending size order
    clusters.sort(key=lambda c : c['size'], reverse=True)

    # FDR-corrected p-values
    fdr_pvalue = empirical_pvalue.gaussian_fdr(zmap)[above_th]

    # Default "nulls"
    if not 'zmax' in nulls:
        nulls['zmax'] = 'bonferroni'
    if not 'smax' in nulls:
        nulls['smax'] = None
    if not 's' in nulls:
        nulls['s'] = None

    # Report significance levels in each cluster
    for c in clusters:
        maxima = c['maxima']
        zscore = zmap_th[maxima]
        pval = sp_stats.norm.sf(zscore)
        # Replace array indices with real coordinates
        c['maxima'] = apply_affine(get_affine(zimg), xyz_th[maxima])
        c['zscore'] = zscore
        c['pvalue'] = pval
        c['fdr_pvalue'] = fdr_pvalue[maxima]

        # Voxel-level corrected p-values
        p = None
        if nulls['zmax'] == 'bonferroni':
            p = bonferroni(pval, nvoxels)
        elif isinstance(nulls['zmax'], np.ndarray):
            p = simulated_pvalue(zscore, nulls['zmax'])
        c['fwer_pvalue'] = p

        # Cluster-level p-values (corrected)
        p = None
        if isinstance(nulls['smax'], np.ndarray):
            p = simulated_pvalue(c['size'], nulls['smax'])
        c['cluster_fwer_pvalue'] = p

        # Cluster-level p-values (uncorrected)
        p = None
        if isinstance(nulls['s'], np.ndarray):
            p = simulated_pvalue(c['size'], nulls['s'])
        c['cluster_pvalue'] = p

    # General info
    info = {'nvoxels': nvoxels,
            'threshold_z': zth,
            'threshold_p': pth,
            'threshold_pcorr': bonferroni(pth, nvoxels)}

    return clusters, info