docs: def neighborhood_cortical_magnification(mesh, coordinates):
      '''
      neighborhood_cortical_magnification(mesh, visual_coordinates) yields a list of neighborhood-
      based cortical magnification values for the vertices in the given mesh if their visual field
      coordinates are given by the visual_coordinates matrix (must be like [x_values, y_values]). If
      either x-value or y-value of a coordinate is either None or numpy.nan, then that cortical
      magnification value is numpy.nan.
      '''
      idcs = _cmag_coord_idcs(coordinates)
      neis = mesh.tess.indexed_neighborhoods
      coords_vis = np.asarray(coordinates if len(coordinates) == 2 else coordinates.T)
      coords_srf = mesh.coordinates
      res = np.full((mesh.vertex_count, 3), np.nan, dtype=np.float)
      res = np.array([row for row in [(np.nan,np.nan,np.nan)] for _ in range(mesh.vertex_count)],
                     dtype=np.float)
      for idx in idcs:
          nei = neis[idx]
          pts_vis = coords_vis[:,nei]
          pts_srf = coords_srf[:,nei]
          x0_vis = coords_vis[:,idx]
          x0_srf = coords_srf[:,idx]
          if any(u is None or np.isnan(u) for pt in pts_vis for u in pt): continue
          # find tangential, radial, and areal magnifications
          x0col_vis = np.asarray([x0_vis]).T
          x0col_srf = np.asarray([x0_srf]).T
          # areal is easy
          voronoi_vis = (pts_vis - x0col_vis) * 0.5 + x0col_vis
          voronoi_srf = (pts_srf - x0col_srf) * 0.5 + x0col_srf
          area_vis = np.sum([geo.triangle_area(x0_vis, a, b)
                             for (a,b) in zip(voronoi_vis.T, np.roll(voronoi_vis, 1, axis=1).T)])
          area_srf = np.sum([geo.triangle_area(x0_srf, a, b)
                             for (a,b) in zip(voronoi_srf.T, np.roll(voronoi_srf, 1, axis=1).T)])
          res[idx,2] = np.inf if np.isclose(area_vis, 0) else area_srf/area_vis
          # radial and tangentual we do together because they are very similar:
          # find the intersection lines then add up their distances along the cortex
          pts_vis = voronoi_vis
          pts_srf = voronoi_srf
          segs_srf = (pts_srf, np.roll(pts_srf, -1, axis=1))
          segs_vis = (pts_vis, np.roll(pts_vis, -1, axis=1))
          segs_vis_t = np.transpose(segs_vis, (2,0,1))
          segs_srf_t = np.transpose(segs_srf, (2,0,1))
          x0norm_vis = npla.norm(x0_vis)
          if not np.isclose(x0norm_vis, 0):
              dirvecs = x0_vis / x0norm_vis
              dirvecs = np.asarray([dirvecs, [-dirvecs[1], dirvecs[0]]])
              for dirno in [0,1]:
                  dirvec = dirvecs[dirno]
                  line = (x0_vis, x0_vis + dirvec)
                  try:
                      isects_vis = np.asarray(geo.line_segment_intersection_2D(line, segs_vis))
                      # okay, these will all be nan but two of them; they are the points we care about
                      isect_idcs = np.unique(np.where(np.logical_not(np.isnan(isects_vis)))[1])
                  except Exception:
                      isect_idcs = []
                  if len(isect_idcs) != 2:
                      res[idx,dirno] = np.nan
                      continue
                  isects_vis = isects_vis[:,isect_idcs].T
                  # we need the distance in visual space
                  len_vis = npla.norm(isects_vis[0] - isects_vis[1])
                  if np.isclose(len_vis, 0): res[idx,dirno] = np.inf
                  else:
                      # we also need the distances on the surface: find the points by projection
                      fsegs_srf = segs_srf_t[isect_idcs]
                      fsegs_vis = segs_vis_t[isect_idcs]
                      s02lens_vis = npla.norm(fsegs_vis[:,0] - fsegs_vis[:,1], axis=1)
                      s01lens_vis = npla.norm(fsegs_vis[:,0] - isects_vis, axis=1)
                      vecs_srf = fsegs_srf[:,1] - fsegs_srf[:,0]
                      s02lens_srf = npla.norm(vecs_srf, axis=1)
                      isects_srf = np.transpose([(s01lens_vis/s02lens_vis)]) * vecs_srf \
                                   + fsegs_srf[:,0]
                      len_srf = np.sum(npla.norm(isects_srf - x0_srf, axis=1))
                      res[idx,dirno] = len_srf / len_vis
      # That's it!
      return res