Loss function to fit unordered points/shape in 2d/3d space

For simplicity I am going start with a toy example first.

Lets suppose we have a set of n points Y in the space, distributed with the shape of the letter M and you want to fit those using 4 lines. These points are not ordered in the sense that the order of the points in the vector Y is shuffled.

So I built a my function f that given 5 points P, it generates 4 lines concatenated joining those 5 points. I want to find the 5 points P that I can use to best describe the data. This is, the lines are close to the points.

So basically my though here is, given values x in [0,1] and my function f(x, P) I can sample points from this function and compute the loss against the Y, Loss(f(x,P), Y). For example if x is torch.linspace(0, 1, 100) I am gong to draw 100 points from my function. I can also draw n points. Whatever it works.

If the points of the ground truth would be ordered and kind of evenly distributed I can use the MSELoss function and I can get the 5 points (I already tried this and it works). However if the Y points are unordered this loss is no longer working.
My question is then, which kind of loss function I can use here to find P?

My problem is a little bit more general and I have a set of points Y that follow a line that curves through the 3d space so I would like to use m concatenated lines given by f(P}) whose distribution best matches Y.

Apparently the The Chamfer distance can work. This distance is computed by summing the squared distances between nearest neighbor correspondences of two point clouds.