docs: def cortex_to_image_interpolation(obj, mask=None, affine=None, method='linear', shape=None):
      '''
      cortex_to_image_interpolation(obj) yields a tuple (indices, interp) where indices is a tuple of
        voxel indices and interp is an interpolation matrix that converts a vector of cortical
        surface vertex values into a vector of of voxel values with the same ordering as those given
        in indices.
  
      The argument obj may be either a subject or a cortex or a tuple of (lh, rh) cortices. If the obj
      specifies two cortices, then the interp matrix will be arranged such that the vertex-value
      vector with which interp should be multiplied must first list the LH values then the RH values;
      i.e., image[indices] = interp.dot(join(lh_values, rh_values)).
  
      The following options are accepted:
        * mask (default: None) specifies the mask of vertices that should be interpolated; if this is
          None, then this will attempt to use the gray_indices of the object if the object is a
          subject, otherwise this is equivalent to 'all'. The special value 'all' indicates that all
          vertices overlapping with the cortex should be interpolated. Mask may be given as a boolean
          mask or as a set or as a tuple/matrix of indices equivalent to numpy.where(binary_mask) or
          its transpose. The value indices that is returned will be identical to mask if mask is given
          in the same format.
        * affine (default: None) specifies the affine transform that is used to align the vertices
          with the voxels. In voxel-space, the voxel with index (i,j,k) is centered at at the point
          (i,j,k) and is 1 unit wide in every direction. If the value None is given, will attempt to
          use the object's vertex_to_voxel_matrix if object is a subject, otherwise will use a
          FreeSurfer-like orientation that places the vertex origin in the center of the image.
        * shape (default: None) specifies the size of the resulting image as a tuple. If None is given
          then the function attempts to deduce the correct shape; if obj is a subject, then its
          image_dimensions tuple is used. Otherwise, the size is deduced from the mask, if possible;
          if it cannot be deduced from the mask, then (256, 256, 256) is used.
        * method (default: 'linear') specifies the method to use. May be 'linear', 'heaviest', or
          'nearest'.
      '''
      # get the min/max values of the coordinates (for deducing sizes, if necessary)
      if mask is None:
          if shape is None and affine is None:
              # we have no way to deduce anything
              shape = (256,256,256)
              affine = ([[-1,0,0],[0,0,-1],[0,1,0]], [128,128,128])
          elif shape is None:
              shape = (256,256,256)
          elif affine is None:
              affine = ([[-1,0,0],[0,0,-1],[0,1,0]], np.asarray(shape) / 2)
          mask = obj.gray_indices if isinstance(obj, Subject) else 'all'
      # okay, having handled that no-arg case, lets parse the argument we have
      if pimms.is_matrix(mask):
          # we take this to be the list; we don't chante its order
          if not is_tuple(mask) or len(mask) != 3:
              mask = np.asarray(mask, dtype=np.int)
              mask = tuple(mask.T if mask.shape[0] != 3 else mask)
      elif isinstance(mask, colls.Set):
          # we have to convert this into a propert mask
          mask = np.asarray(list(mask)).T
          tmp = np.full(shape, False)
          tmp[tuple(mask)] = True
          mask = np.where(tmp)
      elif pimms.is_array(mask, dims=3):
          if shape is None: shape = mask.shape
          mask = np.where(mask)
      elif pimms.is_str(mask) and mask.lower() == 'all':
          if shape is None: shape = (256,256,256)
          mask = np.where(np.ones(shape))
      else:
          raise ValueError('Could not understand mask argument')
      # at this point, we can deduce shape from mask and affine from shape
      if shape is None:  shape  = np.asarray(np.ceil(np.max(mask, axis=0)), dtype=np.int)
      if affine is None:
          if isinstance(obj, Subject):
              affine = obj.vertex_to_voxel_matrix
          else:
              affine = ([[-1,0,0],[0,0,-1],[0,1,0]], np.asarray(shape) / 2)
      affine = to_affine(affine, 3)
      hems = (obj.lh, obj.rh) if isinstance(obj, Subject) else (obj,)
      # all arguments are basically pre-processed; we just need to make the interpolation
      method = 'auto' if method is None else method.lower()
      if method in ['linear', 'auto', 'automatic']:
          interp = [_vertex_to_voxel_linear_interpolation(h, mask, shape, affine) for h in hems]
          if len(interp) == 1: interp = interp[0]
          else: interp = sps.hstack(interp)
      elif method in ['lines', 'line']:
          interp = [_vertex_to_voxel_lines_interpolation(h, mask, shape, affine) for h in hems]
          if len(interp) == 1: interp = interp[0]
          else: interp = sps.hstack(interp)
      elif method in ['heaviest', 'heavy', 'weight', 'weightiest']:
          interp = [_vertex_to_voxel_linear_interpolation(h, mask, shape, affine) for h in hems]
          if len(interp) == 1: interp = interp[0]
          else: interp = sps.hstack(interp)
          # convert to binary matrix:
          (rs,cs) = interp.shape
          argmaxs = np.asarray(interp.argmax(axis=1))[:,0]
          return sps.csr_matrix((np.ones(rs, dtype=np.int), (range(rs), argmaxs)),
                                shape=interp.shape,
                                dtype=np.int)
      elif method in ['nearest', 'near', 'nearest-neighbor', 'nn']:
          aff = npla.inv(affine)
          interp = [_vertex_to_voxel_nearest_interpolation(h, mask, aff) for h in hems]
          if len(interp) == 1: interp = interp[0]
          else: interp = sps.hstack(interp)
      else:
          raise ValueError('unsupported method: %s' % method)
      # That's it; we have the interp matrix and the indices
      return (mask, interp)