code: def density_map(streamlines, affine, vol_dims):
      """Count the number of unique streamlines that pass through each voxel.
  
      Parameters
      ----------
      streamlines : iterable
          A sequence of streamlines.
      affine : array_like (4, 4)
          The mapping from voxel coordinates to streamline points.
          The voxel_to_rasmm matrix, typically from a NIFTI file.
      vol_dims : 3 ints
          The shape of the volume to be returned containing the streamlines
          counts
  
      Returns
      -------
      image_volume : ndarray, shape=vol_dims
          The number of streamline points in each voxel of volume.
  
      Raises
      ------
      IndexError
          When the points of the streamlines lie outside of the return volume.
  
      Notes
      -----
      A streamline can pass through a voxel even if one of the points of the
      streamline does not lie in the voxel. For example a step from [0,0,0] to
      [0,0,2] passes through [0,0,1]. Consider subsegmenting the streamlines when
      the edges of the voxels are smaller than the steps of the streamlines.
  
      """
      lin_T, offset = _mapping_to_voxel(affine)
      counts = np.zeros(vol_dims, "int")
      for sl in streamlines:
          inds = _to_voxel_coordinates(sl, lin_T, offset)
          i, j, k = inds.T
          # this takes advantage of the fact that numpy's += operator only
          # acts once even if there are repeats in inds
          counts[i, j, k] += 1
      return counts
  
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 1 for calling a global function with example usage: density_map(streamlines, affine, vol_dims) and returns: counts