diff --git a/ShoulderCase/CoordinateSystem.m b/ShoulderCase/CoordinateSystem.m
index 8b3bf95..74d483c 100644
--- a/ShoulderCase/CoordinateSystem.m
+++ b/ShoulderCase/CoordinateSystem.m
@@ -1,232 +1,235 @@
 classdef CoordinateSystem < handle
 
 
 
   properties
     origin
     xAxis
     yAxis
     zAxis
   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 points = express(obj,points)
-      % Give the coordinates of points in present coordinate system
-      if (size(points,2) ~= 3)
-        error('Points have to be a Nx3 double array.')
-        return;
-      end
-      points = (points-obj.origin) * obj.rotationMatrix;
+    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