diff --git a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
index aedc83c..3ed8664 100644
--- a/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
+++ b/ShoulderCase/@GlenoidAuto/GlenoidAuto.m
@@ -1,38 +1,38 @@
classdef GlenoidAuto < Glenoid
% To be used with ShoulderAuto data.
% The load() method requires a specific implementation.
methods
function outputArg = load(obj)
% LOADAUTO Load genoid surface from auto segmentation
% Load the points of the glenoid surface from matlab dir
SCase = obj.scapula.shoulder.SCase;
SCaseId4C = SCase.id4C;
matlabDir = SCase.dataMatlabPath;
side = SCase.shoulderAuto.side;
% Import glenoid surface points
- fileName = ['glenoidSurfaceAuto' side '.ply'];
- fileName = [matlabDir '/' fileName];
+ fileName = "glenoidSurfaceAuto" + side + ".ply";
+ fileName = fullfile(matlabDir, fileName);
if exist(fileName,'file') == 2
[points,faces]=loadPly(fileName,1);
obj.surface.points = points;
obj.surface.meanPoint = mean(points);
obj.surface.faces = faces;
else
outputArg = 0;
return
% error(['Could not find file: ' fileName])
end
outputArg = 1; % Should report on loading result
end
end
end
diff --git a/ShoulderCase/@GlenoidManual/GlenoidManual.m b/ShoulderCase/@GlenoidManual/GlenoidManual.m
index cc14619..45b3e7e 100644
--- a/ShoulderCase/@GlenoidManual/GlenoidManual.m
+++ b/ShoulderCase/@GlenoidManual/GlenoidManual.m
@@ -1,39 +1,39 @@
classdef GlenoidManual < Glenoid
% To be used with ShoulderManual data.
% The load() method requires a specific implementation.
methods
function outputArg = load(obj)
% LOAD Summary of this method goes here
% Load the points of the glenoid surface from the amira
% directoty.
SCase = obj.scapula.shoulder.SCase;
SCaseId4C = SCase.id4C;
amiraDir = SCase.dataAmiraPath;
% Import glenoid surface points
- fileName = ['ExtractedSurface' SCaseId4C '.stl'];
- fileName = [amiraDir '/' fileName];
+ fileName = "ExtractedSurface" + SCaseId4C + ".stl";
+ fileName = fullfile(amiraDir, fileName);
if exist(fileName,'file') == 2
[points,faces,~]=loadStl(fileName,1);
obj.surface.points = points;
obj.surface.meanPoint = mean(points);
obj.surface.faces = faces;
else
- warning(['Could not find file: ' fileName]);
+ warning("Could not find file: " + fileName);
outputArg = 0;
return;
% error(['Could not find file: ' fileName])
end
outputArg = 1; % Should report on loading result
end
end
end
diff --git a/ShoulderCase/@Humerus/Humerus.m b/ShoulderCase/@Humerus/Humerus.m
index 906b9d8..337bd21 100644
--- a/ShoulderCase/@Humerus/Humerus.m
+++ b/ShoulderCase/@Humerus/Humerus.m
@@ -1,257 +1,225 @@
classdef Humerus < handle
%HUMERUS Properties and methods associated to the humerus
% Detailed explanation goes here
% Author: Alexandre Terrier, EPFL-LBO
% Creation date: 2018-07-01
% Revision date: 2019-06-29
%
% TO DO:
% Local coordinate system
properties
landmarks % 5 3D points
center % Center of the humeral head (sphere fit on 5 points
radius % Radius of the humeral head (sphere fit on 5 points
jointRadius % Radius of cuvature of the articular surface (todo)
SHSAngle % Scapulo-humeral subluxation angle
SHSPA % Scapulo-humeral subluxation angle in the postero-anterior direction (posterior is negative, as for glenoid version)
SHSIS % Scapulo-humeral subluxation angle in the infero-superior direction
SHSAmpl % Scapulo-humeral subluxation (center ofset / radius)
SHSOrient % Scapulo-humral subluxation orientation (0 degree is posterior)
GHSAmpl % Gleno-humeral subluxation (center ofset / radius)
GHSOrient % Gleno-humral subluxation orientation (0 degree is posterior)
subluxationIndex3D
shoulder
end
methods (Access = ?Shoulder) % Only Shoulder is allowed to construct Humerus
function obj = Humerus(shoulder)
%HUMERUS Construct an instance of this class
% Detailed explanation goes here
obj.landmarks = [];
obj.center = [];
obj.radius = [];
obj.jointRadius = [];
obj.SHSAngle = [];
obj.SHSPA = [];
obj.SHSIS = [];
obj.SHSAmpl = [];
obj.SHSOrient = [];
obj.GHSAmpl = [];
obj.GHSOrient = [];
obj.shoulder = shoulder;
end
end
methods
function outputArg = load(obj)
- %LOAD Load 5 humeral head landmarks
- % TODO: Should be check for 'old' vs 'new' file definition
-
SCase = obj.shoulder.SCase;
- SCaseId4C = SCase.id4C;
- amiraDir = SCase.dataAmiraPath;
-
- % Try 1st with new defiunition of HHLandmarks
- fileName = "newHHLandmarks" + SCaseId4C + ".landmarkAscii";
- fileName = fullfile(amiraDir, fileName);
-
- if exist(fileName,'file') == 2
- importedData = importdata(fileName, ' ', 14);
- else % Check for old version of HHLandmarks
- fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii'];
- fileName = [amiraDir '/' fileName];
- if exist(fileName,'file') == 2
- importedData = importdata(fileName, ' ', 14);
- warning([SCaseId4C ' is using old humeral landmarks definition'])
- else
- outputArg = 0;
- return
- end
+ if exist(fullfile(SCase.dataSlicerPath, "HH_landmarks.mrk.json"), "file")
+ success = obj.loadSlicerLandmarks();
+ elseif not(isempty(dir(fullfile(SCase.dataAmiraPath, "*HH*.landmarksAscii"))))
+ success = obj.loadAmiraLandmarks();
+ else
+ success = 0;
end
- % with other methods
-
- % Check validity of data in file
- % Variable points belos should be renamed landmarks, for consistency
- try
- points = importedData.data;
- catch
- % Cant get landarks from file
- warning([SCaseId4C ' cant load humeral landmarks'])
- outputArg = 0;
- return;
- end
- % We should add a test for validity of points
-
-
- obj.landmarks = points;
- outputArg = 1; % Should report on loading result
+ outputArg = success; % Should report on loading result
end
function outputArg = subluxation(obj)
% Calcul of humeral head subluxation
% The function also get center and radius of the humeral head.
% TODO:
% Calcul would be simpler in the scapular cordinate system, as
% done in Acromion object.
% Radius of curvature of humeral head articular surface
% Might be removed anatomy for concistency
SCase = obj.shoulder.SCase;
scapula = obj.shoulder.scapula;
%% Get data from parent shoulder object
% Shoulder side
side = obj.shoulder.side;
%% Get data from parent scapula object
% Scapula coordinate system
xAxis = scapula.coordSys.PA;
yAxis = scapula.coordSys.IS;
zAxis = scapula.coordSys.ML;
%% Get data from parent glenoid object
% Glenoid center
glenCenter = scapula.glenoid.center;
% Glenoid centerLine from parent glenoid object
glenCenterline = scapula.glenoid.centerLine;
%% Get data from humerus
% humeral head landmarks (5 or ?)
humHead = obj.landmarks; % Humeral head landmarks
%% Anatomical measurements
% Fit sphere on humeral head landmarks
[humHeadCenter,humHeadRadius,HHResiduals, R2HH] = fitSphere(humHead);
humHeadCenter = humHeadCenter'; % transform vertical to horizantal vector
%% Scapulo-humeral subluxation (SHS)
% Vector from glenoid center to humeral head center
SHSvect = humHeadCenter - glenCenter;
% Vector from humeral head center to scapular (z) axis,
% perpendicular to z-axis --> similar to projection in xy plane
% This might be rewritten more clearly
SHSvectxy = SHSvect - (dot(SHSvect,zAxis)*zAxis);
% SHS amplitude (ratio rather than amplitude)
SHSAmpl = norm(SHSvectxy) / humHeadRadius;
% SHS orientation
xProj = dot(SHSvectxy, xAxis);
yProj = dot(SHSvectxy, yAxis);
if xProj > 0 % anterior side
if yProj > 0 % ant-sup quadran
SHSOrient = pi - atan(yProj/xProj);
else % ant-inf quadran
SHSOrient = -pi - atan(yProj/xProj);
end
elseif xProj < 0 % posterior side
SHSOrient = - atan(yProj/xProj);
else % xProj = 0
SHSOrient = pi/2 - pi * (yProj < 0);
end
SHSOrient = rad2deg(SHSOrient);
% SHSAngle: angle between glenHumHead and zAxis
SHSAngle = vrrotvec(SHSvect,zAxis);
SHSAngle = SHSAngle(4);
SHSAngle = rad2deg(SHSAngle);
% SHSAP is angle between zxProjection and zAxis
% Project glenHumHead on yAxis
yProjVect = dot(SHSvectxy, yAxis)*yAxis;
% Remove projection from glenHumHead --> zxProjection
glenHumHeadzx = SHSvect - yProjVect;
a = glenHumHeadzx;
b = zAxis;
rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b
rotVect3 = rotVect4(1:3);
rotAngle = rotVect4(4);
rotSign = sign(dot(rotVect3,yAxis)); % Sign of rotation
SHSPA = rotAngle * rotSign;
SHSPA = -SHSPA; % Convention that posterior is negative
if strcmp(side,'L') % Change signe for left-handed scapula coordinate system
SHSPA = -SHSPA;
end
SHSPA = rad2deg(SHSPA); % Angle in degrees
% SHSIS is angle between yzProjection and zAxis
% Project glenHumHead on xAxis
xProjVect = dot(SHSvectxy, xAxis)*xAxis;
% Remove projection from glenHumHead --> yzProjection
glenHumHeadyz = SHSvect - xProjVect;
a = glenHumHeadyz;
b = zAxis;
rotVect4 = vrrotvec(a,b); % Unit rotation vector to align a to b
rotVect3 = rotVect4(1:3);
rotAngle = rotVect4(4);
rotSign = sign(dot(rotVect3,xAxis)); % Sign of rotation
SHSIS = rotAngle * rotSign;
SHSIS = SHSIS; % Convention that superior is positive
if strcmp(side,'L') % Change signe for left-handed scapula coordinate system
SHSIS = -SHSIS;
end
SHSIS = rad2deg(SHSIS); % Angle in degrees
%% Angle between plane of maximum SH subluxation and CT axis [0,0,1]
% Not used anymore
% SHSPlane = cross(glenHumHead, zAxis);
% G = vrrotvec(SHSPlane, [0 0 1]);
% SHSPlaneAngle = rad2deg(G(4));
% if SHSPlaneAngle > 90
% SHSPlaneAngle = 180 - SHSPlaneAngle;
% end
%% Gleno-humeral subluxation (GHS)
% Vector from glenoid sphere centre to glenoid centerline (perpendicular)
glenCenterlineNorm = glenCenterline/norm(glenCenterline);
GHS = SHSvect - dot(SHSvect, glenCenterlineNorm)*glenCenterlineNorm;
% GHS amplitude
GHSAmpl = norm(GHS) / humHeadRadius;
% GHS orientation
xProj = dot(GHS, -xAxis);
yProj = dot(GHS, yAxis);
if xProj ~= 0
if yProj ~= 0
GHSOrient = atan(yProj/xProj);
else
GHSOrient = pi * (xProj > 0);
end
else
GHSOrient = pi/2 - pi * (yProj > 0);
end
GHSOrient = rad2deg(GHSOrient);
%% Angle between plane of maximum GH subluxation and CT axis [0,0,1]
% Not used anymore
% GHSPlane = cross(glenCenterline, glenHumHead);
% GHSPlaneAngle = vrrotvec(GHSPlane, [0 0 1]);
% GHSPlaneAngle = rad2deg(GHSPlaneAngle(4));
% if GHSPlaneAngle > 90
% GHSPlaneAngle = 180 - GHSPlaneAngle;
% end
%% Set humerus properties
obj.center = humHeadCenter;
obj.radius = humHeadRadius;
obj.SHSAngle = SHSAngle;
obj.SHSPA = SHSPA;
obj.SHSIS = SHSIS;
obj.SHSAmpl = SHSAmpl;
obj.SHSOrient = SHSOrient;
obj.GHSAmpl = GHSAmpl;
obj.GHSOrient = GHSOrient;
% Commented ince not used
% outputArgStruct = struct('HHResiduals',HHResiduals,'R2HH',R2HH); % Useful?
outputArg = 1;
end
end
end
diff --git a/ShoulderCase/@Humerus/loadAmiraLandmarks.m b/ShoulderCase/@Humerus/loadAmiraLandmarks.m
new file mode 100644
index 0000000..302afb3
--- /dev/null
+++ b/ShoulderCase/@Humerus/loadAmiraLandmarks.m
@@ -0,0 +1,44 @@
+function outputArg = loadAmiraLandmarks(obj)
+ %LOAD Load 5 humeral head landmarks
+ % TODO: Should be check for 'old' vs 'new' file definition
+
+ SCase = obj.shoulder.SCase;
+ SCaseId4C = SCase.id4C;
+ amiraDir = SCase.dataAmiraPath;
+
+ % Try 1st with new defiunition of HHLandmarks
+ fileName = "newHHLandmarks" + SCaseId4C + ".landmarkAscii";
+ fileName = fullfile(amiraDir, fileName);
+
+ if exist(fileName,'file') == 2
+ importedData = importdata(fileName, ' ', 14);
+ else % Check for old version of HHLandmarks
+ fileName = ['HHLandmarks' SCaseId4C '.landmarkAscii'];
+ fileName = [amiraDir '/' fileName];
+ if exist(fileName,'file') == 2
+ importedData = importdata(fileName, ' ', 14);
+ warning([SCaseId4C ' is using old humeral landmarks definition'])
+ else
+ outputArg = 0;
+ return
+ end
+ end
+ % with other methods
+
+ % Check validity of data in file
+ % Variable points belos should be renamed landmarks, for consistency
+ try
+ points = importedData.data;
+ catch
+ % Cant get landarks from file
+ warning([SCaseId4C ' cant load humeral landmarks'])
+ outputArg = 0;
+ return;
+ end
+ % We should add a test for validity of points
+
+
+ obj.landmarks = points;
+
+ outputArg = 1; % Should report on loading result
+end
\ No newline at end of file
diff --git a/ShoulderCase/@Humerus/loadSlicerLandmarks.m b/ShoulderCase/@Humerus/loadSlicerLandmarks.m
index 33251ed..9a8e24d 100644
--- a/ShoulderCase/@Humerus/loadSlicerLandmarks.m
+++ b/ShoulderCase/@Humerus/loadSlicerLandmarks.m
@@ -1,7 +1,12 @@
-function loadSlicerLandmarks(obj)
- landmarksFilename = fullfile(...
- obj.shoulder.SCase.dataSlicerPath,...
- "HH_landmarks.mrk.json");
- fileData = jsondecode(fileread(landmarksFilename));
- obj.landmarks = [fileData.markups.controlPoints.position]';
+function output = loadSlicerLandmarks(obj)
+ try
+ landmarksFilename = fullfile(...
+ obj.shoulder.SCase.dataSlicerPath,...
+ "HH_landmarks.mrk.json");
+ fileData = jsondecode(fileread(landmarksFilename));
+ obj.landmarks = [fileData.markups.controlPoints.position]';
+ output = 1;
+ catch
+ output = 0;
+ end
end
\ No newline at end of file
diff --git a/ShoulderCase/@RotatorCuff/RotatorCuff.m b/ShoulderCase/@RotatorCuff/RotatorCuff.m
index 6163329..2e95829 100644
--- a/ShoulderCase/@RotatorCuff/RotatorCuff.m
+++ b/ShoulderCase/@RotatorCuff/RotatorCuff.m
@@ -1,105 +1,110 @@
classdef RotatorCuff < handle
properties
SC
SS
IS
TM
anteroPosteriorImbalance = nan;
shoulder;
end
methods
function obj = RotatorCuff(shoulder)
obj.shoulder = shoulder;
obj.SC = Muscle(obj, "SC");
obj.SS = Muscle(obj, "SS");
obj.IS = Muscle(obj, "IS");
obj.TM = Muscle(obj, "TM");
[~, ~] = mkdir(obj.dataPath);
[~, ~] = mkdir(obj.dataSlicesPath);
end
function output = dataPath(obj)
output = fullfile(obj.shoulder.dataPath, "muscles");
end
function output = dataSlicesPath(obj)
output = fullfile(obj.dataPath, "oblique_slices");
end
function output = summary(obj)
summary = [];
summary = [summary; obj.SC.summary];
summary = [summary; obj.SS.summary];
summary = [summary; obj.IS.summary];
summary = [summary; obj.TM.summary];
output = summary;
end
+ function sliceAndSegment(obj)
+ obj.createRotatorCuffSliceMatthieu();
+ obj.segmentMuscles("rotatorCuffMatthieu", "autoMatthieu");
+ end
+
function output = measure(obj, sliceName, segmentationSet)
try
obj.measureDegenerations(sliceName, segmentationSet);
obj.measureInsertions();
obj.measureAnteroPosteriorImbalance();
end
output = obj.summary;
end
function measureDegenerations(obj, sliceName, segmentationSet)
for muscleName = ["SC", "SS", "IS", "TM"]
obj.(muscleName).measureDegeneration(sliceName, segmentationSet);
end
end
function measureInsertions(obj)
humerus = obj.shoulder.humerus;
humeralHead = Sphere(humerus.center, humerus.radius);
% IS humeral insertion estimation
obj.IS.measureInsertionTangentToSphere(...
humeralHead,...
obj.shoulder.scapula.coordSys.IS,...
-obj.shoulder.scapula.coordSys.PA);
obj.IS.measureForceVector();
% TM humeral insertion estimation
obj.TM.measureInsertionTangentToSphere(...
humeralHead,...
obj.shoulder.scapula.coordSys.IS,...
-obj.shoulder.scapula.coordSys.PA);
obj.TM.measureForceVector();
% SC humeral insertion estimation
obj.SC.measureInsertionTangentToSphere(...
humeralHead,...
obj.shoulder.scapula.coordSys.IS,...
obj.shoulder.scapula.coordSys.PA);
obj.SC.measureForceVector();
end
function measureAnteroPosteriorImbalance(obj)
forceSC = obj.SC.forceVector;
forceIS = obj.IS.forceVector;
forceTM = obj.TM.forceVector;
scapulaCS = obj.shoulder.scapula.coordSys;
forceResultant = Vector(scapulaCS.origin, scapulaCS.origin) + ...
forceSC + forceIS + forceTM;
% Take infero-superior orthogonal component <=> project on scapular axial plane
inferoSuperior = Vector(scapulaCS.origin, scapulaCS.origin + scapulaCS.IS);
imbalanceVector = inferoSuperior.orthogonalComplementToPoint(forceResultant.target);
lateroMedial = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.ML);
anteroPosterior = Vector(scapulaCS.origin, scapulaCS.origin - scapulaCS.PA);
obj.anteroPosteriorImbalance = ...
sign(dot(imbalanceVector, anteroPosterior)) *...
rad2deg(lateroMedial.angle(imbalanceVector));
end
end
end
\ No newline at end of file
diff --git a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
index a9eba51..c02cfdc 100644
--- a/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
+++ b/ShoulderCase/@ScapulaAuto/ScapulaAuto.m
@@ -1,103 +1,102 @@
classdef ScapulaAuto < Scapula
% To be used with ShoulderAuto data.
% The load() method requires a specific implementation.
% Instanciate a GlenoidAuto object.
methods
function createGlenoid(obj)
obj.glenoid = GlenoidAuto(obj);
end
function outputArg = load(obj)
% Load 2 files obtained by SSM auto segmentation
% scapulaSurfaceAuto{L/R}.ply
% scapulaLandmarksAuto{L/R}.mat
outputArg = 1; % To be set to 1 of loading is OK and 0 otherwise
SCase = obj.shoulder.SCase;
SCaseId4C = SCase.id4C;
matlabDir = SCase.dataMatlabPath;
% Set side equal to manual side or 'R' if manual side is not available
- manualSide = obj.shoulder.SCase.shoulderManual.side;
- if isequal(manualSide,'L')
- side = manualSide;
- otherSide = 'R';
- elseif isequal(manualSide,'R')
- side = manualSide;
- otherSide = 'L';
+ if not(isempty(obj.shoulder.side))
+ side = obj.shoulder.side;
+ elseif not(isempty(SCase.shoulderManual.side))
+ side = SCase.shoulderManual.side;
else
- side = 'R';
- otherSide = 'L';
+ side = "R";
end
+ assert(ismember(side, ["R" "L"]), "Shoulder side must be R or L");
+ otherSide = char((side == "R")*'L' + (side == "L")*'R');
+
% Try loading auto segmentation from matlab dir
- fileName = ['scapulaSurfaceAuto' side '.ply'];
+ fileName = "scapulaSurfaceAuto" + side + ".ply";
fileName = fullfile(matlabDir, fileName);
if ~exist(fileName,'file')
side = otherSide;
- fileName = ['scapulaSurfaceAuto' side '.ply'];
+ fileName = "scapulaSurfaceAuto" + side + ".ply";
fileName = fullfile(matlabDir, fileName);
end
%side = SCase.shoulder.side;
if exist(fileName,'file') == 2
try
face_index_start = 1;
[points,faces] = loadPly(fileName,face_index_start); % % Load ply file
% ptCloud = pcread(fileName); % load the pointCloud object in the ply file
% points = ptCloud.Location; % get the array of (x,y,z) points
obj.surface.points = points;
obj.surface.faces = faces;
obj.segmentation = 'A'; % Auto segmentation from matlab
catch
obj.segmentation = 'E'; % Error on loading
outputArg = 0;
end
else
obj.segmentation = 'N'; % No segmentation file
outputArg = 0;
end
% Try loading auto scapula landmarks from matlab dir
- filename = ['scapulaLandmarksAuto' side '.mat'];
- filename = [matlabDir '/' filename];
+ filename = "scapulaLandmarksAuto" + side + ".mat";
+ filename = matlabDir + "/" + filename;
if exist(filename,'file') == 2
try
% Load mat file with scapula landmarks
load(filename, 'ScapulaLandmarks');
% Set loaded landmarks to scapulaAuto (obj)
obj.angulusInferior = ScapulaLandmarks.angulusInferior;
obj.trigonumSpinae = ScapulaLandmarks.trigonumSpinae;
obj.processusCoracoideus = ScapulaLandmarks.processusCoracoideus;
obj.acromioClavicular = ScapulaLandmarks.acromioClavicular;
obj.angulusAcromialis = ScapulaLandmarks.angulusAcromialis;
obj.spinoGlenoidNotch = ScapulaLandmarks.spinoGlenoidNotch;
obj.pillar = ScapulaLandmarks.pillar;
obj.groove = ScapulaLandmarks.groove;
catch
disp(SCaseId4C); % For debug (2 cases)
outputArg = 0;
end
else
outputArg = 0;
end
try
obj.measurePlane;
catch
outputArg = 0;
end
end
end
end
diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m
index c1a7025..31ae096 100644
--- a/ShoulderCase/@Shoulder/Shoulder.m
+++ b/ShoulderCase/@Shoulder/Shoulder.m
@@ -1,24 +1,46 @@
classdef (Abstract) Shoulder < handle
% Contains all the shoulder parts (bones) and their measurements.
%
% Can be used to plot an overview of the case.
properties
side = ""; % R or L
scapula
humerus
rotatorCuff
CTscan = "";
comment = "";
SCase
end
methods
function obj = Shoulder(SCase)
obj.SCase = SCase;
obj.humerus = Humerus(obj);
obj.rotatorCuff = RotatorCuff(obj);
end
+
+ function loadData(obj)
+ obj.scapula.load();
+ obj.scapula.glenoid.load();
+ obj.humerus.load();
+ end
+
+ function measureMorphology(obj)
+ obj.scapula.coordSysSet();
+ obj.scapula.glenoid.morphology();
+ obj.humerus.subluxation();
+ obj.scapula.acromion.morphology();
+ end
+
+ function measureDensity(obj)
+ obj.scapula.genoid.calcDensity();
+ end
+
+ function measureMuscles(obj)
+ obj.scapula.rotatorCuff.sliceAndSegment();
+ obj.scapula.rotatorCuff.measure("rotatorCuffMatthieu", "autoMatthieu");
+ end
end
end