Showing a fitting curve through centroids of white regions in a binary image

I have a binary image where I find the centroids of the white regions, I plotted the centroids on top of the binary image I need to fit a polynomial (Cubic) through the points and show the fitting line on the binary image. what I get is a flipped line (the blue line in the attached image). How can I show the line correctly on the image?

close all
clear
% Image preprocessing
I= imread('Lbend_seg.jpg'); % Read
Gr=rgb2gray(I); %convert to grayscale
BW = imbinarize(Gr,'adaptive','ForegroundPolarity','dark','Sensitivity',0.4); %Convert to binary
BW = imclearborder(BW);%Remove the background
BW = bwareaopen(BW, 50);%Remove the regions which have less than 50 pixels connected
figure 
imshow(BW);
hold on
Kmedian = medfilt2(BW); % denoising the binary image
Kmedian = medfilt2(Kmedian);
[B,L] = bwboundaries(Kmedian,'noholes'); % Detecting boundaries

%Plot the boundaries on top of the binary image
for k = 1:length(B)
   boundary = B{k};
   plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 1.5)
end


%Approximating the boundaries as a polygon
B_tr= cellfun(@transpose,B,'UniformOutput',false); %Transposing each cell inside cell array

for n=1:max(size(B_tr))
     temp=cell2mat(B_tr(n));
     
     reduced{n}=reduce_poly_PR(temp,6);
end
poly_rd= cellfun(@transpose,reduced,'UniformOutput',false); %Transposing each cell inside cell array
poly_rd=poly_rd'; %trasposing the cell

%plotting the polygons on top of the BW image
figure
imshow(BW)
hold on
 for k = 1:length(poly_rd)
   boundary = poly_rd{k};
   boundariesadded=plot(boundary(:,2), boundary(:,1), 'r', 'LineWidth', 2);
 end
 
%Find the centroid of each vertebrae
stats = regionprops('table',BW,'Orientation', 'MajorAxisLength', ...
    'MinorAxisLength', 'Eccentricity', 'Centroid');

%plot the centroid of each vertebrae on top of the polygons
plot(stats.Centroid(:,1),stats.Centroid(:,2),'o')

% Plot the centroids on the blnariy image
figure
imshow(BW)
hold on
plot(stats.Centroid(:,1),stats.Centroid(:,2),'.')
 coefficients=polyfit(stats.Centroid(:,2),stats.Centroid(:,1), 3);
 fittedX = linspace(min(stats.Centroid(:,2)), max(stats.Centroid(:,2)), 500);
fittedY = polyval(coefficients, fittedX);
plot(fittedX, fittedY, 'b-', 'linewidth', 4);
% Visualize the orientation of the ellipse as a sanitary check
% Orientation is the angle between the minor axis and the horizontal
s = regionprops(BW, 'Orientation', 'MajorAxisLength', ...
    'MinorAxisLength', 'Eccentricity', 'Centroid');
figure
imshow(BW)
hold on

phi = linspace(0,2*pi,50);
cosphi = cos(phi);
sinphi = sin(phi);

for k = 1:length(s)
    xbar = s(k).Centroid(1);
    ybar = s(k).Centroid(2);

    a = s(k).MajorAxisLength/2;
    b = s(k).MinorAxisLength/2;

    theta = pi*s(k).Orientation/180;
    R = [ cos(theta)   sin(theta)
         -sin(theta)   cos(theta)];

    xy = [a*cosphi; b*sinphi];
    xy = R*xy;

    x = xy(1,:) + xbar;
    y = xy(2,:) + ybar;

    plot(x,y,'r','LineWidth',2);
    
end
hold off

Answer 1

[Image Analyst]

Wow, you sure do know how to complicate things. It could be a lot simpler. Anyway, you swapped x and y. Try putting them in this order:

plot(fittedY, fittedX, 'b-', 'linewidth', 4);

[Osi fri]

Haha thank you very much Image Analyst, It worked perfectly!