diff --git a/README.md b/README.md new file mode 100755 index 0000000..a36101b --- /dev/null +++ b/README.md @@ -0,0 +1,52 @@ +# Content of segmentation directory + +This README describes the content of the segmentation directory, as the main workflow to perform segmentation of shoulder CT, identification of scapula landmarks, and glenoid surface. + +Author: Alexandre Terrier, EPFL-LBO + +Creation date: 2018-11-27 + +## Updated functions by Alexandre Terrier + +The following functions should be run in this order to get the segmentation of the scapula (SCases in shoulder/data), the bony landmarks, and the glenoid surface. The anatomical measurements are performed from these data, using matlab scripts located in /shoulder/method/matlab/database. + +* segmentScapula.m: function created by Alexandre Terrier (from Elham Taghizadeh code and Bharath Narayanan extension) which calculates the segmentation of all scapulae in the database of CT. It saves the result as a ply file in in case directory. There is a log file in the log directory. + +* findScapulaLandmarksSCase.m: uses findScapulaLandmarks.m for to get bone landmarks o all SCases. There is a log file in the log directory. + +* findGlenoidSurfaceSCase.m: uses findGlenoidSurface.m to get the glenoid surface of all SCases. There is a log file in the log directory. + +* TO DO + * segmentScapula.m should be divided in 2 functions: one for the list of SCase and one for one case, as for findScapulaLandmarksSCase.m and findGlenoidSurfaceSCase.m. + * unused files should be removed + * "SC" functions might be moved to database directory + + +## Initial functions by Elham Taghizadeh + +* howToPrepareData.m + +* howToPrepareData_AT.m + +* “autoSegment.m” is the function that does the segmentation. You should provide the loaded image with some input options (information on the header of functions) and as output it will provide the surface mesh as a struct with the fields of .vertices (coordinates of the nodes) and .faces (connectivity information). + +* “sampleData” contains an example image for segmentation. + +* “functions” folder contains all functions I use for segmentation. + +* “00_vertices” folder contains all preselected nodes of the mesh that are going to be used for automatic landmark detection. + +* “loadModel.m” load the model and all information that is necessary for segmentation. + +* “doMeasurements.m” is the function that gets the segmented bone as input and as output provides the measurements in a struct. + +* To see how the model and images should be prepared, you can find a MATLAB script in the same folder “howToPrepareData.m” that represents how you can call different functions. + +* The files with the name of “tmp” are saved during process of segmentation and are replaced for every call to “autoSegment” function. + +* "Protocol Segmentation Scapula.docx" is a protocol (format word) which explains how to segmented a scapula from the beginning to the end. It is situated in the folder: "shoulder\methods\documentations\protocols" + + +## Codes by Bharath Narayanan + +The directory bharath2remove contains codes developed by Bharath Narayanan, who extended codes of Elham Taghizadeh to work on the shoulder database. These codes have been updated in segmentScapula.m and findScapulaLandmark.m by Alexandre Terrier. diff --git a/findGlenoidSurface.m b/findGlenoidSurface.m new file mode 100755 index 0000000..da13261 --- /dev/null +++ b/findGlenoidSurface.m @@ -0,0 +1,112 @@ +function [GlenoidSurface,message] = findGlenoidSurface(filename) +%FINDGLENOIDSURFACE Finf the glenoid surface of a scapula surface +% Detailed explanation goes here + + +% Input + +% Output + +% Example + +% Author: Alexandre Terrier, EPFL-LBO + +% Date: 2018-09-03 + +% TODO +% options vraiable is missing. It is really needed, more than default +% below? + + + + +% 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 + + +% Add path to segmenation functions folder +path2Add = genpath('./functions'); +addpath(path2Add); + + %Default options +options.side = 'right'; % ensures that the scapula is on the right side by default. +options.curvature_smoothing = 10; +options.verb = 0;% 0: do not print progress bar; 1: print progress bar in the consol + +FV = simplePlyRead(filename); % read the input ply file + +[selectedVertices.GlenoidSurface, ~] = loadOpenFlipperFV(['.' filesep '00_vertices' filesep 'OpenFlipper_GlenoidSurface.ini']); + +Sref = simplePlyRead(['.' filesep '00_vertices' filesep 'ReferenceShape.ply']); +if(strcmp(options.side, 'right')) + reflect = true; +else + reflect = false; +end + +% Prepare sample for measurements +S.FV = FV; +[~, 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; % Vertices and faces +temp_Curvature = S.Curvature; % Curvatures at vertices +temp_v = selectedVertices.GlenoidSurface(temp_Curvature.MaxCurvature(selectedVertices.GlenoidSurface)<1); % Restrict node selection according to curvature +S.indexVertices.GlenoidSurface = temp_v; % indices corresponding to the glenoid surface vertices + + +GlenoidSurface.vertices = FV.vertices(S.indexVertices.GlenoidSurface,:); +GlenoidSurface.faces = getFaces(S.indexVertices.GlenoidSurface,FV); + +message = 'OK'; + +end + +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 \ No newline at end of file diff --git a/findGlenoidSurfaceSCase.m b/findGlenoidSurfaceSCase.m new file mode 100755 index 0000000..cdc2a0a --- /dev/null +++ b/findGlenoidSurfaceSCase.m @@ -0,0 +1,93 @@ +function outputArg = findGlenoidSurfaceSCase(varargin) +%FINDGLENOIDSURFACESCASE Find glenoid surface of segemented SCase +% Detailed explanation goes here + +% The function can be run from (lbovenus) server as follows +% cd /home/shoulder/methods/matlab/segmentation/method +% matlab -nodisplay -nosplash -r "findGlenoidSurfaceSCase;quit" + +% The log file can be checked with shell command +% cd /home/shoulder/methods/matlab/segmentation/method +% tail -f log/findGlenoidSurfaceSCase.log + + +% Input + +% Output + +% Example + +% Author: Alexandre Terrier, EPFL-LBO + +% Date: 2018-09-03 + +% TODO + +tic; % Start stopwatch timer + +% Open log file +filename = 'findGlenoidSurfaceSCase.log'; +logDir = 'log'; +filename = [logDir '/' filename]; +logFileId = fopen(filename, 'w'); +fprintf(logFileId, 'findGlenoidSurfaceSCase log file'); +fprintf(logFileId, ['\n\nDate: ' datestr(datetime)]); + +% Add path to database functions folder +path2Add = '../../database'; +addpath(path2Add); +% Add path to segmenation functions folder +path2Add = genpath('./functions'); +addpath(path2Add); + +% Get list of sCases from varargin +fprintf(logFileId, '\n\nList of SCase'); +if nargin == 0 % findScapulaLandmarksSCase called without argument + SCaseList = listSCase; +elseif nargin == 1 + inputArg = varargin{1}; + SCaseList = listSCase(inputArg); +else + error('inputArg can be empty, N, P, or a SCaseID') +end + +% Loop on list of SCase +nSCase = length(SCaseList); % Number of sCases +for iSCase = 1:nSCase + SCaseID = SCaseList(iSCase).id; + fprintf(logFileId, ['\n\nSCaseID: ' SCaseID]); + currentSCaseDir = SCaseList(iSCase).dir; + matlabDir = [currentSCaseDir '/matlab']; + % Find ply files of scapula surface in matlab dir (there might be 0, 1 or 2) + filenameList = dir([matlabDir '/scapulaSurfaceAuto*.ply']); + if isempty(filenameList) + fprintf(logFileId, '\n No scapulaSurface ply file'); + continue + end + for iFilename = 1:length(filenameList) % Loop on segmented scapula ply files (L/R) + filename = filenameList(iFilename); + filename = [filename.folder '/' filename.name]; + try + [GlenoidSurface,message] = findGlenoidSurface(filename); % Get glenoidSurface + catch + fprintf(logFileId, '\n Error in findGlenoidSurface'); + continue; + end + if isstruct(GlenoidSurface) + % Save glenoid surfaace + if message + fprintf(logFileId, ['\n ' message]); % We might add here R/L/- + end + filename = strrep(filename, 'scapula', 'glenoid'); % Replace Scapula by Glenoid + simplePlyWrite(GlenoidSurface, filename); + fprintf(logFileId, ['\n Save ' filename]); + else + fprintf(logFileId, '\n Error finding landmarks'); + end + end +end % End of SCase loop + +fprintf(logFileId, ['\n\n End of findScapulaLandmarksSCase after ' num2str(toc) ' seconds']); +fclose(logFileId); % Close log file +outputArg = 1; +end diff --git a/findScapulaLandmarks.m b/findScapulaLandmarks.m new file mode 100755 index 0000000..0c675e7 --- /dev/null +++ b/findScapulaLandmarks.m @@ -0,0 +1,330 @@ +function [ScapulaLandmarks, message] = findScapulaLandmarks(fv, options) +%FINDSCAPULALANDMARKS Find scapula landmarks from scapula surface file ply +% Detailed explanation goes here + +% Can only be run from lbovenus to avoid error +% Undefined function or variable 'smoothpatch_curvature_double' +% in function smoothpatch in function findScapulaLandmarks + +% + + +% 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 + + +% The function merge the part of doMeasuremnts (written by Elham +% Taghizadeh) that get the scapula landmarks. The anatomical measurements +% are no more performed here, but in database/measureSCase. The function +% should save the scapular surface, the glenoid surface, and the scapula +% landmarks. + +% This is done by another function findScapulaLandmarks + +% Comment by AT: This should be adapted to the +% ShoulderCase class We should add below the +% identification of the scapula landmarks and the +% glenoid surface. The scapula landmarks are currently +% embeded in function doMeasurements (by ET). Glenoid +% surface is currently obtained as a second step with +% function extractAllGlenoid and extractGlenoidSurface +% (by BN). + +% Variable finalSurface is the argument of function +% doMeasurements to get landmarks and perform +% anatomical measurements. It might be replaced by a +% function getLanmarks. + +% *********************************************************** +% Notes: landmarks.extremities (From Elham T.) +% (1) Bottom of the scapula --> angulusInferior +% (2) Back of the scapula --> trigonumSpinae +% (3) Tip of acromion --> processusCoracoideus ? +% (4) Top coracoid --> acromioClavicular ? +% (5) Back coracoid --> angulusAcromialis ? +% (6) Front scapula --> spinoGlenoidNotch +% ************************************************************ + + +% Perform anatomical data save in this function instead of bharath2remove/extractAllGlenoid.m + + +% 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: + +% TODO +% Check why goove landmarks can be less than 5 + +message = ''; % Added by AT + +% Add path to segmentation functions +%path2Add = 'functions'; +%addpath(path2Add); +addpath(genpath('./functions')); + + +% 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 what 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']); + + +%% "Choose" landmarks based on curvature + +% AT: Choosing "best" vertices might cause error. +% This could be replaced by finding "new" points that best match the curvature criteria + +% AT: S for shape, F for faces and V for vertices +S.FV = fv; +% Line below aleady commented by Elham +% [~, S.FV.vertices] = procrustes(Sref.vertices, fv.vertices, 'scaling', false, 'reflection', reflect); +% AT: procrustes(X,Y) = linear transformation between X and Y, returns the +% tranformed data and the transformation T. +[~, 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 ------ +temp_dist = [12 44 28 20 36]; % AT: Are these numbers distance in mm? It might cause problem for partial scapulae +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) + message = 'The shape is distorted! Could not find enough landmarks on the Groove'; + temp_Landmarks.Groove = temp_vertices; % Added by AT + % return; % Commented by AT +else + try % AT added try/catch because of error in line below + temp_Landmarks.Groove(skipInd, :) = []; + catch + message = 'Less than 5 points remaing in groove after curvature optimisation'; + % Take back initial values + temp_Landmarks.Groove = temp_vertices; + end + +end + +% --------- find the scapula pillar landmarks ------ +temp_dist = [30 58 44 37 51]; % AT: Same as for groove +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)); + try % AT added try/catch because of error in line below + temp_Landmarks.Pillar(iL,:) = temp_vertices(i,:); + catch + message = [message ' Less than 5 points remaing in pillar after curvature optimisation']; + % Take back initial values + temp_Landmarks.Pillar = temp_vertices; + end +end + +% ------------- save everyting (and glenoid surface) in the structure S---------------- +temp_v = selectedVertices.GlenoidSurface(temp_Curvature.MaxCurvature(selectedVertices.GlenoidSurface)<1); +S.indiceVertices.GlenoidSurface = temp_v; +S.Landmarks = temp_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? Transformation obtained by procrustes (T.c = translation) +scapulaExtremities = bsxfun(@minus, out.Landmarks.Extremities, T.c(1, :)) * T.T'; +scapulaGroove = bsxfun(@minus, out.Landmarks.Groove, T.c(1, :)) * T.T'; +scapulaPillar = bsxfun(@minus, out.Landmarks.Pillar, T.c(1, :)) * T.T'; + +% *********************************************************** +% Notes: landmarks.extremities (From Elham T.) +% (1) Bottom of the scapula --> angulusInferior +% (2) Back of the scapula --> trigonumSpinae +% (3) Tip of acromion --> processusCoracoideus ? +% (4) Top coracoid --> acromioClavicular ? +% (5) Back coracoid --> angulusAcromialis ? +% (6) Front scapula --> spinoGlenoidNotch +% ************************************************************ + +ScapulaLandmarks.angulusInferior = scapulaExtremities(1, :); % Angulus Inferior +ScapulaLandmarks.trigonumSpinae = scapulaExtremities(2, :); % Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine +ScapulaLandmarks.processusCoracoideus = scapulaExtremities(3, :); % Processus Coracoideus (most lateral point) +ScapulaLandmarks.acromioClavicular = scapulaExtremities(4, :); % Acromio-clavicular joint (most lateral part on acromion) +ScapulaLandmarks.angulusAcromialis = scapulaExtremities(5, :); % Angulus Acromialis (most laterodorsal point) +ScapulaLandmarks.spinoGlenoidNotch = scapulaExtremities(6, :); % Spino-glenoid notch +ScapulaLandmarks.groove = scapulaGroove; % 5 landmarks on the scapula (supraspinatus) groove +ScapulaLandmarks.pillar = scapulaPillar; % 5 landmarks on the pillar + +if isempty(ScapulaLandmarks.pillar) + ScapulaLandmarks = -1; +end + +end diff --git a/findScapulaLandmarksSCase.m b/findScapulaLandmarksSCase.m new file mode 100755 index 0000000..1d8c9ae --- /dev/null +++ b/findScapulaLandmarksSCase.m @@ -0,0 +1,105 @@ +function outputArg = findScapulaLandmarksSCase(varargin) +%FINDSCAPULALANDMARKS Find scapula landmarks of segmented SCase +% Detailed explanation goes here + +% Can only be run from lbovenus to avoid error +% Undefined function or variable 'smoothpatch_curvature_double' +% in function smoothpatch in function findScapulaLandmarks + +% The function can be run from (lbovenus) server as follows +% cd /home/shoulder/methods/matlab/segmentation/method +% matlab -nodisplay -nosplash -r "findScapulaLandmarksSCase;quit" + +% The log file can be checked with shell command +% cd /home/shoulder/methods/matlab/segmentation/method +% tail -f log/findScapulaLandmarksSCase.log + + +% Input: varargin +% Can be empty, 'N', 'P', or a SCase ('P315') + +% Output: +% log file findScapulaLandmarksSCase.log +% result file scapulaLandmarksAuto{L/R}.mat + +% Associated functions: +% listSCase: +% simplePlyRead: +% findScapulaLandmarks: + +% Author: Alexandre Terrier, EPFL-LBO +% Date: 2018-08-21 + +% TODO + +tic; +% Open log file +filename = 'findScapulaLandmarksSCase.log'; +logDir = 'log'; +filename = [logDir '/' filename]; +logFileId = fopen(filename, 'w'); +fprintf(logFileId, 'findScapulaLandmarksSCase log file'); +fprintf(logFileId, ['\n\nDate: ' datestr(datetime)]); + +% Add path to database functions folder +path2Add = '../../database'; +addpath(path2Add); +% Add path to segmenation functions folder +%path2Add = 'functions'; +%addpath(path2Add); +addpath(genpath('./functions')); + +% Get list of sCases from varargin +fprintf(logFileId, '\n\nList of SCase'); +if nargin == 0 % findScapulaLandmarksSCase called without argument + SCaseList = listSCase; +elseif nargin == 1 + inputArg = varargin{1}; + SCaseList = listSCase(inputArg); +else + error('inputArg can be empty, N, P, or a SCaseID') +end + +% Loop on list of SCase +nSCase = length(SCaseList); % Number of SCase +for iSCase = 1:nSCase + SCaseID = SCaseList(iSCase).id; + fprintf(logFileId, ['\n\nSCaseID: ' SCaseID]); + currentSCaseDir = SCaseList(iSCase).dir; + matlabDir = [currentSCaseDir '/matlab']; + % Find ply files of scapula surface in matlab dir (there might be 0, 1 or 2) + filenameList = dir([matlabDir '/scapulaSurfaceAuto*.ply']); + if isempty(filenameList) + fprintf(logFileId, '\n No ply file'); + continue + end + for iFilename = 1:length(filenameList) % Loop on segmented scapula ply files (L/R) + filename = filenameList(iFilename); + filename = [filename.folder '/' filename.name]; + try + FV = simplePlyRead(filename); % read the input ply file + catch + fprintf(logFileId, '\n Error reading ply file'); + continue; + end + [ScapulaLandmarks, message] = findScapulaLandmarks(FV); % Find scapula landmarks + if isstruct(ScapulaLandmarks) + % Save scapula landmarks + if message + fprintf(logFileId, ['\n ' message]); % We might add here R/L/- + end + filename = strrep(filename, 'Surface', 'Landmarks'); % Replace Surface by Landmarks + filename = strrep(filename, '.ply', '.mat'); % Replace .ply by .mat + save(filename, 'ScapulaLandmarks'); + fprintf(logFileId, ['\n Save ' filename]); + else + fprintf(logFileId, '\n Error finding landmarks'); + end + end +end % End of SCase loop + +fprintf(logFileId, ['\n\n End of findScapulaLandmarksSCase after ' num2str(toc) ' seconds']); + fclose(logFileId); % Close log file + outputArg = 1; +end + diff --git a/segmentScapula.m b/segmentScapula.m new file mode 100755 index 0000000..7c761b1 --- /dev/null +++ b/segmentScapula.m @@ -0,0 +1,453 @@ +function outputArg = scapulaSegment(list) +%SCAPULASEGMENT Semgment scapula in /shoulder/data/{P\N} +% ses statistical shape methods to segments the scapula, and identify +% scapula landmarks and glenoid surface. + +% It can be executed from shell with the command +% matlab -nojvm -nodisplay -nosplash -r "scapulaSegment;quit" > scapulaSegment.log +% For n pathological cases it takes about n minutes + +% Input: either a list of nothing for all + +% Output: + + +% Author: Alexandre Terrier +% Date: 2018-08-24 + +% Comment 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 (scapulaSegment.m) +% Still to be done: +% - Change output directory and file names +% - Obtain and save directly (not afterwards) the scapula landmarks and +% glenoid surface +% - Consistency with ShoulderCase class, fornaming and saving. +% - Use of log file instaed terminal + +% Comment by BN: +% Note: I use 'try' and 'catch' for many of the procedures. This is because +% they use routines that are third party and we do not know why they crash +% or where we might come across problems. In every catch segment, I try to +% store the point of failure so it can be reported to Elham in Bern. + +% TODO: +% Add log file +% Move to database folder +% Perform anatomical data save in this function instead of bharath2remove/extractAllGlenoid.m +% Change main loop to consider N & P, might use listCAses +% Add argument to force segement existing segmentations + + +% Check input arguments +if(nargin < 1) + isList = 0; +else + isList = 1; + listPat = list; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Initialisation + +% Parameters of the segmentation. These parameters are used and you can use +% them as they are: opt.isLeft = true; +opt.isLeft = false; +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]; + + +iSegmented = 1; % counter for segmented scapulas +location = 'P'; % Location - either Normal or Pathological +iSkipped = 1; % counter for skipped scapulas + +% Set paths +segmentationCodePath = '../../segmentation/'; % Segmentation code +dataPath = '../../../../data/'; % Database +modelPath = [segmentationCodePath 'Model/']; % Statistical Shape Models (right & left scapula) + +% Add required functions to matlab path +path = [segmentationCodePath 'method/functions']; +addpath(genpath(path)); + +% 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) +display('Loading statistical shape models...') +tic +disp(' right scapula...'); +modelLeft = loadModel(modelPath, true); % isLeft = true; +disp(' left scapula...'); +modelRight = loadModel(modelPath, false); % isLeft = false; +toc + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Segmentation + +disp('Start segmenting cases...'); +% Loop through data folder, and 2 levels of directories +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; + % disp(folderIn); + disp(['sCase.id: ' iCT]); + % AT comment: P255 is non conventional and should be skipped + if strcmp(iCT,'P255') + disp(' Skip P255 until conform'); + continue; + end + + % Comment AT: We should somehow check (below) that this + % case has been already tried and could not be analyze for + % some reason, and thus avoid loosing 20 minutes in + % re-trying. + + filePly = strcat(folderIn,'matlab/scapulaSurface.ply'); + + % AT: test below should be replaced by listSCase output + % Check if a list of sCase is prodived + if(isList) + % in case of a list being provided, skip files that are + % not in the list + if(~any(strcmp(list,iCT))) + continue; + end + elseif(exist(filePly)) % AT: Should be replaced by ~exist + % otherwise, check if segmentation has already + % been done, and exist loop in this case + % We might force to redo with an addition argument of + % the function. + fprintf('Segmentation already done. \n'); + continue; + end + + % AT: remove the mdir of auto + %{ + % folder for output to be stored + folderOut = [folderIn 'matlab/auto']; + try + if ~exist(folderOut, 'dir') % create folder if it does not exist + mkdir(folderOut); + end + catch + fprintf('error creating the folder \n'); + end + if(folderIn(end)~=filesep) + folderIn = [folderIn filesep]; + end + + if(folderOut(end)~=filesep) + folderOut = [folderOut filesep]; + end + %} + + %% Loading of dicom files + + % 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 + +% disp(['folderIn: ' folderIn]); +% disp(['folderOut: ' folderOut]); +% disp(bones_in); + + % 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) + disp('Loading images...'); + tic; + + 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) + 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('Index exceeds matrix dimensions.'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Index exceeds matrix dimensions. Line 147'; + iSkipped = iSkipped+1; + continue; + end + break; + end + end + end + + fname = boneName; + ind = strfind(fname, '.'); + if(isempty(ind)) + continue; + 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 + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + end + case 'mhd' + try + [imgMHD, info]=read_mhd(boneName); + catch + display([bones_in(b).name ': Problem in reading the file. Check the file and try again later!']); + 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 + display ([ 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)) ')']); + continue; + end + clear imgMHD; + case 'dcm' + dicomDir = strsplit(boneName, filesep); + dicomDir = strjoin(dicomDir(1, 1:end-1), filesep); + [Img, imgInfo] = readDicomStack(dicomDir); + if (isempty(Img)) + display ([ bones_in(b).name ': Could not read the dicom images.']); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Could not read the dicom images.'; + iSkipped = iSkipped+1; + continue; + end + otherwise + continue; + end + fname = bones_in(b).name; + ind = strfind(fname, '.'); + if(isempty(ind)) + ind = length(fname)+1; + else + ind = ind(end); + end + toc; + % End of dicom loading + + %% Identify scapula side (takes about 30 seconds) + disp('Identifying scapula side...'); + try + tic; + [Left, Right] = IdentifyRightAndLeftScapula(Img, imgInfo, modelLeft, modelRight); % Img and imgInfo is similar to autoSegment function. + toc; + catch + fprintf('Scapula identification failed.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Scapula identification failed.'; + iSkipped = iSkipped+1; + continue; + end + + %% Segment scapula (takes about 5-15 minutes) + disp('Segmenting CT...'); + tic; + + if(isempty(Right) && ~isempty(Left)) + disp('Image contains only the left scapula.'); + opt.isLeft = true; + try + finalSurface = autoSegment(Img, imgInfo, modelLeft, opt); + catch + fprintf('Inconsistent data skipped.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Inconsistent Data.'; + iSkipped = iSkipped+1; + continue; + end + end + if(isempty(Left) && ~isempty(Right)) + disp('Image contains only the right scapula.'); + opt.isLeft = false; + try + finalSurface = autoSegment(Img, imgInfo, modelRight, opt); + catch + fprintf('Inconsistent data skipped.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Inconsistent Data.'; + iSkipped = iSkipped+1; + continue; + end + end + if (~isempty(Right) && ~isempty(Left)) + disp('Image contains both left and right scapulae.'); + % IMPORTANT: In case both scapulae are present, right now I choose only the right one. + + % The below lines were used when there were issues in segmenting scapulae with both sides in it. + % Now they can be commented out safely. + + % 4 lines blow commented by AT because of issue + % when both sides are present in CT + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Image contains both left and right scapulae.'; + iSkipped = iSkipped+1; + disp(' skipping'); + continue; + + opt.isLeft = false; + try + finalSurface_Right = autoSegment(Right.img, Right.imgInfo, modelRight, opt); + finalSurface = finalSurface_Right; + catch + fprintf('Inconsistent data skipped.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Inconsistent Data.'; + iSkipped = iSkipped+1; + continue; + end + + % If you want to segment only the left scapula, uncomment the following lines and comment the + % lines above, starting from 'opt.isLeft = False' + +% opt.isLeft = true; +% try +% finalSurface_Left = autoSegment(Left.img, Left.imgInfo, modelLeft, opt); +% finalSurface = finalSurface_Left; +% catch +% fprintf('Inconsistent data skipped.\n'); +% skippedRecord{iSkipped,1} = iCT; +% skippedRecord{iSkipped,2} = 'Inconsistent Data.'; +% iSkipped = iSkipped+1; +% continue; +% end + end + + if (isempty(Right) && isempty(Left)) + disp('Could not recognize, which scapulae it is.'); + end + toc; + % End of segmentation + + %% Write surface file + disp('Writing surface file...'); + try + simplePlyWrite(finalSurface, [folderOut 'scapulaSurfaceAuto.ply']); + catch + fprintf('File writing failed.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'File writing failed.'; + iSkipped = iSkipped+1; + continue; + end + + % Cemment by AT: This should be adapted to the + % ShoulderCase class We should add below the + % identification of the scapula landmarks and the + % glenoid surface. The scapula landmarks are currently + % embeded in function doMeasurements (by ET). Glenoid + % surface is currently obtained as a second step with + % function extractAllGlenoid and extractGlenoidSurface + % (by BN). + + +%{ +% AT: removed the measurement part + % Conduct measurements using Elham's code + disp('Measuring anatomy (scapula landmarks and glenoid surface)...'); + tic; + try + [measurements] = doMeasurements(finalSurface); % Conduct measurements using Elham's code + % From measurements, we need measurements.Landmarks + % Below, we obtain all the measurements using our + % own code. We need to add measurements.Landmarks + % to Result. +% [Result] = scapula_calculation_ply_glenoidManual(subject, directory, location, 0); % Conduct measurements using our code +% Result.landmarks_Auto = measurements.Landmarks; + save( [folderOut 'anatData.mat'],'measurements'); + catch + fprintf('Unable to do measurements.\n'); + skippedRecord{iSkipped,1} = iCT; + skippedRecord{iSkipped,2} = 'Unable to do measurements.'; + iSkipped = iSkipped+1; + continue; + end + + toc; + %} + +% disp(['Processing of sample ' fname(1:ind-1) ' is finished.']) + + segmentedRecord{iSegmented,1} = iCT; + iSegmented = iSegmented + 1; + end + +%{ +% AT: removed this part + % Finally, I save all the skipped and segmented data + % regularly. I do this inside the loop because sometimes, + % when the server connection is lost, I lose all the + % information. + if(iSkipped > 1) + save( 'skippedData-Normal3.mat','skippedRecord'); + end + if(iSegmented > 1) + save( 'segmentedData-Normal3.mat','segmentedRecord'); + end + %} + + + end + end + end +end +outputArg = 1; % Might be extended +end