Here is the MATLAB Code. Basically I'm applying the calculations to four corners of a square which is being transformed.
clear all
close all
clc
dt = 1;
Scale = 3; % Just defined to care of the plotting stuff
Sx = [0 0 1 1 0];
Sy = [0 1 1 0 0];
Sz = [0 0 0 0 0];
P1 = [];
P2 = [];
P3 = [];
P4 = [];
K = [];
X = [];
Y = [];
Z = [];
Radius = [];
Radius_PB = [];
Tprev = 0;
figure
plot(Sx,Sy)
axis(2*[-10 10 -10 10]*Scale);
grid on;
P = 2*Scale*[Sx; Sy; Sz;ones(size(Sx))];
Pprev = P;
for t = 1:dt:80
% R = makehgtform('zrotate',deg2rad(t),'translate',[(t/40) (t/40) 0]);
R = makehgtform('zrotate',deg2rad(t),'translate',[(t/50)^3 (t/50) 0]); % Homogeneous Transformation matrix
% R = makehgtform('zrotate',deg2rad(t),'translate',[2 4 0]);
% R = makehgtform('zrotate',deg2rad(t));
Pnew = R*P; % Evaluation of transformed points
P1 = [P1 Pnew(1:2,1)];
P2 = [P2 Pnew(1:2,2)];
P3 = [P3 Pnew(1:2,3)];
P4 = [P4 Pnew(1:2,4)];
PC = [P1 P2 P3 P4];
%Radius calculation using Frenet-Serret formulas
dr = (Pnew-Pprev); %Calculating dr
Tnew = [dr(:,1)/(sum(dr(:,1).^2))^0.5,...
dr(:,2)/(sum(dr(:,2).^2))^0.5,...
dr(:,3)/(sum(dr(:,3).^2))^0.5,...
dr(:,4)/(sum(dr(:,4).^2))^0.5,...
dr(:,5)/(sum(dr(:,5).^2))^0.5]; %Calculating dr/ds
dT = (Tnew-Tprev);
dTds = [dT(:,1)/(sum(dr(:,1).^2))^0.5,...
dT(:,2)/(sum(dr(:,2).^2))^0.5,...
dT(:,3)/(sum(dr(:,3).^2))^0.5,...
dT(:,4)/(sum(dr(:,4).^2))^0.5,...
dT(:,5)/(sum(dr(:,5).^2))^0.5];
k = [(sum(dTds(:,1).^2))^0.5,...
(sum(dTds(:,2).^2))^0.5,...
(sum(dTds(:,3).^2))^0.5,...
(sum(dTds(:,4).^2))^0.5]
radius = [1/k(1),1/k(2),1/k(3),1/k(4)]
K = [K; k];
Radius = [Radius; radius];
%Radius calculation using Perperndicular bisector method
drxyz = [dr(1:3,1), dr(1:3,2), dr(1:3,3), dr(1:3,4)];
Z = repmat([0; 0; 1],1,4);
Nr = cross(drxyz,Z);
Nr = [Nr(1:2,1)/sum(Nr(1:2,1).^2)^0.5, -Nr(1:2,2)/sum(Nr(1:2,2).^2)^0.5,...
Nr(1:2,3)/sum(Nr(1:2,3).^2)^0.5, -Nr(1:2,4)/sum(Nr(1:2,4).^2)^0.5];
MP = (Pnew+Pprev)/2;
DMP12 = [MP(1:2,2)-MP(1:2,1)];
DMP34 = [MP(1:2,4)-MP(1:2,3)];
lambda = [inv(Nr(:,1:2))*DMP12, inv(Nr(:,3:4))*DMP34]
r1 = MP(1:2,1)+lambda(1)*Nr(:,1)-Pnew(1:2,1); %Calculating the radius for 1st corner
r2 = MP(1:2,2)+lambda(2)*Nr(:,2)-Pnew(1:2,2);
r3 = MP(1:2,3)+lambda(3)*Nr(:,3)-Pnew(1:2,3);
r4 = MP(1:2,4)+lambda(4)*Nr(:,4)-Pnew(1:2,4);
Radius_PB = [Radius_PB; sum(r1.^2)^0.5, sum(r2.^2)^0.5,sum(r3.^2)^0.5, sum(r4.^2)^0.5 ]
%Plotting and annotating
plot(Pnew(1,:),Pnew(2,:),'-k',P1(1,:),P1(2,:),'-r',P2(1,:),P2(2,:),'-b',P3(1,:),...
P3(2,:),'-g',P4(1,:),P4(2,:),'-y');
text(Pnew(1,1),Pnew(2,1),'1');
text(Pnew(1,2),Pnew(2,2),'2');
text(Pnew(1,3),Pnew(2,3),'3');
text(Pnew(1,4),Pnew(2,4),'4');
axis(2*[-10 10 -10 10]*Scale);
grid on;
pause(0.01);
Pprev = Pnew;
Tprev = Tnew;
% pause;
end
figure
s = length(K)
h = plot(1:s,K(1:s,1),'-b',1:s,K(1:s,2),'-r',1:s,K(1:s,3),'-g',1:s,K(1:s,4),'-y');
axis tight
xlabel('Sample No \rightarrow','fontsize',12)
ylabel('\kappa \rightarrow','fontsize',12);
legend(h,'Vertex 1','Vertex 2','Vertex 3','Vertex 4');
grid on
figure
h = plot(1:s,Radius(1:s,1),'-b',1:s,Radius(1:s,2),'-r',1:s,Radius(1:s,3),'-g',1:s,Radius(1:s,4),'-y');
axis tight
xlabel('Sample No \rightarrow','fontsize',12)
ylabel('Radius of Curvature \rightarrow','fontsize',12);
legend(h,'Vertex 1','Vertex 2','Vertex 3','Vertex 4');
grid on
figure
h = plot(1:s,Radius_PB(1:s,1),'-b',1:s,Radius_PB(1:s,2),'-r',1:s,Radius_PB(1:s,3),'-g',1:s,Radius_PB(1:s,4),'-y');
axis tight
Title(' Radius of Curvature from perpendicular bisector','fontsize',12);
xlabel('Sample No \rightarrow','fontsize',12)
ylabel('Radius of Curvature \rightarrow','fontsize',12);
legend(h,'Vertex 1','Vertex 2','Vertex 3','Vertex 4');
grid on
