Homec4science
Diffusion segmentation (master) Browse

segmentation/doMeasurements.m
8222343c77d4master

Last Change
Raw File

doMeasurements.m
View Options

function out = doMeasurements(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
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
if(strcmp(options.side, 'right'))
reflect = true;
else
reflect = false;
end
%% The important vertices needed for measurements
[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
S.FV = fv;
% [~, 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)
display('\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];
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)
display('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
S.GlenoidOrientation = ...
computeGlenoidOrientation(...
S.FV,...
S.indiceVertices.GlenoidSurface,...
S.Landmarks);
out = rmfield(S, 'FV');
out = rmfield(out, 'indiceVertices');
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';
end
%% function to 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 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
c4science · Help