Homec4science
Diffusion segmentation (master) Browse

segmentation/findGlenoidSurface.m
8222343c77d4master

Last Change
Raw File

findGlenoidSurface.m
View Options

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
c4science · Help