Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/anatomy/scapula_calculation.m
3761869a7f6dBashardoust_SystCoord

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
%
%%
% 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 = 2*HHRadius, humeral head diameter (radius*2)
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 = 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);
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
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'
% 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('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