% Run this script to extract severity parameters after completing 
% acetabular parameterisation
% Requires the parameters to be copied into script

close all;
clear all;
clc;

% copy in  acetabular parameters here
% ensure 3 dots included before each new line
fossaloft = 0.45;								
xverts = [...	
	12.389002, 9.025515, -1.83061, -9.06687, -13.501043, ...
	12.592264, 1.154825, -7.479941, -16.7946, -21.330946, -19.217161,...
	11.851671, 0.00632, -14.750265,...
	10.568804, -2.092009, -17.859434, -23.521597, -24.376104, -25.760649,...
	4.764765, -4.891127, -16.631971,-21.305006,-23.364241...
	];	
yverts = [	...
	-2.38927, -13.743106, -25.128605, -18.310436, -14.916698, ...
	1.819653, -1.513599, -7.408069, -16.953693, -12.394756, -12.539417,...
	4.859539, 2.874539, -3.788959,...
	7.775708, 4.038052, 4.101864, 8.900265, 9.63377, 6.78706,...
	14.577237, 18.35484,18.693402,17.661862,11.518841...
	];
zverts = [	...
	22.818437, 18.850099, -4.304259, -13.192351, -13.172758, ...
	23.392809, 24.536729, 22.03325, 5.31563, -7.437163, -10.75511,...
	22.913499, 25.093034, 22.350862,...
	22.574052, 24.229936, 18.014818, 9.539568, 1.836977, -10.031706,...
	20.346742, 15.358934, 4.48929, -2.475462, -8.237913...
	];
					

% define up vector
vectU = [0 1];

% AV angle 
    % blue = anterior edge of acetabulum
    % red =  posterior edge of acetabulum
    % green = other points defining acetabulum
    
figure
hold on
for ii=1:3;
    %figure
    plot(xverts(ii), yverts(ii), 'ro',...
        xverts(5+ii), yverts(5+ii), 'go',...
        xverts(11+ii), yverts(11+ii), 'go',...
        xverts(14+ii), yverts(14+ii), 'go',...
        xverts(20+ii), yverts(20+ii), 'bo',...
        [xverts(ii) xverts(20+ii)], [yverts(ii) yverts(20+ii)],'k-',...
        [xverts(ii) xverts(ii)], [yverts(ii) 0],'k-')
    axis equal  
    title('AV angle calculated from points in transverse plane of pelvis')
    
    % define vector between points
    vectA{ii} = [xverts(20+ii)-xverts(ii) yverts(20+ii)-yverts(ii)];

    % AV angle is angle between up vector and vector between edges of acetabulum 
    costheta = dot(vectU,vectA{ii})/(norm(vectU)*norm(vectA{ii}));
    theta = acos(costheta);
    AV(ii) = radtodeg(theta);
end

for ii =4:5;
    %figure
    plot(xverts(ii), yverts(ii), 'ro',...
        xverts(5+ii), yverts(5+ii), 'go',...
        xverts(14+ii), yverts(14+ii), 'go',...
        xverts(20+ii), yverts(20+ii), 'bo',...
        [xverts(ii) xverts(20+ii)], [yverts(ii) yverts(20+ii)],'k-',...
        [xverts(ii) xverts(ii)], [yverts(ii) 0],'k-')
    axis equal  
    title('AV angle calculated from points in transverse plane of pelvis')    

    % define vector between points
    vectA{ii} = [xverts(20+ii)-xverts(ii) yverts(20+ii)-yverts(ii)];

    % AV angle is angle between up vector and vector between edges of acetabulum 
    costheta = dot(vectU,vectA{ii})/(norm(vectU)*norm(vectA{ii}));
    theta = acos(costheta);
    AV(ii) = radtodeg(theta);
   
end 

% CE angle calculated from: vertical from centre, and centre to covering edge
% of acetabulum

CE = [];

figure 
hold on
for ii=[1,5]
%     figure 
    plot(xverts((ii*5)-4),zverts((ii*5)-4), 'ro',...
        xverts((ii*5)-3),zverts((ii*5)-3), 'go',...
        xverts((ii*5)-2),zverts((ii*5)-2), 'go',...
        xverts((ii*5)-1),zverts((ii*5)-1), 'go',...
        xverts((ii*5)-0),zverts((ii*5)-0), 'go',...
        0,0, 'bo',...
        [0 xverts((ii*5)-4)], [0 zverts((ii*5)-4)],'k-',...
        [0 0], [0 30],'k-')
    axis equal  
    title('CE angle calculated from points in coronal plane of pelvis')

    
    % define vector from centre to top acetabular edge point
    vectC0{ii} = [0-xverts((ii*5)-4), 0-zverts((ii*5)-4)];

    % CE angle is angle between up vector and vector from centre to
    % acetabulum edge
    costheta = dot(vectU,vectC0{ii})/(norm(vectU)*norm(vectC0{ii}));
    theta = acos(costheta);
    CE(ii) = radtodeg(theta);
end 

for ii=2
    % figure
    plot(xverts((ii*5)-4),zverts((ii*5)-4), 'ro',...
        xverts((ii*5)-3),zverts((ii*5)-3), 'go',...
        xverts((ii*5)-2),zverts((ii*5)-2), 'go',...
        xverts((ii*5)-1),zverts((ii*5)-1), 'go',...
        xverts((ii*5)),zverts((ii*5)), 'go',...
        xverts((ii*5)+1),zverts((ii*5)+1), 'go',...
        0,0, 'bo',...
        [0 xverts((ii*5)-4)], [0 zverts((ii*5)-4)],'k-',...
        [0 0], [0 30],'k-')
    axis equal  
    title('CE angle calculated from points in coronal plane of pelvis')

    
    % define vector from centre to top acetabular edge point
    vectC0{ii} = [0-xverts((ii*5)-4), 0-zverts((ii*5)-4)];

    % CE angle is angle between up vector and vector from centre to
    % acetabulum edge
    costheta = dot(vectU,vectC0{ii})/(norm(vectU)*norm(vectC0{ii}));
    theta = acos(costheta);
    CE(ii) = radtodeg(theta);
end 

for ii=4
    % figure
    plot(xverts((ii*5)-5),zverts((ii*5)-5), 'ro',...
        xverts((ii*5)-4),zverts((ii*5)-4), 'go',...
        xverts((ii*5)-3),zverts((ii*5)-3), 'go',...
        xverts((ii*5)-2),zverts((ii*5)-2), 'go',...
        xverts((ii*5)-1),zverts((ii*5)-1), 'go',...
        xverts((ii*5)),zverts((ii*5)), 'go',...
        0,0, 'bo',...
        [0 xverts((ii*5)-5)], [0 zverts((ii*5)-5)],'k-',...
        [0 0], [0 30],'k-')
    axis equal  
    title('CE angle calculated from points in coronal plane of pelvis')

    
    % define vector from centre to top acetabular edge point
    vectC0{ii} = [0-xverts((ii*5)-5), 0-zverts((ii*5)-5)];

    % CE angle is angle between up vector and vector from centre to
    % acetabulum edge
    costheta = dot(vectU,vectC0{ii})/(norm(vectU)*norm(vectC0{ii}));
    theta = acos(costheta);
    CE(ii) = radtodeg(theta);
end 

for ii=3
    % figure
    plot(xverts((ii*5)-3),zverts((ii*5)-3), 'ro',...
        xverts((ii*5)-2),zverts((ii*5)-2), 'go',...
        xverts((ii*5)-1),zverts((ii*5)-1), 'go',...
        0,0, 'bo',...
        [0 xverts((ii*5)-3)], [0 zverts((ii*5)-3)],'k-',...
        [0 0], [0 30],'k-')
    axis equal  
    title('CE angle calculated from points in coronal plane of pelvis')

    
    % define vector from centre to top acetabular edge point
    vectC0{ii} = [0-xverts((ii*5)-3), 0-zverts((ii*5)-3)];

    % CE angle is angle between up vector and vector from centre to
    % acetabulum edge
    costheta = dot(vectU,vectC0{ii})/(norm(vectU)*norm(vectC0{ii}));
    theta = acos(costheta);
    CE(ii) = radtodeg(theta);
end 


% define AVang as the mean anteversion calculated 
% notice that level of error depends on how close z values are between slices 
meanAV = mean(AV);
minAV = min(AV);


% define CEang as the max centre edge calculated
CE = 180 - CE;
meanCE = mean(CE);
maxCE = max(CE);


% pass outputs into a Python file 
fid = fopen('opparams.py', 'w');

fprintf(fid,'CE = ['); 
for j = 1:5
    fprintf(fid,'%.12f, ',CE(j));
end
fprintf(fid,']\n'); 

fprintf(fid,'meanCE = %0.12f\n',meanCE);

fprintf(fid,'maxCE = %0.12f\n',maxCE);

fprintf(fid,'AV = ['); 
for j = 1:5
    fprintf(fid,'%.12f, ',AV(j));
end
fprintf(fid,']\n'); 

fprintf(fid,'meanAV = %0.12f\n',meanAV);

fprintf(fid,'minAV = %0.12f\n',minAV);
