Homec4science
Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (aterrier-plotUpdate) Browse

R8218/ShoulderCase/ShoulderCase.m_bk
98c701372cbdaterrier-plotUpdate

Last Change
Raw File

ShoulderCase.m_bk
View Options

classdef ShoulderCase < handle
% Properties and methods associated to the SoulderCase object.
% Author: Alexandre Terrier, EPFL-LBO
% Creation date: 2018-07-01
% Revision date: 2018-08-11
% TODO:
% Need method to load properties (from amira, matlab, excel, SQL)
% Need method to save properties (in matlab, excel, SQL)
% Remove datapath from constants
%
properties
id % id of the shoulder case, as it appears in the database (Pnnn)
diagnosis
treatment
outcome
patient
shoulder
study
comment
end
properties (Constant, Hidden = true)
dataPath = '../../../data'; % path to data folder from package folder
% This should be done differently. Check where we are, and where we
% should go.
end
properties (Hidden = true)
dataCTPath
dataAmiraPath
dataMatlabPath
id4C % sCaseId with P/N followed by 3 digits --> 4 char
end
methods
function obj = ShoulderCase(sCaseId)
%SHOULDERCASE Construct an instance of this class
% If called with argument sCaseId, set id and path of CT
if nargin > 0
obj.id = sCaseId;
%{
% sCaseId must be P or N followed by 1-3 digit.
if ischar(sCaseId)
if strlength(sCaseId) == 4
obj.id = sCaseId;
else
error('id must have 4 char');
end
else
% error('id must be a string');
end
%}
else
% error('id is required');
end
% Initialsation
obj.diagnosis = '';
obj.treatment = '';
obj.outcome = '';
obj.patient = Patient(obj);
obj.shoulder = Shoulder(obj);
obj.study = '';
obj.comment = '';
obj.dataCTPath = '';
obj.dataAmiraPath = '';
obj.dataMatlabPath = '';
obj.id4C = '';
% Set path of the sCase
% obj.path; To be set only if varargin ~isempty
end
function outputArg = path(obj) % Should be private
% PATH Summary of this method goes here
% Detailed explanation goes here
sCase = obj;
sCaseId = sCase.id;
% The validity of the format should be checked Pnnn or Nnnn.
if (numel(regexp(sCaseId,'^[PN]\d{1,3}$')) == 0)
error(['Invalid format of sCaseId: ' sCaseId{1i} '. CaseID must start with "P" or "N" and be followed by 1 to 3 digits.']);
end
% Set the folder of the sCaseId
sCaseDirLevelPN = sCaseId(1); % Either 'P' or 'N'
strLengthsCaseId = strlength(sCaseId(2:end));
if (strLengthsCaseId < 2)
sCaseDirLevel2 = '0'; % Hunderets
sCaseDirLevel1 = '0'; % Dozent
elseif (strLengthsCaseId < 3)
sCaseDirLevel2 = '0'; % Hunderets
sCaseDirLevel1 = sCaseId(2); % Dozent
else
sCaseDirLevel2 = sCaseId(2); % Hunderets
sCaseDirLevel1 = sCaseId(3); % Dozent
end
% Set sCaseId with fixed 3 digits after P/N (needed when id <
% 10 inloading of data in amira directory.
sCaseId4C = [sCaseDirLevelPN sCaseDirLevel2 sCaseDirLevel1 sCaseId(end)];
obj.id4C = sCaseId4C;
% Check if a (!unique! to be done) directory exists for this sCaseId
dataDir = obj.dataPath;
FindsCaseIdFolder = dir([dataDir '/' sCaseDirLevelPN '/' sCaseDirLevel2 '/' sCaseDirLevel1 '/' sCaseId '*']);
if (isempty(FindsCaseIdFolder)) % No directory for this sCaseId
error(['Missing directory for sCaseId: ' sCaseId]);
end
sCaseDirLevel0 = FindsCaseIdFolder.name;
sCaseDir = [dataDir '/' sCaseDirLevelPN '/' sCaseDirLevel2 '/' sCaseDirLevel1 '/' sCaseDirLevel0];
% Check if this sCaseId has a CT directory with '-1' postfix (preoperative
% sharp CT
FindCTFolder = dir([sCaseDir '/' 'CT*-1']);
if (isempty(FindCTFolder)) % No CT directory for this sCaseId
error(['Missing directory for sCaseId: ' sCaseId]);
% ! should exist script here !
end
CTDir = [sCaseDir '/' FindCTFolder.name];
amiraDir = [CTDir '/amira']; % we should check for existance
matlabDir = [CTDir '/matlab']; % we should check for existance
obj.dataCTPath = CTDir;
obj.dataAmiraPath = amiraDir;
obj.dataMatlabPath = matlabDir;
outputArg = 1; % Should report on directories existence
if exist(amiraDir, 'dir') ~= 7
outputArg = 0;
end
end
function outputArg = output(obj, varargin)
% Below is copy/past from previous version
% This function is used to output the variable Results, used
% by scapula_measurePKG the modified funcion scapula_measure.
% inputArg could be as in scapula_calculation
% 'density', 'References', 'obliqueSlice', 'display'
Result.sCase_id = obj.id; % 1
Result.glenoid_Radius = obj.shoulder.scapula.glenoid.radius; % 2
Result.glenoid_radiusRMSE = obj.shoulder.scapula.glenoid.sphereRMSE; % 3 % 3
%Result.glenoidSphericity = ' ';
%Result.glenoid_biconcave = ' ';
Result.glenoid_depth = obj.shoulder.scapula.glenoid.depth; % 4
Result.glenoid_width = obj.shoulder.scapula.glenoid.width; % 5
Result.glenoid_height = obj.shoulder.scapula.glenoid.height; % 6
Result.glenoid_center_PA = obj.shoulder.scapula.glenoid.centerLocal(1); % 7
Result.glenoid_center_IS = obj.shoulder.scapula.glenoid.centerLocal(2); % 8
Result.glenoid_center_ML = obj.shoulder.scapula.glenoid.centerLocal(3); % 9
Result.glenoid_version_ampl = obj.shoulder.scapula.glenoid.versionAmpl; % 10
Result.glenoid_version_orient = obj.shoulder.scapula.glenoid.versionOrient; % 11
Result.glenoid_Version = obj.shoulder.scapula.glenoid.version; % 12
Result.glenoid_Inclination = obj.shoulder.scapula.glenoid.inclination; % 13
Result.humerus_joint_radius = ' '; % 14
Result.humeral_head_radius = obj.shoulder.humerus.radius; % 15
% Result.humerus_GHsublux_2D = ' ';
% Result.humerus_SHsublux_2D = ' ';
Result.humerus_GHsubluxation_ampl = obj.shoulder.humerus.GHSAmpl; % 16
Result.humerus_GHsubluxation_orient = obj.shoulder.humerus.GHSOrient; % 17
Result.humerus_SHsubluxation_ampl = obj.shoulder.humerus.SHSAmpl; % 18
Result.humerus_SHsubluxation_orient = obj.shoulder.humerus.SHSOrient; % 19
Result.scapula_CSA = obj.shoulder.scapula.acromion.criticalShoulderAngle; % radCSA;
% 5 Lines below should be updated
Result.scapula_CTangle = 0; %CTorientation; %20
Result.scapula_CTangleVersion = 0; % WearPlaneAngle; %21
Result.scapula_CTangleSHS = 0; % SHSPlaneAngle; % 22
Result.scapula_CTangleGHS = 0; % GHSPlaneAngle; % 23
Result.scapula_PlaneRMSE = 0; %PlaneRMSE;%24
Result.scapula_AI = obj.shoulder.scapula.acromion.acromionIndex;
outputArg = Result;
end
function outputArg = saveMatlab(obj)
% Save sCase to matlab file
dir = obj.dataMatlabPath;
% Create dir if not exist
try
if ~exist(dir, 'dir') % create directory if it does not exist
mkdir(dir);
end
catch
fprintf('Error creating the matlab directory \n'); % Should be in log
end
% Save sCase in matlab directoty, in a file named sCaseCNNN.m
filename = 'SCase';
filename = [dir '/' filename '.mat'];
try
sCase = obj;
save(filename, 'sCase');
outputArg = 1;
catch
fprintf('Error creating sCase matlab file \n'); % Should be in log
outputArg = -1;
end
end
function outputArg = saveCSV(~)
% Save sCase to csv file
outputArg = 1;
end
function outputArg = saveExcel(~)
% Save sCase to Excel file
outputArg = 1;
end
function outputArg = saveSQL(~)
% Save sCase to MySQL database
outputArg = 1;
end
function outputArg = plot(obj)
% Plot SCase
figure; % Create figure
hold on; % Superpose objects on the same figure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Scapula
% Plot the scapula surface as points
% If exists, plot the segemented surface (manual)
if strcmp(obj.shoulder.scapula.segmentation, 'M')
points = obj.shoulder.scapula.surface.points;
faces = obj.shoulder.scapula.surface.faces;
x = points(:,1);
y = points(:,2);
z = points(:,3);
trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',1,'EdgeColor','none','FaceLighting','gouraud');
%scatter3(x,y,z,1,'black'); % Point points
end
% If exists, plot the segemented surface (auto)
if strcmp(obj.shoulder.scapulaAuto.segmentation, 'A')
points = obj.shoulder.scapulaAuto.surface.points;
faces = obj.shoulder.scapulaAuto.surface.faces;
x = points(:,1);
y = points(:,2);
z = points(:,3);
trisurf(faces,x,y,z,'Facecolor','red','FaceAlpha',0.25,'EdgeColor','none','FaceLighting','gouraud');
end
% Plot scapula groove (manual)
points = obj.shoulder.scapula.groove;
x = points(:,1);
y = points(:,2);
z = points(:,3);
scatter3(x,y,z,100,'blue');
% If exists, plot segments between manual and auto
% Plot scapula markers
points = [...
obj.shoulder.scapula.angulusInferior;
obj.shoulder.scapula.trigonumSpinae;
obj.shoulder.scapula.processusCoracoideus;
obj.shoulder.scapula.acromioClavicular;
obj.shoulder.scapula.angulusAcromialis;
obj.shoulder.scapula.spinoGlenoidNotch...
];
x = points(:,1);
y = points(:,2);
z = points(:,3);
scatter3(x,y,z,100,'blue','LineWidth',2);
% Plot scapular plane
axisLength = 50;
pt0 = obj.shoulder.scapula.coordSys.origin;
pt1 = pt0 + (obj.shoulder.scapula.coordSys.ML + obj.shoulder.scapula.coordSys.IS) * axisLength ;
pt2 = pt1 - obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]);
pt3 = pt2 - obj.shoulder.scapula.coordSys.ML * pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae] );
pt4 = pt3 + obj.shoulder.scapula.coordSys.IS * pdist([pt1 ; obj.shoulder.scapula.angulusInferior]);
X = [pt1(1) pt2(1) pt3(1) pt4(1)];
Y = [pt1(2) pt2(2) pt3(2) pt4(2)];
Z = [pt1(3) pt2(3) pt3(3) pt4(3)];
patch(X,Y,Z,'black','FaceAlpha',0.3); % Transparent plane
% Plot scapular axis
pt1 = obj.shoulder.scapula.coordSys.origin;
axisLength = 10 + pdist([pt1 ; obj.shoulder.scapula.trigonumSpinae]);
pt2 = pt1 - obj.shoulder.scapula.coordSys.ML * axisLength;
x = [pt1(1), pt2(1)];
y = [pt1(2), pt2(2)];
z = [pt1(3), pt2(3)];
plot3(x,y,z,'Color','blue','LineWidth', 5);
% Plot scapular coordinate system
axisLength = 50; % Length of the cood. sys. axes
pt0 = obj.shoulder.scapula.coordSys.origin;
pt1 = pt0 + obj.shoulder.scapula.coordSys.PA*axisLength; % X
pt2 = pt0 + obj.shoulder.scapula.coordSys.IS*axisLength; % Y
pt3 = pt0 + obj.shoulder.scapula.coordSys.ML*axisLength; % Z
x = [pt0(1), pt1(1)];
y = [pt0(2), pt1(2)];
z = [pt0(3), pt1(3)];
plot3(x,y,z,'Color','red','LineWidth', 5);
x = [pt0(1), pt2(1)];
y = [pt0(2), pt2(2)];
z = [pt0(3), pt2(3)];
plot3(x,y,z,'Color','green','LineWidth', 5);
% Not plotting ML axis to avoid its superposition with scapular axis on glenoid center
% x = [pt0(1), pt3(1)];
% y = [pt0(2), pt3(2)];
% z = [pt0(3), pt3(3)];
% plot3(x,y,z,'Color','blue','LineWidth', 5);
% plot scapular axis on glenoid center
axisLength = 50;
pt1 = obj.shoulder.scapula.glenoid.center;
pt2 = pt1 + obj.shoulder.scapula.coordSys.ML*axisLength;
x = [pt1(1), pt2(1)];
y = [pt1(2), pt2(2)];
z = [pt1(3), pt2(3)];
plot3(x,y,z,'Color','blue','LineWidth', 5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Glenoid
% Plot the glenoid surface
% We might move them in the direction of th glenoid center
% instead of the ML axis
points = obj.shoulder.scapula.glenoid.surface.points +...
obj.shoulder.scapula.coordSys.ML * 0.2;
faces = obj.shoulder.scapula.glenoid.surface.faces;
x = points(:,1);
y = points(:,2);
z = points(:,3);
%scatter3(x,y,z,1,'blach');
trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','yellow','FaceLighting','none');
% plot glenoid centerline
pt1 = obj.shoulder.scapula.glenoid.center;
pt2 = pt1 + obj.shoulder.scapula.glenoid.centerLine;
x = [pt1(1), pt2(1)];
y = [pt1(2), pt2(2)];
z = [pt1(3), pt2(3)];
plot3(x,y,z,'Color','yellow','LineWidth', 5);
% We might add spherical cap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Humerus
% plot humeral head sphere
n = 20; % number of lines
[X,Y,Z] = sphere(n);
X = X*obj.shoulder.humerus.radius;
Y = Y*obj.shoulder.humerus.radius;
Z = Z*obj.shoulder.humerus.radius;
Xc = obj.shoulder.humerus.center(1);
Yc = obj.shoulder.humerus.center(2);
Zc = obj.shoulder.humerus.center(3);
mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
% hidden off; % Transparent sphere
% check 'FaceLighting','gouraud',
% plot humeral head centerline
pt1 = obj.shoulder.scapula.glenoid.center;
pt2 = obj.shoulder.humerus.center;
x = [pt1(1), pt2(1)];
y = [pt1(2), pt2(2)];
z = [pt1(3), pt2(3)];
plot3(x,y,z,'Color','magenta','LineWidth', 5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Graph properties
% axis off; % remove axis
lightPosition = obj.shoulder.scapula.glenoid.center + ...
(...
obj.shoulder.scapula.coordSys.ML + ...
obj.shoulder.scapula.coordSys.IS + ...
obj.shoulder.scapula.coordSys.PA ...
) * 100;
light('Position',lightPosition,'Style','local');
grid on;
xlabel('X (CT)');
ylabel('Y (CT)');
zlabel('Z (CT)');
axis square;
axis vis3d; % fixed limits and scaling
view(obj.shoulder.scapula.coordSys.IS); % view from top
% Align ML axis with window horizontal axis !to do
% obj.shoulder.scapula.coordSys.ML
%angleDegrees = rad2deg(obj.shoulder.scapula.coordSys.ML
% camroll(angleDegrees);
figureHandle = gca;
figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis
figureHandle.Box = 'On';
outputArg = 1;
end
end
end
c4science · Help