diff --git a/ShoulderCase/@ShoulderCase/ShoulderCase.m b/ShoulderCase/@ShoulderCase/ShoulderCase.m index bd6e9e8..a354736 100755 --- a/ShoulderCase/@ShoulderCase/ShoulderCase.m +++ b/ShoulderCase/@ShoulderCase/ShoulderCase.m @@ -1,499 +1,499 @@ 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 shoulder 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 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]; %{ %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 = saveCSV(~) % Save SCase to csv file - logFid = fopen('log/measureSCase.log', 'a'); + logFid = fopen('log/measureSCase.log', 'a'); - dataDir = openConfigFile('config.txt', logFileLoad); - xlsDir = [dataDir '/Excel/xlsFromMatlab']; + dataDir = openConfigFile('config.txt', logFid); + xlsDir = [dataDir '/Excel/xlsFromMatlab']; fid = fopen([xlsDir, '/anatomy.csv'],'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)... % 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