docs: def retinotopy_mesh_field(mesh, mdl,
                            polar_angle=None, eccentricity=None, weight=None,
                            weight_min=0, scale=1, sigma=None, shape=2, suffix=None,
                            max_eccentricity=Ellipsis,
                            max_polar_angle=180,
                            angle_type='both',
                            exclusion_threshold=None):
      '''
      retinotopy_mesh_field(mesh, model) yields a list that can be used with mesh_register as a
        potential term. This should generally be used in a similar fashion to retinotopy_anchors.
  
      Options:
        * polar_angle (default None) specifies that the given data should be used in place of the
          'polar_angle' or 'PRF_polar_angle'  property values. The given argument must be numeric and
          the same length as the the number of vertices in the mesh. If None is given, then the
          property value of the mesh is used; if a list is given and any element is None, then the
          weight for that vertex is treated as a zero. If the option is a string, then the property
          value with the same name isused as the polar_angle data.
        * eccentricity (default None) specifies that the given data should be used in places of the
          'eccentricity' or 'PRF_eccentricity' property values. The eccentricity option is handled 
          virtually identically to the polar_angle option.
        * weight (default None) specifies that the weight or scale of the data; this is handled
          generally like the polar_angle and eccentricity options, but may also be 1, indicating that
          all vertices with polar_angle and eccentricity values defined will be given a weight of 1.
          If weight is left as None, then the function will check for 'weight',
          'variance_explained', 'PRF_variance_explained', and 'retinotopy_weight' values and will use
          the first found (in that order). If none of these is found, then a value of 1 is assumed.
        * weight_min (default 0) specifies that the weight must be higher than the given value inn
          order to be included in the fit; vertices with weights below this value have their weights
          truncated to 0.
        * scale (default 1) specifies a constant by which to multiply all weights for all anchors; the
          value None is interpreted as 1.
        * shape (default 2.0) specifies the exponent in the harmonic function.
        * sigma (default None) specifies, if given as a number, the sigma parameter of the Gaussian
          potential function; if sigma is None, however, the potential function is harmonic.
        * suffix (default None) specifies any additional arguments that should be appended to the 
          potential function description list that is produced by this function; i.e., the 
          retinotopy_anchors function produces a list, and the contents of suffix, if given and not
          None, are appended to that list (see mesh_register).
        * max_eccentricity (default Ellipsis) specifies how the eccentricity portion of the potential
          field should be normalized. Specifically, in order to ensure that polar angle and
          eccentricity contribute roughly equally to the potential, this should be approximately the
          max eccentricity appearing in the data on the mesh. If the argument is the default then the
          actual max eccentricity will be used.
        * max_polar_angle (default: 180) is used the same way as the max_eccentricity function, but if
          Ellipsis is given, the value 180 is always assumed regardless of measured data.
        * exclusion_threshold (default None) specifies that if the initial norm of a vertex's gradient
          is greater than exclusion_threshold * std + median (where std and median are calculated over
          the vertices with non-zero gradients) then its weight is set to 0 and it is not kept as part
          of the potential field.
        * angle_type (default: None) specifies that only one type of angle should be included in the
          mesh; this may be one of 'polar', 'eccen', 'eccentricity', 'angle', or 'polar_angle'. If
          None, then both polar angle and eccentricity are included.
  
      Example:
       # The retinotopy_anchors function is intended for use with mesh_register, as follows:
       # Define our Schira Model:
       model = neuropythy.retinotopy_model()
       # Make sure our mesh has polar angle, eccentricity, and weight data:
       mesh.prop('polar_angle',  polar_angle_vertex_data);
       mesh.prop('eccentricity', eccentricity_vertex_data);
       mesh.prop('weight',       variance_explained_vertex_data);
       # register the mesh using the retinotopy and model:
       registered_mesh = neuropythy.registration.mesh_register(
          mesh,
          ['mesh', retinotopy_mesh_field(mesh, model)],
          max_step_size=0.05,
          max_steps=2000)
      '''
      #TODO: given a 3D mesh and a registered model, we should be able to return a 3D version of the
      # anchors by unprojecting them
      if pimms.is_str(mdl):
          mdl = retinotopy_model(mdl)
      if not isinstance(mdl, RetinotopyMeshModel):
          if isinstance(mdl, RegisteredRetinotopyModel):
              mdl = mdl.model
          if not isinstance(mdl, RetinotopyMeshModel):
              raise RuntimeError('given model is not a RetinotopyMeshModel instance!')
      if not hasattr(mdl, 'data') or 'polar_angle' not in mdl.data or 'eccentricity' not in mdl.data:
          raise ValueError('Retinotopy model does not have polar angle and eccentricity data')
      if not isinstance(mesh, CorticalMesh):
          raise RuntimeError('given mesh is not a CorticalMesh object!')
      n = mesh.vertex_count
      X = mesh.coordinates.T
      if weight_min is None: weight_min = 0
      # make sure we have our polar angle/eccen/weight values:
      # (weight is odd because it might be a single number, so handle that first)
      (polar_angle, eccentricity, weight) = [
          extract_retinotopy_argument(mesh, name, arg, default='empirical')
          for (name, arg) in [
                  ('polar_angle', polar_angle),
                  ('eccentricity', eccentricity),
                  ('weight', [weight for i in range(n)] \
                             if isinstance(weight, Number) or np.issubdtype(type(weight), np.float) \
                             else weight)]]
      # Make sure they contain no None/invalid values
      (polar_angle, eccentricity, weight) = _retinotopy_vectors_to_float(
          polar_angle, eccentricity, weight,
          weight_min=weight_min)
      if np.sum(weight > 0) == 0:
          raise ValueError('No positive weights found')
      idcs = [i for (i,w) in enumerate(weight) if w > 0]
      # Okay, let's get the model data ready
      mdl_1s = np.ones(mdl.forward.coordinates.shape[0])
      mdl_coords = np.dot(mdl.transform, np.vstack((mdl.forward.coordinates.T, mdl_1s)))[:2].T
      mdl_faces  = mdl.forward.triangles
      mdl_data   = np.asarray([mdl.data['polar_angle'], mdl.data['eccentricity']])
      # Get just the relevant weights and the scale
      wgts = weight[idcs] * (1 if scale is None else scale)
      # and the relevant polar angle/eccen data
      msh_data = np.asarray([polar_angle, eccentricity])[:,idcs]
      # format shape correctly
      shape = np.full((len(idcs)), float(shape), dtype=np.float32)
      # Last thing before constructing the field description: normalize both polar angle and eccen to
      # cover a range of 0-1:
      if max_eccentricity is Ellipsis: max_eccentricity = np.max(msh_data[1])
      if max_polar_angle  is Ellipsis: max_polar_angle  = 180
      if max_polar_angle is not None:
          msh_data[0] /= max_polar_angle
          mdl_data[0] /= max_polar_angle
      if max_eccentricity is not None:
          msh_data[1] /= max_eccentricity
          mdl_data[1] /= max_eccentricity
      # Check if we are making an eccentricity-only or a polar-angle-only field:
      if angle_type is not None:
          angle_type = angle_type.lower()
          if angle_type != 'both' and angle_type != 'all':
              convert = {'eccen':'eccen', 'eccentricity':'eccen', 'radius':'eccen',
                         'angle':'angle', 'polar_angle': 'angle', 'polar': 'angle'}
              angle_type = convert[angle_type]
              mdl_data = [mdl_data[0 if angle_type == 'angle' else 1]]
              msh_data = [msh_data[0 if angle_type == 'angle' else 1]]
      # okay, we've partially parsed the data that was given; now we can construct the final list of
      # instructions:
      if sigma is None:
          field_desc = ['mesh-field', 'harmonic', mdl_coords, mdl_faces, mdl_data, idcs, msh_data,
                        'scale', wgts, 'order', shape]
      else:
          if not hasattr(sigma, '__iter__'): sigma = [sigma for _ in wgts]
          field_desc = ['mesh-field', 'gaussian', mdl_coords, mdl_faces, mdl_data, idcs, msh_data,
                        'scale', wgts, 'order', shape, 'sigma', sigma]
      if suffix is not None: field_desc += suffix
      # now, if we want to exclude outliers, we do so here:
      if exclusion_threshold is not None:
          jpe = java_potential_term(mesh, field_desc)
          jcrds = to_java_doubles(mesh.coordinates)
          jgrad = to_java_doubles(np.zeros(mesh.coordinates.shape))
          jpe.calculate(jcrds,jgrad)
          gnorms = np.sum((np.asarray([[x for x in row] for row in jgrad])[:, idcs])**2, axis=0)
          gnorms_pos = gnorms[gnorms > 0]
          mdn = np.median(gnorms_pos)
          std = np.std(gnorms_pos)
          gn_idcs = np.where(gnorms > mdn + std*3.5)[0]
          for i in gn_idcs: wgts[i] = 0;
      return field_desc