diff --git a/functions/readDicomStack.m b/functions/readDicomStack.m
old mode 100644
new mode 100755
index 9cd19b2..23c82f8
--- a/functions/readDicomStack.m
+++ b/functions/readDicomStack.m
@@ -1,346 +1,353 @@
function [img, imgInfo] = readDicomStack(dicomDir, opt)
%
+dicomdict('set','dicom.dic');
% kernelTypes = {'bone', 'lung', 'none'};
if(~isempty(dicomDir) && dicomDir(end) ~= filesep)
dicomDir = [dicomDir filesep];
end
if (nargin == 0)
dicomDir = './';
opt = 'max';
else
if (nargin == 1)
opt = 'max';
end
end
if class(dicomDir)=='char'
filesList = dir([dicomDir '*.dcm']);
if(isempty(filesList))
filesList = dir(dicomDir);
filesList = filesList(3:end);
end
tmp{size(filesList, 1)} = '';
for ii = 1:size(filesList, 1)
tmp{ii} = [dicomDir filesList(ii).name];
end
filesList = tmp;
else
filesList = dicomDir;
end
Headers{size(filesList, 2)} = [];
sliceLoc = zeros(size(filesList, 1), 1);
kernelType = zeros(size(sliceLoc));
for ii = 1:size(filesList, 2)
try
Headers{ii} = dicominfo([filesList{ii}]);
catch
sliceLoc(ii) = inf;
continue;
end
try
sliceLoc(ii) = Headers{ii}.ImagePositionPatient(3);
catch
sliceLoc(ii) = -inf;
end
if (isfield(Headers{ii}, 'ConvolutionKernel'))
switch lower(Headers{ii}.ConvolutionKernel)
case 'bone'
kernelType(ii) = 1;
case 'lung'
kernelType(ii) = 2;
case 'standard'
kernelType(ii) = 2;
case 'soft'
kernelType(ii) = 2;
otherwise
kernelType(ii) = 3;
end
else
kernelType(ii) = 3;
end
end
[kernelType, ind] = sort(kernelType);
sliceLoc = sliceLoc(ind);
Headers = Headers(:, ind);
numStacks = 0;
ii = 1;
while (ii<size(Headers, 2))
jj = ii;
while(jj<=size(Headers, 2) && kernelType(jj) == kernelType(ii)), jj = jj+1; end
tmp = sliceLoc(ii:jj-1);
[sliceLoc(ii:jj-1), ind] = sort(tmp);
tmp = Headers(:, ii:jj-1);
Headers(:, ii:jj-1) = tmp(:, ind);
ii = jj;
numStacks = numStacks+1;
end
if (numStacks == 0)
imgInfo = [];
img = [];
return;
end
+
currentSlice = 1;
img{numStacks} = [];
for ii = 1:numStacks
currentKernel = kernelType(currentSlice);
while(sliceLoc(currentSlice)==-inf), currentSlice = currentSlice+1; end
xMin = 0; xMax = 0; xSize = 0;
yMin = 0; yMax = 0; ySize = 0;
zMin = 0; zMax = 0; zSize = 0;
PixelSpacing = [0, 0];
SliceThickness = 0;
if (isfield(Headers{currentSlice}, 'PixelSpacing'))
PixelSpacing = Headers{currentSlice}.PixelSpacing;
end
% if (isfield(Headers{currentSlice}, 'SliceThickness'))
% SliceThickness = Headers{currentSlice}.SliceThickness;
% end
if (isfield(Headers{currentSlice}, 'ImagePositionPatient'))
xMin = Headers{currentSlice}.ImagePositionPatient(1);
yMin = Headers{currentSlice}.ImagePositionPatient(2);
zMin = Headers{currentSlice}.ImagePositionPatient(3);
end
if (isfield(Headers{currentSlice}, 'Columns'))
xSize = Headers{currentSlice}.Columns;
xMax = xMin + double(xSize - 1) * PixelSpacing(2);
end
if (isfield(Headers{currentSlice}, 'Rows'))
ySize = Headers{currentSlice}.Rows;
yMax = yMin + double(ySize - 1) * PixelSpacing(1);
end
switch currentKernel
case 1
imgInfo(ii).ConvolutionKernel = 'Bone';
case 2
imgInfo(ii).ConvolutionKernel = 'Lung';
otherwise
imgInfo(ii).ConvolutionKernel = 'None';
end
while(currentSlice <= numel(kernelType) && currentKernel == kernelType(currentSlice))
zSize = zSize+1;
try
tmp = dicomread(Headers{currentSlice});
if (isfield(Headers{currentSlice}, 'RescaleSlope') && Headers{currentSlice}.RescaleSlope~=0)
Slope = Headers{currentSlice}.RescaleSlope;
else
Slope = 1;
end
if (isfield(Headers{currentSlice}, 'RescaleIntercept'))
Intercept = Headers{currentSlice}.RescaleIntercept;
else
Intercept = 0;
end
catch
currentSlice = currentSlice + 1;
continue;
end
img{ii}(:, :, zSize) = tmp*Slope + Intercept;
if(SliceThickness == 0 && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
SliceThickness = sliceLoc(currentSlice+1) - sliceLoc(currentSlice);
else
newLoc = zMin+(zSize-1)*SliceThickness;
if(SliceThickness ~= 0 && newLoc ~= sliceLoc(currentSlice) && currentSlice<numel(kernelType) && kernelType(currentSlice+1) == currentKernel)
sliceLoc(currentSlice) = newLoc;
end
end
zMax = sliceLoc(currentSlice);
currentSlice = currentSlice + 1;
end
imgInfo(ii).BoundingBox.xMin = xMin; imgInfo(ii).BoundingBox.xMax = xMax; imgInfo(ii).xSize = xSize;
imgInfo(ii).BoundingBox.yMin = yMin; imgInfo(ii).BoundingBox.yMax = yMax; imgInfo(ii).ySize = ySize;
imgInfo(ii).BoundingBox.zMin = zMin; imgInfo(ii).BoundingBox.zMax = zMax; imgInfo(ii).zSize = zSize;
imgInfo(ii).SliceThickness = (zMax -zMin)/(zSize-1);
imgInfo(ii).PixelSpacing = [(xMax - xMin)/double(xSize-1) (yMax - yMin)/double(ySize-1)];
end
ind = [];
+display(['TEST : length(img) = ' num2str(length(img))]);
+display(['TEST : numStacks = ' num2str(numStacks)]);
+display(['TEST : opt = ' num2str(opt)]);
for ii = 1:numStacks
if(isempty(img{ii}))
ind = [ind, ii];
end
end
+
img(:, ind) = [];
imgInfo(:, ind) = [];
switch opt
case 'all'
case 'first'
img = img{1};
imgInfo = imgInfo(1);
case 'last'
img = img{end};
imgInfo = imgInfo(end);
otherwise
SZs = zeros(size(imgInfo, 2), 1);
for tmp = 1:size(imgInfo, 2)
SZs(tmp) = imgInfo(tmp).zSize;
end
tmp = max(SZs);
ind = find(SZs == tmp);
if(numel(ind)>1)
for ii=1:numel(ind)
if(imgInfo(ind(ii)).ConvolutionKernel == 'Lung')
break;
end
end
if (ii>numel(ind))
ind = ind(1);
else
ind = ind(ii);
end
end
- img = img{ind};
+ display(['TEST : length(img) = ' num2str(length(img))]);
+ display(['TEST : ind = ' num2str(ind)]);
imgInfo = imgInfo(ind);
end
% figure
% for ii = 1:size(img{1}, 3)
% imshow(uint8(img{1}(:, :, ii)/4095*255));
% pause(0.1)
% end
% % % stackNo = 1;
% % %
% % %
% % %
% % % info = dicominfo([dicomDir filesList(1).name]);
% % % tmp = dicomread(info);
% % % imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
% % % imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
% % % imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
% % % imgInfo(stackNo).PixelSpacing = info.PixelSpacing;
% % % imgInfo(stackNo).SliceThickness = 0;
% % % imgInfo(stackNo).xSize = info.Columns;
% % % imgInfo(stackNo).ySize = info.Rows;
% % % imgInfo(stackNo).zSize = 1;
% % % imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + double(imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
% % % imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + double(imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
% % % imgInfo(stackNo).BoundingBox.zMax = double(imgInfo(stackNo).BoundingBox.zMin);
% % % img{stackNo} = tmp;
% % % infoPre = info;
% % %
% % % numSlicesInStack = [];
% % % for ii = 2:size(filesList, 1)
% % % try
% % % info = dicominfo([dicomDir filesList(ii).name]);
% % % catch
% % % continue;
% % % end
% % % if ( ~isequal(info.Columns, infoPre.Columns) || ~isequal(info.Rows, infoPre.Rows) || ~isequal(info.PatientOrientation, infoPre.PatientOrientation))
% % % img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
% % % numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
% % % stackNo = stackNo + 1;
% % % try
% % % imgInfo(stackNo).BoundingBox.xMin = info.ImagePositionPatient(1);
% % % catch
% % % imgInfo(stackNo).BoundingBox.xMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).BoundingBox.yMin = info.ImagePositionPatient(2);
% % % catch
% % % imgInfo(stackNo).BoundingBox.yMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).BoundingBox.zMin = info.ImagePositionPatient(3);
% % % catch
% % % imgInfo(stackNo).BoundingBox.zMin = 0;
% % % end
% % % try
% % % imgInfo(stackNo).PixelSpacing = info.PixelSpacing;
% % % catch
% % % imgInfo(stackNo).PixelSpacing = 1;
% % % end
% % %
% % % imgInfo(stackNo).SliceThickness = 0;
% % % try
% % % imgInfo(stackNo).xSize = info.Columns;
% % % imgInfo(stackNo).ySize = info.Rows;
% % % catch
% % % tmp = dicomread(info);
% % % imgInfo(stackNo).xSize = size(tmp, 2);
% % % imgInfo(stackNo).ySize = size(tmp, 1);
% % % end
% % %
% % %
% % % imgInfo(stackNo).zSize = 1;
% % % imgInfo(stackNo).BoundingBox.xMax = imgInfo(stackNo).BoundingBox.xMin + (imgInfo(stackNo).xSize -1) * imgInfo(stackNo).PixelSpacing(1);
% % % imgInfo(stackNo).BoundingBox.yMax = imgInfo(stackNo).BoundingBox.yMin + (imgInfo(stackNo).ySize -1) * imgInfo(stackNo).PixelSpacing(2);
% % % imgInfo(stackNo).BoundingBox.zMax = imgInfo(stackNo).BoundingBox.zMin;
% % % infoPre = info;
% % % else
% % % try
% % % imgInfo(stackNo).BoundingBox.zMax = info.ImagePositionPatient(3);
% % % catch
% % % imgInfo(stackNo).BoundingBox.zMax = imgInfo(stackNo).BoundingBox.zMax + 1;
% % % end
% % % imgInfo(stackNo).zSize = imgInfo(stackNo).zSize + 1;
% % % imgInfo(stackNo).SliceThickness = (imgInfo(stackNo).BoundingBox.zMax - imgInfo(stackNo).BoundingBox.zMin)/imgInfo(stackNo).zSize;
% % % end
% % % end
% % % img{stackNo} = eval([class(tmp) '(zeros(imgInfo(stackNo).xSize, imgInfo(stackNo).ySize, imgInfo(stackNo).zSize))']);
% % % numSlicesInStack(stackNo) = imgInfo(stackNo).zSize;
% % %
% % % switch opt
% % % case 'all'
% % % f = 1;
% % % for ii = 1:length(numSlicesInStack)
% % % for jj = 1:numSlicesInStack(ii)
% % % info = dicominfo([dicomDir filesList(f).name]);
% % % tmp = dicomread(info);
% % % img{ii}(:, :, jj) = tmp;
% % % f = f + 1;
% % % end
% % % end
% % % case 'first'
% % % for ii = 1:numSlicesInStack(1)
% % % info = dicominfo([dicomDir filesList(ii).name]);
% % % tmp = dicomread(info);
% % % img{1}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(1);
% % % img = img{1};
% % % case 'last'
% % % if numel(numSlicesInStack) == 1
% % % fileInd = 0;
% % % else
% % % fileInd = sum(numSlicesInStack(end-1));
% % % end
% % % for ii = 1:numSlicesInStack(end)
% % % info = dicominfo([dicomDir filesList(ii+fileInd).name]);
% % % tmp = dicomread(info);
% % % img{end}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(end);
% % % img = img{end};
% % % otherwise
% % % [~, ind] = max(numSlicesInStack);
% % % if ind == 1
% % % fileInd = 0;
% % % else
% % % fileInd = sum(numSlicesInStack(1:ind-1));
% % % end
% % % for ii = 1:numSlicesInStack(ind)
% % % info = dicominfo([dicomDir filesList(ii+fileInd).name]);
% % % tmp = dicomread(info);
% % % img{ind}(:, :, ii) = tmp;
% % % end
% % % imgInfo = imgInfo(ind);
% % % img = img{ind};
% % % end
diff --git a/segmentScapula.m b/segmentScapula.m
index e34e1ea..94dddbf 100755
--- a/segmentScapula.m
+++ b/segmentScapula.m
@@ -1,1256 +1,1258 @@
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"
-% WARNING on lbovenus, be sure to run the matlab command with a bash shell
+% 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: varargin
% Can be empty, 'N', 'P', or a SCase ('P315')
% 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
% 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 = '../../../../data/'; % Database
databaseFunctions = '../../database'; % Database functions
% Add database functions to matlab path
addpath(databaseFunctions);
% 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)]);
% 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
% Parameters to redo segmentaions
retryError = 0; % InputArg to add to retry SCase with error
retryAll = 1; % InputArg to add to retry all SCase in SCaseList
% 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;
%% Start loop over sCases in sCaseList
nSCase = length(SCaseList); % Number of sCases
for iSCaseID = 1:nSCase
SCaseID = SCaseList(iSCaseID).id; %
percentProgress = num2str(iSCaseID/nSCase*100, '%3.1f');
fprintf(logFileId, ['\n\nSCaseID: ' SCaseID ' (' percentProgress '%%)']);
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
if (retryAll || (errorFile && retryError) || (~scapulaSurfaceFile && ~errorFile))
% Segment if
% 1) retryAll is true
% 2) Error file and retryError are true
% 3) No scapulaSurface.ply file and no error file
% 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;
continue % 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' scapulaSide 'Err.txt'];
message = 'Scapula side identification error';
fileID = fopen(filename, 'w');
fprintf(fileID, message);
fclose(fileID);
continue; % 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 scapule 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 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;
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 = [matlabDir '/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 = [matlabDir '/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 = [matlabDir '/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 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');
continue;
end
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');
continue;
end
end % End of save surface file(s)
else
if errorFile
fprintf(logFileId, ': Already tried segmentation with error');
else
fprintf(logFileId, ': Already segmented');
end
end % End of SCase to segment
end % End of SCaseList loop
% 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<ST2NM>
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
%}