diff --git a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m index 3b3c273..10a6356 100644 --- a/ShoulderCase/@CoordinateSystem/CoordinateSystem.m +++ b/ShoulderCase/@CoordinateSystem/CoordinateSystem.m @@ -1,177 +1,180 @@ 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 plot(obj,axesColor,center) + 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]; - plot3(X(:,1),X(:,2),X(:,3),'Color',axesColor,'LineWidth', 3); + plotHandle(1) = plot3(X(:,1),X(:,2),X(:,3),'Color',axesColor,'LineWidth', 3); % set the point a bit further to plot axis name X(2,:) = X(2,:) + obj.xAxis*2; - text(X(2,1),X(2,2),X(2,3),'X') + plotHandle(2) = text(X(2,1),X(2,2),X(2,3),'X'); - plot3(Y(:,1),Y(:,2),Y(:,3),'Color',axesColor,'LineWidth', 3); + 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; - text(Y(2,1),Y(2,2),Y(2,3),'Y') + plotHandle(4) = text(Y(2,1),Y(2,2),Y(2,3),'Y'); - plot3(Z(:,1),Z(:,2),Z(:,3),'Color',axesColor,'LineWidth', 3); + 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; - text(Z(2,1),Z(2,2),Z(2,3),'Z') + 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 f97aac4..b9b6105 100644 --- a/ShoulderCase/@Glenoid/plot.m +++ b/ShoulderCase/@Glenoid/plot.m @@ -1,25 +1,27 @@ -function plot(obj) - % plot glenoid surface - points = obj.surface.points; - faces = obj.surface.faces; - x = points(:,1); - y = points(:,2); - z = points(:,3); - - trisurf(faces,x,y,z,... - 'Facecolor','none',... - 'EdgeColor','magenta',... - 'EdgeAlpha',0.1,... - 'FaceLighting','none'); +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'); - % 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)]; + % 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)]; - plot3(x,y,z,... - 'Color','magenta',... - 'LineWidth', 5); + plotHandle(2) = plot3(x,y,z,... + 'Color','magenta',... + 'LineWidth', 5); + + output = plotHandle; end \ No newline at end of file diff --git a/ShoulderCase/@MusclesContainer/plot.m b/ShoulderCase/@MusclesContainer/plot.m index cac30df..9872881 100644 --- a/ShoulderCase/@MusclesContainer/plot.m +++ b/ShoulderCase/@MusclesContainer/plot.m @@ -1,29 +1,30 @@ function output = plot(obj) % For now this function is specific to plotting the rotator cuff segmentation % results SCase = obj.shoulder.SCase; muscleNames = fields(obj.list); rotatorCuffCrossSection = imread(fullfile(obj.slicesDataPath,... [obj.list.(muscleNames{1}).sliceName '_ForSegmentation.png'])); leg = {}; - output = imshow(rotatorCuffCrossSection); + plotHandle(1) = imshow(rotatorCuffCrossSection); hold on colors = {'g','r','b','y'}; for i = 1:length(muscleNames) try segmentedImage = imread(fullfile(obj.list.(muscleNames{i}).maskDataPath,'auto_Mask.png')); catch continue end if (max(segmentedImage,[],'all') > 0) % if segmentation result is not empty leg{end+1} = muscleNames{i}; - contour(segmentedImage,'color',colors{i}); + [~,plotHandle(i+1)] = contour(segmentedImage,'color',colors{i}); end end - legend(leg); + plotHandle(i+2) = legend(leg); + output = plotHandle; end diff --git a/ShoulderCase/@Scapula/plotLandmarksWireframe.m b/ShoulderCase/@Scapula/plotLandmarksWireframe.m index 07622d9..2db7225 100644 --- a/ShoulderCase/@Scapula/plotLandmarksWireframe.m +++ b/ShoulderCase/@Scapula/plotLandmarksWireframe.m @@ -1,35 +1,27 @@ -function plotLandmarksWireframe(obj,color,faceColor) +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]; - clavicle = [SG; PC]; + acromion = [SG; AA; AC; AA; SG]; + clavicle = [SG; PC; SG]; planarScapula = [SG; groove; TS; AI; SG]; - plot3(acromion(:,1),acromion(:,2),acromion(:,3),... - '-o',... - 'Color',color,... - 'MarkerSize',10,... - 'MarkerFaceColor',faceColor); + wireframe = [planarScapula; acromion; clavicle]; - plot3(clavicle(:,1),clavicle(:,2),clavicle(:,3),... + plotHandle = plot3(wireframe(:,1),wireframe(:,2),wireframe(:,3),... '-o',... 'Color',color,... 'MarkerSize',10,... 'MarkerFaceColor',faceColor); - plot3(planarScapula(:,1),planarScapula(:,2),planarScapula(:,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 c23486e..1538d52 100644 --- a/ShoulderCase/@Scapula/plotSurface.m +++ b/ShoulderCase/@Scapula/plotSurface.m @@ -1,16 +1,19 @@ -function plotSurface(obj,color,edgeColor) +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); - trisurf(faces,x,y,z,... + plotHandle = trisurf(faces,x,y,z,... 'Facecolor',color,... 'FaceAlpha',0.1,... 'EdgeColor',edgeColor,... 'EdgeAlpha',0.1,... 'FaceLighting','gouraud'); + + output = plotHandle; end end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m index e427195..0b82305 100755 --- a/ShoulderCase/@ShoulderCase/ShoulderCase.m +++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m @@ -1,217 +1,212 @@ classdef ShoulderCase < handle % Properties and methods associated to the SoulderCase object. % Author: Alexandre Terrier, EPFL-LBO % Matthieu Boubat, EPFL-LBO % Creation date: 2018-07-01 - % Revision date: 2020-06-09 + % Revision date: 2020-07-23 % 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 = []; shoulderManual = []; shoulderAuto = []; study = []; comment = []; end % properties (Constant, Hidden = true) % dataPath = '/Volumes/shoulder/dataDev'; % 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) dataPath = []; dataCTPath = []; dataAmiraPath = []; dataMatlabPath = []; dataDicomPath = []; id4C = []; % SCaseId with P/N followed by 3 digits --> 4 char - - plotFigure = []; - plotOptions = []; - plotLayout = []; - plotHandle = []; end methods (Access = ?ShoulderCaseLoader) function obj = ShoulderCase(SCaseID,dataCTPath) % SCaseID validation rawSCaseID = SCaseIDParser(SCaseID); assert(rawSCaseID.isValidID,'The input argument is not a valid SCaseID.') obj.id = SCaseID; obj.id4C = rawSCaseID.getIDWithNumberOfDigits(4); % path attribution obj.dataCTPath = dataCTPath; % Initialsation obj.patient = Patient(obj); obj.shoulderManual = ShoulderManual(obj); obj.shoulderAuto = ShoulderAuto(obj); obj.propagateDataPath; end end methods function propagateDataPath(obj) % Update current dataPath % Propagate the path to objects in properties obj.dataAmiraPath = fullfile(obj.dataCTPath,'amira'); obj.dataMatlabPath = fullfile(obj.dataCTPath,'matlab'); obj.dataDicomPath = fullfile(obj.dataCTPath,'dicom'); obj.shoulderManual.propagateDataPath; obj.shoulderAuto.propagateDataPath; 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.shoulderManual.scapula.glenoid.radius; % 2 Result.glenoid_radiusRMSE = obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE; % 3 % 3 %Result.glenoidSphericity = ' '; %Result.glenoid_biconcave = ' '; Result.glenoid_depth = obj.shoulderManual.scapula.glenoid.depth; % 4 Result.glenoid_width = obj.shoulderManual.scapula.glenoid.width; % 5 Result.glenoid_height = obj.shoulderManual.scapula.glenoid.height; % 6 Result.glenoid_center_PA = obj.shoulderManual.scapula.glenoid.centerLocal(1); % 7 Result.glenoid_center_IS = obj.shoulderManual.scapula.glenoid.centerLocal(2); % 8 Result.glenoid_center_ML = obj.shoulderManual.scapula.glenoid.centerLocal(3); % 9 Result.glenoid_version_ampl = obj.shoulderManual.scapula.glenoid.versionAmplitude; % 10 Result.glenoid_version_orient = obj.shoulderManual.scapula.glenoid.versionOrientation; % 11 Result.glenoid_Version = obj.shoulderManual.scapula.glenoid.version; % 12 Result.glenoid_Inclination = obj.shoulderManual.scapula.glenoid.inclination; % 13 Result.humerus_joint_radius = ' '; % 14 Result.humeral_head_radius = obj.shoulderManual.humerus.radius; % 15 % Result.humerus_GHsublux_2D = ' '; % Result.humerus_SHsublux_2D = ' '; Result.humerus_GHsubluxation_ampl = obj.shoulderManual.humerus.GHSAmpl; % 16 Result.humerus_GHsubluxation_orient = obj.shoulderManual.humerus.GHSOrient; % 17 Result.humerus_SHsubluxation_ampl = obj.shoulderManual.humerus.SHSAmpl; % 18 Result.humerus_SHsubluxation_orient = obj.shoulderManual.humerus.SHSOrient; % 19 Result.scapula_CSA = obj.shoulderManual.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.shoulderManual.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 = appendToCSV(obj,filename) % Save SCase to csv file logFid = fopen('log/measureSCase.log', 'a'); dataDir = ConfigFileExtractor.getVariable('dataDir'); xlsDir = [dataDir '/Excel/xlsFromMatlab']; fid = fopen([xlsDir '/' filename],'a'); fprintf(fid,[... obj.id ','... % SCase_id obj.shoulderManual.side ','... % shoulder_side num2str(obj.shoulderManual.scapula.glenoid.radius) ','... % glenoid_radius num2str(obj.shoulderManual.scapula.glenoid.fittedSphere.RMSE) ','... % glenoid_sphereRMSE num2str(obj.shoulderManual.scapula.glenoid.depth) ','... % glenoid_depth num2str(obj.shoulderManual.scapula.glenoid.width) ','... % glenoid_width num2str(obj.shoulderManual.scapula.glenoid.height) ','... % glenoid_height num2str(obj.shoulderManual.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA num2str(obj.shoulderManual.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS num2str(obj.shoulderManual.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML num2str(obj.shoulderManual.scapula.glenoid.versionAmplitude) ','... % glenoid_versionAmpl num2str(obj.shoulderManual.scapula.glenoid.versionOrientation) ','... % glenoid_versionOrientation num2str(obj.shoulderManual.scapula.glenoid.version) ','... % glenoid_version num2str(obj.shoulderManual.scapula.glenoid.inclination) ','... % glenoid_inclination num2str(obj.shoulderManual.humerus.jointRadius) ','... % humerus_jointRadius num2str(obj.shoulderManual.humerus.radius) ','... % humerus_headRadius num2str(obj.shoulderManual.humerus.GHSAmpl) ','... % humerus_GHSAmpl num2str(obj.shoulderManual.humerus.GHSOrient) ','... % humerus_GHSOrient num2str(obj.shoulderManual.humerus.SHSAmpl) ','... % humerus_SHSAmpl num2str(obj.shoulderManual.humerus.SHSOrient) ','... % humerus_SHSOrient num2str(obj.shoulderManual.humerus.SHSAngle) ','... % humerus_SHSAgle num2str(obj.shoulderManual.humerus.SHSPA) ','... % humerus_SHSPA num2str(obj.shoulderManual.humerus.SHSIS) ','... % humerus_SHSIS num2str(obj.shoulderManual.scapula.acromion.AI) ','... % acromion_AI num2str(obj.shoulderManual.scapula.acromion.CSA) ','... % acromion_CSA num2str(obj.shoulderManual.scapula.acromion.PSA) ','... % acromion_PSA num2str(obj.shoulderManual.scapula.acromion.AAA) '\n'... % acromion_AAA ]); fclose(fid); fclose(logFid); 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 end end diff --git a/ShoulderCase/@ShoulderCase/plot.m b/ShoulderCase/@ShoulderCase/plot.m index 26fa3cc..14d1827 100644 --- a/ShoulderCase/@ShoulderCase/plot.m +++ b/ShoulderCase/@ShoulderCase/plot.m @@ -1,89 +1,3 @@ function output = plot(obj,varargin) - % Call a ShoulderCase.shoulder.plot function. - % - % input: (''), ('manual'), ('auto') - % input arguments are not case sensitive. - % - % output: call obj.shoulderManual.plot (default) or obj.shoulderAuto.plot. - - if isempty(varargin) - % toggle everything - varargin = {'auto','manual','difference'}; - end - obj.updatePlotOptions(varargin{:}); - - - - % create new figure if necessary - if isempty(findall(obj.plotFigure)) - obj.createNewPlotFigure(); - end - - - - % clear case's figure - clf(obj.plotFigure); - - - - % update plot layout - obj.updatePlotLayout(); - - - - % update plot data - obj.plotChosenData(); - - if not(isempty(findall(obj.plotHandle.scapula))) - - - output = obj.plotOptions; -end - -% % manual/auto setup -% manual = obj.shoulderManual; -% manualColors = {'b','#D9FFFF'}; - -% auto = obj.shoulderAuto; -% autoColors = {'#D95319','#FFFF00'}; - - - -% % plot scapula wireframe -% manual.scapula.plotLandmarksWireframe(manualColors{1},manualColors{2}); -% auto.scapula.plotLandmarksWireframe(autoColors{1},autoColors{2}); - -% % plot scapula surface -% manual.scapula.plotSurface(manualColors{1},manualColors{2}); -% auto.scapula.plotSurface(autoColors{1},autoColors{2}); - -% % plot coordinate systems -% plot(manual.scapula.coordSys,manualColors{1},true); -% plot(auto.scapula.coordSys,autoColors{1},true); - -% % plot glenoid -% plot(manual.scapula.glenoid); -% plot(auto.scapula.glenoid); - -% % plot difference vector -% obj.plotManualAutoDifferences(10); - -% end - - - -% assert(nargin <= 2,'Too many arguments given') -% if nargin == 1 -% output = obj.shoulderManual.plot; -% return -% end -% assert(ischar(varargin{1}),'Given argument must be a char array') -% shoulderToPlot = lower(varargin{1}); % given argument is not case sensitive -% if isequal(shoulderToPlot,'manual') -% output = obj.shoulderManual.plot; -% elseif isequal(shoulderToPlot,'auto') -% output = obj.shoulderAuto.plot; -% else -% error('Wrong type of shoulder given. Input must be ''manual'' or ''auto''.'); -% end -% end \ No newline at end of file + ShoulderCasePlotter(obj); +end \ No newline at end of file diff --git a/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m b/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m index e9f4c45..d604465 100644 --- a/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m +++ b/ShoulderCase/@ShoulderCase/plotManualAutoDifferences.m @@ -1,66 +1,68 @@ -function plotManualAutoDifference(obj,normLimitForOutliers) +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,:)]; - plot3(differenceLine(:,1),differenceLine(:,2),differenceLine(:,3),... + plotHandle(i) = plot3(differenceLine(:,1),differenceLine(:,2),differenceLine(:,3),... ':x',... 'Color',color,... 'LineWidth',3); end + + output = plotHandle; end \ No newline at end of file