diff --git a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
index f760886..2fd2d24 100644
--- a/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
+++ b/ShoulderCase/@MusclesContainer/createRotatorCuffSlicesNathan.m
@@ -1,373 +1,368 @@
 function createRotatorCuffSlicesNathan(obj,sliceName)
   % This function and its results might be reffered as 
   % Nathan algorithm and Nathan slices in some studies.
   %
   %One needs to create two empty files in the current directory :
   %data_degenerescence.mat (a table) and Didnotwork.mat (a
   %string)
   %
   % This script calculates the degeneration of the 4 muscles
   % of the rotator cuff : SS, SC, IS, TM. It is composed of four
   % different steps :
   % - Load CT images (DICOM files) in order to create a 3D matrix
   % consituted with HU (Hounsfield unit) values
   % - Create a slice of this 3D matrix along a patient-specific
   % plane. Do image processing and set saving format usable by
   % automatic segmentation code (Python's code from ARTORG team
   % (UNIBE))
   % - Execute automatic segmentation code on the new created
   % images
   % - Apply degeneration code (from Laboratoire de biom�canique
   % orthop�dique (EPFL)) on new images with automatic
   % segmentation and on image with manual segmentation (from CHUV
   % radiologist). Save degeneration values under the "muscles"
   % properties of each case in the LBO database
   %
   % Some "try-catch" are present in order to catch the errors
   % while keeping running the script on the rest of the database.
   % They are summarized underneath :
   % (DicomReadVolume) : do not manage to load the dicoms
   % Nothing : Automatic segmentation does not work (generally for TM)
   % * : slice function does not work
   % ** : does not want to crop image
   % *** : black image, don't want to save png
 
   % Main code starts here !
 
   %% Load dicom
   SCase = obj.shoulder.SCase;
   dataCTPath = SCase.dataCTPath;
 
   % Create "muscles" folder (and its 4 subfolders) if it does not exist
   folder_muscle=[dataCTPath '/muscles'];
 
   if isfolder(folder_muscle)
   else
   mkdir(dataCTPath,'muscles');
   mkdir(folder_muscle,'IS');
   mkdir(folder_muscle,'SS');
   mkdir(folder_muscle,'SC');
   mkdir(folder_muscle,'TM');
   end
 
   % Use folder "CT-2" (smoothed kernel) if it exists
   dataCTPath2=dataCTPath;
   dataCTPath2(1,end)='2';
 
   if isfolder(dataCTPath2)
       dataCTPath0=dataCTPath2;
   else
       dataCTPath0=dataCTPath;
   end
 
   dicomFolder = [dataCTPath0 '/dicom'];
 
   %Extract the images as a 4D matrix where spatial is the
   %positions of the Pixels in the CT coordinate system
 
   % Return the case's number if it does not manage to load the DICOM files
   try
       [V,spatial,dim] = dicomreadVolume(dicomFolder);
   catch
       outputArg=1;
-      fprintf(logFileID, '\n|      ERROR: dicomreadVolume()');
       return;
   end
 
   %V needs to be 3D matrix containing the value in HU of the
   %corresponding voxel
   V = squeeze(V) ;
 
   %Read the information from the middle dicom image (all image
   %informations are the same for all dicom files) to avoid a problem
   %if other files are located at the beginning or end of the
   %folder
   dicomList = dir(dicomFolder) ;
   size_list = size(dicomList) ;
   dicom_information = dicominfo([dicomFolder '/' dicomList(round(size_list(1)/2)).name]);
 
   %Take into account that the data will need to be linearly transformed
   %from stored space to memory to come back to true HU units afterwards:
   Rescale_slope = dicom_information.RescaleSlope;
   Rescale_intercept = dicom_information.RescaleIntercept;
 
   %Convert to double type to be able to multiply it with a int64
   %matrix afterwards
   V = double(V);
 
   % To filter hard kernel images
 
   %             if isfolder(dataCTPath2)
   %             else
   %             V = imgaussfilt3(V,1) ;
   %             end
 
   size_V = size(V);
 
   %Extract position of each first pixel of each image (upper left corner)
   PatientPositions = spatial.PatientPositions ;
 
   %Get the distance between each pixel in neighboring rows and
   %columns
   PixelSpacings = spatial.PixelSpacings;
 
 
   % Get local coordinate system from parent scapula object(glenoid center)
   origin = obj.shoulder.scapula.coordSys.origin; % do not mix up with origin_image, "origin" is the origin of the Scapula Coordinate System
   xAxis = obj.shoulder.scapula.coordSys.PA;
   yAxis = obj.shoulder.scapula.coordSys.IS;
   zAxis = obj.shoulder.scapula.coordSys.ML;
 
 
   %% Slice from the CT data set
 
   %Get the position of the upper left pixel of the first image
   %and calculate the position of each pixel in the x,y
   %and z direction starting from the "origin" pixel
   origin_image = PatientPositions(1,:); % coordinates of the first pixel
   x_max = origin_image(1)+PixelSpacings(1,1)*(size_V(1)-1);
   y_max = origin_image(2)+PixelSpacings(1,2)*(size_V(2)-1);
   z_max = PatientPositions(end,3);
 
   %Calculate the coefficient of the linear transformation (CT-scan
   %coordinate system (x,y,z) to images coordinates system (i,j,k))
   coefficients_i = polyfit([1 size_V(1)],[origin_image(1) , x_max],1);
   a_i = coefficients_i(1);
   b_i = coefficients_i(2);
 
   coefficients_j = polyfit([1 size_V(2)],[origin_image(2) , y_max],1);
   a_j = coefficients_j(1);
   b_j = coefficients_j(2);
 
   coefficients_k = polyfit([1 size_V(3)],[origin_image(3) , z_max],1);
   a_k = coefficients_k(1);
   b_k = coefficients_k(2);
 
   Pixel_pos_x=(1:size_V(1))*a_i+b_i; % Create the position vector along x-axis
   Pixel_pos_y=(1:size_V(2))*a_j+b_j;
   Pixel_pos_z=(1:size_V(3))*a_k+b_k;
 
   PixSize_z=(z_max-origin_image(3))/(size_V(3)-1); % Space between 2 slices
 
   % The oblique slice to be created is perpendicular to the z
   % axis of the scapula CS and passes through the origin point.
   % We will apply the matlab function "slice", to create the
   % oblique slice. One must first create a slice parallel to 2
   % system's axis. Here we take a plane parallel to the (x,y)
   % plane, passing through the point "origin". Then, we need to
   % rotate this plane so that it becomes normal to the zAxis. V1
   % is the vector normal to the horizontal slice and v2 the
   % vector normal to the final oblique slice, parallel to the
   % zAxis. We need to rotate around a vector v3 and from an angle
   % given by "vrrotvec(v1,v2)".
   v1=[0,0,1];
   v2=zAxis/norm(zAxis);
   rotation_info=vrrotvec(v1,v2);
   v3=rotation_info(1:3);
   theta=rotation_info(4)/pi*180; % Convert it in degrees
 
   % Create the plane to be rotated
   [x,y,z] = meshgrid(Pixel_pos_x,Pixel_pos_y,Pixel_pos_z);
   s=1000; % Extension factor of the slice in order to cut the whole matrix.Set arbitrarily but sufficiently big so that it works for every patient
   Pixel_pos_x_enlarged=[origin_image(1)+(-PixelSpacings(1,1)*s:PixelSpacings(1,1):-PixelSpacings(1,1)) Pixel_pos_x x_max+(PixelSpacings(1,1):PixelSpacings(1,1):PixelSpacings(1,1)*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged=[origin_image(2)+(-PixelSpacings(1,2)*s:PixelSpacings(1,2):-PixelSpacings(1,2)) Pixel_pos_y y_max+(PixelSpacings(1,2):PixelSpacings(1,2):PixelSpacings(1,2)*s)];
 
   figure('visible','off');
   hsp = surf(Pixel_pos_x_enlarged,Pixel_pos_y_enlarged,origin(3)*ones(length(Pixel_pos_y_enlarged),length((Pixel_pos_x_enlarged))));
 
   rotate(hsp,v3,theta,origin);
 
   xd_0 = hsp.XData; % x-coordinates of the plane's points
   yd_0 = hsp.YData; % y-coordinates of the plane's points
   zd_0 = hsp.ZData; % z-coordinates of the plane's points
 
   % Slice the 3D matrix of CT-scan
   try % Return case's number in case slice does not work
   figure('visible','off');
   hss=slice(x,y,z,V,xd_0,yd_0,zd_0,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
   catch
   outputArg=2;
-  fprintf(logFileID,'\n|       ERROR: slice()');
   return
   end
 
   HU_0 = get(hss,'CData'); % HU values at each point of the plane
 
   NaNplaces_0=find(isnan(HU_0));
 
   HU_1=HU_0*Rescale_slope+Rescale_intercept; % getting true HU values
   p_0=10e6; % value outside 3D matrix, set to a high distinguishable number
   HU_1(NaNplaces_0)=p_0;
 
   % If it's a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_1=flipud(HU_1');
   else
   HU_1=HU_1';
   end
 
   % Find the new pixel's sizes necessary to obtain a 1024*1024
   % shoulder slice. The largest diagonal's length is set to 1024
   % pixels
   [left_1,left_2]=find(HU_1~=p_0,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_1,2)~=p_0,1);
   right_r=right_1;
   right_c=size(HU_1,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_1)~=p_0,1);
   up_c=size(HU_1)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_1,3)~=p_0,1);
   down_r=size(HU_1,1)-down_2+1;
   down_c=down_1;
 
   % In case the slice does not cut the whole matrix (which
   % normally never happens), ask to double "s" parameter
 
   if or(or(right_c==size(HU_1,2),left_c==1),or(up_r==1,down_r==size(HU_1,1)))
-    fprintf(logFileID,'\n|       ERROR: the slice doesn''t cut throughout the 3D CT-matrix');
     outputArg = 3;
     return
   end
 
   % Find which diagonal of the matrix' effective slice is the largest
   w=right_c-left_c; % width
   h=down_r-up_r; % height
   if h>w % Set the largest diagonal at 1024 pixels
       l=h;
   else
       l=w;
   end
 
   % Define new pixel's size
   NewPixSize_x=PixelSpacings(1,1)*l/1023;
   NewPixSize_y=PixelSpacings(1,2)*l/1023;
 
   % New plane with the new pixel's sizes
   Pixel_pos_x_enlarged_2=[origin_image(1)+(-NewPixSize_x*s:NewPixSize_x:-NewPixSize_x) origin_image(1):NewPixSize_x:x_max x_max+(NewPixSize_x:NewPixSize_x:NewPixSize_x*s)]; % Create the position vector along x-axis, long enough to fill the diagonal
   Pixel_pos_y_enlarged_2=[origin_image(2)+(-NewPixSize_y*s:NewPixSize_y:-NewPixSize_y) origin_image(2):NewPixSize_y:y_max y_max+(NewPixSize_y:NewPixSize_y:NewPixSize_y*s)];
 
   figure('visible','off');
   hsp2 = surf(Pixel_pos_x_enlarged_2,Pixel_pos_y_enlarged_2,origin(3)*ones(length(Pixel_pos_y_enlarged_2),length((Pixel_pos_x_enlarged_2))));
   rotate(hsp2,v3,theta,origin);
 
   xd = hsp2.XData; % x-coordinates of the rotated plane's points.
   yd = hsp2.YData; % y-coordinates of the rotated plane's points.
   zd = hsp2.ZData; % z-coordinates of the rotated plane's points.
 
   % New slice with the good resolution (1024*1024)
   figure('visible','off');
   hss=slice(x,y,z,V,xd,yd,zd,'Linear');
   % Linear interpolation by default. 'Cubic' or 'Nearest'
 
   HU_2 = get(hss,'CData'); % HU values at each point of the plane
   NaNplaces=find(isnan(HU_2));
   HU_slice=HU_2*Rescale_slope+Rescale_intercept;
   HU_slice_8_bit=HU_slice;
   p_tif=0; % value outside 3D matrix, for the tif image
   p_png=85; % for the png image
   HU_slice(NaNplaces)=p_tif;
 
   % Scaling to fit the training set (for the png only, grayscale)
   HU_interval_png=[-100 160];
   lin_coef=[HU_interval_png(1) 1 ; HU_interval_png(2) 1]^(-1)*[0 255]';
 
   HU_slice_8_bit=lin_coef(1)*HU_slice_8_bit+lin_coef(2);
   HU_slice_8_bit(NaNplaces)=p_png;
 
   % If it a right shoulder or left one, has to be turned in a
   % particular way
   if obj.shoulder.side=='R'
   HU_slice=flipud(HU_slice');
   HU_slice_8_bit=flipud(HU_slice_8_bit');
   else
   HU_slice=HU_slice';
   HU_slice_8_bit=HU_slice_8_bit';
   end
 
   % Determine the location of the 4 corners of the slice
   [left_1,left_2]=find(HU_slice~=p_tif,1);
   left_r=left_1;
   left_c=left_2;
   [right_1,right_2]=find(rot90(HU_slice,2)~=p_tif,1);
   right_r=right_1;
   right_c=size(HU_slice,2)-right_2+1;
   [up_1,up_2]=find(rot90(HU_slice)~=p_tif,1);
   up_c=size(HU_slice,2)-up_1+1;
   up_r=up_2;
   [down_1,down_2]=find(rot90(HU_slice,3)~=p_tif,1);
   down_r=size(HU_slice,1)-down_2+1;
   down_c=down_1;
 
   Folder=[dataCTPath '/matlab'];
 
   % Crop the image
   try % If cropping does not work, return case's number
   HU_slice=HU_slice(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice_8_bit=HU_slice_8_bit(up_r:up_r+1023,left_c:left_c+1023);
 
   HU_slice=int16(HU_slice); % convert into 16-bit
   catch
-  fprintf(logFileID,'\n|       ERROR: crop');
   outputArg=4;
   return
   end
 
   % Create '.tif' file
   file_name_tif=[SCase.id '.tif'];
 
   t = Tiff(fullfile(Folder,file_name_tif),'w');
   setTag(t,'Photometric',Tiff.Photometric.MinIsBlack);
   setTag(t,'Compression',Tiff.Compression.None);
   setTag(t,'BitsPerSample',16);
   setTag(t,'SamplesPerPixel',1);
   setTag(t,'ImageLength',1024);
   setTag(t,'ImageWidth',1024);
   setTag(t,'PlanarConfiguration',Tiff.PlanarConfiguration.Chunky);
   setTag(t,'SampleFormat',Tiff.SampleFormat.Int);
   write(t,HU_slice);
   close(t);
 
   % Create 8 bit '.png' file and save it at the right places
   % (matlab folder and python's input folder)
   file_name_png=[SCase.id '.png'];
 
   HU_slice_8_bit=uint8(HU_slice_8_bit);
   figure('visible','off');
   image(HU_slice_8_bit);
   colormap(gray(256));
 
   try % return case's number if does not manage to write png in the "matlab" folder because it is a black image
   imwrite(HU_slice_8_bit,[Folder '/' file_name_png]);
   catch
-  fprintf(logFileID,'\n|       ERROR: black image');
   outputArg=5;
   return
   end
 
   pythonDir = ConfigFileExtractor.getVariable('pythonDir');
   Folder_input_py= [pythonDir '/input'];
   imwrite(HU_slice_8_bit,[Folder_input_py '/' file_name_png]); % write png file in the python input folder
   % imwrite(HU_slice_8_bit,['./rcseg/' file_name_png]); % write png file in local /rceg folder
 
   % Create .txt file with pixel dimensions
   file_name_info=[SCase.id '.info'];
   text1='# Amira Stacked Slices' ;
   text2=['pixelsize ', num2str(NewPixSize_x),' ',num2str(NewPixSize_y)];
   fileID = fopen([Folder '/' file_name_info],'w');
   fprintf(fileID,'%s\r\n\r\n',text1);
   fprintf(fileID,'%s\n',text2);
   fclose(fileID);
   close all;
 
 
   imwrite(HU_slice_8_bit,fullfile(obj.slicesDataPath,[sliceName '_ForSegmentation.png']));
   imageForMeasurements = HU_slice;
   save(fullfile(obj.slicesDataPath,[sliceName '_ForMeasurements.mat']),'imageForMeasurements');
   imagesPixelSpacings = [NewPixSize_x NewPixSize_y];
   save(fullfile(obj.slicesDataPath,[sliceName '_PixelSpacings.mat']),'imagesPixelSpacings');
 end