code: def classify(self, points, edge=True, col=1, already_sorted=False):
          
          xy = points[:, :2].copy()
          radius = 1e-10 if edge else 0
          n_points = xy.shape[0]
          self.lock() # Precomputes polygons information
          
          # Sorting the points
          if not already_sorted:
              from pyny3d.utils import sort_numpy
              # Sorting by the y column can be faster if set_mesh has been
              # used to generate the points.        
              points, order_back = sort_numpy(xy, col, order_back=True)
          
          # Computing PiP
          sol = np.ones(n_points, dtype=int)*-1
          for i, polygon in enumerate(self):
              pip = polygon.pip(points, sorted_col=col, radius=radius)
              sol[pip] = i
          
          # Return
          if not already_sorted:
              return sol[order_back]
          else:
              return sol
docs: Calculates the belonging relationship between the polygons
          in the Surface and a set of points.
          
          This function enhances the performance of ``Polygon.contains()``
  		when used with multiple non-overlapping polygons (stored in a 
  		Surface) by verifying only the points which are inside the z=0 
  		bounding box of each polygon. To do it fast, it sorts the points
  		and then apply ``Polygon.pip()`` for each Polygon.
          
          :param points: list of (x, y, z) or (x, y) coordinates of the
              points to check. (The z value will not be taken into 
              account).
          :type points: ndarray (shape=(N, 2 or 3))
          :param edge: If True, consider the points in a Polygon's edge
              inside a Polygon.
          :type edge: bool
          :param col: Column to sort or already sorted.
          :type col: int
          :param already_sorted: If True, the method will consider that
              the *points* are already sorted by the column *col*.
          :type already_sorted: bool
          :returns: Index of the Polygon to which each point belongs. 
              -1 if outside the Surface.
          :rtype: ndarray (dtpye=int)