Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (feat_muscle_subdivisions) Browse

R8218/anatomy/scapula_calculation.m
d117fa36fbd6feat_muscle_subdivisions

Last Change
Raw File

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 = 2*HHRadius, humeral head diameter (radius*2)
%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 = Athx2abq*Aisb2thx*Albo2isb*Act2lbo;
%% Output initial position of scapula for abaqus simulation in CT coordinates
% scaleFactor = 24/HHRadius; % =24mm / Humerus radius
ScapulaABQ(1,:) = HC; % Humerus center
ScapulaABQ(2,:) = 10*Act2abq(:,1)'+ScapulaABQ(1,:);
ScapulaABQ(3,:) = 10*Act2abq(:,2)'+ScapulaABQ(1,:);
ScapulaABQ(4,:) = 10*Act2abq(:,3)'+ScapulaABQ(1,:);
% ScapulaABQ(5,1) = scaleFactor;
%% Output coordinates for obliqueSlice
% The ObliqueSlice can be defined as point + (normal OR u-vector and
% v-vector)
FriedmanPlane = [GC 1 0 0 0 1 0];
obliqueSlice = [landmarks(6,:) ScapulaCSYS(2,:) -ScapulaCSYS(3,:)]; % Plane for muscle measurement
obliqueSlice2 = [GC ScapulaCSYS(3,:)]; % Scapula axial plane
obliqueSlice3 = [GC ScapulaCSYS(1,:) GO]; % Max wear plane
obliqueSlice4 = [GC ScapulaCSYS(2,:)]; % Scapular plane
obliqueSlice5 = [GC HO/norm(HO) GHeccentricity/norm(GHeccentricity)]; % Plane for GHS
obliqueSlice6 = [GC HO/norm(HO) SHeccentricity/norm(SHeccentricity)]; % Plane for SHS
if Side == 0
FriedmanPlane(1) = -FriedmanPlane(1);
obliqueSlice(1:3:end) = - obliqueSlice(1:3:end);
obliqueSlice2(1:3:end) = -obliqueSlice2(1:3:end);
obliqueSlice3(1:3:end) = -obliqueSlice3(1:3:end);
obliqueSlice4(1:3:end) = -obliqueSlice4(1:3:end);
obliqueSlice5(1:3:end) = -obliqueSlice5(1:3:end);
obliqueSlice6(1:3:end) = -obliqueSlice6(1:3:end);
end
switch output
case 'density'
%% Data for the calculations of the density in amira
% Needed variables are :
% - CylO: cylinder origin
% - CylA: cylinder alignment point
% - CylR_factor: scale factor for cylinder radius
% - Glenoid_sphere_center
% - Glenoid radius
% - CylO_shift5mm : cylinder origin shifted 5 mm behind
% glenoid surface center
% - CylA2 : cylinder alignment point n?2
% - CylO_shift15mmfront: cylinder origin shifted 15 mm
% frontwards
% - CylA3 : cylinder alignemnt point n?2
% - startcylpoint : origin of the new cylinder starting from
% the AG point
% - endcylpoint: end of cylinder, 40 mm from the startcylpoint
% in the scapula axis
CylO = GC;
Glenoid_sphere_center = GlenoidSphere';
Glenoid_radius = GlenoidRadius;
%% CylA - second point to form a line with glenoid surface
% center point which is parallel to scapula axis
CylA = CylO - 40 * TransverseAxis/norm(TransverseAxis);
%% CylR - distance from glenoid surface center to most
% eccentric point of glenoid surface
%Projection of surface points onto the same plane
PlanePoint = CylO;
PlaneNormal = (CylO-CylA)/norm(CylO-CylA);
ProjectedPoints = zeros(size(p));
DistanceToCylO = zeros(size(p,1),1);
for i = 1:size(p,1)
ProjectedPoints(i,:) = p(i,:) - dot(p(i,:) - PlanePoint, PlaneNormal) * PlaneNormal;
DistanceToCylO(i) = norm(ProjectedPoints(i,:) - CylO);
end
%Selection of the most eccentric point -> radius of cylinder
CylR = max(DistanceToCylO);
%Scale factor of the cylinder
OriginalCylR = 5;
CylR_factor = CylR/OriginalCylR;
% Coordinates of start and end point of the D10 mm cylinder
Deepness = 5; %behind glenoid surface
smallCylR = 5;
CylO_shift5mm = GC - Deepness * TransverseAxis/norm(TransverseAxis);
CylA2 = CylO_shift5mm - smallCylR * TransverseAxis;
% To place big cylinder (finer seleciton, exclusion of distant
% bones)
CylO_shift15mmfront = GC + 15 * TransverseAxis/norm(TransverseAxis);
CylA3 = CylO_shift15mmfront - 40 * TransverseAxis/norm(TransverseAxis);
%%Projection of AG on scapula axis
AG = landmarks(6,:);
p1 = CylA3;
p2 = CylO_shift15mmfront;
v1 = p2-p1;
v2 = AG-p1;
E1 = v1;
E2 = cross(v1,v2);
E3 = cross(E1,E2);
e1 = E1/norm(E1);
e2 = E2/norm(E2);
e3 = E3/norm(E3);
Oc = p1;
Ocx = p1 + e1;
Ocy = p1 + e2;
Ocz = p1 + e3;
O = [0,0,0];
Ox = [1,0,0];
Oy = [0,1,0];
Oz = [0,0,1];
A = [Oc(1),Oc(2),Oc(3),1;
Ocx(1),Ocx(2),Ocx(3),1;
Ocy(1),Ocy(2),Ocy(3),1;
Ocz(1),Ocz(2),Ocz(3),1];
bx = [O(1);
Ox(1);
Oy(1);
Oz(1)];
by = [O(2);
Ox(2);
Oy(2);
Oz(2)];
bz = [O(3);
Ox(3);
Oy(3);
Oz(3)];
B = [bx,by,bz];
X = A\B;
T = [X';
0,0,0,1];
det(T);
AG2 = [AG,1];
AG3=AG2';
diff = T*AG3;
projAGcyl=[diff(1),0,0,1];
startcylpoint=T\projAGcyl';
startcylpoint = [startcylpoint(1),startcylpoint(2),startcylpoint(3)];
endcylpoint = startcylpoint + 40 * TransverseAxis/norm(TransverseAxis);
%% Write data
if Side == 0
CylO(1) = -CylO(1);
CylA(1) = -CylA(1);
Glenoid_sphere_center(1) = -Glenoid_sphere_center(1);
CylO_shift5mm(1) = -CylO_shift5mm(1);
CylA2(1) = -CylA2(1);
CylO_shift15mmfront(1) = -CylO_shift15mmfront(1);
CylA3(1) = -CylA3(1);
startcylpoint(1) = -startcylpoint(1);
endcylpoint(1) = -endcylpoint(1);
end
mkdir('Generated_Amira_TCL_Scripts/CylinderReferences')
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),'w');
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylR_factor, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_sphere_center, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_radius, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift5mm, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA2, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift15mmfront, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA3, '-append','precision',10);
fclose(fid);
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s_AmiraCode.dat',CTn),'w');
fprintf(fid,'set lineset1 [create HxLineSet]\n');
fprintf(fid,'set p1 [$lineset1 addPoint %f %f %f]\n',CylO_shift5mm);
fprintf(fid,'set p2 [$lineset1 addPoint %f %f %f]\n',CylA2);
fprintf(fid,'$lineset1 addLine $p1 $p2\n');
fprintf(fid,'{Line-Set} setIconPosition 20 10\n');
fprintf(fid,'Line-Set fire\n\n');
fprintf(fid,'set lineset2 [create HxLineSet]\n');
fprintf(fid,'set p1 [$lineset2 addPoint 0 0 0]\n');
fprintf(fid,'set p2 [$lineset2 addPoint 5 0 0]\n');
fprintf(fid,'$lineset2 addLine $p1 $p2\n');
fprintf(fid,'{Line-Set2} setIconPosition 20 30\n');
fprintf(fid,'Line-Set2 fire\n\n');
fprintf(fid,'set hideNewModules 0\n');
fprintf(fid,'create {HxCreateSphere} {GlenoidSphereFit}\n');
fprintf(fid,'{GlenoidSphereFit} setIconPosition 20 330\n');
fprintf(fid,'{GlenoidSphereFit} fire\n');
fprintf(fid,'{GlenoidSphereFit} {type} setValue 0\n');
fprintf(fid,'{GlenoidSphereFit} {radius} setMinMax 0 100\n');
fprintf(fid,'{GlenoidSphereFit} {radius} setButtons 0\n');
fprintf(fid,'{GlenoidSphereFit} {radius} setIncrement 0.666667\n');
fprintf(fid,'{GlenoidSphereFit} {radius} setValue %f\n',Glenoid_radius);
fprintf(fid,'{GlenoidSphereFit} {radius} setSubMinMax 0 100\n');
fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{GlenoidSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{GlenoidSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{GlenoidSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
fprintf(fid,'{GlenoidSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
fprintf(fid,'{GlenoidSphereFit} {edgeLength} setButtons 0\n');
fprintf(fid,'{GlenoidSphereFit} {edgeLength} setIncrement 0.332667\n');
fprintf(fid,'{GlenoidSphereFit} {edgeLength} setValue 1.07457\n');
fprintf(fid,'{GlenoidSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setButtons 0\n');
fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setValue 1\n');
fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
fprintf(fid,'GlenoidSphereFit fire\n');
fprintf(fid,'GlenoidSphereFit setViewerMask 16383\n');
fprintf(fid,'GlenoidSphereFit setPickable 1\n\n');
fprintf(fid,'set hideNewModules 0\n');
fprintf(fid,'create {HxCreateSphere} {CorticalSphereFit}\n');
fprintf(fid,'{CorticalSphereFit} setIconPosition 20 350\n');
fprintf(fid,'{CorticalSphereFit} fire\n');
fprintf(fid,'{CorticalSphereFit} {type} setValue 0\n');
fprintf(fid,'{CorticalSphereFit} {radius} setMinMax 0 100\n');
fprintf(fid,'{CorticalSphereFit} {radius} setButtons 0\n');
fprintf(fid,'{CorticalSphereFit} {radius} setIncrement 0.666667\n');
fprintf(fid,'{CorticalSphereFit} {radius} setValue %f\n',Glenoid_radius+3);
fprintf(fid,'{CorticalSphereFit} {radius} setSubMinMax 0 100\n');
fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{CorticalSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{CorticalSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
fprintf(fid,'{CorticalSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
fprintf(fid,'{CorticalSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
fprintf(fid,'{CorticalSphereFit} {edgeLength} setButtons 0\n');
fprintf(fid,'{CorticalSphereFit} {edgeLength} setIncrement 0.332667\n');
fprintf(fid,'{CorticalSphereFit} {edgeLength} setValue 1.07457\n');
fprintf(fid,'{CorticalSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setButtons 0\n');
fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setValue 1\n');
fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
fprintf(fid,'CorticalSphereFit fire\n');
fprintf(fid,'CorticalSphereFit setViewerMask 16383\n');
fprintf(fid,'CorticalSphereFit setPickable 1\n\n');
fprintf(fid,'set lineset3 [create HxLineSet]\n');
fprintf(fid,'set p1 [$lineset3 addPoint %f %f %f]\n',startcylpoint);
fprintf(fid,'set p2 [$lineset3 addPoint %f %f %f]\n',endcylpoint);
fprintf(fid,'$lineset3 addLine $p1 $p2\n');
fprintf(fid,'{Line-Set3} setIconPosition 20 520\n');
fprintf(fid,'Line-Set3 fire\n\n');
fprintf(fid,'set lineset4 [create HxLineSet]\n');
fprintf(fid,'set p1 [$lineset4 addPoint 0 0 0]\n');
fprintf(fid,'set p2 [$lineset4 addPoint 40 0 0]\n');
fprintf(fid,'$lineset4 addLine $p1 $p2\n');
fprintf(fid,'{Line-Set4} setIconPosition 20 540\n');
fprintf(fid,'Line-Set4 fire\n\n');
fprintf(fid,'{Cylinder2.STL} setScaleFactor 1 %f %f\n',CylR_factor,CylR_factor);
fprintf(fid,'{Cylinder2.STL} applyTransform\n\n');
fprintf(fid,'saveProject\n');
fprintf(fid,'clear\n');
fprintf(fid,'echo "Complete center coordinates of ReamingSphere"\n\n');
%fprintf(fid,'
fclose(fid);
case 'References'
mkdir('Generated_Amira_TCL_Scripts/References')
% References data for a study case
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'w');
fprintf(fid,'# Scapula coordinate system\n\n');
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaCSYS(4:7,:), '-append','precision',10);
fclose(fid);
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
fprintf(fid,'\n\n# Abaqus reference\n\n');
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaABQ, '-append','precision',10);
fclose(fid);
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
fprintf(fid,'\n\n# Glenoid centerline\n\n');
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GC, '-append','precision',10);
dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GlenoidSphere', '-append','precision',10);
fclose(fid);
case 'obliqueSlice'
[token, remain] = strtok(directory,'../../..');
directory = ['Z:/' token remain];
pictureName = [CTn '_new'];
mkdir('Generated_Amira_TCL_Scripts/obliqueSlice')
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\obliqueslice\\obliqueSlice_%s.tcl',CTn),'w');
fprintf(fid,'# Create an ObliqueSlice with Friedman plane #\n');
fprintf(fid,'ObliqueSlice orientation hit 0\n');
fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\n',FriedmanPlane);
fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 0 item 2 1 item 3 0 item 4 0\n');
fprintf(fid,'ObliqueSlice fire\n');
fprintf(fid,'set planepos [ObliqueSlice translate getValue]\n');
fprintf(fid,'set maxvalue [ObliqueSlice translate getMaxValue]\n');
fprintf(fid,'set planepos [expr {$maxvalue - int($planepos)}]\n');
fprintf(fid,'ObliqueSlice createImage %s_$planepos\n', CTn);
fprintf(fid,'create HxCastField {CastField}\n');
fprintf(fid,'CastField data connect %s_$planepos\n', CTn);
fprintf(fid,'CastField fire\n');
fprintf(fid,'CastField scaling setValue 0 0.21\n');
fprintf(fid,'CastField scaling setValue 1 200\n');
fprintf(fid,'CastField fire\n');
fprintf(fid,'[ {CastField} create ] setLabel %s_$planepos.jpeg\n', CTn);
if strcmp(location, 'P')
fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" %s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
else
fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" Z:/data/%s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
end
fprintf(fid,'# Get image for muscle measurement #\n');
fprintf(fid,'ObliqueSlice orientation hit 0\n');
fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice fire\n', obliqueSlice); % Muscle evaluation, through notch
fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
fprintf(fid,'ObliqueSlice createImage {%s}\n', pictureName);
if strcmp(location, 'pathologic')
fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
else
fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
end
fprintf(fid,'saveProject\n\n');
fprintf(fid,'ObliqueSlice2 setPlane %f %f %f %f %f %f\nObliqueSlice2 fire\n', obliqueSlice2); % Scapular axial plane
fprintf(fid,'ObliqueSlice2 options setValue 0 1\n');
fprintf(fid,'ObliqueSlice2 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
fprintf(fid,'ObliqueSlice3 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice3 fire\n', obliqueSlice3); % Max wear plane
fprintf(fid,'ObliqueSlice3 options setValue 0 1\n');
fprintf(fid,'ObliqueSlice3 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
fprintf(fid,'ObliqueSlice4 setPlane %f %f %f %f %f %f\nObliqueSlice4 fire\n', obliqueSlice4); % Scapula plane
fprintf(fid,'ObliqueSlice4 options setValue 0 1\n');
fprintf(fid,'ObliqueSlice4 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
fprintf(fid,'ObliqueSlice5 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice5 fire\n', obliqueSlice5); % Scapula plane
fprintf(fid,'ObliqueSlice5 options setValue 0 1\n');
fprintf(fid,'ObliqueSlice5 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
fprintf(fid,'ObliqueSlice6 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice6 fire\n', obliqueSlice6); % Scapula plane
fprintf(fid,'ObliqueSlice6 options setValue 0 1\n');
fprintf(fid,'ObliqueSlice6 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
fprintf(fid,'remove CreateSphere%s HumSphere%s.surf SphereView%s Intersection\n', CTn, CTn, CTn);
fprintf(fid,'remove CreateSphere%s_2 GlenoidSphere%s.surf SphereView%s_2 Intersection2\n', CTn, CTn, CTn);
fprintf(fid,'create HxCreateSphere {CreateSphere%s}\n', CTn);
fprintf(fid,'CreateSphere%s setIconPosition 20 500\n', CTn);
fprintf(fid,'CreateSphere%s radius setMinMax 0 50\n', CTn);
fprintf(fid,'CreateSphere%s radius setValue %f\n', CTn, HHRadius);
if Side == 0
fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, -HC(1));
else
fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, HC(1));
end
fprintf(fid,'CreateSphere%s coords setValue 1 %f\n', CTn, HC(2));
fprintf(fid,'CreateSphere%s coords setValue 2 %f\n', CTn, HC(3));
fprintf(fid,'CreateSphere%s fire\n', CTn);
fprintf(fid,'[ {CreateSphere%s} create\n ] setLabel {HumSphere%s.surf}\n', CTn, CTn);
fprintf(fid,'HumSphere%s.surf setIconPosition 200 500\n', CTn);
fprintf(fid,'HumSphere%s.surf master connect CreateSphere%s\n', CTn, CTn);
fprintf(fid,'HumSphere%s.surf fire\n', CTn);
fprintf(fid,'create HxDisplaySurface {SphereView%s}\n', CTn);
fprintf(fid,'SphereView%s setIconPosition 400 500\n', CTn);
fprintf(fid,'SphereView%s data connect HumSphere%s.surf\n', CTn, CTn);
fprintf(fid,'SphereView%s drawStyle setValue 4\n', CTn);
fprintf(fid,'SphereView%s fire\n', CTn);
fprintf(fid,'create HxOverlayGrid {Intersection}\n');
fprintf(fid,'Intersection module connect ObliqueSlice4\n');
fprintf(fid,'Intersection data connect HumSphere%s.surf\n', CTn);
fprintf(fid,'Intersection lineWidth setValue 2\n');
fprintf(fid,'ObliqueSlice4 setViewerMask 16383\n');
fprintf(fid,'ObliqueSlice4 setPlane %f %f %f 1 0 0 0 1 0\n', HC);
fprintf(fid,'ObliqueSlice4 fire\n\n');
fprintf(fid,'create HxCreateSphere {CreateSphere%s_2}\n', CTn);
fprintf(fid,'CreateSphere%s_2 setIconPosition 20 520\n', CTn);
fprintf(fid,'CreateSphere%s_2 radius setMinMax 0 50\n', CTn);
fprintf(fid,'CreateSphere%s_2 radius setValue %f\n', CTn, GlenoidRadius);
if Side == 0
fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, -GlenoidSphere(1));
else
fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, GlenoidSphere(1));
end
fprintf(fid,'CreateSphere%s_2 coords setValue 1 %f\n', CTn, GlenoidSphere(2));
fprintf(fid,'CreateSphere%s_2 coords setValue 2 %f\n', CTn, GlenoidSphere(3));
fprintf(fid,'CreateSphere%s_2 fire\n', CTn);
fprintf(fid,'[ {CreateSphere%s_2} create\n ] setLabel {GlenoidSphere%s.surf}\n', CTn, CTn);
fprintf(fid,'GlenoidSphere%s.surf setIconPosition 200 520\n', CTn);
fprintf(fid,'GlenoidSphere%s.surf master connect CreateSphere%s_2\n', CTn, CTn);
fprintf(fid,'GlenoidSphere%s.surf fire\n', CTn);
fprintf(fid,'create HxDisplaySurface {SphereView%s_2}\n', CTn);
fprintf(fid,'SphereView%s_2 setIconPosition 400 520\n', CTn);
fprintf(fid,'SphereView%s_2 data connect GlenoidSphere%s.surf\n', CTn, CTn);
fprintf(fid,'SphereView%s_2 drawStyle setValue 4\n', CTn);
fprintf(fid,'SphereView%s_2 fire\n', CTn);
fprintf(fid,'create HxOverlayGrid {Intersection2}\n');
fprintf(fid,'Intersection2 module connect ObliqueSlice5\n');
fprintf(fid,'Intersection2 data connect GlenoidSphere%s.surf\n', CTn);
fprintf(fid,'Intersection2 lineWidth setValue 2\n');
fprintf(fid,'ObliqueSlice5 setViewerMask 16383\n');
fprintf(fid,'ObliqueSlice5 setPlane %f %f %f 1 0 0 0 1 0\n', GlenoidSphere);
fprintf(fid,'ObliqueSlice5 fire\n');
fclose(fid);
case 'display'
mkdir('Generated_Amira_TCL_Scripts/display');
fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\display\\display_%s.tcl',CTn),'w');
% if Side == 0;
% PlaneMean(1) = -PlaneMean(1);
% end
obliqueSlice7 = [PlaneMean ScapulaCSYS(2,:)];
if Side == 0
obliqueSlice7(1:3:end) = -obliqueSlice7(1:3:end);
GlenoidSphere(1) = -GlenoidSphere(1);
end
fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f \nObliqueSlice fire\n', obliqueSlice7); % Muscle evaluation, through notch
fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
fprintf(fid,'ObliqueSlice fire\n\n');
fprintf(fid,'create HxCreateSphere {GlenoidSphere%s}\n', CTn);
fprintf(fid,'GlenoidSphere%s setIconPosition 20 550\n', CTn);
fprintf(fid,'GlenoidSphere%s radius setMinMax 0 50\n', CTn);
fprintf(fid,'GlenoidSphere%s radius setValue %f\n', CTn, GlenoidRadius);
fprintf(fid,'GlenoidSphere%s coords setValue 0 %f\n', CTn, GlenoidSphere(1));
fprintf(fid,'GlenoidSphere%s coords setValue 1 %f\n', CTn, GlenoidSphere(2));
fprintf(fid,'GlenoidSphere%s coords setValue 2 %f\n', CTn, GlenoidSphere(3));
fprintf(fid,'GlenoidSphere%s fire\n\n', CTn);
fprintf(fid,'[ {GlenoidSphere%s} create ] setLabel {HumSphereGlenoid%s.surf}\n', CTn, CTn);
fprintf(fid,'HumSphereGlenoid%s.surf setIconPosition 200 550\n', CTn);
fprintf(fid,'HumSphereGlenoid%s.surf master connect GlenoidSphere%s\n', CTn, CTn);
fprintf(fid,'HumSphereGlenoid%s.surf fire\n\n', CTn);
fprintf(fid,'create HxDisplaySurface {GlenoidSphereView%s}\n', CTn);
fprintf(fid,'GlenoidSphereView%s setIconPosition 400 550\n', CTn);
fprintf(fid,'GlenoidSphereView%s data connect HumSphereGlenoid%s.surf\n', CTn, CTn);
fprintf(fid,'GlenoidSphereView%s drawStyle setValue 4\n', CTn);
fprintf(fid,'GlenoidSphereView%s fire\n', CTn);
fclose(fid);
otherwise
disp('unknown output request')
end
end
c4science · Help