The principle axes you can find with svd(),
XY=[X(:),Y(:)]; %the scattered data
mu=mean(XY);
P=svd(XY-mu,0); %P contains the principle axes
Once you've done the pca, you can rotate/translate your scattered (X,Y) data so that the principle axes are aligned with the coordinate axes and so that the data is centered at the origin.
xy=P.'*(XY-mu); %rotated data
Then, the problem boils down to estimating the dimensions a and b of a non-rotated rectangle. The latter, you might be able to do with fminsearch:
xp=abs(xy(:,1)); yp=abs(xy(:,2));
ab=fminsearch(@(ab) objective(ab,xp,yp), [max(xp),max(yp)] );
function fval = objective(ab,x, y)
a=ab(1)/2; b=ab(2)/2; %half length and width of rectangle
%%% Separate points into 4 regions (positive quadrant only)
in= x<a & y<b; %inside rectangle
out1= x>=a & y<b;
out2= x>=a & y>=b;
out3= x<a & y>=b;
%%% Distances of inner points to boundary
xin=x(in); y=y(in); %points inside
Din=sum( min([abs(xin-a), abs(xin+a), abs(yin-b), abs(yin+b)],[],2) );
%%% Distances of outer points to boundary
Dout1 = sum(x(out1)-a);
Dout2 = sum(vecnorm([x(out2),y(out2)] -[a,b],2,2));
Dout3 = sum(y(out3)-b);
fval=Din+Dout1+Dout2+Dout3;
end
