code: def decompose3d(shape, dim=4):
      """
      Return all (dim-1)-dimensional simplices in a triangulation
      of a cube of a given shape. The vertices in the triangulation
      are indices in a 'flattened' array of the specified shape.
      """
  
      # First do the interior contributions.
      # We first figure out which vertices, edges, triangles, tetrahedra
      # are uniquely associated with an interior voxel
  
      unique = {}
      strides = np.empty(shape, np.bool_).strides
      union = join_complexes(
          *[
              cube_with_strides_center((0, 0, -1), strides),
              cube_with_strides_center((0, -1, 0), strides),
              cube_with_strides_center((0, -1, -1), strides),
              cube_with_strides_center((-1, 0, 0), strides),
              cube_with_strides_center((-1, 0, -1), strides),
              cube_with_strides_center((-1, -1, 0), strides),
              cube_with_strides_center((-1, -1, -1), strides),
          ]
      )
  
      c = cube_with_strides_center((0, 0, 0), strides)
      for i in range(4):
          unique[i + 1] = c[i + 1].difference(union[i + 1])
  
      if dim in unique and dim > 1:
          d = unique[dim]
  
          for i in range(shape[0] - 1):
              for j in range(shape[1] - 1):
                  for k in range(shape[2] - 1):
                      index = (
                          i * strides[0]
                          + j * strides[1]
                          + k * strides[2]
                      )
                      for l in d:
                          yield [index + ii for ii in l]
  
      # There are now contributions from three two-dimensional faces
  
      for _strides, _shape in zip(
          [
              (strides[0], strides[1]),
              (strides[0], strides[2]),
              (strides[1], strides[2]),
          ],
          [
              (shape[0], shape[1]),
              (shape[0], shape[2]),
              (shape[1], shape[2]),
          ],
      ):
          unique = {}
          union = join_complexes(
              *[
                  cube_with_strides_center((0, -1), _strides),
                  cube_with_strides_center((-1, 0), _strides),
                  cube_with_strides_center((-1, -1), _strides),
              ]
          )
  
          c = cube_with_strides_center((0, 0), _strides)
          for i in range(3):
              unique[i + 1] = c[i + 1].difference(union[i + 1])
  
          if dim in unique and dim > 1:
              d = unique[dim]
  
              for i in range(_shape[0] - 1):
                  for j in range(_shape[1] - 1):
                      index = i * _strides[0] + j * _strides[1]
                      for l in d:
                          yield [index + ii for ii in l]
  
      # Finally the one-dimensional faces
  
      for _stride, _shape in zip(strides, shape):
          unique = {}
          union = cube_with_strides_center((-1,), [_stride])
          c = cube_with_strides_center((0,), [_stride])
          for i in range(2):
              unique[i + 1] = c[i + 1].difference(union[i + 1])
  
          if dim in unique and dim > 1:
              d = unique[dim]
  
              for i in range(_shape - 1):
                  index = i * _stride
                  for l in d:
                      yield [index + ii for ii in l]
  
      if dim == 1:
          for i in range(np.prod(shape)):
              yield i
  
file_path: nipy/nipy/algorithms/statistics/utils.py
imports: from itertools import combinations
  import numpy
  from scipy.stats import norm
label_desc: Use image node_id 7 for calling a global function with example usage: decompose3d(shape, dim) without return types