code: 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_fdata() > 0).squeeze())
          zmap = zimg.get_fdata().squeeze()[xyz]
      else:
          xyz = np.where(mask.get_fdata() > 0)
          zmap = zimg.get_fdata()[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 "zmax" not in nulls:
          nulls["zmax"] = "bonferroni"
      if "smax" not in nulls:
          nulls["smax"] = None
      if "s" not 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 isinstance(nulls["zmax"], np.ndarray):
              p = simulated_pvalue(zscore, nulls["zmax"])
          elif nulls["zmax"] == "bonferroni":
              p = bonferroni(pval, nvoxels)
          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
  
imports: import numpy
  import scipy.stats
  from nibabel import Nifti1Image
  from nibabel.affines import apply_affine
  from ..algorithms.graph.field import field_from_graph_and_data
  from ..algorithms.graph.graph import wgraph_from_3d_grid
  from ..algorithms.statistics import empirical_pvalue
  from ..io.nibcompat import get_affine
  from .glm import glm
  from .group.permutation_test import permutation_test_onesample, permutation_test_twosample
label_desc: Use image node_id 3 for calling a global function with example usage: cluster_stats(zimg, mask, height_th, height_control, cluster_th, nulls) and returns: clusters,  info, None,  None