diff --git a/anatomy/scapula_calculation.m b/anatomy/scapula_calculation.m
index c58d372..826febb 100644
--- a/anatomy/scapula_calculation.m
+++ b/anatomy/scapula_calculation.m
@@ -1,1119 +1,1155 @@
 %% scapula_calculation
 %
 % This script makes the anatomic calculation of the glenoid from the Amira files.
 % It is used by the scripts scapula_measure.m and scapula_addmeasure.
 % It should normally not be used directly.
 % Author: EPFL-LBO
 % Date: 2016-11-18
 %
 %%
 
 % Subject is the patient number
 % directory is the location of the datas
 % location it 'Pathologic' or 'Normal'
 % output is 'References', 'obliqueSlice', 'density' and 'display'
 
 
 
 
 
 
 function [Result] = scapula_calculation(subject, directory, location, output)
 
 segmentedScapula = 0;
 
 CTn = strtok(subject,'-');
 CTn = strtok(CTn, location);
 CTn = str2num(CTn); %#ok<ST2NM>
 CTn = sprintf('%s%03d',location,CTn);
 
 
 % Import Data
 
 
 % Import Glenoid Surface
 if exist(sprintf('%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn),'file') == 2
     [p,~,~]=importStl(sprintf('%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn),1);
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find glenoid surface file: %s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn)
 end
 
 % Import Scapula Surface
 if exist(sprintf('%s/%s/CT-%s-1/amira/scapula_%s.stl',directory,subject,subject,CTn),'file') == 2
     [pscapula,~,~]=importStl(sprintf('%s/%s/CT-%s-1/amira/scapula_%s.stl',directory,subject,subject,CTn),1);
     segmentedScapula = 1;
 end
 
 % Import HumeralHead landmarks
 if exist(sprintf('%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData1 = importdata(sprintf('%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     HH = newData1.data;
 else
     warning('Using old humeral landmarks definition')
     if exist(sprintf('%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
         newData1 = importdata(sprintf('%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
         HH = newData1.data;
     else
         error('MATLAB:rmpath:DirNotFound',...
             'Could not find humeral head landmarks file: %s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
     end
 end
 
 % Import Scapula landmarks
 % AI: Angulus Inferior
 % TS: Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine
 % PC: Processus Coracoideus (most lateral point)
-% AC: Acromio-clavicular joint (most lateral part on acromion)
+% AL: Acromion Lateral (most lateral part on acromion) // Before, it was AC: Acromio-clavicular joint (most lateral part on acromion)
 % AA: Angulus Acromialis (most laterodorsal point)
 % AG: Acromion-Glenoid notch
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData2 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     landmarks = newData2.data;
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find scapula landmarks file: %s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 
 end
 
 % Import Scapula Pillar
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData3 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     SP = newData3.data; %Scapula Pillar points
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find axillary border landmarks file: %s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 end
 
 % Import Scapula Groove
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData4 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     SG = newData4.data; %Scapula Groove points
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find supraspinatus fossa landmarks file: %s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 end
 
 % Define anatomical axes
 
 MedLat = landmarks(6,:)-mean(SG);           % Z ?
 AntPos = landmarks(5,:)-landmarks(3,:);     % X ?
 InfSup = landmarks(6,:)-mean(SP);           % Y ?
 
 
 % Find scapula side
 
 A = cross(MedLat,AntPos);
 B = vrrotvec(InfSup,A);    %Comparing Z-axes
 if rad2deg(B(4))>90
     Side = 0; %Right
     p(:,1) = -p(:,1);
     if segmentedScapula == 1
         pscapula(:,1) = -pscapula(:,1); 
     end
     HH(:,1) = -HH(:,1);
     landmarks(:,1) = -landmarks(:,1);
     SG(:,1) = -SG(:,1);
     MedLat = landmarks(6,:)-mean(SG);
     AntPos = landmarks(5,:)-landmarks(3,:);
 else
     Side = 1; %Left
 end
 
 
 %% Fit plane and axis
 
 % Find scapular plane (fit plane to AI and TS and most distal SG)
 % The 10 next line are used to determine the scapulaPlaneNormal regarding
 % the 3 most distal points of the scapulagroove.
 
 
 AI = landmarks(1,:); % Angulus inferioris point
 TS = landmarks(2,:); % Trigonum Spinae point
 EuclidDistance = zeros(5,1);
 for i=1:5 % Calculate the distance between the TS point and each of the 5 groove points
     EuclidDistance(i) = sqrt((landmarks(2,1)-SG(i,1))^2 + (landmarks(2,2)-SG(i,2))^2 + (landmarks(2,3)-SG(i,3))^2);
 end
 
 maxEuclidianDistance = max(EuclidDistance); % Find the largest distance
 pos = EuclidDistance == maxEuclidianDistance;
 GRL = SG(pos,:); % set most lateral point of the groove
 [ScapulaPlaneNormal, PlaneMean, ~, PlaneRMSE, ~] = fitPlane2([AI;TS;GRL]);
 
 
 % Find scapular plane (fit plane to SP and SG)
 %[ScapulaPlaneNormal, PlaneMean, ~, PlaneRMSE, R2Plane] = fitPlane2([SG;SP]);
 
 C = vrrotvec(ScapulaPlaneNormal,AntPos);
 if rad2deg(C(4))>90
     ScapulaPlaneNormal = -ScapulaPlaneNormal;
 end
 
 % First project, then fit
 SG_proj = project2Plane(SG,ScapulaPlaneNormal',[0 0 0],size(SG,1)); % project SG points on scapular plane
 [GrooveAxis, GrooveMean, ~, ~, ~] = fitLine2(SG_proj);
-TransverseAxis = (GrooveAxis/norm(GrooveAxis))';
+TransverseAxis = (GrooveAxis/norm(GrooveAxis))'; % TransverseAxis is the norm of the line that was created by the projection of SG points on scapular plane
 
 D = vrrotvec(TransverseAxis,MedLat);
 if rad2deg(D(4))>90
     TransverseAxis = -TransverseAxis;
 end
 
 % Project GrooveMean to scapular plane
 TransverseAxisMean = project2Plane(GrooveMean,ScapulaPlaneNormal',PlaneMean,1);
 
 % Project AG to transverseAxis
 Origin = TransverseAxisMean + dot((landmarks(6,:)-TransverseAxisMean),TransverseAxis)*TransverseAxis;
 
 % Scapula axial plane normal
 AxialPlane = cross(TransverseAxis, ScapulaPlaneNormal);
 E = vrrotvec(AxialPlane, [0 0 1]);
 CTorientation = rad2deg(E(4));
 if CTorientation > 90
     CTorientation = 180 - CTorientation;
 end
 
-% Closest point
+% Closest point. GC is the closet point between the mean of the extracted
+% surface and all the points of the extracted surface.
+% GC = Glenoid Centre
 GC = mean(p);
 vect = zeros(size(p));
 distance3D = zeros(size(p,1),1);
 for i=1:size(p,1)
     vect(i,:) = p(i,:)-GC;
     distance3D(i) = norm(vect(i,:));
 end
 ClosestNode = distance3D == min(distance3D);
 GC = p(ClosestNode,:);
 
 
 % Fit Sphere on Glenoid Surface
 [GlenoidSphere,GlenoidRadius,GlenoidResiduals, ~] = fitSphere2(p);
 GlenoidRMSE =  norm(GlenoidResiduals)/sqrt(length(GlenoidResiduals));
 
 
 % Fit sphere on Humeral Head landmarks
 [HC,HHRadius, ~, ~] = fitSphere2(HH);
 
 
 % Glenoid version 3D
 GO = GlenoidSphere'-GC;
 
 Version3DRot = vrrotvec(TransverseAxis,GO);
 Version3D = rad2deg(Version3DRot(4));
 
 VersionDirectionRotX = vrrotvec(Version3DRot(1:3),ScapulaPlaneNormal);
 VersionDirectionRotZ = vrrotvec(Version3DRot(1:3),cross(TransverseAxis,ScapulaPlaneNormal));
 
 if rad2deg(VersionDirectionRotX(4)) > 90
     VersionDirection = rad2deg(VersionDirectionRotZ(4));
 else
     VersionDirection = -rad2deg(VersionDirectionRotZ(4));
 end
 if VersionDirection < -180
     VersionDirection = VersionDirection + 360;
 end
 if VersionDirection > 180
     VersionDirection = VersionDirection - 360;
 end
 
 
 % Maximum wear plane
 MaxWearPlane = cross(GO, TransverseAxis);
 F = vrrotvec(MaxWearPlane, [0 0 1]);
 WearPlaneAngle = rad2deg(F(4));
 if WearPlaneAngle > 90
     WearPlaneAngle = 180 - WearPlaneAngle;
 end
 
 
 % Inclination : Angle between GO projected on the scapular plane and the medio-lateral axis
 proj_GO = GO' - dot(GO', ScapulaPlaneNormal) * ScapulaPlaneNormal;
 InclinationRot = vrrotvec(TransverseAxis,proj_GO);
 Inclination = rad2deg(InclinationRot(4));
 if dot(InclinationRot(1:3), ScapulaPlaneNormal) > 0
     Inclination = -Inclination;
 end
 
 
 % Version 2D : Angle between GO projected on the scapular AXIAL plane and the medio-lateral axis
 proj_GO2 = GO' - dot(GO', AxialPlane') * AxialPlane';
 Version2DRot = vrrotvec(TransverseAxis,proj_GO2);
 Version2D = rad2deg(Version2DRot(4));
 if dot(Version2DRot(1:3), AxialPlane') > 0
     Version2D = -Version2D;
 end
 
 
 % Scapulo Humeral Subluxation Index
 HO = HC'-GC;    % glenoid center to humeral center
 SHeccentricity = HO-(dot(HO,TransverseAxis)*TransverseAxis);% vector from HC to transverse axis (perpendicular)
 SHSI = norm(SHeccentricity) / HHRadius; % Mod. on 04.06.2014
 
 SHSIDirectionRotX = vrrotvec(SHeccentricity,ScapulaPlaneNormal);
 SHSIDirectionRotZ = vrrotvec(SHeccentricity,cross(TransverseAxis,ScapulaPlaneNormal));
 
 if rad2deg(SHSIDirectionRotZ(4)) < 90
     SHSIDirection = rad2deg(SHSIDirectionRotX(4));
 else
     SHSIDirection = -rad2deg(SHSIDirectionRotX(4));
 end
 if SHSIDirection < -180
     SHSIDirection = SHSIDirection+360;
 end
 
 SHSPlane = cross(TransverseAxis, HO);
 G = vrrotvec(SHSPlane, [0 0 1]);
 SHSPlaneAngle = rad2deg(G(4));
 if SHSPlaneAngle > 90
     SHSPlaneAngle = 180 - SHSPlaneAngle;
 end
 
 
 % Gleno Humeral Subluxation Index
 GHeccentricity = HO-dot(HO,GO/norm(GO))*(GO/norm(GO));%vector from glenoid sphere centre to glenoid centerline (perpendicular)
 GHSI = norm(GHeccentricity) / HHRadius; % Mod. on 04.06.2014
 
 ScapulaPlaneNormal2Glenoid = project2Plane(ScapulaPlaneNormal',GO,[0 0 0],1);% project y axis on glenoid plane
 
 GHSIDirectionRotX = vrrotvec(GHeccentricity,ScapulaPlaneNormal2Glenoid);
 GHSIDirectionRotZ = vrrotvec(GHeccentricity,cross(GO,ScapulaPlaneNormal2Glenoid));
 
 if rad2deg(GHSIDirectionRotZ(4)) < 90
     GHSIDirection = rad2deg(GHSIDirectionRotX(4));
 else
     GHSIDirection = -rad2deg(GHSIDirectionRotX(4));
 end
 
 if GHSIDirection < -180
     GHSIDirection = GHSIDirection+360;
 end
 
 GHSPlane = cross(GO, HO);
 H = vrrotvec(GHSPlane, [0 0 1]);
 GHSPlaneAngle = rad2deg(H(4));
 if GHSPlaneAngle > 90
     GHSPlaneAngle = 180 - GHSPlaneAngle;
 end
 
 
 % Height of the glenoid cap
 % Project the points of the surface on the centerline
 p_proj = zeros(size(p));
 for i = 1:size(p,1)
     p_proj(i,:) = GC + (dot(GO,(p(i,:)-GC)) / dot(GO,GO)) * GO;
     distance3D(i) = norm(GC-p_proj(i,:));
 end
 
 CapHeight = max(distance3D);
 
 
 % Glenoid height and width
 
 % Scapula coordinate system in CT frame
 ScapulaCSYS(1,:) = TransverseAxis; % X
 ScapulaCSYS(2,:) = ScapulaPlaneNormal; % Y
 ScapulaCSYS(3,:) = AxialPlane; % Z
 ScapulaCSYS(4,:) = Origin;
 ScapulaCSYS(5,:) = Origin + 10*ScapulaCSYS(1,:);
 ScapulaCSYS(6,:) = Origin + 10*ScapulaCSYS(2,:);
 ScapulaCSYS(7,:) = Origin + 10*ScapulaCSYS(3,:);
 
 % Contruction of transformation matrix from CT frame to Scapula
 % frame
 Osc = ScapulaCSYS(4,:);
 Oscx = Origin + ScapulaCSYS(1,:);
 Oscy = Origin + ScapulaCSYS(2,:);
 Oscz = Origin + ScapulaCSYS(3,:);
 
 O = [0,0,0];
 Ox = [1,0,0];
 Oy = [0,1,0];
 Oz = [0,0,1];
 
 A = [Osc(1),Osc(2),Osc(3),1;
     Oscx(1),Oscx(2),Oscx(3),1;
     Oscy(1),Oscy(2),Oscy(3),1;
     Oscz(1),Oscz(2),Oscz(3),1];
 
 bx = [O(1);
     Ox(1);
     Oy(1);
     Oz(1)];
 
 by = [O(2);
     Ox(2);
     Oy(2);
     Oz(2)];
 
 
 bz = [O(3);
     Ox(3);
     Oy(3);
     Oz(3)];
 
 
 B = [bx,by,bz];
 X = A\B;
 T = [X';0,0,0,1]; % Transformation matrix from CT frame to scapula frame
 
 
 psc = zeros(size(p,1),4);
 for i = 1:size(p,1)
     psc(i,:) = T*[p(i,:),1]';
 end
 
 psc_y = psc(:,2);
 psc_z = psc(:,3);
 
 % Difference of extremal values
 i_y_max=1;
 for i = 2:1:size(psc_y)
 
     if psc_y(i) > psc_y(i_y_max)
         i_y_max = i;
     end
 end
 
 i_y_min=1;
 for i = 2:1:size(psc_y)
 
     if psc_y(i) < psc_y(i_y_min)
         i_y_min = i;
     end
 end
 
 i_z_max=1;
 for i = 2:1:size(psc_z)
 
     if psc_z(i) > psc_z(i_z_max)
         i_z_max = i;
     end
 end
 
 i_z_min=1;
 for i = 2:1:size(psc_z)
 
     if psc_z(i) < psc_z(i_z_min)
         i_z_min = i;
     end
 end
 
 scapulaWidth = norm(p(i_y_max,:)-p(i_y_min,:));
 scapulaHeight = norm(p(i_z_max,:)-p(i_z_min,:));
 
 % Acromion Index (AI)
 % AI = GA/GH, where:
-% GA: GlenoAcromial distance in scapula plane = distance from AC to glenoid principal axis
+% GA: GlenoAcromial distance in scapula plane = distance from AL to glenoid principal axis
 % GH: GlenoHumeral distance = 2*HHRadius, humeral head diameter (radius*2) 
 
-AC = landmarks(4,:); % AC: Acromio-clavicular joint (most lateral part on acromion)
+%Calculate the correct AL landmarks if the scapula is segmented
+if segmentedScapula == 1
+    % Transform scapula, glenoid and acromion points to Scapula Coordinate System
+    ex = [1;0;0];
+    ey = [0;1;0];
+    ez = [0;0;1];  
+
+    ev = AxialPlane'; % Y
+    ew = TransverseAxis'; %Z
+    eu = cross(ev, ew);  % X
+
+    R = eu * ex' + ev * ey' + ew * ez';
+    
+    % Transform scapula point to Scapula Coordinate System
+    ScapulaPtsinScapulaCSYS = pscapula*R;
+    % Take the most lateral point (assumption: The most lateral point of the scaplua is always on the acromion
+    [~,ALpositionInScapula] = max(ScapulaPtsinScapulaCSYS(:,3)); 
+    AL = pscapula(ALpositionInScapula,:);% AL: Acromion lateral (most lateral part on acromion) 
+else
+   AL = landmarks(4,:); % AL: Acromion lateral (most lateral part on acromion) 
+end
 
-ACinScapulaPlane = project2Plane(AC,ScapulaPlaneNormal',PlaneMean,size(AC,1));
+% Project AL, glenoid points and GC in scapula plane
+ALinScapulaPlane = project2Plane(AL,ScapulaPlaneNormal',PlaneMean,size(AL,1));
 GlenoidPtsinScapulaPlane = project2Plane(p,ScapulaPlaneNormal',PlaneMean,size(p,1));
 GCinScapulaPlane = project2Plane(GC,ScapulaPlaneNormal',PlaneMean,size(GC,1));
 if segmentedScapula == 1
     ScapulaPtsinScapulaPlane = project2Plane(pscapula,ScapulaPlaneNormal',PlaneMean,size(pscapula,1));
 end
 
 % Figure in 3D to show position of scapula plane, glenoid points and AC
 
 figure;
 d = PlaneMean*ScapulaPlaneNormal;
 if segmentedScapula == 1
     [xx,yy]=ndgrid(min(pscapula(:,1)):max(pscapula(:,1)),min(pscapula(:,2)):max(pscapula(:,2)));
 else
     [xx,yy]=ndgrid((min(p(:,1))-50):(max(p(:,1))+50),(min(p(:,2))-50):(max(p(:,2))+50));
 end    
 zz=(d-ScapulaPlaneNormal(1)*xx-ScapulaPlaneNormal(2)*yy)/ScapulaPlaneNormal(3);
 hold on
 surf(xx,yy,zz,'FaceColor','r','FaceAlpha',0.5, 'EdgeColor', 'none')
 if segmentedScapula == 1
     scatter3(ScapulaPtsinScapulaPlane(:,1),ScapulaPtsinScapulaPlane(:,2),ScapulaPtsinScapulaPlane(:,3),200,[0.5 0.5 0.5],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)    
     hold on;
 end
-scatter3(ACinScapulaPlane(:,1),ACinScapulaPlane(:,2),ACinScapulaPlane(:,3),[],[0 0 1],'filled')
+scatter3(ALinScapulaPlane(:,1),ALinScapulaPlane(:,2),ALinScapulaPlane(:,3),[],[0 0 1],'filled')
 hold on
 scatter3(GlenoidPtsinScapulaPlane(:,1),GlenoidPtsinScapulaPlane(:,2),GlenoidPtsinScapulaPlane(:,3),200,[0 0 1],'Marker','.')
 hold on
 scatter3(GCinScapulaPlane(:,1),GCinScapulaPlane(:,2),GCinScapulaPlane(:,3),[],'r','filled')
 hold on
 daspect([1 1 1])
 xlabel('X (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 zlabel('Z (Scapula CSYS)')
 daspect([1 1 1])
-savefig(sprintf('%s/%s/CT-%s-1/amira/Projection2ScapulaPlane_%s',directory,subject,subject,CTn))
+if segmentedScapula == 1
+    savefig(sprintf('%s/%s/CT-%s-1/amira/Projection2ScapulaPlaneWithSegmentedScapula_%s',directory,subject,subject,CTn))
+else
+    savefig(sprintf('%s/%s/CT-%s-1/amira/Projection2ScapulaPlane_%s',directory,subject,subject,CTn))
+end
 
 % Transform scapula, glenoid and acromion points to Scapula Coordinate System
 
 ex = [1;0;0];
 ey = [0;1;0];
 ez = [0;0;1];  
 
 ev = AxialPlane'; % Y
 ew = TransverseAxis'; %Z
 eu = cross(ev, ew);  % X
 
-R = eu * ex' + ev * ey' + ew * ez';
+R = eu * ex' + ev * ey' + ew * ez'; % Create rotation matrix to transform points to Scapula Coordinate System
 
-ACinScapulaCSYS = AC*R;
+ALinScapulaCSYS = AL*R;
 GlenoidPtsinScapulaCSYS = p*R;
 GCinScapulaCSYS = GC*R;
 TrinScapulaCSYS = TransverseAxis*R;
 ScapulaPlaneNormalinScapulaCSYS = ScapulaPlaneNormal'*R;
 PlaneMeaninScapulaCSYS = PlaneMean*R;
 if segmentedScapula == 1
     ScapulaPtsinScapulaCSYS = pscapula*R;
 end
 
 % Figure in 3D to show position of scapula plane, glenoid points and AC
 
 figure;
 if segmentedScapula == 1
     scatter3(ScapulaPtsinScapulaCSYS(:,1),ScapulaPtsinScapulaCSYS(:,2),ScapulaPtsinScapulaCSYS(:,3),200,[0.5 0.5 0.5],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)    
     hold on;
 end
-scatter3(ACinScapulaCSYS(:,1),ACinScapulaCSYS(:,2),ACinScapulaCSYS(:,3),[],[0 0 1],'filled')
+scatter3(ALinScapulaCSYS(:,1),ALinScapulaCSYS(:,2),ALinScapulaCSYS(:,3),[],[0 0 1],'filled')
 hold on
 scatter3(GlenoidPtsinScapulaCSYS(:,1),GlenoidPtsinScapulaCSYS(:,2),GlenoidPtsinScapulaCSYS(:,3),200,[0 0 1],'Marker','.')
 hold on
 scatter3(GCinScapulaCSYS(:,1),GCinScapulaCSYS(:,2),GCinScapulaCSYS(:,3),[],'r','filled')
 hold on
 x=eu*linspace(0,50);
 plot3(x(:,1)',x(:,2)',x(:,3)','r')
 hold on
 y=ev*linspace(0,50);
 plot3(y(:,1)',y(:,2)',y(:,3)','g')
 hold on
 z=ew*linspace(0,50);
 plot3(z(:,1)',z(:,2)',z(:,3)','b')
 daspect([1 1 1])
 d= ScapulaPlaneNormalinScapulaCSYS*PlaneMeaninScapulaCSYS';
 if segmentedScapula == 1
     [yy,zz]=ndgrid(min(ScapulaPtsinScapulaCSYS(:,2)):max(ScapulaPtsinScapulaCSYS(:,2)),min(ScapulaPtsinScapulaCSYS(:,3)):max(ScapulaPtsinScapulaCSYS(:,3)));
 else
     [yy,zz]=ndgrid((min(GlenoidPtsinScapulaCSYS(:,2))-50):(max(GlenoidPtsinScapulaCSYS(:,2))+50),(min(GlenoidPtsinScapulaCSYS(:,3))-50):(max(GlenoidPtsinScapulaCSYS(:,3))+50));
 end    
 xx = zeros(size(yy));
 xx(:,:) = d/ScapulaPlaneNormalinScapulaCSYS(1); %assuming normal(3)==0 and normal(2)==0
 hold on
 surf(xx,yy,zz,'FaceColor','r','FaceAlpha',0.5, 'EdgeColor', 'none');
 xlabel('X (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 zlabel('Z (Scapula CSYS)')
 daspect([1 1 1])
-savefig(sprintf('%s/%s/CT-%s-1/amira/Transform2ScapulaCSYS_%s',directory,subject,subject,CTn))
+if segmentedScapula == 1
+    savefig(sprintf('%s/%s/CT-%s-1/amira/Transform2ScapulaCSYSwithScapulaSegmented_%s',directory,subject,subject,CTn))
+else
+    savefig(sprintf('%s/%s/CT-%s-1/amira/Transform2ScapulaCSYS_%s',directory,subject,subject,CTn))
+end
 
-% Project scapula, glenoid, AC, GC and scapula axis to scapula plane
+% Project scapula, glenoid, AC, GC and scapula axis to scapula plane in the
+% Scapula Coordinate System
 
-ACinScapulaPlane = project2Plane(ACinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ACinScapulaCSYS,1));
+ALinScapulaPlane = project2Plane(ALinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ALinScapulaCSYS,1));
 GlenoidPtsinScapulaPlane = project2Plane(GlenoidPtsinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(GlenoidPtsinScapulaCSYS,1));
 GCinScapulaPlane = project2Plane(GCinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(GCinScapulaCSYS,1));
 TrinScapulaPlane = project2Plane(TrinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(TrinScapulaCSYS,1));
 if segmentedScapula == 1
     ScapulaPtsinScapulaPlane = project2Plane(ScapulaPtsinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ScapulaPtsinScapulaCSYS,1));
 end
 
 % Find glenoid principal axis (most superior and most inferior projected
 % glenoid points)
 
 GlenoidPrincipalAxis = [GlenoidPtsinScapulaPlane(GlenoidPtsinScapulaPlane(:,2) == min(GlenoidPtsinScapulaPlane(:,2)),:) ;GlenoidPtsinScapulaPlane(GlenoidPtsinScapulaPlane(:,2) == max(GlenoidPtsinScapulaPlane(:,2)),:)];
 
-% Compute GA (distance from AC to glenoid principal axis)
+% Compute GA (distance from AL to glenoid principal axis)
 
-GA = norm(cross(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:),ACinScapulaPlane-GlenoidPrincipalAxis(1,:)))/norm(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:));
+GA = norm(cross(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:),ALinScapulaPlane-GlenoidPrincipalAxis(1,:)))/norm(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:));
 % Compute GH (Humeral head diameter)
 
 GH = 2*HHRadius;
 
 % Compute AI
 
 AI = GA/GH;
 
 % Figure in 2D showing projected glenoid points, glenoid principal axis and
 % AC point
 
 figure;
 daspect([1 1 1])
 if segmentedScapula == 1
     scatter(ScapulaPtsinScapulaPlane(:,3),ScapulaPtsinScapulaPlane(:,2),200,[0.8 0.8 0.8],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)
     hold on
 end
-scatter(ACinScapulaPlane(:,3),ACinScapulaPlane(:,2),[],[0 0 1],'filled')
+scatter(ALinScapulaPlane(:,3),ALinScapulaPlane(:,2),[],[0 0 1],'filled')
 hold on
 scatter(GlenoidPtsinScapulaPlane(:,3),GlenoidPtsinScapulaPlane(:,2),200,[0 0 1],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)
 hold on
 scatter(GlenoidPrincipalAxis(:,3),GlenoidPrincipalAxis(:,2),[],'w','filled','MarkerEdgeColor','k')
 hold on
 scatter(GCinScapulaPlane(:,3),GCinScapulaPlane(:,2),[],'r','filled')
 hold on
 a=TrinScapulaPlane(:,2)/TrinScapulaPlane(:,3);
 b=GCinScapulaPlane(:,2)-a*GCinScapulaPlane(:,3);
 if segmentedScapula == 1
     x=linspace(min(ScapulaPtsinScapulaPlane(:,3)),max(ScapulaPtsinScapulaPlane(:,3)));
 else
     x=linspace((min(GlenoidPtsinScapulaPlane(:,3))-50),(max(GlenoidPtsinScapulaPlane(:,3)))+50);    
 end
 y=a*x+b;
 hold on
 plot(x,y,'--r')
 hold on
 a1=(GlenoidPrincipalAxis(2,2)-GlenoidPrincipalAxis(1,2))/(GlenoidPrincipalAxis(2,3)-GlenoidPrincipalAxis(1,3));
 b1=GlenoidPrincipalAxis(1,2)-a1*GlenoidPrincipalAxis(1,3);
 x=linspace(min(GlenoidPrincipalAxis(:,3))-5,max(GlenoidPrincipalAxis(:,3))+5);
 y1=a1*x+b1;
 plot(x,y1,'--b')
 daspect([1 1 1])
 xlabel('Z (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 title(['Acromion Index = ' num2str(AI) ' (GA = ' num2str(GA) ' GH = ' num2str(GH) ')'])
-saveas(gcf,sprintf('%s/%s/CT-%s-1/amira/AI_%s.png',directory,subject,subject,CTn))
+if segmentedScapula == 1
+    saveas(gcf,sprintf('%s/%s/CT-%s-1/amira/AIwithSegmentedScapula_%s.png',directory,subject,subject,CTn))
+else
+    saveas(gcf,sprintf('%s/%s/CT-%s-1/amira/AI_%s.png',directory,subject,subject,CTn))
+end
 
 close all
 
 % Output
 
 CTtoXYZ = [TransverseAxis' ScapulaPlaneNormal cross(TransverseAxis, ScapulaPlaneNormal)'];
 Cxyz = (GC - Origin) * CTtoXYZ;
 
 % Scapula coordinate system in CT frame
 ScapulaCSYS(3,:) = TransverseAxis; % Z
 ScapulaCSYS(1,:) = ScapulaPlaneNormal; % X
 ScapulaCSYS(2,:) = cross(TransverseAxis,ScapulaPlaneNormal); % Y
 ScapulaCSYS(4,:) = Origin;
 
 Result.sCase_id = strtok(subject,'-');%1
 Result.glenoid_Radius = GlenoidRadius;%2
 Result.glenoid_radiusRMSE = GlenoidRMSE;%3
 %Result.glenoidSpehricity = ' ';
 %Result.glenoid_biconcave = ' ';
 Result.glenoid_depth = CapHeight; %4
 Result.glenoid_width = scapulaWidth;%5
 Result.glenoid_height = scapulaHeight;%6
 Result.glenoid_center_PA = -Cxyz(2); %7
 Result.glenoid_center_IS = Cxyz(3); %8
 Result.glenoid_center_ML = Cxyz(1); %9
 Result.glenoid_version_ampl = Version3D;%10
 Result.glenoid_version_orient = VersionDirection;%11
 Result.glenoid_Version = Version2D; % 12
 Result.glenoid_Inclination = Inclination; %13
 % Result.glenoid_version2D = ' ';
 % Result.glenoid_walch_class = ' ';
 Result.humerus_joint_radius = ' '; %14
 Result.humeral_head_radius = HHRadius;%15
 % Result.humerus_GHsublux_2D = ' ';
 % Result.humerus_SHsublux_2D = ' ';
 Result.humerus_GHsubluxation_ampl = GHSI;%16
 Result.humerus_GHsubluxation_orient = GHSIDirection;%17
 Result.humerus_SHsubluxation_ampl = SHSI;%18
 Result.humerus_SHsubluxation_orient = SHSIDirection;%19
 % Result.scapula_CSA = ' ';
 Result.scapula_CTangle = CTorientation; %20
 Result.scapula_CTangleVersion = WearPlaneAngle; %21
 Result.scapula_CTangleSHS = SHSPlaneAngle; % 22
 Result.scapula_CTangleGHS = GHSPlaneAngle; % 23
 Result.scapula_PlaneRMSE = PlaneRMSE;%24
 Result.scapula_AI = AI;
 
 
 
 %% Additional output
 
 if exist('output','var') == 1
 
     % Scapula coordinate system in CT frame
     ScapulaCSYS(1,:) = TransverseAxis; % X
     ScapulaCSYS(2,:) = ScapulaPlaneNormal; % Y
     ScapulaCSYS(3,:) = AxialPlane; % Z
     ScapulaCSYS(4,:) = Origin;
     ScapulaCSYS(5,:) = Origin + 10*ScapulaCSYS(1,:);
     ScapulaCSYS(6,:) = Origin + 10*ScapulaCSYS(2,:);
     ScapulaCSYS(7,:) = Origin + 10*ScapulaCSYS(3,:);
 
     %% ISB coordinate system in CT frame
 
     Xisb = (landmarks(4,:)-landmarks(2,:))/norm(landmarks(4,:)-landmarks(2,:)); %Med-Lat
     Yisb = cross(Xisb,landmarks(1,:)-landmarks(2,:))/norm(cross(Xisb,landmarks(1,:)-landmarks(2,:)));%Ant-Pos
     Zisb = cross(Xisb,Yisb);%Inf-Sup
 
     %% Locate abq coordinate system in CT frame
 
     Act2isb = [Xisb' Yisb' Zisb'];
     Act2lbo = ScapulaCSYS(1:3,:)';
     Albo2isb = Act2isb/Act2lbo;
     Athx2abq = SpinCalc('EA321toDCM',[40 0 0]);% scapular plane is 40� anterior to frontal plane
     Aisb2thx = SpinCalc('EA321toDCM',[-35.6 -17.8 -2.5]);% initial position McClure et al. 2001
 
     Act2abq = Athx2abq*Aisb2thx*Albo2isb*Act2lbo;
 
     %% Output initial position of scapula for abaqus simulation in CT coordinates
 
     %     scaleFactor = 24/HHRadius;  % =24mm / Humerus radius
     ScapulaABQ(1,:) = HC;   % Humerus center
     ScapulaABQ(2,:) = 10*Act2abq(:,1)'+ScapulaABQ(1,:);
     ScapulaABQ(3,:) = 10*Act2abq(:,2)'+ScapulaABQ(1,:);
     ScapulaABQ(4,:) = 10*Act2abq(:,3)'+ScapulaABQ(1,:);
     %     ScapulaABQ(5,1) = scaleFactor;
 
     %% Output coordinates for obliqueSlice
     % The ObliqueSlice can be defined as point + (normal OR u-vector and
     % v-vector)
     FriedmanPlane = [GC 1 0 0 0 1 0];
     obliqueSlice = [landmarks(6,:) ScapulaCSYS(2,:) -ScapulaCSYS(3,:)];      % Plane for muscle measurement
     obliqueSlice2 = [GC ScapulaCSYS(3,:)];                                  % Scapula axial plane
     obliqueSlice3 = [GC ScapulaCSYS(1,:) GO];                               % Max wear plane
     obliqueSlice4 = [GC ScapulaCSYS(2,:)];                                   % Scapular plane
     obliqueSlice5 = [GC HO/norm(HO) GHeccentricity/norm(GHeccentricity)];   % Plane for GHS
     obliqueSlice6 = [GC HO/norm(HO) SHeccentricity/norm(SHeccentricity)];   % Plane for SHS
 
     if Side == 0
         FriedmanPlane(1) = -FriedmanPlane(1);
         obliqueSlice(1:3:end) = - obliqueSlice(1:3:end);
         obliqueSlice2(1:3:end) = -obliqueSlice2(1:3:end);
         obliqueSlice3(1:3:end) = -obliqueSlice3(1:3:end);
         obliqueSlice4(1:3:end) = -obliqueSlice4(1:3:end);
         obliqueSlice5(1:3:end) = -obliqueSlice5(1:3:end);
         obliqueSlice6(1:3:end) = -obliqueSlice6(1:3:end);
     end
 
 
 
 
     switch output
         case 'density'
             %% Data for the calculations of the density in amira
             % Needed variables are :
             %  - CylO: cylinder origin
             %  - CylA: cylinder alignment point
             %  - CylR_factor: scale factor for cylinder radius
             %  - Glenoid_sphere_center
             %  - Glenoid radius
             %  - CylO_shift5mm : cylinder origin shifted 5 mm behind
             %  glenoid surface center
             %  - CylA2 : cylinder alignment point n�2
             %  - CylO_shift15mmfront: cylinder origin shifted 15 mm
             %  frontwards
             %  - CylA3 : cylinder alignemnt point n�2
             %  - startcylpoint : origin of the new cylinder starting from
             %  the AG point
             %  - endcylpoint: end of cylinder, 40 mm from the startcylpoint
             %  in the scapula axis
 
 
             CylO = GC;
             Glenoid_sphere_center = GlenoidSphere';
             Glenoid_radius = GlenoidRadius;
 
             %% CylA - second point to form a line with glenoid surface
             % center point which is parallel to scapula axis
             CylA = CylO - 40 * TransverseAxis/norm(TransverseAxis);
 
             %% CylR - distance from glenoid surface center to most
             % eccentric point of glenoid surface
 
             %Projection of surface points onto the same plane
             PlanePoint = CylO;
             PlaneNormal = (CylO-CylA)/norm(CylO-CylA);
             
             ProjectedPoints = zeros(size(p));
             DistanceToCylO = zeros(size(p,1),1);
             for i = 1:size(p,1)
                 ProjectedPoints(i,:) = p(i,:) - dot(p(i,:) - PlanePoint, PlaneNormal) * PlaneNormal;
                 DistanceToCylO(i) = norm(ProjectedPoints(i,:) - CylO);
             end
 
             %Selection of the most eccentric point -> radius of cylinder
             CylR = max(DistanceToCylO);
 
             %Scale factor of the cylinder
             OriginalCylR = 5;
             CylR_factor = CylR/OriginalCylR;
 
             % Coordinates of start and end point of the D10 mm cylinder
             Deepness = 5; %behind glenoid surface
             smallCylR = 5;
             CylO_shift5mm = GC - Deepness * TransverseAxis/norm(TransverseAxis);
             CylA2 = CylO_shift5mm - smallCylR * TransverseAxis;
 
             % To place big cylinder (finer seleciton, exclusion of distant
             % bones)
             CylO_shift15mmfront = GC + 15 * TransverseAxis/norm(TransverseAxis);
             CylA3 = CylO_shift15mmfront - 40 * TransverseAxis/norm(TransverseAxis);
 
             %%Projection of AG on scapula axis
             AG = landmarks(6,:);
             p1 = CylA3;
             p2 = CylO_shift15mmfront;
 
             v1 = p2-p1;
             v2 = AG-p1;
 
             E1 = v1;
             E2 = cross(v1,v2);
             E3 = cross(E1,E2);
 
             e1 = E1/norm(E1);
             e2 = E2/norm(E2);
             e3 = E3/norm(E3);
 
             Oc = p1;
             Ocx = p1 + e1;
             Ocy = p1 + e2;
             Ocz = p1 + e3;
 
             O = [0,0,0];
             Ox = [1,0,0];
             Oy = [0,1,0];
             Oz = [0,0,1];
 
             A = [Oc(1),Oc(2),Oc(3),1;
                 Ocx(1),Ocx(2),Ocx(3),1;
                 Ocy(1),Ocy(2),Ocy(3),1;
                 Ocz(1),Ocz(2),Ocz(3),1];
 
             bx = [O(1);
                 Ox(1);
                 Oy(1);
                 Oz(1)];
 
 
             by = [O(2);
                 Ox(2);
                 Oy(2);
                 Oz(2)];
 
 
             bz = [O(3);
                 Ox(3);
                 Oy(3);
                 Oz(3)];
 
 
             B = [bx,by,bz];
             X = A\B;
             T = [X';
                 0,0,0,1];
 
             det(T);
 
             AG2 = [AG,1];
             AG3=AG2';
             diff = T*AG3;
             projAGcyl=[diff(1),0,0,1];
             startcylpoint=T\projAGcyl';
             startcylpoint = [startcylpoint(1),startcylpoint(2),startcylpoint(3)];
             endcylpoint = startcylpoint + 40 * TransverseAxis/norm(TransverseAxis);
 
             %% Write data
             if Side == 0
                 CylO(1) = -CylO(1);
                 CylA(1) = -CylA(1);
                 Glenoid_sphere_center(1) = -Glenoid_sphere_center(1);
                 CylO_shift5mm(1) = -CylO_shift5mm(1);
                 CylA2(1) = -CylA2(1);
                 CylO_shift15mmfront(1) = -CylO_shift15mmfront(1);
                 CylA3(1) = -CylA3(1);
                 startcylpoint(1) = -startcylpoint(1);
                 endcylpoint(1) = -endcylpoint(1);
 
 
 
             end
             mkdir('Generated_Amira_TCL_Scripts/CylinderReferences')
             
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),'w');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylR_factor, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_sphere_center, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_radius, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift5mm, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA2, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift15mmfront, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA3, '-append','precision',10);
             fclose(fid);
 
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s_AmiraCode.dat',CTn),'w');
 
             fprintf(fid,'set lineset1 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset1 addPoint %f %f %f]\n',CylO_shift5mm);
             fprintf(fid,'set p2 [$lineset1 addPoint %f %f %f]\n',CylA2);
             fprintf(fid,'$lineset1 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set} setIconPosition 20 10\n');
             fprintf(fid,'Line-Set fire\n\n');
 
             fprintf(fid,'set lineset2 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset2 addPoint 0 0 0]\n');
             fprintf(fid,'set p2 [$lineset2 addPoint 5 0 0]\n');
             fprintf(fid,'$lineset2 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set2} setIconPosition 20 30\n');
             fprintf(fid,'Line-Set2 fire\n\n');
 
             fprintf(fid,'set hideNewModules 0\n');
             fprintf(fid,'create {HxCreateSphere} {GlenoidSphereFit}\n');
             fprintf(fid,'{GlenoidSphereFit} setIconPosition 20 330\n');
             fprintf(fid,'{GlenoidSphereFit} fire\n');
             fprintf(fid,'{GlenoidSphereFit} {type} setValue 0\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setMinMax 0 100\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setIncrement 0.666667\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setValue %f\n',Glenoid_radius);
             fprintf(fid,'{GlenoidSphereFit} {radius} setSubMinMax 0 100\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setIncrement 0.332667\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setValue 1.07457\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setValue 1\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
             fprintf(fid,'GlenoidSphereFit fire\n');
             fprintf(fid,'GlenoidSphereFit setViewerMask 16383\n');
             fprintf(fid,'GlenoidSphereFit setPickable 1\n\n');
 
             fprintf(fid,'set hideNewModules 0\n');
             fprintf(fid,'create {HxCreateSphere} {CorticalSphereFit}\n');
             fprintf(fid,'{CorticalSphereFit} setIconPosition 20 350\n');
             fprintf(fid,'{CorticalSphereFit} fire\n');
             fprintf(fid,'{CorticalSphereFit} {type} setValue 0\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setMinMax 0 100\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setIncrement 0.666667\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setValue %f\n',Glenoid_radius+3);
             fprintf(fid,'{CorticalSphereFit} {radius} setSubMinMax 0 100\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setIncrement 0.332667\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setValue 1.07457\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setValue 1\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
             fprintf(fid,'CorticalSphereFit fire\n');
             fprintf(fid,'CorticalSphereFit setViewerMask 16383\n');
             fprintf(fid,'CorticalSphereFit setPickable 1\n\n');
 
             fprintf(fid,'set lineset3 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset3 addPoint %f %f %f]\n',startcylpoint);
             fprintf(fid,'set p2 [$lineset3 addPoint %f %f %f]\n',endcylpoint);
             fprintf(fid,'$lineset3 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set3} setIconPosition 20 520\n');
             fprintf(fid,'Line-Set3 fire\n\n');
 
             fprintf(fid,'set lineset4 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset4 addPoint 0 0 0]\n');
             fprintf(fid,'set p2 [$lineset4 addPoint 40 0 0]\n');
             fprintf(fid,'$lineset4 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set4} setIconPosition 20 540\n');
             fprintf(fid,'Line-Set4 fire\n\n');
 
             fprintf(fid,'{Cylinder2.STL} setScaleFactor 1 %f %f\n',CylR_factor,CylR_factor);
             fprintf(fid,'{Cylinder2.STL} applyTransform\n\n');
 
             fprintf(fid,'saveProject\n');
             fprintf(fid,'clear\n');
             fprintf(fid,'echo "Complete center coordinates of ReamingSphere"\n\n');
 
             %fprintf(fid,'
 
             fclose(fid);
 
         case 'References'
             
             mkdir('Generated_Amira_TCL_Scripts/References')
             
             % References data for a study case
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'w');
             fprintf(fid,'# Scapula coordinate system\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaCSYS(4:7,:), '-append','precision',10);
             fclose(fid);
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
             fprintf(fid,'\n\n# Abaqus reference\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaABQ, '-append','precision',10);
             fclose(fid);
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
             fprintf(fid,'\n\n# Glenoid centerline\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GC, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GlenoidSphere', '-append','precision',10);
             fclose(fid);
 
         case 'obliqueSlice'
             
             [token, remain] = strtok(directory,'../../..');
             directory = ['Z:/' token remain];
             pictureName = [CTn '_new'];
 
             mkdir('Generated_Amira_TCL_Scripts/obliqueSlice')
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\obliqueslice\\obliqueSlice_%s.tcl',CTn),'w');
 
             fprintf(fid,'# Create an ObliqueSlice with Friedman plane #\n');
             fprintf(fid,'ObliqueSlice orientation hit 0\n');
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\n',FriedmanPlane);
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 0 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice fire\n');
             fprintf(fid,'set planepos [ObliqueSlice translate getValue]\n');
             fprintf(fid,'set maxvalue [ObliqueSlice translate getMaxValue]\n');
             fprintf(fid,'set planepos [expr {$maxvalue - int($planepos)}]\n');
             fprintf(fid,'ObliqueSlice createImage %s_$planepos\n', CTn);
             fprintf(fid,'create HxCastField {CastField}\n');
             fprintf(fid,'CastField data connect %s_$planepos\n', CTn);
             fprintf(fid,'CastField fire\n');
             fprintf(fid,'CastField scaling setValue 0 0.21\n');
             fprintf(fid,'CastField scaling setValue 1 200\n');
             fprintf(fid,'CastField fire\n');
             fprintf(fid,'[ {CastField} create ] setLabel %s_$planepos.jpeg\n', CTn);
             if strcmp(location, 'P')
                 fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" %s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
             else
                 fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" Z:/data/%s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
             end
 
             fprintf(fid,'# Get image for muscle measurement #\n');
             fprintf(fid,'ObliqueSlice orientation hit 0\n');
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice fire\n', obliqueSlice);    % Muscle evaluation, through notch
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice createImage {%s}\n', pictureName);
             if strcmp(location, 'pathologic')
                 fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
             else
                 fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
             end
             fprintf(fid,'saveProject\n\n');
 
             fprintf(fid,'ObliqueSlice2 setPlane %f %f %f %f %f %f\nObliqueSlice2 fire\n', obliqueSlice2);                          % Scapular axial plane
             fprintf(fid,'ObliqueSlice2 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice2 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice3 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice3 fire\n', obliqueSlice3);                            % Max wear plane
             fprintf(fid,'ObliqueSlice3 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice3 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice4 setPlane %f %f %f %f %f %f\nObliqueSlice4 fire\n', obliqueSlice4);                            % Scapula plane
             fprintf(fid,'ObliqueSlice4 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice4 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice5 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice5 fire\n', obliqueSlice5);                            % Scapula plane
             fprintf(fid,'ObliqueSlice5 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice5 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice6 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice6 fire\n', obliqueSlice6);                            % Scapula plane
             fprintf(fid,'ObliqueSlice6 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice6 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'remove CreateSphere%s HumSphere%s.surf SphereView%s Intersection\n', CTn, CTn, CTn);
             fprintf(fid,'remove CreateSphere%s_2 GlenoidSphere%s.surf SphereView%s_2 Intersection2\n', CTn, CTn, CTn);
             fprintf(fid,'create HxCreateSphere {CreateSphere%s}\n', CTn);
             fprintf(fid,'CreateSphere%s setIconPosition 20 500\n', CTn);
             fprintf(fid,'CreateSphere%s radius setMinMax 0 50\n', CTn);
             fprintf(fid,'CreateSphere%s radius setValue %f\n', CTn, HHRadius);
             if Side == 0
                 fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, -HC(1));
             else
                 fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, HC(1));
             end
             fprintf(fid,'CreateSphere%s coords setValue 1 %f\n', CTn, HC(2));
             fprintf(fid,'CreateSphere%s coords setValue 2 %f\n', CTn, HC(3));
             fprintf(fid,'CreateSphere%s fire\n', CTn);
             fprintf(fid,'[ {CreateSphere%s} create\n ] setLabel {HumSphere%s.surf}\n', CTn, CTn);
             fprintf(fid,'HumSphere%s.surf setIconPosition 200 500\n', CTn);
             fprintf(fid,'HumSphere%s.surf master connect CreateSphere%s\n', CTn, CTn);
             fprintf(fid,'HumSphere%s.surf fire\n', CTn);
             fprintf(fid,'create HxDisplaySurface {SphereView%s}\n', CTn);
             fprintf(fid,'SphereView%s setIconPosition 400 500\n', CTn);
             fprintf(fid,'SphereView%s data connect HumSphere%s.surf\n', CTn, CTn);
             fprintf(fid,'SphereView%s drawStyle setValue 4\n', CTn);
             fprintf(fid,'SphereView%s fire\n', CTn);
             fprintf(fid,'create HxOverlayGrid {Intersection}\n');
             fprintf(fid,'Intersection module connect ObliqueSlice4\n');
             fprintf(fid,'Intersection data connect HumSphere%s.surf\n', CTn);
             fprintf(fid,'Intersection lineWidth setValue 2\n');
             fprintf(fid,'ObliqueSlice4 setViewerMask 16383\n');
             fprintf(fid,'ObliqueSlice4 setPlane %f %f %f 1 0 0 0 1 0\n', HC);
             fprintf(fid,'ObliqueSlice4 fire\n\n');
 
             fprintf(fid,'create HxCreateSphere {CreateSphere%s_2}\n', CTn);
             fprintf(fid,'CreateSphere%s_2 setIconPosition 20 520\n', CTn);
             fprintf(fid,'CreateSphere%s_2 radius setMinMax 0 50\n', CTn);
             fprintf(fid,'CreateSphere%s_2 radius setValue %f\n', CTn, GlenoidRadius);
             if Side == 0
                 fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, -GlenoidSphere(1));
             else
                 fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, GlenoidSphere(1));
             end
             fprintf(fid,'CreateSphere%s_2 coords setValue 1 %f\n', CTn, GlenoidSphere(2));
             fprintf(fid,'CreateSphere%s_2 coords setValue 2 %f\n', CTn, GlenoidSphere(3));
             fprintf(fid,'CreateSphere%s_2 fire\n', CTn);
             fprintf(fid,'[ {CreateSphere%s_2} create\n ] setLabel {GlenoidSphere%s.surf}\n', CTn, CTn);
             fprintf(fid,'GlenoidSphere%s.surf setIconPosition 200 520\n', CTn);
             fprintf(fid,'GlenoidSphere%s.surf master connect CreateSphere%s_2\n', CTn, CTn);
             fprintf(fid,'GlenoidSphere%s.surf fire\n', CTn);
             fprintf(fid,'create HxDisplaySurface {SphereView%s_2}\n', CTn);
             fprintf(fid,'SphereView%s_2 setIconPosition 400 520\n', CTn);
             fprintf(fid,'SphereView%s_2 data connect GlenoidSphere%s.surf\n', CTn, CTn);
             fprintf(fid,'SphereView%s_2 drawStyle setValue 4\n', CTn);
             fprintf(fid,'SphereView%s_2 fire\n', CTn);
             fprintf(fid,'create HxOverlayGrid {Intersection2}\n');
             fprintf(fid,'Intersection2 module connect ObliqueSlice5\n');
             fprintf(fid,'Intersection2 data connect GlenoidSphere%s.surf\n', CTn);
             fprintf(fid,'Intersection2 lineWidth setValue 2\n');
             fprintf(fid,'ObliqueSlice5 setViewerMask 16383\n');
             fprintf(fid,'ObliqueSlice5 setPlane %f %f %f 1 0 0 0 1 0\n', GlenoidSphere);
             fprintf(fid,'ObliqueSlice5 fire\n');
 
             fclose(fid);
 
         case 'display'
 
             mkdir('Generated_Amira_TCL_Scripts/display');
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\display\\display_%s.tcl',CTn),'w');
             %              if Side == 0;
             %                  PlaneMean(1) = -PlaneMean(1);
             %              end
             obliqueSlice7 = [PlaneMean ScapulaCSYS(2,:)];
             if Side == 0
                 obliqueSlice7(1:3:end) = -obliqueSlice7(1:3:end);
                 GlenoidSphere(1) = -GlenoidSphere(1);
             end
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f \nObliqueSlice fire\n', obliqueSlice7);    % Muscle evaluation, through notch
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice fire\n\n');
 
             fprintf(fid,'create HxCreateSphere {GlenoidSphere%s}\n', CTn);
             fprintf(fid,'GlenoidSphere%s setIconPosition 20 550\n', CTn);
             fprintf(fid,'GlenoidSphere%s radius setMinMax 0 50\n', CTn);
             fprintf(fid,'GlenoidSphere%s radius setValue %f\n', CTn, GlenoidRadius);
             fprintf(fid,'GlenoidSphere%s coords setValue 0 %f\n', CTn, GlenoidSphere(1));
             fprintf(fid,'GlenoidSphere%s coords setValue 1 %f\n', CTn, GlenoidSphere(2));
             fprintf(fid,'GlenoidSphere%s coords setValue 2 %f\n', CTn, GlenoidSphere(3));
             fprintf(fid,'GlenoidSphere%s fire\n\n', CTn);
 
             fprintf(fid,'[ {GlenoidSphere%s} create ] setLabel {HumSphereGlenoid%s.surf}\n', CTn, CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf setIconPosition 200 550\n', CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf master connect GlenoidSphere%s\n', CTn, CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf fire\n\n', CTn);
 
             fprintf(fid,'create HxDisplaySurface {GlenoidSphereView%s}\n', CTn);
             fprintf(fid,'GlenoidSphereView%s setIconPosition 400 550\n', CTn);
             fprintf(fid,'GlenoidSphereView%s data connect HumSphereGlenoid%s.surf\n', CTn, CTn);
             fprintf(fid,'GlenoidSphereView%s drawStyle setValue 4\n', CTn);
             fprintf(fid,'GlenoidSphereView%s fire\n', CTn);
             fclose(fid);
         otherwise
             disp('unknown output request')
     end
 end