%% Bone digitisation and geometric properties
% Megan Kenny, 2019, 
% Edited by Matthew Ellison and Hannah Rice
% University of Exeter, UK


% Code to digitise MR image of bone and calculate properties of shape -
% cross sectional area, area moment of inertia, product moment of area
% Outputs include area moment of inertia approximation [Ix, Iy], product moment of area approximation,...
% cross sectional area approximation 

% Import jpeg of MR image of bone cropped to the region of interest
% Basic image processing and edge detection to digitise bone
% Align bone so that x axis is from medial to lateral, y axis is    
% from posterior to anterior                            
% Find coordinates of ant/post/med/lat points  
% Calculate properties of digitised shape in m, m^2, m^4


% Code to rotate into a local coordinate system has been developed
% alongside this, which may be relevant and necessary
                            
clear all
close all
%% Import image

filename = 'Sample_image';   %jpeg of MR image cropped to region of interest

resolution = 0.676*10^-3;               % Resolution of scanner - in metres  i.e. 1 pixel = 0.676mm

tibia = imread([filename, '.jpg']); 

% figure(); imshow(tibia); figure(); imshow(ankle);

tibia_mod = tibia;

%% Edge detection - tibia
% Modify image to help with edge detection - effectively increase the contrast
% Make the denser cortical bone denser - lower pixel value
% Make the muscle lighter
% Make trabecular bone lighter
% Values may need to be adjusted if edge detection not ideal

for i = 1:size(tibia_mod,1)
    for j = 1:size(tibia_mod,2)
        if tibia_mod(i,j)<10        % Modify if necessary
            tibia_mod(i,j) = 0;
        elseif tibia_mod(i,j)>25    % Modify if necessary
            tibia_mod(i,j) = tibia_mod(i,j)*3;
        elseif tibia_mod(i,j)>70    % Modify if necessary
            tibia_mod(i,j) = tibia_mod(i,j)*1.5;
        end 
    end 
end 
% figure(); imshow(tibia_mod)

tibia_edges = edge(tibia_mod(:,:,1), 'canny');     % Edge detection

% ME edit - Remove unwanted edges - Only want the edges that define the inside and outside edge of the cortical bone
% - finds the number of connected regions and removes the smallest until only 2 remain
CC = bwconncomp(tibia_edges); 
while CC.NumObjects >2 
    numPixels = cellfun(@numel,CC.PixelIdxList);
    [~,idx] = min(numPixels);
    tibia_edges(CC.PixelIdxList{idx}) = 0;
    CC = bwconncomp(tibia_edges);
end

%% Confirm the edges that have been detected

figure();
subplot(1,2,1)
imshow(tibia_edges)

answer=questdlg('Digitised ok?','Yes','No');
switch answer
    case 'Yes'
        ;
    case 'No'
        
%% HR edit
        %Identify remaining unwanted edges 
        title('Select unwanted edges')
        [xi,yi] = getpts; 
        % Click the points you want to exclude
        % Double-click/ Shift + click/ Right-click for final point or press enter
        % when finished
        %  Delete or Backspace to remove a point
        yi=round(yi); 
        xi=round(xi); 
            for i=1:length(xi)
            tibia_edges(yi(i),xi(i))=0;
            end
            figure();subplot(1,2,1)
            imshow(tibia_edges)
end

answer=questdlg('Unwanted edges removed?','Yes','No');
switch answer
    case 'Yes'
        ;
    case 'No'       
           title('Select unwanted edges')
        [xi,yi] = getpts; 
        yi=round(yi); 
        xi=round(xi); 
            for i=1:length(xi)
            tibia_edges(yi(i),xi(i))=0;
            end
            figure();subplot(1,2,1)
            imshow(tibia_edges)
end
        %Identify missing edges
        title('Select missing edges')
        [xj,yj] = getpts; 
        yj=round(yj); 
        xj=round(xj); 
        for j=1:length(xj)
        tibia_edges(yj(j),xj(j))=1;
        end
        close all
        figure();subplot(1,2,1)
        imshow(tibia_edges)

% It may be necessary to manually adjust variable tibia-edges to include/exclude all required pixels 


%% Find centroid estimate for tibia

s=regionprops(tibia_edges,'centroid');
centroids = cat(1, s.Centroid);
centroid = [centroids(2,1), centroids(2,2)];    % Take centroid of inside edge of cortical bone as initial approximation of centroid 
centroid_rd = floor(centroid);
% % Plot to check
% figure(); imshow(tibia_edges); hold on ;scatter(centroid(1), centroid(2));

%% ME - Threshold by creating separate images
% Split the bone into two separate images for the inside and outside edge of tibia
CC = bwconncomp(tibia_edges);
bone_edge_outer=tibia_edges;
bone_edge_inner=tibia_edges;

numPixels = cellfun(@numel,CC.PixelIdxList);
[smallest,idx] = min(numPixels);
bone_edge_outer(CC.PixelIdxList{idx}) = 0;
[largest,idx] = max(numPixels);
bone_edge_inner(CC.PixelIdxList{idx}) = 0;

%% Get coordinates of pixels around edge
% Record the row and column number for every cell of array that has a 1 in

count = 1;
for i = 1:size(bone_edge_outer,1)       %iterate along rows
    for j = 1:size(bone_edge_outer,2)   %iterate along columns
        if bone_edge_outer(i,j) == 1
            x_out(count) = i;
            y_out(count) = j;
            count = count+1;
        end
    end
end

count = 1;
for i = 1:size(bone_edge_inner,1)       %iterate along rows
    for j = 1:size(bone_edge_inner,2)   %iterate along columns
        if bone_edge_inner(i,j) == 1
            x_in(count) = i;
            y_in(count) = j;
            count = count+1;
        end
    end
end

%% Sort variables
% Currently every cell is referenced from 1,1, want to reference from centroid
% - so subtract centroid from each
x_out = x_out - centroid_rd(2);
y_out = y_out - centroid_rd(1);
x_in = x_in - centroid_rd(2);
y_in = y_in - centroid_rd(1);

figure(); scatter(x_in,y_in);hold on; scatter(x_out, y_out);

% Rotate coordinates to align with MR image
angle = -pi/2;                                          % Angle in radians
R = [cos(angle), -sin(angle); sin(angle), cos(angle)];  % Counterclockwise rotation matrix with negative angle means clockwise rotation
out = [x_out; y_out];
in = [x_in; y_in];
out_rot = R*out;
in_rot = R*in;

x_in = in_rot(1,:);
y_in = in_rot(2,:);
x_out = out_rot(1,:);
y_out = out_rot(2,:);

% % Plot to check
figure(); scatter(x_in,y_in);hold on; scatter(x_out, y_out);


%% Find angle of points to sort
% Method of splitting bone into triangles requires points to be in order by
% angle made with the x axis

angle_out = atand(y_out./x_out);
angle_in = atand(y_in./x_in);

% Make it so that angles range from 0-359 deg to make easier to sort
for i = 1:length(x_out)                          % Outside points
    if x_out(i) > 0 && y_out(i) >= 0 
                                                 %first quadrant angle 0-90
    elseif x_out(i) <= 0 && y_out(i) > 0 
        angle_out(i) = (90+ angle_out(i)) + 90;  %second quadrant angle 90-180
    elseif x_out(i) < 0 && y_out(i) <= 0 
        angle_out(i) = angle_out(i) + 180;       %third quadrant angle 0-90
    elseif x_out(i) >= 0 && y_out(i) < 0 
        angle_out(i) = (90+angle_out(i)) + 270;  %fourth quadrant angle 0-90
    end 
end 

for i = 1:length(x_in)                           % Inside points
    if x_in(i) > 0 && y_in(i) >= 0 
                                                 %first quadrant angle 0-90
    elseif x_in(i) <= 0 && y_in(i) > 0 
        angle_in(i) = (90+angle_in(i)) + 90;     %second quadrant angle 90-180
    elseif x_in(i) < 0 && y_in(i) <= 0 
        angle_in(i) = angle_in(i) + 180;         %third quadrant angle 0-90
    elseif x_in(i) >= 0 && y_in(i) < 0        
        angle_in(i) = (90+angle_in(i)) + 270;    %fourth quadrant angle 0-90
    end 
end 
     
%% Sort angles
% Sort function at bottom of script

x_out_sorted = [];
y_out_sorted = [];
x_in_sorted = [];
y_in_sorted = [];
while isempty(angle_out) == 0        % Sort outside points
    [x_out_sorted, y_out_sorted, angle_out, x_out, y_out] = find_smallest(angle_out, x_out_sorted, y_out_sorted, x_out, y_out);
end

while isempty(angle_in) == 0         % Sort inside points
    [x_in_sorted, y_in_sorted, angle_in, x_in, y_in] = find_smallest(angle_in, x_in_sorted, y_in_sorted, x_in, y_in);
end

% Plot to check
% figure()
% subplot(1,2,1)
% plot(x_out_sorted, y_out_sorted);
% subplot(1,2,2)
% plot(x_in_sorted, y_in_sorted);

%% Create more points
% Make the estimate more accurate by having more points on circumference to
% create more triangles for approximation

x_out_sorted_2 = zeros(1,2*length(x_out_sorted));      % Create array double length of previous one
x_in_sorted_2 = zeros(1,2*length(x_in_sorted));
y_out_sorted_2 = zeros(1,2*length(x_out_sorted));
y_in_sorted_2 = zeros(1,2*length(x_in_sorted));

for m = 1:2:2*length(x_out_sorted)-2
    x_out_sorted_2(m) = x_out_sorted((m+1)/2);
    x_out_sorted_2(m+1) = (x_out_sorted((m+1)/2)+x_out_sorted(((m+1)/2)+1))/2;
    y_out_sorted_2(m) = y_out_sorted((m+1)/2);
    y_out_sorted_2(m+1) = (y_out_sorted((m+1)/2)+y_out_sorted(((m+1)/2)+1))/2;    
end                % Create point in between each pair of points by averaging both points

for m = 1:2:2*length(x_in_sorted)-2
    x_in_sorted_2(m) = x_in_sorted((m+1)/2);
    x_in_sorted_2(m+1) = (x_in_sorted((m+1)/2)+x_in_sorted(((m+1)/2)+1))/2;
    y_in_sorted_2(m) = y_in_sorted((m+1)/2);
    y_in_sorted_2(m+1) = (y_in_sorted((m+1)/2)+y_in_sorted(((m+1)/2)+1))/2 ;   
end                 % Create point in between each pair of points by averaging both points

% Have to amend method for last points as connect to first points
x_out_sorted_2(length(x_out_sorted_2)-1) = x_out_sorted(length(x_out_sorted));
x_in_sorted_2(length(x_in_sorted_2)-1) = x_in_sorted(length(x_in_sorted));
y_out_sorted_2(length(y_out_sorted_2)-1) = y_out_sorted(length(y_out_sorted));
y_in_sorted_2(length(y_in_sorted_2)-1) = y_in_sorted(length(y_in_sorted));
x_out_sorted_2(length(x_out_sorted_2)) = (x_out_sorted(length(x_out_sorted))+x_out_sorted(1))/2;
x_in_sorted_2(length(x_in_sorted_2)) = (x_in_sorted(length(x_in_sorted))+x_in_sorted(1))/2;
y_out_sorted_2(length(y_out_sorted_2)) = (y_out_sorted(length(y_out_sorted))+y_out_sorted(1))/2;
y_in_sorted_2(length(y_in_sorted_2)) = (y_in_sorted(length(y_in_sorted))+y_in_sorted(1))/2;

% % Plot to check
figure()
subplot(1,2,1)
scatter(x_out_sorted, y_out_sorted); hold on;
scatter(x_in_sorted, y_in_sorted);
subplot(1,2,2)
scatter(x_out_sorted_2, y_out_sorted_2); hold on;
scatter(x_in_sorted_2, y_in_sorted_2);

%% Find coordinates of anterior, posterior, medial, lateral points - specifc to beam bending
% Posterior = (0, y)
% Anterior = (0, y)
% Medial = (x, 0)
% Lateral = (x, 0)

% Find values on axes for outer radii 
y_vals = find(abs(x_out_sorted_2)<=0.5); % y vals when x=0, in case no point exists at 0, there should be point at +/-0.5
x_vals = find(abs(y_out_sorted_2)<=0.5); % x vals when y=0, in case no point exists at 0, there should be point at +/-0.5

ant = y_out_sorted_2(y_vals(1));         % First y value come to is anterior point 
post = y_out_sorted_2(y_vals(end));      % Last y value come to is posterior point

% For medial and lateral if coordinate is positive then lateral, if
% negative then medial

x_points = x_out_sorted_2(x_vals);

count_med = 1;
count_lat = 1;

for i = 1:length(x_points)
    if x_points(i) >= 0    % Lateral
        lat_points(count_lat) = x_points(i);
        count_lat = count_lat+1;
    elseif x_points(i) < 0
        med_points(count_med) = x_points(i);
        count_med = count_med+1;
    end 
end 

med = mode(med_points);
lat = mode(lat_points);

% Plot to check
figure()
plot(x_out_sorted_2, y_out_sorted_2, 'r'); hold on;
scatter(0, post, 'r');
scatter(0, ant, 'b');
scatter(med, 0, 'g');
scatter(lat, 0, 'k');


%% Convert from pixels to m

x_out_m = resolution*x_out_sorted_2;
y_out_m = resolution*y_out_sorted_2;
x_in_m = resolution*x_in_sorted_2;
y_in_m = resolution*y_in_sorted_2;

med = med*resolution;
lat = lat*resolution;
post = post*resolution;
ant = ant*resolution;


%% Run triangle approximation first to find accurate centroid location 

[I_approx_1, I_approx_xy_1, csa_approx_1, number_triangles_tibia_1, radii_in_1, radii_out_1,...
    Cx, Cy] = triangle_approximation_function(x_out_m, y_out_m, x_in_m, y_in_m);

%% Update centroid location 
% First estimate based on the centre of the inside edge of the cortical bone
% New centroid estimate based on the triangles
% Recentre the image on the new calculated centroid and re-do the triangle approximation
    x_out_m = x_out_m - Cx;
    y_out_m = y_out_m - Cy;
    x_in_m = x_in_m - Cx;
    y_in_m = y_in_m - Cy;
    figure; plot(x_out_m, y_out_m); hold on; plot(x_in_m, y_in_m);
    
    [I_approx, I_approx_xy, csa_approx, number_triangles_tibia, radii_in, radii_out,...
        ] = triangle_approximation_function(x_out_m, y_out_m, x_in_m, y_in_m,1);
    

%% Radii ratio for approximation

% In the approximation for area moment of inertia, I = Icm + d^2*A is approximated by I = d^2*A
% For this to be a reasonable approximation the first term Icm has to be
% negligible compared to d^2*A
% Icm will be smaller compared to d^2*A if the triangles are small (more triangles)
% and the bone is distributed far from the centroid (the ratio r_in/(r_out-r_in) is small and the ratio r_out/r_in is large)
% Where radius is used to describe the distance to the centroid of each point digitising
% the edge of the bone
% Preliminary investigation for circular sections suggests values of:
%  r_in/(r_out-r_in) > 1 
%  r_out/r_in < 2
%  number of triangles > 100

radii_ratio = radii_out./radii_in;                      % Ratio r_out/r_in
radii_ratio_average = mean(radii_ratio);                % Average ratio r_out/r_in - ideally < 2
radii_ratio_max = max(radii_ratio);                     % Maximum ratio r_out/r_in - ideally < 2

ratio_width_inner = radii_in./(radii_out-radii_in);     % Ratio r_in/(r_out-r_in)
ratio_width_inner_min = min(ratio_width_inner);         % Minimum ratio r_in/(r_out-r_in) - ideally > 1
ratio_width_inner_average = mean(ratio_width_inner);    % Average ratio r_in/(r_out-r_in) - ideally > 1

% figure()
% subplot(1,2,1)
% plot(radii_ratio)
% legend('Radii ratio - r_out/r_in')
% subplot(1,2,2)
% plot(ratio_width_inner)
% legend('Radii ratio - r_out/r_in')

%% Functions

function [new_smallest_x, new_smallest_y, new_angle_set, new_x, new_y] = find_smallest(angle_set, old_smallest_x, old_smallest_y, x_coords, y_coords)
% While array of angles is not empty
% Find the value of smallest angle(s)
% Find the locations of smallest angle(s)
% Create current array of smallest angles with coords of point that make the smallest
% angle
% New array of smallest is previous (old) smallest plus current smallest
% Remove those elements from set of angles and coordinates
[~, min_loc] = min(angle_set);

count2 = 1;
for i = 1:length(min_loc)
    current_smallest_x(count2) = x_coords(min_loc(i));
    current_smallest_y(count2) = y_coords(min_loc(i));
    count2 = count2 + 1;
end 

new_smallest_x = [old_smallest_x current_smallest_x];
new_smallest_y = [old_smallest_y current_smallest_y];

for i = 1:length(min_loc)
    %delete the element at that location
    new_angle_set = angle_set([1:min_loc(i)-1, min_loc(i)+1:end]);
    new_x = x_coords([1:min_loc(i)-1, min_loc(i)+1:end]);
    new_y = y_coords([1:min_loc(i)-1, min_loc(i)+1:end]);
end 


end