code: def _extract_clusters_from_diam(
      labels, T, XYZ, th, diam, k, nCC, CClabels
  ):
      """
      This recursive function modifies the `labels` input array.
      """
      clust_label = 0
      for i in range(nCC):
          # print "Searching connected component ", i, " out of ", nCC
          I = np.where(CClabels == i)[0]
          extCC = len(I)
          if extCC <= (diam + 1) ** 3:
              diamCC = max_dist(XYZ, I, I)
          else:
              diamCC = diam + 1
          if diamCC <= diam:
              labels[I] = np.zeros(extCC, int) + clust_label
              # print "cluster ", clust_label, ", diam = ", diamCC
              # print "ext = ", len(I), ", diam = ", max_dist(XYZ,I,I)
              clust_label += 1
          else:
              # build the field
              p = len(T[I])
              F = field_from_graph_and_data(
                  wgraph_from_3d_grid(XYZ[:, I].T, k),
                  np.reshape(T[I], (p, 1)),
              )
              # compute the blobs
              idx, parent, label = F.threshold_bifurcations(0, th)
              nidx = np.size(idx)
              height = np.array(
                  [
                      np.ceil(np.sum(label == i) ** (1.0 / 3))
                      for i in np.arange(nidx)
                  ]
              )
              # root = nidx-1
              root = np.where(np.arange(nidx) == parent)[0]
              # Can constraint be met within current region?
              Imin = I[T[I] >= height[root]]
              extmin = len(Imin)
              if extmin <= (diam + 1) ** 3:
                  dmin = max_dist(XYZ, Imin, Imin)
              else:
                  dmin = diam + 1
              if (
                  dmin <= diam
              ):  # If so, search for the largest cluster meeting the constraint
                  Iclust = Imin  # Smallest cluster
                  J = I[T[I] < height[root]]  # Remaining voxels
                  argsortTJ = np.argsort(T[J])[
                      ::-1
                  ]  # Sorted by decreasing T values
                  l = 0
                  L = np.array([J[argsortTJ[l]]], int)
                  diameter = dmin
                  new_diameter = max(dmin, max_dist(XYZ, Iclust, L))
                  while new_diameter <= diam:
                      # print "diameter = " + str(new_diameter)
                      # sys.stdout.flush()
                      Iclust = np.concatenate((Iclust, L))
                      diameter = new_diameter
                      # print "diameter = ", diameter
                      l += 1
                      L = np.array([J[argsortTJ[l]]], int)
                      new_diameter = max(
                          diameter, max_dist(XYZ, Iclust, L)
                      )
                  labels[Iclust] = (
                      np.zeros(len(Iclust), int) + clust_label
                  )
                  # print "cluster ", clust_label, ", diam = ", diameter
                  # print "ext = ", len(Iclust), ", diam = ", max_dist(XYZ,Iclust,Iclust)
                  clust_label += 1
              else:  # If not, search inside sub-regions
                  # print "Searching inside sub-regions "
                  Irest = I[T[I] > height[root]]
                  rest_labels = extract_clusters_from_diam(
                      T[Irest], XYZ[:, Irest], th, diam, k
                  )
                  rest_labels[rest_labels >= 0] += clust_label
                  clust_label = rest_labels.max() + 1
                  labels[Irest] = rest_labels
      return labels
  
imports: import numpy
  import warnings
  from nipy.algorithms.graph import wgraph_from_3d_grid
  from nipy.algorithms.graph.field import Field, field_from_graph_and_data
  from ..utils import zscore
  from .onesample import stat
  from .onesample import stat_mfx
  from .twosample import stat
  from .twosample import stat_mfx
label_desc: Use image node_id 4 for calling a global function with example usage: _extract_clusters_from_diam(labels, T, XYZ, th, diam, k, nCC, CClabels) and returns: labels