diff --git a/segmentScapula.m b/segmentScapula.m index e472a47..e85e97c 100755 --- a/segmentScapula.m +++ b/segmentScapula.m @@ -1,1401 +1,1401 @@ function outputArg = segmentScapula(varargin) %SEGMENTSCAPULA Segment scapula in /shoulder/data/{P\N} % Use statistical shape methods to segment the scapula. % The function can be run from (lbovenus) server as follows % cd /home/shoulder/methods/matlab/segmentation/method % matlab -nodisplay -nosplash -r "segmentScapula;quit" % MBO: WARNING on lbovenus, be sure to run the matlab command with a bash shell % (the result of the command 'echo $SHELL' should be '/bin/bash'). % The default shell on lbovenus accessed by PuTTY from windows is % /bin/tcsh which generates errors when calling: % matlab -nodisplay -nosplash -r "segmentScapula;quit" % To change the default shell from /bin/tcsh to /bin/bash write % "bash" in a ~/.tschrc file. Another way to do this is to call the % command "chsh -s bash" but that did not work for me. % The log file can be checked with shell command % cd /home/shoulder/methods/matlab/segmentation/method % tail -f log/segmentScapula.log % Input arguments can be one or several among: % {Empty, '[NP][1-999]', '[NP]','retryError', 'retryAll', 'prod'} % % Empty: will take into account all the shoulder cases (N and P). % '[NP][1-99]': will take into account these different shoulder cases. % '[NP]': will take into account all the specified shoulder cases (N or P). % 'retryError': will try again to do the segmentation of the scapula of the % chosen shoulder cases that formerly raised an error. % 'retryAll': will try again to do the segmentation of the scapula of all % the chosen shoulder cases. % 'prod': the segementation will be performed in '/home/shoulder/data' % instead of '/home/shoulder/dataDev'. % % Examples: % segmentScapula('prod') % is a valid instruction where every new case is segmented % in '/home/shoulder/data'. % segmentScapula('P400','N520') % is a valid instruction where the cases 'P400' and 'N520' are % segmented only if they have not already been segmented yet (new % cases). % segmentScapula('P400','N520','retryAll') % is a valid instruction where the cases 'P400' and 'N520' are % segmented whether they have already been segmented or not. % segmentScapula('P400','N520','P','retryError') % is a valid instruction where all the 'P' cases are segmented % if a former segmentation of these cases raised an error. % segmentScapula('Z1590','trollingArgumentLol') % is a valid instruction which raises a warning about the % arguments and will result in the same as executing % segmentScapula() i.e. every new case is calculated. % Output: % log file: /home/shoulder/methods/matlab/segmentation/method/log/segmentScapula.log % scapula surface: shoulder/data/{N/P}/{SCaseId}/matlab/scapulaSurfaceAuto.Ply % Associated functions % autoSegment: called to perform the scapula segmentation (original version of Elham Taghizadeh) % Innput: % Output: % IdentifyRightAndLeftScapula: scapula side identification (original version of Elham Taghizadeh) % Innput: % Output: % readDicomStack: in functions directory % read_mhd: in funcions directory % loadOpenFlipperFV: within this file (adapted from doMeasurement version of Elham Taghizadeh) % Innput: % Output: % simplePlyRead: (original version of Elham Taghizadeh, in functions dir) % Innput: % Output: % simplePlyWrite: (original version of Elham Taghizadeh, in functions dir) % Innput: % Output: % loadmodel: load SSM (original version of Elham Taghizadeh, in functions dir) % Innput: % Output: % compute_curvature: functions/toolbox_graph % Innput: % Output: % Author: Alexandre Terrier % Date: 2018-08-24 % Updated by: Matthieu Boubat % Date: 2019-08-27 % Comments by AT: % This code was intially written by Elham Taghizadeh (howToPrepareData), % located in /shoulder/methods/matlab/segmentation. It was then modififed % by Alexandre Terrier (howToPrepareData_AT), Bharath Narayanan % (howToPrepareData_BN), and again by Alexandre Terrier (segmentScapula.m). % The function save the scapular surface. % It takes about 5 day for the P cases % TODO: % Move to database folder % Find better reference to data diretory % Change main loop to consider N & P, might use listCAses % Use matlab dicom load % Add argument to force segement existing segmentations % Save STL too % Remove skippedRecord % princomp will be removed in a future release. Use pca instead. % Check P255, to avoid avoiding it... % Remove glenoid surface and scapula landmarks files if error or re-segmented % Check this warnign for P30 % Warning: keys is empty % > In extractSift3D (line 69) % In autoSegment (line 148) % In segmentScapula (line 433) % Define directories and add matlab paths -dataDir = '/shoulder/dataDev'; % Database databaseFunctions = '../../database'; % Database functions % Add database functions to matlab path addpath(genpath(databaseFunctions)); +dataDir = ConfigFileExtractor.getVariable('dataDir'); % Database % Path to segmentation directory segmentationCodePath = '../../segmentation/'; % Segmentation code modelPath = [segmentationCodePath 'Model/']; % Statistical Shape Models (right & left scapula) % Add segmentation functions to matlab path path = [segmentationCodePath 'method/functions']; addpath(genpath(path)); % Start log file filename = 'segmentScapula.log'; logDir = 'log'; filename = [logDir '/' filename]; logFileId = fopen(filename, 'w'); fprintf(logFileId, 'segmentScapula log file'); fprintf(logFileId, ['\n\nDate: ' datestr(datetime)]); % Parameters to redo segmentaions retryError = false; retryAll = false; specificCases = true; % Get parameters and list of SCases from varargin fprintf(logFileId, '\n\nList of SCase'); cases = {}; for argument = 1:nargin switch varargin{argument} % Test every argument in varargin. case 'prod' dataDir = '/shoulder/data'; case 'retryError' retryError = true; case 'retryAll' retryAll = true; case regexp(varargin{argument},'[NP]','match') specificCases = false; % Detect if argument is 'N' or 'P'. cases = {varargin{argument}}; % ends here. case regexp(varargin{argument},'[NP][1-9][0-9]{0,2}','match') % Test for a correct argument [NP][1-999] if specificCases cases{end+1} = varargin{argument}; % for element = 1:length(cases)-1 % Check if the case has already been added to the case list before if varargin{argument} == cases{element} % cases = cases(1:end-1); % Delete the last element added to cases if it is already in the list end end end otherwise warning(['\n "%s" is not a valid argument and has not been added to the'... ' list of cases.\n segmentScapula arguments format must be "[NP]",'... ' "[NP][1-999]", "retryError", "retryAll", or "prod".'],varargin{argument}) end end SCaseList = listSCase(cases); % Reset counters SCaseCount = 0; % Number of SCase treated (not counting skipped because already segmented) SCaseErrLoadDicomCount = 0; % Number of dicom loading error scapulaCount = 0; % Number of scapula treated (there might be 2 per SCase) scapulaErrSegmentCount = 0; % Number of segmentataion error % Parameters of the segmentation opt.isLeft = false; % AT Why do we keep this here? opt.doSmoothing = true; opt.ParallelProcess = true; opt.numTemplates = 5; opt.threshold = 190; opt.volThr = 35000; opt.ransacMaxIter = 5.0e6; opt.ransacDist = 6; % mm opt.ransacNo = 4; % number of samples for ransac! opt.searchWindowSize = [65, 65, 45]; % Add sift feature extraction toolbox to the path command = [segmentationCodePath 'Libraries/SIFT3D-master/build/lib/wrappers/matlab/setupSift3D.m']; run(command); % Load model and extract the sift features on the meanShape % Takes about 60 seconds fprintf(logFileId, '\n\nLoad statistical shape models'); fprintf(logFileId, '\n Right scapula'); modelLeft = loadModel(modelPath, true); % isLeft = true; fprintf(logFileId, '\n Left scapula'); modelRight = loadModel(modelPath, false); % isLeft = false; %%%% definition of analysis variables % alreadyDone = {}; alreadyDoneLoadDicom = {}; ok = {}; loadDicom = {}; ScapulaSide = {}; rightSegmentation = {}; leftSegmentation = {}; rightPlyFile = {}; leftPlyFile = {}; % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Start loop over sCases in sCaseList nSCase = length(SCaseList); % Number of sCases for iSCaseID = 1:nSCase save('segmentScapula_analysis.mat','ok','loadDicom','ScapulaSide',... 'rightSegmentation','leftSegmentation','rightPlyFile','leftPlyFile',... 'alreadyDone','alreadyDoneLoadDicom'); SCaseID = SCaseList(iSCaseID).id; % percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f'); SCaseDir = SCaseList(iSCaseID).dir; % Set data path of this sCase SCaseCount = SCaseCount+1; % Increment SCaseCount, might be removed since = nSCase matlabDir = [SCaseDir '/matlab']; % Check if this SCaseID is segmented scapulaSurfaceFile = ~isempty(dir([matlabDir '/scapulaSurfaceAuto*.ply'])); % Check if last segmentation produced an error errorFile = ~isempty(dir([matlabDir '/scapulaSurfaceAuto*Err.txt'])); % disp([SCaseID ' ' num2str(surfaceFile) ' ' num2str(errorFile)]); % For debug %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fprintf(logFileId, ['\n\n\n ___________________________________\n|\n| SCaseID: ' SCaseID ' (' percentProgress '%%)']); % Segment if % 1) retryAll is true % 2) Error file and retryError are true % 3) No scapulaSurface.ply file and no error file if ~(retryAll || (errorFile && retryError) || (~scapulaSurfaceFile && ~errorFile)) %% Load dicom CT fprintf(logFileId, '\n| Load CT'); dicomDir = [SCaseDir '/dicom']; [Img, imgInfo] = readDicomStack(dicomDir); % Might be replaced by matlab function if (isempty(Img)) % Write error file filename = [matlabDir '/scapulaSurfaceAutoErr.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); SCaseErrLoadDicomCount = SCaseErrLoadDicomCount+1; disp(SCaseID) alreadyDoneLoadDicom{end+1} = SCaseID; fprintf(logFileId, ': ERROR'); fprintf(logFileId,'\n|___________________________________\n\n') else alreadyDone{end+1} = SCaseID; fprintf(logFileId, ': OK'); fprintf(logFileId,'\n|___________________________________\n\n') end continue end if errorFile fprintf(logFileId, '\n| Already tried segmentation with error'); end while true % Start segmentation % Delete error files associated to the segmentation errorFiles = dir([matlabDir '/scapulaSurfaceAuto*Err.txt']); for iErrorFiles = 1:length(errorFiles) errorFile = [errorFiles(iErrorFiles).folder '/' errorFiles(iErrorFiles).name]; if exist(errorFile, 'file') delete errorFile; end end %% Load dicom CT fprintf(logFileId, '\n| Load CT'); dicomDir = [SCaseDir '/dicom']; [Img, imgInfo] = readDicomStack(dicomDir); % Might be replaced by matlab function if (isempty(Img)) fprintf(logFileId, ': ERROR'); % Write error file filename = [matlabDir '/scapulaSurfaceAutoErr.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); SCaseErrLoadDicomCount = SCaseErrLoadDicomCount+1; disp(SCaseID) loadDicom{end+1} = SCaseID; break % Skip to next SCase if dicom loading error else fprintf(logFileId, ': OK'); end %% Identify scapula side (takes about 30 seconds) fprintf(logFileId, '\n| Identify scapula side'); try [Left, Right] = IdentifyRightAndLeftScapula(Img, imgInfo, modelLeft, modelRight); % Img and imgInfo is similar to autoSegment function. fprintf(logFileId, ': OK'); catch fprintf(logFileId, ': ERROR'); % Write error file filename = [matlabDir '/scapulaSurfaceAuto' 'Err.txt']; message = 'Scapula side identification error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); ScapulaSide{end+1} = SCaseID; break; % Next for loop iteration (next Scase) end % Set scapula side if isempty(Right) if isempty(Left) % Neither right nor left scapulaSide = ''; fprintf(logFileId, '\n| Could not recognize any scapula side'); break; else % Left scapulaSide = 'L'; fprintf(logFileId, '\n| CT contains only the left scapula'); end elseif isempty(Left) % Right scapulaSide = 'R'; fprintf(logFileId, '\n| CT contains only the right scapula'); else % Both right and left scapulaSide = 'B'; fprintf(logFileId, '\n| CT contains both right and left scapulae'); end % End scapula side %% Segment scapula (takes 5-15 minutes) % This switch/case might be replaced by a simple block % We should save after each segmentation saveRight = 0; saveLeft = 0; switch scapulaSide case 'R' % Right fprintf(logFileId, '\n| Segment right scapula '); scapulaCount = scapulaCount+1; opt.isLeft = false; try tic; finalSurfaceRight = autoSegment(Img, imgInfo, modelRight, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveRight = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n| ' ME.identifier]); fprintf(logFileId, ['\n| ' ME.message]); fprintf(logFileId, '\n| Inconsistent data skipped'); % Write error file filename = [matlabDir '/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; rightSegmentation{end+1} = SCaseID; break; % Next iteration of for loop (next Scase) end case 'L' % Left fprintf(logFileId, '\n| Segment left scapula'); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = true; try tic; finalSurfaceLeft = autoSegment(Img, imgInfo, modelLeft, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveLeft = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n| ' ME.identifier]); fprintf(logFileId, ['\n| ' ME.message]); fprintf(logFileId, '\n| Inconsistent data skipped'); % Write error file filename = [matlabDir '/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter leftSegmentation{end+1} = SCaseID; break; % Next iteration of for loop (next Scase) --> Might be removed end case 'B' % Both right and left % Segment the right scapula fprintf(logFileId, '\n| Segment right scapula '); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = false; try tic; finalSurfaceRight = autoSegment(Right.img, Right.imgInfo, modelRight, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveRight = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n| ' ME.identifier]); fprintf(logFileId, ['\n| ' ME.message]); fprintf(logFileId, '\n| Inconsistent data skipped'); % Write error file filename = [matlabDir '/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error (R)'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter rightSegmentation{end+1} = SCaseID; break; end % Segment the left scapula fprintf(logFileId, '\n| Segment left scapula '); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = true; try tic; finalSurfaceLeft = autoSegment(Left.img, Left.imgInfo, modelLeft, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveLeft = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n| ' ME.identifier]); fprintf(logFileId, ['\n| ' ME.message]); fprintf(logFileId, '\n| Inconsistent data skipped'); filename = [matlabDir '/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error (L)'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter leftSegmentation{end+1} = SCaseID; break; end end % End of scapula segmentation %% Save segmented surface file(s) if saveRight % Save the right side filePly = [matlabDir '/scapulaSurfaceAutoR.ply']; fprintf(logFileId, '\n| Save scapula right surface file'); try simplePlyWrite(finalSurfaceRight, filePly); % Write scapula surface % Delete surface error if exists filename = strrep(filePly, '.ply', 'Err.txt'); if exist(filename, 'file') delete(filename); end fprintf(logFileId, ': OK'); catch % Write error file filename = [matlabDir '/scapulaSurfaceAutoRErr.txt']; message = 'Scapula write ply file error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); fprintf(logFileId, ': ERROR'); rightPlyFile{end+1} = SCaseID; break; end break; end if saveLeft % Save the left side filePly = [matlabDir '/scapulaSurfaceAutoL.ply']; fprintf(logFileId, '\n| Save scapula left surface file'); try simplePlyWrite(finalSurfaceLeft, filePly); % Write scapula surface % Delete surface error if exists filename = strrep(filePly, '.ply', 'Err.txt'); if exist(filename, 'file') delete(filename); end fprintf(logFileId, ': OK'); catch fprintf(logFileId,'\n| File writing failed'); % Write error file filename = [matlabDir '/scapulaSurfaceAutoLErr.txt']; message = 'Scapula write ply file error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); fprintf(logFileId, ': ERROR'); leftPlyFile{end+1} = SCaseID; break; end break; end % End of save surface file(s) ok{end+1} = SCaseID; end % End of SCase to segment fprintf(logFileId,'\n|___________________________________\n\n') end % End of SCaseList loop save('segmentScapula_analysis.mat','ok','loadDicom','ScapulaSide',... 'rightSegmentation','leftSegmentation','rightPlyFile','leftPlyFile',... 'alreadyDone','alreadyDoneLoadDicom'); % Message to log file fprintf(logFileId, ['\n\nNumber of SCase treated: ' num2str(SCaseCount)]); fprintf(logFileId, ['\n\nNumber of scapula treated: ' num2str(scapulaCount)]); fprintf(logFileId, ['\n\nNumber of dicom loading error: ' num2str(SCaseErrLoadDicomCount)]); fprintf(logFileId, ['\n\nNumber of segmentation error: ' num2str(scapulaErrSegmentCount)]); fprintf(logFileId, '\n\nEnd of segmentation'); fprintf(logFileId, ['\n\nDate: ' datestr(datetime)]); outputArg = 1; % Might be extended end %{ %% Loop through data folder, and 2 levels of directories % Reset counters SCaseCount = 0; % Number of SCase treated (not counting skipped because already segmented) SCaseErrLoadDicomCount = 0; % Number of dicom loading error scapulaCount = 0; % Number of scapula treated (there might be 2 per SCase) scapulaErrSegmentCount = 0; % Number of segmentataion error % Should be replaced by a list prodided by database/listSCase() location = 'P'; % Location - either Normal or Pathological for dirLevel1 = 0:9 for dirLevel2 = 0:9 directory = strcat(dataPath,location,'/',num2str(dirLevel1),'/',num2str(dirLevel2)); listOfPatients = dir(directory); % get all the folders within the current directory if(length(listOfPatients) > 2) % check if there are any patient data listOfPatients = listOfPatients(3:end); % if there are, store the list of patients nPatients = length(listOfPatients); % Loop through each patient within dirLevel1/dirLevel2 for iPatient = 1:nPatients subject = listOfPatients(iPatient).name; % Extract sCase.id (iCT) from subject iCT = strtok(subject,'-'); iCT = strtok(iCT, location); iCT = str2num(iCT); %#ok iCT = sprintf('%s%03d',location,iCT); dirPatient = strcat(directory,'/',subject,'/CT-',subject,'-1/'); % Should be renamed CTdir folderIn = dirPatient; matlabDir = [dirPatient '/matlab']; fprintf(logFileId, ['\n\nSCaseID: ' iCT]); % AT comment: P255 is non conventional and should be skipped if strcmp(iCT,'P255') fprintf(logFileId, ' Skip P255 until conform'); continue; end SCaseCount = SCaseCount+1; % Increment SCaseCount % Check if this SCaseID is segmented surfaceFile = ~isempty(dir([matlabDir '/scapulaSurfaceAuto*.ply'])); % Check if last segmentation produced an error errorFile = ~isempty(dir([matlabDir '/scapulaSurfaceAuto*Err.txt'])); disp([iCT ' ' num2str(surfaceFile) ' ' num2str(errorFile)]); if (~surfaceFile || (errorFile && retryError) ) % No scapulaSurface file or error during last % segmentation and required to redo % --> Start the segmentation % If errorFile w should delete it !! %% Loading of dicom files % AT: Sould be replaced by regular Matlab dicom load fprintf(logFileId, '\n Load CT'); % List of all images in the folder bones_in = dir([folderIn '*']); % AT: other way to do it in listSCase % Exclude "." and ".." from the list. Put dicom images in a folder! bones_in = bones_in(3:end); % Restrict to folder named "dicom", for use in /home/shoulder/data/ bones_in_idx = find(strcmp({bones_in.name},'dicom') == 1); bones_in = bones_in(bones_in_idx); % We should check it exists % AT comment: The dicom loading below might be repaced by % matlab function (from 2017 version) % Load CT images (takes about 5-15 seconds) for b = 1 : size(bones_in, 1) % Loop on dicom files (AT: directories rather than files) % The above loop might be removed since there is % only on dicom dir per CT fiolder. The dicom % loading below should be replaced by matlab % function % [V,spatial,dim] = dicomreadVolume(source) boneName = [folderIn bones_in(b).name]; % If the folder contains DICOMDIR or each image is in one directory! if (isdir(boneName)) imageFile = dir(boneName); imageFile = imageFile(3:end); for file_counter = 1:size(imageFile, 1) % Loop on dicom list if(strcmp(imageFile(file_counter).name(end-2:end), 'mhd') || strcmp(imageFile(file_counter).name(end-2:end), 'dcm')) try boneName = [boneName, filesep imageFile(3).name]; catch fprintf('\n Index exceeds matrix dimensions.'); skippedRecord{iSkipped,1} = iCT; skippedRecord{iSkipped,2} = 'Index exceeds matrix dimensions. Line 147'; iSkipped = iSkipped+1; continue; % Next iteration of for loop (next image) end break; % Terminate for loop (next image) end end % End of dicom file list end fname = boneName; ind = strfind(fname, '.'); if(isempty(ind)) % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); continue; % Next iteration of for loop (next SCase) end ind = ind(end); fileFormat = fname(ind+1:end); % right now I have code for ".am", ".mhd" and ".dcm" file formats. You % can add more cases here to read data. switch fileFormat case 'am' try [imgInfo, Img] = LoadData_Amira(boneName); catch fprintf(logFileId, ['\n ' bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); end case 'mhd' try [imgMHD, info]=read_mhd(boneName); catch fprintf(logFileId, ['\n ' bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); end eval(['Img = ' class(imgMHD.data) '(zeros(size(imgMHD.data, 2), size(imgMHD.data, 1), size(imgMHD.data, 3)));']); for ii = 1:size(imgMHD.data, 3) Img(:, :, ii) = imgMHD.data(:, :, ii)'; end imgInfo.BoundingBox.xMin = imgMHD.origin(1); imgInfo.BoundingBox.yMin = imgMHD.origin(2); imgInfo.BoundingBox.zMin = imgMHD.origin(3); imgInfo.PixelSpacing(1) = imgMHD.spacing(1); imgInfo.PixelSpacing(2) = imgMHD.spacing(2); imgInfo.SliceThickness = imgMHD.spacing(3); imgInfo.xSize = imgMHD.size(1); imgInfo.ySize = imgMHD.size(2); imgInfo.zSize = imgMHD.size(3); imgInfo.BoundingBox.xMax = imgMHD.origin(1) + (imgMHD.size(1)-1) * imgInfo.PixelSpacing(1); imgInfo.BoundingBox.yMax = imgMHD.origin(2) + (imgMHD.size(2)-1) * imgInfo.PixelSpacing(2); imgInfo.BoundingBox.zMax = imgMHD.origin(3) + (imgMHD.size(3)-1) * imgInfo.SliceThickness; if imgInfo.BoundingBox.zMax - imgInfo.BoundingBox.zMin > 700 fprintf(logFileId, ['\n ' bones_in(b).name ': skipped processing, due to the large size of scan, please crop this image and then try again (' num2str(size(imgMHD.data)) ')']); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error (image size > 700 pixels '; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); continue; % Next iteration of for loop (next SCase) end clear imgMHD; case 'dcm' % AT: this it what we have dicomDir = strsplit(boneName, filesep); dicomDir = strjoin(dicomDir(1, 1:end-1), filesep); [Img, imgInfo] = readDicomStack(dicomDir); if (isempty(Img)) fprintf(logFileId, ['\n ' bones_in(b).name ': Could not read the dicom images.']); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); % Added by AT % source = dicomDir; % [V,spatial,dim] = dicomreadVolume(source); % save([folderIn 'matlab/dicomreadVolume.mat'], V,spatial,dim); continue; % Next iteration of for loop (next SCase) % We should report an error here end otherwise % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Dicom loading error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); SCaseErrLoadDicomCount = SCaseErrLoadDicomCount+1; % Increment counter continue; % Next iteration of for loop (next Scase) end fname = bones_in(b).name; ind = strfind(fname, '.'); if(isempty(ind)) ind = length(fname)+1; else ind = ind(end); end % End of dicom loading % Identify scapula side (takes about 30 seconds) fprintf(logFileId, '\n Identify scapula side'); try [Left, Right] = IdentifyRightAndLeftScapula(Img, imgInfo, modelLeft, modelRight); % Img and imgInfo is similar to autoSegment function. catch fprintf('\n Scapula identification failed.\n'); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula side identification error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); continue; % Next iteration of for loop (next Scase) end % Set scapula side if isempty(Right) if isempty(Left) % Neither right nor left scapulaSide = ''; fprintf(logFileId, '\n Could not recognize scapulae side'); else % Left scapulaSide = 'L'; fprintf(logFileId, '\n CT contains only the left scapula'); end elseif isempty(Left) % Right scapulaSide = 'R'; fprintf(logFileId, '\n CT contains only the right scapula'); else % Both right and left scapulaSide = 'B'; fprintf(logFileId, '\n CT contains both left and right scapulae'); end %% Segment scapula (takes 5-15 minutes) % This switch/case might be replaced by a simple block saveRight = 0; saveLeft = 0; switch scapulaSide case 'R' % Right fprintf(logFileId, '\n Segment right scapula '); scapulaCount = scapulaCount+1; opt.isLeft = false; try tic; finalSurfaceRight = autoSegment(Img, imgInfo, modelRight, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveRight = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n ' ME.identifier]); fprintf(logFileId, ['\n ' ME.message]); fprintf(logFileId, '\n Inconsistent data skipped'); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; continue; % Next iteration of for loop (next Scase) end case 'L' % Left fprintf(logFileId, '\n Segment left scapula'); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = true; try tic; finalSurfaceLeft = autoSegment(Img, imgInfo, modelLeft, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveLeft = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n ' ME.identifier]); fprintf(logFileId, ['\n ' ME.message]); fprintf(logFileId, '\n Inconsistent data skipped'); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter continue; % Next iteration of for loop (next Scase) --> Might be removed end case 'B' % Both right and left % Segment the right scapula fprintf(logFileId, '\n Segment right scapula '); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = false; try tic; finalSurfaceRight = autoSegment(Right.img, Right.imgInfo, modelRight, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveRight = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n ' ME.identifier]); fprintf(logFileId, ['\n ' ME.message]); fprintf(logFileId, '\n Inconsistent data skipped'); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error (R)'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter end % Segment the left scapula fprintf(logFileId, '\n Segment left scapula '); scapulaCount = scapulaCount+1; % Increment counter opt.isLeft = true; try tic; finalSurfaceLeft = autoSegment(Left.img, Left.imgInfo, modelLeft, opt); fprintf(logFileId, [': Done in ' num2str(toc) ' s']); saveLeft = 1; catch ME fprintf(logFileId, ': ERROR'); fprintf(logFileId, ['\n ' ME.identifier]); fprintf(logFileId, ['\n ' ME.message]); fprintf(logFileId, '\n Inconsistent data skipped'); filename = [folderIn 'matlab/scapulaSurfaceAuto' scapulaSide 'Err.txt']; message = 'Scapula segmentation error (L)'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); scapulaErrSegmentCount = scapulaErrSegmentCount+1; % Increment counter end end % End of segmentation %% Write surface file if saveRight % Save the right side filePly = [folderIn 'matlab/scapulaSurfaceAutoR.ply']; fprintf(logFileId, '\n Save scapula right surface file'); try simplePlyWrite(finalSurfaceRight, filePly); % Write scapula surface % Delete surface error if exists filename = strrep(filePly, '.ply', 'Err.txt'); if exist(filename, 'file') delete(filename); end fprintf(logFileId, ': OK'); catch % Write error file filename = [folderIn 'matlab/scapulaSurfaceAutoRErr.txt']; message = 'Scapula write ply file error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); fprintf(logFileId, ': ERROR'); continue; end end if saveLeft % Save the left side filePly = [folderIn 'matlab/scapulaSurfaceAutoL.ply']; fprintf(logFileId, '\n Save scapula left surface file'); try simplePlyWrite(finalSurfaceLeft, filePly); % Write scapula surface % Delete surface error if exists filename = strrep(filePly, '.ply', 'Err.txt'); if exist(filename, 'file') delete(filename); end fprintf(logFileId, ': OK'); catch fprintf(logFileId,'\n File writing failed'); % Write error file filename = [folderIn 'matlab/scapulaSurfaceAutoLErr.txt']; message = 'Scapula write ply file error'; fileID = fopen(filename, 'w'); fprintf(fileID, message); fclose(fileID); fprintf(logFileId, ': ERROR'); continue; end end end % End for loop over dicom directories (in practice there is only one) else fprintf(logFileId, '\n Segmentation already done'); end % End of if segementation exists end % End of SCase loop end % End of if any patient data end % End of dirLevel2 end % End of of dirLevel 1 % Message to log file fprintf(logFileId, ['\n\nNumber of SCase treated: ' num2str(SCaseCount)]); fprintf(logFileId, ['\n\nNumber of scapula treated: ' num2str(scapulaCount)]); fprintf(logFileId, ['\n\nNumber of dicom loading error: ' num2str(SCaseErrLoadDicomCount)]); fprintf(logFileId, ['\n\nNumber of segmetation error: ' num2str(scapulaErrSegmentCount)]); fprintf(logFileId, '\n\nEnd of segmentation'); fprintf(logFileId, ['\n\nDate: ' datestr(datetime)]); outputArg = 1; % Might be extended end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % All below might be removed %% Function to get landmarks, adapted from doMeasurements of Elham Taghizadeh function out = getLandmarks(fv, options) % it works only for the mesh that was segmented using our segmentation % method that provides correspondences to our statistical shape model. % fv.faces: contains the connectivity information of the surface mesh (25000 x 3) % fv.vertices: information on the vertices of the mesh (12502 x 3) % options stores some important information regarding the bone: % .side Can be 'left' or 'right' % .curvature_smoothing Curvature smoothing parameter. Default value is 10 % .verb step-by-step reporting. by default it is set to false % out is a struct contains different measurements %% check the input's validity if(~isfield(fv, 'faces')) error('The provided mesh is not in correct format: the field .faces is not provided'); end if(~isfield(fv, 'vertices')) error('The provided mesh is not in correct format: the field .vertices is not provided'); end sz_F = size(fv.faces); sz_V = size(fv.vertices); if (numel(sz_F)~=2 || sum(sz_F == [25000, 3])~=2 || numel(sz_V)~=2 || sum(sz_V == [12502, 3])~=2) error('The provided mesh is not in correct format.'); end % AT: values of default options if not provided if (nargin == 1) options.side = 'right'; options.curvature_smoothing = 10; options.verb = 0;% 0: do not print progress bar; 1: print progress bar in the consol else if(~isfield(options, 'side')) options.side = 'right'; end if(~isfield(options, 'curvature_smoothing')) options.curvature_smoothing = 10; end if(~isfield(options, 'verb')) options.verb = 0; end end % AT: if side is right, do a relexion, for left? Ask Elham! if(strcmp(options.side, 'right')) reflect = true; else reflect = false; end %% The important vertices needed for measurements % AT: Ask Elham wha is ScapulaToExclude and GlenoidToExclude [selectedVertices.ScapulaToExclude, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_ScapulaToExclude.ini']); [selectedVertices.ScapulaGroove, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Groove.ini']); [selectedVertices.ScapulaPillar, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Pillar.ini']); [selectedVertices.GlenoidToExclude, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_GlenoidToExclude.ini']); [selectedVertices.Spinoglenoidnotch, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_Spinoglenoidnotch.ini']); [selectedVertices.GlenoidSurface, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_GlenoidSurface.ini']); selectedVertices.ScapulaExtremities = load(['.' filesep '00_vertices' filesep 'indiceVertices_ScapulaExtremities.txt']); Sref = simplePlyRead(['.' filesep '00_vertices' filesep 'ReferenceShape.ply']); %% Prepare sample for measurements % AT: Is S for sample, F for faces and V for vertices? S.FV = fv; % Line below aleady commented by Elham %AT: What is T? % [~, S.FV.vertices] = procrustes(Sref.vertices, fv.vertices, 'scaling', false, 'reflection', reflect); [~, S.FV.vertices, T] = procrustes(Sref.vertices, fv.vertices, 'scaling', false, 'reflection', reflect); %try to compute the curvature, if it does not work, we slightly smooth %the surface and retry temp = 1; while(temp) try temp = 0; [S.Curvature.DirectionMinCurvature,... S.Curvature.DirectionMaxCurvature,... S.Curvature.MinCurvature,... S.Curvature.MaxCurvature,... S.Curvature.MeanCurvature,... S.Curvature.GaussianCurvature,... S.Curvature.Normal] = ... compute_curvature(S.FV.vertices, S.FV.faces, options); catch %...if an error occurs, we smooth the face and recompute the curvature S.FV = smoothpatch(S.FV, 1, 1, 0.0001); if(options.verb) disp('\tThe registered shape is smoothed to compute the curvature\n'); end temp = 1; end end %% find Landmarks based on Curvatures temp_FV = S.FV; temp_Curvature = S.Curvature; temp_Landmarks.Extremities = temp_FV.vertices(selectedVertices.ScapulaExtremities,:); % ----------- find the spinoglenoid notch ----------------- %select the vertice which has the smallest gaussian curvature [~,temp] = min(temp_Curvature.GaussianCurvature(selectedVertices.Spinoglenoidnotch)); temp = selectedVertices.Spinoglenoidnotch(temp); temp_Landmarks.Extremities(6,:) = temp_FV.vertices(temp,:); % ------ find the scapula groove landmarks ------ % AT: Choosing "best" vertices might cause error. % This could be replaced by finding "new" points that best match the curvature criteria temp_dist = [12 44 28 20 36]; skipInd = []; for iL=1:5 %get all the candidate vertices temp_indices = selectedVertices.ScapulaGroove; temp_vertices = temp_FV.vertices(temp_indices,:); %get the vector front-scapula->back-scapula [temp_n,~] = fitLineToPoints(temp_vertices); temp_n = temp_n'; temp_Lfront = temp_Landmarks.Extremities(6,:); %from the front-scapula, go X mm in direction of the back, and remove %all candidate vertices behind temp_origine = temp_Lfront + (temp_dist(iL)-2)*temp_n; for i=length(temp_indices):-1:1 temp_p = temp_vertices(i,:)-temp_origine; if(dot(temp_p,temp_n)<0) temp_indices(i)=[]; temp_vertices(i,:)=[]; end end %from the front-scapula, go Y mm in direction of the back, and remove %all candidate vertices ahead temp_origine = temp_Lfront + (temp_dist(iL)+2)*temp_n; for i=length(temp_indices):-1:1 temp_p = temp_vertices(i,:)-temp_origine; if(dot(temp_p,temp_n)>0) temp_indices(i)=[]; temp_vertices(i,:)=[]; end end %get the best vertice based on the min curvature if ~isempty(temp_indices) [temp, i] = min(temp_Curvature.MinCurvature(temp_indices)); temp_Landmarks.Groove(iL,:) = temp_vertices(i,:); else skipInd = [skipInd iL]; end end if(numel(skipInd) == 5) disp('The shape is distorted! Could not find enough landmarks on the Groove'); return; else temp_Landmarks.Groove(skipInd, :) = []; end % --------- find the scapula pillar landmarks ------ temp_dist = [30 58 44 37 51]; for iL=1:5 %get all the candidate vertices temp_indices = selectedVertices.ScapulaPillar; temp_vertices = temp_FV.vertices(temp_indices,:); %get the vector up->down [temp_n,~] = fitLineToPoints(temp_vertices); temp_n = temp_n'; %check the vector point in the right direction temp_up = temp_Landmarks.Extremities(6,:); if(dot(temp_n, temp_FV.vertices(selectedVertices.ScapulaPillar(1,:))-temp_up)<0) temp_n =-temp_n; end %from the front-scapula, go X mm in direction of the back, and remove %all candidate vertices behind temp_origine = temp_up + (temp_dist(iL)-2)*temp_n; for i=length(temp_indices):-1:1 temp_p = temp_vertices(i,:)-temp_origine; if(dot(temp_p,temp_n)<0) temp_indices(i)=[]; temp_vertices(i,:)=[]; end end %from the front-scapula, go Y mm in direction of the back, and remove %all candidate vertices ahead temp_origine = temp_up + (temp_dist(iL)+2)*temp_n; for i=length(temp_indices):-1:1 temp_p = temp_vertices(i,:)-temp_origine; if(dot(temp_p,temp_n)>0) temp_indices(i)=[]; temp_vertices(i,:)=[]; end end %get the best vertice based on the max curvature [temp, i] = max(temp_Curvature.MaxCurvature(temp_indices)); temp_Landmarks.Pillar(iL,:) = temp_vertices(i,:); end % ------------- save everyting (and glenoid surface) in the structure ---------------- temp_v = selectedVertices.GlenoidSurface(temp_Curvature.MaxCurvature(selectedVertices.GlenoidSurface)<1); S.indiceVertices.GlenoidSurface = temp_v; S.Landmarks = temp_Landmarks; %% compute the glenoid orientations % AT: Removed %{ S.GlenoidOrientation = ... computeGlenoidOrientation(... S.FV,... S.indiceVertices.GlenoidSurface,... S.Landmarks); %} % AT: These final lines need to be better understood... out = rmfield(S, 'FV'); % AT: Remove field FV from structure S and put in out out = rmfield(out, 'indiceVertices'); % AT: What is T out.Landmarks.Pillar = bsxfun(@minus, out.Landmarks.Pillar, T.c(1, :)) * T.T'; out.Landmarks.Extremities = bsxfun(@minus, out.Landmarks.Extremities, T.c(1, :)) * T.T'; out.Landmarks.Groove = bsxfun(@minus, out.Landmarks.Groove, T.c(1, :)) * T.T'; % 3 lines below added by AT !!!Error disp('\nStarting the identification of the glenoid surface'); GlenoidSurface.vertices = S.FV.vertices(S.indiceVertices.GlenoidSurface,:); GlenoidSurface.faces = getFaces(S.indiceVertices.GlenoidSurface,S.FV); displ(GlenoidSurface); out.GlenoidSurface = GlenoidSurface; end %% Load vertices numbers function [iVertices, iFaces] = loadOpenFlipperFV(filepath) %% read the file fileID = fopen(filepath,'r'); temp = textscan(fileID,'%s','Delimiter','\n'); fclose(fileID); %% extract vertices %select the line which list all the vertices of interest temp_Vertices = temp{1}{7}; temp_Vertices(1:16) = []; %remove the text 'vertexSelection=' if(isempty(temp_Vertices)==0) temp_Vertices(end) = []; %remove the last character %split the string temp_Vertices = strsplit(temp_Vertices,';')'; %transform the string in nummer iVertices = zeros(length(temp_Vertices),1); for i=1:length(temp_Vertices) iVertices(i,1) = str2num(temp_Vertices{i}) + 1; %add 1 because the vertices start at 0 in OpenFlipper and 1 in Matlab end else iVertices = []; end %% extract faces %select the line which list all the vertices of interest temp_Faces = temp{1}{3}; temp_Faces(1:14) = []; %remove the text 'faceSelection=' if(isempty(temp_Faces)==0) temp_Faces(end) = []; %remove the last character %split the string temp_Faces = strsplit(temp_Faces,';')'; %transform the string in nummer iFaces = zeros(length(temp_Faces),1); for i=1:length(temp_Faces) iFaces(i,1) = str2num(temp_Faces{i}) + 1; %add 1 because the vertices start at 0 in OpenFlipper and 1 in Matlab end else iFaces = []; end end %% Function by Bharath Narayanan function [faceSubSet] = getFaces(indexVertices,FV) % This function outputs all the faces associated with the glenoid % vertices. Note that it outputs the 'local' or 'glenoid' vertex % indices. vertexSubSet = indexVertices; faceSuperSet = FV.faces; faceCount = 1; for i = 1:length(faceSuperSet) face = faceSuperSet(i,:); % Check if all vertices of the existing face are part of the glenoid if(~isempty(find(vertexSubSet == face(1)))) && (~isempty(find(vertexSubSet == face(2)))) && (~isempty(find(vertexSubSet == face(3)))) % If so, determine the local vertex indices of these vertices in the glenoid % frame. index1 = find(vertexSubSet == face(1)); index2 = find(vertexSubSet == face(2)); index3 = find(vertexSubSet == face(3)); faceSubSet(faceCount,:) = [index1,index2,index3]; faceCount = faceCount + 1; end end end %% Glenoid orienation: Could be removed function Results = computeGlenoidOrientation(FV, indiceVertices_GlenoidSurface, Landmarks) % *********************************************************** % Notes: landmarks.extremities % (1) Bottom of the scapula % (2) Back of the scapula % (3) Tip of acromion % (4) Top coracoid % (5) Back coracoid % (6) Front scapula % ************************************************************ %load Data % coordinate of vertices on the Glenoid Surface GlenoidSurface.vertices = FV.vertices(indiceVertices_GlenoidSurface,:); % ----------------------- compute the coordinate system of the scapula ----------------------- % Find the transverse scapular plane (fit plane to landmarks.ScapulaPillar and landmarks.ScapulaGroove) [ScapulaPlane.Normal, ScapulaPlane.Center] = fitPlaneToPoints([Landmarks.Groove; Landmarks.Pillar]); %define the X coord syst CSYS.X = ScapulaPlane.Normal; %check the X axis goes in the right direction temp = Landmarks.Extremities(3,:) - Landmarks.Extremities(4,:); if(sign(dot(CSYS.X,temp))<0) CSYS.X = -CSYS.X; end %define the Z coord syst temp = projectPointsToPlane(Landmarks.Groove,ScapulaPlane.Normal',[0 0 0]); % project landmarks.ScapulaGroove points on scapular plane [CSYS.Z, ~] = fitLineToPoints(temp); % fit a line on the projected landmarks.ScapulaGroove CSYS.Z = CSYS.Z/norm(CSYS.Z); %normalize the vector %check the Z axis goes in the right direction temp = Landmarks.Extremities(6,:) - mean(Landmarks.Groove); if(sign(dot(CSYS.Z,temp))<0) CSYS.Z = -CSYS.Z; end %define the Y coord syst CSYS.Y = cross(CSYS.Z,CSYS.X); % Y % ----------------------- compute the Glenoid center ----------------------- % temp = zeros(size(GlenoidSurface.Vertices)); % temp_distance3D = zeros(size(GlenoidSurface.Vertices,1),1); GlenoidSurface.Center = mean(GlenoidSurface.vertices); %returns the indices of the closest points in X for each point in Y [temp_iSelectedVertice, ~]= dsearchn(GlenoidSurface.vertices,GlenoidSurface.Center); GlenoidSurface.Center = GlenoidSurface.vertices(temp_iSelectedVertice,:); % ----------------------- Fit Sphere on Glenoid Surface ----------------------- [GlenoidSphere.Center,GlenoidSphere.Radius,temp_GlenoidResiduals] = fitSphereToPoints(GlenoidSurface.vertices); GlenoidSphere.RMSE = norm(temp_GlenoidResiduals)/sqrt(length(temp_GlenoidResiduals)); % Glenoid version 3D (spherical coordinates) %centerline Glenoid.Centerline = GlenoidSphere.Center'-GlenoidSurface.Center; %get the centerline vector in the scapula coord system tempX = Glenoid.Centerline * CSYS.X; %project the centerline on the X axis of the scapula tempY = Glenoid.Centerline * CSYS.Y; %project the centerline on the Y axis of the scapula tempZ = Glenoid.Centerline * CSYS.Z; temp_r = norm(Glenoid.Centerline); %length of the centerline %version Glenoid.Version_inDeg = acosd(tempZ/temp_r); %version orientation Glenoid.VersionOrientation_inDeg = 180-rad2deg(atan2(tempY,tempX)); if(Glenoid.VersionOrientation_inDeg > 180) Glenoid.VersionOrientation_inDeg = 360 - Glenoid.VersionOrientation_inDeg; end % ----------------------- Output ----------------------- % Glenoid Results.Glenoid_Center = GlenoidSurface.Center; Results.Glenoid_SphereCenter = GlenoidSphere.Center; Results.Glenoid_SphereRadius = GlenoidSphere.Radius; Results.Glenoid_Centerline = Glenoid.Centerline; Results.Glenoid_VersionInDegree = Glenoid.Version_inDeg; Results.Glenoid_VersionOrientationInDegree = Glenoid.VersionOrientation_inDeg; Results.CoordSys.X = CSYS.X; Results.CoordSys.Y = CSYS.Y; Results.CoordSys.Z = CSYS.Z; end %}