I have a set of points which I need to sort into a set of pareto frontiers (see example diagram here)

I have a method for doing this, but it becomes extremelly slow for a large number of points, 60k+ points.

The method I am using involves iteratively invoking a library function that can find a single pareto front:

1. Start with all points remaining_points: Vec<Point>
2. Create an empty all_frontiers: Vec<Vec<Point>>
3. Filter for points on the pareto front frontier: Vec<Point>, and add them to all_frontiers
4. Remove frontier points from remaining_points
5. Repeat 3-4 until remaining_points.is_empty()

I do not know of (nor could I find) any better approaches to this problem. Can someone give suggestions?