Order 2d points based on distance from each other

Question

Input image:

import cv2
import matplotlib.pyplot as plt
from skimage.morphology import skeletonize

im = cv2.imread('a.png',cv2.IMREAD_GRAYSCALE)
#conts ,_ = cv2.findContours(im, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
#cont = conts[0]
sk = skeletonize(im >0)
plt.imshow(sk)

skeleton_yx = np.argwhere(sk > 0)
pts = np.flip(skeleton_yx, axis=None) #yx -> xy
print(pts.shape) #(774, 2)

I want to order these pts to reproduce the sk image above.

EDIT:
following a comment by Michael Szczesny that refers to this question , the first answer by Imanol Luengo works fine but it's very slow at step 4 (Find the path with smallest cost from all sources):

paths = [list(nx.dfs_preorder_nodes(T, i)) for i in range(len(points))]

some other solutions require knowing the first point to start ordering from, I don't know a way to find it, that generalizes to other problems where there's a single, continuous unsorted line (I guess finding extreme points is possible either by looping over all white pixels or by using filters, will try it).

what I've tried:

def rotational_sort(list_of_xy_coords):
  cx, cy = list_of_xy_coords.mean(0)
  x, y = list_of_xy_coords.T
  angles = np.arctan2(x-cx, y-cy)
  indices = np.argsort(angles)
  return list_of_xy_coords[indices]

im = np.zeros_like(im)
pts2= rotational_sort(pts)
cv2.polylines(im, [pts2], False, 255,2)
plt.imshow(im)

from scipy.signal import savgol_filter
from sklearn.decomposition import PCA
from scipy import interpolate
def XYclean(x,y): 
  xy = np.concatenate((x.reshape(-1,1), y.reshape(-1,1)), axis=1)     
  pca = PCA(2)
  pca.fit(xy)
  #transform into pca space   
  xypca = pca.transform(xy) 
  newx = xypca[:,0]
  newy = xypca[:,1]
  indexSort = np.argsort(x)
  newx = newx[indexSort]
  newy = newy[indexSort]
  #add some more points (optional)
  f = interpolate.interp1d(newx, newy, kind='linear')        
  newX=np.linspace(np.min(newx), np.max(newx), 100)
  newY = f(newX)            
  #return back to old coordinates
  xyclean = pca.inverse_transform(np.concatenate((newX.reshape(-1,1), newY.reshape(-1,1)), axis=1) )
  xc=xyclean[:,0]
  yc = xyclean[:,1]
  return np.hstack((xc.reshape(-1,1),yc.reshape(-1,1))).astype(int)

im = np.zeros_like(im)
X,Y=pts[:,0],pts[:,1]
pts2= XYclean(X,Y)
cv2.polylines(im, [pts2], False, 255,2)
plt.imshow(im)