more comments

STAPLE/algorithms/Kai2014_pelvis.m

@@ -2,15 +2,23 @@
 % reference system of the pelvis described in the following publication: 
 % Kai, Shin, et al. Journal of biomechanics 47.5 (2014):1229-1233. 
 % https://doi.org/10.1016/j.jbiomech.2013.12.013
-% The algorithm creates a reference system from the centre of the femoral 
-% head and the centres of the femoral condyles, after identifying these
-% feature through slicing the bone geometry with planes perpendicular to
-% the longitudinal axis and anterior-posterior axis respectively.
+% The algorithm creates a reference system based on principal axes of
+% inertia and geometrical considerations as described in the original
+% publication. Please note that:
+% 1. the original implementation includes the sacrum, here excluded because
+% difficult to segment from MRI scans.
+% 2. the original publication defines the pelvis reference system
+% differently from the reccomendations of the International Society of
+% Biomechanics.
+% 3. The robustness of the Kai2014 algorithm with respect to different
+% global reference systems is ensured, in our experience, only through the
+% further checks implemented in the pelvis_guess_CS.m function developed
+% for STAPLE_pelvis.
 %
 %   [CS, JCS, pelvisBL] = Kai2014_pelvis(pelvisTri,...
-%                                      side_raw,...
-%                                      result_plots, ...
-%                                      debug_plots, in_mm)
+%                                        side_raw,...
+%                                        result_plots, ...
+%                                        debug_plots, in_mm)
 %
 % Inputs:
 %   pelvisTri - MATLAB triangulation object of the entire pelvic geometry.

STAPLE/algorithms/private/GIBOC_tibia_ProxArtSurf_it1.m

@@ -40,7 +40,7 @@
 end
 
 
-%% Refine and separete medial and lateral AS region
+%% Refine and separate medial and lateral AS region
 % Compute seed points to get a patch of AS on each condyle
 MedPtsInit = mean(ellipsePts) + 2/3*ELP1.b*ELP1.Yel';
 LatPtsInit = mean(ellipsePts) - 2/3*ELP1.b*ELP1.Yel';

STAPLE/algorithms/private/Kai2014_femur_fitSpheres2Condyles.m

@@ -43,7 +43,6 @@
 %  Copyright 2020 Luca Modenese & Jean-Baptiste Renault
 %-------------------------------------------------------------------------%
 
-
 function CS = Kai2014_femur_fitSpheres2Condyles(DistFemTri, CS, debug_plots, debug_prints)
 
 % main plots are in the main method function. 

STAPLE/algorithms/private/pelvis_guess_CS.m

@@ -1,44 +1,58 @@
-%-------------------------------------------------------------------------%
-%  Author:   Jean-Baptiste Renault. 
-%  Copyright 2020 Luca Modenese & Jean-Baptiste Renault
-%-------------------------------------------------------------------------%
-function [RotPseudoISB2Glob, LargestTriangle, BL] = pelvis_guess_CS(Pelvis, debug_plots)
-% Function to test putting back together a correct orientation of the
-% pelvis
+% PELVIS_GUESS_CS Run geometrical checks to correctly estimate the 
+% orientation of an initial pelvis reference system. The convex hull for
+% the pelvis triangulation is computed, and the normal of largest triangle,
+% which connects the ASIS and PUB, identifies the frontal direction. Then, 
+% for the pelvis moved to the reference system defined by its principal
+% axes of inertia, the points that have the largest span along a
+% non-prox-distal axis are identified as the crests of the iliac bones.
+% Using the centroid of the triangulation, a cranial axis is defined and
+% the reference system finalised after the recommendation of the
+% International Society of Biomechanics (ISB).
+%
 % Inputs :
-%           Pelvis : A triangulation of a complete pelvis
-%           debug_plots : A boolean to display plots useful for debugging
+%   pelvisTri - MATLAB triangulation object of the entire pelvic geometry.
+%
+%   debug_plots - enable plots used in debugging. Value: 1 or 0 (default).
 %
 % Output :
-%   BL              A structure containing the pelvis bone landmarks
-%   RotISB2Glob     A rotation matrix containing properly oriented initial
-%                   guess of the X, Y and Z axis of the pelvis ISB CS
-%   CenterVol       Volumetric center of the pelvis
+%   RotPseudoISB2Glob - rotation matrix containing properly oriented initial
+%        guess of the X, Y and Z axis of the pelvis. The axes are not
+%        defined as by ISB definitions, but they are pointing in the same
+%        directions, which is why is named "PseudoISB". This matrix 
+%        represent the body-fixed transformation from this reference system
+%        to ground.
+%
+%   LargestTriangle - MATLAB triangulation that identifies the largest
+%       triangle of the convex hull computed for the pelvis triangulation.
+%
+%   BL - MATLAB structure containing the bony landmarks identified 
+%       on the bone geometries based on the defined reference systems. Each
+%       field is named like a landmark and contain its 3D coordinates.
 %
 % -------------------------------------------------------------------------
-%                           General Idea
-% The largest cross section along the principal inertia axis is located at
-% the front of the pelvis. From that information it's easy to determine the
-% distal to proximal direction. Then the largest triangle of the pelvis
-% convex hull is the one connecting the LASIS RASIS and SYMP.
-% -------------------------------------------------------------------------
+%  Author:   Jean-Baptiste Renault and Luca Modenese. 
+%  Copyright 2020 Luca Modenese & Jean-Baptiste Renault
+%-------------------------------------------------------------------------%
+
+function [RotPseudoISB2Glob, LargestTriangle, BL] = pelvis_guess_CS(pelvisTri,...
+                                                                    debug_plots)
 
 % inputs checks
-if nargin < 2
-    debug_plots = 0;
-end
+if nargin < 2; debug_plots = 0; end
 
 % inertial axes
-[V_all, CenterVol, ~, D ] =  TriInertiaPpties(Pelvis);
+[V_all, CenterVol, ~, D ] =  TriInertiaPpties(pelvisTri);
 
 % smaller moment of inertia is normally the medio/lateral axis. It will be
 % updated anyway. It can be checked using D from TriInertiaPpties.m
 Z0 = V_all(:,1);
 
-%% Get convexHull
-K = convhull(Pelvis.Points);
+% compute convex hull
+K = convhull(pelvisTri.Points);
+
+% transform it in triangulation
 Kold2new(sort(unique(K(:)))) = 1:length(sort(unique(K(:))));
-Pts = Pelvis.Points(sort(unique(K(:))),:);
+Pts = pelvisTri.Points(sort(unique(K(:))),:);
 PelvisConvHull = triangulation(Kold2new(K), Pts);
 
 %% Get the Post-Ant direction by finding the largest triangle of the pelvis
@@ -55,6 +69,9 @@
                                 PelvisConvHull.ConnectivityList(I,:) , :);
 LargestTriangle = triangulation([1 2 3], LargestTriangle.Pts);
 
+% NOTE that we are working using a GIBOC reference system until where the
+% rotation matrix is assembled using ISB conventions(specified in comments)
+
 % vector pointing forward is X
 [~, ind_X] = max(abs(V_all'*LargestTriangle.faceNormal'));
 X0 = V_all(:,ind_X);
@@ -70,15 +87,15 @@
 
 if debug_plots == 1 
     figure()
-    quickPlotTriang(Pelvis)
+    quickPlotTriang(pelvisTri)
     plotArrow(X0, 1, CenterVol, 60, 1, 'r')
     plotArrow(Y0_temp, 1, CenterVol, 60, 1, 'g')
     plotArrow(Z0, 1, CenterVol, 60, 1, 'b')
     title('X0 shold be pointing anteriorly - no interest in other axes')
 end
 
 % transform the pelvis to the new set of inertial axes
-[ PelvisInertia, ~ , ~ ] = TriChangeCS(Pelvis, [X0, Y0_temp, Z0]', CenterVol);
+[ PelvisInertia, ~ , ~ ] = TriChangeCS(pelvisTri, [X0, Y0_temp, Z0]', CenterVol);
 
 if debug_plots == 1 
     figure()
@@ -89,15 +106,18 @@
     title('Check axes of inertia orientation')
 end
 
-% [~, ind_medLat] = max(abs(Y0_temp));
+% get points that could be on iliac crests
 [L1y, ind_P1y] = max(PelvisInertia.Points(:,2));
 [L2y, ind_P2y] = min(PelvisInertia.Points(:,2));
 spanY = abs(L1y)+abs(L2y);
 
+% get points that could be on iliac crests (remaning axis)
 [L1_z, ind_P1z] = max(PelvisInertia.Points(:,3));
 [L2_z, ind_P2z] = min(PelvisInertia.Points(:,3));
 spanZ = abs(L1_z)+abs(L2_z);
 
+% the largest span will identify the iliac crests points and discard other
+% directions
 if spanY>spanZ
     ind_P1 = ind_P1y;
     ind_P2 = ind_P2y;
@@ -116,20 +136,20 @@
 
 % these are the most external points in the iliac wings
 % these are iliac crest tubercles (ICT)
-P1 = Pelvis.Points(ind_P1,:);
-P2 = Pelvis.Points(ind_P2,:);
+P1 = pelvisTri.Points(ind_P1,:);
+P2 = pelvisTri.Points(ind_P2,:);
 P3 = (P1+P2)/2;% midpoint
 
 if debug_plots == 1 
     figure()
-    quickPlotTriang(Pelvis)
+    quickPlotTriang(pelvisTri)
     plotDot(P1, 'k', 7);
     plotDot(P2, 'k', 7);
     plotDot(P3, 'k', 7);
     title('Points should be external points in the iliac wings (glob ref syst)')
 end
 
-%upward vector (perpendicular to X0)
+% upward vector (perpendicular to X0)
 upw_ini = normalizeV(P3-CenterVol');
 upw = upw_ini-(upw_ini'*X0)*X0;
 
@@ -141,6 +161,7 @@
 if debug_plots == 1
     plotArrow(Z0, 1, CenterVol, 60, 1, 'b')
 end
+
 % Until now I have used GIBOC convention, now I build the ISB one!
 % X0 = X0_ISB, Z0 = Y_ISB
 RotPseudoISB2Glob = [X0,  Z0, cross(X0, Z0)];
@@ -149,23 +170,24 @@
 BL.ICT1 = P1;
 BL.ICT2 = P2;
 
-%% Debug Plots
+% debug Plots
 if debug_plots
     figure()
-    hold on
-    axis equal
+    hold on; axis equal
+    % define a ref system to plot
     temp.V = RotPseudoISB2Glob;
     temp.Origin = P3;
+    % plot axes of pelvis (GIBOC)
 %     pl3tVectors(CenterVol, X0, 125);
 %     pl3tVectors(CenterVol, Y0, 175);
 %     pl3tVectors(CenterVol, Z0, 250);
-    trisurf(Pelvis,'facealpha',0.5,'facecolor','b',...
-        'edgecolor','none');
-    trisurf(PelvisConvHull,'facealpha',0.2,'facecolor','c',...
-        'edgecolor',[.3 .3 .3], 'edgealpha', 0.2);
-    trisurf(LargestTriangle,'facealpha',0.8,'facecolor','r',...
-        'edgecolor','k');
+    % plot pelvis, convex hull and largest identified triangle
+    trisurf(pelvisTri,'facealpha',0.5,'facecolor','b', 'edgecolor','none');
+    trisurf(PelvisConvHull,'facealpha',0.2,'facecolor','c','edgecolor',[.3 .3 .3], 'edgealpha', 0.2);
+    trisurf(LargestTriangle,'facealpha',0.8,'facecolor','r','edgecolor','k');
+    % plot axes of pelvis (ISB)
     quickPlotRefSystem(temp)
+    % plot landmarks
     plotDot(P1, 'k', 7);
     plotDot(P2, 'k', 7);
     plotDot(P3, 'k', 7);