Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/scapula_calculation.m
a30b07b1a173Bashardoust_SystCoord

scapula_calculation.m
View Options

			%% 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;

			%subject = char(subject);
			CTn = strtok(subject,'-');
			CTn = strtok(CTn, location);
			CTn = str2num(CTn);
			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
			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

			AG = landmarks(6,:); % We need the inital AG coordinates to use it in "density" part

			% 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)
			% AA: Angulus Acromialis (most laterodorsal point)
			% AG: Acromion-Glenoid notch

			% 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 ?


			%% 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);
			SP(:,1) = -SP(:,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
			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=find(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))';

			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
			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 = find(distance3D == min(distance3D));
			GC = p(ClosestNode,:)

			%% Fit Sphere on Glenoid Surface
			[GlenoidSphere,GlenoidRadius,GlenoidResiduals, ~] = fitSphere2(p);
			GlenoidSphere
			GlenoidRadius
			GlenoidRMSE = norm(GlenoidResiduals)/sqrt(length(GlenoidResiduals));

			%% Fit sphere on Humeral Head landmarks
			[HC,HHRadius,HHResiduals, ~] = fitSphere2(HH);
			HC
			HHRadius
			HHRMSE = norm(HHResiduals)/sqrt(length(HHResiduals));


			%% 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)];

			x1 = A\bx;

			by = [O(2);
			Ox(2);
			Oy(2);
			Oz(2)];

			x2 = A\by;

			bz = [O(3);
			Ox(3);
			Oy(3);
			Oz(3)];

			x3 = A\bz;

			B = [bx,by,bz];
			X = A\B;
			T = [X';0,0,0,1]; % Transformation matrix from CT frame to scapula frame

			[line_p, ~] = size(p);

			for i = 1 : 1 : line_p
			psc(i,:) = T*[p(i,:),1]';
			end

			psc_x = psc(:,1);
			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
			% GH: GlenoHumeral distance = 2HHRadius, humeral head diameter (radius2)

			AC = landmarks(4,:); % AC: Acromio-clavicular joint (most lateral part on acromion)

			ACinScapulaPlane = project2Plane(AC,ScapulaPlaneNormal',PlaneMean,size(AC,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')
			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))

			% 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';

			ACinScapulaCSYS = AC*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')
			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))

			% Project scapula, glenoid, AC, GC and scapula axis to scapula plane

			ACinScapulaPlane = project2Plane(ACinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ACinScapulaCSYS,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)

			GA = norm(cross(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:),ACinScapulaPlane-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')
			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))

			close all

			%% Output

			% Result.SubjectNumber = CTn;%1 A
			% Result.Version3D = Version3D;%2 B
			% Result.VersionDirection = VersionDirection;%3 C
			% Result.SHSubluxationIndex = SHSI;%4 D
			% Result.SHSubluxationDirection = SHSIDirection;%5 E
			% Result.GHSubluxationIndex = GHSI;%6 F
			% Result.GHSubluxationDirection = GHSIDirection;%7 G
			%
			% Result.GlenoidRadius = GlenoidRadius;%8 H
			% Result.HumeralHeadRadius = HHRadius;%9 I
			% Result.RadiusRatio = GlenoidRadius/HHRadius;%10 J
			% Result.GlenoidMedialization = 0;%11 K
			%
			% Result.Side = Side; %12 L 0:Right, 1:Left
			% Result.GrooveRMSE = GrooveRMSE;%13 M
			% Result.PlaneRMSE = PlaneRMSE;%14 N
			% Result.HHRMSE = HHRMSE;%15 O
			% Result.GlenoidRMSE = GlenoidRMSE;%16 P
			%
			% Result.R2Line = R2Line;%17 Q
			% Result.R2Plane = R2Plane;%18 R
			% Result.R2HH = R2HH;%19 S
			% Result.R2Glenoid = R2Glenoid;%20 T
			%
			% Result.CTorientation = CTorientation; %21 U
			% Result.WearPlaneAngle = WearPlaneAngle; %22 V
			% Result.Inclination = Inclination; %23 W
			% Result.Version2D = Version2D; % 24 X
			% Result.SHSPlaneAngle = SHSPlaneAngle; % 25 Y
			% Result.GHSPlaneAngle = GHSPlaneAngle; % 26 Z
			%
			% Result.CapHeight = CapHeight; %27 AB

			% Extra output
			CTtoXYZ = [TransverseAxis' ScapulaPlaneNormal cross(TransverseAxis, ScapulaPlaneNormal)'];
			Cxyz = (GC - Origin) * CTtoXYZ;

			% !! orientation of axes !! Now like Franke (2009) for right scapulae
			% Result.Cx = -Cxyz(2); %28 AA Ant-post
			% Result.Cy = Cxyz(3); %29 AB Sup-inf
			% Result.Cz = Cxyz(1); %30 AC Med-Lat

			% Scapula coordinate system in CT frame
			ScapulaCSYS(3,:) = TransverseAxis; % Z
			ScapulaCSYS(1,:) = ScapulaPlaneNormal; % X
			ScapulaCSYS(2,:) = cross(TransverseAxis,ScapulaPlaneNormal); % Y
			ScapulaCSYS(4,:) = Origin;

			% ISB coordinate system in CT frame

			% Real ISB
			% Zisb = (landmarks(4,:)-landmarks(2,:))/norm(landmarks(4,:)-landmarks(2,:)); %Med-Lat
			% Franke
			Zisb = (GC - landmarks(2,:)) / norm(GC - landmarks(2,:)); % GC - TS
			Xisb = cross(Zisb,landmarks(1,:)-landmarks(2,:))/norm(cross(Zisb,landmarks(1,:)-landmarks(2,:)));% Zisb ^cross^ (AI - TS)
			Yisb = cross(Zisb,Xisb);%Inf-Sup
			%
			% Result.X1 = ScapulaCSYS(1,1); % 30 AD
			% Result.X2 = ScapulaCSYS(1,2);
			% Result.X3 = ScapulaCSYS(1,3);
			%
			% Result.Y1 = ScapulaCSYS(2,1); % 33 AG
			% Result.Y2 = ScapulaCSYS(2,2);
			% Result.Y3 = ScapulaCSYS(2,3);
			%
			% Result.Z1 = ScapulaCSYS(3,1); % 36 AJ
			% Result.Z2 = ScapulaCSYS(3,2);
			% Result.Z3 = ScapulaCSYS(3,3);
			%
			% Result.Xisb1 = Xisb(1); % 39 AM
			% Result.Xisb2 = Xisb(2);
			% Result.Xisb3 = Xisb(3);
			%
			% Result.Yisb1 = Yisb(1); % 42 AP
			% Result.Yisb2 = Yisb(2);
			% Result.Yisb3 = Yisb(3);
			%
			% Result.Zisb1 = Zisb(1); % 45 AS
			% Result.Zisb2 = Zisb(2);
			% Result.Zisb3 = Zisb(3);

			rotX = vrrotvec(ScapulaCSYS(1,:),Xisb);
			% Result.deltaX = rad2deg(rotX(4)); % 48 AV
			rotY = vrrotvec(ScapulaCSYS(2,:),Yisb);
			% Result.deltaY = rad2deg(rotY(4)); % 49 AW
			rotZ = vrrotvec(ScapulaCSYS(3,:),Zisb);
			% Result.deltaZ = rad2deg(rotZ(4)); % 50 AX

			% Result.scapulaWidth = scapulaWidth;
			% Result.scapulaHeight = scapulaHeight;


			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 = Athx2abqAisb2thxAlbo2isb*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);

			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)];

			x1 = A\bx;

			by = [O(2);
			Ox(2);
			Oy(2);
			Oz(2)];

			x2 = A\by;

			bz = [O(3);
			Ox(3);
			Oy(3);
			Oz(3)];

			x3 = A\bz;

			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];
			R=inv(T);
			startcylpoint=R*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
			fid = fopen(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),'w');
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylR_factor, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_sphere_center, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_radius, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift5mm, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA2, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift15mmfront, '-append','precision',10);
			dlmwrite(sprintf('CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA3, '-append','precision',10);
			fclose(fid);

			fid = fopen(sprintf('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'

			% References data for a study case
			fid = fopen(sprintf('References\\References_%s.dat',CTn),'w');
			fprintf(fid,'# Scapula coordinate system\n\n');
			dlmwrite(sprintf('References\\References_%s.dat',CTn),ScapulaCSYS(4:7,:), '-append','precision',10);
			fclose(fid);
			fid = fopen(sprintf('References\\References_%s.dat',CTn),'a');
			fprintf(fid,'\n\n# Abaqus reference\n\n');
			dlmwrite(sprintf('References\\References_%s.dat',CTn),ScapulaABQ, '-append','precision',10);
			fclose(fid);
			fid = fopen(sprintf('References\\References_%s.dat',CTn),'a');
			fprintf(fid,'\n\n# Glenoid centerline\n\n');
			dlmwrite(sprintf('References\\References_%s.dat',CTn),GC, '-append','precision',10);
			dlmwrite(sprintf('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('obliqueSlice')
			fid = fopen(sprintf('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('display');
			fid = fopen(sprintf('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

R8218/scapula_calculation.ma30b07b1a173Bashardoust_SystCoord

scapula_calculation.mView Options

R8218/scapula_calculation.m
a30b07b1a173Bashardoust_SystCoord

scapula_calculation.m
View Options