def classify(self, points, edge=True, col=1, already_sorted=False):
        """
        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)
        """
        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
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)