%Muscle Parameter Default Guess
%Authors: Ashleigh L A Wiseman, James Charles, John R Hutchinson. Simple versus complex muscle modelling in extinct species: A biomechanical case study of the Australopithecus afarensis pelvis and lower extremity

%Code last checked: 03-05-2023
%This code is modified from Bishop et al. 2020: doi:10.1017/pab.2020.46
%Modified code written with assistance by Oliver Demuth.

%This script calculates the 'null hypothesis' estimate of normalized
%muscle mass, fibre length and pennation angle for a given architectural dataset. 
%It does so by fitting an alpha shape to the (scaled) data and then determining 
%the centre of the largest inscried circle for that alpha shape, using the 
%Euclidean distance transform. The mean of the centre of this circle and
%the arithmetic mean of the data is taken as the default guess.
%This more sophisticated approach is taken, as opposed to simply computing
%the centroid of the convex hull of the data, to work better with data that
%occupy a more 'boomerang' spatial distribution.


clear all;close all;clc;

%% Load data for a given muscle, stored in .csv format
[file,path]=uigetfile({'*.csv*'},'Pick .csv containing normallized muscle data');

if file==0
    return
end

Data = csvread([path,file]);
Data = Data(:,1:3);

cd(path);

%% Perform calculations

%Calculate alpha shape, using default alpha radius of 30 (seems to work 
%well on test datasets)
Shape = alphaShape(Data(:,1),Data(:,2),30);
[edges,vertices] = boundaryFacets(Shape);

%Scale results so as to remove effects of (potentially) greatly different
%dimensions in the x and y values.
Min_x = min(Data(:,1));
Min_y = min(Data(:,2));
Min_z = min(Data(:,3));
Range_x = max(Data(:,1))-min(Data(:,1));
Range_y = max(Data(:,2))-min(Data(:,2));
Range_z = max(Data(:,3))-min(Data(:,3));

vertices(:,1) = (vertices(:,1)-Min_x)/Range_x;
vertices(:,2) = (vertices(:,2)-Min_y)/Range_y;
%vertices(:,3) = (vertices(:,3)-Min_z)/Range_z;

x = vertices(:,1);
y = vertices(:,2);
%z = vertices(:,3);

vtx=zeros(size(Data));

vtx(:,1) = (Data(:,1)-Min_x)/Range_x;
vtx(:,2) = (Data(:,2)-Min_y)/Range_y;
vtx(:,3) = (Data(:,3)-Min_z)/Range_z;

AlphaX = vtx(:,1);
AlphaY = vtx(:,2);
AlphaZ = vtx(:,3);

Datashape = alphaShape(AlphaX,AlphaY,AlphaZ,30);



%Convert polygon into a 1000x1000 pixel image
scalingFactor = 1000;
x2 = x*scalingFactor + 1;
y2 = y*scalingFactor + 1;
mask = poly2mask(x2, y2, scalingFactor+1, scalingFactor+1);

%Compute the Euclidean Distance Transform of the mask image
EDT = bwdist(~mask);

%The centre of the largest inscribed circle is located where the EDT is 
%at a maximum
radius = max(EDT(:));
[yC,xC] = find(EDT == radius);
xC = xC(1);
yC = yC(1);

%Scale and shift the center back to the original dimensions of the data.
xC = (xC-1)/scalingFactor;
yC = (yC-1)/scalingFactor;
radius = radius/scalingFactor;
Alpha_xC = xC*Range_x+Min_x;
Alpha_yC = yC*Range_y+Min_y;
radius_X = radius*Range_x;
radius_Y = radius*Range_y;

%Calculate centroid of convex hull of data
[ConHullCentroid,~] = ConvHullCentroid(Data);

%Calculate arithmetic centroid of data
ArithmeticCentroid = mean(Data);

results = zeros(3,3);
resultPlot = zeros(3,3);

results(1,1) = Alpha_xC; 
results(1,2) = Alpha_yC;
results(2,:) = ConHullCentroid;
results(3,:) = ArithmeticCentroid;

resultPlot(:,1) = (results(:,1)-Min_x)/Range_x;
resultPlot(:,2) = (results(:,2)-Min_y)/Range_y;
resultPlot(:,3) = (results(:,3)-Min_z)/Range_z;

%% Report results to user
fprintf(['Alpha shape centroid = ' num2str(Alpha_xC) ', ' num2str(Alpha_yC) '\n'])
fprintf(['Convex hull centroid = ' num2str(ConHullCentroid(1)) ', ' num2str(ConHullCentroid(2)) ', ' num2str(ConHullCentroid(3)) '\n'])
fprintf(['Arithmetic centroid = ' num2str(ArithmeticCentroid(1)) ', ' num2str(ArithmeticCentroid(2)) ', ' num2str(ArithmeticCentroid(3)) '\n'])

%polyVertices = vertices;
%polyVertices(:,3) = [];

polygon = polyshape(vertices);
polygon2 = polyshape(vertices);
polygon.Vertices(:,1) = polygon.Vertices(:,1)*Range_x+Min_x;
polygon.Vertices(:,2) = polygon.Vertices(:,2)*Range_y+Min_y;

figure('Color',[1 1 1]); Shape = plot(Datashape);
set(Shape, 'facealpha', '0.3');
xlabel('Muscle mass/Body mass');
ylabel('Fibre length/Body mass^{0.33}')
zlabel('Pennation Angle (rads)');
hold on
ptc=pointCloud(resultPlot);
[labels,numClusters]=pcsegdist(ptc,0.0001);
s=scatter3(ptc.Location(1,1),ptc.Location(1,2),0,'+','MarkerEdgeColor','r');
t=scatter3(ptc.Location(2,1),ptc.Location(2,2),ptc.Location(2,3),'+','MarkerEdgeColor','k');
u=scatter3(ptc.Location(3,1),ptc.Location(3,2),ptc.Location(3,3),'+','MarkerEdgeColor','b');
s.SizeData = 100;
t.SizeData = 100;
u.SizeData = 100;
pz = patch(polygon2.Vertices(:,1),polygon2.Vertices(:,2),'g');
%h=pcshow(ptc.Location,labels,'MarkerSize', 100)
%colormap(parula(numClusters))
legend('Alpha shape','Alpha centroid','Arithmetic centroid','Convex centroid');
%h=pcshow(resultPlot, 'MarkerSize', 100);

hold off

%Plot results
figure('Color',[1 1 1])
scatter(Data(:,1),Data(:,2),30,'b','filled');
hold on
plot(polygon);
h=scatter(Alpha_xC,Alpha_yC,'+r');
scatter(ArithmeticCentroid(1),ArithmeticCentroid(2),'+k');
scatter(ConHullCentroid(1),ConHullCentroid(2),'+b');
rectangle('Position',[Alpha_xC-radius_X, Alpha_yC-radius_Y, 2*radius_X, 2*radius_Y],'Curvature',[1,1],'LineStyle','--');
hold off
xlabel('Muscle mass/Body mass');
ylabel('Fibre length/Body mass^{0.33}')
legend('Data','Alpha shape','Alpha centroid','Arithmetic centroid','Convex centroid');
ax = ancestor(h, 'axes');
ax.XAxis.Exponent = 0;
%% Helper function to calculate centroid and area of convex hull
function [Centroid,TotalArea] = ConvHullCentroid(Data)

Min_x = min(Data(:,1));
Min_y = min(Data(:,2));
Min_z = min(Data(:,3));
Range_x = max(Data(:,1))-min(Data(:,1));
Range_y = max(Data(:,2))-min(Data(:,2));
Range_z = max(Data(:,3))-min(Data(:,3));

vtx=zeros(size(Data));

vtx(:,1) = (Data(:,1)-Min_x)/Range_x;
vtx(:,2) = (Data(:,2)-Min_y)/Range_y;
vtx(:,3) = (Data(:,3)-Min_z)/Range_z;

x = vtx(:,1);
y = vtx(:,2);
z = vtx(:,3);

Datashape = alphaShape(x,y,z,30);
[tri,pts] = boundaryFacets(Datashape);
TR = triangulation(tri, pts);

numTriangles3D = size(TR,1);

Areas3D = zeros(numTriangles3D,1);
Centres3D = zeros(numTriangles3D,3);
    for i=1:numTriangles3D
        vertex1=[Data(TR.ConnectivityList(i,1),1),Data(TR.ConnectivityList(i,1),2),Data(TR.ConnectivityList(i,1),3)];
        vertex2=[Data(TR.ConnectivityList(i,2),1),Data(TR.ConnectivityList(i,2),2),Data(TR.ConnectivityList(i,2),3)];
        vertex3=[Data(TR.ConnectivityList(i,3),1),Data(TR.ConnectivityList(i,3),2),Data(TR.ConnectivityList(i,3),3)];

        %Compute area of each triangle using cross product
        Areas3D(i) = 0.5*norm(cross(vertex2-vertex1,vertex3-vertex1));

        %Compute centroid of each triangle
        Centres3D(i,:) = (vertex1(1:3)+vertex2(1:3)+vertex3(1:3))/3;
    end

TotalArea = sum(Areas3D);
Centroid = [sum((Centres3D(:,1).*Areas3D))/TotalArea, sum((Centres3D(:,2).*Areas3D))/TotalArea, sum((Centres3D(:,3).*Areas3D))/TotalArea];  
end