The bwdistgeodesic function can help you with the 'cityblock' or ' quasi-euclidean ' method :

For example, if you choose the first edge point(x1,y1) and you compute bwdistgeodesic(bw,y1,x1,'cityblock'), you will get a distance matrix whose values will be in the range [0 : number of connected components -1]. At the location (x1,y1) of your first point , the distance matrix value will be 0, those of the next pixel will be 1, the next 2 , etc... The order gives priority to 4-connectivity, I think it is that you want. Then you just have to sort the values to get your ordered list.

An example :

selected = bwselect(bw,y1,x1);

mat_dist = bwdistgeodesic(selected,y1,x1,'cityblock');

comp = find(selected);

comp(:,2) = mat_dist(comp(:,1));

ordered_list_ind = sortrows(comp,2);

[ordered_list_sub(:,1) ordered_list_sub(:,2)] = ind2sub(size(selected),ordered_list_ind(:,1));

Remarks : The 'cityblock' method works only for 4-connectivity. If you need this for 8-connectivity, replace the 'cityblock' method by the 'quasi-euclidean' one, sorting the values will allow the same result for 8-c.