Efficient algorithm for sorting 2d points into a set of pareto frontiers

I have a set of points which I need to sort into a set of pareto frontiers, example below:

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?

I arrived at a >10x faster algorithm specific to my use case (integer points, 2d, points will often share same x or y value).

1. Find the dimension with the smallest range (x or y)
2. Group points by that dimension into buckets
3. Sort each group by the other dimension
4. Skim the first point of each group, and calculate the pareto frontier
5. For each point in the frontier, remove it from its group in the bucket. If bucket is empty, also remove it
6. Loop 4-5 until buckets are empty