Animating Body Segments during Gait - Rotation Help

I am writing a Matlab script to animate the local coordinate system of the thigh and shank during normal walking motion. I have data from a Vicon motion capture system that gives me the coordinates of points in space, which I have used to extrapolate a transformation matrix from. I have created a transformation matrix by creating orthogonal vectors between a set of trackers in close proximity to each other. Below, I have the code to create the vectors, normalize them, and add them to the transformation matrices (for brevity I only added the right thigh and shank):

% Loading data:
load('markers.mat');
RTIP = markers.RTIP;
RTIA = markers.RTIA;
RTSP = markers.RTSP;
RSIP = markers.RSIP;
RSIA = markers.RSIA;
RSSP = markers.RSSP;

% Length of frames
nFrames = length(RTIP);

% Right Thigh
% RTIP is the origin, RTIA and RTSP are the points next to and above it

% Initialize transformation matrix
Td_thigh=[ones(1,1,nFrames),zeros(1,3,nFrames);zeros(3,4,nFrames)];

% getting r1 and r2
r1_thigh=RTIA-RTIP; % non-normalized y-axis
r2_thigh=RTSP-RTIP;
for i=1:nFrames
    r3_thigh(i,:)=cross(r1_thigh(i,:),r2_thigh(i,:)); % non-normalized x-axis
    r4_thigh(i,:)=cross(r3_thigh(i,:),r1_thigh(i,:)); % non-normalized z-axis
end

% normalizing vectors
for i=1:nFrames
    xhat_thigh(i,:)=r3_thigh(i,:)./norm(r3_thigh(i,:));
    yhat_thigh(i,:)=r1_thigh(i,:)./norm(r1_thigh(i,:));
    zhat_thigh(i,:)=r4_thigh(i,:)./norm(r4_thigh(i,:));
end

% for loop to fill out transformation matrix
for i=1:nFrames
    for j=1:3
        Td_thigh(j+1,1,i)=(RTIP(i,j))';
        Td_thigh(j+1,2,i)=(xhat_thigh(i,j))';
        Td_thigh(j+1,3,i)=(yhat_thigh(i,j))';
        Td_thigh(j+1,4,i)=(zhat_thigh(i,j))';
    end
end

%%

% Right Shank
% RSIP is the origin, RSIA and RSSP are the points next to and above it

% Initialize transformation matrix
Td_shank=[ones(1,1,nFrames),zeros(1,3,nFrames);zeros(3,4,nFrames)];

r1_shank=RSIA-RSIP; % non-normalized y-axis
r2_shank=RSSP-RSIP;
for i=1:nFrames
    r3_shank(i,:)=cross(r1_shank(i,:),r2_shank(i,:)); % non-normalized x-axis
    r4_shank(i,:)=cross(r3_shank(i,:),r1_shank(i,:)); % non-normalized z-axis
end

% normalizing vectors
for i=1:nFrames
    xhat_shank(i,:)=r3_shank(i,:)./norm(r3_shank(i,:));
    yhat_shank(i,:)=r1_shank(i,:)./norm(r1_shank(i,:));
    zhat_shank(i,:)=r4_shank(i,:)./norm(r4_shank(i,:));
end

% for loop to fill out transformation matrix
for i=1:nFrames
    for j=1:3
        Td_shank(j+1,1,i)=(RSIP(i,j))';
        Td_shank(j+1,2,i)=(xhat_shank(i,j))';
        Td_shank(j+1,3,i)=(yhat_shank(i,j))';
        Td_shank(j+1,4,i)=(zhat_shank(i,j))';
    end
end

The transformation matrix that I was taught how to use is in this format:

[ 1 0 0 0 ]

T = [ Xorigin cosXx cosXy cosXz ]

[ Yorigin cosYx cosYy cosYz ]

[ Zorigin cosZx cosZy cosZz ]

(the bottom right 3x3 matrix is the rotation matrix)

Then, I would animate my transformation matrices with this code:

% Create 3D figure
figure;
grid on;
hold on;
axis equal;
xlabel('X');
ylabel('Y');
zlabel('Z');

% Set the viewing angle
azimuth = 135; % Horizontal rotation in degrees
elevation = 30; % Vertical elevation in degrees
view(azimuth, elevation); % Apply the viewing angle
rotate3d on; % allow mouse rotation

% Set axis limits
xlim([0, 800]);
ylim([-4000, 3000]);
zlim([-50, 1400]);

% The length of the coordinate system axes
axisLength = 85;

% Number of frames for the animation
numFrames = size(Td_thigh, 3);

% Initialize lines for the X, Y, Z axes of the coordinate systems
hXAxis_thigh = line([0 0], [0 0], [0 0], 'Color', 'r', 'LineWidth', 2); % X-axis in red
hYAxis_thigh = line([0 0], [0 0], [0 0], 'Color', 'g', 'LineWidth', 2); % Y-axis in green
hZAxis_thigh = line([0 0], [0 0], [0 0], 'Color', 'b', 'LineWidth', 2); % Z-axis in blue

hXAxis_shank = line([0 0], [0 0], [0 0], 'Color', 'r', 'LineWidth', 2); % X-axis in red
hYAxis_shank = line([0 0], [0 0], [0 0], 'Color', 'g', 'LineWidth', 2); % Y-axis in green
hZAxis_shank = line([0 0], [0 0], [0 0], 'Color', 'b', 'LineWidth', 2); % Z-axis in blue

% Animation loop
for i = 1:numFrames

    % Thigh
    % Extract the Rotation matrix for the current frame for the thigh
    R_thigh = Td_thigh(2:4, 2:4, i);

    % Extract the origin of the thigh for the current frame
    origin_thigh = Td_thigh(2:4, 1, i);
    
    % Calculate the end points of the coordinate system axes
    xAxisEndPoint_thigh = origin_thigh + axisLength * R_thigh(1,:)';
    yAxisEndPoint_thigh = origin_thigh + axisLength * R_thigh(2,:)';
    zAxisEndPoint_thigh = origin_thigh + axisLength * R_thigh(3,:)';
    
    % Update the X, Y, Z axes of the coordinate system
    set(hXAxis_thigh, 'XData', [origin_thigh(1), xAxisEndPoint_thigh(1)], 'YData', [origin_thigh(2), xAxisEndPoint_thigh(2)], 'ZData', [origin_thigh(3), xAxisEndPoint_thigh(3)]);
    set(hYAxis_thigh, 'XData', [origin_thigh(1), yAxisEndPoint_thigh(1)], 'YData', [origin_thigh(2), yAxisEndPoint_thigh(2)], 'ZData', [origin_thigh(3), yAxisEndPoint_thigh(3)]);
    set(hZAxis_thigh, 'XData', [origin_thigh(1), zAxisEndPoint_thigh(1)], 'YData', [origin_thigh(2), zAxisEndPoint_thigh(2)], 'ZData', [origin_thigh(3), zAxisEndPoint_thigh(3)]);

    % Shank
    % Extract the rotation matrix for the current frame for the shank
    R_shank = Td_shank(2:4, 2:4, i);

    % Extract the origin of the shank for the current frame
    origin_shank = Td_shank(2:4, 1, i);
    
    % Calculate the end points of the coordinate system axes
    xAxisEndPoint_shank = origin_shank + axisLength * R_shank(1,:)';
    yAxisEndPoint_shank = origin_shank + axisLength * R_shank(2,:)';
    zAxisEndPoint_shank = origin_shank + axisLength * R_shank(3,:)';
    
    % Update the X, Y, Z axes of the coordinate system
    set(hXAxis_shank, 'XData', [origin_shank(1), xAxisEndPoint_shank(1)], 'YData', [origin_shank(2), xAxisEndPoint_shank(2)], 'ZData', [origin_shank(3), xAxisEndPoint_shank(3)]);
    set(hYAxis_shank, 'XData', [origin_shank(1), yAxisEndPoint_shank(1)], 'YData', [origin_shank(2), yAxisEndPoint_shank(2)], 'ZData', [origin_shank(3), yAxisEndPoint_shank(3)]);
    set(hZAxis_shank, 'XData', [origin_shank(1), zAxisEndPoint_shank(1)], 'YData', [origin_shank(2), zAxisEndPoint_shank(2)], 'ZData', [origin_shank(3), zAxisEndPoint_shank(3)]);
    
    % Pause briefly to create an animation effect
    drawnow; % Update the figure window
    pause(0.04);
end

When I run this code, however, my local coordinate systems on my thigh and shank seem to rotate the wrong way. The thigh rotates clockwise when it should rotate counter-clockwise throughout my animation, same with the shank (I'm talking about rotations about the x-axis, which points sideways). I'm completely stumped on where I went wrong. Any help would be greatly appeciated.

[Ezekiel Bibbo]

I have actually partially figured out the answer to my own question.

If I negate the y-axis, like so:

yhat_thigh(i,:)=-r1_thigh(i,:)./norm(r1_thigh(i,:));

and negate the direction the y vector is pointing, like so:

yAxisEndPoint_thigh = origin_thigh + -axisLength * R_thigh(2,:)';

(for both the thigh and the shank) the code runs properly, and the rotations fix themselves...

Although I do not know why this works. Could someone smarter than me explain the reasoning behind this? I still feel as though I made an error.

I attached my updated code for your perusal.