function [I_approx, I_approx_xy, csa_approx, number_triangles, radii_in, radii_out, Cx, Cy] = triangle_approximation_function(x_out, y_out, x_in, y_in, tri_plot) %% Approximating CSA and area moment of inertia using triangles % Calculates the cross sectional area, area moment of inertia about x and y axes and product % moment of area for a shape with an inside and outside edge % Splits the shape into triangles and calculates properties based on these % triangles % CSA is total area of triangles % Area moment of inertia combined area moment of inertia of all triangles % Inputs to function are % - Vectors containing: x coordinates of outside edge % y coordinates of outside edge % x coordinates of inside edge % y coordinates of inside edge % - tri_plot (optional): 1 to include plot of coordinates of triangles % Outputs of function % - Vector of area moment of inertia approximation [Ix, Iy] % - Product moment of area approximation % - Cross sectional area approximation % - Number of triangles used in approximation % - Distance from centre for each point used in approximation - referred to % as radius of each point % - Ratios of radii used in approximation %% Ensure same number of points on inside and outside edges % As required by method of creating triangles diff = length(x_out)-length(x_in); if diff ~= 0 ratio = round(length(x_out)/length(x_in)); % Find ratio of number of outside points to number of inside points x_out = x_out(1:ratio:end); % Take every n of outside points, where n=the ratio y_out = y_out(1:ratio:end); % Take every n of outside points, where n=the ratio diff = length(x_out)-length(x_in); for i = 1:diff % Still any difference then remove every 3rd point until the number of points are equal x_out(3*i) = []; y_out(3*i) = []; end end % figure() % plot(x_out, y_out, 'r'); hold on; % plot(x_in, y_in, 'b'); hold off; %% Create triangles % Create triangles starting from each point around the edge of the bone % Use a given point, the point opposite it and the point after it % Create array of coordinates of triangle vertices tri_coords = zeros(3,2,2*length(x_out)); for i = 1:2:2*length(x_out) %inside triangles if i ~= (2*length(x_out))-1 tri_coords(1,1,i) = x_in((i+1)/2); % x coordinates tri_coords(2,1,i) = x_in(((i+1)/2) +1); tri_coords(3,1,i) = x_out((i+1)/2); tri_coords(1,2,i) = y_in((i+1)/2); % y coordinates tri_coords(2,2,i) = y_in(((i+1)/2) +1); tri_coords(3,2,i) = y_out((i+1)/2); else % last triangle includes points numbered 1 tri_coords(1,1,i) = x_in((i+1)/2); % x coordinates tri_coords(2,1,i) = x_in(1); tri_coords(3,1,i) = x_out((i+1)/2); tri_coords(1,2,i) = y_in((i+1)/2); % y coordinates tri_coords(2,2,i) = y_in(1); tri_coords(3,2,i) = y_out((i+1)/2); end end for j = 2:2:2*length(x_out) %outside triangles if j~=2*length(x_out) tri_coords(1,1,j) = x_out(j/2); % x coordinates tri_coords(2,1,j) = x_out((j/2)+1); tri_coords(3,1,j) = x_in((j/2)+1); tri_coords(1,2,j) = y_out(j/2); % y coordinates tri_coords(2,2,j) = y_out((j/2)+1); tri_coords(3,2,j) = y_in((j/2)+1); else % last triangle includes points numbered 1 tri_coords(1,1,j) = x_out(j/2); % x coordinates tri_coords(2,1,j) = x_out(1); tri_coords(3,1,j) = x_in(1); tri_coords(1,2,j) = y_out(j/2); % y coordinates tri_coords(2,2,j) = y_out(1); tri_coords(3,2,j) = y_in(1); end end % Plot to check if ~exist('tri_plot') tri_plot=0; end if tri_plot==1 figure() for i = 1:size(tri_coords, 3) plot([tri_coords(1,1,i) tri_coords(2,1,i)], [tri_coords(1,2,i) tri_coords(2,2,i)]); hold on; plot([tri_coords(1,1,i) tri_coords(3,1,i)], [tri_coords(1,2,i) tri_coords(3,2,i)]); plot([tri_coords(2,1,i) tri_coords(3,1,i)], [tri_coords(2,2,i) tri_coords(3,2,i)]); end end %% CSA approximation % Herons law % area = sqrt( s(s-a)(s-b)(s-c) ) % s = (a+b+c)/2 % a,b,c side lengths tri_lengths = zeros(3,1,2*length(x_out)); % Create array of side lengths of triangles for i = 1:2*length(x_out) tri_lengths(1,1,i) = sqrt( (tri_coords(1,1,i)-tri_coords(2,1,i))^2 + (tri_coords(1,2,i)-tri_coords(2,2,i))^2 ); tri_lengths(2,1,i) = sqrt( (tri_coords(1,1,i)-tri_coords(3,1,i))^2 + (tri_coords(1,2,i)-tri_coords(3,2,i))^2 ); tri_lengths(3,1,i) = sqrt( (tri_coords(2,1,i)-tri_coords(3,1,i))^2 + (tri_coords(2,2,i)-tri_coords(3,2,i))^2 ); end tri_s = zeros(1,2*length(x_out)); % Create array of half perimeters for i = 1:2*length(x_out) tri_s(1,i) = (tri_lengths(1,1,i)+tri_lengths(2,1,i)+tri_lengths(3,1,i))/2; end tri_area = zeros(1,2*length(x_out)); % Create array of areas of triangles for i = 1:2*length(x_out) tri_area(i) = sqrt( tri_s(i) * (tri_s(i)-tri_lengths(1,1,i)) * (tri_s(i)-tri_lengths(2,1,i)) *(tri_s(i)-tri_lengths(3,1,i))); end % Calculate area tri_area = abs(tri_area); number_triangles = length(tri_area); csa_approx = sum(tri_area); %% Area moment of inertia approximation % Find centroid of each triangle - mean of vertices % Find I for each triangle about x and y axis % Use parallel axis theorem % - Ix = Ix,cm + dx^2*A % - Iy = Iy,cm + dy^2*A % Approximate to I = d^2*A % Distances are coordinates of centroids tri_cent = zeros(2*length(x_out), 2); % Create array of centroids for a = 1:2*length(x_out) tri_cent(a,1) = (tri_coords(1,1,a)+tri_coords(2,1,a)+tri_coords(3,1,a))/3; tri_cent(a,2) = (tri_coords(1,2,a)+tri_coords(2,2,a)+tri_coords(3,2,a))/3; end % Plot to check centroids % figure() % for i = 1:size(tri_coords, 3) %for each triangle % plot([tri_coords(1,1,i) tri_coords(2,1,i)], [tri_coords(1,2,i) tri_coords(2,2,i)]); hold on; % plot([tri_coords(1,1,i) tri_coords(3,1,i)], [tri_coords(1,2,i) tri_coords(3,2,i)]); % plot([tri_coords(2,1,i) tri_coords(3,1,i)], [tri_coords(2,2,i) tri_coords(3,2,i)]); % end % scatter(tri_cent(:,1), tri_cent(:,2)); hold off; tri_cent_sq = tri_cent.^2; tri_I_axis = zeros(2*length(x_out), 2); tri_I_axis(:,1) = tri_area'.*tri_cent_sq(:,2); tri_I_axis(:,2) = tri_area'.*tri_cent_sq(:,1); I_approx(1,1) = sum(tri_I_axis(:,1)); I_approx(1,2) = sum(tri_I_axis(:,2)); %% Product moment of area % Product moment of area including parallel axis theorem % Ixy = Ixy,cm + dx*dy*A % use same approximation as with area moment of inertia % Ixy ~= dx*dy*A tri_I_xy = zeros(2*length(x_out),1); for i = 1:length(tri_I_xy) tri_I_xy(i) = tri_area(i)*tri_cent(i,1)*tri_cent(i,2); end I_approx_xy = sum(tri_I_xy); %% Centroid of tibia % For the calculations to be accurate the centroid should be at 0,0 % Previous centroid based on centre of inside edge of cortical bone % New centroid based on centroid of entire bone within segment slice % Split tibia into i triangles % Cx - x coordinate of centroid, Cy - y coordinate of centroid % Cx = sum(Cix*Ai)/sum(Ai) % Cy = sum(Ciy*Ai)/sum(Ai) Cx_num = 0; Cy_num = 0; for i = 1:2*length(x_out) if tri_area(i)==0 continue; end Cx_num = Cx_num + (tri_cent(i,1)*tri_area(i)); Cy_num = Cy_num + (tri_cent(i,2)*tri_area(i)); end Cx = Cx_num/csa_approx; Cy = Cy_num/csa_approx; %% Radius values for tibia approximation radii_out = zeros(1,length(x_out)); radii_in = zeros(1,length(x_out)); for i = 1:length(x_out) radii_out(1,i) = sqrt( x_out(i)^2 + y_out(i)^2 ); % Outer radius radii_in(1,i) = sqrt( x_in(i)^2 + y_in(i)^2 ); % Inner radius end