Diffusion ShoulderDB_ImportSegmentMeasureAnalyze (Bashardoust_SystCoord) Browse

ShoulderCase.m
8fd68c5307eaBashardoust_SystCoord

ShoulderCase.m
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-12-31

			% 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
			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

			rawSCaseID = SCaseIDParser(SCaseID);
			assert(rawSCaseID.isValidID,'The input argument is not a valid SCaseID.')
			obj.id = SCaseID;

			% Initialsation
			obj.diagnosis = '';
			obj.treatment = '';
			obj.outcome = '';
			obj.patient = Patient(obj);
			obj.shoulderManual = ShoulderManual(obj);
			obj.shoulderAuto = ShoulderAuto(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];


			%{
			%8 lines below to remove when -1 consrain avoided
			% 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];
			%}

			CTdirList=dir([SCaseDir '/CT-*']); % List directories with CT
			iCTdirAmira = 0;
			iCTdir2use = 0;
			for iCTdirList = length(CTdirList):-1:1 % Loop from last to first (4 to 1)
			CTdir = CTdirList(iCTdirList);
			% Check that dir name ends with -1,-2,-3,-4
			dirName = CTdir.name;
			if strcmp(dirName(end-1),'-') % Exclude postoperative 'p' CT
			CTnum = str2num(dirName(end));
			if CTnum <= 4 % Exlude non shoulder (elbow) CT
			% Check that the dir contains a dicom dir
			CTdir = [CTdir.folder '/' CTdir.name];
			dicomDir = [CTdir '/dicom'];
			if exist(dicomDir, 'dir')
			% Chech if amira dir exist
			amiraDir = [CTdir '/amira'];
			if exist(amiraDir, 'dir')
			% This is the CT folder to use
			iCTdirAmira = iCTdirList;
			end
			iCTdir2use = iCTdirList;
			end
			end
			end
			end
			if iCTdir2use == 0
			error(['\n', SCaseID, ' has no valid dicom directory']);
			else
			if iCTdirAmira % If amira dir exist, use it
			iCTdir2use = iCTdirAmira;
			end
			CTdir = CTdirList(iCTdir2use);
			CTdirPath = [CTdir.folder '/' CTdir.name];
			end
			CTDir = CTdirPath;

			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 = appendToCSV(obj,filename)
			% Save SCase to csv file
			logFid = fopen('log/measureSCase.log', 'a');

			dataDir = openConfigFile('config.txt', logFid);
			xlsDir = [dataDir '/Excel/xlsFromMatlab'];

			fid = fopen([xlsDir '/' filename],'a');
			fprintf(fid,[...
			obj.id ','... % SCase_id
			obj.shoulder.side ','... % shoulder_side
			num2str(obj.shoulder.scapula.glenoid.radius) ','... % glenoid_radius
			num2str(obj.shoulder.scapula.glenoid.sphereRMSE) ','... % glenoid_sphereRMSE
			num2str(obj.shoulder.scapula.glenoid.depth) ','... % glenoid_depth
			num2str(obj.shoulder.scapula.glenoid.width) ','... % glenoid_width
			num2str(obj.shoulder.scapula.glenoid.height) ','... % glenoid_height
			num2str(obj.shoulder.scapula.glenoid.centerLocal.x) ','... % glenoid_centerPA
			num2str(obj.shoulder.scapula.glenoid.centerLocal.y) ','... % glenoid_centerIS
			num2str(obj.shoulder.scapula.glenoid.centerLocal.z) ','... % glenoid_centerML
			num2str(obj.shoulder.scapula.glenoid.versionAmpl) ','... % glenoid_versionAmpl
			num2str(obj.shoulder.scapula.glenoid.versionOrient) ','... % glenoid_versionOrient
			num2str(obj.shoulder.scapula.glenoid.version) ','... % glenoid_version
			num2str(obj.shoulder.scapula.glenoid.inclination) ','... % glenoid_inclination
			num2str(obj.shoulder.humerus.jointRadius) ','... % humerus_jointRadius
			num2str(obj.shoulder.humerus.radius) ','... % humerus_headRadius
			num2str(obj.shoulder.humerus.GHSAmpl) ','... % humerus_GHSAmpl
			num2str(obj.shoulder.humerus.GHSOrient) ','... % humerus_GHSOrient
			num2str(obj.shoulder.humerus.SHSAmpl) ','... % humerus_SHSAmpl
			num2str(obj.shoulder.humerus.SHSOrient) ','... % humerus_SHSOrient
			num2str(obj.shoulder.humerus.SHSAngle) ','... % humerus_SHSAgle
			num2str(obj.shoulder.humerus.SHSPA) ','... % humerus_SHSPA
			num2str(obj.shoulder.humerus.SHSIS) ','... % humerus_SHSIS
			num2str(obj.shoulder.scapula.acromion.AI) ','... % acromion_AI
			num2str(obj.shoulder.scapula.acromion.CSA) ','... % acromion_CSA
			num2str(obj.shoulder.scapula.acromion.PSA) ','... % acromion_PSA
			num2str(obj.shoulder.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

			%{
			function outputArg = plot(obj)
			% code moved to file plot.m

			% 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

R8218/ShoulderCase/@ShoulderCase/ShoulderCase.m8fd68c5307eaBashardoust_SystCoord

ShoulderCase.mView Options

R8218/ShoulderCase/@ShoulderCase/ShoulderCase.m
8fd68c5307eaBashardoust_SystCoord

ShoulderCase.m
View Options