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