diff --git a/CT/dicomCopier_gui.fig b/CT/dicomCopier_gui.fig
deleted file mode 100644
index df6de93..0000000
Binary files a/CT/dicomCopier_gui.fig and /dev/null differ
diff --git a/CT/dicomCopier_gui.m b/CT/dicomCopier_gui.m
deleted file mode 100644
index 52b4519..0000000
--- a/CT/dicomCopier_gui.m
+++ /dev/null
@@ -1,448 +0,0 @@
-function varargout = dicomCopier_gui(varargin)
-% DICOMCOPIER_GUI MATLAB code for dicomCopier_gui.fig
-%      DICOMCOPIER_GUI, by itself, creates a new DICOMCOPIER_GUI or raises the existing
-%      singleton*.
-%
-%      H = DICOMCOPIER_GUI returns the handle to a new DICOMCOPIER_GUI or the handle to
-%      the existing singleton*.
-%
-%      DICOMCOPIER_GUI('CALLBACK',hObject,eventData,handles,...) calls the local
-%      function named CALLBACK in DICOMCOPIER_GUI.M with the given input arguments.
-%
-%      DICOMCOPIER_GUI('Property','Value',...) creates a new DICOMCOPIER_GUI or raises the
-%      existing singleton*.  Starting from the left, property value pairs are
-%      applied to the GUI before dicomCopier_gui_OpeningFcn gets called.  An
-%      unrecognized property name or invalid value makes property application
-%      stop.  All inputs are passed to dicomCopier_gui_OpeningFcn via varargin.
-%
-%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
-%      instance to run (singleton)".
-%
-% See also: GUIDE, GUIDATA, GUIHANDLES
-
-% Edit the above text to modify the response to help dicomCopier_gui
-
-% Last Modified by GUIDE v2.5 31-Oct-2013 09:01:18
-
-% Begin initialization code - DO NOT EDIT
-gui_Singleton = 1;
-gui_State = struct('gui_Name',       mfilename, ...
-    'gui_Singleton',  gui_Singleton, ...
-    'gui_OpeningFcn', @dicomCopier_gui_OpeningFcn, ...
-    'gui_OutputFcn',  @dicomCopier_gui_OutputFcn, ...
-    'gui_LayoutFcn',  [] , ...
-    'gui_Callback',   []);
-if nargin && ischar(varargin{1})
-    gui_State.gui_Callback = str2func(varargin{1});
-end
-
-if nargout
-    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
-else
-    gui_mainfcn(gui_State, varargin{:});
-end
-% End initialization code - DO NOT EDIT
-
-
-% --- Executes just before dicomCopier_gui is made visible.
-function dicomCopier_gui_OpeningFcn(hObject, eventdata, handles, varargin)
-% This function has no output args, see OutputFcn.
-% hObject    handle to figure
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-% varargin   command line arguments to dicomCopier_gui (see VARARGIN)
-
-% Choose default command line output for dicomCopier_gui
-handles.output = hObject;
-
-% Update handles structure
-guidata(hObject, handles);
-
-% UIWAIT makes dicomCopier_gui wait for user response (see UIRESUME)
-% uiwait(handles.figure1);
-movegui('center')
-
-
-% --- Outputs from this function are returned to the command line.
-function varargout = dicomCopier_gui_OutputFcn(hObject, eventdata, handles)
-% varargout  cell array for returning output args (see VARARGOUT);
-% hObject    handle to figure
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Get default command line output from handles structure
-varargout{1} = handles.output;
-
-
-% --- Executes on button press in loadbutton.
-function loadbutton_Callback(hObject, eventdata, handles)
-% hObject    handle to loadbutton (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% make clean
-if isfield(handles, 'myinfo');
-    handles = rmfield(handles, 'myinfo');
-end
-
-if isfield(handles, 'mydicom');
-    handles = rmfield(handles, 'mydicom');
-end
-
-%[handles.dicomlist, handles.dicomdir] = uigetfile('D:\DICOM\*.*','Select the images to load','MultiSelect', 'on');
-[handles.dicomlist, handles.dicomdir] = uigetfile('C:\Users\dewarrat\Documents\current projects\Scapula segmentation\raw\Cases 12.2014 n2\*.*','Select the images to load','MultiSelect', 'on');
-if ~handles.dicomdir
-    return
-end
-handles.N = length(handles.dicomlist);
-
-set(handles.slider_slice, 'Max', handles.N)
-set(handles.slider_slice, 'Value', int32(handles.N/2))
-set(handles.slider_slice, 'SliderStep',[1/handles.N 10/handles.N])
-set(handles.slider_slice, 'Min', 1)
-
-set(handles.info_text, 'ForegroundColor', 'black', 'FontWeight', 'normal')
-
-% Read dicom data
-handles.myinfo(1) = dicominfo(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(1))));
-
-set(handles.IPP_text, 'String', sprintf('%s (%s%s)', handles.myinfo(1).PatientID, handles.myinfo(1).PatientName.FamilyName(1), handles.myinfo(1).PatientName.GivenName(1)))
-set(handles.resolution_text, 'String', sprintf('%f mm', handles.myinfo(1).PixelSpacing(1)))
-set(handles.CTdate_text, 'String', sprintf('%s', handles.myinfo(1).AcquisitionDate))
-
-
-if isfield(handles.myinfo(1), 'ConvolutionKernel')
-    if strfind(handles.myinfo(1).ConvolutionKernel, 'BONE')
-        set(handles.radiobutton0, 'Value', 0)
-        set(handles.radiobutton1, 'Value', 1)
-        set(handles.radiobutton2, 'Value', 0)
-    elseif strfind(handles.myinfo(1).ConvolutionKernel, 'STANDARD')
-        set(handles.radiobutton0, 'Value', 0)
-        set(handles.radiobutton1, 'Value', 0)
-        set(handles.radiobutton2, 'Value', 1)
-    else
-        set(handles.radiobutton0, 'Value', 1)
-        set(handles.radiobutton1, 'Value', 0)
-        set(handles.radiobutton2, 'Value', 0)
-    end
-else
-    set(handles.radiobutton0, 'Value', 1)
-    set(handles.radiobutton1, 'Value', 0)
-    set(handles.radiobutton2, 'Value', 0)
-end
-
-% pre-allocate ?
-% handles.mydicom = zeros(handles.myinfo(1).Height, handles.myinfo(1).Width, handles.N,'int16');
-handles.mydicom(:,:,int32(handles.N/2)) = dicomread(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(int32(handles.N/2)))));
-
-handles.slice = int32(handles.N / 2);
-handles.CTmin = -200;
-handles.CTmax = 2000;
-
-handles.bigmax = 3000;
-handles.bigmin = -2000;
-
-set(handles.slider_contrast, 'Max', handles.bigmax - handles.bigmin - 1)
-set(handles.slider_contrast, 'Min', 0)
-set(handles.slider_contrast, 'Value', 900)
-
-set(handles.slider_brightness, 'Max', handles.bigmax)
-set(handles.slider_brightness, 'Min', handles.bigmin)
-set(handles.slider_brightness, 'Value', 1100)
-
-set(handles.slice_num, 'String', handles.slice)
-imshow(handles.mydicom(:,:,handles.slice), [handles.CTmin handles.CTmax])
-guidata(hObject,handles)
-
-
-% --- Executes on slider movement.
-function slider_slice_Callback(hObject, eventdata, handles)
-% hObject    handle to slider_slice (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hints: get(hObject,'Value') returns position of slider
-%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
-
-handles.slice = int32(get(hObject,'Value'));
-set(handles.slice_num, 'String', handles.slice)
-
-handles.mydicom(:,:,handles.slice) = dicomread(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(handles.slice))));
-imshow(handles.mydicom(:,:,handles.slice), [handles.CTmin handles.CTmax])
-
-guidata(hObject,handles)
-
-
-% --- Executes during object creation, after setting all properties.
-function slider_slice_CreateFcn(hObject, eventdata, handles)
-% hObject    handle to slider_slice (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    empty - handles not created until after all CreateFcns called
-
-% Hint: slider controls usually have a light gray background.
-if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
-    set(hObject,'BackgroundColor',[.9 .9 .9]);
-end
-
-
-
-function patient_num_Callback(hObject, eventdata, handles)
-% hObject    handle to patient_num (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hints: get(hObject,'String') returns contents of patient_num as text
-%        str2double(get(hObject,'String')) returns contents of patient_num as a double
-
-
-% --- Executes during object creation, after setting all properties.
-function patient_num_CreateFcn(hObject, eventdata, handles)
-% hObject    handle to patient_num (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    empty - handles not created until after all CreateFcns called
-
-% Hint: edit controls usually have a white background on Windows.
-%       See ISPC and COMPUTER.
-if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
-    set(hObject,'BackgroundColor','white');
-end
-
-
-
-
-% --- Executes on slider movement.
-function slider_contrast_Callback(hObject, eventdata, handles)
-% hObject    handle to slider_contrast (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hints: get(hObject,'Value') returns position of slider
-%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
-
-if ~isfield(handles, 'mydicom');
-    return
-end
-
-handles.CTmax = (get(handles.slider_brightness,'Max') + get(handles.slider_brightness,'Min') - get(handles.slider_brightness,'Value')) + (get(hObject,'Max') - get(hObject,'Value'))/2;
-handles.CTmin = (get(handles.slider_brightness,'Max') + get(handles.slider_brightness,'Min') - get(handles.slider_brightness,'Value')) - (get(hObject,'Max') - get(hObject,'Value'))/2;
-
-if handles.CTmax > handles.bigmax
-    handles.CTmax = handles.bigmax;
-end
-
-if handles.CTmin < handles.bigmin
-    handles.CTmin = handles.bigmin;
-end
-
-imshow(handles.mydicom(:,:,handles.slice), [handles.CTmin handles.CTmax])
-
-guidata(hObject,handles)
-
-% --- Executes during object creation, after setting all properties.
-function slider_contrast_CreateFcn(hObject, eventdata, handles)
-% hObject    handle to slider_contrast (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    empty - handles not created until after all CreateFcns called
-
-% Hint: slider controls usually have a light gray background.
-if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
-    set(hObject,'BackgroundColor',[.9 .9 .9]);
-end
-
-
-
-% --- Executes on slider movement.
-function slider_brightness_Callback(hObject, eventdata, handles)
-% hObject    handle to slider_brightness (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hints: get(hObject,'Value') returns position of slider
-%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
-
-if ~isfield(handles, 'mydicom');
-    return
-end
-
-handles.CTmax = (get(hObject,'Max') + get(hObject,'Min') - get(hObject,'Value')) + (get(handles.slider_contrast,'Max') - get(handles.slider_contrast,'Value'))/2;
-handles.CTmin = (get(hObject,'Max') + get(hObject,'Min') - get(hObject,'Value')) - (get(handles.slider_contrast,'Max') - get(handles.slider_contrast,'Value'))/2;
-
-if handles.CTmax > handles.bigmax
-    handles.CTmax = handles.bigmax;
-end
-
-if handles.CTmin < handles.bigmin
-    handles.CTmin = handles.bigmin;
-end
-
-imshow(handles.mydicom(:,:,handles.slice), [handles.CTmin handles.CTmax])
-
-guidata(hObject,handles)
-
-
-% --- Executes during object creation, after setting all properties.
-function slider_brightness_CreateFcn(hObject, eventdata, handles)
-% hObject    handle to slider_brightness (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    empty - handles not created until after all CreateFcns called
-
-% Hint: slider controls usually have a light gray background.
-if isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
-    set(hObject,'BackgroundColor',[.9 .9 .9]);
-end
-
-
-% --- Executes on button press in copy_button.
-function copy_button_Callback(hObject, eventdata, handles)
-% hObject    handle to copy_button (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-sortPictures = 1;
-
-if isempty(get(handles.patient_num,'String'))
-    set(handles.info_text, 'ForegroundColor', 'red', 'FontWeight', 'bold', 'String', 'Please give a patient number')
-    return
-else
-    patient = get(handles.patient_num,'String');
-end
-
-if isempty(get(handles.output_dir,'String'))
-    outputdir = uigetdir;
-    if outputdir
-        set(handles.output_dir,'String', outputdir)
-    end
-else
-    outputdir = get(handles.output_dir,'String');
-end
-
-if ~outputdir
-    set(handles.info_text, 'ForegroundColor', 'red', 'FontWeight', 'bold', 'String', 'Please select an output directoy')
-    return
-end
-
-set(handles.info_text, 'ForegroundColor', 'black', 'FontWeight', 'normal')
-
-% Create name string
-newname = [patient, '-', handles.myinfo(1).PatientID, '_', handles.myinfo(1).PatientName.FamilyName(1), handles.myinfo(1).PatientName.GivenName(1)];
-if get(handles.radiobutton1,'Value')
-    newname = [newname, '-bone.'];
-elseif get(handles.radiobutton2,'Value')
-    newname = [newname, '-soft.'];
-else
-    newname = [newname, '.'];
-end
-
-% Create output directory
-outputdir = [outputdir, '\', patient, '-', handles.myinfo(1).PatientID, '\', 'CT-', patient, '-', handles.myinfo(1).PatientID, '-xxx\dicom\'];
-mkdir(outputdir)
-
-% Sort pictures in correct order
-if sortPictures == 1
-    for h = 1:handles.N
-        infoTemp = dicominfo(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(h))));
-        SliceLoc(h,1) = h;
-        SliceLoc(h,2) = infoTemp.SliceLocation;
-    end
-    SliceSorted = sortrows(SliceLoc,2);
-end
-
-% Read all data
-for i = 1:handles.N
-    str = sprintf('Copying dicom...\n%d / %d\n', i, handles.N);
-    set(handles.info_text, 'String', str)
-    drawnow
-    % Anonymous
-    if sortPictures == 1
-        caseID = SliceSorted(i,1);
-    else
-        caseID = i;
-    end
-     
-    handles.myinfo(i) = dicominfo(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(caseID))));
-    
-    if get(handles.radiobutton1,'Value')
-        handles.myinfo(i).PatientName.GivenName = [handles.myinfo(i).PatientName.FamilyName(1), handles.myinfo(i).PatientName.GivenName(1), '_bone'];
-    elseif get(handles.radiobutton2,'Value')
-        handles.myinfo(i).PatientName.GivenName = [handles.myinfo(i).PatientName.FamilyName(1), handles.myinfo(i).PatientName.GivenName(1), '_soft'];
-    else
-        handles.myinfo(i).PatientName.GivenName = [handles.myinfo(i).PatientName.FamilyName(1), handles.myinfo(i).PatientName.GivenName(1)];
-    end
-    
-    handles.myinfo(i).PatientName.FamilyName = [patient, '_', handles.myinfo(1).PatientID];
-
-    % Read and copy
-    dicomwrite(dicomread(sprintf('%s%s', handles.dicomdir, char(handles.dicomlist(caseID)))),[outputdir, newname, sprintf('%04d',i), '.dcm'], handles.myinfo(i));
-end
-
-set(handles.info_text, 'String', 'Copied')
-
-
-function output_dir_Callback(hObject, eventdata, handles)
-% hObject    handle to output_dir (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hints: get(hObject,'String') returns contents of output_dir as text
-%        str2double(get(hObject,'String')) returns contents of output_dir as a double
-
-
-% --- Executes during object creation, after setting all properties.
-function output_dir_CreateFcn(hObject, eventdata, handles)
-% hObject    handle to output_dir (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    empty - handles not created until after all CreateFcns called
-
-% Hint: edit controls usually have a white background on Windows.
-%       See ISPC and COMPUTER.
-if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor'))
-    set(hObject,'BackgroundColor','white');
-end
-
-
-% --- Executes on button press in radiobutton1.
-function radiobutton1_Callback(hObject, eventdata, handles)
-% hObject    handle to radiobutton1 (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hint: get(hObject,'Value') returns toggle state of radiobutton1
-
-if get(hObject,'Value')
-    set(handles.radiobutton0,'Value',0)
-    set(handles.radiobutton2,'Value',0)
-else
-    set(hObject,'Value',1)
-end
-
-
-% --- Executes on button press in radiobutton2.
-function radiobutton2_Callback(hObject, eventdata, handles)
-% hObject    handle to radiobutton2 (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hint: get(hObject,'Value') returns toggle state of radiobutton2
-
-if get(hObject,'Value')
-    set(handles.radiobutton0,'Value',0)
-    set(handles.radiobutton1,'Value',0)
-else
-    set(hObject,'Value',1)
-end
-
-
-% --- Executes on button press in radiobutton0.
-function radiobutton0_Callback(hObject, eventdata, handles)
-% hObject    handle to radiobutton0 (see GCBO)
-% eventdata  reserved - to be defined in a future version of MATLAB
-% handles    structure with handles and user data (see GUIDATA)
-
-% Hint: get(hObject,'Value') returns toggle state of radiobutton0
-
-if get(hObject,'Value')
-    set(handles.radiobutton1,'Value',0)
-    set(handles.radiobutton2,'Value',0)
-else
-    set(hObject,'Value',1)
-end
diff --git a/CT/readme.txt b/CT/readme.txt
index 02077e5..fe20383 100755
--- a/CT/readme.txt
+++ b/CT/readme.txt
@@ -1 +1,11 @@
-The matlab function infoCT in this folder only works with matlab R2015a on a Windows system.
+infoCT:
+ - The matlab function infoCT in this folder only works with matlab R2015a on a Windows system.
+ 
+ - Must be run from /shoulder/methods/matlab/database including folder CT/ to the path.
+
+infoCT-->CTdicomInfo:
+ - CTdicomInfo replaces infoCT to be used in latest Matlab version 2017b
+
+dicomCopier:
+ - dicomCopier_v001: working version used before 2018
+ - dicomCopier_v002: updated in fev2018
\ No newline at end of file
diff --git a/XLS_MySQL_DB/patientFromExcel.m b/XLS_MySQL_DB/patientFromExcel.m
index 6564400..0699074 100644
--- a/XLS_MySQL_DB/patientFromExcel.m
+++ b/XLS_MySQL_DB/patientFromExcel.m
@@ -1,106 +1,107 @@
 function [patient] = patientFromExcel(rawExcel)
 %PATIENT builds the structure array patients from
 %   Detailed explanation goes here
 
 fprintf('\nGet patient from Excel... ');
 
 % define the patient structure array
 patient = struct(...
     'patient_id', [], ...
     'IPP',        [], ... % to be later move to another table for anonymous constrain
     'initials',   [], ... % to be later move to another table for anonymous constrain
     'anonym_birth_date',   [], ... % to be later move to another table for anonymous constrain
     'gender',     [], ...
     'birth_date', [], ...
     'height',     [], ...
     'weight',     [],...  % should we put weight here since it can vary?
     'comment',     []...  
     );
 
 % get column of variables
 header = rawExcel(1,:);
 sCase_id_col   = find(strcmp([header],'sCase.id'));
 IPP_col        = find(strcmp([header],'anonymity.IPP'));
 initials_col   = find(strcmp([header],'anonymity.initials'));
 anonym_birth_date_col   = find(strcmp([header],'anonymity.birth_date'));
 gender_col     = find(strcmp([header],'patient.gender'));
 birth_date_col = find(strcmp([header],'patient.birth_date'));
 height_col     = find(strcmp([header],'patient.height'));
 weight_col     = find(strcmp([header],'patient.weight'));
 comment_col    = find(strcmp([header],'patient.comment'));
 % loop over rows of Excel table to build patient structure
 % patients are not included if caseExcel, IPP, initials, gender, birth_date are not
 % defined
 [rowN, ~] = size(rawExcel);
 patient_id = 0;
 for row_idx = 2:rowN
+    row_idx
     row = rawExcel(row_idx, :); % get the entire row
     
     % check that there is a sCase for this row
     sCase_id  = row{sCase_id_col}; % sCase_id from Excel
     if ~isempty(sCase_id)
         
         % check validity of data
         IPP        = row{IPP_col};
         initials   = row{initials_col};
         gender        = row{gender_col};
         anonym_birth_date = row{anonym_birth_date_col};
         anonym_birth_date = datetime(anonym_birth_date,'InputFormat','dd.MM.yyyy','ConvertFrom','excel');
         anonym_birth_date.Format = 'yyyy-MM-dd';
         birth_date = row{birth_date_col};
         birth_date = datetime(birth_date,'InputFormat','dd.MM.yyyy','ConvertFrom','excel');
         birth_date.Format = 'yyyy-MM-dd';
         % get muscle comment
         comment = row{comment_col};
         if isnan(comment)
             comment = '';
         end
         
         if ~isnan(IPP) & ...
                 ~isnan(initials) & ...
                 ~isnan(gender) & ...
                 ~isnat(birth_date) & ...
                 ~isnat(anonym_birth_date)
             
             % Check that patient in row sCase is not already in structure
             % We asume here that IPP uniquely identifies a patient
             sameIPP = find([patient.IPP] == IPP);
             if ~isempty(sameIPP)
                 % IPP is allready in patient array, extra-checks to avoid false same patient
                 if birth_date ~= patient(sameIPP).birth_date % not same birthdate
                     fprintf('\nSame IPP (%i) but different birth date in case %s !', IPP, sCase_id);
                 end
                 if gender ~= patient(sameIPP).gender % not same gender
                     fprintf('\nSame IPP (%i) but different gender in case %s !', IPP, sCase_id);
                 end
                 if ~strcmp(initials, patient(sameIPP).initials) % not same initials
                     fprintf('\nSame IPP (%i) but different initials in case %s !', IPP, sCase_id);
                 end
             else % patient is not in patient structure array, so add it
                 patient_id = patient_id + 1;
                 patient_idx = patient_id;
                 
                 height = row{height_col};
                 if height == ' ' % to avoid strange behavior of the excel import
                     height = NaN;
                 end
                 weight = row{weight_col};
                 
                 patient(patient_idx).patient_id = patient_id;
                 patient(patient_idx).IPP        = IPP;
                 patient(patient_idx).initials   = initials;
                 patient(patient_idx).gender     = gender;
                 patient(patient_idx).birth_date = birth_date;
                 patient(patient_idx).height     = height;
                 patient(patient_idx).weight     = weight;
                 patient(patient_idx).anonym_birth_date = anonym_birth_date;
                 patient(patient_idx).comment = comment;
             end
         end
     end
 end
 
 fprintf('Done\n');
 
 end
 
diff --git a/anatomy/HumeruSphere.m b/anatomy/HumeruSphere.m
index 0892dfd..e2aa57c 100644
--- a/anatomy/HumeruSphere.m
+++ b/anatomy/HumeruSphere.m
@@ -1,118 +1,118 @@
 %% Humerusphere
 %
 % This script creates a tcl text file named "SphereScript.tcl" to display the humerus sphere in
 % Amira.
 % This script is used in the "Anatomical measurements with Amira".
 
 % Subject is the patient number
 % Type is 'normal' or 'pathologic'
 
 
 
 function HumeruSphere(subject,type)
 
-directory = '../../../data'; % Define directory
+directory = '../../../../data'; % Define directory
 
 switch type
     
     
     case 'pathologic'
         type = 'P';
-        folder = sprintf('P%03d',subject); % format the patient name to P###
+        folder = sprintf('P%03d',subject) % format the patient name to P###
         subjectString = num2str(subject);
         levelDir1 = str2num(subjectString(1));
         levelDir2 = str2num(subjectString(2));
         
         % Loop over the list of patients and extract the complete name of the specific patient
-        listing = dir(sprintf('%s/%s/%d/%d', directory, type, levelDir1, levelDir2));
+        listing = dir(sprintf('%s/%s/%d/%d', directory, type, levelDir1, levelDir2))%;
         for i = 1:1:numel(listing)
             name = listing(i).name;
             t = strfind (name, folder);
             if t == 1;
                 patientName = listing(i).name;
             end
         end
         
         % Import HumeralHead landmarks
         if exist(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder),'file') == 2
             newData1 = importdata(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder), ' ', 14);
             HH = newData1.data;
             
         elseif exist(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder),'file') == 2
             warning('Using old humeral landmarks definition')
             newData1 = importdata(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder), ' ', 14);
             HH = newData1.data;
         else
             error('MATLAB:rmpath:DirNotFound',...
-                'Could not find humeral head landmarks file: %s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder)
+                'Could not find humeral head landmarks file: %s/%s/%d/%d/%s/CT-%s-1b/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder)
         end
         
         
         
     case 'normal'
         type = 'N';
         folder = sprintf('N%03d',subject); % format the patient name to N###
         subjectString = num2str(subject);
         levelDir1 = str2num(subjectString(1));
         levelDir2 = str2num(subjectString(2));
         
         % Loop over the list of patients and extract the complete name of the specific patient
         listing = dir(sprintf('%s/%s/%d/%d', directory, type, levelDir1, levelDir2));
         for i = 1:1:numel(listing)
             name = listing(i).name;
             t = strfind (name, folder);
             if t == 1;
                 patientName = listing(i).name;
             end
         end
         
         if exist(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder),'file') == 2
             newData1 = importdata(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder), ' ', 14);
             HH = newData1.data;
             
         elseif exist(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder),'file') == 2
             warning('Using old humeral landmarks definition')
             newData1 = importdata(sprintf('%s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder), ' ', 14);
             HH = newData1.data;
         else
             error('MATLAB:rmpath:DirNotFound',...
-                'Could not find humeral head landmarks file: %s/%s/%d/%d/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder)
+                'Could not find humeral head landmarks file: %s/%s/%d/%d/%s/CT-%s-1b/amira/HHLandmarks%s.landmarkAscii',directory, type, levelDir1, levelDir2, patientName,patientName,folder)
         end
         
     otherwise
         error('Unexpected shoulder type.');
 end
 
 %% Fit sphere on Humeral Head landmarks
 % [HC,HHRadius,HHResiduals] = fitSphere(HH);
 [HC,HHRadius,~, ~] = fitSphere2(HH);
 
 mkdir('Generated_Amira_TCL_Scripts/HumerusSphereScripts')
 
 
 fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\HumerusSphereScripts\\SphereScript_%s.tcl',folder),'w');
 fprintf(fid,'remove CreateSphere%s HumSphere%s.surf SphereView%s Intersect\n', folder, folder, folder);
 fprintf(fid,'create HxCreateSphere {CreateSphere%s}\n', folder);
 fprintf(fid,'CreateSphere%s setIconPosition 20 520\n', folder);
 fprintf(fid,'CreateSphere%s radius setMinMax 0 50\n', folder);
 fprintf(fid,'CreateSphere%s radius setValue %f\n', folder, HHRadius);
 fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', folder, HC(1));
 fprintf(fid,'CreateSphere%s coords setValue 1 %f\n', folder, HC(2));
 fprintf(fid,'CreateSphere%s coords setValue 2 %f\n', folder, HC(3));
 fprintf(fid,'CreateSphere%s fire\n', folder);
 fprintf(fid,'[ {CreateSphere%s} create\n ] setLabel {HumSphere%s.surf}\n', folder, folder);
 fprintf(fid,'HumSphere%s.surf setIconPosition 200 520\n', folder);
 fprintf(fid,'HumSphere%s.surf master connect CreateSphere%s\n', folder, folder);
 fprintf(fid,'HumSphere%s.surf fire\n', folder);
 fprintf(fid,'create HxDisplaySurface {SphereView%s}\n', folder);
 fprintf(fid,'SphereView%s setIconPosition 400 520\n', folder);
 fprintf(fid,'SphereView%s data connect HumSphere%s.surf\n', folder, folder);
 fprintf(fid,'SphereView%s drawStyle setValue 4\n', folder);
 fprintf(fid,'SphereView%s fire\n', folder);
 fprintf(fid,'create HxOverlayGrid Intersect\n');
 fprintf(fid,'{Intersect} {data} connect HumSphere%s.surf\n', folder);
 fprintf(fid,'{Intersect} {module} connect ObliqueSlice4\n');
 fprintf(fid,'{Intersect} {lineWidth} setIndex 0 2\n');
 fprintf(fid,'ObliqueSlice4 setPlane %f %f %f 1 0 0 0 1 0\n', HC);
 fprintf(fid,'ObliqueSlice4 fire\n');
 fclose(fid);
\ No newline at end of file
diff --git a/anatomy/scapula_addmeasure.m b/anatomy/scapula_addmeasure.m
index b9ffeb4..47e2e9f 100644
--- a/anatomy/scapula_addmeasure.m
+++ b/anatomy/scapula_addmeasure.m
@@ -1,91 +1,93 @@
+
 %% scapula_measure
 %
 % This script returns additional measurements from the Amira files.
 % 3 differents output are possible:
 
 %  - 'References' creates a text file in References folder containing coordinates to built in Solidworks 2 Coordinates
 %     Systems. One for Abaqus and one for the scapula. There are two
 %     additional coordinates to built the glenoid centerline.
 %
 %  - 'density' creates a textfile in CylinderReferences folder containing
 %     an Amira script used in the density measurements protocol. It is needed to adjust the Cylinder at the correct position and with the correct size.
 %
 %  - 'obliqueSlice' creates a textfile in obliqueSlice folder containing an
 %     Amira script used in Amira to obtain an image of the Friedman plane (2D
 %     glenoid orientation measurement) and an image for muscles measurement.
 %
 %  - 'display' creates a textfile in "display" folder containing an Amira
 %     script used to display the glenoid sphere and the scapula plane.
 
 
 % Author: EPFL-LBO
 % Date: 2016-09-05
 %
 %%
 
 % Output is 'References' or 'obliqueSlice' 'display' or 'density'
 % Type is 'normal' or 'pathologic'
 % Subject is the patient number
 
 
 function scapula_addmeasure(output,type,subject)
 
-directory = '../../../data';  % Define directory
+directory = 'Y:/data';
+%'../../../data';  % Define directory
 
 switch type
     case 'pathologic'
         location = 'P';
         CTn = sprintf('P%d',str2num(subject));   % format the patient name to P###
         levelDir1 = str2num(subject(1));
         levelDir2 = str2num(subject(2));
         
         
-        finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
+        finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2)
         
         listing = dir(sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2));
         for i = 1:1:numel(listing)
             name = listing(i).name;
             t = strfind(name, CTn);
             if t == 1;
                 patientName = listing(i).name;
             end
         end
         
         
     case 'normal'
         location = 'N';
         CTn = sprintf('N%03d',subject);   % format the patient name to N###
         patient = num2str(subject);
         levelDir1 = str2num(patient(1));
         levelDir2 = str2num(patient(2));
         finalDirectory = sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2);
         
         listing = dir(sprintf('%s/%s/%d/%d', directory, location, levelDir1, levelDir2));
         for i = 1:1:numel(listing)
             name = listing(i).name;
             t = strfind(name, CTn);
             if t == 1;
                 patientName = listing(i).name;
             end
         end
 end
 
 switch output
     
     case 'References'
         output = 'References';
         scapula_calculation(patientName, finalDirectory, location, output);
         
     case 'obliqueSlice'
         output = 'obliqueSlice';
         scapula_calculation(patientName, finalDirectory, location, output);
         
     case 'density'
         output = 'density';
         scapula_calculation(patientName, finalDirectory, location, output);
         
     case 'display'
         output = 'display';
         scapula_calculation(patientName, finalDirectory, location, output);
 end
 end
\ No newline at end of file
diff --git a/anatomy/scapula_calculation.m b/anatomy/scapula_calculation.m
index 444e3fb..a2ae17c 100644
--- a/anatomy/scapula_calculation.m
+++ b/anatomy/scapula_calculation.m
@@ -1,1155 +1,1152 @@
 %% scapula_calculation
 %
 % This script makes the anatomic calculation of the glenoid from the Amira files.
 % It is used by the scripts scapula_measure.m and scapula_addmeasure.
 % It should normally not be used directly.
 % Author: EPFL-LBO
 % Date: 2016-11-18
 %
 %%
-
-% Subject is the patient number
+% INPUTS
+% Subject: is the patient number
 % directory is the location of the datas
 % location it 'Pathologic' or 'Normal'
 % output is 'References', 'obliqueSlice', 'density' and 'display'
 
 
-
-
-
-
 function [Result] = scapula_calculation(subject, directory, location, output)
 
 segmentedScapula = 0;
-
-CTn = strtok(subject,'-');
-CTn = strtok(CTn, location);
-CTn = str2num(CTn); %#ok<ST2NM>
-CTn = sprintf('%s%03d',location,CTn);
+subject
+CTn = strtok(subject,'-')
+CTn = strtok(CTn, location)
+CTn = str2num(CTn) %#ok<ST2NM>
+CTn = sprintf('%s%03d',location,CTn)
 
 
 % Import Data
-
+ls
 
 % Import Glenoid Surface
 if exist(sprintf('%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn),'file') == 2
     [p,~,~]=importStl(sprintf('%s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn),1);
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find glenoid surface file: %s/%s/CT-%s-1/amira/ExtractedSurface%s.stl',directory,subject,subject,CTn)
 end
 
 % Import Scapula Surface
 if exist(sprintf('%s/%s/CT-%s-1/amira/scapula_%s.stl',directory,subject,subject,CTn),'file') == 2
     [pscapula,~,~]=importStl(sprintf('%s/%s/CT-%s-1/amira/scapula_%s.stl',directory,subject,subject,CTn),1);
     segmentedScapula = 1;
 end
 
 % Import HumeralHead landmarks
 if exist(sprintf('%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData1 = importdata(sprintf('%s/%s/CT-%s-1/amira/newHHLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     HH = newData1.data;
 else
     warning('Using old humeral landmarks definition')
     if exist(sprintf('%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
         newData1 = importdata(sprintf('%s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
         HH = newData1.data;
     else
         error('MATLAB:rmpath:DirNotFound',...
             'Could not find humeral head landmarks file: %s/%s/CT-%s-1/amira/HHLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
     end
 end
 
 % Import Scapula landmarks
 % AI: Angulus Inferior
 % TS: Trigonum Spinae Scapulae (the midpoint of the triangular surface on the medial border of the scapula in line with the scaoular spine
 % PC: Processus Coracoideus (most lateral point)
 % AL: Acromion Lateral (most lateral part on acromion) // Before, it was AC: Acromio-clavicular joint (most lateral part on acromion)
 % AA: Angulus Acromialis (most laterodorsal point)
 % AG: Acromion-Glenoid notch
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData2 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     landmarks = newData2.data;
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find scapula landmarks file: %s/%s/CT-%s-1/amira/ScapulaLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 
 end
 
 % Import Scapula Pillar
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData3 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     SP = newData3.data; %Scapula Pillar points
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find axillary border landmarks file: %s/%s/CT-%s-1/amira/ScapulaPillarLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 end
 
 % Import Scapula Groove
 if exist(sprintf('%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn),'file') == 2
     newData4 = importdata(sprintf('%s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn), ' ', 14);
     SG = newData4.data; %Scapula Groove points
 else
     error('MATLAB:rmpath:DirNotFound',...
         'Could not find supraspinatus fossa landmarks file: %s/%s/CT-%s-1/amira/ScapulaGrooveLandmarks%s.landmarkAscii',directory,subject,subject,CTn)
 end
 
 % Define anatomical axes
 
 MedLat = landmarks(6,:)-mean(SG);           % Z ?
 AntPos = landmarks(5,:)-landmarks(3,:);     % X ?
 InfSup = landmarks(6,:)-mean(SP);           % Y ?
 
 
 % Find scapula side
 
 A = cross(MedLat,AntPos);
 B = vrrotvec(InfSup,A);    %Comparing Z-axes
 if rad2deg(B(4))>90
     Side = 0; %Right
     p(:,1) = -p(:,1);
     if segmentedScapula == 1
         pscapula(:,1) = -pscapula(:,1); 
     end
     HH(:,1) = -HH(:,1);
     landmarks(:,1) = -landmarks(:,1);
     SG(:,1) = -SG(:,1);
     MedLat = landmarks(6,:)-mean(SG);
     AntPos = landmarks(5,:)-landmarks(3,:);
 else
     Side = 1; %Left
 end
 
 
 %% Fit plane and axis
 
 % Find scapular plane (fit plane to AI and TS and most distal SG)
 % The 10 next line are used to determine the scapulaPlaneNormal regarding
 % the 3 most distal points of the scapulagroove.
 
 
 AI = landmarks(1,:); % Angulus inferioris point
 TS = landmarks(2,:); % Trigonum Spinae point
 EuclidDistance = zeros(5,1);
 for i=1:5 % Calculate the distance between the TS point and each of the 5 groove points
     EuclidDistance(i) = sqrt((landmarks(2,1)-SG(i,1))^2 + (landmarks(2,2)-SG(i,2))^2 + (landmarks(2,3)-SG(i,3))^2);
 end
 
 maxEuclidianDistance = max(EuclidDistance); % Find the largest distance
 pos = EuclidDistance == maxEuclidianDistance;
 GRL = SG(pos,:); % set most lateral point of the groove
 [ScapulaPlaneNormal, PlaneMean, ~, PlaneRMSE, ~] = fitPlane2([AI;TS;GRL]);
 
 
 % Find scapular plane (fit plane to SP and SG)
 %[ScapulaPlaneNormal, PlaneMean, ~, PlaneRMSE, R2Plane] = fitPlane2([SG;SP]);
 
 C = vrrotvec(ScapulaPlaneNormal,AntPos);
 if rad2deg(C(4))>90
     ScapulaPlaneNormal = -ScapulaPlaneNormal;
 end
 
 % First project, then fit
 SG_proj = project2Plane(SG,ScapulaPlaneNormal',[0 0 0],size(SG,1)); % project SG points on scapular plane
 [GrooveAxis, GrooveMean, ~, ~, ~] = fitLine2(SG_proj);
 TransverseAxis = (GrooveAxis/norm(GrooveAxis))'; % TransverseAxis is the norm of the line that was created by the projection of SG points on scapular plane
 
 D = vrrotvec(TransverseAxis,MedLat);
 if rad2deg(D(4))>90
-    TransverseAxis = -TransverseAxis;
+    TransverseAxis = -TransverseAxis
 end
 
 % Project GrooveMean to scapular plane
 TransverseAxisMean = project2Plane(GrooveMean,ScapulaPlaneNormal',PlaneMean,1);
 
 % Project AG to transverseAxis
 Origin = TransverseAxisMean + dot((landmarks(6,:)-TransverseAxisMean),TransverseAxis)*TransverseAxis;
 
 % Scapula axial plane normal
 AxialPlane = cross(TransverseAxis, ScapulaPlaneNormal);
 E = vrrotvec(AxialPlane, [0 0 1]);
 CTorientation = rad2deg(E(4));
 if CTorientation > 90
     CTorientation = 180 - CTorientation;
 end
 
 % Closest point. GC is the closet point between the mean of the extracted
 % surface and all the points of the extracted surface.
 % GC = Glenoid Centre
 GC = mean(p);
 vect = zeros(size(p));
 distance3D = zeros(size(p,1),1);
 for i=1:size(p,1)
     vect(i,:) = p(i,:)-GC;
     distance3D(i) = norm(vect(i,:));
 end
 ClosestNode = distance3D == min(distance3D);
 GC = p(ClosestNode,:);
 
 
 % Fit Sphere on Glenoid Surface
 [GlenoidSphere,GlenoidRadius,GlenoidResiduals, ~] = fitSphere2(p);
 GlenoidRMSE =  norm(GlenoidResiduals)/sqrt(length(GlenoidResiduals));
 
 
 % Fit sphere on Humeral Head landmarks
 [HC,HHRadius, ~, ~] = fitSphere2(HH);
 
 
 % Glenoid version 3D
 GO = GlenoidSphere'-GC;
 
 Version3DRot = vrrotvec(TransverseAxis,GO);
 Version3D = rad2deg(Version3DRot(4));
 
 VersionDirectionRotX = vrrotvec(Version3DRot(1:3),ScapulaPlaneNormal);
 VersionDirectionRotZ = vrrotvec(Version3DRot(1:3),cross(TransverseAxis,ScapulaPlaneNormal));
 
 if rad2deg(VersionDirectionRotX(4)) > 90
     VersionDirection = rad2deg(VersionDirectionRotZ(4));
 else
     VersionDirection = -rad2deg(VersionDirectionRotZ(4));
 end
 if VersionDirection < -180
     VersionDirection = VersionDirection + 360;
 end
 if VersionDirection > 180
     VersionDirection = VersionDirection - 360;
 end
 
 
 % Maximum wear plane
 MaxWearPlane = cross(GO, TransverseAxis);
 F = vrrotvec(MaxWearPlane, [0 0 1]);
 WearPlaneAngle = rad2deg(F(4));
 if WearPlaneAngle > 90
     WearPlaneAngle = 180 - WearPlaneAngle;
 end
 
 
 % Inclination : Angle between GO projected on the scapular plane and the medio-lateral axis
 proj_GO = GO' - dot(GO', ScapulaPlaneNormal) * ScapulaPlaneNormal;
 InclinationRot = vrrotvec(TransverseAxis,proj_GO);
 Inclination = rad2deg(InclinationRot(4));
 if dot(InclinationRot(1:3), ScapulaPlaneNormal) > 0
     Inclination = -Inclination;
 end
 
 
 % Version 2D : Angle between GO projected on the scapular AXIAL plane and the medio-lateral axis
 proj_GO2 = GO' - dot(GO', AxialPlane') * AxialPlane';
 Version2DRot = vrrotvec(TransverseAxis,proj_GO2);
 Version2D = rad2deg(Version2DRot(4));
 if dot(Version2DRot(1:3), AxialPlane') > 0
     Version2D = -Version2D;
 end
 
 
 % Scapulo Humeral Subluxation Index
 HO = HC'-GC;    % glenoid center to humeral center
 SHeccentricity = HO-(dot(HO,TransverseAxis)*TransverseAxis);% vector from HC to transverse axis (perpendicular)
 SHSI = norm(SHeccentricity) / HHRadius; % Mod. on 04.06.2014
 
 SHSIDirectionRotX = vrrotvec(SHeccentricity,ScapulaPlaneNormal);
 SHSIDirectionRotZ = vrrotvec(SHeccentricity,cross(TransverseAxis,ScapulaPlaneNormal));
 
 if rad2deg(SHSIDirectionRotZ(4)) < 90
     SHSIDirection = rad2deg(SHSIDirectionRotX(4));
 else
     SHSIDirection = -rad2deg(SHSIDirectionRotX(4));
 end
 if SHSIDirection < -180
     SHSIDirection = SHSIDirection+360;
 end
 
 SHSPlane = cross(TransverseAxis, HO);
 G = vrrotvec(SHSPlane, [0 0 1]);
 SHSPlaneAngle = rad2deg(G(4));
 if SHSPlaneAngle > 90
     SHSPlaneAngle = 180 - SHSPlaneAngle;
 end
 
 
 % Gleno Humeral Subluxation Index
 GHeccentricity = HO-dot(HO,GO/norm(GO))*(GO/norm(GO));%vector from glenoid sphere centre to glenoid centerline (perpendicular)
 GHSI = norm(GHeccentricity) / HHRadius; % Mod. on 04.06.2014
 
 ScapulaPlaneNormal2Glenoid = project2Plane(ScapulaPlaneNormal',GO,[0 0 0],1);% project y axis on glenoid plane
 
 GHSIDirectionRotX = vrrotvec(GHeccentricity,ScapulaPlaneNormal2Glenoid);
 GHSIDirectionRotZ = vrrotvec(GHeccentricity,cross(GO,ScapulaPlaneNormal2Glenoid));
 
 if rad2deg(GHSIDirectionRotZ(4)) < 90
     GHSIDirection = rad2deg(GHSIDirectionRotX(4));
 else
     GHSIDirection = -rad2deg(GHSIDirectionRotX(4));
 end
 
 if GHSIDirection < -180
     GHSIDirection = GHSIDirection+360;
 end
 
 GHSPlane = cross(GO, HO);
 H = vrrotvec(GHSPlane, [0 0 1]);
 GHSPlaneAngle = rad2deg(H(4));
 if GHSPlaneAngle > 90
     GHSPlaneAngle = 180 - GHSPlaneAngle;
 end
 
 
 % Height of the glenoid cap
 % Project the points of the surface on the centerline
 p_proj = zeros(size(p));
 for i = 1:size(p,1)
     p_proj(i,:) = GC + (dot(GO,(p(i,:)-GC)) / dot(GO,GO)) * GO;
     distance3D(i) = norm(GC-p_proj(i,:));
 end
 
 CapHeight = max(distance3D);
 
 
 % Glenoid height and width
 
 % Scapula coordinate system in CT frame
 ScapulaCSYS(1,:) = TransverseAxis; % X
 ScapulaCSYS(2,:) = ScapulaPlaneNormal; % Y
 ScapulaCSYS(3,:) = AxialPlane; % Z
 ScapulaCSYS(4,:) = Origin;
 ScapulaCSYS(5,:) = Origin + 10*ScapulaCSYS(1,:);
 ScapulaCSYS(6,:) = Origin + 10*ScapulaCSYS(2,:);
 ScapulaCSYS(7,:) = Origin + 10*ScapulaCSYS(3,:);
 
 % Contruction of transformation matrix from CT frame to Scapula
 % frame
 Osc = ScapulaCSYS(4,:);
 Oscx = Origin + ScapulaCSYS(1,:);
 Oscy = Origin + ScapulaCSYS(2,:);
 Oscz = Origin + ScapulaCSYS(3,:);
 
 O = [0,0,0];
 Ox = [1,0,0];
 Oy = [0,1,0];
 Oz = [0,0,1];
 
 A = [Osc(1),Osc(2),Osc(3),1;
     Oscx(1),Oscx(2),Oscx(3),1;
     Oscy(1),Oscy(2),Oscy(3),1;
     Oscz(1),Oscz(2),Oscz(3),1];
 
 bx = [O(1);
     Ox(1);
     Oy(1);
     Oz(1)];
 
 by = [O(2);
     Ox(2);
     Oy(2);
     Oz(2)];
 
 
 bz = [O(3);
     Ox(3);
     Oy(3);
     Oz(3)];
 
 
 B = [bx,by,bz];
 X = A\B;
 T = [X';0,0,0,1]; % Transformation matrix from CT frame to scapula frame
 
 
 psc = zeros(size(p,1),4);
 for i = 1:size(p,1)
     psc(i,:) = T*[p(i,:),1]';
 end
 
 psc_y = psc(:,2);
 psc_z = psc(:,3);
 
 % Difference of extremal values
 i_y_max=1;
 for i = 2:1:size(psc_y)
 
     if psc_y(i) > psc_y(i_y_max)
         i_y_max = i;
     end
 end
 
 i_y_min=1;
 for i = 2:1:size(psc_y)
 
     if psc_y(i) < psc_y(i_y_min)
         i_y_min = i;
     end
 end
 
 i_z_max=1;
 for i = 2:1:size(psc_z)
 
     if psc_z(i) > psc_z(i_z_max)
         i_z_max = i;
     end
 end
 
 i_z_min=1;
 for i = 2:1:size(psc_z)
 
     if psc_z(i) < psc_z(i_z_min)
         i_z_min = i;
     end
 end
 
 scapulaWidth = norm(p(i_y_max,:)-p(i_y_min,:));
 scapulaHeight = norm(p(i_z_max,:)-p(i_z_min,:));
 
 % Acromion Index (AI)
 % AI = GA/GH, where:
 % GA: GlenoAcromial distance in scapula plane = distance from AL to glenoid principal axis
 % GH: GlenoHumeral distance = 2*HHRadius, humeral head diameter (radius*2) 
 
 %Calculate the correct AL landmarks if the scapula is segmented
 if segmentedScapula == 1
     % Transform scapula, glenoid and acromion points to Scapula Coordinate System
     ex = [1;0;0];
     ey = [0;1;0];
     ez = [0;0;1];  
 
     ev = AxialPlane'; % Y
     ew = TransverseAxis'; %Z
     eu = cross(ev, ew);  % X
 
     R = eu * ex' + ev * ey' + ew * ez';
     
     % Transform scapula point to Scapula Coordinate System
     ScapulaPtsinScapulaCSYS = pscapula*R;
     % Take the most lateral point (assumption: The most lateral point of the scaplua is always on the acromion
     [~,ALpositionInScapula] = max(ScapulaPtsinScapulaCSYS(:,3)); 
     AL = pscapula(ALpositionInScapula,:);% AL: Acromion lateral (most lateral part on acromion) 
 else
    AL = landmarks(4,:); % AL: Acromion lateral (most lateral part on acromion) 
 end
 
 % Project AL, glenoid points and GC in scapula plane
 ALinScapulaPlane = project2Plane(AL,ScapulaPlaneNormal',PlaneMean,size(AL,1));
 GlenoidPtsinScapulaPlane = project2Plane(p,ScapulaPlaneNormal',PlaneMean,size(p,1));
 GCinScapulaPlane = project2Plane(GC,ScapulaPlaneNormal',PlaneMean,size(GC,1));
 if segmentedScapula == 1
     ScapulaPtsinScapulaPlane = project2Plane(pscapula,ScapulaPlaneNormal',PlaneMean,size(pscapula,1));
 end
 
+% CSA calculation (work done by Bharath
+[radCSA, degCSA] = computeCSA(GlenoidPtsinScapulaPlane, ALinScapulaPlane);
+
 % Figure in 3D to show position of scapula plane, glenoid points and AC
 
 figure;
 d = PlaneMean*ScapulaPlaneNormal;
 if segmentedScapula == 1
     [xx,yy]=ndgrid(min(pscapula(:,1)):max(pscapula(:,1)),min(pscapula(:,2)):max(pscapula(:,2)));
 else
     [xx,yy]=ndgrid((min(p(:,1))-50):(max(p(:,1))+50),(min(p(:,2))-50):(max(p(:,2))+50));
 end    
 zz=(d-ScapulaPlaneNormal(1)*xx-ScapulaPlaneNormal(2)*yy)/ScapulaPlaneNormal(3);
 hold on
 surf(xx,yy,zz,'FaceColor','r','FaceAlpha',0.5, 'EdgeColor', 'none')
 if segmentedScapula == 1
     scatter3(ScapulaPtsinScapulaPlane(:,1),ScapulaPtsinScapulaPlane(:,2),ScapulaPtsinScapulaPlane(:,3),200,[0.5 0.5 0.5],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)    
     hold on;
 end
 scatter3(ALinScapulaPlane(:,1),ALinScapulaPlane(:,2),ALinScapulaPlane(:,3),[],[0 0 1],'filled')
 hold on
 scatter3(GlenoidPtsinScapulaPlane(:,1),GlenoidPtsinScapulaPlane(:,2),GlenoidPtsinScapulaPlane(:,3),200,[0 0 1],'Marker','.')
 hold on
 scatter3(GCinScapulaPlane(:,1),GCinScapulaPlane(:,2),GCinScapulaPlane(:,3),[],'r','filled')
 hold on
 daspect([1 1 1])
 xlabel('X (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 zlabel('Z (Scapula CSYS)')
 daspect([1 1 1])
 if segmentedScapula == 1
     savefig(sprintf('%s/%s/CT-%s-1/amira/Projection2ScapulaPlaneWithSegmentedScapula_%s',directory,subject,subject,CTn))
 else
     savefig(sprintf('%s/%s/CT-%s-1/amira/Projection2ScapulaPlane_%s',directory,subject,subject,CTn))
 end
 
 % Transform scapula, glenoid and acromion points to Scapula Coordinate System
 
 ex = [1;0;0];
 ey = [0;1;0];
 ez = [0;0;1];  
 
 ev = AxialPlane'; % Y
 ew = TransverseAxis'; %Z
 eu = cross(ev, ew);  % X
 
 R = eu * ex' + ev * ey' + ew * ez'; % Create rotation matrix to transform points to Scapula Coordinate System
 
 ALinScapulaCSYS = AL*R;
 GlenoidPtsinScapulaCSYS = p*R;
 GCinScapulaCSYS = GC*R;
 TrinScapulaCSYS = TransverseAxis*R;
 ScapulaPlaneNormalinScapulaCSYS = ScapulaPlaneNormal'*R;
 PlaneMeaninScapulaCSYS = PlaneMean*R;
 if segmentedScapula == 1
     ScapulaPtsinScapulaCSYS = pscapula*R;
 end
 
 % Figure in 3D to show position of scapula plane, glenoid points and AC
 
 figure;
 if segmentedScapula == 1
     scatter3(ScapulaPtsinScapulaCSYS(:,1),ScapulaPtsinScapulaCSYS(:,2),ScapulaPtsinScapulaCSYS(:,3),200,[0.5 0.5 0.5],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)    
     hold on;
 end
 scatter3(ALinScapulaCSYS(:,1),ALinScapulaCSYS(:,2),ALinScapulaCSYS(:,3),[],[0 0 1],'filled')
 hold on
 scatter3(GlenoidPtsinScapulaCSYS(:,1),GlenoidPtsinScapulaCSYS(:,2),GlenoidPtsinScapulaCSYS(:,3),200,[0 0 1],'Marker','.')
 hold on
 scatter3(GCinScapulaCSYS(:,1),GCinScapulaCSYS(:,2),GCinScapulaCSYS(:,3),[],'r','filled')
 hold on
 x=eu*linspace(0,50);
 plot3(x(:,1)',x(:,2)',x(:,3)','r')
 hold on
 y=ev*linspace(0,50);
 plot3(y(:,1)',y(:,2)',y(:,3)','g')
 hold on
 z=ew*linspace(0,50);
 plot3(z(:,1)',z(:,2)',z(:,3)','b')
 daspect([1 1 1])
 d= ScapulaPlaneNormalinScapulaCSYS*PlaneMeaninScapulaCSYS';
 if segmentedScapula == 1
     [yy,zz]=ndgrid(min(ScapulaPtsinScapulaCSYS(:,2)):max(ScapulaPtsinScapulaCSYS(:,2)),min(ScapulaPtsinScapulaCSYS(:,3)):max(ScapulaPtsinScapulaCSYS(:,3)));
 else
     [yy,zz]=ndgrid((min(GlenoidPtsinScapulaCSYS(:,2))-50):(max(GlenoidPtsinScapulaCSYS(:,2))+50),(min(GlenoidPtsinScapulaCSYS(:,3))-50):(max(GlenoidPtsinScapulaCSYS(:,3))+50));
 end    
 xx = zeros(size(yy));
 xx(:,:) = d/ScapulaPlaneNormalinScapulaCSYS(1); %assuming normal(3)==0 and normal(2)==0
 hold on
 surf(xx,yy,zz,'FaceColor','r','FaceAlpha',0.5, 'EdgeColor', 'none');
 xlabel('X (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 zlabel('Z (Scapula CSYS)')
 daspect([1 1 1])
 if segmentedScapula == 1
     savefig(sprintf('%s/%s/CT-%s-1/amira/Transform2ScapulaCSYSwithScapulaSegmented_%s',directory,subject,subject,CTn))
 else
     savefig(sprintf('%s/%s/CT-%s-1/amira/Transform2ScapulaCSYS_%s',directory,subject,subject,CTn))
 end
 
 % Project scapula, glenoid, AC, GC and scapula axis to scapula plane in the
 % Scapula Coordinate System
 
 ALinScapulaPlane = project2Plane(ALinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ALinScapulaCSYS,1));
 GlenoidPtsinScapulaPlane = project2Plane(GlenoidPtsinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(GlenoidPtsinScapulaCSYS,1));
 GCinScapulaPlane = project2Plane(GCinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(GCinScapulaCSYS,1));
 TrinScapulaPlane = project2Plane(TrinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(TrinScapulaCSYS,1));
 if segmentedScapula == 1
     ScapulaPtsinScapulaPlane = project2Plane(ScapulaPtsinScapulaCSYS,ScapulaPlaneNormalinScapulaCSYS,PlaneMeaninScapulaCSYS,size(ScapulaPtsinScapulaCSYS,1));
 end
 
 % Find glenoid principal axis (most superior and most inferior projected
 % glenoid points)
 
 GlenoidPrincipalAxis = [GlenoidPtsinScapulaPlane(GlenoidPtsinScapulaPlane(:,2) == min(GlenoidPtsinScapulaPlane(:,2)),:) ;GlenoidPtsinScapulaPlane(GlenoidPtsinScapulaPlane(:,2) == max(GlenoidPtsinScapulaPlane(:,2)),:)];
 
 % Compute GA (distance from AL to glenoid principal axis)
-
 GA = norm(cross(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:),ALinScapulaPlane-GlenoidPrincipalAxis(1,:)))/norm(GlenoidPrincipalAxis(2,:)-GlenoidPrincipalAxis(1,:));
-% Compute GH (Humeral head diameter)
 
+% Compute GH (Humeral head diameter)
 GH = 2*HHRadius;
 
 % Compute AI
-
 AI = GA/GH;
 
 % Figure in 2D showing projected glenoid points, glenoid principal axis and
 % AC point
-
 figure;
 daspect([1 1 1])
 if segmentedScapula == 1
     scatter(ScapulaPtsinScapulaPlane(:,3),ScapulaPtsinScapulaPlane(:,2),200,[0.8 0.8 0.8],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)
     hold on
 end
 scatter(ALinScapulaPlane(:,3),ALinScapulaPlane(:,2),[],[0 0 1],'filled')
 hold on
 scatter(GlenoidPtsinScapulaPlane(:,3),GlenoidPtsinScapulaPlane(:,2),200,[0 0 1],'Marker','.','MarkerFaceAlpha',0.5,'MarkerEdgeAlpha',0.5)
 hold on
 scatter(GlenoidPrincipalAxis(:,3),GlenoidPrincipalAxis(:,2),[],'w','filled','MarkerEdgeColor','k')
 hold on
 scatter(GCinScapulaPlane(:,3),GCinScapulaPlane(:,2),[],'r','filled')
 hold on
 a=TrinScapulaPlane(:,2)/TrinScapulaPlane(:,3);
 b=GCinScapulaPlane(:,2)-a*GCinScapulaPlane(:,3);
 if segmentedScapula == 1
     x=linspace(min(ScapulaPtsinScapulaPlane(:,3)),max(ScapulaPtsinScapulaPlane(:,3)));
 else
     x=linspace((min(GlenoidPtsinScapulaPlane(:,3))-50),(max(GlenoidPtsinScapulaPlane(:,3)))+50);    
 end
 y=a*x+b;
 hold on
 plot(x,y,'--r')
 hold on
 a1=(GlenoidPrincipalAxis(2,2)-GlenoidPrincipalAxis(1,2))/(GlenoidPrincipalAxis(2,3)-GlenoidPrincipalAxis(1,3));
 b1=GlenoidPrincipalAxis(1,2)-a1*GlenoidPrincipalAxis(1,3);
 x=linspace(min(GlenoidPrincipalAxis(:,3))-5,max(GlenoidPrincipalAxis(:,3))+5);
 y1=a1*x+b1;
 plot(x,y1,'--b')
 daspect([1 1 1])
 xlabel('Z (Scapula CSYS)')
 ylabel('Y (Scapula CSYS)')
 title(['Acromion Index = ' num2str(AI) ' (GA = ' num2str(GA) ' GH = ' num2str(GH) ')'])
 if segmentedScapula == 1
     saveas(gcf,sprintf('%s/%s/CT-%s-1/amira/AIwithSegmentedScapula_%s.png',directory,subject,subject,CTn))
 else
     saveas(gcf,sprintf('%s/%s/CT-%s-1/amira/AI_%s.png',directory,subject,subject,CTn))
 end
 
 close all
 
 % Output
 
 CTtoXYZ = [TransverseAxis' ScapulaPlaneNormal cross(TransverseAxis, ScapulaPlaneNormal)'];
 Cxyz = (GC - Origin) * CTtoXYZ;
 
 % Scapula coordinate system in CT frame
 ScapulaCSYS(3,:) = TransverseAxis; % Z
 ScapulaCSYS(1,:) = ScapulaPlaneNormal; % X
 ScapulaCSYS(2,:) = cross(TransverseAxis,ScapulaPlaneNormal); % Y
 ScapulaCSYS(4,:) = Origin;
 
 Result.sCase_id = strtok(subject,'-');%1
 Result.glenoid_Radius = GlenoidRadius;%2
 Result.glenoid_radiusRMSE = GlenoidRMSE;%3
-%Result.glenoidSpehricity = ' ';
+%Result.glenoidSphericity = ' ';
 %Result.glenoid_biconcave = ' ';
 Result.glenoid_depth = CapHeight; %4
 Result.glenoid_width = scapulaWidth;%5
 Result.glenoid_height = scapulaHeight;%6
 Result.glenoid_center_PA = -Cxyz(2); %7
 Result.glenoid_center_IS = Cxyz(3); %8
 Result.glenoid_center_ML = Cxyz(1); %9
 Result.glenoid_version_ampl = Version3D;%10
 Result.glenoid_version_orient = VersionDirection;%11
 Result.glenoid_Version = Version2D; % 12
 Result.glenoid_Inclination = Inclination; %13
 % Result.glenoid_version2D = ' ';
 % Result.glenoid_walch_class = ' ';
 Result.humerus_joint_radius = ' '; %14
 Result.humeral_head_radius = HHRadius;%15
 % Result.humerus_GHsublux_2D = ' ';
 % Result.humerus_SHsublux_2D = ' ';
 Result.humerus_GHsubluxation_ampl = GHSI;%16
 Result.humerus_GHsubluxation_orient = GHSIDirection;%17
 Result.humerus_SHsubluxation_ampl = SHSI;%18
 Result.humerus_SHsubluxation_orient = SHSIDirection;%19
-% Result.scapula_CSA = ' ';
+Result.scapula_CSA = degCSA % radCSA;
 Result.scapula_CTangle = CTorientation; %20
 Result.scapula_CTangleVersion = WearPlaneAngle; %21
 Result.scapula_CTangleSHS = SHSPlaneAngle; % 22
 Result.scapula_CTangleGHS = GHSPlaneAngle; % 23
 Result.scapula_PlaneRMSE = PlaneRMSE;%24
 Result.scapula_AI = AI;
 
+Result
 
 
 %% Additional output
 
 if exist('output','var') == 1
 
     % Scapula coordinate system in CT frame
     ScapulaCSYS(1,:) = TransverseAxis; % X
     ScapulaCSYS(2,:) = ScapulaPlaneNormal; % Y
     ScapulaCSYS(3,:) = AxialPlane; % Z
     ScapulaCSYS(4,:) = Origin;
     ScapulaCSYS(5,:) = Origin + 10*ScapulaCSYS(1,:);
     ScapulaCSYS(6,:) = Origin + 10*ScapulaCSYS(2,:);
     ScapulaCSYS(7,:) = Origin + 10*ScapulaCSYS(3,:);
 
     %% ISB coordinate system in CT frame
 
     Xisb = (landmarks(4,:)-landmarks(2,:))/norm(landmarks(4,:)-landmarks(2,:)); %Med-Lat
     Yisb = cross(Xisb,landmarks(1,:)-landmarks(2,:))/norm(cross(Xisb,landmarks(1,:)-landmarks(2,:)));%Ant-Pos
     Zisb = cross(Xisb,Yisb);%Inf-Sup
 
     %% Locate abq coordinate system in CT frame
 
     Act2isb = [Xisb' Yisb' Zisb'];
     Act2lbo = ScapulaCSYS(1:3,:)';
     Albo2isb = Act2isb/Act2lbo;
-    Athx2abq = SpinCalc('EA321toDCM',[40 0 0]);% scapular plane is 40� anterior to frontal plane
+    Athx2abq = SpinCalc('EA321toDCM',[40 0 0]);% scapular plane is 40? anterior to frontal plane
     Aisb2thx = SpinCalc('EA321toDCM',[-35.6 -17.8 -2.5]);% initial position McClure et al. 2001
 
     Act2abq = Athx2abq*Aisb2thx*Albo2isb*Act2lbo;
 
     %% Output initial position of scapula for abaqus simulation in CT coordinates
 
     %     scaleFactor = 24/HHRadius;  % =24mm / Humerus radius
     ScapulaABQ(1,:) = HC;   % Humerus center
     ScapulaABQ(2,:) = 10*Act2abq(:,1)'+ScapulaABQ(1,:);
     ScapulaABQ(3,:) = 10*Act2abq(:,2)'+ScapulaABQ(1,:);
     ScapulaABQ(4,:) = 10*Act2abq(:,3)'+ScapulaABQ(1,:);
     %     ScapulaABQ(5,1) = scaleFactor;
 
     %% Output coordinates for obliqueSlice
     % The ObliqueSlice can be defined as point + (normal OR u-vector and
     % v-vector)
     FriedmanPlane = [GC 1 0 0 0 1 0];
     obliqueSlice = [landmarks(6,:) ScapulaCSYS(2,:) -ScapulaCSYS(3,:)];      % Plane for muscle measurement
     obliqueSlice2 = [GC ScapulaCSYS(3,:)];                                  % Scapula axial plane
     obliqueSlice3 = [GC ScapulaCSYS(1,:) GO];                               % Max wear plane
     obliqueSlice4 = [GC ScapulaCSYS(2,:)];                                   % Scapular plane
     obliqueSlice5 = [GC HO/norm(HO) GHeccentricity/norm(GHeccentricity)];   % Plane for GHS
     obliqueSlice6 = [GC HO/norm(HO) SHeccentricity/norm(SHeccentricity)];   % Plane for SHS
 
     if Side == 0
         FriedmanPlane(1) = -FriedmanPlane(1);
         obliqueSlice(1:3:end) = - obliqueSlice(1:3:end);
         obliqueSlice2(1:3:end) = -obliqueSlice2(1:3:end);
         obliqueSlice3(1:3:end) = -obliqueSlice3(1:3:end);
         obliqueSlice4(1:3:end) = -obliqueSlice4(1:3:end);
         obliqueSlice5(1:3:end) = -obliqueSlice5(1:3:end);
         obliqueSlice6(1:3:end) = -obliqueSlice6(1:3:end);
     end
 
 
 
 
     switch output
         case 'density'
             %% Data for the calculations of the density in amira
             % Needed variables are :
             %  - CylO: cylinder origin
             %  - CylA: cylinder alignment point
             %  - CylR_factor: scale factor for cylinder radius
             %  - Glenoid_sphere_center
             %  - Glenoid radius
             %  - CylO_shift5mm : cylinder origin shifted 5 mm behind
             %  glenoid surface center
-            %  - CylA2 : cylinder alignment point n�2
+            %  - CylA2 : cylinder alignment point n?2
             %  - CylO_shift15mmfront: cylinder origin shifted 15 mm
             %  frontwards
-            %  - CylA3 : cylinder alignemnt point n�2
+            %  - CylA3 : cylinder alignemnt point n?2
             %  - startcylpoint : origin of the new cylinder starting from
             %  the AG point
             %  - endcylpoint: end of cylinder, 40 mm from the startcylpoint
             %  in the scapula axis
 
 
             CylO = GC;
             Glenoid_sphere_center = GlenoidSphere';
             Glenoid_radius = GlenoidRadius;
 
             %% CylA - second point to form a line with glenoid surface
             % center point which is parallel to scapula axis
             CylA = CylO - 40 * TransverseAxis/norm(TransverseAxis);
 
             %% CylR - distance from glenoid surface center to most
             % eccentric point of glenoid surface
 
             %Projection of surface points onto the same plane
             PlanePoint = CylO;
             PlaneNormal = (CylO-CylA)/norm(CylO-CylA);
             
             ProjectedPoints = zeros(size(p));
             DistanceToCylO = zeros(size(p,1),1);
             for i = 1:size(p,1)
                 ProjectedPoints(i,:) = p(i,:) - dot(p(i,:) - PlanePoint, PlaneNormal) * PlaneNormal;
                 DistanceToCylO(i) = norm(ProjectedPoints(i,:) - CylO);
             end
 
             %Selection of the most eccentric point -> radius of cylinder
             CylR = max(DistanceToCylO);
 
             %Scale factor of the cylinder
             OriginalCylR = 5;
             CylR_factor = CylR/OriginalCylR;
 
             % Coordinates of start and end point of the D10 mm cylinder
             Deepness = 5; %behind glenoid surface
             smallCylR = 5;
             CylO_shift5mm = GC - Deepness * TransverseAxis/norm(TransverseAxis);
             CylA2 = CylO_shift5mm - smallCylR * TransverseAxis;
 
             % To place big cylinder (finer seleciton, exclusion of distant
             % bones)
             CylO_shift15mmfront = GC + 15 * TransverseAxis/norm(TransverseAxis);
             CylA3 = CylO_shift15mmfront - 40 * TransverseAxis/norm(TransverseAxis);
 
             %%Projection of AG on scapula axis
             AG = landmarks(6,:);
             p1 = CylA3;
             p2 = CylO_shift15mmfront;
 
             v1 = p2-p1;
             v2 = AG-p1;
 
             E1 = v1;
             E2 = cross(v1,v2);
             E3 = cross(E1,E2);
 
             e1 = E1/norm(E1);
             e2 = E2/norm(E2);
             e3 = E3/norm(E3);
 
             Oc = p1;
             Ocx = p1 + e1;
             Ocy = p1 + e2;
             Ocz = p1 + e3;
 
             O = [0,0,0];
             Ox = [1,0,0];
             Oy = [0,1,0];
             Oz = [0,0,1];
 
             A = [Oc(1),Oc(2),Oc(3),1;
                 Ocx(1),Ocx(2),Ocx(3),1;
                 Ocy(1),Ocy(2),Ocy(3),1;
                 Ocz(1),Ocz(2),Ocz(3),1];
 
             bx = [O(1);
                 Ox(1);
                 Oy(1);
                 Oz(1)];
 
 
             by = [O(2);
                 Ox(2);
                 Oy(2);
                 Oz(2)];
 
 
             bz = [O(3);
                 Ox(3);
                 Oy(3);
                 Oz(3)];
 
 
             B = [bx,by,bz];
             X = A\B;
             T = [X';
                 0,0,0,1];
 
             det(T);
 
             AG2 = [AG,1];
             AG3=AG2';
             diff = T*AG3;
             projAGcyl=[diff(1),0,0,1];
             startcylpoint=T\projAGcyl';
             startcylpoint = [startcylpoint(1),startcylpoint(2),startcylpoint(3)];
             endcylpoint = startcylpoint + 40 * TransverseAxis/norm(TransverseAxis);
 
             %% Write data
             if Side == 0
                 CylO(1) = -CylO(1);
                 CylA(1) = -CylA(1);
                 Glenoid_sphere_center(1) = -Glenoid_sphere_center(1);
                 CylO_shift5mm(1) = -CylO_shift5mm(1);
                 CylA2(1) = -CylA2(1);
                 CylO_shift15mmfront(1) = -CylO_shift15mmfront(1);
                 CylA3(1) = -CylA3(1);
                 startcylpoint(1) = -startcylpoint(1);
                 endcylpoint(1) = -endcylpoint(1);
 
 
 
             end
             mkdir('Generated_Amira_TCL_Scripts/CylinderReferences')
             
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),'w');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylR_factor, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_sphere_center, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),Glenoid_radius, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift5mm, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA2, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylO_shift15mmfront, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s.dat',CTn),CylA3, '-append','precision',10);
             fclose(fid);
 
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\CylinderReferences\\CylinderReferences_%s_AmiraCode.dat',CTn),'w');
 
             fprintf(fid,'set lineset1 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset1 addPoint %f %f %f]\n',CylO_shift5mm);
             fprintf(fid,'set p2 [$lineset1 addPoint %f %f %f]\n',CylA2);
             fprintf(fid,'$lineset1 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set} setIconPosition 20 10\n');
             fprintf(fid,'Line-Set fire\n\n');
 
             fprintf(fid,'set lineset2 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset2 addPoint 0 0 0]\n');
             fprintf(fid,'set p2 [$lineset2 addPoint 5 0 0]\n');
             fprintf(fid,'$lineset2 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set2} setIconPosition 20 30\n');
             fprintf(fid,'Line-Set2 fire\n\n');
 
             fprintf(fid,'set hideNewModules 0\n');
             fprintf(fid,'create {HxCreateSphere} {GlenoidSphereFit}\n');
             fprintf(fid,'{GlenoidSphereFit} setIconPosition 20 330\n');
             fprintf(fid,'{GlenoidSphereFit} fire\n');
             fprintf(fid,'{GlenoidSphereFit} {type} setValue 0\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setMinMax 0 100\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setIncrement 0.666667\n');
             fprintf(fid,'{GlenoidSphereFit} {radius} setValue %f\n',Glenoid_radius);
             fprintf(fid,'{GlenoidSphereFit} {radius} setSubMinMax 0 100\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
             fprintf(fid,'{GlenoidSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{GlenoidSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setIncrement 0.332667\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setValue 1.07457\n');
             fprintf(fid,'{GlenoidSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setButtons 0\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setValue 1\n');
             fprintf(fid,'{GlenoidSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
             fprintf(fid,'GlenoidSphereFit fire\n');
             fprintf(fid,'GlenoidSphereFit setViewerMask 16383\n');
             fprintf(fid,'GlenoidSphereFit setPickable 1\n\n');
 
             fprintf(fid,'set hideNewModules 0\n');
             fprintf(fid,'create {HxCreateSphere} {CorticalSphereFit}\n');
             fprintf(fid,'{CorticalSphereFit} setIconPosition 20 350\n');
             fprintf(fid,'{CorticalSphereFit} fire\n');
             fprintf(fid,'{CorticalSphereFit} {type} setValue 0\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setMinMax 0 100\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setIncrement 0.666667\n');
             fprintf(fid,'{CorticalSphereFit} {radius} setValue %f\n',Glenoid_radius+3);
             fprintf(fid,'{CorticalSphereFit} {radius} setSubMinMax 0 100\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 0 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 0 %f\n',Glenoid_sphere_center(1));
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 1 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 1 %f\n',Glenoid_sphere_center(2));
             fprintf(fid,'{CorticalSphereFit} {coords} setMinMax 2 -3.40282346638529e+038 3.40282346638529e+038\n');
             fprintf(fid,'{CorticalSphereFit} {coords} setValue 2 %f\n',Glenoid_sphere_center(3));
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setIncrement 0.332667\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setValue 1.07457\n');
             fprintf(fid,'{CorticalSphereFit} {edgeLength} setSubMinMax 0.00999999977648258 5\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setMinMax 0.100000001490116 20\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setButtons 0\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setIncrement 1.32667\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setValue 1\n');
             fprintf(fid,'{CorticalSphereFit} {nopPerUnit2} setSubMinMax 0.100000001490116 20\n');
             fprintf(fid,'CorticalSphereFit fire\n');
             fprintf(fid,'CorticalSphereFit setViewerMask 16383\n');
             fprintf(fid,'CorticalSphereFit setPickable 1\n\n');
 
             fprintf(fid,'set lineset3 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset3 addPoint %f %f %f]\n',startcylpoint);
             fprintf(fid,'set p2 [$lineset3 addPoint %f %f %f]\n',endcylpoint);
             fprintf(fid,'$lineset3 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set3} setIconPosition 20 520\n');
             fprintf(fid,'Line-Set3 fire\n\n');
 
             fprintf(fid,'set lineset4 [create HxLineSet]\n');
             fprintf(fid,'set p1 [$lineset4 addPoint 0 0 0]\n');
             fprintf(fid,'set p2 [$lineset4 addPoint 40 0 0]\n');
             fprintf(fid,'$lineset4 addLine $p1 $p2\n');
             fprintf(fid,'{Line-Set4} setIconPosition 20 540\n');
             fprintf(fid,'Line-Set4 fire\n\n');
 
             fprintf(fid,'{Cylinder2.STL} setScaleFactor 1 %f %f\n',CylR_factor,CylR_factor);
             fprintf(fid,'{Cylinder2.STL} applyTransform\n\n');
 
             fprintf(fid,'saveProject\n');
             fprintf(fid,'clear\n');
             fprintf(fid,'echo "Complete center coordinates of ReamingSphere"\n\n');
 
             %fprintf(fid,'
 
             fclose(fid);
 
         case 'References'
             
             mkdir('Generated_Amira_TCL_Scripts/References')
             
             % References data for a study case
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'w');
             fprintf(fid,'# Scapula coordinate system\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaCSYS(4:7,:), '-append','precision',10);
             fclose(fid);
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
             fprintf(fid,'\n\n# Abaqus reference\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),ScapulaABQ, '-append','precision',10);
             fclose(fid);
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),'a');
             fprintf(fid,'\n\n# Glenoid centerline\n\n');
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GC, '-append','precision',10);
             dlmwrite(sprintf('Generated_Amira_TCL_Scripts\\References\\References_%s.dat',CTn),GlenoidSphere', '-append','precision',10);
             fclose(fid);
 
         case 'obliqueSlice'
             
             [token, remain] = strtok(directory,'../../..');
             directory = ['Z:/' token remain];
             pictureName = [CTn '_new'];
 
             mkdir('Generated_Amira_TCL_Scripts/obliqueSlice')
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\obliqueslice\\obliqueSlice_%s.tcl',CTn),'w');
 
             fprintf(fid,'# Create an ObliqueSlice with Friedman plane #\n');
             fprintf(fid,'ObliqueSlice orientation hit 0\n');
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\n',FriedmanPlane);
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 0 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice fire\n');
             fprintf(fid,'set planepos [ObliqueSlice translate getValue]\n');
             fprintf(fid,'set maxvalue [ObliqueSlice translate getMaxValue]\n');
             fprintf(fid,'set planepos [expr {$maxvalue - int($planepos)}]\n');
             fprintf(fid,'ObliqueSlice createImage %s_$planepos\n', CTn);
             fprintf(fid,'create HxCastField {CastField}\n');
             fprintf(fid,'CastField data connect %s_$planepos\n', CTn);
             fprintf(fid,'CastField fire\n');
             fprintf(fid,'CastField scaling setValue 0 0.21\n');
             fprintf(fid,'CastField scaling setValue 1 200\n');
             fprintf(fid,'CastField fire\n');
             fprintf(fid,'[ {CastField} create ] setLabel %s_$planepos.jpeg\n', CTn);
             if strcmp(location, 'P')
                 fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" %s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
             else
                 fprintf(fid,'%s_$planepos.jpeg exportData "JPEG" Z:/data/%s/%s/CT-%s-1/amira/%s_$planepos.jpeg\n\n', CTn, directory, subject, subject, CTn );
             end
 
             fprintf(fid,'# Get image for muscle measurement #\n');
             fprintf(fid,'ObliqueSlice orientation hit 0\n');
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice fire\n', obliqueSlice);    % Muscle evaluation, through notch
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice createImage {%s}\n', pictureName);
             if strcmp(location, 'pathologic')
                 fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
             else
                 fprintf(fid,'%s exportData "2D Tiff" %s/%s/CT-%s-1/amira/%s.tif\n', pictureName, directory, subject, subject, pictureName );
             end
             fprintf(fid,'saveProject\n\n');
 
             fprintf(fid,'ObliqueSlice2 setPlane %f %f %f %f %f %f\nObliqueSlice2 fire\n', obliqueSlice2);                          % Scapular axial plane
             fprintf(fid,'ObliqueSlice2 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice2 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice3 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice3 fire\n', obliqueSlice3);                            % Max wear plane
             fprintf(fid,'ObliqueSlice3 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice3 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice4 setPlane %f %f %f %f %f %f\nObliqueSlice4 fire\n', obliqueSlice4);                            % Scapula plane
             fprintf(fid,'ObliqueSlice4 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice4 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice5 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice5 fire\n', obliqueSlice5);                            % Scapula plane
             fprintf(fid,'ObliqueSlice5 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice5 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'ObliqueSlice6 setPlane %f %f %f %f %f %f %f %f %f\nObliqueSlice6 fire\n', obliqueSlice6);                            % Scapula plane
             fprintf(fid,'ObliqueSlice6 options setValue 0 1\n');
             fprintf(fid,'ObliqueSlice6 sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n\n');
 
             fprintf(fid,'remove CreateSphere%s HumSphere%s.surf SphereView%s Intersection\n', CTn, CTn, CTn);
             fprintf(fid,'remove CreateSphere%s_2 GlenoidSphere%s.surf SphereView%s_2 Intersection2\n', CTn, CTn, CTn);
             fprintf(fid,'create HxCreateSphere {CreateSphere%s}\n', CTn);
             fprintf(fid,'CreateSphere%s setIconPosition 20 500\n', CTn);
             fprintf(fid,'CreateSphere%s radius setMinMax 0 50\n', CTn);
             fprintf(fid,'CreateSphere%s radius setValue %f\n', CTn, HHRadius);
             if Side == 0
                 fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, -HC(1));
             else
                 fprintf(fid,'CreateSphere%s coords setValue 0 %f\n', CTn, HC(1));
             end
             fprintf(fid,'CreateSphere%s coords setValue 1 %f\n', CTn, HC(2));
             fprintf(fid,'CreateSphere%s coords setValue 2 %f\n', CTn, HC(3));
             fprintf(fid,'CreateSphere%s fire\n', CTn);
             fprintf(fid,'[ {CreateSphere%s} create\n ] setLabel {HumSphere%s.surf}\n', CTn, CTn);
             fprintf(fid,'HumSphere%s.surf setIconPosition 200 500\n', CTn);
             fprintf(fid,'HumSphere%s.surf master connect CreateSphere%s\n', CTn, CTn);
             fprintf(fid,'HumSphere%s.surf fire\n', CTn);
             fprintf(fid,'create HxDisplaySurface {SphereView%s}\n', CTn);
             fprintf(fid,'SphereView%s setIconPosition 400 500\n', CTn);
             fprintf(fid,'SphereView%s data connect HumSphere%s.surf\n', CTn, CTn);
             fprintf(fid,'SphereView%s drawStyle setValue 4\n', CTn);
             fprintf(fid,'SphereView%s fire\n', CTn);
             fprintf(fid,'create HxOverlayGrid {Intersection}\n');
             fprintf(fid,'Intersection module connect ObliqueSlice4\n');
             fprintf(fid,'Intersection data connect HumSphere%s.surf\n', CTn);
             fprintf(fid,'Intersection lineWidth setValue 2\n');
             fprintf(fid,'ObliqueSlice4 setViewerMask 16383\n');
             fprintf(fid,'ObliqueSlice4 setPlane %f %f %f 1 0 0 0 1 0\n', HC);
             fprintf(fid,'ObliqueSlice4 fire\n\n');
 
             fprintf(fid,'create HxCreateSphere {CreateSphere%s_2}\n', CTn);
             fprintf(fid,'CreateSphere%s_2 setIconPosition 20 520\n', CTn);
             fprintf(fid,'CreateSphere%s_2 radius setMinMax 0 50\n', CTn);
             fprintf(fid,'CreateSphere%s_2 radius setValue %f\n', CTn, GlenoidRadius);
             if Side == 0
                 fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, -GlenoidSphere(1));
             else
                 fprintf(fid,'CreateSphere%s_2 coords setValue 0 %f\n', CTn, GlenoidSphere(1));
             end
             fprintf(fid,'CreateSphere%s_2 coords setValue 1 %f\n', CTn, GlenoidSphere(2));
             fprintf(fid,'CreateSphere%s_2 coords setValue 2 %f\n', CTn, GlenoidSphere(3));
             fprintf(fid,'CreateSphere%s_2 fire\n', CTn);
             fprintf(fid,'[ {CreateSphere%s_2} create\n ] setLabel {GlenoidSphere%s.surf}\n', CTn, CTn);
             fprintf(fid,'GlenoidSphere%s.surf setIconPosition 200 520\n', CTn);
             fprintf(fid,'GlenoidSphere%s.surf master connect CreateSphere%s_2\n', CTn, CTn);
             fprintf(fid,'GlenoidSphere%s.surf fire\n', CTn);
             fprintf(fid,'create HxDisplaySurface {SphereView%s_2}\n', CTn);
             fprintf(fid,'SphereView%s_2 setIconPosition 400 520\n', CTn);
             fprintf(fid,'SphereView%s_2 data connect GlenoidSphere%s.surf\n', CTn, CTn);
             fprintf(fid,'SphereView%s_2 drawStyle setValue 4\n', CTn);
             fprintf(fid,'SphereView%s_2 fire\n', CTn);
             fprintf(fid,'create HxOverlayGrid {Intersection2}\n');
             fprintf(fid,'Intersection2 module connect ObliqueSlice5\n');
             fprintf(fid,'Intersection2 data connect GlenoidSphere%s.surf\n', CTn);
             fprintf(fid,'Intersection2 lineWidth setValue 2\n');
             fprintf(fid,'ObliqueSlice5 setViewerMask 16383\n');
             fprintf(fid,'ObliqueSlice5 setPlane %f %f %f 1 0 0 0 1 0\n', GlenoidSphere);
             fprintf(fid,'ObliqueSlice5 fire\n');
 
             fclose(fid);
 
         case 'display'
 
             mkdir('Generated_Amira_TCL_Scripts/display');
             fid = fopen(sprintf('Generated_Amira_TCL_Scripts\\display\\display_%s.tcl',CTn),'w');
             %              if Side == 0;
             %                  PlaneMean(1) = -PlaneMean(1);
             %              end
             obliqueSlice7 = [PlaneMean ScapulaCSYS(2,:)];
             if Side == 0
                 obliqueSlice7(1:3:end) = -obliqueSlice7(1:3:end);
                 GlenoidSphere(1) = -GlenoidSphere(1);
             end
             fprintf(fid,'ObliqueSlice setPlane %f %f %f %f %f %f \nObliqueSlice fire\n', obliqueSlice7);    % Muscle evaluation, through notch
             fprintf(fid,'ObliqueSlice sampling setState item 0 3 item 1 1 item 2 1 item 3 0 item 4 0\n');
             fprintf(fid,'ObliqueSlice fire\n\n');
 
             fprintf(fid,'create HxCreateSphere {GlenoidSphere%s}\n', CTn);
             fprintf(fid,'GlenoidSphere%s setIconPosition 20 550\n', CTn);
             fprintf(fid,'GlenoidSphere%s radius setMinMax 0 50\n', CTn);
             fprintf(fid,'GlenoidSphere%s radius setValue %f\n', CTn, GlenoidRadius);
             fprintf(fid,'GlenoidSphere%s coords setValue 0 %f\n', CTn, GlenoidSphere(1));
             fprintf(fid,'GlenoidSphere%s coords setValue 1 %f\n', CTn, GlenoidSphere(2));
             fprintf(fid,'GlenoidSphere%s coords setValue 2 %f\n', CTn, GlenoidSphere(3));
             fprintf(fid,'GlenoidSphere%s fire\n\n', CTn);
 
             fprintf(fid,'[ {GlenoidSphere%s} create ] setLabel {HumSphereGlenoid%s.surf}\n', CTn, CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf setIconPosition 200 550\n', CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf master connect GlenoidSphere%s\n', CTn, CTn);
             fprintf(fid,'HumSphereGlenoid%s.surf fire\n\n', CTn);
 
             fprintf(fid,'create HxDisplaySurface {GlenoidSphereView%s}\n', CTn);
             fprintf(fid,'GlenoidSphereView%s setIconPosition 400 550\n', CTn);
             fprintf(fid,'GlenoidSphereView%s data connect HumSphereGlenoid%s.surf\n', CTn, CTn);
             fprintf(fid,'GlenoidSphereView%s drawStyle setValue 4\n', CTn);
             fprintf(fid,'GlenoidSphereView%s fire\n', CTn);
             fclose(fid);
         otherwise
             disp('unknown output request')
     end
 end
diff --git a/anatomy/scapula_measure.m b/anatomy/scapula_measure.m
index a0be7f6..1188452 100644
--- a/anatomy/scapula_measure.m
+++ b/anatomy/scapula_measure.m
@@ -1,133 +1,140 @@
 function scapula_measure(caseID)
 %%SCAPULA_MEASURE fills the Shoulder database with scapula anatomical data
 %
 % This function fills the Shoulder database with the scapula anatomical
 % data of the case given as input. The anatomical measurements are stored
 % simultaneously in a XLS file and in the MySQL database.
 %
 % USE: scapula_measure(caseID)
 %
 % INPUT caseID: Case IDs for which you want to compute scapula anatomical
 % data. Cell array that stores the case IDs as a char in the form 'P###' or
 % 'N###' (starts with "P" or "N" followed by 1 to 3 digits).
 %
 % REMARKS The anatomical data are saved in the ScapulaAnatomicalData.xls
 % file.
 %
 % created with MATLAB ver.: 8.6.0.267246 (R2015b) on Windows 7
 % Author: EPFL-LBO-VMC
 % Date: 20-Apr-2017
 %
+DBaccess=0; % Variable for database access deactivation due to  driver problems
     CTDatabaseLocation = '../../../data'; % Location of the CT database
     XLSShoulderDatabaseLocation = '../../../data/Excel/'; % Location of the XLS ShoulderDatabase
 
     % Check validity of input argument (cell array)
     validateattributes(caseID,{'cell'},{'nonempty'});
 
     % Check that imput case IDs are unique
     [caseIDNames,~,uniqueCaseID] = unique(caseID);
     countOfCaseID = hist(uniqueCaseID,unique(uniqueCaseID))';
     freqCaseIDs = struct('name',caseIDNames,'freq',num2cell(countOfCaseID));
     repeatedValues = find([freqCaseIDs.freq] > 1);
 
     if(~isempty(repeatedValues))
         warning('Input contains repeated Case ID: %s. Repeated occurences of that Case ID will be omitted. \n', freqCaseIDs(repeatedValues).name)
         caseID = caseIDNames;
     end
 
     %  open mySQL connection
-    conn = openSQL();
-
+    if DBaccess
+        conn = openSQL();
+    end
     % Get the list of exisiting cases in XLS database
+    exist([XLSShoulderDatabaseLocation 'xlsFromMatlab/anatomy.xls']) % add by JSM for debbuging
     [~,~,currDatabase] = xlsread([XLSShoulderDatabaseLocation 'xlsFromMatlab/anatomy.xls']);
     currDatabaseCaseIDlist = currDatabase(2:end,1);
 
     % --Get CT directory adress from CT database location and case ID--
 
     % Compute scapula anatomical data and fill database
     for i = 1:length(caseID)
 
         % Check validity of input arguments (char and format 'P###' or 'N###')
         validateattributes(caseID{i},{'char'},{'nonempty'});
 
         if (numel(regexp(caseID{i},'^[PN]\d{1,3}$')) == 0)
             error(['Invalid format of CaseID: ' caseID{i} '. CaseID must start with "P" or "N" and be followed by 1 to 3 digits.']);
         end
 
         levelDir1 = caseID{i}(1);
 
         if (length(caseID{i}(2:end)) < 2)
             levelDir2 = '0';
             levelDir3 = '0';
         elseif (length(caseID{i}(2:end)) < 3)
             levelDir2 = '0';
             levelDir3 = caseID{i}(2);
         else
             levelDir2 = caseID{i}(2);
             levelDir3 = caseID{i}(3);
         end
-
-        FindCaseCTFolder = dir([CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' caseID{i} '*']);
+        
+%         [CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' caseID{i} '*']
+        FindCaseCTFolder = dir([CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' caseID{i} '*'])
         if (isempty(FindCaseCTFolder))
             error(['Missing CT directory for CaseID: ' caseID{i}]);
         end
 
         CaseID_IPP = FindCaseCTFolder.name;
         CaseCTFolder = [CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' CaseID_IPP];
 
         % Compute scapula measurements and store in the XLS database
         if exist([CaseCTFolder '/CT-' CaseID_IPP '-1/amira'],'dir') == 7 %if a folder "amira" exist, the patient name is added to the "measured" column
-
             finalDirectory = [CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3] ;
             anatomy = scapula_calculation(CaseID_IPP, finalDirectory, levelDir1); % Function that calculates the anatomical data
 
             % Write to XLS file
             % Replace line if already existing or append
             if(any(strcmp(currDatabaseCaseIDlist, caseID{i})))
                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/anatomy.xls'],struct2cell(anatomy)','Feuil1',['A' int2str(find(strcmp(currDatabaseCaseIDlist, caseID{i}))+1)]);
             else
                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/anatomy.xls'],struct2cell(anatomy)','Feuil1',['A' int2str(length(currDatabaseCaseIDlist)+1+i)]);
             end
+            struct2cell(anatomy)'
 
             % Write to MySQL database
             % Replace line if already existing or issue warning message
-            sqlquery = ['SELECT CT_id FROM CT WHERE shoulder_id IN (SELECT shoulder_id FROM sCase WHERE folder_name = "' caseID{i} '")'];
-            curs = exec(conn,sqlquery);
-            if ~isempty(curs.Message)
-                fprintf('\nMessage: %s\n', curs.Message);
-            end
-
-            curs=fetch(curs);
-            sqlresult = curs.Data;
-
-            if isnumeric(sqlresult{1})
-
-                % Update data in table 'glenoid'
-                data = {sqlresult{1} anatomy.glenoid_Radius anatomy.glenoid_Radius/anatomy.glenoid_radiusRMSE '' anatomy.glenoid_depth anatomy.glenoid_width anatomy.glenoid_height anatomy.glenoid_version_ampl anatomy.glenoid_version_orient anatomy.glenoid_center_PA anatomy.glenoid_center_IS anatomy.glenoid_center_ML anatomy.glenoid_Version anatomy.glenoid_Inclination 1 ''};
-                fieldNames = {'CT_id','radius','sphericity','biconcave','depth','width','height','version_ampl','version_orient','center_PA','center_IS','center_ML','version', 'inclination', 'version_2D', 'walch_class'};
-                upsert(conn, 'glenoid', fieldNames, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0], data);
-
-                % Update data in table 'humerus'
-                data = {sqlresult{1} 0.0 anatomy.humeral_head_radius anatomy.humerus_SHsubluxation_ampl anatomy.humerus_SHsubluxation_orient anatomy.humerus_GHsubluxation_ampl anatomy.humerus_GHsubluxation_orient 0.0 0.0};
-                fieldNames = {'CT_id','joint_radius','head_radius','SHsublux_ampl','SHsublux_orient','GHsublux_ampl','GHsublux_orient','SHsublux_2D','GHsublux_2D'};
-                upsert(conn, 'humerus', fieldNames, [1 0 0 0 0 0 0 0 0], data);
-
-                % Update data in table 'scapula'
-                data = {sqlresult{1} anatomy.scapula_CTangle anatomy.scapula_AI};
-                fieldNames = {'CT_id','CT_angle','AI'};
-                upsert(conn, 'scapula', fieldNames, [1 0 0], data);
-
-            else
-                error(['No entry for caseID ' caseID{i} ' was found in the ' conn.Instance ' MySQL database.']);
+            if DBaccess
+                sqlquery = ['SELECT CT_id FROM CT WHERE shoulder_id IN (SELECT shoulder_id FROM sCase WHERE folder_name = "' caseID{i} '")'];
+                curs = exec(conn,sqlquery);
+                if ~isempty(curs.Message)
+                    fprintf('\nMessage: %s\n', curs.Message);
+                end
+                        
+                curs=fetch(curs);
+                sqlresult = curs.Data;
+
+                if isnumeric(sqlresult{1})
+
+                    % Update data in table 'glenoid'
+                    data = {sqlresult{1} anatomy.glenoid_Radius anatomy.glenoid_Radius/anatomy.glenoid_radiusRMSE '' anatomy.glenoid_depth anatomy.glenoid_width anatomy.glenoid_height anatomy.glenoid_version_ampl anatomy.glenoid_version_orient anatomy.glenoid_center_PA anatomy.glenoid_center_IS anatomy.glenoid_center_ML anatomy.glenoid_Version anatomy.glenoid_Inclination 1 ''};
+                    fieldNames = {'CT_id','radius','sphericity','biconcave','depth','width','height','version_ampl','version_orient','center_PA','center_IS','center_ML','version', 'inclination', 'version_2D', 'walch_class'};
+                    upsert(conn, 'glenoid', fieldNames, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0], data);
+
+                    % Update data in table 'humerus'
+                    data = {sqlresult{1} 0.0 anatomy.humeral_head_radius anatomy.humerus_SHsubluxation_ampl anatomy.humerus_SHsubluxation_orient anatomy.humerus_GHsubluxation_ampl anatomy.humerus_GHsubluxation_orient 0.0 0.0};
+                    fieldNames = {'CT_id','joint_radius','head_radius','SHsublux_ampl','SHsublux_orient','GHsublux_ampl','GHsublux_orient','SHsublux_2D','GHsublux_2D'};
+                    upsert(conn, 'humerus', fieldNames, [1 0 0 0 0 0 0 0 0], data);
+    
+                    % Update data in table 'scapula'
+                    data = {sqlresult{1} anatomy.scapula_CTangle anatomy.scapula_AI};
+                    fieldNames = {'CT_id','CT_angle','AI'};
+                    upsert(conn, 'scapula', fieldNames, [1 0 0], data);
+
+                else
+                    error(['No entry for caseID ' caseID{i} ' was found in the ' conn.Instance ' MySQL database.']);
+                end
             end
-
         else
             error(['Amira files missing for CaseID: ' caseID{i}]);
         end
 
     end
 
     %  Close mySQL database
-    closeSQL(conn);
+    if DBaccess
+        closeSQL(conn);
+    end
 
 end
diff --git a/listSCase.m b/listSCase.m
new file mode 100755
index 0000000..bbcf5f6
--- /dev/null
+++ b/listSCase.m
@@ -0,0 +1,106 @@
+function SCaseList = listSCase(varargin)
+%LISTSCASE output list of sCases (id and path), accordind to inputList and filter
+%arguments. Each sCase directory should at least contain at directory with a
+%name staring with 'CT' and ending with '1', containing a 'dicom'
+%subdirectory.
+%
+%   The optional input argument inputList if an array of sCases. If empty, sCaseList
+%   will search for all sCases. The optional input argument filter can restrict le list. 
+
+% Syntax: sCaseList = listsCase(varargin)
+% 
+% Input: 
+%   inputArg: can be nothing (=all), 'N', 'P', or a sCaseId
+%
+% Output:
+%   outputArg: 1 if ok and 0 otherwise
+% 
+% Example: listSCase, listSCase('P315')
+% Other M-files required:
+% Subfunctions:
+% See also:
+% 
+% Author: Alexandre Terrier
+% EPFL-LBO
+% Creation date: 2018-07-01
+% Revision date: 2018-08-09
+% 
+% TO DO:
+% When varargin = a sCaseId
+
+
+%% Initialisation of paths
+% Should be replaced by a function
+dataDir = '../../../data'; % location of data
+logDir = 'log'; % log file in xls folder
+
+%% Open log file
+
+logFile = [logDir '/listSCase.log'];
+logFileId = fopen(logFile, 'w');
+fprintf(logFileId, 'listSCase log file');
+fprintf(logFileId, '\nDate: ');
+fprintf(logFileId, datestr(datetime));
+fprintf(logFileId, '\nGet list of all SCase with a valid dicom folder');
+
+%%
+% Struct contains id and associated directory
+SCaseList = struct(...
+    'id'       ,[],...
+    'dir'      ,[],...
+    'comment'  ,[]...   % Might be used later
+    );
+sCaseListIdx = 0;
+
+%% Loop over all cases
+
+% At the root of dataDir, two directories should be considered:
+% N (normal) and P (pathological).
+dirL0Array = ['N';'P'];
+if strcmp(varargin,'N')
+    dirL0Array = 'N';
+elseif strcmp(varargin,'P')
+    dirL0Array = 'P';
+end
+
+for dirL0idx = 1:length(dirL0Array) % N or P
+    for dirL1idx = 0:9 % Hunderts
+        for dirL2idx = 0:9 % Tens
+            dirL0Str = dirL0Array(dirL0idx);
+            dirL1Str = num2str(dirL1idx);
+            dirL2Str = num2str(dirL2idx);
+            dirTens  = [dataDir '/' dirL0Str '/' dirL1Str '/' dirL2Str];
+            dirTensList  = dir(dirTens); % List all directories in tens directory
+            dirTensList  = dirTensList(~ismember({dirTensList.name},{'.','..'})); % Remove . and ..
+            % We might add here a filter (only diretories, {N/P} followed
+            % by digits
+            for dirL3idx = 1:numel(dirTensList)
+                sCaseDirName = dirTensList(dirL3idx).name; % Root directory of this case
+                sCaseDirPath = dirTensList(dirL3idx).folder;
+                % Check that this sCase directory contains a directory with the
+                % name [CT?*-1] and with a dicom directory in it.
+                sCaseID = strtok(sCaseDirName, '-'); % Get chars before '-'
+                % Listing of (CT) folders in sCaseDirPath staring with CT
+                % and ending with '-1'
+                CTdir = dir([sCaseDirPath '/' sCaseDirName '/CT*-1']);
+                if ~isempty(CTdir)
+                    CTdirPath = [CTdir.folder '/' CTdir.name];
+                    dicomDirName = [CTdirPath '/dicom'];
+                    if exist(dicomDirName, 'dir')
+                        sCaseListIdx = sCaseListIdx + 1;
+                        SCaseList(sCaseListIdx).id = sCaseID;
+                        SCaseList(sCaseListIdx).dir = CTdirPath;
+                        % We might get here the readme.txt if present
+                    else
+                        fprintf(logFileId, ['\n', sCaseID, ' has no valid dicom directory']);
+                    end
+                else
+                    fprintf(logFileId, ['\n', sCaseID, ' has no CT*-1 directory']);
+                end
+            end
+        end
+    end
+end
+fprintf(logFileId, ['\nNumber of SCase: ' num2str(length(SCaseList))]);
+end
+
diff --git a/muscles/degenerationAll.m b/muscles/degenerationAll.m
index 3dccfb7..033bf2d 100644
--- a/muscles/degenerationAll.m
+++ b/muscles/degenerationAll.m
@@ -1,158 +1,166 @@
 function muscles_results = degenerationAll(caseID, musclesList)
 %%DEGENERATIONALL returns a structure array containing muscle 
 % measurements for each case ID entered as input and fills the Shoulder 
 % database with muscle degeneration measurements
 %
 % This function fills the Shoulder database with the muscle degeneration
 % values of the case given as input. The measurements are stored in a XLS
 % file and in the MySQL database.
 %
 % USE: degenerationAll(caseID, musclesList)
 %
 % INPUT     caseID: Case IDs for which you want to compute muscle degeneration.
 %           Cell array that stores the case IDs as a char in the form 'P###' or
 %           'N###' (starts with "P" or "N" followed by 1 to 3 digits).
 %
 %           musclesList: List of muscles which you want to compute muscle
 %           degeneration. Cell array that stores each msucle name as a char.
 %
 % OUTPUT    muscles_results: structure array containing fields 'Case_ID', 
 %           and fields "muscle name" with muscle measurements for the 
-%           current patient: 'Stot', 'Satrophy', 'Sinfiltration', 
-%           'Sosteochondroma' and 'Sdegeneration', 'atrophy',  
+%           current patient: areas : 'Stot', 'Satrophy', 'Sinfiltration', 
+%           'Sosteochondroma' and 'Sdegeneration', and ratios: 'atrophy',  
 %           'infiltration', 'osteochondroma' and 'degeneration'
 %
 % REMARKS   The muscle values (ratios) are saved in the muscles.xls file  
 %           and in the shoulder MySQL database.
 %
 % created with MATLAB ver.: 9.2.0.556344 (R2017a) on Windows 7
 % Author: EPFL-LBO-VMC
 % Date: 27-Jun-2017
 %
 
     CTDatabaseLocation = '../../../data'; % Location of the CT database
     XLSShoulderDatabaseLocation = '../../../data/Excel/'; % Location of the XLS ShoulderDatabase
 
     HUThresholdMuscle = 0; %#ok<NASGU>
     HUThresholdFat = 30; %#ok<NASGU>
     HUThresholdOsteochondroma = 166; %#ok<NASGU>
     
     % Check validity of input argument (cell array)
     validateattributes(caseID,{'cell'},{'nonempty'});
 
     % Check that imput case IDs are unique
     [caseIDNames,~,uniqueCaseID] = unique(caseID);
     countOfCaseID = hist(uniqueCaseID,unique(uniqueCaseID))';
     freqCaseIDs = struct('name',caseIDNames,'freq',num2cell(countOfCaseID));
     repeatedValues = find([freqCaseIDs.freq] > 1);
 
     if(~isempty(repeatedValues))
         warning('Input contains repeated Case ID: %s. Repeated occurences of that Case ID will be omitted. \n', freqCaseIDs(repeatedValues).name)
         caseID = caseIDNames';
     end
 
     %  open mySQL connection
     conn = openSQL();
 
     % Get the list of exisiting cases in XLS database
     [~,~,currDatabase] = xlsread([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls']);
     currDatabaseCaseIDlist = currDatabase(2:end,1);
     
     %Initialize structure array for results
     muscles_results = [];
     muscles_results = setfield(muscles_results,{length(caseID),1},'Case_ID',[]);
     
     % Compute muscle degeneration values and fill database
     for i = 1:length(caseID)
 
         % Check validity of input arguments (char and format 'P###' or 'N###')
         validateattributes(caseID{i},{'char'},{'nonempty'});
 
         if (numel(regexp(caseID{i},'^[PN]\d{1,3}$')) == 0)
             error(['Invalid format of CaseID: ' caseID{i} '. CaseID must start with "P" or "N" and be followed by 1 to 3 digits.']);
         end
 
         % Find CT folder for the given Case ID
         levelDir1 = caseID{i}(1);
 
         if (length(caseID{i}(2:end)) < 2)
             levelDir2 = '0';
             levelDir3 = '0';
         elseif (length(caseID{i}(2:end)) < 3)
             levelDir2 = '0';
             levelDir3 = caseID{i}(2);
         else
             levelDir2 = caseID{i}(2);
             levelDir3 = caseID{i}(3);
         end
 
         FindCaseCTFolder = dir([CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' caseID{i} '*']);
         if (isempty(FindCaseCTFolder))
             error(['Missing CT directory for CaseID: ' caseID{i}]);
         end
 
         CaseID_IPP = FindCaseCTFolder.name;
         CaseCTFolder = [CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' CaseID_IPP];
 
         % Compute musle degeneration values and store in the XLS database
         if exist([CaseCTFolder '/CT-' CaseID_IPP '-1/muscles'],'dir') == 7 %if a folder "muscles" exist, the patient name is added to the "measured" column
 
             muscleDirectory = [CaseCTFolder '/CT-' CaseID_IPP '-1/muscles']; %#ok<NASGU>
 
             muscles_results(i).Case_ID = caseID{i}; 
 
             for j=1:length(musclesList)
                 
                 disp(caseID{i});
                 % Compute muscle degeneration values
                 eval([musclesList{j} '= muscleDegeneration(caseID{i}, musclesList{j},HUThresholdMuscle,HUThresholdFat,HUThresholdOsteochondroma,muscleDirectory);']);
                 eval(['muscles_results = setfield(muscles_results,{i,1}, musclesList{j}, ' musclesList{j} ');']);
 
             end
 
             muscles = cellstr(muscles_results(i).Case_ID)';
             for j = 1:length(musclesList)
                 eval(['muscles = horzcat(muscles, transpose(struct2cell(muscles_results(i).' musclesList{j} ')));']);
             end
 
             % Write to XLS file
             % Replace line if already existing or append
             if(any(strcmp(currDatabaseCaseIDlist, caseID{i})))
-                xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles(:,[1,7:10,16:19,25:28,34:37]),'Feuil1',['A' int2str(find(strcmp(currDatabaseCaseIDlist, caseID{i}))+1)]);
+                xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles(:,[1,8:11,18:21,28:31,38:41]),'Feuil1',['A' int2str(find(strcmp(currDatabaseCaseIDlist, caseID{i}))+1)]);
             else
-                xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles(:,[1,7:10,16:19,25:28,34:37]),'Feuil1',['A' int2str(length(currDatabaseCaseIDlist)+1+i)]);
+                xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles(:,[1,8:11,18:21,28:31,38:41]),'Feuil1',['A' int2str(length(currDatabaseCaseIDlist)+1+i)]);
             end
+            
+            %ORIGINAL but weird results coz shifted - correction by Yasmine in Aug2018
+%             if(any(strcmp(currDatabaseCaseIDlist, caseID{i})))
+%                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles_test.xls'],muscles(:,[1,7:10,16:19,25:28,34:37]),'Feuil1',['A' int2str(find(strcmp(currDatabaseCaseIDlist, caseID{i}))+1)]);
+%             else
+%                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles_test.xls'],muscles(:,[1,7:10,16:19,25:28,34:37]),'Feuil1',['A' int2str(length(currDatabaseCaseIDlist)+1+i)]);
+%             end
+            
             % Write to MySQL database
             % Replace line if already existing or issue warning message
             sqlquery = ['SELECT CT_id FROM CT WHERE shoulder_id IN (SELECT shoulder_id FROM sCase WHERE folder_name = "' caseID{i} '")'];
             curs = exec(conn,sqlquery);
             if ~isempty(curs.Message)
                 fprintf('\nMessage: %s\n', curs.Message);
             end
 
             curs=fetch(curs);
             sqlresult = curs.Data;
 
             if isnumeric(sqlresult{1})
 
                 % Update data in table 'muscle'
                 data = {sqlresult{1} ... 
                     muscles_results(i).SS.atrophy muscles_results(i).SS.infiltration muscles_results(i).SS.osteochondroma muscles_results(i).SS.degeneration ...
                     muscles_results(i).IS.atrophy muscles_results(i).IS.infiltration muscles_results(i).IS.osteochondroma muscles_results(i).IS.degeneration ...
                     muscles_results(i).SC.atrophy muscles_results(i).SC.infiltration muscles_results(i).SC.osteochondroma muscles_results(i).SC.degeneration ...
                     muscles_results(i).TM.atrophy muscles_results(i).TM.infiltration muscles_results(i).TM.osteochondroma muscles_results(i).TM.degeneration};
                 fieldNames = {'CT_id','SSA','SSI','SSO','SSD','ISA','ISI','ISO','ISD','SCA','SCI','SCO','SCD', 'TMA', 'TMI', 'TMO', 'TMD'};
                 upsert(conn, 'muscle', fieldNames, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0], data);
 
             else
                 error(['No entry for caseID ' caseID{i} ' was found in the ' conn.Instance ' MySQL database.']);
             end
                         
         else
             error(['"muscles" folder missing for CaseID: ' caseID{i}]);
         end
 
     end
 
    
 end
diff --git a/muscles/degenerationAll3D.m b/muscles/degenerationAll3D.m
index 3bc8704..a849f8f 100644
--- a/muscles/degenerationAll3D.m
+++ b/muscles/degenerationAll3D.m
@@ -1,192 +1,193 @@
 function muscles_results = degenerationAll3D(caseID, musclesList)
 %%DEGENERATIONALL3D returns a structure array containing muscle
 % measurements for each case ID entered as input and fills the Shoulder
 % database with muscle degeneration measurements
 %
 % This function fills the Shoulder database with the muscle degeneration
 % values of the case given as input, for all CT slices available.
 % The measurements are stored in a XLS file and in the MySQL database.
 %
 % USE: degenerationAll3D(caseID, musclesList)
 %
 % INPUT     caseID: Case IDs for which you want to compute muscle degeneration.
 %           Cell array that stores the case IDs as a char in the form 'P###' or
 %           'N###' (starts with "P" or "N" followed by 1 to 3 digits).
 %
 %           musclesList: List of muscles which you want to compute muscle
 %           degeneration. Cell array that stores each msucle name as a char.
 %
 % OUTPUT    muscles_results: structure array containing fields 'Case_ID',
 %           and structure array field "muscle name" with muscle
 %           measurements for each slice for the current patient: 'Stot',
 %           'Satrophy', 'Sinfiltration', 'Sosteochondroma' and
 %           'Sdegeneration', 'atrophy', 'infiltration', 'osteochondroma'
 %           and 'degeneration'
 %
 % REMARKS   The muscle values (ratios) for all available CT slices are
 %           saved in the muscles.xls file and in the shoulder MySQL
 %           database.
 %
 % created with MATLAB ver.: 9.2.0.556344 (R2017a) on Windows 7
 % Author: EPFL-LBO-VMC
 % Date: 29-kun-2017
 %
 
     CTDatabaseLocation = '../../../data'; % Location of the CT database
-%     XLSShoulderDatabaseLocation = 'C:/Users/vmalfroy/Dropbox/shoulderDatabase/'; % Location of the XLS ShoulderDatabase
+
+    %     XLSShoulderDatabaseLocation = 'C:/Users/vmalfroy/Dropbox/shoulderDatabase/'; % Location of the XLS ShoulderDatabase
     XLSShoulderDatabaseLocation = '../../../data/Excel/'; % Location of the XLS ShoulderDatabase
     HUThresholdMuscle = 0; %#ok<NASGU>
     HUThresholdFat = 30; %#ok<NASGU>
     HUThresholdOsteochondroma = 166; %#ok<NASGU>
 
     % Check validity of input argument (cell array)
     validateattributes(caseID,{'cell'},{'nonempty'});
 
     % Check that imput case IDs are unique
     [caseIDNames,~,uniqueCaseID] = unique(caseID);
     countOfCaseID = hist(uniqueCaseID,unique(uniqueCaseID))';
     freqCaseIDs = struct('name',caseIDNames,'freq',num2cell(countOfCaseID));
     repeatedValues = find([freqCaseIDs.freq] > 1);
 
     if(~isempty(repeatedValues))
         warning('Input contains repeated Case ID: %s. Repeated occurences of that Case ID will be omitted. \n', freqCaseIDs(repeatedValues).name)
         caseID = caseIDNames';
     end
 
     %  open mySQL connection
     conn = openSQL();
 
     % Get the list of exisiting cases in XLS database
     [~,~,currDatabase] = xlsread([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls']);
     currDatabaseCaseIDlist = currDatabase(2:end,1);
 
     %Initialize structure array for results
     muscles_results = [];
     muscles_results = setfield(muscles_results,{length(caseID),1},'Case_ID',[]);
 
     % Compute muscle degeneration values and fill database
     for i = 1:length(caseID)
 
         % Check validity of input arguments (char and format 'P###' or 'N###')
         validateattributes(caseID{i},{'char'},{'nonempty'});
 
         if (numel(regexp(caseID{i},'^[PN]\d{1,3}$')) == 0)
             error(['Invalid format of CaseID: ' caseID{i} '. CaseID must start with "P" or "N" and be followed by 1 to 3 digits.']);
         end
 
         % Find CT folder for the given Case ID
         levelDir1 = caseID{i}(1);
 
         if (length(caseID{i}(2:end)) < 2)
             levelDir2 = '0';
             levelDir3 = '0';
         elseif (length(caseID{i}(2:end)) < 3)
             levelDir2 = '0';
             levelDir3 = caseID{i}(2);
         else
             levelDir2 = caseID{i}(2);
             levelDir3 = caseID{i}(3);
         end
 
         FindCaseCTFolder = dir([CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' caseID{i} '*']);
         if (isempty(FindCaseCTFolder))
             error(['Missing CT directory for CaseID: ' caseID{i}]);
         end
 
         CaseID_IPP = FindCaseCTFolder.name;
         CaseCTFolder = [CTDatabaseLocation '/' levelDir1 '/' levelDir2 '/' levelDir3 '/' CaseID_IPP];
 
         % Compute musle degeneration values and store in the XLS database
         if exist([CaseCTFolder '/CT-' CaseID_IPP '-1/muscles'],'dir') == 7 %if a folder "muscles" exist, the patient name is added to the "measured" column
 
             muscleDirectory = [CaseCTFolder '/CT-' CaseID_IPP '-1/muscles'];
 
             slicesList = dir(muscleDirectory);
             slicesList = slicesList(3:end);
 
             muscles_results(i).Case_ID = caseID{i};
 
             field2sort = 'SliceNumber'; %#ok<NASGU>
 
             for j=1:length(slicesList)
 
                 for k=1:length(musclesList)
 
                     disp(caseID{i});
                     
                     % Compute muscle degeneration values
                     eval([musclesList{k} '= muscleDegeneration3D(caseID{i}, musclesList{k},HUThresholdMuscle,HUThresholdFat,HUThresholdOsteochondroma,muscleDirectory);']);
                     eval(['muscles_results = setfield(muscles_results,{i,1}, musclesList{k}, ' musclesList{k} ');']);
                     
                     % Sort muscle field by slice number
                     eval(['[~,I] = sort(arrayfun (@(x) x.(field2sort), muscles_results(i,1).' musclesList{k} '));']);
                     eval(['muscles_results(i,1).' musclesList{k} ' =  muscles_results(i,1).' musclesList{k} '(I);']);
                 
                 end
 
             end
 
             muscles = cellstr(muscles_results(i).CT_id)';
             
             
             for k = 1:length(musclesList)               
                 eval(['muscles = horzcat(muscles, transpose(struct2cell(muscles_results.' musclesList{k} ')));']);
             end
 
             % Write to XLS file
             % Replace line if already existing or append
             if(any(strcmp(currDatabaseCaseIDlist, caseID{i})))
                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles,'Feuil1',['A' int2str(find(strcmp(currDatabaseCaseIDlist, caseID{i}))+1)]);
             else
                 xlswrite([XLSShoulderDatabaseLocation 'xlsFromMatlab/muscles.xls'],muscles,'Feuil1',['A' int2str(length(currDatabaseCaseIDlist)+1+i)]);
             end
             % Write to MySQL database
             % Replace line if already existing or issue warning message
             sqlquery = ['SELECT CT_id FROM CT WHERE shoulder_id IN (SELECT shoulder_id FROM sCase WHERE folder_name = "' caseID{i} '")'];
             curs = exec(conn,sqlquery);
             if ~isempty(curs.Message)
                 fprintf('\nMessage: %s\n', curs.Message);
             end
 
             curs=fetch(curs);
             sqlresult = curs.Data;
 
             if isnumeric(sqlresult{1})
 
                 SS_atrophy = sum([muscles_results(i).SS(:).Satrophy])/sum([muscles_results(i).SS(:).Stot]);
                 SS_infiltration = sum([muscles_results(i).SS(:).Sinfiltration])/sum([muscles_results(i).SS(:).Stot]);
                 SS_osteochondroma = sum([muscles_results(i).SS(:).Sosteochondroma])/sum([muscles_results(i).SS(:).Stot]);
                 SS_degeneration = (sum([muscles_results(i).SS(:).Stot]) - sum([muscles_results(i).SS(:).Sdegeneration]))/sum([muscles_results(i).SS(:).Stot]);
   
                 IS_atrophy = sum([muscles_results(i).IS(:).Satrophy])/sum([muscles_results(i).IS(:).Stot]);
                 IS_infiltration = sum([muscles_results(i).IS(:).Sinfiltration])/sum([muscles_results(i).IS(:).Stot]);
                 IS_osteochondroma = sum([muscles_results(i).IS(:).Sosteochondroma])/sum([muscles_results(i).IS(:).Stot]);
                 IS_degeneration = (sum([muscles_results(i).IS(:).Stot]) - sum([muscles_results(i).IS(:).Sdegeneration]))/sum([muscles_results(i).IS(:).Stot]);
                 
                 SC_atrophy = sum([muscles_results(i).SC(:).Satrophy])/sum([muscles_results(i).SC(:).Stot]);
                 SC_infiltration = sum([muscles_results(i).SC(:).Sinfiltration])/sum([muscles_results(i).SC(:).Stot]);
                 SC_osteochondroma = sum([muscles_results(i).SC(:).Sosteochondroma])/sum([muscles_results(i).SC(:).Stot]);
                 SC_degeneration = (sum([muscles_results(i).SC(:).Stot]) - sum([muscles_results(i).SC(:).Sdegeneration]))/sum([muscles_results(i).SC(:).Stot]);
                 
                 TM_atrophy = sum([muscles_results(i).TM(:).Satrophy])/sum([muscles_results(i).TM(:).Stot]);
                 TM_infiltration = sum([muscles_results(i).TM(:).Sinfiltration])/sum([muscles_results(i).TM(:).Stot]);
                 TM_osteochondroma = sum([muscles_results(i).TM(:).Sosteochondroma])/sum([muscles_results(i).TM(:).Stot]);
                 TM_degeneration = (sum([muscles_results(i).TM(:).Stot]) - sum([muscles_results(i).TM(:).Sdegeneration]))/sum([muscles_results(i).TM(:).Stot]);
        
                 % Update data in table 'muscle'
                 data = {sqlresult{1} SS_atrophy SS_infiltration SS_osteochondroma SS_degeneration IS_atrophy IS_infiltration IS_osteochondroma IS_degeneration SC_atrophy SC_infiltration SC_osteochondroma SC_degeneration TM_atrophy TM_infiltration TM_osteochondroma TM_degeneration};
                 fieldNames = {'CT_id','SSA','SSI','SSO','SSD','ISA','ISI','ISO','ISD','SCA','SCI','SCO','SCD', 'TMA', 'TMI', 'TMO', 'TMD'};
                 upsert(conn, 'muscle', fieldNames, [1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0], data);
 
             else
                 error(['No entry for caseID ' caseID{i} ' was found in the ' conn.Instance ' MySQL database.']);
             end
 
         else
             error(['"muscles" folder missing for CaseID: ' caseID{i}]);
         end
 
     end
 
 
 end
diff --git a/muscles/muscleContours.m b/muscles/muscleContours.m
index 0e4fdcd..e4e18f8 100644
--- a/muscles/muscleContours.m
+++ b/muscles/muscleContours.m
@@ -1,148 +1,148 @@
 %% muscleContours
 
 % This script allows to select the contours on Tiff images of the four  
 % cuff rotator muscles: Supraspinatus (SS), Infraspinatus (IS), Subscapularis (SC)
 % and Teres Minor (TM).
 % 
 
 % INPUT: 
 % subject: ### subject number
 % type: 'normal' or 'pathologic'
 
 % Author: EPFL-LBO
 % Date: 2016-09-05
 %%
 
 function [] = muscleContours(subject,type)
 
 directory = '../../../data'; % Location of the CT database
 location = 'P';
-CTn = sprintf('P%03d',subject);
+CTn = sprintf('P%03d',subject)
 
-listDirInCurrDir = dir([directory '/' CTn(1) '/' CTn(2) '/' CTn(3)]);
+listDirInCurrDir = dir([directory '/' CTn(1) '/' CTn(2) '/' CTn(3)])
 listDirInCurrDir = listDirInCurrDir(3:end);
 
 for i=1:length(listDirInCurrDir)
         if ~isempty(regexp(listDirInCurrDir(i).name,['^[PN]\d{0,2}[' CTn(4) '][-]\w*'], 'once'))
               
             listDirInCaseFolder = dir([listDirInCurrDir(i).folder '\' listDirInCurrDir(i).name]);
             listDirInCaseFolder = listDirInCaseFolder(3:end);
             
             % In all CT folders in the current Case Folder, find CT #1
             for j=1:length(listDirInCaseFolder)
                 if ~isempty(regexp(listDirInCaseFolder(j).name,'^[C][T][-][PN]\d{1,3}[-]\d*[-][1]', 'once'))
                     outputPath = [listDirInCaseFolder(j).folder '\' listDirInCaseFolder(j).name];
                 end
             end            
         end
 end
 
 
 pathologic_list = dir(sprintf('%s/%s/P*', directory, location)); %Create the case list
 for i = 1:1:numel(pathologic_list)  % Loop over the list of patients and extract the complete name of the specific patient
     name = pathologic_list(i).name;
     t = strfind(name,CTn);
     if t == 1;
         patient_name = name;
     end
 end
 
 finaldirectory = [outputPath '\amira\*tif'];
 muscleName = cellstr(['SS';'IS';'SC';'TM']);
 
 close all; % Close all figure windows except those created by imtool.
 imtool close all; % Close all figure windows created by imtool.
 workspace; % Make sure the workspace panel is showing.
 fontSize = 16;
 imageWindow = [-100 200];
 
 % Choose file
 [inputFile, inputPath] = uigetfile(finaldirectory,'Select the image for muscle atrophy measurement');
 whereisp = strfind(inputFile, 'P');
 if isempty(whereisp)
     whereisp = strfind(inputFile, 'N');
     if isempty(whereisp)
         return
     end
 end
 subject = inputFile(whereisp:whereisp+3);
 
 % Read gray scale image.
 inputImage = imread(sprintf('%s%s',inputPath,inputFile));
 
 if size(inputImage,3) == 3
 % Calculate the monochrome luminance by combining the RGB values according to the NTSC standard
     inputImage = 0.2989*inputImage(:,:,1)...
         +0.5870*inputImage(:,:,2)...
         +0.1140*inputImage(:,:,3);
 end
 
 for i=1:4;
     
     finished = 'Redo';
     while strcmp('Redo',finished)
         msg = 0; % 0 = no message
         %         clc; % Clear command window.
         close all; % Close all figure windows except those created by imtool.
         imtool close all; % Close all figure windows created by imtool.
         imshow(inputImage,imageWindow);
         imagetitle = sprintf('Select %s ideal section', char(muscleName(i)));
         title(imagetitle, 'FontSize', fontSize);
         set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
         
         if msg == 1
             message = sprintf('Select the ideal section of your muscle in the image. \nLeft click to place points and draw a polygon. Click again on the first point to close the shape.\nYou can move the points once the line is closed. Double-click or right-click and "Create Mask" to finish.');
             msg = 0;
             uiwait(msgbox(message));
         end
         
         binaryImage1 = roipoly();
         
         if isempty(binaryImage1)
             error('Invalid selection or interface closed')
         end
         
         binaryImage1 = bwmorph(binaryImage1,'majority');
         % Get coordinates of the boundary of the polygonal drawn region.
         structBoundaries1 = bwboundaries(binaryImage1);
         xy1 = structBoundaries1{1}; % Get n by 2 array of x,y coordinates.
         x1 = xy1(:, 2); % Columns.
         y1 = xy1(:, 1); % Rows.
         hold on
         plot(x1, y1, 'color','green');
         
         % check contours
         finished = questdlg('Please verify the contours. Press "Save" if they are satisfactory or "Redo" if you would like to redo the selection. ',...
             'Contour verification','Save','Redo','Save');
     end
     
     %% Report results
     % mkdir('./results',subject);
     if exist([outputPath '\muscles'],'dir') == 0
         mkdir(outputPath,'muscles');
     end
     
     outputPath2 = [outputPath '\muscles\' char(muscleName(i)) '\'];
     mkdir(outputPath2);
     outputFile = sprintf('ssContours%s',subject);
     
     % Check for overwriting
     overwrite = exist(sprintf('%s%s.mat',outputPath2,outputFile), 'file');
     switch overwrite
         case 2
             isoverwrite = questdlg('Overwrite existing file ?', 'Warning', 'No', 'Yes', 'No');
             switch isoverwrite
                 case 'Yes'
                     save(sprintf('%s%s.mat',outputPath2,outputFile), 'inputImage','binaryImage1')
                     saveas(gcf,sprintf('%s%s.tif',outputPath2,outputFile))
                     fprintf('%s%s.mat\n%s%s.tif\n',outputPath2,outputFile,outputPath2,outputFile);
             end
             
         case 0
             save(sprintf('%s%s.mat',outputPath2,outputFile), 'inputImage','binaryImage1')
             saveas(gcf,sprintf('%s%s.tif',outputPath2,outputFile))
             fprintf('%s%s.mat\n%s%s.tif\n',outputPath2,outputFile,outputPath2,outputFile);
     end
 end
 close all;
 end
\ No newline at end of file