scapula_calculation.m
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Jul 27, 14:22

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
	%
	%%
	% INPUTS
	% 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
	CTn = strtok(subject,'-')
	CTn = strtok(CTn, location)
	CTn = str2num(CTn) %#ok<ST2NM>
	CTn = sprintf('%s%03d',location,CTn)


	% Import Data
	ls

	% 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)
	% 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 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. 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 AL to glenoid principal axis
	% GH: GlenoHumeral distance = 2HHRadius, humeral head diameter (radius2)

	%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

	% 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

	% CSA calculation (work done by Bharath
	[radCSA, degCSA] = computeCSA(GlenoidPtsinScapulaPlane, ALinScapulaPlane);

	% 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(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])
	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'; % Create rotation matrix to transform points to Scapula Coordinate System

	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(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])
	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 in the
	% Scapula Coordinate System

	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 AL to glenoid principal axis)
	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(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) ')'])
	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.glenoidSphericity = ' ';
	%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 = degCSA % radCSA;
	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;

	Result


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

	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

scapula_calculation.mNo OneTemporaryActions

File Metadata

scapula_calculation.mView Options

Event Timeline

scapula_calculation.m
No OneTemporary
Actions

scapula_calculation.m
View Options