code: def connectivity_matrix(
      streamlines,
      affine,
      label_volume,
      inclusive=False,
      symmetric=True,
      return_mapping=False,
      mapping_as_streamlines=False,
  ):
      """Count the streamlines that start and end at each label pair.
  
      Parameters
      ----------
      streamlines : sequence
          A sequence of streamlines.
      affine : array_like (4, 4)
          The mapping from voxel coordinates to streamline coordinates.
          The voxel_to_rasmm matrix, typically from a NIFTI file.
      label_volume : ndarray
          An image volume with an integer data type, where the intensities in the
          volume map to anatomical structures.
      inclusive: bool
          Whether to analyze the entire streamline, as opposed to just the
          endpoints. False by default.
      symmetric : bool, True by default
          Symmetric means we don't distinguish between start and end points. If
          symmetric is True, ``matrix[i, j] == matrix[j, i]``.
      return_mapping : bool, False by default
          If True, a mapping is returned which maps matrix indices to
          streamlines.
      mapping_as_streamlines : bool, False by default
          If True voxel indices map to lists of streamline objects. Otherwise
          voxel indices map to lists of integers.
  
      Returns
      -------
      matrix : ndarray
          The number of connection between each pair of regions in
          `label_volume`.
      mapping : defaultdict(list)
          ``mapping[i, j]`` returns all the streamlines that connect region `i`
          to region `j`. If `symmetric` is True mapping will only have one key
          for each start end pair such that if ``i < j`` mapping will have key
          ``(i, j)`` but not key ``(j, i)``.
      """
  
      # Error checking on label_volume
      kind = label_volume.dtype.kind
      labels_positive = (kind == "u") or (
          (kind == "i") and (label_volume.min() >= 0)
      )
      valid_label_volume = labels_positive and label_volume.ndim == 3
      if not valid_label_volume:
          raise ValueError(
              "label_volume must be a 3d integer array with"
              "non-negative label values"
          )
  
      matrix = np.zeros(
          (np.max(label_volume) + 1, np.max(label_volume) + 1),
          dtype=np.int64,
      )
  
      mapping = defaultdict(list)
      lin_T, offset = _mapping_to_voxel(affine)
  
      if inclusive:
          for i, sl in enumerate(streamlines):
              sl = _to_voxel_coordinates(sl, lin_T, offset)
              x, y, z = sl.T
              if symmetric:
                  crossed_labels = np.unique(label_volume[x, y, z])
              else:
                  crossed_labels = np.unique(
                      label_volume[x, y, z], return_index=True
                  )
                  crossed_labels = crossed_labels[0][
                      np.argsort(crossed_labels[1])
                  ]
  
              for comb in combinations(crossed_labels, 2):
                  matrix[comb] += 1
  
                  if return_mapping:
                      if mapping_as_streamlines:
                          mapping[comb].append(streamlines[i])
                      else:
                          mapping[comb].append(i)
  
      else:
          streamlines_end = np.array(
              [sl[0 :: len(sl) - 1] for sl in streamlines]
          )
          streamlines_end = _to_voxel_coordinates(
              streamlines_end, lin_T, offset
          )
          x, y, z = streamlines_end.T
          if symmetric:
              end_labels = np.sort(label_volume[x, y, z], axis=0)
          else:
              end_labels = label_volume[x, y, z]
          np.add.at(matrix, (end_labels[0].T, end_labels[1].T), 1)
  
          if return_mapping:
              if mapping_as_streamlines:
                  for i, (a, b) in enumerate(end_labels.T):
                      mapping[a, b].append(streamlines[i])
              else:
                  for i, (a, b) in enumerate(end_labels.T):
                      mapping[a, b].append(i)
  
      if symmetric:
          matrix = np.maximum(matrix, matrix.T)
  
      if return_mapping:
          return (matrix, mapping)
      else:
          return matrix
  
imports: from collections import defaultdict, OrderedDict
  from functools import wraps
  from itertools import combinations, groupby
  from warnings import warn
  import numpy
  from nibabel.affines import apply_affine
  from scipy.spatial.distance import cdist
  from dipy.core.geometry import dist_to_corner
  from dipy.tracking import metrics
  from dipy.tracking.vox2track import _streamlines_in_mask
  from dipy.tracking._utils import _mapping_to_voxel, _to_voxel_coordinates
label_desc: Use image node_id 2 for calling a global function with example usage: connectivity_matrix(streamlines, affine, label_volume, inclusive, symmetric, return_mapping, mapping_as_streamlines) and returns: matrix,  mapping, matrix