diff --git a/ShoulderCase/@Shoulder/Shoulder.m b/ShoulderCase/@Shoulder/Shoulder.m
index 9fad905..137fff2 100644
--- a/ShoulderCase/@Shoulder/Shoulder.m
+++ b/ShoulderCase/@Shoulder/Shoulder.m
@@ -1,221 +1,221 @@
 classdef (Abstract) Shoulder < handle
 
   properties
     side % R or L
     scapula
     humerus
     muscles
     CTscan
     comment
   end
 
 
 
   properties (Hidden = true)
     SCase
   end
 
 
 
   methods (Abstract, Access = protected)
     createScapula(obj);
   end
 
 
 
   methods
     function obj = Shoulder(SCase)
       obj.side    = '';
       obj.humerus = Humerus(obj);
       obj.muscles = Muscles(obj);
       obj.CTscan  = '';
       obj.comment = '';
 
       obj.SCase   = SCase;
 
       obj.createScapula;
     end
 
     function output = plot(obj)
       % Create a figure and plot inside:
       % the scapula surface, coordSys, groove points and landmarks;
       % the glenoid surface, sphere and centerLine;
       % the humerus sphere and centerLine.
 
       figure;
       hold on;
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       %% Scapula
       % Plot the segmented surface
       if not(isempty(obj.scapula.surface))
           points = obj.scapula.surface.points;
           faces  = obj.scapula.surface.faces;
           x = points(:,1);
           y = points(:,2);
           z = points(:,3);
           trisurf(faces,x,y,z,'Facecolor','yellow','FaceAlpha',0.5,'EdgeColor','none','FaceLighting','gouraud');
       end
 
       % Plot scapula groove
       if not(isempty(obj.scapula.groove))
         points = obj.scapula.groove;
         x = points(:,1);
         y = points(:,2);
         z = points(:,3);
         scatter3(x,y,z,100,'blue');
       end
 
       % Plot scapula markers
       points = [...
           obj.scapula.angulusInferior;
           obj.scapula.trigonumSpinae;
           obj.scapula.processusCoracoideus;
           obj.scapula.acromioClavicular;
           obj.scapula.angulusAcromialis;
           obj.scapula.spinoGlenoidNotch...
           ];
       x = points(:,1);
       y = points(:,2);
       z = points(:,3);
       scatter3(x,y,z,100,'blue','LineWidth',2);
 
 
-      % Plot scapular plane
-      plot(obj.scapula.plane,100);
-
-      % Plot scapular axis
-      pt1 = obj.scapula.coordSys.origin;
-      axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]);
-      pt2 = pt1 - obj.scapula.coordSys.ML * axisLength;
-      x = [pt1(1), pt2(1)];
-      y = [pt1(2), pt2(2)];
-      z = [pt1(3), pt2(3)];
-      plot3(x,y,z,'Color','blue','LineWidth', 5);
-
-      % Plot scapular coordinate system
-      axisLength = 50; % Length of the coord. sys. axes
-      pt0 = obj.scapula.coordSys.origin;
-      pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X
-      pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y
-      pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z
-      x = [pt0(1), pt1(1)];
-      y = [pt0(2), pt1(2)];
-      z = [pt0(3), pt1(3)];
-      plot3(x,y,z,'Color','red','LineWidth', 5);
-      x = [pt0(1), pt2(1)];
-      y = [pt0(2), pt2(2)];
-      z = [pt0(3), pt2(3)];
-      plot3(x,y,z,'Color','green','LineWidth', 5);
-      x = [pt0(1), pt3(1)];
-      y = [pt0(2), pt3(2)];
-      z = [pt0(3), pt3(3)];
-      plot3(x,y,z,'Color','blue','LineWidth', 5);
-
-      % plot scapular axis on glenoid center
-      axisLength = 50;
-      pt1 = obj.scapula.glenoid.center;
-      pt2 = pt1 + obj.scapula.coordSys.ML*axisLength;
-      x = [pt1(1), pt2(1)];
-      y = [pt1(2), pt2(2)];
-      z = [pt1(3), pt2(3)];
-      plot3(x,y,z,'Color','magenta','LineWidth', 5);
+%       % Plot scapular plane
+%       %plot(obj.scapula.plane,100);
+% 
+%       % Plot scapular axis
+%       pt1 = obj.scapula.coordSys.origin;
+%       axisLength = 10 + pdist([pt1 ; obj.scapula.trigonumSpinae]);
+%       pt2 = pt1 - obj.scapula.coordSys.ML * axisLength;
+%       x = [pt1(1), pt2(1)];
+%       y = [pt1(2), pt2(2)];
+%       z = [pt1(3), pt2(3)];
+%       plot3(x,y,z,'Color','blue','LineWidth', 5);
+% 
+%       % Plot scapular coordinate system
+%       axisLength = 50; % Length of the coord. sys. axes
+%       pt0 = obj.scapula.coordSys.origin;
+%       pt1 = pt0 + obj.scapula.coordSys.PA*axisLength; % X
+%       pt2 = pt0 + obj.scapula.coordSys.IS*axisLength; % Y
+%       pt3 = pt0 + obj.scapula.coordSys.ML*axisLength; % Z
+%       x = [pt0(1), pt1(1)];
+%       y = [pt0(2), pt1(2)];
+%       z = [pt0(3), pt1(3)];
+%       plot3(x,y,z,'Color','red','LineWidth', 5);
+%       x = [pt0(1), pt2(1)];
+%       y = [pt0(2), pt2(2)];
+%       z = [pt0(3), pt2(3)];
+%       plot3(x,y,z,'Color','green','LineWidth', 5);
+%       x = [pt0(1), pt3(1)];
+%       y = [pt0(2), pt3(2)];
+%       z = [pt0(3), pt3(3)];
+%       plot3(x,y,z,'Color','blue','LineWidth', 5);
+% 
+%       % plot scapular axis on glenoid center
+%       axisLength = 50;
+%       pt1 = obj.scapula.glenoid.center;
+%       pt2 = pt1 + obj.scapula.coordSys.ML*axisLength;
+%       x = [pt1(1), pt2(1)];
+%       y = [pt1(2), pt2(2)];
+%       z = [pt1(3), pt2(3)];
+%       plot3(x,y,z,'Color','magenta','LineWidth', 5);
+% 
+% 
+%       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%       % Glenoid
+% 
+%       % Plot the glenoid surface
+%       % We might move them in the direction of th glenoid center
+%       % instead of the ML axis
+% 
+%       points = obj.scapula.glenoid.surface.points +...
+%           obj.scapula.coordSys.ML * 0.2;
+%       faces  = obj.scapula.glenoid.surface.faces;
+%       x = points(:,1);
+%       y = points(:,2);
+%       z = points(:,3);
+%       trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none');
+% 
+% 
+%       % plot glenoid centerline
+%       pt1 = obj.scapula.glenoid.center;
+%       pt2 = pt1 + obj.scapula.glenoid.centerLine;
+%       x = [pt1(1), pt2(1)];
+%       y = [pt1(2), pt2(2)];
+%       z = [pt1(3), pt2(3)];
+%       plot3(x,y,z,'Color','magenta','LineWidth', 5);
+% 
+%       % We might add spherical cap
+
+%       % plot glenoid sphere
+%       n = 20; % number of lines
+%       [X,Y,Z] = sphere(n);
+%       X = X*obj.scapula.glenoid.radius;
+%       Y = Y*obj.scapula.glenoid.radius;
+%       Z = Z*obj.scapula.glenoid.radius;
+%       Xc = obj.scapula.glenoid.center(1);
+%       Yc = obj.scapula.glenoid.center(2);
+%       Zc = obj.scapula.glenoid.center(3);
+%       mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      % Glenoid
-
-      % Plot the glenoid surface
-      % We might move them in the direction of th glenoid center
-      % instead of the ML axis
-
-      points = obj.scapula.glenoid.surface.points +...
-          obj.scapula.coordSys.ML * 0.2;
-      faces  = obj.scapula.glenoid.surface.faces;
-      x = points(:,1);
-      y = points(:,2);
-      z = points(:,3);
-      trisurf(faces,x,y,z,'Facecolor','none','FaceAlpha',1,'EdgeColor','magenta','FaceLighting','none');
-
-
-      % plot glenoid centerline
-      pt1 = obj.scapula.glenoid.center;
-      pt2 = pt1 + obj.scapula.glenoid.centerLine;
-      x = [pt1(1), pt2(1)];
-      y = [pt1(2), pt2(2)];
-      z = [pt1(3), pt2(3)];
-      plot3(x,y,z,'Color','magenta','LineWidth', 5);
-
-      % We might add spherical cap
-
-      % plot glenoid sphere
-      n = 20; % number of lines
-      [X,Y,Z] = sphere(n);
-      X = X*obj.scapula.glenoid.radius;
-      Y = Y*obj.scapula.glenoid.radius;
-      Z = Z*obj.scapula.glenoid.radius;
-      Xc = obj.scapula.glenoid.center(1);
-      Yc = obj.scapula.glenoid.center(2);
-      Zc = obj.scapula.glenoid.center(3);
-      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
-
-
-      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      % Humerus
-
-      % plot humeral head sphere
-      n = 20; % number of lines
-      [X,Y,Z] = sphere(n);
-      X = X*obj.humerus.radius;
-      Y = Y*obj.humerus.radius;
-      Z = Z*obj.humerus.radius;
-      Xc = obj.humerus.center(1);
-      Yc = obj.humerus.center(2);
-      Zc = obj.humerus.center(3);
-      mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
-
-
-      % plot humeral head centerline
-      pt1 = obj.scapula.glenoid.center;
-      pt2 = obj.humerus.center;
-      x = [pt1(1), pt2(1)];
-      y = [pt1(2), pt2(2)];
-      z = [pt1(3), pt2(3)];
-      plot3(x,y,z,'Color','magenta','LineWidth', 5);
+%       % Humerus
+% 
+%       % plot humeral head sphere
+%       n = 20; % number of lines
+%       [X,Y,Z] = sphere(n);
+%       X = X*obj.humerus.radius;
+%       Y = Y*obj.humerus.radius;
+%       Z = Z*obj.humerus.radius;
+%       Xc = obj.humerus.center(1);
+%       Yc = obj.humerus.center(2);
+%       Zc = obj.humerus.center(3);
+%       mesh(X+Xc,Y+Yc,Z+Zc,'LineStyle','none','FaceColor','magenta','FaceAlpha',0.3,'FaceLighting','gouraud');
+% 
+% 
+%       % plot humeral head centerline
+%       pt1 = obj.scapula.glenoid.center;
+%       pt2 = obj.humerus.center;
+%       x = [pt1(1), pt2(1)];
+%       y = [pt1(2), pt2(2)];
+%       z = [pt1(3), pt2(3)];
+%       plot3(x,y,z,'Color','magenta','LineWidth', 5);
 
 
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Graph properties
 
       axis off; % remove axis
       lightPosition = obj.scapula.glenoid.center + ...
           (...
           obj.scapula.coordSys.ML + ...
           obj.scapula.coordSys.IS + ...
           obj.scapula.coordSys.PA ...
           ) * 100;
       light('Position',lightPosition,'Style','local');
       grid on;
       xlabel('X (CT)');
       ylabel('Y (CT)');
       zlabel('Z (CT)');
       axis square;
       axis vis3d; % fixed limits and scaling
       viewPoint = obj.scapula.coordSys.ML + ...
                   obj.scapula.coordSys.IS + ...
                   obj.scapula.coordSys.PA;
       viewPoint = obj.scapula.coordSys.PA;
       view(viewPoint); % view from lateral side
       % Align IS axis with window vertical axis
       camup(obj.scapula.coordSys.IS);
 
       figureHandle = gca;
       figureHandle.DataAspectRatio = [1 1 1]; % Same relative length of data units along each axis
       figureHandle.Box = 'On';
       set(gcf,'color','w'); % Set background color to white
 
       output = 1;
     end
   end
 end
diff --git a/ShoulderCase/CoordinateSystem.m b/ShoulderCase/CoordinateSystem.m
index abf03f6..28a53db 100644
--- a/ShoulderCase/CoordinateSystem.m
+++ b/ShoulderCase/CoordinateSystem.m
@@ -1,236 +1,236 @@
 classdef CoordinateSystem < handle
 
 
 
   properties
     origin
     xAxis
     yAxis
     zAxis
-    deviate
+    
   end
 
 
 
   properties (Hidden = true)
     rotationMatrix
 
     % enable retro-compatibility with codes that uses PA, IS, and ML axes
     PA
     IS
     ML
   end
 
 
 
   methods
 
     function obj = CoordinateSystem()
 %       obj.origin    = [0 0 0];
 %       obj.xAxis     = [1 0 0];
 %       obj.yAxis     = [0 1 0];
 %       obj.zAxis     = [0 0 1];
 
     end
 
     function set.xAxis(obj,value)
       % enable retro-compatibility with codes that uses a PA axis
       obj.xAxis               = value;
       obj.PA                  = value;
       obj.rotationMatrix(:,1) = value';
     end
 
     function set.yAxis(obj,value)
       % enable retro-compatibility with codes that uses an IS axis
       obj.yAxis               = value;
       obj.IS                  = value;
       obj.rotationMatrix(:,2) = value';
     end
 
     function set.zAxis(obj,value)
       % enable retro-compatibility with codes that uses a ML axis
       obj.zAxis               = value;
       obj.ML                  = value;
       obj.rotationMatrix(:,3) = value';
     end
 
     function output = getRotationMatrix(obj)
       output = obj.rotationMatrix;
     end
 
     function localPoints = express(obj,globalPoints)
       % Give the coordinates of global cartesian points in present coordinate system
       assert(size(globalPoints,2) == 3,'Points have to be a Nx3 double array.')
       localPoints = (globalPoints-obj.origin) * obj.rotationMatrix;
     end
 
     function globalPoints = getGlobalCoordinatesOfLocalPoints(obj,localPoints)
         % Give the coordinates of local points in global cartesian coordinate system
         assert(size(localPoints,2) == 3,'Points have to be a Nx3 double array.')
         globalPoints = (localPoints * obj.rotationMatrix') + obj.origin;
     end
 
     function projected = projectOnXYPlane(obj,points)
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.xAxis);(obj.origin + obj.yAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function projected = projectOnYZPlane(obj,points)
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.yAxis);(obj.origin + obj.zAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function projected = projectOnZXPlane(obj,points)
       projectionPlane = Plane;
       projectionPlane.fit([obj.origin;(obj.origin + obj.zAxis);(obj.origin + obj.xAxis)]);
       projected = projectionPlane.projectOnPlane(points);
     end
 
     function output = isOrthogonal(obj)
       % The axes comparison can't be done looking for strict equality due to
       % rounded values. Therefore the values must be evaluated with a given
       % tolerance
       minAngle = (pi/2)-(10^-6); % Arbitrarily defined as a deviation of
                                  % a millionth of a radian
       if dot(obj.xAxis,obj.yAxis) > norm(obj.xAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
       if dot(obj.xAxis,obj.zAxis) > norm(obj.xAxis)*norm(obj.zAxis)*cos(minAngle)
         output = false;
         return
       end
       if dot(obj.zAxis,obj.yAxis) > norm(obj.zAxis)*norm(obj.yAxis)*cos(minAngle)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isRightHanded(obj)
       % The comparison between a theoretical right-handed axis and the actual
       % z Axis can't be done looking for strict vector equality due to rounded
       % values. Therefore the norm of the sum of the two vectors is evaluated.
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       rightHandedAxis = cross(obj.xAxis,obj.yAxis);
       if norm(rightHandedAxis+obj.zAxis) < norm(rightHandedAxis)
         output = false;
         return
       end
       output = true;
     end
 
     function output = isLeftHanded(obj)
       if not(obj.isOrthogonal)
         output = false;
         return
       end
       if obj.isRightHanded
         output = false;
         return
       end
       output = true;
     end
 
     function setSystemWithScapulaLandmarks(obj,scapula)
       % Calculate the (EPFL) scapular coordinate system
       % The EPFL scapular coordinate system is defined for a right scapula see
       % paper (10.1302/0301-620X.96B4.32641). The X axis is
       % antero-posterior, the Y axis is infero-superior, the Z axis
       % is meddio-lateral. This system is right-handed.
       % This orientation corresponds approximatively to the ISB
       % recommendadtion (doi:10.1016/j.jbiomech.2004.05.042).
       % For left scapula we keep the x axis as postero-anterior and
       % get a left-handed coordinate system.
 
       if (length(scapula.groove) < 2)
         error(['Can''t set scapular coordinate system with actual',...
                ' scapula.groove']);
         return;
       end
 
       % The four following methods have to be called in this specific order
       obj.setXAxisWithScapulaLandmarks(scapula);
       obj.setZAxisWithScapulaLandmarks(scapula);
       obj.setYAxisWithScapulaLandmarks(scapula);
       obj.setOriginWithScapulaLandmarks(scapula);
     end
 
     function setSystemWith(obj)
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       % New method to set a coordinate system can be written here
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
       %
     end
 
   end
 
 
 
   methods (Access = private, Hidden = true)
 
     function setOriginWithScapulaLandmarks(obj,scapula)
       % Set the origin of the coordinate system to the spino-gemoid notch
       % projected on the scapular axis
       spinoGlenoid         = scapula.spinoGlenoidNotch;
       grooveMeanProjection = mean(scapula.plane.projectOnPlane(scapula.groove));
       obj.origin = grooveMeanProjection + ...
                    dot((spinoGlenoid - grooveMeanProjection),obj.zAxis)*obj.zAxis;
     end
 
     function setXAxisWithScapulaLandmarks(obj,scapula)
       % Set the postero-anterior axis to be the scapula plane normal
       obj.xAxis = scapula.plane.normal;
     end
 
     function setYAxisWithScapulaLandmarks(obj,scapula)
       % Set the infero-superior axis to be orthogonal with the two other axes
       superior    = scapula.spinoGlenoidNotch;
       inferior    = scapula.angulusInferior;
       obj.yAxis   = cross(obj.xAxis,obj.zAxis);
       obj.yAxis   = orientVectorToward(obj.yAxis,(superior-inferior));
     end
 
     function setZAxisWithScapulaLandmarks(obj,scapula)
       % Set the media-lateral axis to be the line fitted to the projection of
       % the groove points on the scapula plane
       lateral     = scapula.spinoGlenoidNotch;
       medial      = scapula.trigonumSpinae;
       groovePointsProjection = scapula.plane.projectOnPlane(scapula.groove);
       grooveAxis  = fitLine(groovePointsProjection);
       grooveAxis  = (grooveAxis/norm(grooveAxis))';
       obj.zAxis   = orientVectorToward(grooveAxis,(lateral-medial));
     end
 
   end
 
 
 
 end