diff --git a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m index 10a6356..f9dd441 100644 --- a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m +++ b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m @@ -1,180 +1,181 @@ classdef CoordinateSystem < handle % Useful to define, use, and test a coordinate system. % % Can be tested on orthogonality, right or left handedness. % % Can be used to express global points in local coordinates % and vice versa. % Can be used to project points on planes spanned by the % CoordinateSystem axes. properties origin = []; end properties (Hidden = true) xAxis = []; yAxis = []; zAxis = []; rotationMatrix = []; end methods function obj = CoordinateSystem() end function output = isEmpty(obj) output = isempty(obj.origin) ||... isempty(obj.xAxis) ||... isempty(obj.yAxis) ||... isempty(obj.zAxis); end function set.xAxis(obj,value) obj.xAxis = value; obj.rotationMatrix(:,1) = value'/norm(value); end function set.yAxis(obj,value) obj.yAxis = value; obj.rotationMatrix(:,2) = value'/norm(value); end function set.zAxis(obj,value) obj.zAxis = value; obj.rotationMatrix(:,3) = value'/norm(value); end function output = getRotationMatrix(obj) assert(obj.isOrthogonal,"CoordinateSystem is not orthogonal. Can't get its rotation matrix"); output = obj.rotationMatrix; end function localPoints = express(obj,globalPoints) % Give the coordinates of global cartesian points in present coordinate system assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); assert(size(globalPoints,2) == 3,'Points have to be a Nx3 double array.') localPoints = (globalPoints-obj.origin) * obj.rotationMatrix; end function globalPoints = getGlobalCoordinatesOfLocalPoints(obj,localPoints) % Give the coordinates of local points in global cartesian coordinate system assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); assert(size(localPoints,2) == 3,'Points have to be a Nx3 double array.') globalPoints = (localPoints * obj.rotationMatrix') + obj.origin; end function projected = projectOnXYPlane(obj,points) assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.xAxis);(obj.origin + obj.yAxis)]); projected = projectionPlane.projectOnPlane(points); end function projected = projectOnYZPlane(obj,points) assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.yAxis);(obj.origin + obj.zAxis)]); projected = projectionPlane.projectOnPlane(points); end function projected = projectOnZXPlane(obj,points) assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); projectionPlane = Plane; projectionPlane.fit([obj.origin;(obj.origin + obj.zAxis);(obj.origin + obj.xAxis)]); projected = projectionPlane.projectOnPlane(points); end function output = isOrthogonal(obj) % The axes comparison can't be done looking for strict equality due to % rounded values. Therefore the values must be evaluated with a given % tolerance assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of % a millionth of a radian if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle) output = false; return end if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle) output = false; return end if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle) output = false; return end output = true; end function output = isRightHanded(obj) % The comparison between a theoretical right-handed axis and the actual % z Axis can't be done looking for strict vector equality due to rounded % values. Therefore the norm of the sum of the two vectors is evaluated. assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); if not(obj.isOrthogonal) output = false; return end rightHandedAxis = cross(obj.xAxis,obj.yAxis); if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis) output = false; return end output = true; end function output = isLeftHanded(obj) assert(not(obj.isEmpty),'Coordinate System has not been initialized yet.'); if not(obj.isOrthogonal) output = false; return end if obj.isRightHanded output = false; return end output = true; end function output = plot(obj,axesColor,center) axisLength = 50; if (center) origin = [0 0 0]; else origin = obj.origin; end X = [origin; origin + obj.xAxis*axisLength]; Y = [origin; origin + obj.yAxis*axisLength]; Z = [origin; origin + obj.zAxis*axisLength]; plotHandle(1) = plot3(X(:,1),X(:,2),X(:,3),'Color',axesColor,'LineWidth', 3); + hold on % set the point a bit further to plot axis name X(2,:) = X(2,:) + obj.xAxis*2; plotHandle(2) = text(X(2,1),X(2,2),X(2,3),'X'); plotHandle(3) = plot3(Y(:,1),Y(:,2),Y(:,3),'Color',axesColor,'LineWidth', 3); % set the point a bit further to plot axis name Y(2,:) = Y(2,:) + obj.yAxis*2; plotHandle(4) = text(Y(2,1),Y(2,2),Y(2,3),'Y'); plotHandle(5) = plot3(Z(:,1),Z(:,2),Z(:,3),'Color',axesColor,'LineWidth', 3); % set the point a bit further to plot axis name Z(2,:) = Z(2,:) + obj.zAxis*2; plotHandle(6) = text(Z(2,1),Z(2,2),Z(2,3),'Z'); output = plotHandle; end end end diff --git a/ShoulderCase/@Glenoid/plot.m b/ShoulderCase/@Glenoid/plot.m index b9b6105..5c17b29 100644 --- a/ShoulderCase/@Glenoid/plot.m +++ b/ShoulderCase/@Glenoid/plot.m @@ -1,27 +1,29 @@ function output = plot(obj) % plot glenoid surface points = obj.surface.points; faces = obj.surface.faces; x = points(:,1); y = points(:,2); z = points(:,3); plotHandle(1) = trisurf(faces,x,y,z,... 'Facecolor','none',... 'EdgeColor','magenta',... 'EdgeAlpha',0.1,... 'FaceLighting','none'); + hold on + % plot glenoid centerline pt1 = obj.center; pt2 = pt1 + obj.centerLine; x = [pt1(1), pt2(1)]; y = [pt1(2), pt2(2)]; z = [pt1(3), pt2(3)]; plotHandle(2) = plot3(x,y,z,... 'Color','magenta',... 'LineWidth', 5); output = plotHandle; end \ No newline at end of file diff --git a/ShoulderCase/@Scapula/plotLandmarks.m b/ShoulderCase/@Scapula/plotLandmarks.m new file mode 100644 index 0000000..6a2bf9e --- /dev/null +++ b/ShoulderCase/@Scapula/plotLandmarks.m @@ -0,0 +1,54 @@ +function output = plotLandmarks(obj,color,faceColor) + % get landmarks + AI = obj.angulusInferior; + TS = obj.trigonumSpinae; + PC = obj.processusCoracoideus; + AC = obj.acromioClavicular; + AA = obj.angulusAcromialis; + SG = obj.spinoGlenoidNotch; + if not(isempty(obj.groove)) + % sort groove from most lateral to most medial + rawGroove = obj.groove; + groove = []; + for i = 1:size(rawGroove,1) + nextGrooveIndex = findLongest3DVector(TS-rawGroove); + groove(i,:) = rawGroove(nextGrooveIndex(1),:); + rawGroove(nextGrooveIndex,:) = []; + end + end + + % plot wireframe + acromion = [SG; AA; AC]; + clavicle = [SG; PC]; + scapulaPlane = [groove(1,:); TS; AI; groove(1,:); SG]; + + plotHandle = []; + + plotHandle(end+1) = plot3(acromion(:,1),acromion(:,2),acromion(:,3),... + '-o',... + 'Color',color,... + 'MarkerSize',10,... + 'MarkerFaceColor',faceColor); + + hold on + + plotHandle(end+1) = plot3(clavicle(:,1),clavicle(:,2),clavicle(:,3),... + '-o',... + 'Color',color,... + 'MarkerSize',10,... + 'MarkerFaceColor',faceColor); + + plotHandle(end+1) = plot3(scapulaPlane(:,1),scapulaPlane(:,2),scapulaPlane(:,3),... + '-o',... + 'Color',color,... + 'MarkerSize',10,... + 'MarkerFaceColor',faceColor); + + plotHandle(end+1) = plot3(groove(:,1),groove(:,2),groove(:,3),... + '-o',... + 'Color',color,... + 'MarkerSize',10,... + 'MarkerFaceColor',faceColor); + + output = plotHandle; +end \ No newline at end of file diff --git a/ShoulderCase/@Scapula/plotLandmarksWireframe.m b/ShoulderCase/@Scapula/plotLandmarksWireframe.m deleted file mode 100644 index 2db7225..0000000 --- a/ShoulderCase/@Scapula/plotLandmarksWireframe.m +++ /dev/null @@ -1,27 +0,0 @@ -function output = plotLandmarksWireframe(obj,color,faceColor) - % get landmarks - AI = obj.angulusInferior; - TS = obj.trigonumSpinae; - PC = obj.processusCoracoideus; - AC = obj.acromioClavicular; - AA = obj.angulusAcromialis; - SG = obj.spinoGlenoidNotch; - if not(isempty(obj.groove)) - groove = obj.groove; - end - - % plot wireframe - acromion = [SG; AA; AC; AA; SG]; - clavicle = [SG; PC; SG]; - planarScapula = [SG; groove; TS; AI; SG]; - - wireframe = [planarScapula; acromion; clavicle]; - - plotHandle = plot3(wireframe(:,1),wireframe(:,2),wireframe(:,3),... - '-o',... - 'Color',color,... - 'MarkerSize',10,... - 'MarkerFaceColor',faceColor); - - output = plotHandle; -end \ No newline at end of file diff --git a/ShoulderCase/@Scapula/plotSurface.m b/ShoulderCase/@Scapula/plotSurface.m index 1538d52..8ee89f8 100644 --- a/ShoulderCase/@Scapula/plotSurface.m +++ b/ShoulderCase/@Scapula/plotSurface.m @@ -1,19 +1,21 @@ function output = plotSurface(obj,color,edgeColor) output = []; if not(isempty(obj.surface)) points = obj.surface.points; faces = obj.surface.faces; x = points(:,1); y = points(:,2); z = points(:,3); plotHandle = trisurf(faces,x,y,z,... 'Facecolor',color,... 'FaceAlpha',0.1,... 'EdgeColor',edgeColor,... 'EdgeAlpha',0.1,... 'FaceLighting','gouraud'); + hold on + output = plotHandle; end end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m b/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m index d604465..b646aec 100644 --- a/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m +++ b/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m @@ -1,68 +1,70 @@ function output = plotManualAutoDifference(obj,normLimitForOutliers) % Plot the vectors between auto and manual points. % % input: normLimitForOutliers % % The vector is green if below the value given to normLimitForOutliers, % otherwise it is red. % initialisation manual = obj.shoulderManual; manualPoints = []; auto = obj.shoulderAuto; autoPoints = []; % find maximum number of common groove points minGroovePoints = min(size(manual.scapula.groove),size(auto.scapula.groove)); if (minGroovePoints > 0) manualPoints = [manual.scapula.groove(1:minGroovePoints,:)]; autoPoints = [auto.scapula.groove(1:minGroovePoints,:)]; end % create points arrays manualPoints = [manualPoints;... manual.scapula.angulusInferior;... manual.scapula.trigonumSpinae;... manual.scapula.processusCoracoideus;... manual.scapula.acromioClavicular;... manual.scapula.angulusAcromialis;... manual.scapula.spinoGlenoidNotch;... manual.scapula.coordSys.origin;... manual.scapula.glenoid.center;... ]; autoPoints = [autoPoints;... auto.scapula.angulusInferior;... auto.scapula.trigonumSpinae;... auto.scapula.processusCoracoideus;... auto.scapula.acromioClavicular;... auto.scapula.angulusAcromialis;... auto.scapula.spinoGlenoidNotch;... auto.scapula.coordSys.origin;... auto.scapula.glenoid.center;... ]; % measure norm difference manualAutoDifference = vecnorm(autoPoints-manualPoints,2,2); % plot the difference vectors for i = 1:size(manualAutoDifference,1) if (manualAutoDifference(i) > normLimitForOutliers) % outlier color = 'red'; else % valid point color = 'green'; end differenceLine = [manualPoints(i,:);... autoPoints(i,:)]; plotHandle(i) = plot3(differenceLine(:,1),differenceLine(:,2),differenceLine(:,3),... ':x',... 'Color',color,... 'LineWidth',3); + + hold on end output = plotHandle; end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCasePlotter/ShoulderCasePlotter.m b/ShoulderCase/@ShoulderCasePlotter/ShoulderCasePlotter.m index 44fb4c6..be9f0b5 100644 --- a/ShoulderCase/@ShoulderCasePlotter/ShoulderCasePlotter.m +++ b/ShoulderCase/@ShoulderCasePlotter/ShoulderCasePlotter.m @@ -1,53 +1,53 @@ classdef ShoulderCasePlotter < handle % Create a figure with ShoulderCase's plots % and buttons to toggle the data visualisation. % % Used by the ShoulderCase.plot() method. properties (Access = private) SCase fig axesHandle plotHandle buttonHandle - options = {'wireframe',... + options = {'landmarks',... 'scapula surface',... 'glenoid',... 'coordinate system',... 'centered coordinate system',... 'rotator cuff',... 'difference'}; end methods (Access = ?ShoulderCase) function obj = ShoulderCasePlotter(SCase) obj.SCase = SCase; obj.fig = figure('Name',SCase.id,... 'NumberTitle','off',... 'units','normalized',... 'outerposition',[0 0 1 1]); obj.axesHandle = containers.Map; obj.plotHandle = {}; obj.buttonHandle = {}; obj.initializeLayout(); obj.plot(); obj.initializeButton(); end end methods (Access = private) initializeButton(obj); initializeLayout(obj); plot(obj); end end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCasePlotter/initializeButton.m b/ShoulderCase/@ShoulderCasePlotter/initializeButton.m index 2963299..f077d1c 100644 --- a/ShoulderCase/@ShoulderCasePlotter/initializeButton.m +++ b/ShoulderCase/@ShoulderCasePlotter/initializeButton.m @@ -1,48 +1,48 @@ function initializeButton(obj) autoPanel = uipanel(obj.fig,... - 'Title','Auto',... + 'Title','Auto (in blue)',... 'units','pixels',... 'Position',[10 320 180 150 ]); for i = 1:length(obj.options)-1 try obj.buttonHandle.auto(i) = uicontrol(autoPanel,... 'Style','checkbox',... 'Value',1,... 'String',obj.options{i},... 'position',[10,130-i*20,150,20],... 'Callback',{@setDataVisibility,obj.plotHandle.auto(obj.options{i})}); end end manualPanel = uipanel(obj.fig,... - 'Title','Manual',... + 'Title','Manual (in red)',... 'units','pixels',... 'Position',[10 160 180 150 ]); for i = 1:length(obj.options)-1 try obj.buttonHandle.manual(i) = uicontrol(manualPanel,... 'Style','checkbox',... 'Value',1,... 'String',obj.options{i},... 'position',[10,130-i*20,150,20],... 'Callback',{@setDataVisibility,obj.plotHandle.manual(obj.options{i})}); end end try obj.buttonHandle.difference = uicontrol(obj.fig,... 'Style','checkbox',... 'Value',1,... 'String','auto/manual difference',... 'position',[20,130,150,20],... 'Callback',{@setDataVisibility,obj.plotHandle.difference('difference')}); end end function setDataVisibility(button,event,dataHandle) for i = 1:length(dataHandle) set(dataHandle(i),'Visible',logical(button.Value)); end end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCasePlotter/initializeLayout.m b/ShoulderCase/@ShoulderCasePlotter/initializeLayout.m index a95c7f5..911f4c2 100644 --- a/ShoulderCase/@ShoulderCasePlotter/initializeLayout.m +++ b/ShoulderCase/@ShoulderCasePlotter/initializeLayout.m @@ -1,6 +1,8 @@ function initializeLayout(obj) obj.axesHandle('scapula') = subplot(2,4,[1 2 5 6]); + set(obj.axesHandle('scapula'),'Visible',false); obj.axesHandle('centered coordinate system') = subplot(2,4,[3 4]); + set(obj.axesHandle('centered coordinate system'),'Visible',false); obj.axesHandle('manual muscles') = subplot(2,4,7); obj.axesHandle('auto muscles') = subplot(2,4,8); end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCasePlotter/plot.m b/ShoulderCase/@ShoulderCasePlotter/plot.m index c3d0501..d85dd52 100644 --- a/ShoulderCase/@ShoulderCasePlotter/plot.m +++ b/ShoulderCase/@ShoulderCasePlotter/plot.m @@ -1,91 +1,94 @@ function plot(obj) % set method colors autoColors = {'b','#D9FFFF'}; - manualColors = {'#D95319','#FFFF00'}; + manualColors = {'#D95319','#D9FFFF'}; % plot scapula data axes(obj.axesHandle('scapula')); title('Auto is blue, manual is red'); hold on; grid on axis vis3d rotate3d on auto = containers.Map; try - auto('wireframe') = obj.SCase.shoulderAuto.scapula.plotLandmarksWireframe(autoColors{1},autoColors{2}); + auto('landmarks') = obj.SCase.shoulderAuto.scapula.plotLandmarks(autoColors{1},autoColors{1}); end try - auto('scapula surface') = obj.SCase.shoulderAuto.scapula.plotSurface(autoColors{1},autoColors{2}); + auto('scapula surface') = obj.SCase.shoulderAuto.scapula.plotSurface(autoColors{2},autoColors{1}); end try auto('glenoid') = obj.SCase.shoulderAuto.scapula.glenoid.plot(); end try auto('coordinate system') = obj.SCase.shoulderAuto.scapula.coordSys.plot(autoColors{1},false); end manual = containers.Map; try - manual('wireframe') = obj.SCase.shoulderManual.scapula.plotLandmarksWireframe(manualColors{1},manualColors{2}); + manual('landmarks') = obj.SCase.shoulderManual.scapula.plotLandmarks(manualColors{1},manualColors{1}); end try - manual('scapula surface') = obj.SCase.shoulderManual.scapula.plotSurface(manualColors{1},manualColors{2}); + manual('scapula surface') = obj.SCase.shoulderManual.scapula.plotSurface(manualColors{2},manualColors{1}); end try manual('glenoid') = obj.SCase.shoulderManual.scapula.glenoid.plot(); end try manual('coordinate system') = obj.SCase.shoulderManual.scapula.coordSys.plot(manualColors{1},false); end outliersNormLimit = 10; difference = containers.Map; try difference('difference') = obj.SCase.plotManualAutoDifferences(outliersNormLimit); end - + view(obj.SCase.shoulderAuto.scapula.coordSys.ML); + + % plot coordinateSystem data axes(obj.axesHandle('centered coordinate system')); title('Auto is blue, manual is red'); hold on; grid on axis vis3d rotate3d on try auto('centered coordinate system') = obj.SCase.shoulderAuto.scapula.coordSys.plot(autoColors{1},true); end try manual('centered coordinate system') = obj.SCase.shoulderManual.scapula.coordSys.plot(manualColors{1},true); end - + view(obj.SCase.shoulderAuto.scapula.coordSys.ML); + % plot muscles segmentation axes(obj.axesHandle('auto muscles')); - title('Auto'); + title('Auto segmentation with auto landmarks'); hold on; try auto('rotator cuff') = obj.SCase.shoulderAuto.muscles.plot(); end axes(obj.axesHandle('manual muscles')); - title('Manual'); + title('Auto segmentation with manual landmarks'); hold on; try manual('rotator cuff') = obj.SCase.shoulderManual.muscles.plot(); end % store the data handles obj.plotHandle.auto = auto; obj.plotHandle.manual = manual; obj.plotHandle.difference = difference; end diff --git a/ShoulderCase/findLongest3DVector.m b/ShoulderCase/findLongest3DVector.m index 7f2f01d..dc0366c 100644 --- a/ShoulderCase/findLongest3DVector.m +++ b/ShoulderCase/findLongest3DVector.m @@ -1,11 +1,11 @@ function index = findLongest3DVector(vectors) - dimension = find(size(vectors) == 3); - if (dimension ~= 2) + if (size(vectors,2) ~= 3) error('You should provide an Nx3 vector array'); index = 0; return end type = 2; + dimension = 2; norms = vecnorm(vectors,type,dimension); index = find(norms == max(norms)); end