Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

R8218/scapula_measure.m
56d95df3d151Bashardoust_SystCoord

Last Change
Raw File

scapula_measure.m
View Options

%% scapula_measure
%
% This script returns the anatomy measurements from the Amira files in an
% excel file.
% The file name is 'OutputAll' for type 'all', 'Output' for type 'normal'
% or 'pathologic'. For the 'specific' measure, the file name is 'Specific_scapula_measure.xls'
% The script makes measurements for patients from the "unmeasured-cases" list or for specific patients.
% Author: EPFL-LBO
% Date: 2016-09-05
%
%%
% Source is 'list' or 'specific' (for only one case)
% Type is 'normal' or 'pathologic' or 'all'
% Subject is the patient number
function scapula_measure(source,type,subject)
% Define the directory and inmitialize n and m
directory = '../../../data';
n=1;
m=1;
switch source
case 'list'
switch type
%% Creation of the case list according to the type: 'normal, pathologic or all
case 'pathologic'
location = 'P';
l=1;
rootDir = '../../../data/P';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/P*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
sCases{l,1} = folderInfo(i).name;
sCases{l,2} = levelDir1;
sCases{l,3} = levelDir2;
l=l+1;
end
end
end
case 'normal'
location = 'N';
l=1;
rootDir = '../../../data/N';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/N*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
sCases{l,1} = folderInfo(i).name;
sCases{l,2} = levelDir1;
sCases{l,3} = levelDir2;
l=l+1;
end
end
end
case 'all'
location = 'N';
l=1;
rootDir = '../../../data/N';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/N*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
sCases_normal{l,1} = folderInfo(i).name;
sCases_normal{l,2} = levelDir1;
sCases_normal{l,3} = levelDir2;
sCases_normal{l,4} = 'N';
l=l+1;
end
end
end
location = 'P';
l=1;
rootDir = '../../../data/P';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/P*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
sCases_pathologic{l,1} = folderInfo(i).name;
sCases_pathologic{l,2} = levelDir1;
sCases_pathologic{l,3} = levelDir2;
sCases_pathologic{l,4} = 'P';
l=l+1;
end
end
end
sCases = [sCases_normal;sCases_pathologic];
otherwise
error('Unexpected shoulder type: "normal" or "pathologic"');
end
%% Loop over the all cases
if strcmp(type,'all') == 1
%Loop over the normal cases
for k = 2 : length(sCases)+1
patient = char(sCases(k-1,1));
levelDir1 = cell2mat(sCases(k-1,2));
levelDir2 = cell2mat(sCases(k-1,3));
location = char(sCases(k-1,4));
finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
if exist(sprintf('%s/%s/CT-%s-1/amira', finalDirectory,patient,patient),'dir') == 7 %if a folder "amira" exist, the patient name is added to the "measured" column
result = scapula_calculation(patient, finalDirectory, location); % Function that calculates the anatomical datas
vars = fieldnames(result);
if k==2
for m=1:length(vars) %Loop to write the fieldname in the first line of the excel file
Output{k-1,m}= vars{m};
end
end
for m = 1:length(vars)
Output{k,m} = result.(vars{m}); %Loop to write the data in the excel file
end
else
Output{k,1} = strtok(patient,'-');
end
Progress = 100*k/length(sCases); %Display the percentage progress
clc
fprintf('All %3.0f%%\n ',Progress)
fprintf('All %d on %d\n ',k-1 , length(sCases))
end
% %Loop over the pathologic cases
% for k = length(case_list_normal)+2 : length(case_list)+1
% result = scapula_calculation(case_list(k-1), directory, 'Pathologic'); % Function that calculates the anatomical datas
% vars = fieldnames(result);
%
% for m = 1:length(vars)
% Output{k,m} = result.(vars{m}); %Loop to write the data in the excel file
% end
%
% PProgress = 100*k/length(case_list); %Display the percentage progress
% clc
% fprintf('Pathologic %3.0f%%\n ',PProgress)
% fprintf('Pathologic %d on %d',k ,length(case_list_pathologic))
% end
delete('anatomy.xls');
xlswrite('anatomy.xls',Output);
end
%% Loop over the Normal cases OR the pathologic cases according to the type
if strcmp(type,'all') == 0
for j = 2:1:length(sCases)+1 % Loop over the list
patient = char(sCases(j-1,1));
levelDir1 = cell2mat(sCases(j-1,2));
levelDir2 = cell2mat(sCases(j-1,3));
finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
if exist(sprintf('%s/%s/CT-%s-1/amira', finalDirectory,patient,patient),'dir') == 7 %if a folder "amira" exist, the patient name is added to the "measured" column
result = scapula_calculation(patient, finalDirectory, location); % Function that calculates the anatomical datas
vars = fieldnames(result);
if j==2
for m=1:length(vars) %Loop to write the fieldname in the first line of the excel file
Output{j-1,m}= vars{m};
end
end
for m = 1:length(vars) %Loop to write the data in the excel file
Output{j,m} = result.(vars{m});
end
else
Output{j,1} = strtok(patient,'-');
end
if strcmp(type,'normal') == 1
Nprogress = 100*j/length(sCases); %Display the percentage progress for normal cases
clc
fprintf('Normal %3.0f%%\n',Nprogress)
fprintf('Normal %d on %d ',j-1 , length(sCases))
end
if strcmp(type,'pathologic') == 1
Pprogress = 100*j/length(sCases); %Display the percentage progress for pathologic cases
clc
fprintf('Pathologic %3.0f%%\n' ,Pprogress)
fprintf('Pathologic %d on %d',j-1 ,length(sCases))
end
end
delete('Output.xls');
xlswrite('Output.xls',Output); %Create the excel file
end
case 'specific'
directory = '../../../data';
switch type
case 'pathologic'
location = 'P';
l=1;
rootDir = '../../../data/P';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/P*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
if exist(sprintf('%s/%d/%d/%s/CT-%s-1/amira', rootDir,levelDir1,levelDir2,name,name),'dir') == 7 %if a folder "amira" exist, the patient name is added to the "measured" column
sCases{l,1} = folderInfo(i).name;
sCases{l,2} = levelDir1;
sCases{l,3} = levelDir2;
l=l+1;
end
end
end
end
CTn = sprintf('P%d-',subject); % format the patient name to P###
for i = 1:1:numel(sCases(:,1)) % Loop over the list of patients and extract the complete name of the specific patient
name = char(sCases(i,1));
t = strfind(name,CTn);
if t == 1;
patient = name;
levelDir1 = cell2mat(sCases(i,2));
levelDir2 = cell2mat(sCases(i,3));
end
end
finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
result = scapula_calculation(patient, finalDirectory, location); % Function that calculates the anatomical datas
vars = fieldnames(result);
for m=1:length(vars) %Loop to write the fieldname in the first line of the excel file
Output{1,m}= vars{m};
end
for n = 1:length(vars) %Loop to write the data in the excel file
Output{2,n} = result.(vars{n});
end
delete ('Specific_scapula_measure.xls');
xlswrite('Specific_scapula_measure.xls',Output); %Create the excel file
case 'normal'
location = 'N';
l=1;
rootDir = '../../../data/N';
for levelDir1 = 0:9
for levelDir2 = 0:9
curentDir = sprintf('%s/%d/%d/N*', rootDir, levelDir1, levelDir2);
folderInfo = dir(curentDir);
for i=1:1:numel(folderInfo)
name = folderInfo(i).name;
if exist(sprintf('%s/%d/%d/%s/CT-%s-1/amira', rootDir,levelDir1,levelDir2,name,name),'dir') == 7 %if a folder "amira" exist, the patient name is added to the "measured" column
sCases{l,1} = folderInfo(i).name;
sCases{l,2} = levelDir1;
sCases{l,3} = levelDir2;
l=l+1;
end
end
end
end
CTn = sprintf('N%d-',subject); % format the patient name to P###
%CTn = sprintf('N%03d',subject); % format the patient name to P###
% patient = sprintf('N%03d',subject);
% subject = subject+1000; % format the patient name to N###
%normal_list = dir(sprintf('%s/%s/N*', directory, location));
for i = 1:1:numel(sCases(:,1)) % Loop over the list of patients and extract the complete name of the specific patient
name = char(sCases(i,1));
t = strfind(name,CTn);
if t == 1;
patient = name;
levelDir1 = cell2mat(sCases(i,2));
levelDir2 = cell2mat(sCases(i,3));
end
end
finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
result = scapula_calculation(patient, finalDirectory, location); % Function that calculates the anatomical datas
vars = fieldnames(result);
for m=1:length(vars) %Loop to write the fieldname in the first line of the excel file
Output{1,m}= vars{m};
end
for n = 1:length(vars) %Loop to write the data in the excel file
Output{2,n} = result.(vars{n});
end
delete ('Specific_scapula_measure.xls');
xlswrite('Specific_scapula_measure.xls',Output); %Create the excel file
end
end
%% This function returns all the measurments
% function [Result] = anatomy(subject,directory, location)
%
% patient_num = char(subject);
% CTn = strtok(patient_num,'-');
%
% %% Import Data
%
%
% % Import Glenoid Surface
% if exist(sprintf('%s/%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,location,patient_num,patient_num,CTn),'file') == 2
% [p,~,~]=importStl(sprintf('%s/%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,location,patient_num,patient_num,CTn),1);
% else
% error('MATLAB:rmpath:DirNotFound',...
% 'Could not find glenoid surface file: %s/%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,location,patient_num,patient_num,CTn)
% end
%
%
%
% %%Import HumeralHead landmarks
% if exist(sprintf('%s/%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn),'file') == 2
% newData1 = importdata(sprintf('%s/%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn), ' ', 14);
% HH = newData1.data;
% else
% warning('Using old humeral landmarks definition')
% if exist(sprintf('%s/%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn),'file') == 2
% newData1 = importdata(sprintf('%s/%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn), ' ', 14);
% HH = newData1.data;
% else
% error('MATLAB:rmpath:DirNotFound',...
% 'Could not find humeral head landmarks file: %s/%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn)
% end
% end
%
%
% % Import Scapula landmarks
% if exist(sprintf('%s/%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn),'file') == 2
% newData2 = importdata(sprintf('%s/%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn), ' ', 14);
% landmarks = newData2.data;
% else
% error('MATLAB:rmpath:DirNotFound',...
% 'Could not find scapula landmarks file: %s/%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,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/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn),'file') == 2
% newData3 = importdata(sprintf('%s/%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn), ' ', 14);
% SP = newData3.data; %Scapula Pillar points
% else
% error('MATLAB:rmpath:DirNotFound',...
% 'Could not find axillary border landmarks file: %s/%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn)
% end
%
% % Import Scapula Groove
% if exist(sprintf('%s/%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn),'file') == 2
% newData4 = importdata(sprintf('%s/%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,CTn), ' ', 14);
% SG = newData4.data; %Scapula Groove points
% else
% error('MATLAB:rmpath:DirNotFound',...
% 'Could not find supraspinatus fossa landmarks file: %s/%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,location,patient_num,patient_num,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);
% 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 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, ~, GrooveRMSE, R2Line] = 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, R2Glenoid] = fitSphere2(p);
% GlenoidSphere
% GlenoidRadius
% GlenoidRMSE = norm(GlenoidResiduals)/sqrt(length(GlenoidResiduals));
%
% %% Fit sphere on Humeral Head landmarks
% [HC,HHRadius,HHResiduals, R2HH] = 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, column_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,:));
%
% %% 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;
%
% end
end
c4science · Help